TABLE OF CONTENTS


m_paw_correlations/calc_ubare [ Functions ]

[ Top ] [ m_paw_correlations ] [ Functions ]

NAME

 calc_ubare

FUNCTION

 Calculate the bare interaction on atomic orbitals

INPUTS

  itypatcor = value of itypat for correlated species
  lpawu = angular momentum for correlated species
  pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data:
  pawang
     %lmax=Maximum value of angular momentum l+1
     %gntselect((2*l_max-1)**2,l_max**2,l_max**2)=
                     selection rules for Gaunt coefficients
  pawrad(ntypat) <type(pawrad_type)>=paw radial mesh and related data:
     %mesh_size=Dimension of radial mesh
     %rad(mesh_size)=The coordinates of all the points of the radial mesh
     %radfact(mesh_size)=Factor used to compute radial integrals

OUTPUT

NOTES

PARENTS

      pawpuxinit

CHILDREN

      poisson,simp_gen,wrtout

SOURCE

2637  subroutine calc_ubare(itypatcor,lpawu,pawang,pawrad,pawtab,rmax)
2638 
2639 
2640 !This section has been created automatically by the script Abilint (TD).
2641 !Do not modify the following lines by hand.
2642 #undef ABI_FUNC
2643 #define ABI_FUNC 'calc_ubare'
2644 !End of the abilint section
2645 
2646  implicit none
2647 !Arguments ------------------------------------
2648  integer, intent(in)   :: itypatcor,lpawu
2649  type(pawang_type),intent(in) :: pawang
2650  type(pawrad_type),intent(in) :: pawrad
2651  type(pawtab_type),target,intent(in) :: pawtab
2652  real(dp), optional, intent(in) :: rmax
2653 
2654 !Local variables ------------------------------
2655 !scalars
2656  integer :: ilmn,ilmn1,iln,iln1,isel,isel1,itypat,jlmn,jlmn1,jln,jln1
2657  integer :: klm,klm1,klmn,klmn1,ll,lm0
2658  integer :: lmin,lmax,lmn2_size,mesh_size,meshsz,mm
2659  real(dp) :: norm,r_for_intg,rg,rg1,ubare,uint,uint_tmp
2660  character(len=800) :: message
2661 !arrays
2662  real(dp),allocatable :: ff(:),gg(:),phiphj(:),phiphj1(:)
2663 
2664 !************************************************************************
2665 
2666  itypat=itypatcor
2667  write(message,'(11a,f12.4,2a,i7,2a,f12.4,2a,i7,2a,f12.4)') &
2668 & ch10," =======================================================================",ch10, &
2669 & "  == Calculation of diagonal bare Coulomb interaction on ATOMIC orbitals ",ch10, &
2670 & "     (it is assumed that the wavefunction for the first reference ",ch10, &
2671 & "             energy in PAW atomic data is an atomic eigenvalue)",ch10,ch10, &
2672 & " Max value of the radius in atomic data file   =", pawrad%rmax ,ch10, &
2673 & " Max value of the mesh   in atomic data file   =", pawrad%mesh_size,ch10, &
2674 & " PAW radius is                                 =", pawtab%rpaw,ch10, &
2675 & " PAW value of the mesh for integration is      =", pawrad%int_meshsz,ch10, &
2676 & " Integral of atomic wavefunction until rpaw    =", pawtab%ph0phiint(1)
2677  if(.not.present(rmax)) then
2678    call wrtout(ab_out,message,'COLL')
2679    call wrtout(std_out,message,'COLL')
2680  end if
2681 
2682  mesh_size=pawrad%mesh_size
2683 
2684 !  Definition of the mesh used for integration.
2685  if(present(rmax)) then
2686    if(rmax>pawrad%rmax)  then
2687      write(message, '(a)' ) 'calc_ubare: the radius cannot be larger than the maximum radius of the mesh'
2688      MSG_ERROR(message)
2689    end if
2690    meshsz=pawrad_ifromr(pawrad,rmax)+5
2691    r_for_intg=rmax
2692  else
2693    meshsz=pawtab%partialwave_mesh_size
2694    r_for_intg=pawrad%rad(meshsz)  ! (we could use r_for_intg=-1)
2695  end if
2696 
2697  lmn2_size=pawtab%lmn2_size
2698  ABI_ALLOCATE(ff,(mesh_size))
2699  ABI_ALLOCATE(gg,(mesh_size))
2700  ABI_ALLOCATE(phiphj,(mesh_size))
2701  ABI_ALLOCATE(phiphj1,(mesh_size))
2702  do klmn=1,lmn2_size
2703    ilmn=pawtab%indklmn(7,klmn);jlmn=pawtab%indklmn(8,klmn)
2704    ! Select lpawu and first projectors il=jl=lpawu and first proj only
2705    if (( pawtab%indklmn(3,klmn)+pawtab%indklmn(4,klmn)==2*lpawu).and. &
2706 &   (-pawtab%indklmn(3,klmn)+pawtab%indklmn(4,klmn)==2*lpawu).and. &
2707 &   (pawtab%indlmn(3,ilmn)==1).and.(pawtab%indlmn(3,jlmn)==1) ) then
2708      klm=pawtab%indklmn(1,klmn);iln=pawtab%indlmn(5,ilmn);jln=pawtab%indlmn(5,jlmn)
2709      lmin=pawtab%indklmn(3,klmn);lmax=pawtab%indklmn(4,klmn)
2710      phiphj(1:meshsz)=pawtab%phi(1:meshsz,iln)*pawtab%phi(1:meshsz,jln)
2711      !write(6,*) "A",klmn,pawtab%klmntomn(1,klmn),pawtab%klmntomn(2,klmn),&
2712      !&pawtab%indklmn(7,klmn),pawtab%indklmn(8,klmn),pawtab%klmntomn(3,klmn),pawtab%klmntomn(4,klmn)
2713      do ll=lmin,lmin,2
2714        lm0=ll*ll+ll+1
2715        ff(1:meshsz)=phiphj(1:meshsz)
2716        call simp_gen(norm,ff,pawrad,r_for_intg=r_for_intg)
2717        call poisson(ff,ll,pawrad,gg)
2718        do klmn1=klmn,lmn2_size
2719          ilmn1=pawtab%indklmn(7,klmn);jlmn1=pawtab%indklmn(8,klmn)
2720          ! Select lpawu and first projectors il=jl=lpawu and first proj only
2721          if (( pawtab%indklmn(3,klmn1)+pawtab%indklmn(4,klmn1)==2*lpawu).and. &
2722 &         (-pawtab%indklmn(3,klmn1)+pawtab%indklmn(4,klmn1)==2*lpawu).and. &
2723 &         (pawtab%indlmn(3,ilmn1)==1).and.(pawtab%indlmn(3,jlmn1)==1) ) then
2724            !write(6,*) "A1",klmn1,pawtab%klmntomn(1,klmn1),pawtab%klmntomn(2,klmn1),&
2725            !&pawtab%indklmn(7,klmn1),pawtab%indklmn(8,klmn1),pawtab%klmntomn(3,klmn1),pawtab%klmntomn(4,klmn1)
2726            klm1=pawtab%indklmn(1,klmn1);iln1=pawtab%indlmn(5,ilmn1);jln1=pawtab%indlmn(5,jlmn1)
2727            phiphj1(1:meshsz)=pawtab%phi(1:meshsz,iln1)*pawtab%phi(1:meshsz,jln1)
2728            uint_tmp=zero
2729            if ((ll==lmin)) then
2730              ff(1)=zero
2731              ff(2:meshsz)=phiphj1(2:meshsz)*gg(2:meshsz)*two/pawrad%rad(2:meshsz)
2732              call simp_gen(uint_tmp,ff,pawrad,r_for_intg=r_for_intg)
2733            end if
2734            uint=zero
2735            do mm=-ll,ll
2736              isel =pawang%gntselect(lm0+mm,klm)
2737              isel1=pawang%gntselect(lm0+mm,klm1)
2738              if (isel>0.and.isel1>0) then
2739                rg =pawang%realgnt(isel)
2740                rg1=pawang%realgnt(isel1)
2741                uint=uint+uint_tmp*rg*rg1*two_pi
2742              end if
2743            end do
2744            if((pawtab%indklmn(5,klmn)==pawtab%indklmn(6,klmn)).and.&
2745 &           (pawtab%indklmn(5,klmn1)==pawtab%indklmn(6,klmn1)).and.&
2746 &           (pawtab%indklmn(5,klmn)==pawtab%indklmn(5,klmn1))) then
2747              ubare=uint*Ha_eV
2748            end if
2749          end if
2750        end do
2751      end do
2752    end if
2753  end do
2754  ABI_DEALLOCATE(gg)
2755  ABI_DEALLOCATE(ff)
2756  ABI_DEALLOCATE(phiphj)
2757  ABI_DEALLOCATE(phiphj1)
2758 
2759  write(message,'(a,3(a,f12.4,a),2a,f12.4,a)') ch10," For an atomic wfn truncated at rmax =",r_for_intg,ch10,&
2760 & "     The norm of the wfn is                    =",norm,ch10,&
2761 & "     The bare interaction (no renormalization) =",ubare," eV",ch10,&
2762 & "     The bare interaction (for a renorm. wfn ) =",ubare/norm/norm," eV"
2763  call wrtout(ab_out,message,'COLL')
2764  call wrtout(std_out,message,'COLL')
2765  if(r_for_intg < 10_dp .and. .not.present(rmax)) then
2766    write(message,'(a,f6.2,4a)') '   ( WARNING: The radial mesh in the atomic data file is cut at',r_for_intg,ch10,&
2767 &   '   Use XML atomic data files to compute the bare Coulomb interaction',ch10,&
2768 &   '   on a true normalized atomic wavefunction )'
2769    call wrtout(ab_out,message,'COLL')
2770    call wrtout(std_out,message,'COLL')
2771  end if
2772  if(present(rmax)) then
2773    write(message,'(2a)')  " =======================================================================",ch10
2774    call wrtout(ab_out,message,'COLL')
2775    call wrtout(std_out,message,'COLL')
2776  end if
2777 
2778  end subroutine calc_ubare

m_paw_correlations/m_paw_correlations [ Modules ]

[ Top ] [ Modules ]

NAME

  m_paw_correlations

FUNCTION

  This module contains several routines related to the treatment of electronic
    correlations in the PAW approach (DFT+U, exact-exchange, ...).

COPYRIGHT

 Copyright (C) 2018-2018 ABINIT group (BA,FJ,MT)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .

SOURCE

18 #if defined HAVE_CONFIG_H
19 #include "config.h"
20 #endif
21 
22 #include "abi_common.h"
23 
24 MODULE m_paw_correlations
25 
26  use defs_basis
27  use defs_abitypes
28  use m_errors
29  use m_abicore
30  use m_xmpi
31 
32  use m_linalg_interfaces
33  use m_special_funcs
34  use m_io_tools, only : open_file
35 
36  use m_pawang,      only : pawang_type,pawang_init,pawang_free
37  use m_pawrad,      only : pawrad_type,simp_gen,nderiv_gen,pawrad_ifromr,poisson
38  use m_pawtab,      only : pawtab_type
39  use m_pawrhoij,    only : pawrhoij_type
40  use m_paw_ij,      only : paw_ij_type
41  use m_paw_sphharm, only : mat_mlms2jmj,mat_slm2ylm
42  use m_paw_io,      only : pawio_print_ij
43  use m_paral_atom,  only : get_my_atmtab,free_my_atmtab
44 
45  implicit none
46 
47  private
48 
49 !public procedures.
50  public :: pawpuxinit   ! Initialize some data for PAW+U/PAW+LocalExactExchange/PAW+DMFT
51  public :: pawuenergy   ! Compute contributions to energy for PAW+U
52  public :: pawxenergy   ! Compute contributions to energy for PAW+[local exact exchange]
53  public :: setnoccmmp   ! Compute LDA+U density matrix nocc_{m,m_prime} or impose it
54  public :: setrhoijpbe0 ! Impose value of rhoij for using an auxiliairy file (PBE0 only)
55  public :: calc_ubare   ! Calculate the bare interaction on atomic orbitals
56 
57 CONTAINS  !========================================================================================

m_paw_correlations/pawpuxinit [ Functions ]

[ Top ] [ m_paw_correlations ] [ Functions ]

NAME

 pawpuxinit

FUNCTION

 Initialize some starting values of several arrays used in
 PAW+U/+DMFT or local exact-exchange calculations

 A-define useful indices for LDA+U/local exact-exchange
 B-Compute overlap between atomic wavefunction
 C-Compute matrix elements of coulomb interaction (see PRB vol.52 5467) [[cite:Liechenstein1995]]
    (angular part computed from Gaunt coefficients)

INPUTS

  dmatpuopt= select expression for the density matrix
  exchmix= mixing factor for local exact-exchange
  jpawu(ntypat)= value of J
  llexexch(ntypat)= value of l on which local exact-exchange applies
  llpawu(ntypat)= value of l on which LDA+U applies
  ntypat=number of types of atoms in unit cell.
  pawang <type(pawang_type)>=paw angular mesh and related data
     %lmax=Maximum value of angular momentum l+1
     %gntselect((2*l_max-1)**2,l_max**2,l_max**2)=
                     selection rules for Gaunt coefficients
  pawprtvol=output printing level for PAW
  pawrad(ntypat) <type(pawrad_type)>=paw radial mesh and related data
  pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data:
  upawu(ntypat)= value of U
  use_dmft = 0 no PAW+DMFT, =1 PAW+DMFT
  useexexch= 0 if no local exact-exchange; 1 if local exact-exchange
  usepawu= 0 if no LDA+U; 1 if LDA+U

OUTPUT

  pawtab <type(pawtab_type)>=paw tabulated data read at start:
     %euijkl=(2,2,lmn2_size,lmn2_size) array for computing LDA+U terms without occupancies
     %ij_proj=nproj*(nproju+1)/2
     %klmntomn(4,lmn2_size) = Array giving im, jm ,in, and jn for each klmn=(ilmn,jlmn)
     %lnproju(nproj)= value of ln for projectors on which paw+u/local exact-exchange acts.
     %nproju=number of projectors for orbitals on which paw+u/local exact-exchange acts.
     %phiphjint(pawtabitypat%ij_proj)=Integral of Phi(:,i)*Phi(:,j) for correlated orbitals.
     %usepawu=0 if no LDA+U; 1 if LDA+U
     %useexexch=0 if no local exact-exchange; 1 if local exact-exchange
     === if usepawu>0
     %jpawu= value of J
     %upawu= value of U
     %vee(2*lpawu+1,:,:,:)=matrix of the screened interaction for correlated orbitals
     === if useexexch>0
     %fk
     %vex(2*lpawu+1,:,:,:)=matrix of the screened interaction for correlated orbitals

PARENTS

      bethe_salpeter,gstate,m_entropyDMFT,respfn,screening,sigma

CHILDREN

      calc_ubare,poisson,simp_gen,wrtout

SOURCE

