Line data    Source code

       1             : !--------------------------------------------------------------------------------
       2             : ! Copyright (c) 2020 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
       3             : ! This file is part of FLEUR and available as free software under the conditions
       4             : ! of the MIT license as expressed in the LICENSE file in more detail.
       5             : !--------------------------------------------------------------------------------
       6             : MODULE m_forceb8
       7             :    !-----------------------------------------------------------------------------
       8             :    ! Implements the surface contribution to the force. Madsen Eq.(B8)
       9             :    !
      10             :    ! FZJ 15/3-01 GMadsen
      11             :    !-----------------------------------------------------------------------------
      12             : CONTAINS
      13           0 :    SUBROUTINE force_b8(atoms,ecwk,stars, sym,cell, jspin, force,f_b8)
      14             : 
      15             :       USE m_constants, ONLY : tpi_const
      16             :       USE m_sphbes
      17             :       USE m_stern
      18             :       USE m_types
      19             : 
      20             :       IMPLICIT NONE
      21             : 
      22             :       TYPE(t_sym),   INTENT(IN) :: sym
      23             :       TYPE(t_stars), INTENT(IN) :: stars
      24             :       TYPE(t_cell),  INTENT(IN) :: cell
      25             :       TYPE(t_atoms), INTENT(IN) :: atoms
      26             : 
      27             :       INTEGER, INTENT(IN)    :: jspin
      28             :       COMPLEX, INTENT(IN)    :: ecwk(stars%ng3)
      29             :       COMPLEX, INTENT(INOUT) :: f_b8(3,atoms%ntype)
      30             :       REAL,    INTENT(INOUT) :: force(:,:,:)
      31             : 
      32           0 :       INTEGER g(3),nst,stg(3,sym%nop),ia,istr,i,j,jj,iType
      33           0 :       REAL    fj(0:atoms%lmaxd),rotkzz(3),rstg(3,sym%nop)
      34             :       REAL    frmt,gl,pha,s
      35           0 :       COMPLEX taup(sym%nop),factor,fact,fstar(3),fsur2(3)
      36             :    
      37           0 :       DO iType=1,atoms%ntype
      38           0 :          ia = atoms%firstAtom(iType)
      39           0 :          frmt = 2.0*tpi_const*atoms%rmt(iType)**2
      40           0 :          fsur2(1:3) = cmplx(0.0,0.0)         
      41             : 
      42             :          ! Skip G=(0,0,0) [no contribution to ekin]
      43           0 :          DO istr=2,stars%ng3_fft 
      44           0 :             g(:)     = stars%kv3(:,istr)
      45           0 :             fstar(:) = cmplx(0.0,0.0)
      46           0 :             CALL stern(sym,cell,g, nst,stg,taup,gl,rstg)
      47             : 
      48           0 :             CALL sphbes(atoms%lmax(iType),atoms%rmt(iType)*gl,fj)
      49           0 :             fact = ecwk(istr) * fj(1) / gl
      50             : 
      51           0 :             DO jj=1,nst
      52             :                pha=(atoms%taual(1,ia)*stg(1,jj)+atoms%taual(2,ia)*stg(2,jj)&
      53           0 :                    +atoms%taual(3,ia)*stg(3,jj))*tpi_const
      54             : 
      55             :                ! Swapped sin and cos because there's an i in the equation
      56           0 :                factor = fact * cmplx(-sin(pha),cos(pha)) * taup(jj)
      57           0 :                DO i=1,3
      58           0 :                   fstar(i) = fstar(i) + factor*rstg(i,jj)
      59             :                END DO
      60             :             END DO
      61           0 :             DO i=1,3
      62           0 :                fsur2(i)=fsur2(i)+fstar(i)*frmt
      63             :             END DO
      64             :          END DO
      65           0 :          DO i=1,3
      66           0 :             f_b8(i,iType) = f_b8(i,iType) + fsur2(i)
      67           0 :             force(i,iType,jspin) = force(i,iType,jspin) + real(fsur2(i))
      68             :          END DO
      69             :       END DO
      70             : 
      71           0 :    END SUBROUTINE force_b8
      72             : END MODULE m_forceb8