Line data    Source code

       1             : MODULE m_force_a3
       2             : CONTAINS
       3         334 :    SUBROUTINE force_a3(atoms,sym,sphhar,input,rho,vr,force)
       4             :       !--------------------------------------------------------------------------
       5             :       ! Hellman-Feynman force contribution à la Rici et al.
       6             :       ! 
       7             :       ! Equation A3, Phys. Rev. B 43, 6411
       8             :       !--------------------------------------------------------------------------
       9             :       USE m_intgr, ONLY : intgr3
      10             :       USE m_constants
      11             :       USE m_types
      12             : 
      13             :       IMPLICIT NONE
      14             : 
      15             :       TYPE(t_atoms),  INTENT(IN) :: atoms
      16             :       TYPE(t_sym),    INTENT(IN) :: sym
      17             :       TYPE(t_sphhar), INTENT(IN) :: sphhar
      18             :       TYPE(t_input),  INTENT(IN) :: input
      19             : 
      20             :       REAL, INTENT (IN)    :: rho(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,input%jspins)
      21             :       REAL, INTENT (IN)    ::  vr(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,input%jspins)
      22             :       REAL, INTENT (INOUT) :: force(3,atoms%ntype,input%jspins)
      23             : 
      24             :       ! Local scalars
      25             :       COMPLEX, PARAMETER :: czero=(0.0,0.0)
      26             :       REAL    a3_1, a3_2, s34, s38, w
      27             :       INTEGER i, ir, jsp, lh, lindex, mem, mindex, n, nd, na
      28             : 
      29             :       ! Local arrays
      30             :       COMPLEX forc_a3(3), grd1(3,-1:1), gv(3)
      31         334 :       REAL    rhoaux(atoms%jmtd)
      32             : 
      33             :       ! Set components of gradient in terms of Ylms.
      34         334 :       s34 = SQRT(3.0/(4.0*pi_const))
      35         334 :       s38 = SQRT(3.0/(8.0*pi_const))
      36         334 :       grd1(1,0) = czero
      37         334 :       grd1(2,0) = czero
      38         334 :       grd1(3,0) = CMPLX(s34,0.0)
      39         334 :       grd1(1,-1) = CMPLX(s38,0.0)
      40         334 :       grd1(2,-1) = CMPLX(0.0,-s38)
      41         334 :       grd1(3,-1) = czero
      42         334 :       grd1(1,1) = CMPLX(-s38,0.0)
      43         334 :       grd1(2,1) = CMPLX(0.0,-s38)
      44         334 :       grd1(3,1) = czero
      45             : 
      46         334 :       CALL timestart("force_a3")
      47             : 
      48         334 :       WRITE  (oUnit,*)
      49             : 
      50         863 :       DO jsp = 1, input%jspins
      51        1771 :          DO n = 1, atoms%ntype
      52         908 :             na = atoms%firstAtom(n)
      53        1437 :             IF (atoms%l_geo(n)) THEN
      54         904 :                nd = sym%ntypsy(na)
      55             : 
      56        3616 :                DO i = 1, 3
      57        3616 :                   forc_a3(i) = czero
      58             :                END DO
      59             : 
      60        1600 :                DO lh = 1, sphhar%nlh(nd)
      61        1600 :                   lindex = sphhar%llh(lh,nd)
      62             : 
      63        1600 :                   IF (lindex.GT.1) EXIT
      64             : 
      65        2784 :                   DO i = 1, 3
      66        2784 :                      gv(i) = czero
      67             :                   END DO
      68             : 
      69             :                   ! Sum over all m for a particular lattice harmonic.
      70        1880 :                   DO mem = 1, sphhar%nmem(lh,nd)
      71        1184 :                      mindex = sphhar%mlh(mem,lh,nd)
      72        5432 :                      DO i = 1, 3
      73        4736 :                         gv(i) = gv(i) + sphhar%clnu(mem,lh,nd)*grd1(i,mindex)
      74             :                      END DO
      75             :                   END DO
      76             : 
      77      474908 :                   DO ir = 1, atoms%jri(n)
      78      474908 :                      rhoaux(ir) = rho(ir,lh,n,jsp)/ (atoms%rmsh(ir,n)**2)*(1.0- (atoms%rmsh(ir,n)/atoms%rmt(n))**3)
      79             :                   END DO
      80             : 
      81         696 :                   CALL intgr3(rhoaux,atoms%rmsh(1,n),atoms%dx(n),atoms%jri(n),w)
      82             : 
      83         696 :                   a3_1 = 4.0*pi_const/3.0*w
      84             : 
      85             :                   !+Gu
      86         696 :                   a3_2 = vr(atoms%jri(n),lh,n,jsp) / (input%jspins*atoms%rmt(n))
      87             :                   !-Gu
      88             : 
      89        3688 :                   DO i = 1, 3
      90        2784 :                      forc_a3(i) = forc_a3(i) + (a3_1+a3_2)*gv(i)*atoms%zatom(n)
      91             :                   END DO
      92             : 
      93             :                END DO ! lh (0:sphhar%nlh(nd))
      94             : 
      95             :                ! Add onto existing forces.
      96        3616 :                DO i = 1, 3
      97        3616 :                   force(i,n,jsp) = force(i,n,jsp) + REAL(forc_a3(i))
      98             :                END DO
      99             : 
     100             :                ! Write out result.
     101         904 :                WRITE (oUnit,FMT=8010) n
     102         904 :                WRITE (oUnit,FMT=8020) (forc_a3(i),i=1,3)
     103             : 8010           FORMAT (' FORCES: EQUATION A3 FOR ATOM TYPE',i4)
     104             : 8020           FORMAT (' FX_A3=',2f10.6,' FY_A3=',2f10.6,' FZ_A3=',2f10.6)
     105             : 
     106             :             END IF ! atoms%l_geo(n)
     107             :          END DO ! n (1:atoms%ntype)
     108             :       END DO ! jsp (1:input%jpins)
     109             : 
     110         334 :       CALL timestop("force_a3")
     111             : 
     112         334 :    END SUBROUTINE force_a3
     113             : END MODULE m_force_a3