120  subroutine pawpuxinit(dmatpuopt,exchmix,f4of2_sla,f6of2_sla,jpawu,llexexch,llpawu,&
121 &           ntypat,pawang,pawprtvol,pawrad,pawtab,upawu,use_dmft,useexexch,usepawu,&
122 &           ucrpa) ! optional argument
123 
124 
125 !This section has been created automatically by the script Abilint (TD).
126 !Do not modify the following lines by hand.
127 #undef ABI_FUNC
128 #define ABI_FUNC 'pawpuxinit'
129 !End of the abilint section
130 
131  implicit none
132 
133 !Arguments ---------------------------------------------
134 !scalars
135  integer,intent(in) :: dmatpuopt,ntypat,pawprtvol,use_dmft,useexexch,usepawu
136  real(dp),intent(in) :: exchmix
137  type(pawang_type), intent(in) :: pawang
138  integer,optional, intent(in) :: ucrpa
139 !arrays
140  integer,intent(in) :: llexexch(ntypat),llpawu(ntypat)
141  real(dp),intent(in) :: jpawu(ntypat),upawu(ntypat)
142  real(dp),intent(in) :: f4of2_sla(ntypat),f6of2_sla(ntypat)
143  type(pawrad_type),intent(inout) :: pawrad(ntypat)
144  type(pawtab_type),target,intent(inout) :: pawtab(ntypat)
145 
146 !Local variables ---------------------------------------
147 !scalars
148  integer :: icount,icountp,il,ilmn,ilmnp,isela,iselb,itemp,itypat,iu,iup,j0lmn,jl,jlmn,jlmnp,ju,jup
149  integer :: klm0u,klm0x,klma,klmb,klmn,klmna,klmnb,klmn0,klmn0p,kln,kln1,kln2,kyc,lcur,lexexch,lkyc,ll,ll1
150  integer :: lmexexch,lmkyc,lmn_size,lmn2_size,lmpawu,lpawu
151  integer :: m1,m11,m2,m21,m3,m31,m4,m41
152  integer :: mesh_size,int_meshsz,mkyc,sig,sigp,sz1
153  logical :: compute_euijkl,compute_euij_fll
154  real(dp) :: ak,fact,f4of2,f6of2,int1,intg,phiint_ij,phiint_ipjp,vee1,vee2
155  character(len=500) :: message
156 !arrays
157  integer,ABI_CONTIGUOUS pointer :: indlmn(:,:)
158  real(dp) :: euijkl_temp(2,2),euijkl_temp2(2,2),euijkl_dc(2,2)
159  real(dp),allocatable :: ff(:),fk(:),gg(:)
160 
161 ! *************************************************************************
162 
163  DBG_ENTER("COLL")
164 
165  if(useexexch==0.and.usepawu==0.and.use_dmft==0) return
166 
167 !PAW+U and local exact-exchange restriction
168  if(useexexch>0.and.usepawu>0)then
169    do itypat=1,ntypat
170      if (llpawu(itypat)/=llexexch(itypat).and.llpawu(itypat)/=-1.and.llexexch(itypat)/=-1) then
171        write(message, '(5a,i2,3a)' )&
172 &       '  When PAW+U (usepawu>0) and local exact-exchange (exexch>0)',ch10,&
173 &       '  are selected together, they must apply on the same',ch10,&
174 &       '  angular momentum (lpawu/=lexexch forbidden, here for typat=',itypat,') !',ch10,&
175 &       '  Action: correct your input file.'
176        MSG_ERROR(message)
177      end if
178    end do
179  end if
180 
181 !Print title
182  if((usepawu>=1.and.usepawu<=6).or.useexexch>0) write(message, '(3a)' ) ch10,ch10," ******************************************"
183  if(usepawu==1) then
184    write(message, '(3a)' ) trim(message),ch10," LDA+U Method used: FLL"
185  else if(usepawu==2) then
186    write(message, '(3a)' ) trim(message),ch10," LDA+U Method used: AMF"
187  else if(usepawu==3) then
188    write(message, '(3a)' ) trim(message),ch10," LDA+U Method used: AMF (alternative)"
189  else if(usepawu==4) then
190    write(message, '(3a)' ) trim(message),ch10," LDA+U Method used: FLL with no spin polarization in the xc functional"
191  else if(usepawu==5) then
192    write(message, '(3a)' ) trim(message),ch10," LDA+U Method used: FLL (no use of occupation matrix)"
193  else if(usepawu==6) then
194    write(message, '(3a)' ) trim(message),ch10," LDA+U Method used: AMF (no use of occupation matrix)"
195  end if
196  if(useexexch>0) write(message, '(3a)' ) trim(message),ch10," PAW Local Exact exchange: PBE0"
197  if((usepawu>=1.and.usepawu<=6).or.useexexch>0) &
198  write(message, '(3a)' ) trim(message),ch10," ******************************************"
199  if(use_dmft==0) then
200    call wrtout(ab_out,message,'COLL')
201    call wrtout(std_out,  message,'COLL')
202  end if
203 !if(use_dmft>0) then
204 !write(message, '(3a)' ) ch10, " (see DMFT data in log file) "
205 !call wrtout(ab_out,message,'COLL')
206 !endif
207 
208 !Loop on atom types
209  do itypat=1,ntypat
210    indlmn => pawtab(itypat)%indlmn
211    lmn_size=pawtab(itypat)%lmn_size
212    lmn2_size=pawtab(itypat)%lmn2_size
213    mesh_size=pawtab(itypat)%mesh_size
214    int_meshsz=pawrad(itypat)%int_meshsz
215    lcur=-1
216 
217 !  PAW+U data
218    if (usepawu>0.or.use_dmft>0) then
219      lcur=llpawu(itypat)
220      pawtab(itypat)%lpawu=lcur
221      if(lcur/=-1) then
222        pawtab(itypat)%usepawu=usepawu
223        pawtab(itypat)%upawu=upawu(itypat)
224        pawtab(itypat)%jpawu=jpawu(itypat)
225        pawtab(itypat)%f4of2_sla=f4of2_sla(itypat)
226        pawtab(itypat)%f6of2_sla=f6of2_sla(itypat)
227      else
228        pawtab(itypat)%usepawu=0
229        pawtab(itypat)%upawu=zero
230        pawtab(itypat)%jpawu=zero
231        pawtab(itypat)%f4of2_sla=zero
232        pawtab(itypat)%f6of2_sla=zero
233      end if
234    end if
235 
236 !  Local exact-echange data
237    if (useexexch>0) then
238      lcur=llexexch(itypat)
239      pawtab(itypat)%lexexch=lcur
240      pawtab(itypat)%exchmix=exchmix
241      if(pawtab(itypat)%lexexch==-1) pawtab(itypat)%useexexch=0
242      if(pawtab(itypat)%lexexch/=-1) pawtab(itypat)%useexexch=useexexch
243    end if
244 
245 !  Select only atoms with +U
246    if(lcur/=-1) then
247 
248 !    Compute number of projectors for LDA+U/local exact-exchange/LDA+DMFT
249      icount=count(indlmn(1,1:lmn_size)==lcur)
250      pawtab(itypat)%nproju=icount/(2*lcur+1)
251      if(useexexch>0.and.pawtab(itypat)%nproju>2)  then
252        write(message, '(a,a,a)' )&
253 &       '  Error on the number of projectors ',ch10,&
254 &       '  more than 2 projectors is not allowed for local exact-exchange'
255        MSG_ERROR(message)
256      end if
257      if(pawtab(itypat)%nproju*(2*lcur+1)/=icount)  then
258        message = 'pawpuxinit: Error on the number of projectors '
259        MSG_BUG(message)
260      end if
261      write(message, '(a,a,i4,a,a,i4)' ) ch10,&
262 &     ' pawpuxinit : for species ',itypat,ch10,&
263 &     '   number of projectors is',pawtab(itypat)%nproju
264      call wrtout(std_out,message,'COLL')
265 
266      pawtab(itypat)%ij_proj=pawtab(itypat)%nproju*(pawtab(itypat)%nproju+1)/2
267 
268 !    ==================================================
269 !    A-define useful indexes
270 !    --------------------------------------------------
271      if (allocated(pawtab(itypat)%lnproju)) then
272        ABI_DEALLOCATE(pawtab(itypat)%lnproju)
273      end if
274      ABI_ALLOCATE(pawtab(itypat)%lnproju,(pawtab(itypat)%nproju))
275      icount=0
276      do ilmn=1,lmn_size
277        if(indlmn(1,ilmn)==lcur) then
278          itemp=icount/(2*lcur+1)
279          if (itemp*(2*lcur+1)==icount) then
280            pawtab(itypat)%lnproju(itemp+1)=indlmn(5,ilmn)
281          end if
282          icount=icount+1
283        end if
284      end do
285 
286      if (allocated(pawtab(itypat)%klmntomn)) then
287        ABI_DEALLOCATE(pawtab(itypat)%klmntomn)
288      end if
289      ABI_ALLOCATE(pawtab(itypat)%klmntomn,(4,lmn2_size))
290      do jlmn=1,lmn_size
291        jl= indlmn(1,jlmn)
292        j0lmn=jlmn*(jlmn-1)/2
293        do ilmn=1,jlmn
294          il= indlmn(1,ilmn)
295          klmn=j0lmn+ilmn
296          pawtab(itypat)%klmntomn(1,klmn)=indlmn(2,ilmn)+il+1
297          pawtab(itypat)%klmntomn(2,klmn)=indlmn(2,jlmn)+jl+1
298          pawtab(itypat)%klmntomn(3,klmn)=indlmn(3,ilmn)
299          pawtab(itypat)%klmntomn(4,klmn)=indlmn(3,jlmn)
300        end do
301      end do
302 
303 !    ==================================================
304 !    B-PAW+U: overlap between atomic wavefunctions
305 !    --------------------------------------------------
306 !    if (usepawu>0) then
307      if(dmatpuopt==1) then
308        write(message, '(4a)' ) ch10,&
309 &       ' pawpuxinit : dmatpuopt=1 ',ch10,&
310 &       '   PAW+U: dens. mat. constructed by projection on atomic wfn inside PAW augm. region(s)'
311        call wrtout(std_out,message,'COLL')
312        write(message, '(8a)' ) ch10,&
313 &       ' pawpuxinit: WARNING: Check that the first partial wave for lpawu:', ch10, &
314 &       '                      - Is an atomic eigenfunction  ',ch10, &
315 &       '                      - Is normalized ',ch10, &
316 &       '                      In other cases, choose dmatpuopt=2'
317        call wrtout(std_out,message,'COLL')
318      else if(dmatpuopt==2) then
319        write(message, '(6a)' ) ch10,&
320 &       ' pawpuxinit : dmatpuopt=2 ',ch10,&
321 &       '   PAW+U: dens. mat. constructed by selecting contribution',ch10,&
322 &       '          for each angular momentum to the density (inside PAW augm. region(s))'
323        call wrtout(std_out,message,'COLL')
324      else if(dmatpuopt==3) then
325        write(message, '(a,a,a,a,a,a)' ) ch10,&
326 &       ' pawpuxinit : dmatpuopt=3 ',ch10,&
327 &       '    PAW+U: dens. mat. constructed by projection on atomic wfn inside PAW augm. region(s)',ch10,&
328 &       '           and normalized inside PAW augm. region(s)'
329        call wrtout(std_out,message,'COLL')
330        write(message, '(6a)' ) ch10,&
331 &       ' pawpuxinit: WARNING: Check that the first partial wave for lpawu:', ch10, &
332 &       '                     is an atomic eigenfunction',ch10, &
333 &       '                     In the other case, choose dmatpuopt=2'
334        call wrtout(std_out,message,'COLL')
335      end if
336 
337      ABI_ALLOCATE(ff,(mesh_size))
338      ff(:)=zero
339 
340      if (allocated(pawtab(itypat)%ph0phiint)) then
341        ABI_DEALLOCATE(pawtab(itypat)%ph0phiint)
342      end if
343      if (allocated(pawtab(itypat)%zioneff)) then
344        ABI_DEALLOCATE(pawtab(itypat)%zioneff)
345      end if
346      ABI_ALLOCATE(pawtab(itypat)%ph0phiint,(pawtab(itypat)%nproju))
347      ABI_ALLOCATE(pawtab(itypat)%zioneff,(pawtab(itypat)%nproju))
348 
349      icount=0
350      do iu=1,pawtab(itypat)%nproju
351 !      write(std_out,*)'DJA iu',iu,' mesh_size',pawtab(itypat)%mesh_size
352 !      do ju=2,pawtab(itypat)%mesh_size
353 !      ff(ju)=pawtab(itypat)%phi(ju,pawtab(itypat)%lnproju(iu))/pawrad(itypat)%rad(ju)
354 !      write(std_out,fmt='(i5,3e15.5)')ju,pawrad(itypat)%rad(ju),ff(ju),&
355 !      &         RadFnH(pawrad(itypat)%rad(ju),4,3,15.0_dp)
356 !      end do
357 !      ff(1:mesh_size)=pawtab(itypat)%phi(1:mesh_size,pawtab(itypat)%lnproju(iu))**2
358 !      call simp_gen(int1,ff,pawrad(itypat))
359 !      write(std_out,*)'DJA iu',iu,'int1 ',int1
360 !      write(std_out,*)'DJA int1,IRadFnH',int1,IRadFnH(0.0_dp,pawrad(itypat)%rmax,4,3,12)
361 !      Calculation of zioneff
362        ju=pawtab(itypat)%mesh_size-1
363        ak=pawtab(itypat)%phi(ju,pawtab(itypat)%lnproju(iu))/pawtab(itypat)%phi(ju+1,pawtab(itypat)%lnproju(iu))
364        ak=ak*(pawrad(itypat)%rad(ju+1)/pawrad(itypat)%rad(ju))**(pawtab(itypat)%lpawu-1)
365        pawtab(itypat)%zioneff(iu)=log(ak)/(pawrad(itypat)%rad(ju+1)-pawrad(itypat)%rad(ju))
366 !      Calculation of ph0phiint
367        ff(1:mesh_size)=pawtab(itypat)%phi(1:mesh_size,pawtab(itypat)%lnproju(1))&
368 &       *pawtab(itypat)%phi(1:mesh_size,pawtab(itypat)%lnproju(iu))
369        call simp_gen(int1,ff,pawrad(itypat))
370        pawtab(itypat)%ph0phiint(iu)=int1
371      end do
372      if(abs(pawprtvol)>=2) then
373        do icount=1,pawtab(itypat)%nproju
374          write(message, '(a,a,i2,f9.5,a)' ) ch10,&
375 &         '  pawpuxinit: icount, ph0phiint(icount)=',icount,pawtab(itypat)%ph0phiint(icount)
376          call wrtout(std_out,message,'COLL')
377          write(message, '(a,f15.5)' ) &
378 &         '  pawpuxinit: zioneff=',pawtab(itypat)%zioneff(icount)
379          call wrtout(std_out,message,'COLL')
380        end do
381        write(message, '(a)' ) ch10
382        call wrtout(std_out,message,'COLL')
383      end if
384 
385      if (allocated(pawtab(itypat)%phiphjint)) then
386        ABI_DEALLOCATE(pawtab(itypat)%phiphjint)
387      end if
388      ABI_ALLOCATE(pawtab(itypat)%phiphjint,(pawtab(itypat)%ij_proj))
389 
390      icount=0
391      do ju=1,pawtab(itypat)%nproju
392        do iu=1,ju
393          icount=icount+1
394          if ((dmatpuopt==1).and.(useexexch==0)) then
395            pawtab(itypat)%phiphjint(icount)=pawtab(itypat)%ph0phiint(iu)*&
396 &           pawtab(itypat)%ph0phiint(ju)
397          else if((dmatpuopt==2).or.(useexexch>0)) then
398            ff(1:mesh_size)=pawtab(itypat)%phi(1:mesh_size,pawtab(itypat)%lnproju(iu))&
399 &           *pawtab(itypat)%phi(1:mesh_size,pawtab(itypat)%lnproju(ju))
400            call simp_gen(int1,ff,pawrad(itypat))
401            pawtab(itypat)%phiphjint(icount)=int1
402          else if((dmatpuopt>=3).and.(useexexch==0)) then
403            pawtab(itypat)%phiphjint(icount)=pawtab(itypat)%ph0phiint(iu)* &
404 &           pawtab(itypat)%ph0phiint(ju)/pawtab(itypat)%ph0phiint(1)**(dmatpuopt-2)
405          else
406            write(message, '(3a)' )&
407 &           '  PAW+U: dmatpuopt has a wrong value !',ch10,&
408 &           '  Action : change value in input file'
409            MSG_ERROR(message)
410          end if
411        end do
412      end do
413      if(pawtab(itypat)%ij_proj/=icount)  then
414        message = ' Error in the loop for calculating phiphjint '
415        MSG_ERROR(message)
416      end if
417      ABI_DEALLOCATE(ff)
418      if(abs(pawprtvol)>=2) then
419        do icount=1,pawtab(itypat)%ij_proj
420          write(message, '(a,a,i2,f9.5,a)' ) ch10,&
421 &         '  PAW+U: icount, phiphjint(icount)=',icount,pawtab(itypat)%phiphjint(icount)
422          call wrtout(std_out,message,'COLL')
423        end do
424      end if
425 !    end if
426 
427 !    ======================================================================
428 !    C-PAW+U: Matrix elements of coulomb interaction (see PRB vol.52 5467) [[cite:Liechenstein1995]]
429 !    1. angular part computed from Gaunt coefficients
430 !    --------------------------------------------------------------------
431      if (usepawu>0) then
432        lpawu=lcur
433 
434 !      a. compute F(k)
435 !      ---------------------------------------------
436        ABI_ALLOCATE(fk,(lpawu+1))
437        fk(1)=pawtab(itypat)%upawu
438 !      cf Slater Physical Review 165, p 665 (1968) [[cite:Slater1958]]
439 !      write(std_out,*) "f4of2_sla",pawtab(itypat)%f4of2_sla
440        if(lpawu==0) then
441          fk(1)=fk(1)
442        else if(lpawu==1) then
443          fk(2)=pawtab(itypat)%jpawu*5._dp
444        else if(lpawu==2) then
445 !        f4of2=0._dp
446          if(pawtab(itypat)%f4of2_sla<-0.1_dp)  then
447            f4of2=0.625_dp
448            pawtab(itypat)%f4of2_sla=f4of2
449          else
450            f4of2=pawtab(itypat)%f4of2_sla
451          end if
452          fk(2)=pawtab(itypat)%jpawu*14._dp/(One+f4of2)
453          fk(3)=fk(2)*f4of2
454 !        if(abs(pawprtvol)>=2) then
455          write(message,'(a,3x,a,f9.4,f9.4,f9.4,f9.4)') ch10,&
456 &         "Slater parameters F^0, F^2, F^4 are",fk(1),fk(2),fk(3)
457          call wrtout(std_out,message,'COLL')
458 !        end if
459        else if(lpawu==3) then
460          f4of2=0.6681_dp
461          f6of2=0.4943_dp
462          if(pawtab(itypat)%f4of2_sla<-0.1_dp)  then
463            f4of2=0.6681_dp
464            pawtab(itypat)%f4of2_sla=f4of2
465          else
466            f4of2=pawtab(itypat)%f4of2_sla
467          end if
468          if(pawtab(itypat)%f6of2_sla<-0.1_dp)  then
469            f6of2=0.4943_dp
470            pawtab(itypat)%f6of2_sla=f6of2
471          else
472            f6of2=pawtab(itypat)%f6of2_sla
473          end if
474          fk(2)=pawtab(itypat)%jpawu*6435._dp/(286._dp+195._dp*f4of2+250._dp*f6of2)
475          fk(3)=fk(2)*f4of2
476          fk(4)=fk(2)*f6of2
477          write(std_out,'(a,3x,a,f9.4,f9.4,f9.4,f9.4)') ch10,&
478 &         "Slater parameters F^0, F^2, F^4, F^6 are",fk(1),fk(2),fk(3),fk(4)
479        else
480          write(message, '(a,i0,2a)' )&
481 &         ' lpawu=',lpawu,ch10,&
482 &         ' lpawu not equal to 0 ,1 ,2 or 3 is not allowed'
483          MSG_ERROR(message)
484        end if
485 
486 !      b. Compute ak and vee.
487 !      ---------------------------------------------
488        if (allocated(pawtab(itypat)%vee)) then
489          ABI_DEALLOCATE(pawtab(itypat)%vee)
490        end if
491        sz1=2*lpawu+1
492        ABI_ALLOCATE(pawtab(itypat)%vee,(sz1,sz1,sz1,sz1))
493        pawtab(itypat)%vee=zero
494        lmpawu=(lpawu-1)**2+2*(lpawu-1)+1  ! number of m value below correlated orbitals
495        klm0u=lmpawu*(lmpawu+1)/2          ! value of klmn just below correlated orbitals
496 !      --------- 4 loops for interaction matrix
497        do m1=-lpawu,lpawu
498          m11=m1+lpawu+1
499          do m2=-lpawu,m1
500            m21=m2+lpawu+1
501 !          klma= number of pair before correlated orbitals +
502 !          number of pair for m1 lower than correlated orbitals
503 !          (m1+lpawu+1)*(lpawu-1) + number of pairs for correlated orbitals
504 !          before (m1,m2) + number of pair for m2 lower than current value
505            klma=klm0u+m11*lmpawu+(m11-1)*m11/2+m21
506            do m3=-lpawu,lpawu
507              m31=m3+lpawu+1
508              do m4=-lpawu,m3
509                m41=m4+lpawu+1
510                klmb=klm0u+m31*lmpawu+(m31-1)*m31/2+m41
511 !              --------- loop on k=1,2,3 (4 if f orbitals)
512                do kyc=1,2*lpawu+1,2
513                  lkyc=kyc-1
514                  lmkyc=(lkyc+1)*(lkyc)+1
515                  ak=zero
516                  do mkyc=-lkyc,lkyc,1
517                    isela=pawang%gntselect(lmkyc+mkyc,klma)
518                    iselb=pawang%gntselect(lmkyc+mkyc,klmb)
519                    if (isela>0.and.iselb>0) ak=ak +pawang%realgnt(isela)*pawang%realgnt(iselb)
520                  end do
521 !                ----- end loop on k=1,2,3 (4 if f orbitals)
522                  ak=ak/(two*dble(lkyc)+one)
523                  pawtab(itypat)%vee(m11,m31,m21,m41)=ak*fk(lkyc/2+1)+pawtab(itypat)%vee(m11,m31,m21,m41)
524                end do  !kyc
525                pawtab(itypat)%vee(m11,m31,m21,m41)=pawtab(itypat)%vee(m11,m31,m21,m41)*four_pi
526                pawtab(itypat)%vee(m21,m31,m11,m41)=pawtab(itypat)%vee(m11,m31,m21,m41)
527                pawtab(itypat)%vee(m11,m41,m21,m31)=pawtab(itypat)%vee(m11,m31,m21,m41)
528                pawtab(itypat)%vee(m21,m41,m11,m31)=pawtab(itypat)%vee(m11,m31,m21,m41)
529              end do
530            end do
531          end do
532        end do
533        ABI_DEALLOCATE(fk)
534 
535      !  testu=0
536      !  write(std_out,*) " Matrix of interaction vee(m1,m2,m1,m2)"
537      !  do m1=1,2*lpawu+1
538      !    write(std_out,'(2x,14(f12.6,2x))') (pawtab(itypat)%vee(m1,m2,m1,m2),m2=1,2*lpawu+1)
539      !    do m2=1,2*lpawu+1
540      !      testu=testu+ pawtab(itypat)%vee(m1,m2,m1,m2)
541      !   enddo
542      !  enddo
543      !  testu=testu/((two*lpawu+one)**2)
544      !  write(std_out,*) "------------------------"
545      !  write(std_out,'(a,f12.6)') " U=", testu
546      !  write(std_out,*) "------------------------"
547      !  write(std_out,*) " Matrix of interaction vee(m1,m2,m1,m2)-vee(m1,m2,m2,m1)"
548      !  do m1=1,2*lpawu+1
549      !    write(std_out,'(2x,14(f12.6,2x))') ((pawtab(itypat)%vee(m1,m2,m1,m2)-pawtab(itypat)%vee(m1,m2,m2,m1)),m2=1,2*lpawu+1)
550      !    do m2=1,2*lpawu+1
551      !    if(m1/=m2) testumj=testumj+ pawtab(itypat)%vee(m1,m2,m1,m2)-pawtab(itypat)%vee(m1,m2,m2,m1)
552      !   enddo
553      !  enddo
554      !  testumj=testumj/((two*lpawu)*(two*lpawu+one))
555      !  write(std_out,*) "------------------------"
556      !  write(std_out,'(a,f12.6)') " U-J=", testumj
557      !  write(std_out,*) "------------------------"
558      !  write(std_out,*) "------------------------"
559      !  write(std_out,'(a,f12.6)')  " J=", testu-testumj
560      !  write(std_out,*) "------------------------"
561 
562 !      c. For usepawu == 5 or 6, compute euijkl
563 !      ---------------------------------------------
564        compute_euijkl=(usepawu==5.or.usepawu==6)
565        if (compute_euijkl) then
566          if (allocated(pawtab(itypat)%euijkl)) then
567            ABI_DEALLOCATE(pawtab(itypat)%euijkl)
568          end if
569          ABI_ALLOCATE(pawtab(itypat)%euijkl,(2,2,lmn_size,lmn_size,lmn_size,lmn_size))
570          pawtab(itypat)%euijkl = zero
571          compute_euij_fll = .false.
572          euijkl_temp2=zero
573          if (usepawu==5) then ! Only for FLL
574            if (allocated(pawtab(itypat)%euij_fll)) then ! allocate euij_fll for FLL
575              ABI_DEALLOCATE(pawtab(itypat)%euij_fll)
576            end if
577            ABI_ALLOCATE(pawtab(itypat)%euij_fll,(lmn2_size))
578            pawtab(itypat)%euij_fll = zero
579            compute_euij_fll = .true.
580          end if
581 
582 !        loop on i,j
583          do klmna=1,lmn2_size
584            ilmn=pawtab(itypat)%indklmn(7,klmna) ! i
585            jlmn=pawtab(itypat)%indklmn(8,klmna) ! j
586            if (pawtab(itypat)%indlmn(1,ilmn)==lpawu.and.pawtab(itypat)%indlmn(1,jlmn)==lpawu) then ! only correlated orbitals
587              iu = pawtab(itypat)%indlmn(3,ilmn) ! ni
588              ju = pawtab(itypat)%indlmn(3,jlmn) ! nj
589              phiint_ij = pawtab(itypat)%phiphjint(iu+(ju*(ju-1))/2) ! iu <= ju by construction (ilmn<=jlmn)
590              m2 = pawtab(itypat)%indlmn(2,ilmn) ! mi
591              m21=m2+lpawu+1
592              m1 = pawtab(itypat)%indlmn(2,jlmn) ! mj
593              m11=m1+lpawu+1
594 
595              if (compute_euij_fll.and.m1==m2) then ! FLL
596                pawtab(itypat)%euij_fll(klmna) = - half * phiint_ij * ( pawtab(itypat)%jpawu - pawtab(itypat)%upawu )
597              end if
598 
599 !            loop on ip,jp (=k,l)
600              do klmnb=1,lmn2_size
601                ilmnp=pawtab(itypat)%indklmn(7,klmnb) ! ip (=k)
602                jlmnp=pawtab(itypat)%indklmn(8,klmnb) ! jp (=l)
603                if (pawtab(itypat)%indlmn(1,ilmnp)==lpawu.and.pawtab(itypat)%indlmn(1,jlmnp)==lpawu) then ! correlated orbitals
604                  iup = pawtab(itypat)%indlmn(3,ilmnp) ! nip
605                  jup = pawtab(itypat)%indlmn(3,jlmnp) ! njp
606                  phiint_ipjp = pawtab(itypat)%phiphjint(iup+(jup*(jup-1))/2) ! iup <= jup by construction (ilmnp<=jlmnp)
607                  m4 = pawtab(itypat)%indlmn(2,ilmnp) ! mip
608                  m41=m4+lpawu+1
609                  m3 = pawtab(itypat)%indlmn(2,jlmnp) ! mjp
610                  m31=m3+lpawu+1
611 
612                  euijkl_dc = zero
613 !                Compute the double-counting part of euijkl (invariant when exchanging i<-->j or ip<-->jp)
614                  if (m1==m2.and.m3==m4) then ! In that case, we have to add the double-counting term
615 
616                      do sig=1,2
617                        do sigp=1,2
618 
619                          if (usepawu==1.or.usepawu==5) then ! FLL
620 
621                            if (sig==sigp) then
622                              euijkl_dc(sig,sigp) = &
623 &                             phiint_ij * phiint_ipjp * ( pawtab(itypat)%upawu - pawtab(itypat)%jpawu )
624                            else
625                              euijkl_dc(sig,sigp) = &
626 &                             phiint_ij * phiint_ipjp * pawtab(itypat)%upawu
627                            end if
628 
629                          else if (usepawu==2.or.usepawu==6) then ! AMF
630 
631                            if (sig==sigp) then
632                              euijkl_dc(sig,sigp) = &
633 &                             two*lpawu/(two*lpawu+one) * phiint_ij * phiint_ipjp * ( pawtab(itypat)%upawu - pawtab(itypat)%jpawu )
634                            else
635                              euijkl_dc(sig,sigp) = &
636 &                             phiint_ij * phiint_ipjp * pawtab(itypat)%upawu
637                            end if
638 
639                          end if
640 
641                        end do ! sigp
642                      end do ! sig
643 
644                  end if ! double-counting term
645 
646                  vee1 = pawtab(itypat)%vee(m11,m31,m21,m41)
647 !                Note : vee(13|24) = vee(23|14) ( so : i    <--> j     )
648 !                       vee(13|24) = vee(14|23) ( so : ip   <--> jp    )
649 !                       vee(13|24) = vee(24|13) ( so : i,ip <--> j,jp  )
650 !                Also : vee(13|24) = vee(31|42) ( so : i,j  <--> ip,jp )
651 !                ==> vee1 is invariant with respect to the permutations i <--> j , ip <--> jp and i,ip <--> j,jp
652 !                ( The term 'phiint_ij * phiint_ipjp' has the same properties)
653                  do sig=1,2
654                    do sigp=1,2
655                      euijkl_temp(sig,sigp) = phiint_ij * phiint_ipjp * vee1
656                    end do
657                  end do
658 
659                  vee2 = pawtab(itypat)%vee(m11,m31,m41,m21)
660 !                Note : vee(13|42) = vee(43|12) ( so : ip   <--> j     )
661 !                       vee(13|42) = vee(12|43) ( so : i    <--> jp    )
662 !                       vee(13|42) = vee(42|13) ( so : i,ip <--> jp,j  )
663 !                Also : vee(13|42) = vee(31|24) ( so : i,j  <--> ip,jp )
664 !                Combining the third and fourth rule we get:
665 !                       vee(13|42) = vee(42|13) = vee(24|31) ( so : i,ip  <--> j,jp )
666 !                ==> vee2 is invariant only with respect to the permutation i,ip <--> j,jp
667 
668 !                Terms i,j,ip,jp (m2,m1,m4,m3) and j,i,jp,ip (m1,m2,m3,m4)
669                  do sig=1,2
670                    euijkl_temp2(sig,sig) = phiint_ij * phiint_ipjp * vee2
671                  end do
672                  pawtab(itypat)%euijkl(:,:,ilmn,jlmn,ilmnp,jlmnp) = euijkl_temp(:,:) - euijkl_temp2(:,:) - euijkl_dc(:,:)
673                  pawtab(itypat)%euijkl(:,:,jlmn,ilmn,jlmnp,ilmnp) = pawtab(itypat)%euijkl(:,:,ilmn,jlmn,ilmnp,jlmnp)
674 
675 !                Term j,i,ip,jp (m1,m2,m4,m3)
676                  vee2 = pawtab(itypat)%vee(m21,m31,m41,m11)
677                  do sig=1,2
678                    euijkl_temp2(sig,sig) = phiint_ij * phiint_ipjp * vee2
679                  end do
680                  pawtab(itypat)%euijkl(:,:,jlmn,ilmn,ilmnp,jlmnp) = euijkl_temp(:,:) - euijkl_temp2(:,:) - euijkl_dc(:,:)
681 
682 !                Term i,j,jp,ip (m2,m1,m3,m4)
683                  vee2 = pawtab(itypat)%vee(m11,m41,m31,m21)
684                  do sig=1,2
685                    euijkl_temp2(sig,sig) = phiint_ij * phiint_ipjp * vee2
686                  end do
687                  pawtab(itypat)%euijkl(:,:,ilmn,jlmn,jlmnp,ilmnp) = euijkl_temp(:,:) - euijkl_temp2(:,:) - euijkl_dc(:,:)
688 
689                end if ! correlated orbitals
690              end do ! klmnb
691            end if ! correlated orbitals
692          end do ! klmna
693 
694        end if ! compute_euijkl
695      end if ! usepawu
696 
697 !    ======================================================================
698 !    D-Local ex-exchange: Matrix elements of coulomb interaction and Fk
699 !    ----------------------------------------------------------------------
700      if (useexexch>0) then
701        lexexch=lcur
702 
703 !      a. compute F(k)
704 !      ---------------------------------------------
705        if (allocated(pawtab(itypat)%fk)) then
706          ABI_DEALLOCATE(pawtab(itypat)%fk)
707        end if
708        ABI_ALLOCATE(pawtab(itypat)%fk,(6,4))
709        pawtab(itypat)%fk=zero
710        ABI_ALLOCATE(ff,(mesh_size))
711        ABI_ALLOCATE(gg,(mesh_size))
712        ff(:)=zero;gg(:)=zero
713        kln=(pawtab(itypat)%lnproju(1)*( pawtab(itypat)%lnproju(1)+1)/2)
714        do ll=1,lexexch+1
715          ll1=2*ll-2
716          if (int_meshsz<mesh_size) ff(int_meshsz+1:mesh_size)=zero
717          ff(1:int_meshsz)=pawtab(itypat)%phiphj(1:int_meshsz,kln)
718          call poisson(ff,ll1,pawrad(itypat),gg)
719          ff(1)=zero
720          ff(2:mesh_size)=(pawtab(itypat)%phiphj(2:mesh_size,kln)*gg(2:mesh_size))&
721 &         /pawrad(itypat)%rad(2:mesh_size)
722          call simp_gen(intg,ff,pawrad(itypat))
723          pawtab(itypat)%fk(1,ll)=intg*(two*ll1+one)
724        end do
725        if (pawtab(itypat)%nproju==2) then
726          kln1=kln+pawtab(itypat)%lnproju(1)
727          kln2=kln1+1
728          do ll=1,lexexch+1
729            ll1=2*ll-2
730            if (int_meshsz<mesh_size) ff(int_meshsz+1:mesh_size)=zero
731            ff(1:int_meshsz)=pawtab(itypat)%phiphj(1:int_meshsz,kln1)
732            call poisson(ff,ll1,pawrad(itypat),gg)
733            ff(1)=zero
734            ff(2:mesh_size)=(pawtab(itypat)%phiphj(2:mesh_size,kln1)*gg(2:mesh_size))&
735 &           /pawrad(itypat)%rad(2:mesh_size)
736            call simp_gen(intg,ff,pawrad(itypat))
737            pawtab(itypat)%fk(2,ll)=intg*(two*ll1+one)
738          end do
739          do ll=1,lexexch+1
740            ll1=2*ll-2
741            if (int_meshsz<mesh_size) ff(int_meshsz+1:mesh_size)=zero
742            ff(1:int_meshsz)=pawtab(itypat)%phiphj(1:int_meshsz,kln2)
743            call poisson(ff,ll1,pawrad(itypat),gg)
744            ff(1)=zero
745            ff(2:mesh_size)=(pawtab(itypat)%phiphj(2:mesh_size,kln2)*gg(2:mesh_size))&
746 &           /pawrad(itypat)%rad(2:mesh_size)
747            call simp_gen(intg,ff,pawrad(itypat))
748            pawtab(itypat)%fk(3,ll)=intg*(two*ll1+one)
749          end do
750          do ll=1,lexexch+1
751            ll1=2*ll-2
752            if (int_meshsz<mesh_size) ff(int_meshsz+1:mesh_size)=zero
753            ff(1:int_meshsz)=pawtab(itypat)%phiphj(1:int_meshsz,kln)
754            call poisson(ff,ll1,pawrad(itypat),gg)
755            ff(1)=zero
756            ff(2:mesh_size)=(pawtab(itypat)%phiphj(2:mesh_size,kln1)*gg(2:mesh_size))&
757 &           /pawrad(itypat)%rad(2:mesh_size)
758            call simp_gen(intg,ff,pawrad(itypat))
759            pawtab(itypat)%fk(4,ll)=intg*(two*ll1+one)
760          end do
761          do ll=1,lexexch+1
762            ll1=2*ll-2
763            if (int_meshsz<mesh_size) ff(int_meshsz+1:mesh_size)=zero
764            ff(1:int_meshsz)=pawtab(itypat)%phiphj(1:int_meshsz,kln)
765            call poisson(ff,ll1,pawrad(itypat),gg)
766            ff(1)=zero
767            ff(2:mesh_size)=(pawtab(itypat)%phiphj(2:mesh_size,kln2)*gg(2:mesh_size))&
768 &           /pawrad(itypat)%rad(2:mesh_size)
769            call simp_gen(intg,ff,pawrad(itypat))
770            pawtab(itypat)%fk(5,ll)=intg*(two*ll1+one)
771          end do
772          do ll=1,lexexch+1
773            ll1=2*ll-2
774            if (int_meshsz<mesh_size) ff(int_meshsz+1:mesh_size)=zero
775            ff(1:int_meshsz)=pawtab(itypat)%phiphj(1:int_meshsz,kln1)
776            call poisson(ff,ll1,pawrad(itypat),gg)
777            ff(1)=zero
778            ff(2:mesh_size)=(pawtab(itypat)%phiphj(2:mesh_size,kln2)*gg(2:mesh_size))&
779 &           /pawrad(itypat)%rad(2:mesh_size)
780            call simp_gen(intg,ff,pawrad(itypat))
781            pawtab(itypat)%fk(6,ll)=intg*(two*ll1+one)
782          end do
783          f4of2=0.6681_dp
784          f6of2=0.4943_dp
785        end if
786        ABI_DEALLOCATE(ff)
787        ABI_DEALLOCATE(gg)
788 
789 !      b. Compute vex.
790 !      ---------------------------------------------
791        if (allocated(pawtab(itypat)%vex)) then
792          ABI_DEALLOCATE(pawtab(itypat)%vex)
793        end if
794        sz1=2*lexexch+1
795        ABI_ALLOCATE(pawtab(itypat)%vex,(sz1,sz1,sz1,sz1,4))
796        pawtab(itypat)%vex=zero
797        lmexexch=(lexexch-1)**2+2*(lexexch-1)+1  ! number of m value below correlated orbitals
798        klm0x=lmexexch*(lmexexch+1)/2            ! value of klmn just below correlated orbitals
799 !      --------- 4 loops for interaction matrix
800        do m1=-lexexch,lexexch
801          m11=m1+lexexch+1
802          do m2=-lexexch,m1
803            m21=m2+lexexch+1
804 !          klma= number of pair before correlated orbitals +
805 !          number of pair for m1 lower than correlated orbitals
806 !          (m1+lexexch+1)*(lexexch-1) + number of pairs for correlated orbitals
807 !          before (m1,m2) + number of pair for m2 lower than current value
808            klma=klm0x+m11*lmexexch+(m11-1)*m11/2+m21
809            do m3=-lexexch,lexexch
810              m31=m3+lexexch+1
811              do m4=-lexexch,m3
812                m41=m4+lexexch+1
813                klmb=klm0x+m31*lmexexch+(m31-1)*m31/2+m41
814 !              --------- loop on k=1,2,3 (4 if f orbitals)
815                do kyc=1,2*lexexch+1,2
816                  lkyc=kyc-1
817                  ll=(kyc+1)/2
818                  lmkyc=(lkyc+1)*(lkyc)+1
819                  ak=zero
820                  do mkyc=-lkyc,lkyc,1
821                    isela=pawang%gntselect(lmkyc+mkyc,klma)
822                    iselb=pawang%gntselect(lmkyc+mkyc,klmb)
823                    if (isela>0.and.iselb>0) ak=ak +pawang%realgnt(isela)*pawang%realgnt(iselb)
824                  end do
825 !                ----- end loop on k=1,2,3 (4 if f orbitals)
826                  pawtab(itypat)%vex(m11,m31,m21,m41,ll)=ak/(two*dble(lkyc)+one)
827                end do  !kyc
828                do ll=1,4
829                  pawtab(itypat)%vex(m11,m31,m21,m41,ll)=pawtab(itypat)%vex(m11,m31,m21,m41,ll)*four_pi
830                  pawtab(itypat)%vex(m21,m31,m11,m41,ll)=pawtab(itypat)%vex(m11,m31,m21,m41,ll)
831                  pawtab(itypat)%vex(m11,m41,m21,m31,ll)=pawtab(itypat)%vex(m11,m31,m21,m41,ll)
832                  pawtab(itypat)%vex(m21,m41,m11,m31,ll)=pawtab(itypat)%vex(m11,m31,m21,m41,ll)
833                end do
834              end do
835            end do
836          end do
837        end do
838 
839      end if !useexexch>0
840 
841      if (present(ucrpa)) then
842        if (ucrpa>=1) then
843          call calc_ubare(itypat,lcur,pawang,pawrad(itypat),pawtab(itypat))
844          call calc_ubare(itypat,lcur,pawang,pawrad(itypat),pawtab(itypat),pawtab(itypat)%rpaw)
845        end if
846      end if
847    end if !lcur/=-1
848  end do !end loop on typat
849 
850  DBG_EXIT("COLL")
851 
852  end subroutine pawpuxinit

