force_a8.F90 13.1 KB
Newer Older
1 2 3 4 5 6
MODULE m_forcea8
  ! ************************************************************
  ! Pulay 1st term  force contribution a la Rici et al., eq. A8
  !
  ! ************************************************************
CONTAINS
7
  SUBROUTINE force_a8(input,atoms,sphhar,&
8 9 10 11 12 13 14 15 16
       jsp, vr,rho, f_a12,f_a21,f_b4,f_b8,&
       force)
    !
    USE m_intgr, ONLY : intgr3
    USE m_constants,ONLY: pi_const,sfp_const
    USE m_gaunt, ONLY :gaunt1
    USE m_differentiate,ONLY: difcub
    USE m_types
    IMPLICIT NONE
17
    TYPE(t_input),INTENT(IN)       :: input
18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 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 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 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 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261
    TYPE(t_sphhar),INTENT(IN)      :: sphhar
    TYPE(t_atoms),INTENT(IN)       :: atoms
    !     ..
    !     .. Scalar Arguments ..
    INTEGER, INTENT (IN) :: jsp 
    !     ..
    !     .. Array Arguments ..
    COMPLEX, INTENT (IN) :: f_a12(3,atoms%ntypd),f_a21(3,atoms%ntypd)
    COMPLEX, INTENT (IN) :: f_b4(3,atoms%ntypd),f_b8(3,atoms%ntypd)
    REAL,    INTENT (IN) :: vr(atoms%jmtd,0:sphhar%nlhd,atoms%ntypd)
    REAL,    INTENT (IN) :: rho(:,0:,:,:)!(atoms%jmtd,0:sphhar%nlhd,atoms%ntypd,DIMENSION%jspd)
    REAL,    INTENT (INOUT) :: force(:,:,:) !(3,ntypd,jspd)
    !     ..
    !     .. Local Scalars ..
    COMPLEX aaa,bbb,ccc,ddd,eee,fff
    REAL a8_1,a8_2,qval,rr,truerho,truev,xi
    INTEGER i,ir,j,l,l1,l2,lh1,lh2 ,m1,m2,mem1,mem2,n,nd,na,m
    !     ..
    !     .. Local Arrays ..
    COMPLEX forc_a8(3),gv(3),f_sum(3)
    REAL rhoaux(atoms%jmtd),rhodif(atoms%jmtd)
    !     ..
    !     .. Statement Functions ..
    REAL   alpha,beta,delta,epslon,gamma,phi 
    INTEGER krondel
    !     ..
    !     .. Data statements ..
    COMPLEX,PARAMETER:: czero=CMPLX(0.000,0.000)
    COMPLEX,PARAMETER:: ci = CMPLX(0.0,1.0)
    !
    !  inline functions:
    !
    ! Kronecker delta for arguments >=0 AND <0
    !
    !
    krondel(i,j) = MIN(ABS(i)+1,ABS(j)+1)/MAX(ABS(i)+1,ABS(j)+1)* (1+SIGN(1,i)*SIGN(1,j))/2
    alpha(l,m) = (l+1)*0.5e0*SQRT(REAL((l-m)* (l-m-1))/ REAL((2*l-1)* (2*l+1)))
    beta(l,m) = l*0.5e0*SQRT(REAL((l+m+2)* (l+m+1))/ REAL((2*l+1)* (2*l+3)))
    GAMMA(l,m) = (l+1)*0.5e0*SQRT(REAL((l+m)* (l+m-1))/ REAL((2*l-1)* (2*l+1)))
    delta(l,m) = l*0.5e0*SQRT(REAL((l-m+2)* (l-m+1))/ REAL((2*l+1)* (2*l+3)))
    epslon(l,m) = (l+1)*SQRT(REAL((l-m)* (l+m))/ REAL((2*l-1)* (2*l+1)))
    phi(l,m) = l*SQRT(REAL((l-m+1)* (l+m+1))/REAL((2*l+1)* (2*l+3)))
    !     ..
    WRITE  (6,*)
    WRITE (16,*)
    !
    na = 1
    DO  n = 1,atoms%ntype
       IF (atoms%l_geo(n)) THEN
          DO i = 1,3
             forc_a8(i) = czero
          END DO
          !
          !
          nd = atoms%ntypsy(na)
          !
          CALL intgr3(rho(:,0,n,jsp),atoms%rmsh(:,n),atoms%dx(n),atoms%jri(n),qval)
          !
          !     check if l=0 density is correct;
          !     note that in general also all l>0
          !     components of the density have been multiplied by r**2
          !     (see for example subr. checkdop which constructs true MT-density at Rmt)
          !     factor sqrt(4pi) comes from Y_00 * \int d\Omega = 1/sqrt(4pi) * 4pi
          !     write(16,1616) qval*sfp
