LCOV - code coverage report
Current view: top level - ldau - v_mmp.F90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 57 62 91.9 %
Date: 2019-09-08 04:53:50 Functions: 1 1 100.0 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------
       2             : ! Copyright (c) 2016 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             : 
       7             : MODULE m_vmmp
       8             :   !     ************************************************************
       9             :   !     This subroutine calculates the potential matrix v^{s}_{m,m'}
      10             :   !     for a given atom 'n' and l-quantum number 'l'. The l,u,j's for
      11             :   !     all atoms are stored in lda_u, the density matrix is ns_mmp,
      12             :   !     and the e-e- interaction potential is u(m1,m2,m3,m4,n).
      13             :   !     For details see Eq.(16) of Shick et al. PRB 60, 10765 (1999)
      14             :   !
      15             :   !     Additionally, the total energy contribution of LDA+U (Eq.24)
      16             :   !     is calculated (e_ldau).
      17             :   !     Part of the LDA+U package                   G.B., Oct. 2000
      18             :   !     
      19             :   !     Extension to multiple U per atom type  G.M. 2017
      20             :   !     ************************************************************
      21             : CONTAINS
      22          28 :   SUBROUTINE v_mmp(sym,atoms,jspins,ns_mmp,u,f0,f2, vs_mmp,results)
      23             : 
      24             :     USE m_types
      25             :     USE m_constants
      26             :     IMPLICIT NONE
      27             :     TYPE(t_sym),INTENT(IN)          :: sym
      28             :     TYPE(t_results),INTENT(INOUT)   :: results
      29             :     TYPE(t_atoms),INTENT(IN)        :: atoms
      30             :     !
      31             :     ! ..  Arguments ..
      32             :     INTEGER, INTENT(IN)    :: jspins 
      33             :     REAL,    INTENT(IN)    :: u(-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const,&
      34             :                                 -lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const,atoms%n_u)
      35             :     REAL,    INTENT(IN)    :: f0(atoms%n_u),f2(atoms%n_u)
      36             :     COMPLEX, INTENT(OUT)   :: vs_mmp(-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const,atoms%n_u,jspins)
      37             : 
      38             :     COMPLEX, INTENT(INOUT) :: ns_mmp(-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const,atoms%n_u,jspins)
      39             : 
      40             :     ! ..  Local Variables ..
      41             :     INTEGER ispin,jspin,l ,mp,p,q,itype,m,i_u
      42             :     REAL rho_tot,u_htr,j_htr,e_ee,ns_sum,spin_deg,e_dc,e_dcc
      43          56 :     REAL rho_sig(jspins),v_diag(jspins),eta(0:jspins)
      44             :     !
      45             :     ! Use around-mean-field limit (true) of atomic limit (false)
      46             :     !
      47             :     !
      48             :     ! Loop over atoms
      49             :     !
      50          28 :     spin_deg = 1.0 / (3 - jspins)
      51          28 :     results%e_ldau = 0.0
      52             : 
      53         140 :     DO i_u = 1, atoms%n_u
      54         112 :        iType = atoms%lda_u(i_u)%atomType
      55         112 :        l = atoms%lda_u(i_u)%l
      56         112 :        u_htr = atoms%lda_u(i_u)%u / hartree_to_ev_const
      57         112 :        j_htr = atoms%lda_u(i_u)%j / hartree_to_ev_const
      58         112 :        u_htr = f0(i_u)/hartree_to_ev_const
      59         112 :        IF (l.EQ.1) THEN
      60          56 :           j_htr = f2(i_u)/(5*hartree_to_ev_const)
      61          56 :        ELSE IF (l.EQ.2) THEN
      62          56 :           j_htr = 1.625*f2(i_u)/(14*hartree_to_ev_const)
      63           0 :        ELSE IF (l.EQ.3) THEN
      64           0 :           j_htr = (286.+195*451/675+250*1001/2025)*f2(i_u)/(6435*hartree_to_ev_const)
      65             :        END IF
      66             :        !
      67             :        ! calculate spin-density 'rho_sig' and total density 'rho_tot'
      68             :        !
      69         112 :        rho_tot = 0.0
      70         224 :        DO ispin = 1,jspins
      71         112 :           rho_sig(ispin) = 0.0
      72         560 :           DO m = -l,l
      73         560 :              rho_sig(ispin) = rho_sig(ispin) + REAL(ns_mmp(m,m,i_u,ispin))
      74             :           END DO
      75         224 :           rho_tot = rho_tot + rho_sig(ispin)
      76             :        END DO
      77         112 :        rho_sig(1) = rho_sig(1) * spin_deg  ! if jspins = 1, divide by 2
      78         112 :        IF (atoms%lda_u(i_u)%l_amf) THEN
      79           0 :           eta(1) = rho_sig(1) / (2*l + 1) 
      80           0 :           eta(jspins) = rho_sig(jspins) / (2*l + 1) 
      81           0 :           eta(0) = (eta(1) + eta(jspins) ) / 2
      82             :        ELSE
      83         112 :           eta(0) = 1.0
      84         112 :           eta(1) = 1.0
      85         112 :           eta(jspins) = 1.0
      86             :        END IF
      87             :        !
      88             :        !--------------------------------------------------------------------------------------------+
      89             :        !  s       --                                        s'                    1        s   1    |
      90             :        ! V     =  >  ( <m,p|V|m',q> - <m,p|V|q,m'> d     ) n     + d    ( -U (n - -) + J (n  - -) ) |
      91             :        !  m,m'    --                                s,s'    p,q     m,m'          2            2    |
      92             :        !        p,q,s'                                                                              |
      93             :        !--------------------------------------------------------------------------------------------+     
      94             :        ! initialise vs_mmp
      95             :        !
      96             :        !IF (sym%invs) THEN
      97             :        !   vs_mmp(:,:,i_u,:) = ns_mmp(:,:,i_u,:)
      98             :        !   DO ispin = 1,jspins
      99             :        !      DO m = -l,l
     100             :        !         DO mp = -l,l
     101             :        !            ns_mmp(m,mp,i_u,ispin) = vs_mmp(-m,-mp,i_u,ispin)
     102             :        !         END DO
     103             :        !      END DO
     104             :        !   END DO
     105             :        !END IF
     106         112 :        vs_mmp(:,:,i_u,:) = CMPLX(0.0,0.0)
     107             :        !
     108             :        ! outer spin loop - set up v_mmp
     109             :        !
     110         224 :        DO ispin = 1,jspins
     111         560 :           DO m = -l,l
     112        4368 :              DO mp =-l,l 
     113        6160 :                 DO jspin = 1,jspins
     114        1904 :                    IF (ispin.EQ.jspin) THEN
     115       18928 :                       DO p = -l,l
     116       89488 :                          DO q = -l,l
     117             :                             vs_mmp(m,mp,i_u,ispin) = vs_mmp(m,mp,i_u,ispin) +  &
     118       48048 :                                ns_mmp(p, q,i_u,jspin) * ( u(m,p,mp,q,i_u) - u(m,p,q,mp,i_u) ) 
     119             :                          END DO
     120             :                       END DO
     121             :                    END IF
     122        3808 :                    IF ((ispin.NE.jspin).OR.(jspins.EQ.1)) THEN
     123       18928 :                       DO p = -l,l
     124       89488 :                          DO q = -l,l
     125             :                             vs_mmp(m,mp,i_u,ispin) = vs_mmp(m,mp,i_u,ispin) +  &
     126       48048 :                                u(m,p,mp,q,i_u) * ns_mmp(p, q,i_u,jspin) 
     127             :                          END DO
     128             :                       END DO
     129             :                    END IF
     130             :                 END DO
     131             : 
     132             :              END DO ! m' loop
     133             :           END DO   ! m  loop
     134             :        END DO      ! outer spin loop
     135             :        !
     136             :        !  set diagonal terms and correct for non-spin-polarised case
     137             :        !
     138         336 :        DO ispin = 1,jspins
     139         112 :           v_diag(ispin) = - u_htr * ( rho_tot - 0.5*eta(0) ) + j_htr * ( rho_sig(ispin) - 0.5*eta(ispin) )
     140         672 :           DO m = -l,l
     141        4256 :              DO mp = -l,l
     142        2352 :                 vs_mmp(m,mp,i_u,ispin) = vs_mmp(m,mp,i_u,ispin) * spin_deg
     143             :              END DO
     144         560 :              vs_mmp(m,m,i_u,ispin) = vs_mmp(m,m,i_u,ispin) + v_diag(ispin)
     145             :           END DO
     146             :        END DO
     147             : 
     148             :        !----------------------------------------------------------------------+
     149             :        !              s                                                       !
     150             :        !  ee      1  ---   s        s                     1        s  1       !
     151             :        ! E  (n) = -  >    n      ( V     + d     ( U (n - -) - J (n - -) ))   !
     152             :        !          2  ---   m,m'     m,m'    m,m'          2           2       !
     153             :        !             m,m'                                                     !
     154             :        !----------------------------------------------------------------------+
     155             : 
     156             :        e_ee = 0.0
     157         336 :        DO ispin = 1,jspins
     158         560 :           DO m = -l,l
     159        4256 :              DO mp =-l,l
     160        2352 :                 e_ee=e_ee+REAL(vs_mmp(m,mp,i_u,ispin)*ns_mmp(m,mp,i_u,ispin))
     161             :              END DO
     162         560 :              e_ee = e_ee - v_diag(ispin) * REAL( ns_mmp(m,m,i_u,ispin) )
     163             :           END DO
     164             :        END DO
     165             : 
     166             :        !----------------------------------------------------------------------+
     167             :        !   dc       ee      U           J  --   s   s       1                 |
     168             :        !  E      = E  (n) - - n (n-1) + -  >   n  (n -1)  - - (U-J) n         |
     169             :        !   LDA+U            2           2  --               2                 |
     170             :        !                                    s                                 |
     171             :        !----------------------------------------------------------------------+
     172             : 
     173             :        ns_sum = 0.0
     174         336 :        DO ispin = 1,jspins
     175         224 :           ns_sum = ns_sum + rho_sig(ispin) * (rho_sig(ispin) - eta(ispin))
     176             :        END DO
     177         112 :        e_dc = u_htr * rho_tot * ( rho_tot - eta(0) ) - j_htr * ns_sum
     178         112 :        e_dcc = (u_htr - j_htr) * rho_tot
     179             : 
     180         112 :        ns_sum = ns_sum / spin_deg
     181             :        !       e_ldau = e_ldau + (e_ee -  u_htr * rho_tot * ( rho_tot - 1. ) 
     182             :        !    +    + j_htr * ns_sum  - (u_htr - j_htr) * rho_tot) * neq(itype)
     183             :        !       write(*,*) e_ldau
     184         140 :        results%e_ldau = results%e_ldau + ( e_ee - e_dc - e_dcc) * atoms%neq(itype)
     185             :        !       write(*,*) e_ldau
     186             : 
     187             :     END DO ! loop over U parameters
     188             : 
     189          28 :     results%e_ldau = results%e_ldau / 2
     190             : 
     191          28 :   END SUBROUTINE v_mmp
     192             : END MODULE m_vmmp

Generated by: LCOV version 1.13