m_paw_correlations/pawuenergy [ Functions ]

[ Top ] [ m_paw_correlations ] [ Functions ]

NAME

 pawuenergy

FUNCTION

 Compute contributions to energy for PAW+U calculations

INPUTS

  iatom=index of current atom (absolute index, the index on current proc)
  noccmmp(2*lpawu+1,2*lpawu+1,nspden)=density matrix in the PAW augm. region
  nocctot(nspden)=number of electrons in the correlated subspace
  pawprtvol=control print volume and debugging output for PAW
  pawtab <type(pawtab_type)>=paw tabulated starting data:
     %lpawu=l used for lda+u
     %vee(2*lpawu+1*4)=screened coulomb matrix
  dmft_dc,e_ee,e_dc,e_dcdc,u_dmft,j_dmft= optional arguments for DMFT

OUTPUT

  eldaumdc= PAW+U contribution to total energy
  eldaumdcdc= PAW+U contribution to double-counting total energy

PARENTS

      m_energy,pawdenpot

CHILDREN

      wrtout

SOURCE

 886  subroutine pawuenergy(iatom,eldaumdc,eldaumdcdc,noccmmp,nocctot,pawprtvol,pawtab,&
 887  &                     dmft_dc,e_ee,e_dc,e_dcdc,u_dmft,j_dmft) ! optional arguments (DMFT)
 888 
 889 
 890 !This section has been created automatically by the script Abilint (TD).
 891 !Do not modify the following lines by hand.
 892 #undef ABI_FUNC
 893 #define ABI_FUNC 'pawuenergy'
 894 !End of the abilint section
 895 
 896  implicit none
 897 
 898 !Arguments ---------------------------------------------
 899 !scalars
 900  integer,intent(in) :: iatom,pawprtvol
 901  integer,optional,intent(in) :: dmft_dc
 902  real(dp),intent(in) :: noccmmp(:,:,:,:),nocctot(:)
 903  real(dp),intent(inout) :: eldaumdc,eldaumdcdc
 904  real(dp),optional,intent(inout) :: e_ee,e_dc,e_dcdc
 905  real(dp),optional,intent(in) :: j_dmft,u_dmft
 906  type(pawtab_type),intent(in) :: pawtab
 907 
 908 !Local variables ---------------------------------------
 909 !scalars
 910 !Option for interaction energy in case of non-collinear magnetism:
 911 !           1: E_int=-U/4.N.(N-2)
 912 !           2: E_int=-U/2.(Nup.(Nup-1)+Ndn.(Ndn-1))
 913  integer,parameter :: option_interaction=1
 914  integer :: cplex_occ,dmftdc,ispden,jspden,lpawu,m1,m11,m2,m21,m3,m31,m4,m41,nspden
 915  real(dp) :: eks_opt3,edcdc_opt3,edcdctemp,edctemp,eldautemp,jpawu,jpawu_dc,mnorm,mx,my,mz
 916  real(dp) :: n_sig,n_sigs,n_msig,n_msigs,n_dndn,n_tot,n_upup
 917  real(dp) :: n12_ud_im,n12_du_im
 918  real(dp) :: n12_ud_re,n12_du_re
 919  real(dp) :: n34_ud_im,n34_du_im
 920  real(dp) :: n34_ud_re,n34_du_re
 921  real(dp) :: upawu
 922  real(dp),allocatable :: n12_sig(:),n34_msig(:),n34_sig(:)
 923  character(len=500) :: message
 924 
 925 ! *****************************************************
 926 
 927  nspden=size(nocctot)
 928  cplex_occ=size(noccmmp,1)
 929 
 930  if (size(noccmmp,4)/=nspden) then
 931    message='size of nocctot and noccmmp are inconsistent!'
 932    MSG_BUG(message)
 933  end if
 934  if(present(dmft_dc))  then
 935    dmftdc=dmft_dc
 936    if(pawtab%usepawu<10) then
 937      write(message,'(a,i5)') "usepawu should be =10 if dmft_dc is present ",pawtab%usepawu
 938      MSG_BUG(message)
 939    end if
 940  else
 941    dmftdc=0
 942  end if
 943 
 944  DBG_ENTER("COLL")
 945 
 946  lpawu=pawtab%lpawu
 947  upawu=pawtab%upawu;if(present(u_dmft)) upawu=u_dmft
 948  jpawu=pawtab%jpawu;if(present(j_dmft)) jpawu=j_dmft
 949 
 950 !======================================================
 951 !Compute LDA+U Energy
 952 !-----------------------------------------------------
 953 
 954  eldautemp=zero
 955  edcdc_opt3=zero
 956  eks_opt3=zero
 957 
 958  if (pawtab%usepawu/=5.and.pawtab%usepawu/=6) then
 959 
 960    ABI_ALLOCATE(n12_sig,(cplex_occ))
 961    ABI_ALLOCATE(n34_msig,(cplex_occ))
 962    ABI_ALLOCATE(n34_sig,(cplex_occ))
 963    do ispden=1,min(nspden,2)
 964      jspden=min(nspden,2)-ispden+1
 965 
 966 !    Compute n_sigs and n_msigs for pawtab%usepawu=3
 967      if (nspden<=2) then
 968        n_sig =nocctot(ispden)
 969        n_msig=nocctot(jspden)
 970        n_tot=n_sig+n_msig
 971      else
 972        n_tot=nocctot(1)
 973        mx=nocctot(2)
 974        my=nocctot(3)
 975        mz=nocctot(4)
 976        mnorm=sqrt(mx*mx+my*my+mz*mz)
 977        if (ispden==1) then
 978 !        n_sig =half*(n_tot+mnorm)
 979 !        n_msig=half*(n_tot-mnorm)
 980          n_sig =half*(n_tot+sign(mnorm,mz))
 981          n_msig=half*(n_tot-sign(mnorm,mz))
 982        else
 983 !        n_sig =half*(n_tot-mnorm)
 984 !        n_msig=half*(n_tot+mnorm)
 985          n_sig =half*(n_tot-sign(mnorm,mz))
 986          n_msig=half*(n_tot+sign(mnorm,mz))
 987        end if
 988      end if
 989      n_sigs =n_sig/(float(2*lpawu+1))
 990      n_msigs =n_msig/(float(2*lpawu+1))
 991 !    if(pawtab%usepawu==3) then
 992 !    write(message,fmt=12) "noccmmp11 ",ispden,noccmmp(1,1,1,ispden)
 993 !    call wrtout(std_out,message,'COLL')
 994 !    write(message,fmt=12) "noccmmp11 ",jspden,noccmmp(1,1,1,jspden)
 995 !    call wrtout(std_out,message,'COLL')
 996 !    write(message,fmt=12) "n_sig      ",ispden,n_sig
 997 !    call wrtout(std_out,message,'COLL')
 998 !    write(message,fmt=12) "n_msig     ",jspden,n_msig
 999 !    call wrtout(std_out,message,'COLL')