8000      FORMAT (' FORCE_A8: valence charge=',1p,e16.8)
          !
          !    PART I of FORCE_A8
          !
          DO  lh1 = 0,sphhar%nlh(nd)
             l1 = sphhar%llh(lh1,nd)
             DO  lh2 = 0,sphhar%nlh(nd)
                l2 = sphhar%llh(lh2,nd)
                !
                DO i = 1,3
                   gv(i) = czero
                END DO
                !
                !    sum over all m for particular symm. harmonic
                DO mem1 = 1,sphhar%nmem(lh1,nd)
                   m1 = sphhar%mlh(mem1,lh1,nd)
                   DO mem2 = 1,sphhar%nmem(lh2,nd)
                      m2 = sphhar%mlh(mem2,lh2,nd)
                      gv(1) = gv(1) + SQRT(2.e0*pi_const/3.e0)*&
                           &                       sphhar%clnu(mem1,lh1,nd)*sphhar%clnu(mem2,lh2,nd)*&
                           &                       (gaunt1(1,l1,l2,-1,m1,m2,atoms%lmaxd)-&
                           &                       gaunt1(1,l1,l2,1,m1,m2,atoms%lmaxd))
                      gv(2) = gv(2) - ci*SQRT(2.e0*pi_const/3.e0)*&
                           &                       sphhar%clnu(mem1,lh1,nd)*sphhar%clnu(mem2,lh2,nd)*&
                           &                       (gaunt1(1,l1,l2,-1,m1,m2,atoms%lmaxd)+&
                           &                       gaunt1(1,l1,l2,1,m1,m2,atoms%lmaxd))
                      gv(3) = gv(3) + SQRT(4.e0*pi_const/3.e0)*&
                           &                       sphhar%clnu(mem1,lh1,nd)*sphhar%clnu(mem2,lh2,nd)*&
                           &                       gaunt1(1,l1,l2,0,m1,m2,atoms%lmaxd)
                   END DO
                END DO
                !
                !     note that in general also all l>0
                !     components of the density have been multiplied by r**2
                !     here we need the true radial denisity for performing the derivative
                !     therefore we divide by r**2
                DO ir = 1,atoms%jri(n)
                   rhoaux(ir) = rho(ir,lh2,n,jsp)/ (atoms%rmsh(ir,n)**2)
                END DO
                !
                !
                ! NOTE: we should have here: vr  = vtrue
                !     difcub performs analytic derivative of Lagrangian of 3rd order
                !
                xi = atoms%rmsh(1,n)
                rr = xi*xi
                rhodif(1) = difcub(atoms%rmsh(1,n),rhoaux(1),xi)*&
                     &                     vr(1,lh1,n)*rr
                DO ir = 2,atoms%jri(n) - 2
                   xi = atoms%rmsh(ir,n)
                   rr = xi*xi
                   rhodif(ir) = difcub(atoms%rmsh(ir-1,n),rhoaux(ir-1),xi)*&
                        &                         vr(ir,lh1,n)*rr
                END DO
                !
                xi = atoms%rmsh(atoms%jri(n)-1,n)
                rr = xi*xi
                rhodif(atoms%jri(n)-1) = difcub(atoms%rmsh(atoms%jri(n)-3,n),&
                     &                            rhoaux(atoms%jri(n)-3),xi)*&
                     &                            vr(atoms%jri(n)-1,lh1,n)*rr
                !
                xi = atoms%rmsh(atoms%jri(n),n)
                rr = xi*xi
                rhodif(atoms%jri(n)) = difcub(atoms%rmsh(atoms%jri(n)-3,n),&
                     &                          rhoaux(atoms%jri(n)-3),xi)*&
                     &                          vr(atoms%jri(n),lh1,n)*rr
                ! NOTE: vr(l=0) is EXPLICITELY multiplied by r/sqrt(4pi) to be the TRUE
                ! r*V which is needed for the radial Schr. equ.
                ! here, we need the l=0 component, v_l=0(r) which will be multiplied by
                ! Y_00 in the lm expansion; therefore we MUST recorrect vr(l=0) by
                ! the incerse factor sqrt(4pi)/r. We do the correction for the product
                ! array rhodif
                IF (lh1.EQ.0) THEN
                   DO ir = 1,atoms%jri(n)
                      rhodif(ir) = rhodif(ir)/atoms%rmsh(ir,n)*sfp_const
                   END DO
                END IF
                CALL intgr3(rhodif,atoms%rmsh(1,n),atoms%dx(n),atoms%jri(n),a8_1)
                !
                !
                DO i = 1,3
                   forc_a8(i) = forc_a8(i) + a8_1*gv(i)
                END DO
                !
                !  lh1,lh2 loops end
             ENDDO
          ENDDO
          !     END OF PART I
          !
          !    PART II of FORCE_A8
          !
          DO  lh1 = 0,sphhar%nlh(nd)
             l1 = sphhar%llh(lh1,nd)
             DO  lh2 = 0,sphhar%nlh(nd)
                l2 = sphhar%llh(lh2,nd)
                !
                DO ir = 1,atoms%jri(n)
                   truev = vr(ir,lh1,n)
                   truerho = rho(ir,lh2,n,jsp)/ (atoms%rmsh(ir,n)**2)
                   rhoaux(ir) = truev*truerho*atoms%rmsh(ir,n)
                END DO
                IF (lh1.EQ.0) THEN
                   DO ir = 1,atoms%jri(n)
                      rhoaux(ir) = rhoaux(ir)/atoms%rmsh(ir,n)*sfp_const
                   END DO
                END IF
                CALL intgr3(rhoaux,atoms%rmsh(1,n),atoms%dx(n),atoms%jri(n),a8_2)
                !
                DO i = 1,3
                   gv(i) = (0.0,0.0)
                END DO
                !
                !    sum over all m for particular sym. harmonic
                DO mem1 = 1,sphhar%nmem(lh1,nd)
                   m1 = sphhar%mlh(mem1,lh1,nd)
                   DO mem2 = 1,sphhar%nmem(lh2,nd)
                      m2 = sphhar%mlh(mem2,lh2,nd)
                      !  NOTE: delta(-m,m')(-1)^m was applied because we have integrand
                      !        Y nabla Y instead of Y* nabla Y.
                      !        Then we know that  Y*(lm) = (-1)^m Y(l,-m).
                      !                      and  Y(lm)  = (-1)^m Y*(l,-m).
                      !        Therefore  (-1)^m delta (-m,m') appears
                      !
                      aaa = alpha(l2,m2)*sphhar%clnu(mem1,lh1,nd)*&
                           &                     sphhar%clnu(mem2,lh2,nd)* (-1)**m1*krondel(l1,l2-1)*&
                           &                     krondel(-m1,m2+1)
                      bbb = beta(l2,m2)*sphhar%clnu(mem1,lh1,nd)*&
                           &                     sphhar%clnu(mem2,lh2,nd)* (-1)**m1*krondel(l1,l2+1)*&
                           &                     krondel(-m1,m2+1)
                      ccc = GAMMA(l2,m2)*sphhar%clnu(mem1,lh1,nd)*&
                           &                     sphhar%clnu(mem2,lh2,nd)* (-1)**m1*krondel(l1,l2-1)*&
                           &                     krondel(-m1,m2-1)
                      ddd = delta(l2,m2)*sphhar%clnu(mem1,lh1,nd)*&
                           &                     sphhar%clnu(mem2,lh2,nd)* (-1)**m1*krondel(l1,l2+1)*&
                           &                     krondel(-m1,m2-1)
                      eee = epslon(l2,m2)*sphhar%clnu(mem1,lh1,nd)*&
                           &                     sphhar%clnu(mem2,lh2,nd)* (-1)**m1*krondel(l1,l2-1)*&
                           &                     krondel(-m1,m2)
                      fff = phi(l2,m2)*sphhar%clnu(mem1,lh1,nd)*&
                           &                     sphhar%clnu(mem2,lh2,nd)* (-1)**m1*krondel(l1,l2+1)*&
                           &                     krondel(-m1,m2)
                      !
                      gv(1) = gv(1) + aaa + bbb - ccc - ddd
                      gv(2) = gv(2) - ci* (aaa+bbb+ccc+ddd)
                      gv(3) = gv(3) + eee - fff
                      !
                      !  end of summation m1,m2
                   END DO
                END DO
                !
                DO i = 1,3
                   forc_a8(i) = forc_a8(i) + a8_2*gv(i)
                END DO
                !
                !  lh1,lh2 loops end
             ENDDO
          ENDDO
          !
          !     sum to existing forces
          !
          DO i = 1,3
             force(i,n,jsp) = force(i,n,jsp) + REAL(forc_a8(i))
          END DO
          !
          !     write result
          !
          WRITE (6,FMT=8010) n
          WRITE (16,FMT=8010) n
          WRITE (6,FMT=8020) (forc_a8(i),i=1,3)
          WRITE (16,FMT=8020) (forc_a8(i),i=1,3)
8010      FORMAT (' FORCES: EQUATION A8 FOR ATOM TYPE',i4)
8020      FORMAT (' FX_A8=',2f10.6,' FY_A8=',2f10.6,' FZ_A8=',2f10.6)

       ENDIF
       na = na + atoms%neq(n)
    ENDDO
    !
    !
    !     write now also result of a12 & a21 ( also b4 and b8 )
    !     and b4
262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301
    IF (.NOT.input%l_useapw) THEN

       WRITE  (6,*)
       WRITE (16,*)
       DO n=1,atoms%ntype
          IF (atoms%l_geo(n)) THEN
             WRITE  (6,FMT=8030) n
             WRITE (16,FMT=8030) n
             WRITE  (6,FMT=8040) (f_a12(i,n),i=1,3)
             WRITE (16,FMT=8040) (f_a12(i,n),i=1,3)
          ENDIF
8030      FORMAT (' FORCES: EQUATION A12 FOR ATOM TYPE',i4)
8040      FORMAT (' FX_A12=',2f10.6,' FY_A12=',2f10.6,' FZ_A12=',2f10.6)
       ENDDO
    ELSE
       WRITE  (6,*)
       WRITE (16,*)
       DO n=1,atoms%ntype
          IF (atoms%l_geo(n)) THEN
             WRITE  (6,FMT=8070) n
             WRITE (16,FMT=8070) n
             WRITE  (6,FMT=8080) (f_b4(i,n),i=1,3)
             WRITE (16,FMT=8080) (f_b4(i,n),i=1,3)
          ENDIF
8070      FORMAT (' FORCES: EQUATION B4 FOR ATOM TYPE',i4)
8080      FORMAT (' FX_B4=',2f10.6,' FY_B4=',2f10.6,' FZ_B4=',2f10.6)
       ENDDO
       WRITE  (6,*)
       WRITE (16,*)
       DO n=1,atoms%ntype
          IF (atoms%l_geo(n)) THEN
             WRITE  (6,FMT=8090) n
             WRITE (16,FMT=8090) n
             WRITE  (6,FMT=8100) (f_b8(i,n),i=1,3)
             WRITE (16,FMT=8100) (f_b8(i,n),i=1,3)
          ENDIF
8090      FORMAT (' FORCES: EQUATION B8 FOR ATOM TYPE',i4)
8100      FORMAT (' FX_B8=',2f10.6,' FY_B8=',2f10.6,' FZ_B8=',2f10.6)
       ENDDO
    ENDIF
302 303 304 305 306 307 308 309 310 311 312 313 314 315 316
    WRITE  (6,*)
    WRITE (16,*)
    DO n=1,atoms%ntype
       IF (atoms%l_geo(n)) THEN
          WRITE  (6,FMT=8050) n
          WRITE (16,FMT=8050) n
          WRITE  (6,FMT=8060) (f_a21(i,n),i=1,3)
          WRITE (16,FMT=8060) (f_a21(i,n),i=1,3)
       ENDIF
8050   FORMAT (' FORCES: EQUATION A21 FOR ATOM TYPE',i4)
8060   FORMAT (' FX_A21=',2f10.6,' FY_A21=',2f10.6,' FZ_A21=',2f10.6)
    ENDDO

  END SUBROUTINE force_a8
END MODULE m_forcea8