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

          Line data    Source code
       1             : MODULE m_force_a3
       2             : CONTAINS
       3         162 :   SUBROUTINE force_a3(atoms,sphhar,&
       4             :        &                    input,&
       5         162 :        &                    rho,vr,&
       6         162 :        &                    force)
       7             :     ! ************************************************************
       8             :     ! Hellman-Feynman force contribution a la Rici et al.
       9             :     ! ************************************************************
      10             :     !
      11             :     USE m_intgr, ONLY : intgr3
      12             :     USE m_constants, ONLY : pi_const
      13             :     USE m_types
      14             :     IMPLICIT NONE
      15             :     TYPE(t_input),INTENT(IN)   :: input
      16             :     TYPE(t_sphhar),INTENT(IN)   :: sphhar
      17             :     TYPE(t_atoms),INTENT(IN)   :: atoms
      18             :     !     .. Array Arguments ..
      19             :     REAL,    INTENT (IN) ::  vr(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,input%jspins)
      20             :     REAL,    INTENT (IN) :: rho(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,input%jspins)
      21             :     REAL,    INTENT (INOUT) :: force(3,atoms%ntype,input%jspins)
      22             :     !     ..
      23             :     !     .. Local Scalars ..
      24             :     COMPLEX,PARAMETER:: czero=(0.0,0.0)
      25             :     REAL a3_1,a3_2,s34,s38,w
      26             :     INTEGER i,ir,jsp,lh,lindex,mem,mindex,n,nd,na
      27             :     !     ..
      28             :     !     .. Local Arrays ..
      29             :     COMPLEX forc_a3(3),grd1(3,-1:1),gv(3)
      30         324 :     REAL rhoaux(atoms%jmtd)
      31             :     !     ..
      32             :     !     set components of gradient in terms of Ylms
      33             :     !
      34         162 :     s34 = SQRT(3.0/ (4.0*pi_const))
      35         162 :     s38 = SQRT(3.0/ (8.0*pi_const))
      36         162 :     grd1(1,0) = czero
      37         162 :     grd1(2,0) = czero
      38         162 :     grd1(3,0) = CMPLX(s34,0.0)
      39         162 :     grd1(1,-1) = CMPLX(s38,0.0)
      40         162 :     grd1(2,-1) = CMPLX(0.0,-s38)
      41         162 :     grd1(3,-1) = czero
      42         162 :     grd1(1,1) = CMPLX(-s38,0.0)
      43         162 :     grd1(2,1) = CMPLX(0.0,-s38)
      44         162 :     grd1(3,1) = czero
      45             :     !
      46         162 :     WRITE  (6,*)
      47         457 :     spin: DO jsp = 1,input%jspins
      48         295 :        na = 1
      49        1020 :        DO n = 1,atoms%ntype
      50         563 :           IF (atoms%l_geo(n)) THEN
      51             : 
      52        3941 :              DO i = 1,3
      53        2252 :                 forc_a3(i) = czero
      54             :              END DO
      55             :              !
      56         563 :              nd = atoms%ntypsy(na)
      57             :              !
      58         762 :              lat_har: DO lh = 1,sphhar%nlh(nd)
      59         762 :                 lindex = sphhar%llh(lh,nd)
      60         762 :                 IF (lindex.GT.1) EXIT
      61        1393 :                 DO i = 1,3
      62         796 :                    gv(i) = czero
      63             :                 END DO
      64             :                 !    sum over all m for particular symm. harmonic
      65         542 :                 DO mem = 1,sphhar%nmem(lh,nd)
      66         343 :                    mindex = sphhar%mlh(mem,lh,nd)
      67        1571 :                    DO i = 1,3
      68        1372 :                       gv(i) = gv(i) + sphhar%clnu(mem,lh,nd)*grd1(i,mindex)
      69             :                    END DO
      70             :                 END DO
      71      120874 :                 DO ir = 1,atoms%jri(n)
      72             :                    rhoaux(ir) = rho(ir,lh,n,jsp)/ (atoms%rmsh(ir,n)**2)*&
      73      120874 :                         &                         (1.0- (atoms%rmsh(ir,n)/atoms%rmt(n))**3)
      74             :                 END DO
      75         199 :                 CALL intgr3(rhoaux,atoms%rmsh(1,n),atoms%dx(n),atoms%jri(n),w)
      76         199 :                 a3_1 = 4.0*pi_const/3.0*w
      77             :                 !+Gu
      78         199 :                 a3_2 = vr(atoms%jri(n),lh,n,jsp)/(input%jspins*atoms%rmt(n))
      79             :                 !-Gu
      80        1359 :                 DO i = 1,3
      81         796 :                    forc_a3(i) = forc_a3(i) + (a3_1+a3_2)*gv(i)*atoms%zatom(n)
      82             :                 END DO
      83             :              ENDDO lat_har
      84             : 
      85        3941 :              DO i = 1,3
      86        2252 :                 force(i,n,jsp) = force(i,n,jsp) + REAL(forc_a3(i))
      87             :              END DO
      88             :              !
      89             :              !     write result
      90         563 :              WRITE (6,FMT=8010) n
      91         563 :              WRITE (6,FMT=8020) (forc_a3(i),i=1,3)
      92             : 8010         FORMAT (' FORCES: EQUATION A3 FOR ATOM TYPE',i4)
      93             : 8020         FORMAT (' FX_A3=',2f10.6,' FY_A3=',2f10.6,' FZ_A3=',2f10.6)
      94             : 
      95             :           ENDIF               ! atoms%l_geo(n) = .true.
      96         858 :           na = na + atoms%neq(n)
      97             :        ENDDO
      98             : 
      99             :     ENDDO spin
     100             : 
     101         162 :   END SUBROUTINE force_a3
     102             : END MODULE m_force_a3

Generated by: LCOV version 1.13