1000 !    write(message,fmt=12) "n_sigs     ",ispden,n_sigs
1001 !    call wrtout(std_out,message,'COLL')
1002 !    write(message,fmt=12) "n_msigs    ",jspden,n_msigs
1003 !    call wrtout(std_out,message,'COLL')
1004 !    endif
1005 !    12 format(a,i4,e20.10)
1006 
1007 !    Compute interaction energy E_{ee}
1008      do m1=-lpawu,lpawu
1009        m11=m1+lpawu+1
1010        do m2=-lpawu,lpawu
1011          m21=m2+lpawu+1
1012          n12_sig(:)=noccmmp(:,m11,m21,ispden)
1013          if(m21==m11.and.(pawtab%usepawu==3.or.dmftdc==3)) n12_sig(1)=n12_sig(1)-n_sigs
1014          do m3=-lpawu,lpawu
1015            m31=m3+lpawu+1
1016            do m4=-lpawu,lpawu
1017              m41=m4+lpawu+1
1018              n34_sig(:) =noccmmp(:,m31,m41,ispden)
1019              n34_msig(:)=noccmmp(:,m31,m41,jspden)
1020              if(m31==m41.and.(pawtab%usepawu==3.or.dmftdc==3)) then
1021                n34_sig(1)= n34_sig(1) - n_sigs
1022                n34_msig(1)= n34_msig(1) - n_msigs
1023              end if
1024              eldautemp=eldautemp &
1025 &             + n12_sig(1)*n34_msig(1)*pawtab%vee(m11,m31,m21,m41) &
1026 &             + n12_sig(1)*n34_sig(1) *(pawtab%vee(m11,m31,m21,m41)-pawtab%vee(m11,m31,m41,m21))
1027              if(cplex_occ==2) then
1028                eldautemp=eldautemp &
1029 &               - n12_sig(2)*n34_msig(2)*pawtab%vee(m11,m31,m21,m41) &
1030 &               - n12_sig(2)*n34_sig(2) *(pawtab%vee(m11,m31,m21,m41)-pawtab%vee(m11,m31,m41,m21))
1031              end if
1032              if (pawtab%usepawu==3.or.dmftdc==3) then
1033                edcdc_opt3=edcdc_opt3 &
1034 &               + n_sigs*n34_msig(1)*pawtab%vee(m11,m31,m21,m41) &
1035 &               + n_sigs*n34_sig(1) *(pawtab%vee(m11,m31,m21,m41)-pawtab%vee(m11,m31,m41,m21))
1036                eks_opt3=eks_opt3 &
1037 &               + noccmmp(1,m11,m21,ispden)*n34_msig(1)*pawtab%vee(m11,m31,m21,m41) &
1038 &               + noccmmp(1,m11,m21,ispden)*n34_sig(1) *(pawtab%vee(m11,m31,m21,m41)-pawtab%vee(m11,m31,m41,m21))
1039                if(cplex_occ==2) then
1040                  eks_opt3=eks_opt3 &
1041 &                 - noccmmp(2,m11,m21,ispden)*n34_msig(2)*pawtab%vee(m11,m31,m21,m41) &
1042 &                 - noccmmp(2,m11,m21,ispden)*n34_sig(2) *(pawtab%vee(m11,m31,m21,m41)-pawtab%vee(m11,m31,m41,m21))
1043                end if
1044              end if
1045            end do ! m4
1046          end do ! m3
1047        end do ! m2
1048      end do ! m1
1049 
1050    end do ! ispden
1051    if (nspden==1) eldautemp=two*eldautemp ! Non-magn. system: sum up and dn energies
1052    ABI_DEALLOCATE(n12_sig)
1053    ABI_DEALLOCATE(n34_msig)
1054    ABI_DEALLOCATE(n34_sig)
1055 
1056  end if ! usepawu/=5 or 6
1057 
1058 !Non-collinear magnetism: add non-diagonal term; see (Eq 3) in PRB 72, 024458 (2005) [[cite:Shurikov2005]]
1059  if (nspden==4) then
1060    do m1=-lpawu,lpawu
1061      m11=m1+lpawu+1
1062      do m2=-lpawu,lpawu
1063        m21=m2+lpawu+1
1064        n12_ud_re=noccmmp(1,m11,m21,3) ! updn
1065        n12_ud_im=noccmmp(2,m11,m21,3) ! updn
1066        n12_du_re=noccmmp(1,m11,m21,4) ! dnup
1067        n12_du_im=noccmmp(2,m11,m21,4) ! dnup
1068        do m3=-lpawu,lpawu
1069          m31=m3+lpawu+1
1070          do m4=-lpawu,lpawu
1071            m41=m4+lpawu+1
1072            n34_ud_re=noccmmp(1,m31,m41,3)  ! updn
1073            n34_ud_im=noccmmp(2,m31,m41,3)  ! updn
1074            n34_du_re=noccmmp(1,m31,m41,4)  ! dnup
1075            n34_du_im=noccmmp(2,m31,m41,4)  ! dnup
1076            eldautemp=eldautemp-pawtab%vee(m11,m31,m41,m21) &
1077 &           *(n12_ud_re*n34_du_re-n12_ud_im*n34_du_im &
1078 &           +n12_du_re*n34_ud_re-n12_du_im*n34_ud_im)
1079            if (pawtab%usepawu==3.or.dmftdc==3) then
1080              eks_opt3=eks_opt3-pawtab%vee(m11,m31,m41,m21) &
1081 &             *(n12_ud_re*n34_du_re-n12_ud_im*n34_du_im &
1082 &             +n12_du_re*n34_ud_re-n12_du_im*n34_ud_im)
1083            end if
1084          end do ! m4
1085        end do ! m3
1086      end do ! m2
1087    end do ! m1
1088  end if
1089 
1090 !Divide eldautemp by 2; see (Eq 1) in PRB 77, 155104 (2008) [[cite:Amadon2008a]]
1091  eldautemp=half*eldautemp
1092 
1093 !if (nspden==1) then
1094 !n_tot=two*nocctot(1)
1095 !n_upup=nocctot(1)
1096 !n_dndn=nocctot(1)
1097 !else if (nspden==2) then
1098 !n_tot=nocctot(1)+nocctot(2)
1099 !n_upup=nocctot(1)
1100 !n_dndn=nocctot(2)
1101 !else if (nspden==4) then
1102 !n_tot=nocctot(1)
1103 !mx=nocctot(2)
1104 !my=nocctot(3)
1105 !mz=nocctot(4)
1106 !mnorm=sqrt(mx*mx+my*my+mz*mz)
1107 !n_upup=half*(n_tot+mnorm)
1108 !n_dndn=half*(n_tot-mnorm)
1109 !end if
1110  n_upup=n_sig
1111  n_dndn=n_msig
1112 
1113  edcdctemp=zero;edctemp=zero
1114 
1115 !Full localized limit
1116  if((pawtab%usepawu==1.or.pawtab%usepawu==4).or.(dmftdc==1.or.dmftdc==4.or.dmftdc==5)) then
1117    jpawu_dc=jpawu
1118    if(dmftdc==4)  then
1119      jpawu_dc=zero
1120    end if
1121    edcdctemp=edcdctemp-half*upawu*n_tot**2
1122    edctemp  =edctemp  +half*upawu*(n_tot*(n_tot-one))
1123    if (nspden/=4.or.option_interaction==2) then
1124      if(dmftdc/=5.and.pawtab%usepawu/=4) then
1125        edcdctemp=edcdctemp+half*jpawu_dc*(n_upup**2+n_dndn**2)
1126        edctemp  =edctemp  -half*jpawu_dc*(n_upup*(n_upup-one)+n_dndn*(n_dndn-one))
1127      else if(dmftdc==5.or.pawtab%usepawu==4)  then
1128        edcdctemp=edcdctemp+quarter*jpawu_dc*n_tot**2
1129        edctemp  =edctemp  -quarter*jpawu_dc*(n_tot*(n_tot-two))
1130      end if
1131    else if (nspden==4.and.option_interaction==1) then
1132 !    write(message,'(a)') "  warning: option_interaction==1 for test         "
1133 !    call wrtout(std_out,message,'COLL')
1134      edcdctemp=edcdctemp+quarter*jpawu_dc*n_tot**2
1135      edctemp  =edctemp  -quarter*jpawu_dc*(n_tot*(n_tot-two))
1136    else if (nspden==4.and.option_interaction==3) then
1137 !    edcdctemp= \frac{J}/{4}[ N(N) + \vect{m}.\vect{m}]
1138      edcdctemp=edcdctemp+quarter*jpawu_dc*(n_tot**2 + &
1139 &     mx**2+my**2+mz**2)  ! +\frac{J}/{4}\vect{m}.\vect{m}
1140 !    edctemp= -\frac{J}/{4}[ N(N-2) + \vect{m}.\vect{m}]
1141      edctemp  =edctemp  -quarter*jpawu_dc*(  &
1142 &     (n_tot*(n_tot-two)) +   &
1143 &     mx**2+my**2+mz**2)  ! -\frac{J}/{4}\vect{m}.\vect{m}
1144    end if
1145 
1146 !  Around mean field
1147  else if(pawtab%usepawu==2.or.dmftdc==2) then
1148    edctemp=edctemp+upawu*(n_upup*n_dndn)&
1149 &   +half*(upawu-jpawu)*(n_upup**2+n_dndn**2) &
1150 &   *(dble(2*lpawu)/dble(2*lpawu+1))
1151    edcdctemp=-edctemp
1152  else if(pawtab%usepawu==3.or.dmftdc==3) then
1153    edcdctemp=edcdc_opt3
1154    if(abs(pawprtvol)>=3) then
1155      write(message,fmt=11) "edcdc_opt3          ",edcdc_opt3
1156      call wrtout(std_out,message,'COLL')
1157      write(message,fmt=11) "eks_opt3            ",eks_opt3
1158      call wrtout(std_out,message,'COLL')
1159      write(message,fmt=11) "eks+edcdc_opt3      ",eks_opt3+edcdc_opt3
1160      call wrtout(std_out,message,'COLL')
1161      write(message,fmt=11) "(eks+edcdc_opt3)/2  ",(eks_opt3+edcdc_opt3)/2.d0
1162      call wrtout(std_out,message,'COLL')
1163    end if
1164  end if
1165 
1166 
1167  eldaumdc  =eldaumdc  +eldautemp-edctemp
1168  eldaumdcdc=eldaumdcdc-eldautemp-edcdctemp
1169 
1170 !if(pawtab%usepawu/=10.or.pawprtvol>=3) then
1171  if(abs(pawprtvol)>=3) then
1172    if(pawtab%usepawu<10) then
1173      write(message, '(5a,i4)')ch10,'======= LDA+U Energy terms (in Hartree) ====',ch10,&
1174 &     ch10,' For Atom ',iatom
1175    else if (pawtab%usepawu >= 10) then
1176      write(message, '(5a,i4)')ch10,'  ===   LDA+U Energy terms for the DMFT occupation matrix ==',ch10,&
1177 &     ch10,' For Atom ',iatom
1178    end if
1179 
1180    call wrtout(std_out,message,'COLL')
1181    write(message, '(a)' )"   Contributions to the direct expression of energy:"
1182    call wrtout(std_out,  message,'COLL')
1183    write(message,fmt=11) "     Double counting  correction   =",edctemp
1184    call wrtout(std_out,  message,'COLL')
1185    write(message,fmt=11) "     Interaction energy            =",eldautemp
1186    call wrtout(std_out,  message,'COLL')
1187    write(message,fmt=11) "     Total LDA+U Contribution      =",eldautemp-edctemp
1188    call wrtout(std_out,  message,'COLL')
1189    write(message, '(a)' )' '
1190    call wrtout(std_out,  message,'COLL')
1191    write(message, '(a)' )"   For the ""Double-counting"" decomposition:"
1192    call wrtout(std_out,  message,'COLL')
1193    write(message,fmt=11) "     LDA+U Contribution            =",-eldautemp-edcdctemp
1194    call wrtout(std_out,  message,'COLL')
1195    11 format(a,e20.10)
1196    if(abs(pawprtvol)>=2) then
1197      write(message,fmt=11)"     edcdctemp                     =",edcdctemp
1198      call wrtout(std_out,  message,'COLL')
1199      write(message,fmt=11)"     eldaumdcdc for current atom   =",-eldautemp-edcdctemp
1200      call wrtout(std_out,  message,'COLL')
1201      write(message, '(a)' )' '
1202      call wrtout(std_out,  message,'COLL')
1203      write(message,fmt=11)"   pawuenergy: -VUKS pred          =",eldaumdcdc-eldaumdc
1204      call wrtout(std_out,  message,'COLL')
1205    end if
1206    write(message, '(a)' )' '
1207    call wrtout(std_out,  message,'COLL')
1208  end if
1209 
1210 !For DMFT calculation
1211  if(present(e_ee))   e_ee=e_ee+eldautemp
1212  if(present(e_dc))   e_dc=e_dc+edctemp
1213  if(present(e_dcdc)) e_dcdc=e_dcdc+edcdctemp
1214 
1215  DBG_EXIT("COLL")
1216 
1217  end subroutine pawuenergy

m_paw_correlations/pawxenergy [ Functions ]

[ Top ] [ m_paw_correlations ] [ Functions ]

NAME

 pawxenergy

FUNCTION

 Compute contributions to energy for PAW+ local exact exchange calculations

INPUTS

  pawprtvol=control print volume and debugging output for PAW
  pawrhoij <type(pawrhoij_type)>= paw rhoij occupancies and related data
  pawtab <type(pawtab_type)>=paw tabulated starting data:
     %lexexch=l used for local exact-exchange
     %vex(2*lexexch+1*4)=screened coulomb matrix

SIDE EFFECTS

  eexex=energy is updated with the contribution of the cuyrrent atom

PARENTS

      pawdenpot

CHILDREN

      wrtout

SOURCE

