Line data    Source code

       1             : MODULE m_dftUPotential
       2             : 
       3             :    USE m_constants
       4             :    USE m_juDFT
       5             :    USE m_types
       6             :    USE m_uj2f
       7             :    USE m_doubleCounting
       8             :    USE m_coulombPotential
       9             : 
      10             :    IMPLICIT NONE
      11             : 
      12             :    CONTAINS
      13             : 
      14         160 :    SUBROUTINE dftUPotential(density, ldau, jspins, equivalentAtoms, l_spinoffd, potential, ldaUEnergy, spinavg_dc)
      15             : 
      16             :       !This subroutine calculates the DFT+U potential matrix
      17             : 
      18             :       COMPLEX,          INTENT(IN)     :: density(-lmaxU_const:,-lmaxU_const:,:)
      19             :       TYPE(t_utype),    INTENT(IN)     :: ldau
      20             :       INTEGER,          INTENT(IN)     :: jspins
      21             :       INTEGER,          INTENT(IN)     :: equivalentAtoms
      22             :       LOGICAL,          INTENT(IN)     :: l_spinoffd
      23             :       COMPLEX,          INTENT(INOUT)  :: potential(-lmaxU_const:,-lmaxU_const:,:)
      24             :       REAL,             INTENT(INOUT)  :: ldaUEnergy
      25             :       LOGICAL, OPTIONAL,INTENT(IN)     :: spinavg_dc
      26             : 
      27             : 
      28             :       LOGICAL :: spinavg_dc_local
      29             :       INTEGER :: m,ispin,jspin
      30             :       REAL    :: u_htr,j_htr,f0,f2,f4,f6,energy_contribution,total_charge, double_counting
      31             :       COMPLEX, ALLOCATABLE :: Vdc(:,:,:)
      32             : 
      33         160 :       spinavg_dc_local = .false.
      34         320 :       if(present(spinavg_dc)) spinavg_dc_local = spinavg_dc .and. jspins==2
      35             : 
      36         160 :       call coulombPotential(density, ldau, jspins, l_spinoffd, potential, energy_contribution)
      37             : 
      38         160 :       CALL uj2f(jspins,ldau,f0,f2,f4,f6)
      39         160 :       u_htr = f0/hartree_to_ev_const
      40         160 :       IF (ldau%l.EQ.1) THEN
      41          52 :          j_htr = f2/(5*hartree_to_ev_const)
      42         108 :       ELSE IF (ldau%l.EQ.2) THEN
      43          96 :          j_htr = 1.625*f2/(14*hartree_to_ev_const)
      44          12 :       ELSE IF (ldau%l.EQ.3) THEN
      45          12 :          j_htr = (286.+195.*451./675.+250.*1001./2025.)*f2/(6435*hartree_to_ev_const)
      46             :       END IF
      47             : 
      48             :       !Add double counting terms
      49             :       Vdc = doubleCountingPot(density, ldau, u_htr, j_htr, l_spinoffd,&
      50         160 :                               .FALSE.,spinavg_dc_local,0.0)
      51       12472 :       potential = potential - Vdc
      52             : 
      53             :       double_counting = doubleCountingEnergy(density, ldau, u_htr, j_htr, l_spinoffd,&
      54         160 :                                              .FALSE.,spinavg_dc_local,0.0)
      55             : 
      56         160 :       if(ldau%l_amf) then
      57           0 :          energy_contribution = energy_contribution - double_counting
      58             :       else
      59         160 :          energy_contribution = energy_contribution - double_counting
      60             :       endif
      61             :       
      62         160 :       total_charge = 0.0
      63         376 :       DO ispin = 1, MIN(2,SIZE(density,3))
      64        1400 :          DO m = -ldau%l, ldau%l
      65        1240 :             total_charge = total_charge +  REAL(density(m,m,ispin))
      66             :          ENDDO
      67             :       ENDDO
      68         160 :       energy_contribution = energy_contribution - (u_htr-j_htr)/2.0 * total_charge
      69             :       
      70         160 :       ldaUEnergy = ldaUEnergy + energy_contribution * equivalentAtoms
      71             : 
      72         160 :    END SUBROUTINE dftUPotential
      73             : 
      74             : END MODULE m_dftUPotential