LCOV - code coverage report
Current view: top level - force - force_a8.F90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 92 102 90.2 %
Date: 2019-09-08 04:53:50 Functions: 1 1 100.0 %

          Line data    Source code
       1             : MODULE m_forcea8
       2             :   ! ************************************************************
       3             :   ! Pulay 1st term  force contribution a la Rici et al., eq. A8
       4             :   !
       5             :   ! ************************************************************
       6             : CONTAINS
       7           2 :   SUBROUTINE force_a8(input,atoms,sphhar,jsp,vr,rho,force,results)
       8             :     !
       9             :     USE m_intgr, ONLY : intgr3
      10             :     USE m_constants, ONLY: pi_const, sfp_const, ImagUnit
      11             :     USE m_gaunt, ONLY :gaunt1
      12             :     USE m_differentiate,ONLY: difcub
      13             :     USE m_types
      14             :     USE m_juDFT
      15             :     IMPLICIT NONE
      16             :     TYPE(t_input),INTENT(IN)       :: input
      17             :     TYPE(t_sphhar),INTENT(IN)      :: sphhar
      18             :     TYPE(t_atoms),INTENT(IN)       :: atoms
      19             :     TYPE(t_force),INTENT(IN)       :: force
      20             :     TYPE(t_results),INTENT(INOUT)  :: results
      21             :     !     ..
      22             :     !     .. Scalar Arguments ..
      23             :     INTEGER, INTENT (IN) :: jsp 
      24             :     !     ..
      25             :     !     .. Array Arguments ..
      26             :     REAL,    INTENT (IN) :: vr(atoms%jmtd,0:sphhar%nlhd,atoms%ntype)
      27             :     REAL,    INTENT (IN) :: rho(:,0:,:,:)!(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,input%jspins)
      28             :     !     ..
      29             :     !     .. Local Scalars ..
      30             :     COMPLEX aaa,bbb,ccc,ddd,eee,fff
      31             :     REAL a8_1,a8_2,qval,rr,truerho,truev,xi
      32             :     INTEGER i,ir,j,l,l1,l2,lh1,lh2 ,m1,m2,mem1,mem2,n,nd,na,m
      33             :     !     ..
      34             :     !     .. Local Arrays ..
      35             :     COMPLEX forc_a8(3),gv(3),f_sum(3)
      36           4 :     REAL rhoaux(atoms%jmtd),rhodif(atoms%jmtd)
      37             :     !     ..
      38             :     !     .. Statement Functions ..
      39             :     REAL   alpha,beta,delta,epslon,gamma,phi 
      40             :     INTEGER krondel
      41             :     !     ..
      42             :     !     .. Data statements ..
      43             :     COMPLEX,PARAMETER:: czero=CMPLX(0.000,0.000)
      44             : 
      45             :     ! Kronecker delta for arguments >=0 AND <0
      46             :     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
      47             :     alpha(l,m) = (l+1)*0.5e0*SQRT(REAL((l-m)* (l-m-1))/ REAL((2*l-1)* (2*l+1)))
      48             :     beta(l,m) = l*0.5e0*SQRT(REAL((l+m+2)* (l+m+1))/ REAL((2*l+1)* (2*l+3)))
      49             :     GAMMA(l,m) = (l+1)*0.5e0*SQRT(REAL((l+m)* (l+m-1))/ REAL((2*l-1)* (2*l+1)))
      50             :     delta(l,m) = l*0.5e0*SQRT(REAL((l-m+2)* (l-m+1))/ REAL((2*l+1)* (2*l+3)))
      51             :     epslon(l,m) = (l+1)*SQRT(REAL((l-m)* (l+m))/ REAL((2*l-1)* (2*l+1)))
      52             :     phi(l,m) = l*SQRT(REAL((l-m+1)* (l+m+1))/REAL((2*l+1)* (2*l+3)))
      53             : 
      54           2 :     CALL timestart("force_a8")
      55             : 
      56           2 :     WRITE  (6,*)
      57             :  
      58           2 :     na = 1
      59           6 :     DO  n = 1,atoms%ntype
      60           4 :        IF (atoms%l_geo(n)) THEN
      61          28 :           DO i = 1,3
      62          16 :              forc_a8(i) = czero
      63             :           END DO
      64             :           !
      65             :           !
      66           4 :           nd = atoms%ntypsy(na)
      67             :           !
      68           4 :           CALL intgr3(rho(:,0,n,jsp),atoms%rmsh(:,n),atoms%dx(n),atoms%jri(n),qval)
      69             : 
      70             :           !     check if l=0 density is correct;
      71             :           !     note that in general also all l>0
      72             :           !     components of the density have been multiplied by r**2
      73             :           !     (see for example subr. checkdop which constructs true MT-density at Rmt)
      74             :           !     factor sqrt(4pi) comes from Y_00 * \int d\Omega = 1/sqrt(4pi) * 4pi
      75             : 8000      FORMAT (' FORCE_A8: valence charge=',1p,e16.8)
      76             : 
      77             :           !    PART I of FORCE_A8
      78          76 :           DO  lh1 = 0,sphhar%nlh(nd)
      79          72 :              l1 = sphhar%llh(lh1,nd)
      80        1372 :              DO  lh2 = 0,sphhar%nlh(nd)
      81        1296 :                 l2 = sphhar%llh(lh2,nd)
      82             :                 !
      83        5184 :                 DO i = 1,3
      84        5184 :                    gv(i) = czero
      85             :                 END DO
      86             :                 !
      87             :                 !    sum over all m for particular symm. harmonic
      88        9216 :                 DO mem1 = 1,sphhar%nmem(lh1,nd)
      89        7920 :                    m1 = sphhar%mlh(mem1,lh1,nd)
      90       57616 :                    DO mem2 = 1,sphhar%nmem(lh2,nd)
      91       48400 :                       m2 = sphhar%mlh(mem2,lh2,nd)
      92             :                       gv(1) = gv(1) + SQRT(2.e0*pi_const/3.e0)*&
      93             :                            &                       sphhar%clnu(mem1,lh1,nd)*sphhar%clnu(mem2,lh2,nd)*&
      94             :                            &                       (gaunt1(1,l1,l2,-1,m1,m2,atoms%lmaxd)-&
      95       48400 :                            &                       gaunt1(1,l1,l2,1,m1,m2,atoms%lmaxd))
      96             :                       gv(2) = gv(2) - ImagUnit*SQRT(2.e0*pi_const/3.e0)*&
      97             :                            &                       sphhar%clnu(mem1,lh1,nd)*sphhar%clnu(mem2,lh2,nd)*&
      98             :                            &                       (gaunt1(1,l1,l2,-1,m1,m2,atoms%lmaxd)+&
      99       48400 :                            &                       gaunt1(1,l1,l2,1,m1,m2,atoms%lmaxd))
     100             :                       gv(3) = gv(3) + SQRT(4.e0*pi_const/3.e0)*&
     101             :                            &                       sphhar%clnu(mem1,lh1,nd)*sphhar%clnu(mem2,lh2,nd)*&
     102       56320 :                            &                       gaunt1(1,l1,l2,0,m1,m2,atoms%lmaxd)
     103             :                    END DO
     104             :                 END DO
     105             :                 !
     106             :                 !     note that in general also all l>0
     107             :                 !     components of the density have been multiplied by r**2
     108             :                 !     here we need the true radial denisity for performing the derivative
     109             :                 !     therefore we divide by r**2
     110      948672 :                 DO ir = 1,atoms%jri(n)
     111      948672 :                    rhoaux(ir) = rho(ir,lh2,n,jsp)/ (atoms%rmsh(ir,n)**2)
     112             :                 END DO
     113             :                 !
     114             :                 !
     115             :                 ! NOTE: we should have here: vr  = vtrue
     116             :                 !     difcub performs analytic derivative of Lagrangian of 3rd order
     117             :                 !
     118        1296 :                 xi = atoms%rmsh(1,n)
     119        1296 :                 rr = xi*xi
     120             :                 rhodif(1) = difcub(atoms%rmsh(1,n),rhoaux(1),xi)*&
     121        1296 :                      &                     vr(1,lh1,n)*rr
     122      944784 :                 DO ir = 2,atoms%jri(n) - 2
     123      943488 :                    xi = atoms%rmsh(ir,n)
     124      943488 :                    rr = xi*xi
     125             :                    rhodif(ir) = difcub(atoms%rmsh(ir-1,n),rhoaux(ir-1),xi)*&
     126      944784 :                         &                         vr(ir,lh1,n)*rr
     127             :                 END DO
     128             :                 !
     129        1296 :                 xi = atoms%rmsh(atoms%jri(n)-1,n)
     130        1296 :                 rr = xi*xi
     131             :                 rhodif(atoms%jri(n)-1) = difcub(atoms%rmsh(atoms%jri(n)-3,n),&
     132             :                      &                            rhoaux(atoms%jri(n)-3),xi)*&
     133        1296 :                      &                            vr(atoms%jri(n)-1,lh1,n)*rr
     134             :                 !
     135        1296 :                 xi = atoms%rmsh(atoms%jri(n),n)
     136        1296 :                 rr = xi*xi
     137             :                 rhodif(atoms%jri(n)) = difcub(atoms%rmsh(atoms%jri(n)-3,n),&
     138             :                      &                          rhoaux(atoms%jri(n)-3),xi)*&
     139        1296 :                      &                          vr(atoms%jri(n),lh1,n)*rr
     140             :                 ! NOTE: vr(l=0) is EXPLICITELY multiplied by r/sqrt(4pi) to be the TRUE
     141             :                 ! r*V which is needed for the radial Schr. equ.
     142             :                 ! here, we need the l=0 component, v_l=0(r) which will be multiplied by
     143             :                 ! Y_00 in the lm expansion; therefore we MUST recorrect vr(l=0) by
     144             :                 ! the incerse factor sqrt(4pi)/r. We do the correction for the product
     145             :                 ! array rhodif
     146        1296 :                 IF (lh1.EQ.0) THEN
     147      105336 :                    DO ir = 1,atoms%jri(n)
     148       52704 :                       rhodif(ir) = rhodif(ir)/atoms%rmsh(ir,n)*sfp_const
     149             :                    END DO
     150             :                 END IF
     151        1296 :                 CALL intgr3(rhodif,atoms%rmsh(1,n),atoms%dx(n),atoms%jri(n),a8_1)
     152             :                 !
     153             :                 !
     154        5256 :                 DO i = 1,3
     155        5184 :                    forc_a8(i) = forc_a8(i) + a8_1*gv(i)
     156             :                 END DO
     157             :                 !
     158             :                 !  lh1,lh2 loops end
     159             :              ENDDO
     160             :           ENDDO
     161             :           !     END OF PART I
     162             :           !
     163             :           !    PART II of FORCE_A8
     164             :           !
     165          76 :           DO  lh1 = 0,sphhar%nlh(nd)
     166          72 :              l1 = sphhar%llh(lh1,nd)
     167        1372 :              DO  lh2 = 0,sphhar%nlh(nd)
     168        1296 :                 l2 = sphhar%llh(lh2,nd)
     169             :                 !
     170      948672 :                 DO ir = 1,atoms%jri(n)
     171      947376 :                    truev = vr(ir,lh1,n)
     172      947376 :                    truerho = rho(ir,lh2,n,jsp)/ (atoms%rmsh(ir,n)**2)
     173      948672 :                    rhoaux(ir) = truev*truerho*atoms%rmsh(ir,n)
     174             :                 END DO
     175        1296 :                 IF (lh1.EQ.0) THEN
     176      105336 :                    DO ir = 1,atoms%jri(n)
     177       52704 :                       rhoaux(ir) = rhoaux(ir)/atoms%rmsh(ir,n)*sfp_const
     178             :                    END DO
     179             :                 END IF
     180        1296 :                 CALL intgr3(rhoaux,atoms%rmsh(1,n),atoms%dx(n),atoms%jri(n),a8_2)
     181             :                 !
     182        5184 :                 DO i = 1,3
     183        5184 :                    gv(i) = (0.0,0.0)
     184             :                 END DO
     185             :                 !
     186             :                 !    sum over all m for particular sym. harmonic
     187        9216 :                 DO mem1 = 1,sphhar%nmem(lh1,nd)
     188        7920 :                    m1 = sphhar%mlh(mem1,lh1,nd)
     189       57616 :                    DO mem2 = 1,sphhar%nmem(lh2,nd)
     190       48400 :                       m2 = sphhar%mlh(mem2,lh2,nd)
     191             :                       !  NOTE: delta(-m,m')(-1)^m was applied because we have integrand
     192             :                       !        Y nabla Y instead of Y* nabla Y.
     193             :                       !        Then we know that  Y*(lm) = (-1)^m Y(l,-m).
     194             :                       !                      and  Y(lm)  = (-1)^m Y*(l,-m).
     195             :                       !        Therefore  (-1)^m delta (-m,m') appears
     196             :                       !
     197             :                       aaa = alpha(l2,m2)*sphhar%clnu(mem1,lh1,nd)*&
     198             :                            &                     sphhar%clnu(mem2,lh2,nd)* (-1)**m1*krondel(l1,l2-1)*&
     199       48400 :                            &                     krondel(-m1,m2+1)
     200             :                       bbb = beta(l2,m2)*sphhar%clnu(mem1,lh1,nd)*&
     201             :                            &                     sphhar%clnu(mem2,lh2,nd)* (-1)**m1*krondel(l1,l2+1)*&
     202       48400 :                            &                     krondel(-m1,m2+1)
     203             :                       ccc = GAMMA(l2,m2)*sphhar%clnu(mem1,lh1,nd)*&
     204             :                            &                     sphhar%clnu(mem2,lh2,nd)* (-1)**m1*krondel(l1,l2-1)*&
     205       48400 :                            &                     krondel(-m1,m2-1)
     206             :                       ddd = delta(l2,m2)*sphhar%clnu(mem1,lh1,nd)*&
     207             :                            &                     sphhar%clnu(mem2,lh2,nd)* (-1)**m1*krondel(l1,l2+1)*&
     208       48400 :                            &                     krondel(-m1,m2-1)
     209             :                       eee = epslon(l2,m2)*sphhar%clnu(mem1,lh1,nd)*&
     210             :                            &                     sphhar%clnu(mem2,lh2,nd)* (-1)**m1*krondel(l1,l2-1)*&
     211       48400 :                            &                     krondel(-m1,m2)
     212             :                       fff = phi(l2,m2)*sphhar%clnu(mem1,lh1,nd)*&
     213             :                            &                     sphhar%clnu(mem2,lh2,nd)* (-1)**m1*krondel(l1,l2+1)*&
     214       48400 :                            &                     krondel(-m1,m2)
     215             :                       !
     216       48400 :                       gv(1) = gv(1) + aaa + bbb - ccc - ddd
     217       48400 :                       gv(2) = gv(2) - ImagUnit* (aaa+bbb+ccc+ddd)
     218       56320 :                       gv(3) = gv(3) + eee - fff
     219             :                       !
     220             :                       !  end of summation m1,m2
     221             :                    END DO
     222             :                 END DO
     223             :                 !
     224        9144 :                 DO i = 1,3
     225        5184 :                    forc_a8(i) = forc_a8(i) + a8_2*gv(i)
     226             :                 END DO
     227             :                 !
     228             :                 !  lh1,lh2 loops end
     229             :              ENDDO
     230             :           ENDDO
     231             :           !
     232             :           !     sum to existing forces
     233             :           !
     234          28 :           DO i = 1,3
     235          16 :              results%force(i,n,jsp) = results%force(i,n,jsp) + REAL(forc_a8(i))
     236             :           END DO
     237             :           !
     238             :           !     write result
     239             :           !
     240           4 :           WRITE (6,FMT=8010) n
     241           4 :           WRITE (6,FMT=8020) (forc_a8(i),i=1,3)
     242             : 8010      FORMAT (' FORCES: EQUATION A8 FOR ATOM TYPE',i4)
     243             : 8020      FORMAT (' FX_A8=',2f10.6,' FY_A8=',2f10.6,' FZ_A8=',2f10.6)
     244             : 
     245             :        ENDIF
     246           6 :        na = na + atoms%neq(n)
     247             :     ENDDO
     248             :     !
     249             :     !
     250             :     !     write now also result of a12 & a21 ( also b4 and b8 )
     251             :     !     and b4
     252           2 :     IF (.NOT.input%l_useapw) THEN
     253             : 
     254           2 :        WRITE  (6,*)
     255           6 :        DO n=1,atoms%ntype
     256           6 :           IF (atoms%l_geo(n)) THEN
     257           4 :              WRITE  (6,FMT=8030) n
     258           4 :              WRITE  (6,FMT=8040) (force%f_a12(i,n),i=1,3)
     259             :           ENDIF
     260             : 8030      FORMAT (' FORCES: EQUATION A12 FOR ATOM TYPE',i4)
     261             : 8040      FORMAT (' FX_A12=',2f10.6,' FY_A12=',2f10.6,' FZ_A12=',2f10.6)
     262             :        ENDDO
     263             :     ELSE
     264           0 :        WRITE  (6,*)
     265           0 :        DO n=1,atoms%ntype
     266           0 :           IF (atoms%l_geo(n)) THEN
     267           0 :              WRITE  (6,FMT=8070) n
     268           0 :              WRITE  (6,FMT=8080) (force%f_b4(i,n),i=1,3)
     269             :           ENDIF
     270             : 8070      FORMAT (' FORCES: EQUATION B4 FOR ATOM TYPE',i4)
     271             : 8080      FORMAT (' FX_B4=',2f10.6,' FY_B4=',2f10.6,' FZ_B4=',2f10.6)
     272             :        ENDDO
     273           0 :        WRITE  (6,*)
     274           0 :        DO n=1,atoms%ntype
     275           0 :           IF (atoms%l_geo(n)) THEN
     276           0 :              WRITE  (6,FMT=8090) n
     277           0 :              WRITE  (6,FMT=8100) (force%f_b8(i,n),i=1,3)
     278             :           ENDIF
     279             : 8090      FORMAT (' FORCES: EQUATION B8 FOR ATOM TYPE',i4)
     280             : 8100      FORMAT (' FX_B8=',2f10.6,' FY_B8=',2f10.6,' FZ_B8=',2f10.6)
     281             :        ENDDO
     282             :     ENDIF
     283           2 :     WRITE  (6,*)
     284           6 :     DO n=1,atoms%ntype
     285           6 :        IF (atoms%l_geo(n)) THEN
     286           4 :           WRITE  (6,FMT=8050) n
     287           4 :           WRITE  (6,FMT=8060) (force%f_a21(i,n),i=1,3)
     288             :        ENDIF
     289             : 8050   FORMAT (' FORCES: EQUATION A21 FOR ATOM TYPE',i4)
     290             : 8060   FORMAT (' FX_A21=',2f10.6,' FY_A21=',2f10.6,' FZ_A21=',2f10.6)
     291             :     ENDDO
     292             : 
     293           2 :     CALL timestop("force_a8")
     294             : 
     295           2 :   END SUBROUTINE force_a8
     296             : END MODULE m_forcea8

Generated by: LCOV version 1.13