1247  subroutine pawxenergy(eexex,pawprtvol,pawrhoij,pawtab)
1248 
1249 
1250 !This section has been created automatically by the script Abilint (TD).
1251 !Do not modify the following lines by hand.
1252 #undef ABI_FUNC
1253 #define ABI_FUNC 'pawxenergy'
1254 !End of the abilint section
1255 
1256  implicit none
1257 
1258 !Arguments ---------------------------------------------
1259 !scalars
1260  integer,intent(in) :: pawprtvol
1261  real(dp),intent(inout) :: eexex
1262  type(pawrhoij_type),intent(in) :: pawrhoij
1263  type(pawtab_type),intent(in) :: pawtab
1264 
1265 !Local variables ---------------------------------------
1266 !scalars
1267  integer :: irhoij,irhoij1,ispden,jrhoij,jrhoij1,klmn,klmn1,lexexch,ll,m11,m21,m31,m41,n1
1268  integer :: n2,n3,n4,nk,nn1,nn2
1269  real(dp) :: eexextemp
1270  character(len=500) :: message
1271 !arrays
1272  integer :: indn(3,3)
1273  real(dp) :: factnk(6)
1274 
1275 ! *****************************************************
1276 
1277  DBG_ENTER("COLL")
1278 
1279  lexexch=pawtab%lexexch
1280  if (pawtab%nproju==1) nk=1
1281  if (pawtab%nproju==2) nk=6
1282  factnk(1)=one;factnk(2)=one;factnk(3)=one
1283  factnk(4)=two;factnk(5)=two;factnk(6)=two
1284  indn(1,1)=1;indn(1,2)=4;indn(1,3)=5
1285  indn(2,1)=4;indn(2,2)=2;indn(2,3)=6
1286  indn(3,1)=5;indn(3,2)=6;indn(3,3)=3
1287 
1288 !======================================================
1289 !Compute local exact exchange Energy
1290 !-----------------------------------------------------
1291  eexextemp=zero
1292 
1293  do ispden=1,pawrhoij%nspden
1294    jrhoij=1
1295    do irhoij=1,pawrhoij%nrhoijsel
1296      klmn=pawrhoij%rhoijselect(irhoij)
1297      if(pawtab%indklmn(3,klmn)==0.and.pawtab%indklmn(4,klmn)==2*lexexch) then
1298        m11=pawtab%klmntomn(1,klmn);m21=pawtab%klmntomn(2,klmn)
1299        n1=pawtab%klmntomn(3,klmn);n2=pawtab%klmntomn(4,klmn)
1300        nn1=(n1*n2)/2+1
1301        jrhoij1=1
1302        do irhoij1=1,pawrhoij%nrhoijsel
1303          klmn1=pawrhoij%rhoijselect(irhoij1)
1304          if(pawtab%indklmn(3,klmn1)==0.and.pawtab%indklmn(4,klmn1)==2*lexexch) then
1305            m31=pawtab%klmntomn(1,klmn1);m41=pawtab%klmntomn(2,klmn1)
1306            n3=pawtab%klmntomn(3,klmn1);n4=pawtab%klmntomn(4,klmn1)
1307            nn2=(n3*n4)/2+1
1308            do ll=1,lexexch+1
1309 !            eexextemp=eexextemp-pawtab%vex(m11,m31,m41,m21,ll)*factnk(indn(nn1,nn2))*pawtab%fk(indn(nn1,nn2),ll)&
1310 !            &        *pawrhoij%rhoijp(jrhoij,ispden)*pawrhoij%rhoijp(jrhoij1,ispden)
1311              eexextemp=eexextemp-pawtab%vex(m11,m31,m41,m21,ll)*pawtab%dltij(klmn)*pawtab%fk(indn(nn1,nn2),ll)&
1312 &             *pawtab%dltij(klmn1)*pawrhoij%rhoijp(jrhoij,ispden)*pawrhoij%rhoijp(jrhoij1,ispden)
1313            end do
1314          end if
1315          jrhoij1=jrhoij1+pawrhoij%cplex
1316        end do !irhoij1
1317      end if
1318      jrhoij=jrhoij+pawrhoij%cplex
1319    end do !irhoij
1320  end do ! ispden
1321  eexextemp=eexextemp/two
1322  eexex=eexex+eexextemp*pawtab%exchmix
1323 
1324  if (abs(pawprtvol)>=2) then
1325    write(message, '(a)' )"   Contributions to the direct expression of energy:"
1326    call wrtout(std_out,message,'COLL')
1327    write(message,fmt='(a,e20.10,a)') "     HF exchange energy  =",eexextemp,ch10
1328    call wrtout(std_out,message,'COLL')
1329  end if
1330 
1331  DBG_EXIT("COLL")
1332 
1333  end subroutine pawxenergy

m_paw_correlations/setnoccmmp [ Functions ]

[ Top ] [ m_paw_correlations ] [ Functions ]

NAME

 setnoccmmp

