force_a4.f90 5.27 KB
Newer Older
1 2
MODULE m_force_a4
CONTAINS
3
  SUBROUTINE force_a4(atoms,sphhar,input,dimension,&
4 5 6 7 8 9 10 11 12 13 14 15
       &                    vr,&
       &                    force)
    !
    ! ************************************************************
    ! rhocore force contribution a la Rici et al.
    !
    ! ************************************************************
    !
    USE m_intgr, ONLY : intgr0,intgr3
    USE m_constants, ONLY : sfp_const
    USE m_differentiate,ONLY: difcub
    USE m_types
16
    USE m_cdn_io
17
    IMPLICIT NONE
18 19 20 21
    TYPE(t_input),INTENT(IN)     :: input
    TYPE(t_sphhar),INTENT(IN)    :: sphhar
    TYPE(t_atoms),INTENT(IN)     :: atoms
    TYPE(t_dimension),INTENT(IN) :: dimension
22 23
    !     ..
    !     .. Array Arguments ..
Daniel Wortmann's avatar
Daniel Wortmann committed
24 25
    REAL,    INTENT (IN) :: vr(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,input%jspins)
    REAL,    INTENT (INOUT) :: force(3,atoms%ntype,input%jspins)
26 27 28 29 30 31 32
    !     ..
    !     .. Local Scalars ..
    REAL a4_1,a4_2,qcore,s13,s23,w,xi
    INTEGER i,ir,jsp,lh,lindex,mem,mindex,n,nd,na
    !     ..
    !     .. Local Arrays ..
    COMPLEX forc_a4(3),gv(3),ycomp1(3,-1:1)
33 34
    REAL rhoaux(atoms%jmtd),rhoc(atoms%jmtd,atoms%ntype,input%jspins)
    REAL tec(atoms%ntype,input%jspins),qintc(atoms%ntype,input%jspins)
35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52
    !     ..
    !     .. Data statements ..
    COMPLEX,PARAMETER:: czero=(0.000,0.000)
    !     ..
    !
    !     set ylm related components
    !
    s13 = SQRT(1.0/3.0)
    s23 = SQRT(2.0/3.0)
    ycomp1(1,0) = czero
    ycomp1(2,0) = czero
    ycomp1(3,0) = CMPLX(2.0*s13,0.0)
    ycomp1(1,-1) = CMPLX(s23,0.0)
    ycomp1(2,-1) = CMPLX(0.0,-s23)
    ycomp1(3,-1) = czero
    ycomp1(1,1) = CMPLX(-s23,0.0)
    ycomp1(2,1) = CMPLX(0.0,-s23)
    ycomp1(3,1) = czero
53 54 55 56
    !     --->    read in core density
    CALL readCoreDensity(input,atoms,dimension,rhoc,tec,qintc)

    DO jsp = 1,input%jspins
57
       na = 1
58
       DO n = 1,atoms%ntype
59 60 61 62 63
          IF (atoms%l_geo(n)) THEN
             DO i = 1,3
                forc_a4(i) = czero
             END DO
             !
64
             CALL intgr0(rhoc(:,n,jsp),atoms%rmsh(1,n),atoms%dx(n),atoms%jri(n),qcore)
65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91
8000         FORMAT (' FORCE_A4: core charge=',1p,e16.8)
             !
             !
             nd = atoms%ntypsy(na)
             !
             !
             lh_loop: DO  lh = 1,sphhar%nlh(nd)
                lindex = sphhar%llh(lh,nd)
                IF (lindex.GT.1) EXIT lh_loop
                DO i = 1,3
                   gv(i) = czero
                END DO
                !
                !    sum over all m for particular symm. harmonic
                DO mem = 1,sphhar%nmem(lh,nd)
                   mindex = sphhar%mlh(mem,lh,nd)
                   DO i = 1,3
                      gv(i) = gv(i) + sphhar%clnu(mem,lh,nd)*ycomp1(i,mindex)
                   END DO
                END DO
                !
                !
                !     construct integrand rhocore*d/dr(v)*r**2
                !     note: rhocore is already multiplied by r**2 and srt(4.*pi)
                !     difcub performs analytic derivative of Lagrangian of 3rd order
                !
                xi = atoms%rmsh(1,n)
92
                rhoaux(1) = difcub(atoms%rmsh(1,n),vr(1,lh,n,jsp),xi)*rhoc(1,n,jsp)
93 94 95
                DO ir = 2,atoms%jri(n) - 2
                   xi = atoms%rmsh(ir,n)
                   rhoaux(ir) = difcub(atoms%rmsh(ir-1,n),&
96
                        &                                vr(ir-1,lh,n,jsp),xi) * rhoc(ir,n,jsp)
97 98 99 100 101
                END DO
                !
                ir = atoms%jri(n) - 1
                xi = atoms%rmsh(ir,n)
                rhoaux(ir) = difcub(atoms%rmsh(atoms%jri(n)-3,n),&
102
                     &                              vr(atoms%jri(n)-3,lh,n,jsp),xi)*rhoc(ir,n,jsp)
103 104 105 106
                !
                ir = atoms%jri(n)
                xi = atoms%rmsh(ir,n)
                rhoaux(ir) = difcub(atoms%rmsh(atoms%jri(n)-3,n),&
107
                     &                              vr(atoms%jri(n)-3,lh,n,jsp),xi)*rhoc(ir,n,jsp)
108 109 110 111 112 113 114
                CALL intgr3(rhoaux,atoms%rmsh(1,n),atoms%dx(n),atoms%jri(n),w)
                a4_1 = 0.5*w/sfp_const
                !
                !     construct integrand rhocore*v*r
                !     note: rhocore is already multiplied by r**2 and srt(4.*pi)
                !
                DO ir = 1,atoms%jri(n)
115
                   rhoaux(ir) = rhoc(ir,n,jsp)/atoms%rmsh(ir,n)*vr(ir,lh,n,jsp)
116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147
                END DO
                !
                CALL intgr3(rhoaux,atoms%rmsh(1,n),atoms%dx(n),atoms%jri(n),w)
                a4_2 = w/sfp_const
                !
                DO i = 1,3
                   forc_a4(i) = forc_a4(i) - (a4_1+a4_2)*gv(i)
                END DO
                !
                !  lh loop ends
             ENDDO lh_loop
             !
             !     sum to existing forces
             !
             DO i = 1,3
                force(i,n,jsp) = force(i,n,jsp) + REAL(forc_a4(i))
             END DO
             !
             !     write result
             !
             WRITE (6,FMT=8010) n
             WRITE (6,FMT=8020) (forc_a4(i),i=1,3)
8010         FORMAT (' FORCES: EQUATION A4 FOR ATOM TYPE',i4)
8020         FORMAT (' FX_A4=',2f10.6,' FY_A4=',2f10.6,' FZ_A4=',2f10.6)
             ! type loop ends
          ENDIF
          na = na + atoms%neq(n)
       ENDDO
       ! spin loop ends
    ENDDO
  END SUBROUTINE force_a4
END MODULE m_force_a4