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

          Line data    Source code
       1             : MODULE m_force_a4
       2             : CONTAINS
       3           2 :   SUBROUTINE force_a4(atoms,sphhar,input,dimension,&
       4           2 :        &                    vr,&
       5           2 :        &                    force)
       6             :     !
       7             :     ! ************************************************************
       8             :     ! rhocore force contribution a la Rici et al.
       9             :     !
      10             :     ! ************************************************************
      11             :     !
      12             :     USE m_intgr, ONLY : intgr0,intgr3
      13             :     USE m_constants, ONLY : sfp_const
      14             :     USE m_differentiate,ONLY: difcub
      15             :     USE m_types
      16             :     USE m_cdn_io
      17             :     IMPLICIT NONE
      18             :     TYPE(t_input),INTENT(IN)     :: input
      19             :     TYPE(t_sphhar),INTENT(IN)    :: sphhar
      20             :     TYPE(t_atoms),INTENT(IN)     :: atoms
      21             :     TYPE(t_dimension),INTENT(IN) :: dimension
      22             :     !     ..
      23             :     !     .. Array Arguments ..
      24             :     REAL,    INTENT (IN) :: vr(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,input%jspins)
      25             :     REAL,    INTENT (INOUT) :: force(3,atoms%ntype,input%jspins)
      26             :     !     ..
      27             :     !     .. Local Scalars ..
      28             :     REAL a4_1,a4_2,qcore,s13,s23,w,xi
      29             :     INTEGER i,ir,jsp,lh,lindex,mem,mindex,n,nd,na
      30             :     !     ..
      31             :     !     .. Local Arrays ..
      32             :     COMPLEX forc_a4(3),gv(3),ycomp1(3,-1:1)
      33           4 :     REAL rhoaux(atoms%jmtd),rhoc(atoms%jmtd,atoms%ntype,input%jspins)
      34           8 :     REAL tec(atoms%ntype,input%jspins),qintc(atoms%ntype,input%jspins)
      35             :     !     ..
      36             :     !     .. Data statements ..
      37             :     COMPLEX,PARAMETER:: czero=(0.000,0.000)
      38             :     !     ..
      39             :     !
      40             :     !     set ylm related components
      41             :     !
      42           2 :     s13 = SQRT(1.0/3.0)
      43           2 :     s23 = SQRT(2.0/3.0)
      44           2 :     ycomp1(1,0) = czero
      45           2 :     ycomp1(2,0) = czero
      46           2 :     ycomp1(3,0) = CMPLX(2.0*s13,0.0)
      47           2 :     ycomp1(1,-1) = CMPLX(s23,0.0)
      48           2 :     ycomp1(2,-1) = CMPLX(0.0,-s23)
      49           2 :     ycomp1(3,-1) = czero
      50           2 :     ycomp1(1,1) = CMPLX(-s23,0.0)
      51           2 :     ycomp1(2,1) = CMPLX(0.0,-s23)
      52           2 :     ycomp1(3,1) = czero
      53             :     !     --->    read in core density
      54           2 :     CALL readCoreDensity(input,atoms,dimension,rhoc,tec,qintc)
      55             : 
      56           4 :     DO jsp = 1,input%jspins
      57           2 :        na = 1
      58           8 :        DO n = 1,atoms%ntype
      59           4 :           IF (atoms%l_geo(n)) THEN
      60          28 :              DO i = 1,3
      61          16 :                 forc_a4(i) = czero
      62             :              END DO
      63             :              !
      64           4 :              CALL intgr0(rhoc(:,n,jsp),atoms%rmsh(1,n),atoms%dx(n),atoms%jri(n),qcore)
      65             : 8000         FORMAT (' FORCE_A4: core charge=',1p,e16.8)
      66             :              !
      67             :              !
      68           4 :              nd = atoms%ntypsy(na)
      69             :              !
      70             :              !
      71           8 :              lh_loop: DO  lh = 1,sphhar%nlh(nd)
      72           8 :                 lindex = sphhar%llh(lh,nd)
      73           8 :                 IF (lindex.GT.1) EXIT lh_loop
      74          28 :                 DO i = 1,3
      75          16 :                    gv(i) = czero
      76             :                 END DO
      77             :                 !
      78             :                 !    sum over all m for particular symm. harmonic
      79          16 :                 DO mem = 1,sphhar%nmem(lh,nd)
      80          12 :                    mindex = sphhar%mlh(mem,lh,nd)
      81          52 :                    DO i = 1,3
      82          48 :                       gv(i) = gv(i) + sphhar%clnu(mem,lh,nd)*ycomp1(i,mindex)
      83             :                    END DO
      84             :                 END DO
      85             :                 !
      86             :                 !
      87             :                 !     construct integrand rhocore*d/dr(v)*r**2
      88             :                 !     note: rhocore is already multiplied by r**2 and srt(4.*pi)
      89             :                 !     difcub performs analytic derivative of Lagrangian of 3rd order
      90             :                 !
      91           4 :                 xi = atoms%rmsh(1,n)
      92           4 :                 rhoaux(1) = difcub(atoms%rmsh(1,n),vr(1,lh,n,jsp),xi)*rhoc(1,n,jsp)
      93        2916 :                 DO ir = 2,atoms%jri(n) - 2
      94        2912 :                    xi = atoms%rmsh(ir,n)
      95             :                    rhoaux(ir) = difcub(atoms%rmsh(ir-1,n),&
      96        2916 :                         &                                vr(ir-1,lh,n,jsp),xi) * rhoc(ir,n,jsp)
      97             :                 END DO
      98             :                 !
      99           4 :                 ir = atoms%jri(n) - 1
     100           4 :                 xi = atoms%rmsh(ir,n)
     101             :                 rhoaux(ir) = difcub(atoms%rmsh(atoms%jri(n)-3,n),&
     102           4 :                      &                              vr(atoms%jri(n)-3,lh,n,jsp),xi)*rhoc(ir,n,jsp)
     103             :                 !
     104           4 :                 ir = atoms%jri(n)
     105           4 :                 xi = atoms%rmsh(ir,n)
     106             :                 rhoaux(ir) = difcub(atoms%rmsh(atoms%jri(n)-3,n),&
     107           4 :                      &                              vr(atoms%jri(n)-3,lh,n,jsp),xi)*rhoc(ir,n,jsp)
     108           4 :                 CALL intgr3(rhoaux,atoms%rmsh(1,n),atoms%dx(n),atoms%jri(n),w)
     109           4 :                 a4_1 = 0.5*w/sfp_const
     110             :                 !
     111             :                 !     construct integrand rhocore*v*r
     112             :                 !     note: rhocore is already multiplied by r**2 and srt(4.*pi)
     113             :                 !
     114        2928 :                 DO ir = 1,atoms%jri(n)
     115        2928 :                    rhoaux(ir) = rhoc(ir,n,jsp)/atoms%rmsh(ir,n)*vr(ir,lh,n,jsp)
     116             :                 END DO
     117             :                 !
     118           4 :                 CALL intgr3(rhoaux,atoms%rmsh(1,n),atoms%dx(n),atoms%jri(n),w)
     119           4 :                 a4_2 = w/sfp_const
     120             :                 !
     121          20 :                 DO i = 1,3
     122          16 :                    forc_a4(i) = forc_a4(i) - (a4_1+a4_2)*gv(i)
     123             :                 END DO
     124             :                 !
     125             :                 !  lh loop ends
     126             :              ENDDO lh_loop
     127             :              !
     128             :              !     sum to existing forces
     129             :              !
     130          28 :              DO i = 1,3
     131          16 :                 force(i,n,jsp) = force(i,n,jsp) + REAL(forc_a4(i))
     132             :              END DO
     133             :              !
     134             :              !     write result
     135             :              !
     136           4 :              WRITE (6,FMT=8010) n
     137           4 :              WRITE (6,FMT=8020) (forc_a4(i),i=1,3)
     138             : 8010         FORMAT (' FORCES: EQUATION A4 FOR ATOM TYPE',i4)
     139             : 8020         FORMAT (' FX_A4=',2f10.6,' FY_A4=',2f10.6,' FZ_A4=',2f10.6)
     140             :              ! type loop ends
     141             :           ENDIF
     142           6 :           na = na + atoms%neq(n)
     143             :        ENDDO
     144             :        ! spin loop ends
     145             :     ENDDO
     146           2 :   END SUBROUTINE force_a4
     147             : END MODULE m_force_a4

Generated by: LCOV version 1.13