FUNCTION

 PAW+U only:
   Compute density matrix nocc_{m,m_prime}
   or
   Impose value of density matrix using dmatpawu input array, then symetrize it.

 noccmmp^{\sigma}_{m,m'}=\sum_{ni,nj}[\rho^{\sigma}_{ni,nj}*phiphjint_{ni,nj}]

INPUTS

  compute_dmat= flag: if 1, nocc_{m,mp} is computed
  dimdmat=first dimension of dmatpawu array
  dmatpawu(dimdmat,dimdmat,nsppol*nspinor,natpawu)=input density matrix to be copied into noccmpp
  dmatudiag= flag controlling the use of diagonalization:
             0: no diagonalization of nocc_{m,mp}
             1: diagonalized nocc_{m,mp} matrix is printed
             2: dmatpawu matrix is expressed in the basis where nocc_(m,mp} is diagonal
  impose_dmat= flag: if 1, nocc_{m,mp} is replaced by dmatpawu
  indsym(4,nsym,natom)=indirect indexing array for atom labels
  mpi_atmtab(:)=--optional-- indexes of the atoms treated by current proc
  comm_atom=--optional-- MPI communicator over atoms
  my_natom=number of atoms treated by current processor
  natom=number of atoms in cell
  natpawu=number of atoms on which PAW+U is applied
  nspinor=number of spinorial components of the wavefunctions
  nsppol=number of independant spin components
  nsym=number of symmetry elements in space group
  ntypat=number of atom types
  paw_ij(my_natom) <type(paw_ij_type)>=paw arrays given on (i,j) channels
  pawang <type(pawang_type)>=paw angular mesh and related data
  pawrhoij(my_natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data
  pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data
  spinat(3,matom)=initial spin of each atom, in unit of hbar/2
  symafm(nsym)=(anti)ferromagnetic part of symmetry operations
  typat(natom)=type for each atom
  useexexch=1 if local-exact-exchange is activated
  usepawu=1 if PAW+U is activated

OUTPUT

   paw_ij(natom)%noccmmp(cplex_dij,2*lpawu+1,2*lpawu+1,nsppol or ndij)=density matrix

NOTES

 For non-collinear magnetism,
 - nocc_{m,mp} is computed as:noccmmp(:,:,:,1)=   n{m,mp}
                              noccmmp(:,:,:,2)=   m_x{m,mp}
                              noccmmp(:,:,:,3)=   m_y{m,mp}
                              noccmmp(:,:,:,4)=   m_z{m,mp}
 - but nocc_{m,mp} is stored as: noccmmp(:,:,:,1)=   n^{up,up}_{m,mp}
                                 noccmmp(:,:,:,2)=   n^{dn,dn}_{m,mp}
                                 noccmmp(:,:,:,3)=   n^{up,dn}_{m,mp}
                                 noccmmp(:,:,:,4)=   n^{dn,up}_{m,mp}
   We choose to have noccmmp complex when ndij=4 (ie nspinor=2)
    If ndij=4 and pawspnorb=0, one could keep noccmmp real
    with the n11, n22, Re(n12), Im(n21) representation, but it would
    less clear to change the representation when pawspnorb is activated.
   If ndij=4, nocc_{m,mp} is transformed to the Ylm basis
    and then to the J, M_J basis (if cplex_dij==2)

  Note that n_{m,mp}=<mp|hat(n)|m> because rhoij=<p_j|...|p_i>

PARENTS

      afterscfloop,pawdenpot,pawprt,scfcv,vtorho

CHILDREN

      dgemm,dsyev,free_my_atmtab,get_my_atmtab,mat_mlms2jmj,mat_slm2ylm
      wrtout,zgemm,zheev

SOURCE

1410 subroutine setnoccmmp(compute_dmat,dimdmat,dmatpawu,dmatudiag,impose_dmat,indsym,my_natom,natom,&
1411 &                     natpawu,nspinor,nsppol,nsym,ntypat,paw_ij,pawang,pawprtvol,pawrhoij,pawtab,&
1412 &                     spinat,symafm,typat,useexexch,usepawu, &
1413 &                     mpi_atmtab,comm_atom) ! optional arguments (parallelism)
1414 
1415 
1416 !This section has been created automatically by the script Abilint (TD).
1417 !Do not modify the following lines by hand.
1418 #undef ABI_FUNC
1419 #define ABI_FUNC 'setnoccmmp'
1420 !End of the abilint section
1421 
1422  implicit none
1423 
1424 !Arguments ---------------------------------------------
1425 !scalars
1426  integer,intent(in) :: compute_dmat,dimdmat,dmatudiag,impose_dmat,my_natom,natom,natpawu
1427  integer,intent(in) :: nspinor,nsppol,nsym,ntypat,useexexch,usepawu
1428  integer,optional,intent(in) :: comm_atom
1429  type(pawang_type),intent(in) :: pawang
1430  integer,intent(in) :: pawprtvol
1431 !arrays
1432  integer,intent(in) :: indsym(4,nsym,natom),symafm(nsym),typat(natom)
1433  integer,optional,target,intent(in) :: mpi_atmtab(:)
1434  real(dp),intent(in) :: dmatpawu(dimdmat,dimdmat,nspinor*nsppol,natpawu*impose_dmat)
1435  !real(dp),intent(in) :: dmatpawu(:,:,:,:)
1436  real(dp),intent(in) :: spinat(3,natom)
1437  type(paw_ij_type),intent(inout) :: paw_ij(my_natom)
1438  type(pawrhoij_type),intent(in) :: pawrhoij(my_natom)
1439  type(pawtab_type),intent(in) :: pawtab(ntypat)
1440 
1441 !Local variables ---------------------------------------
1442 !scalars
1443  integer,parameter :: limp=0 ! could become an input variable
1444  integer :: at_indx,cplex_dij,dmatudiag_loc,iafm,iatom,iatom_tot,iatpawu,icount
1445  integer :: ilm,im1,im2,in1,in2,info,iplex,irot,ispden, irhoij,itypat,jlm,jrhoij
1446  integer :: jspden,klmn,kspden,lcur,ldim,lmax,lmin,lpawu,lwork,my_comm_atom,ndij,nmat,nspden,nsploop
1447  logical,parameter :: afm_noncoll=.true.  ! TRUE if antiferro symmetries are used with non-collinear magnetism
1448  logical :: antiferro,my_atmtab_allocated,noccsym_error,paral_atom,use_afm
1449  real(dp),parameter :: invsqrt2=one/sqrt2
1450  real(dp) :: factafm,mnorm,mx,my,mz,ntot,nup,ndn,snorm,sx,sy,szp,szm
1451  character(len=4) :: wrt_mode
1452  character(len=500) :: message
1453 !arrays
1454  integer :: nsym_used(2)
1455  integer,pointer :: my_atmtab(:)
1456  real(dp) :: sumocc(2)
1457  real(dp),allocatable :: eig(:),hdp(:,:,:),hdp2(:,:),noccmmptemp(:,:,:,:),noccmmp_tmp(:,:,:,:)
1458  real(dp),allocatable :: rwork(:),ro(:),noccmmp2(:,:,:,:),nocctot2(:)
1459  complex(dpc),allocatable :: noccmmp_ylm(:,:,:),noccmmp_jmj(:,:),noccmmp_slm(:,:,:)
1460  complex(dpc),allocatable :: zhdp(:,:),zhdp2(:,:),znoccmmp_tmp(:,:),zwork(:)
1461  character(len=9),parameter :: dspin(6)=  (/"up       ","down     ","up-up    ","down-down","Re[up-dn]","Im[up-dn]"/)
1462  character(len=9),parameter :: dspinc(6)= (/"up       ","down     ","up-up    ","down-down","up-dn    ","dn-up    "/)
1463  character(len=9),parameter :: dspinc2(6)=(/"up       ","down     ","dn-dn    ","up-up    ","dn-up    ","up-dn    "/)
1464  character(len=9),parameter :: dspinm(6)= (/"dn       ","up i     ","n        ","mx       ","my       ","mz       "/)
1465  type(coeff4_type),allocatable :: tmp_noccmmp(:)
1466 
1467 !*********************************************************************
1468 
1469  DBG_ENTER("COLL")
1470 
1471 !Tests
1472  if (my_natom>0) then
1473    if (nsppol/=paw_ij(1)%nsppol) then
1474      message=' inconsistent values for nsppol !'
1475      MSG_BUG(message)
1476    end if
1477    if (compute_dmat>0) then
1478      if (pawrhoij(1)%nspden/=paw_ij(1)%nspden.and.&
1479 &     pawrhoij(1)%nspden/=4.and.paw_ij(1)%nspden/=1) then
1480        message=' inconsistent values for nspden !'
1481        MSG_BUG(message)
1482      end if
1483    end if
1484  end if
1485  if (usepawu>0.and.useexexch>0) then
1486    message=' usepawu>0 and useexexch>0 not allowed !'
1487    MSG_BUG(message)
1488  end if
1489  if (impose_dmat/=0.and.dimdmat==0) then
1490    message=' dmatpawu must be allocated when impose_dmat/=0 !'
1491    MSG_BUG(message)
1492  end if
1493  if (usepawu>0.and.compute_dmat/=0.and.impose_dmat/=0.and.pawang%nsym==0) then
1494    message=' pawang%zarot must be allocated !'
1495    MSG_BUG(message)
1496  end if
1497 
1498 !Some inits
1499  if (usepawu==0.and.useexexch==0) return
1500  nspden=1;ndij=1;cplex_dij=1
1501  if (my_natom>0) then
1502    nspden=paw_ij(1)%nspden
1503    ndij=paw_ij(1)%ndij
1504    cplex_dij=paw_ij(1)%cplex_dij
1505  end if
1506  antiferro=(nspden==2.and.nsppol==1)
1507  use_afm=((antiferro).or.((nspden==4).and.afm_noncoll))
1508  dmatudiag_loc=dmatudiag
1509  if (dmatudiag==2.and.(dimdmat==0.or.impose_dmat==0)) dmatudiag_loc=1
1510 
1511 !Set up parallelism over atoms
1512  paral_atom=(present(comm_atom).and.(my_natom/=natom))
1513  nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
1514  my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
1515  call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,my_natom_ref=my_natom) !vz_d
1516  wrt_mode='COLL';if (paral_atom) wrt_mode='PERS'
1517 
1518 !If needed, store dmatpu in suitable format in tmp_noccmmp
1519  if (usepawu>0.and.impose_dmat/=0) then
1520    iatpawu=0
1521    ABI_DATATYPE_ALLOCATE(tmp_noccmmp,(natom))
1522    do iatom_tot=1,natom
1523      itypat=typat(iatom_tot)
1524      lpawu=pawtab(itypat)%lpawu
1525      if (lpawu/=-1) then
1526        iatpawu=iatpawu+1
1527        if (ndij/=4) then
1528          ABI_ALLOCATE(tmp_noccmmp(iatom_tot)%value,(cplex_dij,2*lpawu+1,2*lpawu+1,nsppol))
1529          tmp_noccmmp(iatom_tot)%value(1,1:2*lpawu+1,1:2*lpawu+1,1:nsppol)=&
1530 &         dmatpawu(1:2*lpawu+1,1:2*lpawu+1,1:nsppol,iatpawu)
1531        else
1532          ABI_ALLOCATE(tmp_noccmmp(iatom_tot)%value,(cplex_dij,2*lpawu+1,2*lpawu+1,ndij))
1533          tmp_noccmmp(iatom_tot)%value=zero
1534          if(limp==0) then ! default reading
1535            snorm=sqrt(spinat(1,iatom_tot)**2+spinat(1,iatom_tot)**2+spinat(3,iatom_tot)**2)
1536            if (snorm>tol12) then
1537              sx=half*spinat(1,iatom_tot)/snorm
1538              sy=half*spinat(2,iatom_tot)/snorm
1539              szp=half*(one+spinat(3,iatom_tot)/snorm)
1540              szm=half*(one-spinat(3,iatom_tot)/snorm)
1541            else
1542              sx=zero;sy=zero
1543              szp=one;szm=zero
1544            end if
1545            do im2=1,2*lpawu+1
1546              do im1=1,2*lpawu+1
1547                nup=dmatpawu(im1,im2,1,iatpawu);ndn=dmatpawu(im1,im2,2,iatpawu)
1548                tmp_noccmmp(iatom_tot)%value(1,im1,im2,1)=nup*szp+ndn*szm
1549                tmp_noccmmp(iatom_tot)%value(1,im1,im2,2)=nup*szm+ndn*szp
1550                tmp_noccmmp(iatom_tot)%value(1,im1,im2,3)=(nup-ndn)*sx
1551                tmp_noccmmp(iatom_tot)%value(1,im1,im2,4)=(ndn-nup)*sy
1552              end do
1553            end do
1554          else if(limp>=1) then
1555            ABI_ALLOCATE(noccmmp_ylm,(2*lpawu+1,2*lpawu+1,ndij))
1556            noccmmp_ylm=czero
1557            ABI_ALLOCATE(noccmmp_slm,(2*lpawu+1,2*lpawu+1,ndij))
1558            noccmmp_slm=czero
1559            ABI_ALLOCATE(noccmmp_jmj,(2*(2*lpawu+1),2*(2*lpawu+1)))
1560            noccmmp_jmj=czero
1561            if(limp==1) then ! read input matrix in J,M_J basis (l-1/2, then l+1/2)
1562              noccmmp_jmj=czero
1563              do im1=1,2*lpawu+1
1564                noccmmp_jmj(im1,im1)=cmplx(dmatpawu(im1,im1,1,iatpawu),zero,kind=dp)
1565                noccmmp_jmj(im1+lpawu,im1+lpawu)=cmplx(dmatpawu(im1+lpawu,im1+lpawu,2,iatpawu),zero,kind=dp)
1566              end do
1567              write(message,'(2a,i3,a)') ch10,'== Atom ',iatom_tot,&
1568 &             ' == Imposed occupation matrix (in the J M_J basis: L-1/2 and L+1/2 states)'
1569              call wrtout(std_out,message,wrt_mode)
1570              call mat_mlms2jmj(lpawu,noccmmp_ylm,noccmmp_jmj,ndij,&
1571 &             2,2,pawprtvol,std_out,wrt_mode) !  optspin=1: up spin are first
1572            end if
1573            if(limp==2) then ! read input matrix in Ylm basis
1574              noccmmp_ylm=czero
1575              do im1=1,2*lpawu+1
1576                noccmmp_ylm(im1,im1,1)=cmplx(dmatpawu(im1,im1,1,iatpawu),zero,kind=dp)
1577                noccmmp_ylm(im1,im1,2)=cmplx(dmatpawu(im1,im1,2,iatpawu),zero,kind=dp)
1578              end do
1579              write(message,'(2a,i3,a)') ch10,'== Atom ',iatom_tot,&
1580 &             ' == Imposed occupation matrix (in the Ylm basis), for dn and up spin'
1581              call wrtout(std_out,message,wrt_mode)
1582            end if
1583            call mat_slm2ylm(lpawu,noccmmp_ylm,noccmmp_slm,ndij,&
1584 &           2,2,pawprtvol,std_out,wrt_mode) ! optspin=1 because up spin are first
1585 !          interchange upup and dndn
1586            if(limp>=1) then
1587              tmp_noccmmp(iatom_tot)%value(1,:,:,1)=real(noccmmp_slm(:,:,2))
1588              tmp_noccmmp(iatom_tot)%value(2,:,:,1)=aimag(noccmmp_slm(:,:,2))
1589              tmp_noccmmp(iatom_tot)%value(1,:,:,2)=real(noccmmp_slm(:,:,1))
1590              tmp_noccmmp(iatom_tot)%value(2,:,:,2)=aimag(noccmmp_slm(:,:,1))
1591              tmp_noccmmp(iatom_tot)%value(1,:,:,3)=real(noccmmp_slm(:,:,4))
1592              tmp_noccmmp(iatom_tot)%value(2,:,:,3)=aimag(noccmmp_slm(:,:,4))
1593              tmp_noccmmp(iatom_tot)%value(1,:,:,4)=real(noccmmp_slm(:,:,3))
1594              tmp_noccmmp(iatom_tot)%value(2,:,:,4)=aimag(noccmmp_slm(:,:,3))
1595            end if
1596            if(abs(pawprtvol)>2) then
1597              write(message, '(2a)' ) ch10,&
1598 &             " Check Imposed density matrix in different basis"
1599              call wrtout(std_out,message,wrt_mode)
1600              call mat_slm2ylm(lpawu,noccmmp_slm,noccmmp_ylm,ndij,&
1601 &             1,2,pawprtvol,std_out,wrt_mode) ! optspin=1 because up spin are first
1602              call mat_mlms2jmj(lpawu,noccmmp_ylm,noccmmp_jmj,ndij,1,2,&
1603 &             pawprtvol,std_out,wrt_mode) !  optspin=1: up spin are first
1604            end if
1605            ABI_DEALLOCATE(noccmmp_ylm)
1606            ABI_DEALLOCATE(noccmmp_jmj)
1607            ABI_DEALLOCATE(noccmmp_slm)
1608          end if
1609        end if
1610      end if
1611    end do
1612  end if  ! impose_dmat/=0
1613 
1614 !Print message
1615  if (usepawu>0.and.impose_dmat/=0) then
1616    if (dmatudiag_loc/=2) then
1617      write(message,'(6a)') ch10,'Occupation matrix for correlated orbitals is kept constant',ch10,&
1618 &     'and equal to dmatpawu from input file !',ch10,&
1619 &     '----------------------------------------------------------'
1620    else
1621      write(message,'(6a)') ch10,'Occupation matrix for correlated orbitals is imposed',ch10,&
1622 &     'and equal to dmatpawu in the diagonal basis !',ch10,&
1623 &     '----------------------------------------------------------'
1624    end if
1625    call wrtout(std_out,message,'COLL')
1626  end if
1627 
1628  if (usepawu>0.and.dmatudiag_loc/=0.and.my_natom>0) then
1629    write(message,'(4a)') ch10,'Diagonalized occupation matrix "noccmmp" is printed !',ch10,&
1630 &   '-------------------------------------------------------------'
1631    call wrtout(std_out,message,wrt_mode)
1632  end if
1633 
1634 !Loops over atoms
1635  do iatom=1,my_natom
1636    iatom_tot=iatom;if (paral_atom) iatom_tot=my_atmtab(iatom)
1637    itypat=pawrhoij(iatom)%itypat
1638 
1639    if (useexexch>0) then
1640      lcur=pawtab(itypat)%lexexch
1641    else if (usepawu>0) then
1642      lcur=pawtab(itypat)%lpawu
1643    end if
1644    if (lcur/=-1) then
1645 
1646 !    ########################################################################################
1647 !    # Compute nocc_mmp
1648 !    ########################################################################################
1649      if ((usepawu>0.and.compute_dmat/=0).or.useexexch>0) then
1650 
1651 
1652        paw_ij(iatom)%noccmmp(:,:,:,:)=zero
1653 
1654 !      Loop over spin components
1655        ABI_ALLOCATE(noccmmptemp,(cplex_dij,2*lcur+1,2*lcur+1,ndij))
1656        noccmmptemp(:,:,:,:)=zero
1657        if(ndij==4)  then
1658          ABI_ALLOCATE(noccmmp2,(cplex_dij,2*lcur+1,2*lcur+1,ndij))
1659        end if
1660        if(ndij==4)  then
1661          ABI_ALLOCATE(nocctot2,(ndij))
1662        end if
1663        do ispden=1,ndij
1664          jrhoij=1
1665          do irhoij=1,pawrhoij(iatom)%nrhoijsel
1666            klmn=pawrhoij(iatom)%rhoijselect(irhoij)
1667            im1=pawtab(itypat)%klmntomn(1,klmn)
1668            im2=pawtab(itypat)%klmntomn(2,klmn)
1669            in1=pawtab(itypat)%klmntomn(3,klmn)
1670            in2=pawtab(itypat)%klmntomn(4,klmn)
1671            lmin=pawtab(itypat)%indklmn(3,klmn)
1672            lmax=pawtab(itypat)%indklmn(4,klmn)
1673            ABI_ALLOCATE(ro,(cplex_dij))
1674            if (ndij==1) then
1675              ro(1)=half*pawrhoij(iatom)%rhoijp(jrhoij,1)
1676            else if (ndij==2) then
1677              ro(1)=pawrhoij(iatom)%rhoijp(jrhoij,ispden)
1678            else  ! ndij==4
1679 !            Non-collinear magnetism: transfer rhoij to ro_c (keep n, m storage because
1680 !            it is easier for the computation of noccmmp from rhoij)
1681 !            cplex_dij has to be used here, because it is the dimension of rhoijp
1682              ro(1:cplex_dij)=pawrhoij(iatom)%rhoijp(jrhoij:jrhoij-1+cplex_dij,ispden)
1683            end if
1684            if(lmin==0.and.lmax==2*lcur) then
1685              icount=in1+(in2*(in2-1))/2
1686              if(pawtab(itypat)%ij_proj<icount)  then
1687                message=' PAW+U: Problem in the loop for calculating noccmmp !'
1688                MSG_BUG(message)
1689              end if
1690              if(in1/=in2) then
1691                if(im2<=im1) then
1692                  noccmmptemp(:,im1,im2,ispden)=noccmmptemp(:,im1,im2,ispden)+ro(:)*pawtab(itypat)%phiphjint(icount)
1693                end if
1694              end if
1695              if(im2>=im1) then
1696 
1697                paw_ij(iatom)%noccmmp(:,im1,im2,ispden)=paw_ij(iatom)%noccmmp(:,im1,im2,ispden) &
1698 &               +ro(:)*pawtab(itypat)%phiphjint(icount)
1699              end if
1700            end if
1701            jrhoij=jrhoij+pawrhoij(iatom)%cplex
1702            ABI_DEALLOCATE(ro)
1703          end do ! irhoij
1704          do im2=1,2*lcur+1
1705            do im1=1,im2
1706              paw_ij(iatom)%noccmmp(1,im1,im2,ispden)=paw_ij(iatom)%noccmmp(1,im1,im2,ispden) &
1707 &             +noccmmptemp(1,im2,im1,ispden)
1708              if(cplex_dij==2) paw_ij(iatom)%noccmmp(2,im1,im2,ispden)=paw_ij(iatom)%noccmmp(2,im1,im2,ispden) &
1709 &             -noccmmptemp(2,im2,im1,ispden)
1710            end do
1711          end do
1712          do im1=1,2*lcur+1
1713            do im2=1,im1
1714              paw_ij(iatom)%noccmmp(1,im1,im2,ispden)=paw_ij(iatom)%noccmmp(1,im2,im1,ispden)
1715              if(cplex_dij==2) paw_ij(iatom)%noccmmp(2,im1,im2,ispden)=-paw_ij(iatom)%noccmmp(2,im2,im1,ispden)
1716            end do
1717          end do
1718        end do ! ispden
1719        ABI_DEALLOCATE(noccmmptemp)
1720 !      Compute noccmmp2, occupation matrix in the spin basis (upup, dndn, updn, dnup)
1721        if(ndij==4) then
1722          noccmmp2(:,:,:,:)=zero
1723          do im1=1,2*lcur+1
1724            do im2=1,2*lcur+1
1725              noccmmp2(1,im1,im2,1)=half*(paw_ij(iatom)%noccmmp(1,im1,im2,1)+paw_ij(iatom)%noccmmp(1,im1,im2,4))
1726              noccmmp2(2,im1,im2,1)=half*(paw_ij(iatom)%noccmmp(2,im1,im2,1)+paw_ij(iatom)%noccmmp(2,im1,im2,4))
1727              noccmmp2(1,im1,im2,2)=half*(paw_ij(iatom)%noccmmp(1,im1,im2,1)-paw_ij(iatom)%noccmmp(1,im1,im2,4))
1728              noccmmp2(2,im1,im2,2)=half*(paw_ij(iatom)%noccmmp(2,im1,im2,1)-paw_ij(iatom)%noccmmp(2,im1,im2,4))
1729              noccmmp2(1,im1,im2,3)=half*(paw_ij(iatom)%noccmmp(1,im1,im2,2)+paw_ij(iatom)%noccmmp(2,im1,im2,3))
1730              noccmmp2(2,im1,im2,3)=half*(paw_ij(iatom)%noccmmp(2,im1,im2,2)-paw_ij(iatom)%noccmmp(1,im1,im2,3))
1731              noccmmp2(1,im1,im2,4)=half*(paw_ij(iatom)%noccmmp(1,im1,im2,2)-paw_ij(iatom)%noccmmp(2,im1,im2,3))
1732              noccmmp2(2,im1,im2,4)=half*(paw_ij(iatom)%noccmmp(2,im1,im2,2)+paw_ij(iatom)%noccmmp(1,im1,im2,3))
1733            end do
1734          end do
1735          if(abs(pawprtvol)>=1) then
1736            write(message,'(2a)') ch10,"== Calculated occupation matrix for correlated orbitals in the n, m basis :"
1737            call wrtout(std_out,message,wrt_mode)
1738            do ispden=1,ndij
1739              write(message,'(3a)') ch10,"Calculated occupation matrix for component ",trim(dspinm(ispden+2*(ndij/4)))
1740              call wrtout(std_out,message,wrt_mode)
1741              do im1=1,lcur*2+1  ! ( order of indices in noccmmp is exchanged in order to have the same convention as rhoij: transposition is done after )
1742                if(cplex_dij==1)&
1743 &               write(message,'(12(1x,9(1x,f10.5)))')&
1744 &               (paw_ij(iatom)%noccmmp(1,im2,im1,ispden),im2=1,lcur*2+1)
1745                if(cplex_dij==2)&
1746 !              &               write(message,'(12(1x,9(1x,"(",f7.3,",",f7.3,")")))')&
1747 &               write(message,'(12(1x,9(1x,"(",f10.5,",",f10.5,")")))')&
1748 &               (paw_ij(iatom)%noccmmp(:,im2,im1,ispden),im2=1,lcur*2+1)
1749                call wrtout(std_out,message,wrt_mode)
1750              end do
1751            end do
1752          end if ! pawprtvol >=1
1753        end if
1754 
1755 !      Compute total number of electrons per spin
1756        paw_ij(iatom)%nocctot(:)=zero ! contains nmmp in the n m representation
1757        if(ndij==4) nocctot2(:)=zero ! contains nmmp in the upup dndn updn dnup  representation
1758        do ispden=1,ndij
1759          do im1=1,2*lcur+1
1760            if(ndij==4) then
1761              paw_ij(iatom)%nocctot(ispden)=paw_ij(iatom)%nocctot(ispden)+paw_ij(iatom)%noccmmp(1,im1,im1,ispden)
1762              nocctot2(ispden)=nocctot2(ispden)+noccmmp2(1,im1,im1,ispden)
1763            else
1764              paw_ij(iatom)%nocctot(ispden)=paw_ij(iatom)%nocctot(ispden)+paw_ij(iatom)%noccmmp(1,im1,im1,ispden)
1765            end if
1766          end do
1767        end do
1768 !      noccmmp will now be in the up up , dn dn... representation and now n_mmp=<m|n|mp> instead of <mp|n|m> !
1769        if(ndij==4) then
1770          do ispden=1,ndij
1771            do iplex=1,cplex_dij
1772              do im1=1,2*lcur+1
1773                do im2=1,2*lcur+1
1774                  paw_ij(iatom)%noccmmp(iplex,im1,im2,ispden)=noccmmp2(iplex,im2,im1,ispden) ! now, noccmmp is in the upup dndn updn dnup representation
1775                end do
1776              end do
1777            end do
1778          end do
1779          ABI_DEALLOCATE(noccmmp2)
1780        end if
1781 !      Printing of new nocc_mmp
1782        if ((usepawu>0.and.usepawu<10).or.(usepawu>=10.and.pawprtvol>=3)) &
1783        write(message, '(2a)' )  ch10,'========== LDA+U DATA =================================================== '
1784        if (useexexch>0) write(message, '(2a)' )ch10,'======= Local ex-exchange (PBE0) DATA =================================== '
1785        if(((usepawu>0.and.usepawu<10).or.(usepawu>=10.and.pawprtvol>=3)).or.useexexch>0) call wrtout(std_out,message,wrt_mode)
1786        if (usepawu>=10.and.pawprtvol>=3) then
1787          write(message, '(6a)' )  ch10,'    ( A DFT+DMFT calculation is carried out                              ',&
1788          ch10,'      Thus, the following LDA+U occupation matrices are not physical     ',&
1789          ch10,'      and just informative )'
1790          call wrtout(std_out,message,wrt_mode)
1791        end if
1792        if(usepawu<10.or.pawprtvol>=3) then ! Always write except if DMFT and pawprtvol low
1793          write(message,'(2a,i5,a,i4,a)') ch10,"====== For Atom", iatom_tot,&
1794 &         ", occupations for correlated orbitals. l =",lcur,ch10
1795          call wrtout(std_out,message,wrt_mode)
1796          if(ndij==2) then
1797            do ispden=1,2
1798              write(message,'(a,i4,3a,f10.5)') "Atom", iatom_tot,". Occupations for spin ",&
1799 &             trim(dspin(ispden))," =",paw_ij(iatom)%nocctot(ispden)
1800              call wrtout(std_out,message,wrt_mode)
1801            end do
1802            write(message,'(a,i4,a,2x,e16.8)') "=> On atom",iatom_tot,", local Mag. is  ",&
1803 &           paw_ij(iatom)%nocctot(2)-paw_ij(iatom)%nocctot(1)
1804            call wrtout(std_out,message,wrt_mode)
1805          end if
1806          if(ndij==4) then
1807            ntot=paw_ij(iatom)%nocctot(1)
1808            mx=paw_ij(iatom)%nocctot(2)
1809            my=paw_ij(iatom)%nocctot(3)
1810            mz=paw_ij(iatom)%nocctot(4)
1811            mnorm=sqrt(mx*mx+my*my+mz*mz)
1812            nup=nocctot2(1)
1813            ndn=nocctot2(2)
1814            write(message,'(a,i4,a,2x,e16.8)') "=> On atom",iatom_tot,", local Mag. x is ",mx
1815            call wrtout(std_out,message,wrt_mode)
1816            write(message,'(14x,a,2x,e16.8)') "  local Mag. y is ",my
1817            call wrtout(std_out,message,wrt_mode)
1818            write(message,'(14x,a,2x,e16.8)') "  local Mag. z is ",mz
1819            call wrtout(std_out,message,wrt_mode)
1820            write(message,'(14x,a,2x,e16.8)') "  norm of Mag. is ",mnorm
1821            call wrtout(std_out,message,wrt_mode)
1822            write(message,'(14x,a,2x,f10.5)') "  occ. of majority spin is ",half*(ntot+mnorm)  ! to be checked versus direct calc from noccmmp
1823            call wrtout(std_out,message,wrt_mode)
1824            if(abs(pawprtvol)>=1) write(message,'(14x,a,2x,f10.5)') "  occ. for spin up (along z) ",nup
1825            if(abs(pawprtvol)>=1) then
1826              call wrtout(std_out,message,wrt_mode)
1827            end if
1828            write(message,'(14x,a,2x,f10.5)') "  occ. of minority spin is ",half*(ntot-mnorm)
1829            call wrtout(std_out,message,wrt_mode)
1830            if(abs(pawprtvol)>=1) write(message,'(14x,a,2x,f10.5)') "  occ. for spin dn (along z) ",ndn
1831            if(abs(pawprtvol)>=1) then
1832              call wrtout(std_out,message,wrt_mode)
1833            end if
1834            if(ndij==4)  then
1835              ABI_DEALLOCATE(nocctot2)
1836            end if
1837          end if
1838          write(message,'(2a)') ch10,"== Calculated occupation matrix for correlated orbitals:"
1839          call wrtout(std_out,message,wrt_mode)
1840          do ispden=1,ndij
1841            write(message,'(3a)') ch10,"Calculated occupation matrix for component ",trim(dspinc(ispden+2*(ndij/4)))
1842            call wrtout(std_out,message,wrt_mode)
1843            do im1=1,lcur*2+1
1844              if(cplex_dij==1)&
1845 &             write(message,'(12(1x,9(1x,f10.5)))')&
1846 &             (paw_ij(iatom)%noccmmp(1,im1,im2,ispden),im2=1,lcur*2+1)
1847              if(cplex_dij==2)&
1848 &             write(message,'(12(1x,9(1x,"(",f7.3,",",f7.3,")")))')&
1849 &             (paw_ij(iatom)%noccmmp(:,im1,im2,ispden),im2=1,lcur*2+1)
1850              call wrtout(std_out,message,wrt_mode)
1851            end do
1852          end do
1853        end if
1854 
1855 !      Transformation matrices: real->complex spherical harmonics (for test)
1856        if(ndij==4.and.abs(pawprtvol)>=0) then
1857          ABI_ALLOCATE(noccmmp_ylm,(2*lcur+1,2*lcur+1,ndij))
1858          noccmmp_ylm=czero
1859          ABI_ALLOCATE(noccmmp_slm,(2*lcur+1,2*lcur+1,ndij))
1860          noccmmp_slm=czero
1861          ABI_ALLOCATE(noccmmp_jmj,(2*(2*lcur+1),2*(2*lcur+1)))
1862          noccmmp_jmj=czero
1863 !        go from real notation for complex noccmmp to complex notation in noccmmp_slm
1864          noccmmp_slm(:,:,:)=cmplx(paw_ij(iatom)%noccmmp(1,:,:,:)&
1865 &         ,paw_ij(iatom)%noccmmp(2,:,:,:),kind=dp)
1866          call mat_slm2ylm(lcur,noccmmp_slm,noccmmp_ylm,ndij,1,1,pawprtvol,std_out,wrt_mode) ! optspin=1: up spin are first
1867 
1868          do ispden=1,ndij
1869            write(message,'(3a)') ch10,"Calculated Ylm occupation matrix for component ",trim(dspinc(ispden+2*(ndij/4)))
1870            call wrtout(std_out,message,wrt_mode)
1871            do im1=1,lcur*2+1
1872              write(message,'(12(1x,9(1x,"(",f9.5,",",f9.5,")")))') (noccmmp_ylm(im1,im2,ispden),im2=1,lcur*2+1)
1873              call wrtout(std_out,message,wrt_mode)
1874            end do
1875          end do
1876          call mat_mlms2jmj(lcur,noccmmp_ylm,noccmmp_jmj,ndij,1,1,pawprtvol,std_out,wrt_mode) !  optspin=1: up spin are first
1877          ABI_DEALLOCATE(noccmmp_ylm)
1878          ABI_DEALLOCATE(noccmmp_jmj)
1879          ABI_DEALLOCATE(noccmmp_slm)
1880        end if !ndij==4
1881 
1882      end if ! impose_dmat==0
1883 
1884 !    ########################################################################################
1885 !    # Diagonalize nocc_mmp
1886 !    ########################################################################################
1887      if(usepawu>0.and.dmatudiag_loc>0) then
1888 
1889        lpawu=lcur;ldim=2*lpawu+1
1890        ABI_ALLOCATE(noccmmp_tmp,(1,ldim,ldim,ndij))
1891        if (ndij==4)  then
1892          ABI_ALLOCATE(znoccmmp_tmp,(2*ldim,2*ldim))
1893        end if
1894 
1895 !      Select noccmmp for this atom
1896        do ispden=1,ndij
1897          noccmmp_tmp(1,:,:,ispden)=paw_ij(iatom)%noccmmp(1,:,:,ispden)
1898        end do
1899        if (ndij==4) then
1900          do im2=1,ldim
1901            do im1=1,ldim
1902              znoccmmp_tmp(im1     ,     im2)=cmplx(paw_ij(iatom)%noccmmp(1,im1,im2,1)&
1903 &             ,paw_ij(iatom)%noccmmp(2,im1,im2,1),kind=dp)
1904              znoccmmp_tmp(ldim+im1,ldim+im2)=cmplx(paw_ij(iatom)%noccmmp(1,im1,im2,2)&
1905 &             ,paw_ij(iatom)%noccmmp(2,im1,im2,2),kind=dp)
1906              znoccmmp_tmp(     im1,ldim+im2)=cmplx(paw_ij(iatom)%noccmmp(1,im1,im2,3)&
1907 &             ,paw_ij(iatom)%noccmmp(2,im1,im2,3),kind=dp)
1908              znoccmmp_tmp(ldim+im1,     im2)=cmplx(paw_ij(iatom)%noccmmp(1,im1,im2,4)&
1909 &             ,paw_ij(iatom)%noccmmp(2,im1,im2,4),kind=dp)
1910            end do
1911          end do
1912        end if
1913 
1914 !      Diagonalize nocc_mmp
1915        if (ndij/=4) then
1916          ABI_ALLOCATE(hdp,(ldim,ldim,ndij))
1917          hdp=zero
1918          lwork=3*ldim-1
1919          ABI_ALLOCATE(rwork,(lwork))
1920          ABI_ALLOCATE(eig,(ldim))
1921          do ispden=1,ndij
1922            call dsyev('v','u',ldim,noccmmp_tmp(1,:,:,ispden),ldim,eig,rwork,lwork,info)
1923            if(info/=0) then
1924              message=' Error in diagonalization of noccmmp (DSYEV)!'
1925              MSG_ERROR(message)
1926            end if
1927            do ilm=1,ldim
1928              hdp(ilm,ilm,ispden)=eig(ilm)
1929            end do
1930          end do ! ispden
1931          ABI_DEALLOCATE(rwork)
1932          ABI_DEALLOCATE(eig)
1933        else
1934          ABI_ALLOCATE(hdp,(2*ldim,2*ldim,1))
1935          hdp=zero
1936          lwork=4*ldim-1
1937          ABI_ALLOCATE(rwork,(6*ldim-2))
1938          ABI_ALLOCATE(zwork,(lwork))
1939          ABI_ALLOCATE(eig,(2*ldim))
1940          call zheev('v','u',2*ldim,znoccmmp_tmp,2*ldim,eig,zwork,lwork,rwork,info)
1941          if(info/=0) then
1942            message=' Error in diagonalization of znoccmmp_tmp (zheev) !'
1943            MSG_ERROR(message)
1944          end if
1945          do ilm=1,2*ldim
1946            hdp(ilm,ilm,1)=eig(ilm)
1947          end do
1948          ABI_DEALLOCATE(rwork)
1949          ABI_DEALLOCATE(zwork)
1950          ABI_DEALLOCATE(eig)
1951        end if
1952 
1953 !      Print diagonalized matrix and eigenvectors
1954        do ispden=1,size(hdp,3)
1955          write(message,'(2a,i3,a)') ch10,'== Atom ',iatom_tot,' == Diagonalized Occupation matrix'
1956          if (ndij==1) write(message,fmt='(2a)')     trim(message)," for spin up =="
1957          if (ndij==2) write(message,fmt='(2a,i3,a)')trim(message)," for spin ",ispden," =="
1958          if (ndij==4) write(message,fmt='(2a,i3,a)')trim(message)," =="
1959          call wrtout(std_out,message,wrt_mode)
1960          do ilm=1,size(hdp,1)
1961            write(message,'(12(1x,9(1x,f10.5)))') (hdp(ilm,jlm,ispden),jlm=1,size(hdp,2))
1962            call wrtout(std_out,message,wrt_mode)
1963          end do
1964        end do ! ispden
1965        if(abs(pawprtvol)>=1) then
1966          if (ndij/=4) then
1967            do ispden=1,ndij
1968              write(message,'(2a,i3,a)') ch10,'== Atom ',iatom_tot,' == Eigenvectors'
1969              if (ndij==1) write(message,fmt='(2a)')     trim(message),' for spin up =='
1970              if (ndij==2) write(message,fmt='(2a,i3,a)')trim(message),' for spin ',ispden,' =='
1971              call wrtout(std_out,message,wrt_mode)
1972              do ilm=1,ldim
1973                write(message,'(12(1x,9(1x,f10.5)))') (noccmmp_tmp(1,ilm,jlm,ispden),jlm=1,ldim)
1974                call wrtout(std_out,message,wrt_mode)
1975              end do
1976            end do
1977          else
1978            write(message,'(2a,i3,a)') ch10,'== Atom ',iatom_tot,' == Eigenvectors (spinors) in the real harmonics basis =='
1979            call wrtout(std_out,message,wrt_mode)
1980            do ilm=1,2*ldim
1981              write(message,'(12(1x,9(1x,"(",f7.3,",",f7.3,")")))') (znoccmmp_tmp(ilm,jlm),jlm=1,2*ldim)
1982              call wrtout(std_out,message,wrt_mode)
1983            end do
1984          end if
1985        end if
1986 
1987 !      Back rotation of diagonalized matrix and printing
1988        if(abs(pawprtvol)>=1) then
1989          if (ndij/=4) then
1990            ABI_ALLOCATE(hdp2,(ldim,ldim))
1991            do ispden=1,ndij
1992              call dgemm('n','t',ldim,ldim,ldim,one,hdp(:,:,ispden),ldim,noccmmp_tmp(1,:,:,ispden),ldim,zero,hdp2,ldim)
1993              call dgemm('n','n',ldim,ldim,ldim,one,noccmmp_tmp(1,:,:,ispden),ldim,hdp2,ldim,zero,hdp(:,:,ispden),ldim)
1994              noccmmp_tmp(1,:,:,ispden)=hdp(:,:,ispden)
1995            end do ! ispden
1996            ABI_DEALLOCATE(hdp2)
1997          else
1998            ABI_ALLOCATE(zhdp,(2*ldim,2*ldim))
1999            ABI_ALLOCATE(zhdp2,(2*ldim,2*ldim))
2000            zhdp(:,:)=cmplx(hdp(:,:,1),zero,kind=dp)
2001            zhdp2(:,:)=cmplx(zero,zero,kind=dp)
2002            call zgemm('n','c',2*ldim,2*ldim,2*ldim,cone,zhdp,2*ldim,znoccmmp_tmp,2*ldim,czero,zhdp2,2*ldim)
2003            zhdp(:,:)=cmplx(zero,zero,kind=dp)
2004            call zgemm('n','n',2*ldim,2*ldim,2*ldim,cone,znoccmmp_tmp,2*ldim,zhdp2,2*ldim,czero,zhdp,2*ldim)
2005            znoccmmp_tmp=zhdp
2006            ABI_DEALLOCATE(zhdp)
2007            ABI_DEALLOCATE(zhdp2)
2008          end if
2009          nmat=ndij ; if(ndij==4.and.cplex_dij==2) nmat=1
2010          do ispden=1,nmat
2011            write(message,'(2a,i3,a)') ch10,'== Atom ',iatom_tot,&
2012 &           ' == Rotated back diagonalized matrix'
2013            if (ndij==1) write(message,fmt='(2a)')     trim(message)," for spin up =="
2014            if (ndij==2) write(message,fmt='(2a,i3,a)')trim(message)," for spin ",ispden," =="
2015            if (ndij==4.and.cplex_dij==2) write(message,fmt='(4a)')     trim(message)," for all component "
2016            call wrtout(std_out,message,wrt_mode)
2017            do ilm=1,ldim*cplex_dij
2018              if(ndij==1.or.ndij==2)&
2019 &             write(message,'(12(1x,9(1x,f10.5)))')&
2020 &             (noccmmp_tmp(1,ilm,jlm,ispden),jlm=1,ldim)
2021              if(ndij==4.and.cplex_dij==2)&
2022 &             write(message,'(12(1x,18(1x,"(",f7.3,",",f7.3,")")))')&
2023 &             (znoccmmp_tmp(ilm,jlm),jlm=1,ldim*cplex_dij)
2024              call wrtout(std_out,message,wrt_mode)
2025            end do
2026          end do ! ispden
2027        end if
2028        ABI_DEALLOCATE(hdp)
2029 
2030      end if ! dmatudiag_loc
2031 
2032 !    ########################################################################################
2033 !    # Impose value of nocc_mmp from dmatpu; symetrize it
2034 !    ########################################################################################
2035      if (usepawu>0.and.impose_dmat/=0) then
2036 
2037        lpawu=lcur
2038        nsploop=nsppol;if (ndij==4) nsploop=4
2039        noccsym_error=.false.
2040 
2041 !      Loop over spin components
2042        do ispden=1,nsploop
2043          if (ndij/=4) then
2044            jspden=min(3-ispden,paw_ij(iatom)%nsppol)
2045          else if (ispden<=2) then
2046            jspden=3-ispden
2047          else
2048            jspden=ispden
2049          end if
2050 
2051 !        Loops over components of nocc_mmp
2052          do jlm=1,2*lpawu+1
2053            do ilm=1,2*lpawu+1
2054 
2055              if(nsym>1.and.ndij<4) then
2056 
2057                nsym_used(1:2)=0
2058                sumocc(1:2)=zero
2059 
2060 !              Accumulate values of nocc_mmp over symmetries
2061                do irot=1,nsym
2062                  if ((symafm(irot)/=1).and.(.not.use_afm)) cycle
2063                  kspden=ispden;if (symafm(irot)==-1) kspden=jspden
2064                  factafm=one;if (ispden>3) factafm=dble(symafm(irot))
2065                  iafm=1;if ((antiferro).and.(symafm(irot)==-1)) iafm=2
2066                  nsym_used(iafm)=nsym_used(iafm)+1
2067                  at_indx=indsym(4,irot,iatom_tot)
2068                  do im2=1,2*lpawu+1
2069                    do im1=1,2*lpawu+1
2070 !                    Be careful: use here R_rel^-1 in term of spherical harmonics
2071 !                    which is tR_rec in term of spherical harmonics
2072 !                    so, use transpose[zarot]
2073                      sumocc(iafm)=sumocc(iafm)+factafm*tmp_noccmmp(at_indx)%value(1,im1,im2,kspden) &
2074 &                     *pawang%zarot(im1,ilm,lpawu+1,irot)&
2075 &                     *pawang%zarot(im2,jlm,lpawu+1,irot)
2076 !                    sumocc(iafm)=sumocc(iafm)+factafm*tmp_noccmmp(at_indx)%value(im1,im2,kspden) &
2077 !                    &                     *pawang%zarot(ilm,im1,lpawu+1,irot)&
2078 !                    &                     *pawang%zarot(jlm,im2,lpawu+1,irot)
2079                    end do
2080                  end do
2081                end do ! End loop over symmetries
2082 
2083 !              Store new values of nocc_mmp
2084                paw_ij(iatom)%noccmmp(1,ilm,jlm,ispden)=sumocc(1)/nsym_used(1)
2085                if (.not.noccsym_error)&
2086 &               noccsym_error=(abs(paw_ij(iatom)%noccmmp(1,ilm,jlm,ispden) &
2087 &               -tmp_noccmmp(iatom_tot)%value(1,ilm,jlm,ispden))>tol5)
2088 
2089 !              Antiferromagnetic case: has to fill up "down" component of nocc_mmp
2090                if (antiferro.and.nsym_used(2)>0) paw_ij(iatom)%noccmmp(1,ilm,jlm,2)=sumocc(2)/nsym_used(2)
2091 
2092              else  ! nsym=1
2093 
2094 !              Case without symetries
2095                paw_ij(iatom)%noccmmp(:,ilm,jlm,ispden)= tmp_noccmmp(iatom_tot)%value(:,ilm,jlm,ispden)
2096              end if
2097 
2098            end do !ilm
2099          end do !jlm
2100        end do ! ispden
2101        do ispden=1,nsploop
2102          paw_ij(iatom)%nocctot(ispden)=zero ! contains nmmp in the n m representation
2103          do im1=1,2*lcur+1
2104            if(ndij==4.and.ispden==1) then
2105 !            in this case, on computes total number or electron for double counting correction
2106              paw_ij(iatom)%nocctot(ispden)=paw_ij(iatom)%nocctot(ispden)+&
2107 &             paw_ij(iatom)%noccmmp(1,im1,im1,1)+paw_ij(iatom)%noccmmp(1,im1,im1,2)
2108            else if(ndij==4.and.ispden==2) then
2109              paw_ij(iatom)%nocctot(ispden)=paw_ij(iatom)%nocctot(ispden)+&
2110 &             paw_ij(iatom)%noccmmp(1,im1,im1,3)+paw_ij(iatom)%noccmmp(1,im1,im1,4)
2111            else if(ndij==4.and.ispden==3) then
2112              paw_ij(iatom)%nocctot(ispden)=paw_ij(iatom)%nocctot(ispden)-&
2113 &             paw_ij(iatom)%noccmmp(2,im1,im1,3)+paw_ij(iatom)%noccmmp(2,im1,im1,4)
2114            else if(ndij==4.and.ispden==4) then
2115              paw_ij(iatom)%nocctot(ispden)=paw_ij(iatom)%nocctot(ispden)+&
2116 &             paw_ij(iatom)%noccmmp(2,im1,im1,1)-paw_ij(iatom)%noccmmp(2,im1,im1,2)
2117            else
2118              paw_ij(iatom)%nocctot(ispden)=paw_ij(iatom)%nocctot(ispden)+&
2119 &             paw_ij(iatom)%noccmmp(1,im1,im1,ispden)
2120            end if
2121          end do
2122        end do ! ispden
2123 
2124 !      Printing of new nocc_mmp
2125        do ispden=1,ndij
2126          if(dmatudiag_loc==2) then
2127            write(message,'(2a,i3,a)') ch10,'== Atom ',iatom_tot,&
2128 &           ' == Imposed occupation matrix (in the basis of diagonalization!!)'
2129          else
2130            write(message,'(2a,i3,a)') ch10,'== Atom ',iatom_tot,&
2131 &           ' == Imposed occupation matrix'
2132          end if
2133          if (ndij==1) write(message,fmt='(2a)')     trim(message)," for spin up =="
2134          if (ndij==2) write(message,fmt='(2a,i3,a)')trim(message)," for spin ",ispden," =="
2135          if (ndij==4) write(message,fmt='(4a)')     trim(message)," for component ", &
2136 &         trim(dspinc(ispden+2*(ndij/4)))," =="
2137          call wrtout(std_out,message,wrt_mode)
2138          do ilm=1,2*lpawu+1
2139            if(cplex_dij==1)&
2140 &           write(message,'(12(1x,9(1x,f10.5)))')&
2141 &           (paw_ij(iatom)%noccmmp(1,ilm,jlm,ispden),jlm=1,2*lpawu+1)
2142            if(cplex_dij==2)&
2143 &           write(message,'(12(1x,9(1x,"(",f7.3,",",f7.3,")")))')&
2144 &           (paw_ij(iatom)%noccmmp(:,ilm,jlm,ispden),jlm=1,2*lpawu+1)
2145            call wrtout(std_out,message,wrt_mode)
2146          end do
2147        end do
2148 
2149 !      WARNING if symmetrization changes the matrix
2150        if (noccsym_error) then
2151          write(message, '(a,i4,6a)' ) &
2152          '   After symmetrization, imposed occupation matrix for atom ',iatom_tot,ch10,&
2153 &         '   is different from dmatpawu value set in input file !',ch10,&
2154 &         '   It is likely that dmatpawu does not match the symmetry operations of the system.',ch10,&
2155 &         '   Action: change dmatpawu in input file or increase precision until 0.00001'
2156          MSG_WARNING(message)
2157        end if
2158 
2159      end if ! impose_dmat/=0
2160 
2161 !    ########################################################################################
2162 !    # Rotate imposed occupation matrix in the non-diagonal basis
2163 !    ########################################################################################
2164      if (usepawu>0.and.impose_dmat/=0.and.dmatudiag_loc==2) then
2165 
2166        lpawu=lcur;ldim=2*lpawu+1
2167 
2168 !      Rotation of imposed nocc_mmp
2169        if (ndij/=4) then
2170          ABI_ALLOCATE(hdp2,(ldim,ldim))
2171          do ispden=1,ndij
2172            call dgemm('n','t',ldim,ldim,ldim,one,&
2173 &           paw_ij(iatom)%noccmmp(1,:,:,ispden),ldim,noccmmp_tmp(1,:,:,ispden),ldim,zero,hdp2,ldim)
2174            call dgemm('n','n',ldim,ldim,ldim,one,&
2175 &           noccmmp_tmp(1,:,:,ispden),ldim,hdp2,ldim,zero,paw_ij(iatom)%noccmmp(1,:,:,ispden),ldim)
2176          end do ! ispden
2177          ABI_DEALLOCATE(hdp2)
2178        else
2179          ABI_ALLOCATE(zhdp,(2*ldim,2*ldim))
2180          ABI_ALLOCATE(zhdp2,(2*ldim,2*ldim))
2181          do im2=1,ldim
2182            do im1=1,ldim
2183              zhdp(     im1,     im2)=cmplx(paw_ij(iatom)%noccmmp(1,im1,im2,1),zero,kind=dp)  ! to be checked
2184              zhdp(ldim+im1,ldim+im2)=cmplx(paw_ij(iatom)%noccmmp(1,im1,im2,2),zero,kind=dp)  ! to be checked
2185              zhdp(     im1,ldim+im2)=&
2186 &             cmplx(paw_ij(iatom)%noccmmp(1,im1,im2,3),+paw_ij(iatom)%noccmmp(1,im1,im2,4),kind=dp)  ! to be checked
2187              zhdp(ldim+im1,     im2)=&
2188 &             cmplx(paw_ij(iatom)%noccmmp(1,im2,im1,3),-paw_ij(iatom)%noccmmp(1,im2,im1,4),kind=dp)  ! to be checked
2189            end do
2190          end do
2191          call zgemm('n','c',2*ldim,2*ldim,2*ldim,cone,zhdp,2*ldim,znoccmmp_tmp,2*ldim,czero,zhdp2,2*ldim)
2192          call zgemm('n','n',2*ldim,2*ldim,2*ldim,cone,znoccmmp_tmp,2*ldim,zhdp2,2*ldim,czero,zhdp,2*ldim)
2193          do jlm=1,ldim
2194            do ilm=1,ldim
2195              paw_ij(iatom)%noccmmp(1,ilm,jlm,1)= real(znoccmmp_tmp(     ilm,     jlm))  ! to be checked
2196              paw_ij(iatom)%noccmmp(1,ilm,jlm,2)= real(znoccmmp_tmp(ldim+ilm,ldim+jlm))  ! to be checked
2197              paw_ij(iatom)%noccmmp(1,ilm,jlm,3)= real(znoccmmp_tmp(     ilm,ldim+jlm))  ! to be checked
2198              paw_ij(iatom)%noccmmp(1,ilm,jlm,4)=aimag(znoccmmp_tmp(     ilm,ldim+jlm))  ! to be checked
2199            end do
2200          end do
2201          ABI_DEALLOCATE(zhdp)
2202          ABI_DEALLOCATE(zhdp2)
2203        end if
2204 
2205 !      Printing of rotated imposed matrix
2206        do ispden=1,ndij
2207          write(message,'(2a,i3,a)') ch10,'== Atom ',iatom_tot,&
2208 &         ' == Imposed density matrix in original basis'
2209          if (ndij==1) write(message,fmt='(2a)')     trim(message)," for spin up =="
2210          if (ndij==2) write(message,fmt='(2a,i3,a)')trim(message)," for spin ",ispden," =="
2211          if (ndij==4) write(message,fmt='(4a)')     trim(message)," for component ", &
2212 &         trim(dspin(ispden+2*(ndij/4)))," =="
2213          call wrtout(std_out,message,wrt_mode)
2214          do ilm=1,2*lpawu+1
2215            write(message,'(12(1x,9(1x,f10.5)))') (paw_ij(iatom)%noccmmp(1,ilm,jlm,ispden),jlm=1,2*lpawu+1)  ! to be checked
2216            call wrtout(std_out,message,wrt_mode)
2217          end do
2218        end do ! ispden
2219 
2220      end if ! dmatudiag_loc==2
2221 
2222      if (usepawu>0.and.dmatudiag_loc>0) then
2223        ABI_DEALLOCATE(noccmmp_tmp)
2224        if (ndij==4)  then
2225          ABI_DEALLOCATE(znoccmmp_tmp)
2226        end if
2227      end if
2228 
2229      paw_ij(iatom)%has_pawu_occ=2
2230 
2231    end if ! lcur
2232  end do ! iatom
2233 
2234 !Memory deallocation
2235  if (usepawu>0.and.impose_dmat/=0) then
2236    do iatom_tot=1,natom
2237      lpawu=pawtab(typat(iatom_tot))%lpawu
2238      if (lpawu/=-1)  then
2239        ABI_DEALLOCATE(tmp_noccmmp(iatom_tot)%value)
2240      end if
2241    end do
2242    ABI_DATATYPE_DEALLOCATE(tmp_noccmmp)
2243  end if
2244 
2245 !Destroy atom table used for parallelism
2246  call free_my_atmtab(my_atmtab,my_atmtab_allocated)
2247 
2248  DBG_EXIT("COLL")
2249 
2250 end subroutine setnoccmmp

m_paw_correlations/setrhoijpbe0 [ Functions ]

[ Top ] [ m_paw_correlations ] [ Functions ]

NAME

 setrhoijpbe0

FUNCTION

 PAW local exact exchange only:
 Impose value of rhoij for f electrons using an auxiliairy file

INPUTS

  dtset <type(dataset_type)>=all input variables for this dataset
  initialized= if 0, the initialization of the gstate run is not yet finished
  istep=index of the number of steps in the routine scfcv
  istep_mix=index of the number of steps for the SCF mixing (can be <istep)
  mpi_atmtab(:)=--optional-- indexes of the atoms treated by current proc
  comm_atom=--optional-- MPI communicator over atoms
  mpi_comm_read=MPI communicator containing all the processes reading the PBE0 file
  my_natom=number of atoms treated by current processor
  natom=number of atoms in cell
  ntypat=number of types of atoms in unit cell
  pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
  typat(natom)=type integer for each atom in cell

SIDE EFFECTS

  istep_mix=index of the number of steps for the SCF mixing (can be <istep)
  pawrhoij(natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data

NOTES

  Only valid for f electrons !!!

PARENTS

      scfcv

CHILDREN

      free_my_atmtab,get_my_atmtab,pawio_print_ij,wrtout,xmpi_sum

SOURCE

2294 subroutine setrhoijpbe0(dtset,initialized,istep,istep_mix,&
2295 &                       mpi_comm_read,my_natom,natom,ntypat,pawrhoij,pawtab,typat,&
2296 &                       mpi_atmtab,comm_atom) ! optional arguments (parallelism)
2297 
2298 
2299 !This section has been created automatically by the script Abilint (TD).
2300 !Do not modify the following lines by hand.
2301 #undef ABI_FUNC
2302 #define ABI_FUNC 'setrhoijpbe0'
2303 !End of the abilint section
2304 
2305  implicit none
2306 
2307 !Arguments ---------------------------------------------
2308 !scalars
2309  integer,intent(in) :: initialized,istep,mpi_comm_read,my_natom,natom,ntypat
2310  integer,intent(inout) :: istep_mix
2311  integer,optional,intent(in) :: comm_atom
2312  type(dataset_type),intent(in) :: dtset
2313 !arrays
2314  integer,intent(in) :: typat(natom)
2315  integer,optional,target,intent(in) :: mpi_atmtab(:)
2316  type(pawrhoij_type),intent(inout) :: pawrhoij(my_natom)
2317  type(pawtab_type),intent(in) :: pawtab(ntypat)
2318 
2319 !Local variables ---------------------------------------
2320 !scalars
2321  integer,parameter :: ll=3
2322  integer :: iatom,iatom_tot,ierr,ii,ios,iread,irhoij,ispden,itypat,jj,klmn,my_comm_atom,my_rank,nselect
2323  integer :: nstep1,nstep1_abs,rhoijshft,rhoijsz
2324  logical :: my_atmtab_allocated,paral_atom,test0
2325  character(len=9),parameter :: filnam='rhoijpbe0'
2326  character(len=9),parameter :: dspin(6)=(/"up       ","down     ","up-up    ","down-down","Re[up-dn]","Im[up-dn]"/)
2327  character(len=500) :: strg, message
2328 !arrays
2329  integer, allocatable :: nspden_tmp(:)
2330  integer,pointer :: my_atmtab(:)
2331  real(dp),allocatable :: rhoijtmp(:,:),rhoijtmp1(:,:),rhoijtmp2(:,:,:,:)
2332 
2333 ! *********************************************************************
2334 
2335  DBG_ENTER("COLL")
2336 
2337 !Test existence of file and open it
2338  inquire(file=filnam,iostat=ios,exist=test0)
2339  if(.not.test0) return
2340 
2341 !Look for parallelisation over atomic sites
2342  paral_atom=(present(comm_atom).and.(my_natom/=natom))
2343 
2344 !Test if exact-exch. is on f electrons
2345  test0=.false.
2346  do itypat=1,ntypat
2347    if (pawtab(itypat)%useexexch>0.and.pawtab(itypat)%lexexch/=ll) test0=.true.
2348  end do
2349  if (test0) then
2350    write(message, '(3a,i1,a)' ) &
2351 &   ' Local exact exchange: occ. matrix can only be imposed for l=',ll,' !'
2352    MSG_ERROR(message)
2353  end if
2354 
2355 !============================================================
2356 !===== First case: no parallelisation over atomic sites =====
2357 !============================================================
2358 
2359  if (.not.paral_atom) then
2360 
2361 !  Open file
2362    if (open_file(filnam,message,unit=77,form='formatted') /= 0) then
2363      MSG_ERROR(message)
2364    end if
2365 
2366 !  Read step number and eventually exit
2367    nstep1=0;test0=.false.
2368    do while (.not.test0)
2369      read(77,'(A)') strg
2370      test0=(strg(1:1)/="#")
2371      if (test0) read(unit=strg,fmt=*) nstep1
2372    end do
2373    nstep1_abs=abs(nstep1)
2374    if (nstep1_abs==0.or.istep>nstep1_abs.or.(nstep1>0.and.initialized/=0)) then
2375      close(77)
2376 !    Reinitalize mixing when rhoij is allowed to change; for experimental purpose...
2377      if (dtset%userib==1234.and.istep==1+nstep1_abs.and.(nstep1<0.or.initialized==0)) istep_mix=1
2378      return
2379    end if
2380 
2381 !  Loop on atoms
2382    do iatom=1,natom
2383      itypat=typat(iatom)
2384      if (pawtab(itypat)%useexexch>0) then
2385 
2386 !      Set sizes depending on ll
2387        rhoijsz=4*ll+2
2388        rhoijshft=2*ll*ll
2389 
2390 !      Uncompress rhoij
2391        ABI_ALLOCATE(rhoijtmp,(pawrhoij(iatom)%lmn2_size,pawrhoij(iatom)%nspden))
2392        do ispden=1,pawrhoij(iatom)%nspden
2393          rhoijtmp=zero
2394          do irhoij=1,pawrhoij(iatom)%nrhoijsel
2395            klmn=pawrhoij(iatom)%rhoijselect(irhoij)
2396            rhoijtmp(klmn,ispden)=pawrhoij(iatom)%rhoijp(irhoij,ispden)
2397          end do
2398        end do
2399 !      Read rhoij from file
2400        ABI_ALLOCATE(rhoijtmp1,(rhoijsz,rhoijsz))
2401        do ispden=1,pawrhoij(iatom)%nspden
2402          do ii=1,rhoijsz
2403            test0=.false.
2404            do while (.not.test0)
2405              read(77,'(A)') strg
2406              test0=(strg(1:1)/="#")
2407              if (test0)  read(unit=strg,fmt=*) (rhoijtmp1(ii,jj), jj=1,rhoijsz)
2408            end do
2409          end do
2410 
2411 !        Impose rhoij
2412          do jj=1,rhoijsz
2413            do ii=1,jj
2414              rhoijtmp((jj+rhoijshft)*((jj+rhoijshft)-1)/2+ii+rhoijshft,ispden)=rhoijtmp1(ii,jj)
2415            end do
2416          end do
2417 
2418        end do
2419        ABI_DEALLOCATE(rhoijtmp1)
2420 
2421 !      Compress rhoij
2422        nselect=0
2423        do klmn=1,pawrhoij(iatom)%lmn2_size
2424          if (any(abs(rhoijtmp(klmn,:))>tol10)) then
2425            nselect=nselect+1
2426            do ispden=1,pawrhoij(iatom)%nspden
2427              pawrhoij(iatom)%rhoijp(nselect,ispden)=rhoijtmp(klmn,ispden)
2428            end do
2429            pawrhoij(iatom)%rhoijselect(nselect)=klmn
2430          end if
2431        end do
2432        pawrhoij(iatom)%nrhoijsel=nselect
2433        ABI_DEALLOCATE(rhoijtmp)
2434 
2435 !      Print new rhoij
2436        do ispden=1,pawrhoij(iatom)%nspden
2437          write(message,'(2a,i3,a)') ch10,'== Atom ',iatom,&
2438 &         ' == Imposed occupation matrix'
2439          if (pawrhoij(iatom)%nspden==1) write(message,fmt='(2a)')     trim(message)," for spin up =="
2440          if (pawrhoij(iatom)%nspden==2) write(message,fmt='(2a,i3,a)')trim(message)," for spin ",ispden," =="
2441          if (pawrhoij(iatom)%nspden==4) write(message,fmt='(4a)')     trim(message)," for component ", &
2442 &         trim(dspin(ispden+2*(pawrhoij(iatom)%nspden/4)))," =="
2443          call wrtout(std_out,message,'COLL')
2444          call pawio_print_ij(std_out,pawrhoij(iatom)%rhoijp(:,ispden),pawrhoij(iatom)%nrhoijsel,&
2445 &         pawrhoij(iatom)%cplex,pawrhoij(iatom)%lmn_size,ll,&
2446 &         pawtab(itypat)%indlmn(1,1:pawtab(itypat)%lmn_size),&
2447 &         1,-1,pawrhoij(iatom)%rhoijselect(:),-1.d0,1,mode_paral='COLL')
2448        end do
2449 
2450 !      End loop on atoms
2451      end if
2452    end do
2453 
2454 !  Close file
2455    close (77)
2456 
2457  else
2458 
2459 !  ============================================================
2460 !  ====== 2nd case: no parallelisation over atomic sites =====
2461 !  ============================================================
2462 
2463    my_rank=xmpi_comm_rank(mpi_comm_read)
2464 
2465 !  Read step number and eventually exit
2466    iread=0
2467    if (my_rank==0) then
2468      if (open_file(filnam,message,unit=77,form='formatted') /=0 ) then
2469        MSG_ERROR(message)
2470      end if
2471      nstep1=0;test0=.false.
2472      do while (.not.test0)
2473        read(77,'(A)') strg
2474        test0=(strg(1:1)/="#")
2475        if (test0) read(unit=strg,fmt=*) nstep1
2476      end do
2477      nstep1_abs=abs(nstep1)
2478      if (nstep1_abs==0.or.istep>nstep1_abs.or.(nstep1>0.and.initialized/=0)) then
2479        close(77)
2480 !      Reinitalize mixing when rhoij is allowed to change; for experimental purpose...
2481        if (dtset%userib==1234.and.istep==1+nstep1_abs.and.(nstep1<0.or.initialized==0)) istep_mix=1
2482        iread=1
2483      end if
2484    end if
2485    call xmpi_sum(iread,mpi_comm_read,ierr)
2486    if (iread/=0) return
2487 
2488 !  Set up parallelism over atoms
2489    nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
2490    my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
2491    call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,my_natom_ref=my_natom)
2492 
2493 !  Store number of component for rhoij
2494    ABI_ALLOCATE(nspden_tmp,(natom))
2495    nspden_tmp(:)=zero
2496    do iatom=1,my_natom
2497      iatom_tot=my_atmtab(iatom)
2498      nspden_tmp(iatom_tot)=pawrhoij(iatom)%nspden
2499    end do
2500    call xmpi_sum(nspden_tmp,mpi_comm_read,ierr)
2501 
2502 !  To be improve if too much memory
2503    ABI_ALLOCATE(rhoijtmp2,(natom,4,rhoijsz,rhoijsz))
2504    rhoijtmp2=zero
2505 
2506 !  Read rhoij from file
2507    if (my_rank==0) then
2508      do iatom=1,natom
2509        itypat=typat(iatom)
2510        if (pawtab(itypat)%useexexch>0) then
2511          rhoijsz=4*ll+2
2512          do ispden=1,nspden_tmp(iatom)
2513            do ii=1,rhoijsz
2514              test0=.false.
2515              do while (.not.test0)
2516                read(77,'(A)') strg
2517                test0=(strg(1:1)/="#")
2518                if (test0)  read(unit=strg,fmt=*) (rhoijtmp2(iatom,ispden,ii,jj),jj=1,rhoijsz)
2519              end do
2520            end do
2521          end do
2522        end if
2523      end do
2524    end if
2525    call xmpi_sum(rhoijtmp2,mpi_comm_read,ierr)
2526 
2527 !  Now, distribute rhoij
2528    do iatom=1,my_natom
2529      iatom_tot=my_atmtab(iatom)
2530      itypat=pawrhoij(iatom)%itypat
2531 
2532      if (pawtab(itypat)%useexexch>0) then
2533 
2534 !      Set sizes depending on ll
2535        rhoijsz=4*ll+2
2536        rhoijshft=2*ll*ll
2537 
2538 !      Uncompress rhoij
2539        ABI_ALLOCATE(rhoijtmp,(pawrhoij(iatom)%lmn2_size,pawrhoij(iatom)%nspden))
2540        do ispden=1,pawrhoij(iatom)%nspden
2541          rhoijtmp=zero
2542          do irhoij=1,pawrhoij(iatom)%nrhoijsel
2543            klmn=pawrhoij(iatom)%rhoijselect(irhoij)
2544            rhoijtmp(klmn,ispden)=pawrhoij(iatom)%rhoijp(irhoij,ispden)
2545          end do
2546 
2547 !        Impose rhoij
2548          do jj=1,rhoijsz
2549            do ii=1,jj
2550              rhoijtmp((jj+rhoijshft)*((jj+rhoijshft)-1)/2+ii+rhoijshft,ispden)=rhoijtmp2(iatom_tot,ispden,ii,jj)
2551            end do
2552          end do
2553 
2554        end do
2555 
2556 !      Compress rhoij
2557        nselect=0
2558        do klmn=1,pawrhoij(iatom)%lmn2_size
2559          if (any(abs(rhoijtmp(klmn,:))>tol10)) then
2560            nselect=nselect+1
2561            do ispden=1,pawrhoij(iatom)%nspden
2562              pawrhoij(iatom)%rhoijp(nselect,ispden)=rhoijtmp(klmn,ispden)
2563            end do
2564            pawrhoij(iatom)%rhoijselect(nselect)=klmn
2565          end if
2566        end do
2567        pawrhoij(iatom)%nrhoijsel=nselect
2568        ABI_DEALLOCATE(rhoijtmp)
2569 
2570      end if ! useexexch>0
2571 
2572 !    Print new rhoij
2573      do ispden=1,pawrhoij(iatom)%nspden
2574        write(message,'(2a,i3,a)') ch10,'== Atom ',iatom,' == Imposed occupation matrix'
2575        if (pawrhoij(iatom)%nspden==1) write(message,fmt='(2a)')     trim(message)," for spin up =="
2576        if (pawrhoij(iatom)%nspden==2) write(message,fmt='(2a,i3,a)')trim(message)," for spin ",ispden," =="
2577        if (pawrhoij(iatom)%nspden==4) write(message,fmt='(4a)')     trim(message)," for component ", &
2578 &       trim(dspin(ispden+2*(pawrhoij(iatom)%nspden/4)))," =="
2579        call wrtout(std_out,message,'PERS')
2580        call pawio_print_ij(std_out,pawrhoij(iatom)%rhoijp(:,ispden),pawrhoij(iatom)%nrhoijsel,&
2581 &       pawrhoij(iatom)%cplex,pawrhoij(iatom)%lmn_size,ll,&
2582 &       pawtab(itypat)%indlmn(1,1:pawtab(itypat)%lmn_size),&
2583 &       1,-1,pawrhoij(iatom)%rhoijselect(:),-1.d0,1,mode_paral='PERS')
2584      end do
2585 
2586 !    end loop on atoms
2587    end do
2588 
2589    ABI_DEALLOCATE(nspden_tmp)
2590    ABI_DEALLOCATE(rhoijtmp2)
2591 
2592 !  Destroy atom table used for parallelism
2593    call free_my_atmtab(my_atmtab,my_atmtab_allocated)
2594 
2595 !  ============================================================
2596  end if ! paral_atom
2597 
2598  DBG_EXIT("COLL")
2599 
2600 end subroutine setrhoijpbe0