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_forcea21U
       7             : CONTAINS
       8          12 :    SUBROUTINE force_a21_U(atoms,itype,isp,we,ne,usdus,v_mmp,eigVecCoeffs,aveccof,bveccof,cveccof,a21)
       9             :       !--------------------------------------------------------------------------
      10             :       ! This subroutine calculates the lda+U contribution to the HF forces,
      11             :       ! similar to the A21 term, according to eqn. (22) of F. Tran et al.
      12             :       ! Comp.Phys.Comm. 179 (2008) 784-790
      13             :       !--------------------------------------------------------------------------
      14             : 
      15             :       USE m_constants
      16             :       USE m_types_setup
      17             :       USE m_types_usdus
      18             :       USE m_types_cdnval
      19             : 
      20             :       IMPLICIT NONE
      21             : 
      22             :       TYPE(t_usdus),        INTENT(IN) :: usdus
      23             :       TYPE(t_atoms),        INTENT(IN) :: atoms
      24             :       TYPE(t_eigVecCoeffs), INTENT(IN) :: eigVecCoeffs
      25             : 
      26             :       INTEGER, INTENT (IN)    :: itype,isp,ne
      27             : 
      28             :       REAL,    INTENT(IN)    :: we(:) !(ne)
      29             :       COMPLEX, INTENT(IN)    :: v_mmp(-lmaxU_const:,-lmaxU_const:,:) !(-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const,atoms%n_u+atoms%n_hia)
      30             :       COMPLEX, INTENT(IN)    :: aveccof(:,:,0:,:) !(3,ne,0:atoms%lmaxd*(atoms%lmaxd+2),atoms%nat)
      31             :       COMPLEX, INTENT(IN)    :: bveccof(:,:,0:,:) !(3,ne,0:atoms%lmaxd*(atoms%lmaxd+2),atoms%nat)
      32             :       COMPLEX, INTENT(IN)    :: cveccof(:,-atoms%llod:,:,:,:) !(3,-atoms%llod:atoms%llod,ne,atoms%nlod,atoms%nat)
      33             :       REAL,    INTENT(INOUT) :: a21(:,:) !(3,atoms%nat)
      34             : 
      35             :       COMPLEX v_a, v_b, v_c, p1, p2, p3
      36             :       INTEGER lo, lop, l, lp, mp, lm, lmp, iatom, ie, i, m, i_u
      37             : 
      38             :       !--- ABBREVIATIONS --------------------------------------------------------
      39             :       ! ccof       : coefficient of the local orbital function (u_lo*Y_lm)
      40             :       ! cveccof    : is defined equivalently to aveccof, but with the LO-fct.
      41             :       ! for information on nlo,llo,uulon,dulon, and uloulopn see
      42             :       ! comments in setlomap.
      43             :       !--------------------------------------------------------------------------
      44             : 
      45          60 :       DO i_u = 1, atoms%n_u+atoms%n_hia
      46             : 
      47          48 :          IF(atoms%lda_u(i_u)%atomType/=itype) CYCLE
      48             : 
      49          24 :          l = atoms%lda_u(i_u)%l
      50             : 
      51             :          ! Add contribution for the regular LAPWs (like force_a21, but with
      52             :          ! the potential matrix, v_mmp, instead of the tuu, tdd ...)
      53             : 
      54         120 :          DO m = -l,l
      55          96 :             lm = l* (l+1) + m
      56         528 :             DO mp = -l,l
      57         408 :                lmp = l* (l+1) + mp
      58         408 :                v_a = v_mmp(m,mp,i_u)
      59         408 :                v_b = v_mmp(m,mp,i_u) * usdus%ddn(l,itype,isp)
      60         912 :                DO iatom = atoms%firstAtom(itype), atoms%firstAtom(itype) + atoms%neq(itype) - 1
      61        6528 :                   DO ie = 1,ne
      62       23256 :                      DO i = 1,3
      63       17136 :                         p1 = (CONJG(eigVecCoeffs%abcof(ie,lm,0,iatom,isp)) * v_a) * aveccof(i,ie,lmp,iatom)
      64       17136 :                         p2 = (CONJG(eigVecCoeffs%abcof(ie,lm,1,iatom,isp)) * v_b) * bveccof(i,ie,lmp,iatom)
      65       22848 :                         a21(i,iatom) = a21(i,iatom) + 2.0*AIMAG(p1 + p2) * we(ie)/atoms%neq(itype)
      66             :                      END DO
      67             :                   END DO
      68             :                END DO
      69             :             END DO ! mp
      70             :          END DO   ! m
      71             : 
      72             :          ! If there are also LOs on this atom, with the same l as
      73             :          ! the one of LDA+U, add another few terms
      74             : 
      75          36 :           DO lo = 1,atoms%nlo(itype)
      76          48 :             IF (l == atoms%llo(lo,itype)) THEN
      77           0 :                DO m = -l,l
      78           0 :                   lm = l* (l+1) + m
      79           0 :                   DO mp = -l,l
      80           0 :                      lmp = l* (l+1) + mp
      81           0 :                      v_a = v_mmp(m,mp,i_u)
      82           0 :                      v_b = v_mmp(m,mp,i_u) * usdus%uulon(lo,itype,isp)
      83           0 :                      v_c = v_mmp(m,mp,i_u) * usdus%dulon(lo,itype,isp)
      84           0 :                      DO iatom =  atoms%firstAtom(itype), atoms%firstAtom(itype) + atoms%neq(itype) - 1
      85           0 :                         DO ie = 1,ne
      86           0 :                            DO i = 1,3
      87           0 :                               p1 = v_a * (CONJG(eigVecCoeffs%ccof(m,ie,lo,iatom,isp)) * cveccof(i,mp,ie,lo,iatom))
      88             :                               p2 = v_b * (CONJG(eigVecCoeffs%abcof(ie,lm,0,iatom,isp)) * cveccof(i,mp,ie,lo,iatom) + &
      89           0 :                                           CONJG(eigVecCoeffs%ccof(m,ie,lo,iatom,isp)) * aveccof(i,ie,lmp,iatom))
      90             :                               p3 = v_c * (CONJG(eigVecCoeffs%abcof(ie,lm,1,iatom,isp)) * cveccof(i,mp,ie,lo,iatom) + &
      91           0 :                                           CONJG(eigVecCoeffs%ccof(m,ie,lo,iatom,isp)) * bveccof(i,ie,lmp,iatom))
      92           0 :                               a21(i,iatom) = a21(i,iatom) + 2.0*AIMAG(p1 + p2 + p3)*we(ie)/atoms%neq(itype)
      93             :                            END DO
      94             :                         END DO
      95             :                      END DO
      96             :                   END DO
      97             :                END DO
      98             :             END IF   ! l == atoms%llo(lo,itype)
      99             :          END DO     ! lo = 1,atoms%nlo
     100             : 
     101             :       END DO
     102          12 :    END SUBROUTINE force_a21_U
     103             : END MODULE m_forcea21U