Line data Source code
1 : !!****m* m_paw_correlations/m_paw_correlations
2 : !! NAME
3 : !! m_paw_correlations
4 : !!
5 : !! FUNCTION
6 : !! This module contains several routines related to the treatment of electronic
7 : !! correlations in the PAW approach (DFT+U, exact-exchange, ...).
8 : !!
9 : !! COPYRIGHT
10 : !! Copyright (C) 2018-2026 ABINIT group (BA,FJ,MT)
11 : !! This file is distributed under the terms of the
12 : !! GNU General Public License, see ~abinit/COPYING
13 : !! or http://www.gnu.org/copyleft/gpl.txt .
14 : !!
15 : !! SOURCE
16 :
17 : #if defined HAVE_CONFIG_H
18 : #include "config.h"
19 : #endif
20 :
21 : #include "abi_common.h"
22 :
23 : MODULE m_paw_correlations
24 :
25 : use defs_basis
26 : use m_errors
27 : use m_abicore
28 : use m_xmpi
29 : use m_dtset
30 : use m_linalg_interfaces
31 : use m_special_funcs
32 : use m_fstrings, only : int2char4
33 : use m_io_tools, only : get_unit,open_file
34 : use m_paw_dmft, only : paw_dmft_type
35 : use m_pawang, only : pawang_type,pawang_init,pawang_free
36 : use m_pawrad, only : pawrad_free,pawrad_init,pawrad_type,simp_gen,nderiv_gen,pawrad_ifromr,poisson
37 : use m_pawtab, only : pawtab_type,pawtab_nullify,pawtab_free,pawtab_set_flags
38 : use m_pawrhoij, only : pawrhoij_type,pawrhoij_gather, pawrhoij_nullify, pawrhoij_free
39 : use m_paw_ij, only : paw_ij_type,paw_ij_gather, paw_ij_free, paw_ij_nullify
40 : use m_paw_sphharm, only : mat_mlms2jmj,mat_slm2ylm,slxyzs
41 : use m_paw_io, only : pawio_print_ij
42 : use m_paw_yukawa, only : compute_slater,get_lambda
43 : use m_paral_atom, only : get_my_atmtab,free_my_atmtab
44 : use m_copy, only : alloc_copy
45 : use m_geometry, only : vcart2ylm
46 : use m_atomdata
47 :
48 : implicit none
49 :
50 : private
51 :
52 : !public procedures.
53 : public :: pawpuxinit ! Initialize some data for PAW+U/PAW+LocalExactExchange/PAW+DMFT
54 : public :: calc_vee ! Compute vee for DFT+U
55 : public :: pawuenergy ! Compute contributions to energy for PAW+U
56 : public :: pawxenergy ! Compute contributions to energy for PAW+[local exact exchange]
57 : public :: setnoccmmp ! Compute DFT+U density matrix nocc_{m,m_prime} or impose it
58 : public :: setrhoijpbe0 ! Impose value of rhoij for using an auxiliairy file (PBE0 only)
59 : public :: calc_ubare ! Calculate the bare interaction on atomic orbitals
60 : public :: loc_orbmom_cal ! calculate local orbital magnetic moments
61 : CONTAINS !========================================================================================
62 : !!***
63 :
64 : !----------------------------------------------------------------------
65 :
66 : !!****f* m_paw_correlations/pawpuxinit
67 : !! NAME
68 : !! pawpuxinit
69 : !!
70 : !! FUNCTION
71 : !! Initialize some starting values of several arrays used in
72 : !! PAW+U/+DMFT or local exact-exchange calculations
73 : !!
74 : !! A-define useful indices for DFT+U/local exact-exchange
75 : !! B-Compute overlap between atomic wavefunction
76 : !! C-Compute matrix elements of coulomb interaction (see PRB vol.52 5467) [[cite:Liechenstein1995]]
77 : !! (angular part computed from Gaunt coefficients)
78 : !!
79 : !! INPUTS
80 : !! dmatpuopt= select expression for the density matrix
81 : !! dmft_dc= option for the double-counting scheme in DMFT
82 : !! dmft_orbital(ntypat)= option for the choice of the DMFT radial orbital
83 : !! dmft_orbital_filepath= name of the DMFT orbital file
84 : !! exchmix= mixing factor for local exact-exchange
85 : !! is_dfpt=true if we are running a DFPT calculation
86 : !! jpawu(ntypat)= value of J
87 : !! llexexch(ntypat)= value of l on which local exact-exchange applies
88 : !! llpawu(ntypat)= value of l on which DFT+U applies
89 : !! ntypat=number of types of atoms in unit cell.
90 : !! pawang <type(pawang_type)>=paw angular mesh and related data
91 : !! %lmax=Maximum value of angular momentum l+1
92 : !! %gntselect((2*l_max-1)**2,l_max**2,l_max**2)=
93 : !! selection rules for Gaunt coefficients
94 : !! pawprtvol=output printing level for PAW
95 : !! pawrad(ntypat) <type(pawrad_type)>=paw radial mesh and related data
96 : !! pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data:
97 : !! upawu(ntypat)= value of U
98 : !! use_dmft = 0 no PAW+DMFT, =1 PAW+DMFT
99 : !! useexexch= 0 if no local exact-exchange; 1 if local exact-exchange
100 : !! usepawu= 0 if no DFT+U; /=0 if DFT+U
101 : !!
102 : !! OUTPUT
103 : !! pawtab <type(pawtab_type)>=paw tabulated data read at start:
104 : !! %euijkl=(3,lmn2_size,lmn2_size)= array for computing DFT+U terms without occupancies
105 : !! %ij_proj= nproj*(nproju+1)/2
106 : !! %klmntomn(4,lmn2_size)= Array giving im, jm ,in, and jn for each klmn=(ilmn,jlmn)
107 : !! %lnproju(nproj)= value of ln for projectors on which paw+u/local exact-exchange acts.
108 : !! %nproju=number of projectors for orbitals on which paw+u/local exact-exchange acts.
109 : !! %phiphjint(pawtabitypat%ij_proj)=Integral of Phi(:,i)*Phi(:,j) for correlated orbitals.
110 : !! %usepawu=0 if no DFT+U; /=0 if DFT+U
111 : !! %useexexch=0 if no local exact-exchange; 1 if local exact-exchange
112 : !! === if usepawu/=0
113 : !! %jpawu= value of J
114 : !! %upawu= value of U
115 : !! %vee(2*lpawu+1,:,:,:)=matrix of the screened interaction for correlated orbitals
116 : !! === if useexexch/=0
117 : !! %fk
118 : !! %vex(2*lpawu+1,:,:,:)=matrix of the screened interaction for correlated orbitals
119 : !!
120 : !! SOURCE
121 :
122 1494 : subroutine pawpuxinit(dmatpuopt,exchmix,f4of2_sla,f6of2_sla,is_dfpt,jpawu,llexexch,llpawu,&
123 1494 : & nspinor,ntypat,option_interaction,pawang,pawprtvol,pawrad,pawtab,upawu,use_dmft,&
124 : & useexexch,usepawu,&
125 : & ucrpa,lmagCalc,dmft_orbital,dmft_dc,dmft_orbital_filepath,& ! optional argument
126 : & dmft_yukawa_param,dmft_yukawa_lambda,dmft_yukawa_epsilon) ! optional argument
127 :
128 : !Arguments ---------------------------------------------
129 : !scalars
130 : integer,intent(in) :: dmatpuopt,nspinor,ntypat,pawprtvol,use_dmft,useexexch,usepawu
131 : !Option for interaction energy in case of non-collinear magnetism:
132 : ! 1: E_int=-J/4.N.(N-2)
133 : ! 2: E_int=-J/2.(Nup.(Nup-1)+Ndn.(Ndn-1)) (Nup and Ndn are ill-defined)
134 : ! 3: E_int=-J/4.( N.(N-2) + mx^2 + my^2 + mz^2 )
135 : ! Default is 3
136 : integer,intent(in) :: option_interaction
137 : logical :: is_dfpt
138 : real(dp),intent(in) :: exchmix
139 : type(pawang_type), intent(in) :: pawang
140 : integer,optional, intent(in) :: dmft_dc,ucrpa
141 : !arrays
142 : integer,intent(in) :: llexexch(ntypat),llpawu(ntypat)
143 : real(dp),intent(in) :: jpawu(ntypat),upawu(ntypat)
144 : real(dp),intent(in) :: f4of2_sla(ntypat),f6of2_sla(ntypat)
145 : type(pawrad_type),intent(inout) :: pawrad(ntypat)
146 : type(pawtab_type),target,intent(inout) :: pawtab(ntypat)
147 : logical,optional,intent(in) :: lmagCalc
148 : integer,optional,intent(in) :: dmft_yukawa_param
149 : integer,optional,intent(in) :: dmft_orbital(ntypat)
150 : real(dp),optional,intent(in) :: dmft_yukawa_epsilon,dmft_yukawa_lambda
151 : character(len=fnlen),optional,intent(in) :: dmft_orbital_filepath
152 : !Local variables ---------------------------------------
153 : !scalars
154 : integer :: icount,ierr,il,ilmn,ilmnp,ir,isela,iselb,itemp,itypat,iu,iup,j0lmn,jl,jlmn,jlmnp,ju,jup
155 : integer :: klm0x,klma,klmb,klmn,klmna,klmnb,kln,kln1,kln2,kyc,lcur,lexexch,lkyc,ll,ll1
156 : integer :: lmexexch,lmkyc,lmn_size,lmn2_size,lpawu
157 : integer :: m1,m11,m2,m21,m3,m31,m4,m41
158 : integer :: me,mesh_size,mesh_type,meshsz,int_meshsz,mkyc,unt,sz1
159 : integer :: option_interaction_, Loc_prtvol
160 : logical :: compute_euijkl,compute_euij_fll,lexist
161 : real(dp) :: ak,eps,f4of2,f6of2,int1,intg,jh,lambda,lstep,phiint_ij,phiint_ipjp,rstep,uh,vee1,vee2
162 : character(len=4) :: tag,tag2
163 : character(len=500) :: message,tmpfil
164 : logical :: lmagCalc_
165 : !arrays
166 1494 : integer,ABI_CONTIGUOUS pointer :: indlmn(:,:)
167 : real(dp) :: euijkl_temp(3),euijkl_temp2(3),euijkl_dc(3)
168 1494 : real(dp),allocatable :: ff(:),fk(:),gg(:)
169 1494 : type(pawrad_type) :: pawrad_tmp
170 : ! *************************************************************************
171 :
172 : DBG_ENTER("COLL")
173 1494 : Loc_prtvol = 3
174 1494 : lmagCalc_ = .False.
175 1494 : if (present(lmagCalc)) then
176 72 : if (lmagCalc .eqv. .True.) lmagCalc_ = .True.
177 72 : Loc_prtvol = 0
178 : end if
179 :
180 : !No correlations= nothing to do
181 1494 : if(useexexch==0.and.usepawu==0.and.use_dmft==0) then
182 3059 : do itypat=1,ntypat
183 3059 : pawtab(itypat)%usepawu=0;pawtab(itypat)%useexexch=0;pawtab(itypat)%exchmix=zero
184 : end do
185 : return
186 : end if
187 :
188 : !PAW+U and local exact-exchange restriction
189 235 : if(useexexch/=0.and.usepawu/=0)then
190 0 : do itypat=1,ntypat
191 0 : if (llpawu(itypat)/=llexexch(itypat).and.llpawu(itypat)/=-1.and.llexexch(itypat)/=-1) then
192 : write(message, '(5a,i2,3a)' )&
193 0 : & ' When PAW+U (usepawu/=0) and local exact-exchange (exexch/=0)',ch10,&
194 0 : & ' are selected together, they must apply on the same',ch10,&
195 0 : & ' angular momentum (lpawu/=lexexch forbidden, here for typat=',itypat,') !',ch10,&
196 0 : & ' Action: correct your input file.'
197 0 : ABI_ERROR(message)
198 : end if
199 : end do
200 : end if
201 :
202 : !Print title
203 53 : if((abs(usepawu)>=1.and.abs(usepawu)<=4).or.useexexch/=0.and.(.not.lmagCalc_)) &
204 236 : & write(message, '(3a)' ) ch10,ch10," ******************************************"
205 : if(usepawu==1) then
206 143 : write(message, '(3a)' ) trim(message),ch10," DFT+U Method used: FLL"
207 : else if(usepawu==2) then
208 4 : write(message, '(3a)' ) trim(message),ch10," DFT+U Method used: AMF"
209 : else if(usepawu==3) then
210 0 : write(message, '(3a)' ) trim(message),ch10," DFT+U Method used: AMF (alternative)"
211 : else if(usepawu==4) then
212 24 : write(message, '(3a)' ) trim(message),ch10," DFT+U Method used: FLL with no spin polarization in the xc functional"
213 : else if(usepawu==-1) then
214 5 : write(message, '(3a)' ) trim(message),ch10," DFT+U Method used: FLL (no use of occupation matrix) - experimental"
215 : else if(usepawu==-2) then
216 2 : write(message, '(3a)' ) trim(message),ch10," DFT+U Method used: AMF (no use of occupation matrix) - experimental"
217 : else if(usepawu==-4) then
218 4 : write(message, '(3a)' ) trim(message),ch10," DFT+U Method used: FLL with no spin polarization in the xc functional &
219 8 : & (no use of occupation matrix) - experimental"
220 : end if
221 235 : if(useexexch/=0) write(message, '(3a)' ) trim(message),ch10," PAW Local Exact exchange: PBE0"
222 235 : if((abs(usepawu)>=1.and.abs(usepawu)<=4).or.useexexch/=0 .and.(.not.lmagCalc_)) then
223 183 : if (nspinor==2) then
224 99 : write(message, '(3a,i1)' ) trim(message),ch10," Magnetic DC : option_interaction = ",option_interaction
225 : end if
226 183 : write(message, '(3a)' ) trim(message),ch10," ******************************************"
227 : end if
228 235 : if(use_dmft==0 .and. abs(usepawu)<=4 .and.(.not.lmagCalc_)) then
229 111 : call wrtout(ab_out,message,'COLL')
230 111 : call wrtout(std_out, message,'COLL')
231 : end if
232 : !if(use_dmft>0) then
233 : !write(message, '(3a)' ) ch10, " (see DMFT data in log file) "
234 : !call wrtout(ab_out,message,'COLL')
235 : !endif
236 235 : option_interaction_ = option_interaction
237 235 : if(abs(usepawu)>=10.and.nspinor==2.and.option_interaction/=1 .and.(.not.lmagCalc_)) then
238 1 : option_interaction_ = 1
239 1 : write(message, '(a)' ) "When usepawu>=10, option_interaction for DC is set to 1"
240 1 : call wrtout(std_out,message,'COLL')
241 : end if
242 235 : if(usepawu<0.and.nspinor==2.and.option_interaction_==2 .and.(.not.lmagCalc_)) then
243 0 : write(message, '(a)' ) "option_interaction=2 is not implemented for usepawu<0. Change 'usepawu' or 'optdcmagpawu' in the input."
244 0 : ABI_ERROR(message)
245 : end if
246 :
247 : !Loop on atom types
248 635 : do itypat=1,ntypat
249 400 : indlmn => pawtab(itypat)%indlmn
250 400 : lmn_size=pawtab(itypat)%lmn_size
251 400 : lmn2_size=pawtab(itypat)%lmn2_size
252 400 : mesh_size=pawtab(itypat)%mesh_size
253 400 : int_meshsz=pawrad(itypat)%int_meshsz
254 400 : lcur=-1
255 :
256 400 : if (use_dmft > 0) then
257 : if (dmft_dc == 8 .and. (f4of2_sla(itypat) >= -0.1_dp .or. &
258 89 : & f6of2_sla(itypat) >= -0.1_dp) .and. dmft_yukawa_param == 1) then
259 0 : message = "dmft_dc=8 and dmft_yukawa_param=1 not compatible with custom f4of2 and f6of2"
260 0 : ABI_ERROR(message)
261 : end if
262 : end if
263 :
264 : ! PAW+U data
265 400 : if (usepawu/=0.or.use_dmft>0) then
266 398 : lcur=llpawu(itypat)
267 398 : pawtab(itypat)%lpawu=lcur
268 398 : if(lcur/=-1) then
269 300 : pawtab(itypat)%usepawu=usepawu
270 300 : pawtab(itypat)%upawu=upawu(itypat)
271 300 : pawtab(itypat)%jpawu=jpawu(itypat)
272 300 : pawtab(itypat)%f4of2_sla=f4of2_sla(itypat)
273 300 : pawtab(itypat)%f6of2_sla=f6of2_sla(itypat)
274 300 : pawtab(itypat)%option_interaction_pawu=option_interaction_
275 : else
276 98 : pawtab(itypat)%usepawu=0
277 98 : pawtab(itypat)%upawu=zero
278 98 : pawtab(itypat)%jpawu=zero
279 98 : pawtab(itypat)%f4of2_sla=zero
280 98 : pawtab(itypat)%f6of2_sla=zero
281 98 : pawtab(itypat)%option_interaction_pawu=option_interaction_
282 : end if
283 : end if
284 :
285 : ! Local exact-echange data
286 400 : if (useexexch/=0) then
287 2 : lcur=llexexch(itypat)
288 2 : pawtab(itypat)%lexexch=lcur
289 2 : pawtab(itypat)%exchmix=exchmix
290 2 : if(pawtab(itypat)%lexexch==-1) pawtab(itypat)%useexexch=0
291 2 : if(pawtab(itypat)%lexexch/=-1) pawtab(itypat)%useexexch=useexexch
292 : end if
293 :
294 : ! Select only atoms with +U
295 635 : if(lcur/=-1) then
296 :
297 : ! Compute number of projectors for DFT+U/local exact-exchange/DFT+DMFT
298 5516 : icount=count(indlmn(1,1:lmn_size)==lcur)
299 301 : pawtab(itypat)%nproju=icount/(2*lcur+1)
300 301 : if(useexexch/=0.and.pawtab(itypat)%nproju>2) then
301 : write(message, '(a,a,a)' )&
302 0 : & ' Error on the number of projectors ',ch10,&
303 0 : & ' more than 2 projectors is not allowed for local exact-exchange'
304 0 : ABI_ERROR(message)
305 : end if
306 301 : if(pawtab(itypat)%nproju*(2*lcur+1)/=icount) then
307 0 : message = 'pawpuxinit: Error on the number of projectors '
308 0 : ABI_BUG(message)
309 : end if
310 301 : if ((.not.lmagCalc_)) then
311 175 : write(message, '(a,a,i4,a,a,i4)' ) ch10,&
312 175 : & ' pawpuxinit : for species ',itypat,ch10,&
313 350 : & ' number of projectors is',pawtab(itypat)%nproju
314 175 : call wrtout(std_out,message,'COLL')
315 :
316 : end if
317 301 : pawtab(itypat)%ij_proj=pawtab(itypat)%nproju*(pawtab(itypat)%nproju+1)/2
318 :
319 : ! ==================================================
320 : ! A-define useful indexes
321 : ! --------------------------------------------------
322 301 : if (allocated(pawtab(itypat)%lnproju)) then
323 203 : ABI_FREE(pawtab(itypat)%lnproju)
324 : end if
325 903 : ABI_MALLOC(pawtab(itypat)%lnproju,(pawtab(itypat)%nproju))
326 301 : icount=0
327 5516 : do ilmn=1,lmn_size
328 5516 : if(indlmn(1,ilmn)==lcur) then
329 2351 : itemp=icount/(2*lcur+1)
330 2351 : if (itemp*(2*lcur+1)==icount) then
331 543 : pawtab(itypat)%lnproju(itemp+1)=indlmn(5,ilmn)
332 : end if
333 2351 : icount=icount+1
334 : end if
335 : end do
336 :
337 301 : if (allocated(pawtab(itypat)%klmntomn)) then
338 203 : ABI_FREE(pawtab(itypat)%klmntomn)
339 : end if
340 903 : ABI_MALLOC(pawtab(itypat)%klmntomn,(4,lmn2_size))
341 5516 : do jlmn=1,lmn_size
342 5215 : jl= indlmn(1,jlmn)
343 5215 : j0lmn=jlmn*(jlmn-1)/2
344 57758 : do ilmn=1,jlmn
345 52242 : il= indlmn(1,ilmn)
346 52242 : klmn=j0lmn+ilmn
347 52242 : pawtab(itypat)%klmntomn(1,klmn)=indlmn(2,ilmn)+il+1
348 52242 : pawtab(itypat)%klmntomn(2,klmn)=indlmn(2,jlmn)+jl+1
349 52242 : pawtab(itypat)%klmntomn(3,klmn)=indlmn(3,ilmn)
350 57457 : pawtab(itypat)%klmntomn(4,klmn)=indlmn(3,jlmn)
351 : end do
352 : end do
353 :
354 : ! ==================================================
355 : ! B-PAW+U: overlap between atomic wavefunctions
356 : ! --------------------------------------------------
357 301 : if(dmatpuopt==1 .and.(.not.lmagCalc_)) then
358 68 : write(message, '(4a)' ) ch10,&
359 68 : & ' pawpuxinit : dmatpuopt=1 ',ch10,&
360 136 : & ' PAW+U: dens. mat. constructed by projection on atomic wfn inside PAW augm. region(s)'
361 68 : call wrtout(std_out,message,'COLL')
362 68 : write(message, '(8a)' ) ch10,&
363 68 : & ' pawpuxinit: WARNING: Check that the first partial wave for lpawu:', ch10, &
364 68 : & ' - Is an atomic eigenfunction ',ch10, &
365 68 : & ' - Is normalized ',ch10, &
366 136 : & ' In other cases, choose dmatpuopt=2'
367 68 : call wrtout(std_out,message,'COLL')
368 233 : else if(dmatpuopt==2 .and.(.not.lmagCalc_)) then
369 83 : write(message, '(6a)' ) ch10,&
370 83 : & ' pawpuxinit : dmatpuopt=2 ',ch10,&
371 83 : & ' PAW+U: dens. mat. constructed by selecting contribution',ch10,&
372 166 : & ' for each angular momentum to the density (inside PAW augm. region(s))'
373 83 : call wrtout(std_out,message,'COLL')
374 150 : else if(dmatpuopt==3 .and.(.not.lmagCalc_)) then
375 23 : write(message, '(a,a,a,a,a,a)' ) ch10,&
376 23 : & ' pawpuxinit : dmatpuopt=3 ',ch10,&
377 23 : & ' PAW+U: dens. mat. constructed by projection on atomic wfn inside PAW augm. region(s)',ch10,&
378 46 : & ' and normalized inside PAW augm. region(s)'
379 23 : call wrtout(std_out,message,'COLL')
380 23 : write(message, '(6a)' ) ch10,&
381 23 : & ' pawpuxinit: WARNING: Check that the first partial wave for lpawu:', ch10, &
382 23 : & ' is an atomic eigenfunction',ch10, &
383 46 : & ' In the other case, choose dmatpuopt=2'
384 23 : call wrtout(std_out,message,'COLL')
385 : end if
386 :
387 903 : ABI_MALLOC(ff,(mesh_size))
388 191194 : ff(:)=zero
389 :
390 301 : if (allocated(pawtab(itypat)%ph0phiint)) then
391 203 : ABI_FREE(pawtab(itypat)%ph0phiint)
392 : end if
393 301 : if (allocated(pawtab(itypat)%zioneff)) then
394 203 : ABI_FREE(pawtab(itypat)%zioneff)
395 : end if
396 903 : ABI_MALLOC(pawtab(itypat)%ph0phiint,(pawtab(itypat)%nproju))
397 602 : ABI_MALLOC(pawtab(itypat)%zioneff,(pawtab(itypat)%nproju))
398 :
399 301 : icount=0
400 844 : do iu=1,pawtab(itypat)%nproju
401 : ! write(std_out,*)'DJA iu',iu,' mesh_size',pawtab(itypat)%mesh_size
402 : ! do ju=2,pawtab(itypat)%mesh_size
403 : ! ff(ju)=pawtab(itypat)%phi(ju,pawtab(itypat)%lnproju(iu))/pawrad(itypat)%rad(ju)
404 : ! write(std_out,fmt='(i5,3e15.5)')ju,pawrad(itypat)%rad(ju),ff(ju),&
405 : ! & RadFnH(pawrad(itypat)%rad(ju),4,3,15.0_dp)
406 : ! end do
407 : ! ff(1:mesh_size)=pawtab(itypat)%phi(1:mesh_size,pawtab(itypat)%lnproju(iu))**2
408 : ! call simp_gen(int1,ff,pawrad(itypat))
409 : ! write(std_out,*)'DJA iu',iu,'int1 ',int1
410 : ! write(std_out,*)'DJA int1,IRadFnH',int1,IRadFnH(0.0_dp,pawrad(itypat)%rmax,4,3,12)
411 : ! Calculation of zioneff
412 543 : ju=pawtab(itypat)%mesh_size-1
413 543 : ak=pawtab(itypat)%phi(ju,pawtab(itypat)%lnproju(iu))/pawtab(itypat)%phi(ju+1,pawtab(itypat)%lnproju(iu))
414 543 : ak=ak*(pawrad(itypat)%rad(ju+1)/pawrad(itypat)%rad(ju))**(pawtab(itypat)%lpawu-1)
415 543 : pawtab(itypat)%zioneff(iu)=log(ak)/(pawrad(itypat)%rad(ju+1)-pawrad(itypat)%rad(ju))
416 : ! Calculation of ph0phiint
417 : ff(1:mesh_size)=pawtab(itypat)%phi(1:mesh_size,pawtab(itypat)%lnproju(1))&
418 361679 : & *pawtab(itypat)%phi(1:mesh_size,pawtab(itypat)%lnproju(iu))
419 543 : call simp_gen(int1,ff,pawrad(itypat))
420 844 : pawtab(itypat)%ph0phiint(iu)=int1
421 : end do
422 301 : if(abs(pawprtvol)>=2) then
423 287 : do icount=1,pawtab(itypat)%nproju
424 188 : write(message, '(a,a,i2,f9.5,a)' ) ch10,&
425 376 : & ' pawpuxinit: icount, ph0phiint(icount)=',icount,pawtab(itypat)%ph0phiint(icount)
426 188 : call wrtout(std_out,message,'COLL')
427 : write(message, '(a,f15.5)' ) &
428 188 : & ' pawpuxinit: zioneff=',pawtab(itypat)%zioneff(icount)
429 287 : call wrtout(std_out,message,'COLL')
430 : end do
431 99 : write(message, '(a)' ) ch10
432 99 : call wrtout(std_out,message,'COLL')
433 : end if
434 :
435 301 : if (allocated(pawtab(itypat)%phiphjint)) then
436 203 : ABI_FREE(pawtab(itypat)%phiphjint)
437 : end if
438 903 : ABI_MALLOC(pawtab(itypat)%phiphjint,(pawtab(itypat)%ij_proj))
439 :
440 301 : icount=0
441 844 : do ju=1,pawtab(itypat)%nproju
442 1629 : do iu=1,ju
443 785 : icount=icount+1
444 1328 : if ((dmatpuopt==1).and.(useexexch==0)) then
445 : pawtab(itypat)%phiphjint(icount)=pawtab(itypat)%ph0phiint(iu)*&
446 144 : & pawtab(itypat)%ph0phiint(ju)
447 641 : else if((dmatpuopt==2).or.(useexexch/=0)) then
448 : ff(1:mesh_size)=pawtab(itypat)%phi(1:mesh_size,pawtab(itypat)%lnproju(iu))&
449 334164 : & *pawtab(itypat)%phi(1:mesh_size,pawtab(itypat)%lnproju(ju))
450 569 : call simp_gen(int1,ff,pawrad(itypat))
451 569 : pawtab(itypat)%phiphjint(icount)=int1
452 72 : else if((dmatpuopt>=3).and.(useexexch==0)) then
453 : pawtab(itypat)%phiphjint(icount)=pawtab(itypat)%ph0phiint(iu)* &
454 72 : & pawtab(itypat)%ph0phiint(ju)/pawtab(itypat)%ph0phiint(1)**(dmatpuopt-2)
455 : else
456 : write(message, '(3a)' )&
457 0 : & ' PAW+U: dmatpuopt has a wrong value !',ch10,&
458 0 : & ' Action : change value in input file'
459 0 : ABI_ERROR(message)
460 : end if
461 : end do
462 : end do
463 301 : if(pawtab(itypat)%ij_proj/=icount) then
464 0 : message = ' Error in the loop for calculating phiphjint '
465 0 : ABI_ERROR(message)
466 : end if
467 301 : ABI_FREE(ff)
468 301 : if(abs(pawprtvol)>=2) then
469 376 : do icount=1,pawtab(itypat)%ij_proj
470 277 : write(message, '(a,a,i2,f9.5,a)' ) ch10,&
471 554 : & ' PAW+U: icount, phiphjint(icount)=',icount,pawtab(itypat)%phiphjint(icount)
472 376 : call wrtout(std_out,message,'COLL')
473 : end do
474 : end if
475 : ! end if
476 :
477 : ! ======================================================================
478 : ! C-PAW+U: Matrix elements of coulomb interaction (see PRB vol.52 5467) [[cite:Liechenstein1995]]
479 : ! 1. angular part computed from Gaunt coefficients
480 : ! --------------------------------------------------------------------
481 301 : if (usepawu/=0) then
482 300 : lpawu=lcur
483 :
484 300 : if (allocated(pawtab(itypat)%vee)) then
485 203 : ABI_FREE(pawtab(itypat)%vee)
486 : end if
487 300 : sz1=2*lpawu+1
488 1800 : ABI_MALLOC(pawtab(itypat)%vee,(sz1,sz1,sz1,sz1))
489 :
490 : call calc_vee(pawtab(itypat)%f4of2_sla,pawtab(itypat)%f6of2_sla,pawtab(itypat)%jpawu,&
491 300 : & pawtab(itypat)%lpawu,pawang,pawtab(itypat)%upawu,pawtab(itypat)%vee,Loc_prtvol)
492 :
493 : ! testu=0
494 : ! write(std_out,*) " Matrix of interaction vee(m1,m2,m1,m2)"
495 : ! do m1=1,2*lpawu+1
496 : ! write(std_out,'(2x,14(f12.6,2x))') (pawtab(itypat)%vee(m1,m2,m1,m2),m2=1,2*lpawu+1)
497 : ! do m2=1,2*lpawu+1
498 : ! testu=testu+ pawtab(itypat)%vee(m1,m2,m1,m2)
499 : ! enddo
500 : ! enddo
501 300 : if (.not.lmagCalc_) then
502 174 : write(message,'(a)') ch10
503 174 : call wrtout(std_out,message,'COLL')
504 174 : write(message,'(a)') " Matrix of interaction vee(m1,m2,m1,m2)"
505 174 : call wrtout(std_out,message,'COLL')
506 1078 : do m1=1,2*lpawu+1
507 5710 : write(message,'(2x,14(f20.14,2x))') (pawtab(itypat)%vee(m1,m2,m1,m2)*Ha_eV,m2=1,2*lpawu+1)
508 1078 : call wrtout(std_out,message,'COLL')
509 : ! do m2=1,2*lpawu+1
510 : ! testu=testu+ pawtab(itypat)%vee(m1,m2,m1,m2)
511 : ! enddo
512 : enddo
513 174 : write(message,'(a)') ch10
514 174 : call wrtout(std_out,message,'COLL')
515 : end if
516 :
517 : ! testu=testu/((two*lpawu+one)**2)
518 : ! write(std_out,*) "------------------------"
519 : ! write(std_out,'(a,f12.6)') " U=", testu
520 : ! write(std_out,*) "------------------------"
521 : ! write(std_out,*) " Matrix of interaction vee(m1,m2,m1,m2)-vee(m1,m2,m2,m1)"
522 : ! do m1=1,2*lpawu+1
523 : ! 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)
524 : ! do m2=1,2*lpawu+1
525 : ! if(m1/=m2) testumj=testumj+ pawtab(itypat)%vee(m1,m2,m1,m2)-pawtab(itypat)%vee(m1,m2,m2,m1)
526 : ! enddo
527 : ! enddo
528 : ! testumj=testumj/((two*lpawu)*(two*lpawu+one))
529 : ! write(std_out,*) "------------------------"
530 : ! write(std_out,'(a,f12.6)') " U-J=", testumj
531 : ! write(std_out,*) "------------------------"
532 : ! write(std_out,*) "------------------------"
533 : ! write(std_out,'(a,f12.6)') " J=", testu-testumj
534 : ! write(std_out,*) "------------------------"
535 :
536 : ! c. For DFPT (or with exp. values usepawu=-1,-2 or -4), compute euijkl
537 : ! ---------------------------------------------
538 300 : compute_euijkl=(is_dfpt.or.usepawu<0)
539 : if (compute_euijkl) then
540 11 : if (allocated(pawtab(itypat)%euijkl)) then
541 7 : ABI_FREE(pawtab(itypat)%euijkl)
542 : end if
543 66 : ABI_MALLOC(pawtab(itypat)%euijkl,(3,lmn_size,lmn_size,lmn_size,lmn_size))
544 4377414 : pawtab(itypat)%euijkl = zero
545 11 : compute_euij_fll = .false.
546 11 : euijkl_temp2=zero
547 11 : if (abs(usepawu)==1.or.abs(usepawu)==4) then ! Only for FLL
548 9 : if (allocated(pawtab(itypat)%euij_fll)) then ! allocate euij_fll for FLL
549 5 : ABI_FREE(pawtab(itypat)%euij_fll)
550 : end if
551 27 : ABI_MALLOC(pawtab(itypat)%euij_fll,(lmn2_size))
552 1468 : pawtab(itypat)%euij_fll = zero
553 : compute_euij_fll = .true.
554 : end if
555 :
556 : ! loop on i,j
557 1812 : do klmna=1,lmn2_size
558 1801 : ilmn=pawtab(itypat)%indklmn(7,klmna) ! i
559 1801 : jlmn=pawtab(itypat)%indklmn(8,klmna) ! j
560 1812 : if (pawtab(itypat)%indlmn(1,ilmn)==lpawu.and.pawtab(itypat)%indlmn(1,jlmn)==lpawu) then ! only correlated orbitals
561 565 : iu = pawtab(itypat)%indlmn(3,ilmn) ! ni
562 565 : ju = pawtab(itypat)%indlmn(3,jlmn) ! nj
563 565 : phiint_ij = pawtab(itypat)%phiphjint(iu+(ju*(ju-1))/2) ! iu <= ju by construction (ilmn<=jlmn)
564 565 : m2 = pawtab(itypat)%indlmn(2,ilmn) ! mi
565 565 : m21=m2+lpawu+1
566 565 : m1 = pawtab(itypat)%indlmn(2,jlmn) ! mj
567 565 : m11=m1+lpawu+1
568 :
569 565 : if (compute_euij_fll.and.m1==m2) then ! FLL
570 125 : pawtab(itypat)%euij_fll(klmna) = - half * phiint_ij * ( pawtab(itypat)%jpawu - pawtab(itypat)%upawu )
571 : end if
572 :
573 : ! loop on ip,jp (=k,l)
574 95980 : do klmnb=1,lmn2_size
575 95415 : ilmnp=pawtab(itypat)%indklmn(7,klmnb) ! ip (=k)
576 95415 : jlmnp=pawtab(itypat)%indklmn(8,klmnb) ! jp (=l)
577 95980 : if (pawtab(itypat)%indlmn(1,ilmnp)==lpawu.and.pawtab(itypat)%indlmn(1,jlmnp)==lpawu) then ! correlated orbitals
578 30475 : iup = pawtab(itypat)%indlmn(3,ilmnp) ! nip
579 30475 : jup = pawtab(itypat)%indlmn(3,jlmnp) ! njp
580 30475 : phiint_ipjp = pawtab(itypat)%phiphjint(iup+(jup*(jup-1))/2) ! iup <= jup by construction (ilmnp<=jlmnp)
581 30475 : m4 = pawtab(itypat)%indlmn(2,ilmnp) ! mip
582 30475 : m41=m4+lpawu+1
583 30475 : m3 = pawtab(itypat)%indlmn(2,jlmnp) ! mjp
584 30475 : m31=m3+lpawu+1
585 :
586 30475 : euijkl_dc(:) = zero
587 : ! Compute the double-counting part of euijkl (invariant when exchanging i<-->j or ip<-->jp)
588 : ! Must be consistent with pawuenergy and pawpupot
589 30475 : if (m1==m2.and.m3==m4) then ! In that case, we have to add the double-counting term
590 :
591 2275 : if (abs(usepawu)==1.and.nspinor==1) then ! FLL
592 :
593 : euijkl_dc(1) = &
594 450 : & phiint_ij * phiint_ipjp * ( pawtab(itypat)%upawu - pawtab(itypat)%jpawu )
595 : euijkl_dc(2) = &
596 450 : & phiint_ij * phiint_ipjp * pawtab(itypat)%upawu
597 :
598 1825 : else if (abs(usepawu)==2.and.nspinor==1) then ! AMF
599 :
600 : euijkl_dc(1) = &
601 450 : & two*lpawu/(two*lpawu+one) * phiint_ij * phiint_ipjp * ( pawtab(itypat)%upawu - pawtab(itypat)%jpawu )
602 : euijkl_dc(2) = &
603 450 : & phiint_ij * phiint_ipjp * pawtab(itypat)%upawu
604 :
605 1375 : else if (abs(usepawu)==4.or.nspinor>1) then ! FLL without polarization in XC or nspinor>1
606 :
607 : euijkl_dc(1:2) = &
608 4125 : & phiint_ij * phiint_ipjp * ( pawtab(itypat)%upawu - half*pawtab(itypat)%jpawu )
609 :
610 : ! Add term taking into account global magnetization
611 1375 : if (abs(usepawu)/=4.and.pawtab(itypat)%option_interaction_pawu==3) then
612 :
613 250 : euijkl_dc(1) = euijkl_dc(1) - half * phiint_ij * phiint_ipjp * pawtab(itypat)%jpawu
614 250 : euijkl_dc(2) = euijkl_dc(2) + half * phiint_ij * phiint_ipjp * pawtab(itypat)%jpawu
615 250 : euijkl_dc(3) = eUijkl_dc(3) - phiint_ij * phiint_ipjp * pawtab(itypat)%jpawu
616 :
617 : end if
618 :
619 : end if
620 :
621 : end if ! double-counting term
622 :
623 : ! Array of size 3:
624 : ! 1st element : coupled with up/up or down/down terms
625 : ! 2nd element : coupled with up/down or down/up terms (collinear part)
626 : ! 3rd element : coupled with up/down or down/up terms (non-collinear part)
627 30475 : euijkl_temp(:) = zero
628 30475 : euijkl_temp2(:) = zero
629 :
630 30475 : vee1 = pawtab(itypat)%vee(m11,m31,m21,m41)
631 : ! Note : vee(13|24) = vee(23|14) ( so : i <--> j )
632 : ! vee(13|24) = vee(14|23) ( so : ip <--> jp )
633 : ! vee(13|24) = vee(24|13) ( so : i,ip <--> j,jp )
634 : ! Also : vee(13|24) = vee(31|42) ( so : i,j <--> ip,jp )
635 : ! ==> vee1 is invariant with respect to the permutations i <--> j , ip <--> jp and i,ip <--> j,jp
636 : ! ( The term 'phiint_ij * phiint_ipjp' has the same properties)
637 91425 : euijkl_temp(1:2) = phiint_ij * phiint_ipjp * vee1
638 :
639 30475 : vee2 = pawtab(itypat)%vee(m11,m31,m41,m21)
640 : ! Note : vee(13|42) = vee(43|12) ( so : ip <--> j )
641 : ! vee(13|42) = vee(12|43) ( so : i <--> jp )
642 : ! vee(13|42) = vee(42|13) ( so : i,ip <--> jp,j )
643 : ! Also : vee(13|42) = vee(31|24) ( so : i,j <--> ip,jp )
644 : ! Combining the third and fourth rule we get:
645 : ! vee(13|42) = vee(42|13) = vee(24|31) ( so : i,ip <--> j,jp )
646 : ! ==> vee2 is invariant only with respect to the permutation i,ip <--> j,jp
647 :
648 : ! Terms i,j,ip,jp (m2,m1,m4,m3) and j,i,jp,ip (m1,m2,m3,m4)
649 30475 : euijkl_temp2(1) = phiint_ij * phiint_ipjp * vee2
650 30475 : euijkl_temp2(3) = phiint_ij * phiint_ipjp * vee2
651 121900 : pawtab(itypat)%euijkl(:,ilmn,jlmn,ilmnp,jlmnp) = euijkl_temp(:) - euijkl_temp2(:) - euijkl_dc(:)
652 121900 : pawtab(itypat)%euijkl(:,jlmn,ilmn,jlmnp,ilmnp) = pawtab(itypat)%euijkl(:,ilmn,jlmn,ilmnp,jlmnp)
653 :
654 : ! Term j,i,ip,jp (m1,m2,m4,m3)
655 30475 : vee2 = pawtab(itypat)%vee(m21,m31,m41,m11)
656 30475 : euijkl_temp2(1) = phiint_ij * phiint_ipjp * vee2
657 30475 : euijkl_temp2(3) = phiint_ij * phiint_ipjp * vee2
658 121900 : pawtab(itypat)%euijkl(:,jlmn,ilmn,ilmnp,jlmnp) = euijkl_temp(:) - euijkl_temp2(:) - euijkl_dc(:)
659 :
660 : ! Term i,j,jp,ip (m2,m1,m3,m4)
661 30475 : vee2 = pawtab(itypat)%vee(m11,m41,m31,m21)
662 30475 : euijkl_temp2(1) = phiint_ij * phiint_ipjp * vee2
663 30475 : euijkl_temp2(3) = phiint_ij * phiint_ipjp * vee2
664 121900 : pawtab(itypat)%euijkl(:,ilmn,jlmn,jlmnp,ilmnp) = euijkl_temp(:) - euijkl_temp2(:) - euijkl_dc(:)
665 :
666 : end if ! correlated orbitals
667 : end do ! klmnb
668 : end if ! correlated orbitals
669 : end do ! klmna
670 :
671 : end if ! compute_euijkl
672 : end if ! usepawu
673 :
674 : ! ======================================================================
675 : ! D-Local ex-exchange: Matrix elements of coulomb interaction and Fk
676 : ! ----------------------------------------------------------------------
677 301 : if (useexexch/=0) then
678 1 : lexexch=lcur
679 :
680 : ! a. compute F(k)
681 : ! ---------------------------------------------
682 1 : if (allocated(pawtab(itypat)%fk)) then
683 0 : ABI_FREE(pawtab(itypat)%fk)
684 : end if
685 1 : ABI_MALLOC(pawtab(itypat)%fk,(6,4))
686 29 : pawtab(itypat)%fk=zero
687 2 : ABI_MALLOC(ff,(mesh_size))
688 2 : ABI_MALLOC(gg,(mesh_size))
689 2161 : ff(:)=zero;gg(:)=zero
690 1 : kln=(pawtab(itypat)%lnproju(1)*( pawtab(itypat)%lnproju(1)+1)/2)
691 4 : do ll=1,lexexch+1
692 3 : ll1=2*ll-2
693 18 : if (int_meshsz<mesh_size) ff(int_meshsz+1:mesh_size)=zero
694 3228 : ff(1:int_meshsz)=pawtab(itypat)%phiphj(1:int_meshsz,kln)
695 3 : call poisson(ff,ll1,pawrad(itypat),gg)
696 3 : ff(1)=zero
697 : ff(2:mesh_size)=(pawtab(itypat)%phiphj(2:mesh_size,kln)*gg(2:mesh_size))&
698 3240 : & /pawrad(itypat)%rad(2:mesh_size)
699 3 : call simp_gen(intg,ff,pawrad(itypat))
700 4 : pawtab(itypat)%fk(1,ll)=intg*(two*ll1+one)
701 : end do
702 1 : if (pawtab(itypat)%nproju==2) then
703 1 : kln1=kln+pawtab(itypat)%lnproju(1)
704 1 : kln2=kln1+1
705 4 : do ll=1,lexexch+1
706 3 : ll1=2*ll-2
707 18 : if (int_meshsz<mesh_size) ff(int_meshsz+1:mesh_size)=zero
708 3228 : ff(1:int_meshsz)=pawtab(itypat)%phiphj(1:int_meshsz,kln1)
709 3 : call poisson(ff,ll1,pawrad(itypat),gg)
710 3 : ff(1)=zero
711 : ff(2:mesh_size)=(pawtab(itypat)%phiphj(2:mesh_size,kln1)*gg(2:mesh_size))&
712 3240 : & /pawrad(itypat)%rad(2:mesh_size)
713 3 : call simp_gen(intg,ff,pawrad(itypat))
714 4 : pawtab(itypat)%fk(2,ll)=intg*(two*ll1+one)
715 : end do
716 4 : do ll=1,lexexch+1
717 3 : ll1=2*ll-2
718 18 : if (int_meshsz<mesh_size) ff(int_meshsz+1:mesh_size)=zero
719 3228 : ff(1:int_meshsz)=pawtab(itypat)%phiphj(1:int_meshsz,kln2)
720 3 : call poisson(ff,ll1,pawrad(itypat),gg)
721 3 : ff(1)=zero
722 : ff(2:mesh_size)=(pawtab(itypat)%phiphj(2:mesh_size,kln2)*gg(2:mesh_size))&
723 3240 : & /pawrad(itypat)%rad(2:mesh_size)
724 3 : call simp_gen(intg,ff,pawrad(itypat))
725 4 : pawtab(itypat)%fk(3,ll)=intg*(two*ll1+one)
726 : end do
727 4 : do ll=1,lexexch+1
728 3 : ll1=2*ll-2
729 18 : if (int_meshsz<mesh_size) ff(int_meshsz+1:mesh_size)=zero
730 3228 : ff(1:int_meshsz)=pawtab(itypat)%phiphj(1:int_meshsz,kln)
731 3 : call poisson(ff,ll1,pawrad(itypat),gg)
732 3 : ff(1)=zero
733 : ff(2:mesh_size)=(pawtab(itypat)%phiphj(2:mesh_size,kln1)*gg(2:mesh_size))&
734 3240 : & /pawrad(itypat)%rad(2:mesh_size)
735 3 : call simp_gen(intg,ff,pawrad(itypat))
736 4 : pawtab(itypat)%fk(4,ll)=intg*(two*ll1+one)
737 : end do
738 4 : do ll=1,lexexch+1
739 3 : ll1=2*ll-2
740 18 : if (int_meshsz<mesh_size) ff(int_meshsz+1:mesh_size)=zero
741 3228 : ff(1:int_meshsz)=pawtab(itypat)%phiphj(1:int_meshsz,kln)
742 3 : call poisson(ff,ll1,pawrad(itypat),gg)
743 3 : ff(1)=zero
744 : ff(2:mesh_size)=(pawtab(itypat)%phiphj(2:mesh_size,kln2)*gg(2:mesh_size))&
745 3240 : & /pawrad(itypat)%rad(2:mesh_size)
746 3 : call simp_gen(intg,ff,pawrad(itypat))
747 4 : pawtab(itypat)%fk(5,ll)=intg*(two*ll1+one)
748 : end do
749 4 : do ll=1,lexexch+1
750 3 : ll1=2*ll-2
751 18 : if (int_meshsz<mesh_size) ff(int_meshsz+1:mesh_size)=zero
752 3228 : ff(1:int_meshsz)=pawtab(itypat)%phiphj(1:int_meshsz,kln1)
753 3 : call poisson(ff,ll1,pawrad(itypat),gg)
754 3 : ff(1)=zero
755 : ff(2:mesh_size)=(pawtab(itypat)%phiphj(2:mesh_size,kln2)*gg(2:mesh_size))&
756 3240 : & /pawrad(itypat)%rad(2:mesh_size)
757 3 : call simp_gen(intg,ff,pawrad(itypat))
758 4 : pawtab(itypat)%fk(6,ll)=intg*(two*ll1+one)
759 : end do
760 1 : f4of2=0.6681_dp
761 1 : f6of2=0.4943_dp
762 : end if
763 1 : ABI_FREE(ff)
764 1 : ABI_FREE(gg)
765 :
766 : ! b. Compute vex.
767 : ! ---------------------------------------------
768 1 : if (allocated(pawtab(itypat)%vex)) then
769 0 : ABI_FREE(pawtab(itypat)%vex)
770 : end if
771 1 : sz1=2*lexexch+1
772 7 : ABI_MALLOC(pawtab(itypat)%vex,(sz1,sz1,sz1,sz1,4))
773 3125 : pawtab(itypat)%vex=zero
774 1 : lmexexch=(lexexch-1)**2+2*(lexexch-1)+1 ! number of m value below correlated orbitals
775 1 : klm0x=lmexexch*(lmexexch+1)/2 ! value of klmn just below correlated orbitals
776 : ! --------- 4 loops for interaction matrix
777 6 : do m1=-lexexch,lexexch
778 5 : m11=m1+lexexch+1
779 21 : do m2=-lexexch,m1
780 15 : m21=m2+lexexch+1
781 : ! klma= number of pair before correlated orbitals +
782 : ! number of pair for m1 lower than correlated orbitals
783 : ! (m1+lexexch+1)*(lexexch-1) + number of pairs for correlated orbitals
784 : ! before (m1,m2) + number of pair for m2 lower than current value
785 15 : klma=klm0x+m11*lmexexch+(m11-1)*m11/2+m21
786 95 : do m3=-lexexch,lexexch
787 75 : m31=m3+lexexch+1
788 315 : do m4=-lexexch,m3
789 225 : m41=m4+lexexch+1
790 225 : klmb=klm0x+m31*lmexexch+(m31-1)*m31/2+m41
791 : ! --------- loop on k=1,2,3 (4 if f orbitals)
792 225 : do kyc=1,2*lexexch+1,2
793 675 : lkyc=kyc-1
794 675 : ll=(kyc+1)/2
795 675 : lmkyc=(lkyc+1)*(lkyc)+1
796 675 : ak=zero
797 4050 : do mkyc=-lkyc,lkyc,1
798 3375 : isela=pawang%gntselect(lmkyc+mkyc,klma)
799 3375 : iselb=pawang%gntselect(lmkyc+mkyc,klmb)
800 4050 : if (isela>0.and.iselb>0) ak=ak +pawang%realgnt(isela)*pawang%realgnt(iselb)
801 : end do
802 : ! ----- end loop on k=1,2,3 (4 if f orbitals)
803 675 : pawtab(itypat)%vex(m11,m31,m21,m41,ll)=ak/(two*dble(lkyc)+one)
804 : end do !kyc
805 1200 : do ll=1,4
806 900 : pawtab(itypat)%vex(m11,m31,m21,m41,ll)=pawtab(itypat)%vex(m11,m31,m21,m41,ll)*four_pi
807 900 : pawtab(itypat)%vex(m21,m31,m11,m41,ll)=pawtab(itypat)%vex(m11,m31,m21,m41,ll)
808 900 : pawtab(itypat)%vex(m11,m41,m21,m31,ll)=pawtab(itypat)%vex(m11,m31,m21,m41,ll)
809 1125 : pawtab(itypat)%vex(m21,m41,m11,m31,ll)=pawtab(itypat)%vex(m11,m31,m21,m41,ll)
810 : end do
811 : end do
812 : end do
813 : end do
814 : end do
815 :
816 : end if !useexexch/=0
817 :
818 301 : if (present(ucrpa)) then
819 292 : if (ucrpa>=1) then
820 0 : call calc_ubare(itypat,lcur,pawang,pawrad(itypat),pawtab(itypat))
821 0 : call calc_ubare(itypat,lcur,pawang,pawrad(itypat),pawtab(itypat),pawtab(itypat)%rpaw)
822 : end if
823 : end if
824 :
825 301 : if (use_dmft > 0) then
826 :
827 52 : write(tag,'(i4)') itypat
828 : write(message,'(4a)') &
829 52 : & ch10,' ===== Build DMFT radial orbital for atom type ',trim(adjustl(tag)),' ========'
830 52 : call wrtout(std_out,message,"COLL")
831 :
832 52 : ABI_SFREE(pawtab(itypat)%proj)
833 52 : ABI_SFREE(pawtab(itypat)%proj2)
834 :
835 52 : me = xmpi_comm_rank(xmpi_world)
836 :
837 52 : if (dmft_orbital(itypat) > 0) then ! use atomic orbital from PAW dataset
838 52 : if (dmft_orbital(itypat) > pawtab(itypat)%nproju) then
839 0 : write(message,*) "For atom type:",itypat,"you need to set dmft_orbital to a value", &
840 0 : & " lower than",pawtab(itypat)%nproju
841 0 : ABI_ERROR(message)
842 : end if
843 52 : write(message,'(2a,i1,a)') ch10," Using atomic orbital number ",dmft_orbital(itypat)," from PAW dataset"
844 52 : call wrtout(std_out,message,"COLL")
845 52 : meshsz = pawrad(itypat)%int_meshsz
846 156 : ABI_MALLOC(pawtab(itypat)%proj,(meshsz))
847 40827 : pawtab(itypat)%proj(:) = pawtab(itypat)%phi(1:meshsz,pawtab(itypat)%lnproju(dmft_orbital(itypat)))
848 : else ! read orbital from file
849 0 : call int2char4(itypat,tag2)
850 0 : tmpfil = trim(adjustl(dmft_orbital_filepath)) // '_' // tag2
851 0 : write(message,'(3a)') ch10," Using wavefunction from file ",trim(tmpfil)
852 0 : call wrtout(std_out,message,"COLL")
853 0 : inquire(file=trim(tmpfil),exist=lexist)
854 0 : if (.not. lexist) ABI_ERROR("File "//trim(tmpfil)//" does not exist !")
855 0 : if (me == 0) then
856 0 : unt = get_unit()
857 0 : open(unit=unt,file=trim(tmpfil),status='unknown',form='formatted')
858 0 : read(unt,*,iostat=ierr) meshsz
859 : end if ! me=0
860 0 : call xmpi_bcast(meshsz,0,xmpi_world,ierr)
861 0 : ABI_MALLOC(pawtab(itypat)%proj,(meshsz))
862 0 : if (me == 0) then
863 0 : do ir=1,meshsz
864 0 : read(unt,*,iostat=ierr) pawtab(itypat)%proj(ir)
865 : end do ! ir
866 0 : close(unt)
867 : end if ! me=0
868 0 : call xmpi_bcast(ierr,0,xmpi_world,ir)
869 0 : if (ierr /= 0) ABI_ERROR("Error when reading file "//trim(tmpfil))
870 0 : call xmpi_bcast(pawtab(itypat)%proj(:),0,xmpi_world,ierr)
871 : end if ! dmft_orbital
872 :
873 52 : mesh_type = pawrad(itypat)%mesh_type
874 52 : lstep = pawrad(itypat)%lstep
875 52 : rstep = pawrad(itypat)%rstep
876 :
877 52 : call pawrad_init(pawrad_tmp,meshsz,mesh_type,rstep,lstep)
878 40827 : call simp_gen(int1,pawtab(itypat)%proj(1:meshsz)**2,pawrad_tmp)
879 :
880 52 : write(message,'(a,f6.4)') " Squared norm of the DMFT orbital: ",int1
881 52 : call wrtout(std_out,message,"COLL")
882 :
883 52 : int1 = sqrt(int1)
884 :
885 52 : if (dmft_dc == 8) then
886 :
887 0 : if (dmft_orbital(itypat) > 0) then
888 0 : write(message ,'(7a)') "WARNING: You are using dmft_dc=8 while using an atomic orbital from ", &
889 0 : & "the PAW dataset. In our current implementation, we assume that ", &
890 0 : & "the projection of the orbital on [dmftbandi,dmftbandf] is the same ", &
891 0 : & "as the orbital itself, and this can hardly be the case with a truncated atomic ", &
892 0 : & "orbital. Please compute the projection of the atomic orbital with dmft_prtwan=1, ", &
893 0 : & "and then use this projection as your DMFT orbital with dmft_orbital=-1. This is explained ",&
894 0 : & "in the tutorial."
895 0 : ABI_WARNING(message)
896 : end if
897 :
898 0 : ABI_MALLOC(pawtab(itypat)%proj2,(meshsz))
899 0 : pawtab(itypat)%proj2(:) = (pawtab(itypat)%proj(:)/int1)**2
900 :
901 0 : ABI_MALLOC(fk,(lcur+1))
902 :
903 0 : if (dmft_yukawa_param <= 2) then
904 : ! Get correspondence U,J <-> lambda,epsilon or U <-> lambda depending on the value of dmft_yukawa_param
905 : call get_lambda(lcur,pawrad_tmp,pawtab(itypat)%proj2(:),meshsz, &
906 0 : & pawtab(itypat)%upawu,pawtab(itypat)%jpawu,lambda,eps,dmft_yukawa_param)
907 0 : else if (dmft_yukawa_param == 3) then
908 0 : call compute_slater(lcur,pawrad_tmp,pawtab(itypat)%proj2(:),meshsz,zero,one,fk(:))
909 0 : lambda = zero
910 0 : eps = fk(1) / pawtab(itypat)%upawu
911 0 : fk(:) = fk(:) / eps
912 0 : else if (dmft_yukawa_param == 4) then
913 0 : lambda = dmft_yukawa_lambda
914 0 : eps = dmft_yukawa_epsilon
915 : end if
916 :
917 0 : pawtab(itypat)%lambda = lambda
918 0 : pawtab(itypat)%eps = eps
919 :
920 : ! Recompute Slater integrals
921 0 : if (dmft_yukawa_param /= 3) then
922 0 : call compute_slater(lcur,pawrad_tmp,pawtab(itypat)%proj2(:),meshsz,lambda,eps,fk(:))
923 : end if
924 :
925 0 : write(message,'(3a)') ch10," Yukawa parameters for atom type: ",adjustl(tag)
926 0 : call wrtout(std_out,message,"COLL")
927 0 : write(message,'(a,f9.4)') " Lambda: ",lambda
928 0 : call wrtout(std_out,message,"COLL")
929 0 : write(message,'(a,f9.4)') " Epsilon:",eps
930 0 : call wrtout(std_out,message,"COLL")
931 :
932 : ! Recompute U tensor with new Slater integrals
933 0 : f4of2 = - one
934 0 : f6of2 = - one
935 0 : uh = fk(1)
936 :
937 0 : if (lcur == 0) then
938 0 : jh = zero
939 0 : else if (lcur == 1) then
940 0 : jh = fk(2) / dble(5.)
941 0 : else if (lcur == 2) then
942 0 : f4of2 = fk(3) / fk(2)
943 0 : jh = fk(2) * (one + f4of2) / dble(14.)
944 0 : else if (lcur == 3) then
945 0 : f6of2 = fk(4) / fk(2)
946 0 : f4of2 = fk(3) / fk(2)
947 0 : jh = fk(2) * (dble(286.)+dble(195.)*f4of2+dble(250.)*f6of2) / dble(6435.)
948 : else
949 0 : write(message,'(a,i0,2a)') ' lpawu=',lpawu,ch10, &
950 0 : & ' lpawu not equal to 0, 1, 2 or 3 is not allowed'
951 0 : ABI_ERROR(message)
952 : end if
953 :
954 0 : write(message,'(a,6x,f9.4)') " U:",uh
955 0 : call wrtout(std_out,message,"COLL")
956 0 : write(message,'(a,6x,f9.4)') " J:",jh
957 0 : call wrtout(std_out,message,"COLL")
958 :
959 0 : pawtab(itypat)%upawu = uh
960 0 : pawtab(itypat)%jpawu = jh
961 :
962 0 : call calc_vee(f4of2,f6of2,jh,lcur,pawang,uh,pawtab(itypat)%vee(:,:,:,:),Loc_prtvol)
963 :
964 0 : ABI_FREE(fk)
965 :
966 : end if ! dmft_dc=8
967 :
968 52 : call pawrad_free(pawrad_tmp)
969 :
970 : end if ! use_dmft > 0
971 :
972 : end if !lcur/=-1
973 : end do !end loop on typat
974 :
975 : DBG_EXIT("COLL")
976 :
977 1729 : end subroutine pawpuxinit
978 : !!***
979 :
980 : !----------------------------------------------------------------------
981 :
982 : !!****f* m_paw_correlations/calc_vee
983 : !! NAME
984 : !! calc_vee
985 : !!
986 : !! FUNCTION
987 : !!
988 : !! Compute matrix elements of coulomb interaction (see PRB vol.52 5467) [[cite:Liechenstein1995]]
989 : !! (angular part computed from Gaunt coefficients)
990 : !!
991 : !! INPUTS
992 : !1 f4of2_sla= Ratio of Slater integrals.
993 : !1 f6of2_sla= Ratio of Slater integrals.
994 : !! jpawu= value of J
995 : !! lpawu= value of l on which DFT+U applies
996 : !! upawu= value of U
997 : !!
998 : !! OUTPUT
999 : !! vee(2*lpawu+1,:,:,:)=matrix of the screened interaction for correlated orbitals
1000 : !!
1001 : !! SOURCE
1002 :
1003 300 : subroutine calc_vee(f4of2_sla,f6of2_sla,jpawu,lpawu,pawang,upawu,vee,prtvol)
1004 :
1005 : !Arguments ---------------------------------------------
1006 : !scalars
1007 : integer,intent(in) :: lpawu
1008 : integer,optional,intent(in) :: prtvol
1009 : real(dp),intent(in) :: upawu,jpawu
1010 : real(dp),intent(inout) :: f4of2_sla,f6of2_sla
1011 : type(pawang_type), intent(in) :: pawang
1012 : !arrays
1013 : real(dp),intent(out) :: vee(2*lpawu+1,2*lpawu+1,2*lpawu+1,2*lpawu+1)
1014 :
1015 : !Local variables ---------------------------------------
1016 : !scalars
1017 : integer :: isela,iselb
1018 : integer :: klm0u,klma,klmb,kyc,lkyc
1019 : integer :: lmkyc,lmpawu
1020 : integer :: m1,m11,m2,m21,m3,m31,m4,m41,prtvol_
1021 : integer :: mkyc,sz1
1022 : real(dp) :: ak,f4of2,f6of2
1023 : character(len=500) :: message
1024 : !arrays
1025 300 : real(dp),allocatable :: fk(:)
1026 :
1027 : ! *************************************************************************
1028 :
1029 : DBG_ENTER("COLL")
1030 :
1031 :
1032 300 : prtvol_ = 3
1033 300 : if (present(prtvol)) then
1034 300 : prtvol_ = prtvol
1035 : end if
1036 : ! Select only atoms with +U
1037 300 : if(lpawu/=-1) then
1038 :
1039 : ! ======================================================================
1040 : ! C-PAW+U: Matrix elements of coulomb interaction (see PRB vol.52 5467) [[cite:Liechenstein1995]]
1041 : ! 1. angular part computed from Gaunt coefficients
1042 : ! --------------------------------------------------------------------
1043 : ! a. compute F(k)
1044 : ! ---------------------------------------------
1045 900 : ABI_MALLOC(fk,(lpawu+1))
1046 300 : fk(1)=upawu
1047 : ! cf Slater Physical Review 165, p 665 (1968) [[cite:Slater1958]]
1048 : ! write(std_out,*) "f4of2_sla",pawtab(itypat)%f4of2_sla
1049 300 : if(lpawu==0) then
1050 : fk(1)=fk(1)
1051 300 : else if(lpawu==1) then
1052 114 : fk(2)=jpawu*5._dp
1053 186 : else if(lpawu==2) then
1054 : ! f4of2=0._dp
1055 163 : if(f4of2_sla<-0.1_dp) then
1056 145 : f4of2=0.625_dp
1057 145 : f4of2_sla=f4of2
1058 : else
1059 : f4of2=f4of2_sla
1060 : end if
1061 163 : fk(2)=jpawu*14._dp/(One+f4of2)
1062 163 : fk(3)=fk(2)*f4of2
1063 163 : if(abs(prtvol_)>=2) then
1064 145 : write(message,'(a,3x,a,f9.4,f9.4,f9.4,f9.4)') ch10,&
1065 290 : & "Slater parameters F^0, F^2, F^4 are",fk(1),fk(2),fk(3)
1066 145 : call wrtout(std_out,message,'COLL')
1067 : end if
1068 23 : else if(lpawu==3) then
1069 23 : f4of2=0.6681_dp
1070 23 : f6of2=0.4943_dp
1071 23 : if(f4of2_sla<-0.1_dp) then
1072 21 : f4of2=0.6681_dp
1073 21 : f4of2_sla=f4of2
1074 : else
1075 : f4of2=f4of2_sla
1076 : end if
1077 23 : if(f6of2_sla<-0.1_dp) then
1078 21 : f6of2=0.4943_dp
1079 21 : f6of2_sla=f6of2
1080 : else
1081 : f6of2=f6of2_sla
1082 : end if
1083 23 : fk(2)=jpawu*6435._dp/(286._dp+195._dp*f4of2+250._dp*f6of2)
1084 23 : fk(3)=fk(2)*f4of2
1085 23 : fk(4)=fk(2)*f6of2
1086 23 : if(abs(prtvol_)>=2) then
1087 23 : write(std_out,'(a,3x,a,f9.4,f9.4,f9.4,f9.4)') ch10,&
1088 46 : & "Slater parameters F^0, F^2, F^4, F^6 are",fk(1),fk(2),fk(3),fk(4)
1089 : end if
1090 : else
1091 : write(message, '(a,i0,2a)' )&
1092 0 : & ' lpawu=',lpawu,ch10,&
1093 0 : & ' lpawu not equal to 0 ,1 ,2 or 3 is not allowed'
1094 0 : ABI_ERROR(message)
1095 : end if
1096 :
1097 : ! b. Compute ak and vee.
1098 : ! ---------------------------------------------
1099 : ! if (allocated(vee)) then
1100 : ! ABI_DEALLOCATE(vee)
1101 : ! end if
1102 : sz1=2*lpawu+1
1103 : ! ABI_ALLOCATE(vee,(sz1,sz1,sz1,sz1))
1104 205520 : vee=zero
1105 300 : lmpawu=(lpawu-1)**2+2*(lpawu-1)+1 ! number of m value below correlated orbitals
1106 300 : klm0u=lmpawu*(lmpawu+1)/2 ! value of klmn just below correlated orbitals
1107 : ! --------- 4 loops for interaction matrix
1108 1618 : do m1=-lpawu,lpawu
1109 1318 : m11=m1+lpawu+1
1110 5391 : do m2=-lpawu,m1
1111 3773 : m21=m2+lpawu+1
1112 : ! klma= number of pair before correlated orbitals +
1113 : ! number of pair for m1 lower than correlated orbitals
1114 : ! (m1+lpawu+1)*(lpawu-1) + number of pairs for correlated orbitals
1115 : ! before (m1,m2) + number of pair for m2 lower than current value
1116 3773 : klma=klm0u+m11*lmpawu+(m11-1)*m11/2+m21
1117 23876 : do m3=-lpawu,lpawu
1118 18785 : m31=m3+lpawu+1
1119 81369 : do m4=-lpawu,m3
1120 58811 : m41=m4+lpawu+1
1121 58811 : klmb=klm0u+m31*lmpawu+(m31-1)*m31/2+m41
1122 : ! --------- loop on k=1,2,3 (4 if f orbitals)
1123 58811 : do kyc=1,2*lpawu+1,2
1124 190361 : lkyc=kyc-1
1125 190361 : lmkyc=(lkyc+1)*(lkyc)+1
1126 190361 : ak=zero
1127 1270006 : do mkyc=-lkyc,lkyc,1
1128 1079645 : isela=pawang%gntselect(lmkyc+mkyc,klma)
1129 1079645 : iselb=pawang%gntselect(lmkyc+mkyc,klmb)
1130 1270006 : if (isela>0.and.iselb>0) ak=ak +pawang%realgnt(isela)*pawang%realgnt(iselb)
1131 : end do
1132 : ! ----- end loop on k=1,2,3 (4 if f orbitals)
1133 190361 : ak=ak/(two*dble(lkyc)+one)
1134 190361 : vee(m11,m31,m21,m41)=ak*fk(lkyc/2+1)+vee(m11,m31,m21,m41)
1135 : end do !kyc
1136 58811 : vee(m11,m31,m21,m41)=vee(m11,m31,m21,m41)*four_pi
1137 58811 : vee(m21,m31,m11,m41)=vee(m11,m31,m21,m41)
1138 58811 : vee(m11,m41,m21,m31)=vee(m11,m31,m21,m41)
1139 77596 : vee(m21,m41,m11,m31)=vee(m11,m31,m21,m41)
1140 : end do
1141 : end do
1142 : end do
1143 : end do
1144 300 : ABI_FREE(fk)
1145 : endif
1146 :
1147 300 : end subroutine calc_vee
1148 : !!***
1149 :
1150 : !----------------------------------------------------------------------
1151 : !!****f* m_paw_correlations/pawuenergy
1152 : !! NAME
1153 : !! pawuenergy
1154 : !!
1155 : !! FUNCTION
1156 : !! Compute contributions to energy for PAW+U calculations
1157 : !!
1158 : !! INPUTS
1159 : !! iatom=index of current atom (absolute index, the index on current proc)
1160 : !! noccmmp(2*lpawu+1,2*lpawu+1,nspden)=density matrix in the PAW augm. region
1161 : !! nocctot(nspden)=number of electrons in the correlated subspace
1162 : !! pawprtvol=control print volume and debugging output for PAW
1163 : !! pawtab <type(pawtab_type)>=paw tabulated starting data:
1164 : !! %lpawu=l used for dft+u
1165 : !! %vee(2*lpawu+1*4)=screened coulomb matrix
1166 : !! dmft_dc,e_ee,e_dc,e_dcdc,u_dmft,j_dmft= optional arguments for DMFT
1167 : !!
1168 : !! OUTPUT
1169 : !! edftumdc= PAW+U contribution to total energy
1170 : !! edftumdcdc= PAW+U contribution to double-counting total energy
1171 : !!
1172 : !! SOURCE
1173 :
1174 2707 : subroutine pawuenergy(iatom,edftumdc,edftumdcdc,noccmmp,nocctot,pawprtvol,pawtab,&
1175 : & dmft_dc,e_ee,e_dc,e_dcdc,u_dmft,j_dmft,paw_dmft) ! optional arguments (DMFT)
1176 :
1177 : !Arguments ---------------------------------------------
1178 : !scalars
1179 : integer,intent(in) :: iatom,pawprtvol
1180 : integer,optional,intent(in) :: dmft_dc
1181 : real(dp),intent(in) :: noccmmp(:,:,:,:),nocctot(:)
1182 : real(dp),intent(inout) :: edftumdc,edftumdcdc
1183 : real(dp),optional,intent(inout) :: e_ee,e_dc,e_dcdc
1184 : real(dp),optional,intent(in) :: j_dmft,u_dmft
1185 : type(paw_dmft_type),optional,intent(in) :: paw_dmft
1186 : type(pawtab_type),intent(in) :: pawtab
1187 : !Local variables ---------------------------------------
1188 : !scalars
1189 : integer :: cplex_occ,dmftdc,ispden,jspden,lpawu,m1,m11,m2,m21,m3,m31,m4,m41,nspden
1190 : real(dp) :: eks_opt3,edcdc_opt3,edcdctemp,edctemp,edftutemp,jpawu,jpawu_dc,mnorm,mx,my,mz
1191 : real(dp) :: n_sig,n_sigs,n_msig,n_msigs,n_dndn,n_tot,n_upup
1192 : real(dp) :: n12_ud_im,n12_du_im
1193 : real(dp) :: n12_ud_re,n12_du_re
1194 : real(dp) :: n34_ud_im,n34_du_im
1195 : real(dp) :: n34_ud_re,n34_du_re
1196 : real(dp) :: upawu
1197 2707 : real(dp),allocatable :: n12_sig(:),n34_msig(:),n34_sig(:)
1198 : character(len=4) :: tag
1199 : character(len=500) :: message
1200 : ! *****************************************************
1201 :
1202 2707 : nspden=size(nocctot)
1203 2707 : cplex_occ=size(noccmmp,1)
1204 :
1205 2707 : if (size(noccmmp,4)/=nspden) then
1206 0 : message='size of nocctot and noccmmp are inconsistent!'
1207 0 : ABI_BUG(message)
1208 : end if
1209 2707 : if (pawtab%usepawu<0) then
1210 0 : message='not allowed for usepawu<0!'
1211 0 : ABI_BUG(message)
1212 : end if
1213 2707 : if(present(dmft_dc)) then
1214 244 : dmftdc=dmft_dc
1215 244 : if(pawtab%usepawu<10) then
1216 0 : write(message,'(a,i5)') "usepawu should be =10 if dmft_dc is present ",pawtab%usepawu
1217 0 : ABI_BUG(message)
1218 : end if
1219 : else
1220 : dmftdc=0
1221 : end if
1222 :
1223 : DBG_ENTER("COLL")
1224 :
1225 2707 : lpawu=pawtab%lpawu
1226 2707 : upawu=pawtab%upawu;if(present(u_dmft)) upawu=u_dmft
1227 2707 : jpawu=pawtab%jpawu;if(present(j_dmft)) jpawu=j_dmft
1228 :
1229 : !======================================================
1230 : !Compute DFT+U Energy
1231 : !-----------------------------------------------------
1232 :
1233 2707 : edftutemp=zero
1234 2707 : edcdc_opt3=zero
1235 2707 : eks_opt3=zero
1236 :
1237 8121 : ABI_MALLOC(n12_sig,(cplex_occ))
1238 5414 : ABI_MALLOC(n34_msig,(cplex_occ))
1239 5414 : ABI_MALLOC(n34_sig,(cplex_occ))
1240 7771 : do ispden=1,min(nspden,2)
1241 5064 : jspden=min(nspden,2)-ispden+1
1242 :
1243 : ! Compute n_sigs and n_msigs for pawtab%usepawu=3
1244 5064 : if (nspden<=2) then
1245 3794 : n_sig =nocctot(ispden)
1246 3794 : n_msig=nocctot(jspden)
1247 3794 : n_tot=n_sig+n_msig
1248 : else
1249 1270 : n_tot=nocctot(1)
1250 1270 : mx=nocctot(2)
1251 1270 : my=nocctot(3)
1252 1270 : mz=nocctot(4)
1253 1270 : mnorm=sqrt(mx*mx+my*my+mz*mz)
1254 1270 : if (ispden==1) then
1255 : ! n_sig =half*(n_tot+mnorm)
1256 : ! n_msig=half*(n_tot-mnorm)
1257 635 : n_sig =half*(n_tot+sign(mnorm,mz))
1258 635 : n_msig=half*(n_tot-sign(mnorm,mz))
1259 : else
1260 : ! n_sig =half*(n_tot-mnorm)
1261 : ! n_msig=half*(n_tot+mnorm)
1262 635 : n_sig =half*(n_tot-sign(mnorm,mz))
1263 635 : n_msig=half*(n_tot+sign(mnorm,mz))
1264 : end if
1265 : end if
1266 5064 : n_sigs =n_sig/(float(2*lpawu+1))
1267 5064 : n_msigs =n_msig/(float(2*lpawu+1))
1268 : ! if(pawtab%usepawu==3) then
1269 : ! write(message,fmt=12) "noccmmp11 ",ispden,noccmmp(1,1,1,ispden)
1270 : ! call wrtout(std_out,message,'COLL')
1271 : ! write(message,fmt=12) "noccmmp11 ",jspden,noccmmp(1,1,1,jspden)
1272 : ! call wrtout(std_out,message,'COLL')
1273 : ! write(message,fmt=12) "n_sig ",ispden,n_sig
1274 : ! call wrtout(std_out,message,'COLL')
1275 : ! write(message,fmt=12) "n_msig ",jspden,n_msig
1276 : ! call wrtout(std_out,message,'COLL')
1277 : ! write(message,fmt=12) "n_sigs ",ispden,n_sigs
1278 : ! call wrtout(std_out,message,'COLL')
1279 : ! write(message,fmt=12) "n_msigs ",jspden,n_msigs
1280 : ! call wrtout(std_out,message,'COLL')
1281 : ! endif
1282 : ! 12 format(a,i4,e20.10)
1283 :
1284 : ! Compute interaction energy E_{ee}
1285 33077 : do m1=-lpawu,lpawu
1286 25306 : m11=m1+lpawu+1
1287 160386 : do m2=-lpawu,lpawu
1288 130016 : m21=m2+lpawu+1
1289 303481 : n12_sig(:)=noccmmp(:,m11,m21,ispden)
1290 130016 : if(m21==m11.and.(pawtab%usepawu==3.or.dmftdc==3)) n12_sig(1)=n12_sig(1)-n_sigs
1291 840556 : do m3=-lpawu,lpawu
1292 685234 : m31=m3+lpawu+1
1293 4519754 : do m4=-lpawu,lpawu
1294 3704504 : m41=m4+lpawu+1
1295 8614009 : n34_sig(:) =noccmmp(:,m31,m41,ispden)
1296 8614009 : n34_msig(:)=noccmmp(:,m31,m41,jspden)
1297 3704504 : if(m31==m41.and.(pawtab%usepawu==3.or.dmftdc==3)) then
1298 0 : n34_sig(1)= n34_sig(1) - n_sigs
1299 0 : n34_msig(1)= n34_msig(1) - n_msigs
1300 : end if
1301 : edftutemp=edftutemp &
1302 : & + n12_sig(1)*n34_msig(1)*pawtab%vee(m11,m31,m21,m41) &
1303 3704504 : & + n12_sig(1)*n34_sig(1) *(pawtab%vee(m11,m31,m21,m41)-pawtab%vee(m11,m31,m41,m21))
1304 3704504 : if(cplex_occ==2) then
1305 : edftutemp=edftutemp &
1306 : & - n12_sig(2)*n34_msig(2)*pawtab%vee(m11,m31,m21,m41) &
1307 1205001 : & - n12_sig(2)*n34_sig(2) *(pawtab%vee(m11,m31,m21,m41)-pawtab%vee(m11,m31,m41,m21))
1308 : end if
1309 4389738 : if (pawtab%usepawu==3.or.dmftdc==3) then
1310 : edcdc_opt3=edcdc_opt3 &
1311 : & + n_sigs*n34_msig(1)*pawtab%vee(m11,m31,m21,m41) &
1312 0 : & + n_sigs*n34_sig(1) *(pawtab%vee(m11,m31,m21,m41)-pawtab%vee(m11,m31,m41,m21))
1313 : eks_opt3=eks_opt3 &
1314 : & + noccmmp(1,m11,m21,ispden)*n34_msig(1)*pawtab%vee(m11,m31,m21,m41) &
1315 0 : & + noccmmp(1,m11,m21,ispden)*n34_sig(1) *(pawtab%vee(m11,m31,m21,m41)-pawtab%vee(m11,m31,m41,m21))
1316 0 : if(cplex_occ==2) then
1317 : eks_opt3=eks_opt3 &
1318 : & - noccmmp(2,m11,m21,ispden)*n34_msig(2)*pawtab%vee(m11,m31,m21,m41) &
1319 0 : & - noccmmp(2,m11,m21,ispden)*n34_sig(2) *(pawtab%vee(m11,m31,m21,m41)-pawtab%vee(m11,m31,m41,m21))
1320 : end if
1321 : end if
1322 : end do ! m4
1323 : end do ! m3
1324 : end do ! m2
1325 : end do ! m1
1326 :
1327 : end do ! ispden
1328 2707 : if (nspden==1) edftutemp=two*edftutemp ! Non-magn. system: sum up and dn energies
1329 2707 : ABI_FREE(n12_sig)
1330 2707 : ABI_FREE(n34_msig)
1331 2707 : ABI_FREE(n34_sig)
1332 :
1333 : !Non-collinear magnetism: add non-diagonal term; see (Eq 3) in PRB 72, 024458 (2005) [[cite:Shurikov2005]]
1334 2707 : if (nspden==4) then
1335 3872 : do m1=-lpawu,lpawu
1336 3237 : m11=m1+lpawu+1
1337 20491 : do m2=-lpawu,lpawu
1338 16619 : m21=m2+lpawu+1
1339 16619 : n12_ud_re=noccmmp(1,m11,m21,3) ! updn
1340 16619 : n12_ud_im=noccmmp(2,m11,m21,3) ! updn
1341 16619 : n12_du_re=noccmmp(1,m11,m21,4) ! dnup
1342 16619 : n12_du_im=noccmmp(2,m11,m21,4) ! dnup
1343 105989 : do m3=-lpawu,lpawu
1344 86133 : m31=m3+lpawu+1
1345 554683 : do m4=-lpawu,lpawu
1346 451931 : m41=m4+lpawu+1
1347 451931 : n34_ud_re=noccmmp(1,m31,m41,3) ! updn
1348 451931 : n34_ud_im=noccmmp(2,m31,m41,3) ! updn
1349 451931 : n34_du_re=noccmmp(1,m31,m41,4) ! dnup
1350 451931 : n34_du_im=noccmmp(2,m31,m41,4) ! dnup
1351 : edftutemp=edftutemp-pawtab%vee(m11,m31,m41,m21) &
1352 : & *(n12_ud_re*n34_du_re-n12_ud_im*n34_du_im &
1353 451931 : & +n12_du_re*n34_ud_re-n12_du_im*n34_ud_im)
1354 538064 : if (pawtab%usepawu==3.or.dmftdc==3) then
1355 : eks_opt3=eks_opt3-pawtab%vee(m11,m31,m41,m21) &
1356 : & *(n12_ud_re*n34_du_re-n12_ud_im*n34_du_im &
1357 0 : & +n12_du_re*n34_ud_re-n12_du_im*n34_ud_im)
1358 : end if
1359 : end do ! m4
1360 : end do ! m3
1361 : end do ! m2
1362 : end do ! m1
1363 : end if
1364 :
1365 : !Divide edftutemp by 2; see (Eq 1) in PRB 77, 155104 (2008) [[cite:Amadon2008a]]
1366 2707 : edftutemp=half*edftutemp
1367 :
1368 : !if (nspden==1) then
1369 : !n_tot=two*nocctot(1)
1370 : !n_upup=nocctot(1)
1371 : !n_dndn=nocctot(1)
1372 : !else if (nspden==2) then
1373 : !n_tot=nocctot(1)+nocctot(2)
1374 : !n_upup=nocctot(1)
1375 : !n_dndn=nocctot(2)
1376 : !else if (nspden==4) then
1377 : !n_tot=nocctot(1)
1378 : !mx=nocctot(2)
1379 : !my=nocctot(3)
1380 : !mz=nocctot(4)
1381 : !mnorm=sqrt(mx*mx+my*my+mz*mz)
1382 : !n_upup=half*(n_tot+mnorm)
1383 : !n_dndn=half*(n_tot-mnorm)
1384 : !end if
1385 2707 : n_upup=n_sig
1386 2707 : n_dndn=n_msig
1387 :
1388 2707 : edcdctemp=zero;edctemp=zero
1389 :
1390 : !Full localized limit
1391 2707 : if((pawtab%usepawu==1.or.pawtab%usepawu==4).or.(dmftdc==1.or.dmftdc==4.or.dmftdc==5)) then
1392 2611 : jpawu_dc=jpawu
1393 2611 : if(dmftdc==4) then
1394 2637 : jpawu_dc=zero
1395 : end if
1396 2637 : edcdctemp=edcdctemp-half*upawu*n_tot**2
1397 2637 : edctemp =edctemp +half*upawu*(n_tot*(n_tot-one))
1398 2637 : if (nspden/=4.or.pawtab%option_interaction_pawu==2) then
1399 2002 : if(dmftdc/=5.and.pawtab%usepawu/=4) then
1400 1806 : edcdctemp=edcdctemp+half*jpawu_dc*(n_upup**2+n_dndn**2)
1401 1806 : edctemp =edctemp -half*jpawu_dc*(n_upup*(n_upup-one)+n_dndn*(n_dndn-one))
1402 : else if(dmftdc==5.or.pawtab%usepawu==4) then
1403 196 : edcdctemp=edcdctemp+quarter*jpawu_dc*n_tot**2
1404 196 : edctemp =edctemp -quarter*jpawu_dc*(n_tot*(n_tot-two))
1405 : end if
1406 635 : else if (nspden==4.and.(pawtab%usepawu==4.or.pawtab%option_interaction_pawu==1)) then
1407 : ! write(message,'(a)') " warning: option_interaction==1 for test "
1408 : ! call wrtout(std_out,message,'COLL')
1409 263 : edcdctemp=edcdctemp+quarter*jpawu_dc*n_tot**2
1410 263 : edctemp =edctemp -quarter*jpawu_dc*(n_tot*(n_tot-two))
1411 372 : else if (nspden==4.and.pawtab%option_interaction_pawu==3) then
1412 : ! edcdctemp= \frac{J}/{4}[ N(N) + \vect{m}.\vect{m}]
1413 : edcdctemp=edcdctemp+quarter*jpawu_dc*(n_tot**2 + &
1414 372 : & mx**2+my**2+mz**2) ! +\frac{J}/{4}\vect{m}.\vect{m}
1415 : ! edctemp= -\frac{J}/{4}[ N(N-2) + \vect{m}.\vect{m}]
1416 : edctemp =edctemp -quarter*jpawu_dc*( &
1417 : & (n_tot*(n_tot-two)) + &
1418 372 : & mx**2+my**2+mz**2) ! -\frac{J}/{4}\vect{m}.\vect{m}
1419 : end if
1420 :
1421 : ! Around mean field
1422 70 : else if(pawtab%usepawu==2.or.dmftdc==2) then
1423 : edctemp=edctemp+upawu*(n_upup*n_dndn)&
1424 : & +half*(upawu-jpawu)*(n_upup**2+n_dndn**2) &
1425 70 : & *(dble(2*lpawu)/dble(2*lpawu+1))
1426 70 : edcdctemp=-edctemp
1427 0 : else if(pawtab%usepawu==6.or.dmftdc==6) then
1428 : edctemp=edctemp+upawu*(n_tot*n_tot/4_dp)&
1429 : & +half*(upawu-jpawu)*(n_tot**2+n_tot**2)/4_dp &
1430 0 : & *(dble(2*lpawu)/dble(2*lpawu+1))
1431 0 : edcdctemp=-edctemp
1432 0 : else if(pawtab%usepawu==3.or.dmftdc==3) then
1433 0 : edcdctemp=edcdc_opt3
1434 0 : if(abs(pawprtvol)>=3) then
1435 0 : write(message,fmt=11) "edcdc_opt3 ",edcdc_opt3
1436 0 : call wrtout(std_out,message,'COLL')
1437 0 : write(message,fmt=11) "eks_opt3 ",eks_opt3
1438 0 : call wrtout(std_out,message,'COLL')
1439 0 : write(message,fmt=11) "eks+edcdc_opt3 ",eks_opt3+edcdc_opt3
1440 0 : call wrtout(std_out,message,'COLL')
1441 0 : write(message,fmt=11) "(eks+edcdc_opt3)/2 ",(eks_opt3+edcdc_opt3)/2.d0
1442 0 : call wrtout(std_out,message,'COLL')
1443 : end if
1444 0 : else if(dmftdc==7) then
1445 0 : edctemp=upawu*(dble(paw_dmft%dmft_nominal(iatom))-half)*n_tot-half*jpawu*(dble(paw_dmft%dmft_nominal(iatom))-one)*n_tot
1446 : edcdctemp=zero
1447 0 : else if(dmftdc==8) then
1448 0 : edctemp=paw_dmft%edc(iatom)
1449 0 : edcdctemp=paw_dmft%edc(iatom)-paw_dmft%edcdc(iatom)
1450 : end if
1451 :
1452 2707 : edftumdc =edftumdc +edftutemp-edctemp
1453 2707 : edftumdcdc=edftumdcdc-edftutemp-edcdctemp
1454 :
1455 : !if(pawtab%usepawu/=10.or.pawprtvol>=3) then
1456 2707 : if(abs(pawprtvol)>=3) then
1457 1361 : write(tag,'(i4)') iatom
1458 1361 : tag = trim(adjustl(tag))
1459 1361 : if(pawtab%usepawu<10) then
1460 1126 : write(message, '(6a)') ch10,'======= DFT+U Energy terms (in Hartree) ====',ch10,&
1461 2252 : & ch10,' For Atom ',tag
1462 : else if (pawtab%usepawu >= 10) then
1463 235 : write(message, '(6a)') ch10,' === DFT+U Energy terms from the DMFT occupation matrix ==',ch10,&
1464 470 : & ch10,' For Atom ',tag
1465 : end if
1466 :
1467 1361 : call wrtout(std_out,message,'COLL')
1468 1361 : write(message, '(a)' )" Contributions to the direct expression of energy:"
1469 1361 : call wrtout(std_out, message,'COLL')
1470 1361 : write(message,fmt=11) " Double counting correction =",edctemp
1471 1361 : call wrtout(std_out, message,'COLL')
1472 1361 : write(message,fmt=11) " Interaction energy =",edftutemp
1473 1361 : call wrtout(std_out, message,'COLL')
1474 1361 : write(message,fmt=11) " Total DFT+U Contribution =",edftutemp-edctemp
1475 1361 : call wrtout(std_out, message,'COLL')
1476 1361 : write(message, '(a)' )' '
1477 1361 : call wrtout(std_out, message,'COLL')
1478 1361 : write(message, '(a)' )" For the ""Double-counting"" decomposition:"
1479 1361 : call wrtout(std_out, message,'COLL')
1480 1361 : write(message,fmt=11) " DFT+U Contribution =",-edftutemp-edcdctemp
1481 1361 : call wrtout(std_out, message,'COLL')
1482 : 11 format(a,e20.10)
1483 : if(abs(pawprtvol)>=2) then
1484 1361 : write(message,fmt=11)" edcdctemp =",edcdctemp
1485 1361 : call wrtout(std_out, message,'COLL')
1486 1361 : write(message,fmt=11)" edftumdcdc for current atom =",-edftutemp-edcdctemp
1487 1361 : call wrtout(std_out, message,'COLL')
1488 1361 : write(message, '(a)' )' '
1489 1361 : call wrtout(std_out, message,'COLL')
1490 1361 : write(message,fmt=11)" pawuenergy: -VUKS pred =",edftumdcdc-edftumdc
1491 1361 : call wrtout(std_out, message,'COLL')
1492 : end if
1493 1361 : write(message, '(a)' )' '
1494 1361 : call wrtout(std_out, message,'COLL')
1495 : end if
1496 :
1497 : !For DMFT calculation
1498 2707 : if(present(e_ee)) e_ee=e_ee+edftutemp
1499 2707 : if(present(e_dc)) e_dc=e_dc+edctemp
1500 2707 : if(present(e_dcdc)) e_dcdc=e_dcdc+edcdctemp
1501 :
1502 : DBG_EXIT("COLL")
1503 :
1504 2707 : end subroutine pawuenergy
1505 : !!***
1506 :
1507 : !----------------------------------------------------------------------
1508 :
1509 : !!****f* m_paw_correlations/pawxenergy
1510 : !! NAME
1511 : !! pawxenergy
1512 : !!
1513 : !! FUNCTION
1514 : !! Compute contributions to energy for PAW+ local exact exchange calculations
1515 : !!
1516 : !! INPUTS
1517 : !! pawprtvol=control print volume and debugging output for PAW
1518 : !! pawrhoij <type(pawrhoij_type)>= paw rhoij occupancies and related data
1519 : !! pawtab <type(pawtab_type)>=paw tabulated starting data:
1520 : !! %lexexch=l used for local exact-exchange
1521 : !! %vex(2*lexexch+1*4)=screened coulomb matrix
1522 : !!
1523 : !! SIDE EFFECTS
1524 : !! eexex=energy is updated with the contribution of the cuyrrent atom
1525 : !!
1526 : !! SOURCE
1527 :
1528 42 : subroutine pawxenergy(eexex,pawprtvol,pawrhoij,pawtab)
1529 :
1530 : !Arguments ---------------------------------------------
1531 : !scalars
1532 : integer,intent(in) :: pawprtvol
1533 : real(dp),intent(inout) :: eexex
1534 : type(pawrhoij_type),intent(in) :: pawrhoij
1535 : type(pawtab_type),intent(in) :: pawtab
1536 :
1537 : !Local variables ---------------------------------------
1538 : !scalars
1539 : integer :: irhoij,irhoij1,ispden,jrhoij,jrhoij1,klmn,klmn1,lexexch,ll,m11,m21,m31,m41,n1
1540 : integer :: n2,n3,n4,nk,nn1,nn2
1541 : real(dp) :: eexextemp
1542 : character(len=500) :: message
1543 : !arrays
1544 : integer :: indn(3,3)
1545 : real(dp) :: factnk(6)
1546 :
1547 : ! *****************************************************
1548 :
1549 : DBG_ENTER("COLL")
1550 :
1551 42 : if (pawrhoij%qphase==2) then
1552 0 : message='pawxenergy: local exact-exchange not compatible with qphase=2!'
1553 0 : ABI_ERROR(message)
1554 : end if
1555 :
1556 42 : lexexch=pawtab%lexexch
1557 : if (pawtab%nproju==1) nk=1
1558 42 : if (pawtab%nproju==2) nk=6
1559 42 : factnk(1)=one;factnk(2)=one;factnk(3)=one
1560 42 : factnk(4)=two;factnk(5)=two;factnk(6)=two
1561 42 : indn(1,1)=1;indn(1,2)=4;indn(1,3)=5
1562 42 : indn(2,1)=4;indn(2,2)=2;indn(2,3)=6
1563 42 : indn(3,1)=5;indn(3,2)=6;indn(3,3)=3
1564 :
1565 : !======================================================
1566 : !Compute local exact exchange Energy
1567 : !-----------------------------------------------------
1568 42 : eexextemp=zero
1569 :
1570 126 : do ispden=1,pawrhoij%nspden
1571 84 : jrhoij=1
1572 6790 : do irhoij=1,pawrhoij%nrhoijsel
1573 6664 : klmn=pawrhoij%rhoijselect(irhoij)
1574 6664 : if(pawtab%indklmn(3,klmn)==0.and.pawtab%indklmn(4,klmn)==2*lexexch) then
1575 3780 : m11=pawtab%klmntomn(1,klmn);m21=pawtab%klmntomn(2,klmn)
1576 3780 : n1=pawtab%klmntomn(3,klmn);n2=pawtab%klmntomn(4,klmn)
1577 3780 : nn1=(n1*n2)/2+1
1578 3780 : jrhoij1=1
1579 315980 : do irhoij1=1,pawrhoij%nrhoijsel
1580 312200 : klmn1=pawrhoij%rhoijselect(irhoij1)
1581 312200 : if(pawtab%indklmn(3,klmn1)==0.and.pawtab%indklmn(4,klmn1)==2*lexexch) then
1582 176820 : m31=pawtab%klmntomn(1,klmn1);m41=pawtab%klmntomn(2,klmn1)
1583 176820 : n3=pawtab%klmntomn(3,klmn1);n4=pawtab%klmntomn(4,klmn1)
1584 176820 : nn2=(n3*n4)/2+1
1585 707280 : do ll=1,lexexch+1
1586 : eexextemp=eexextemp-pawtab%vex(m11,m31,m41,m21,ll)*pawtab%dltij(klmn)*pawtab%fk(indn(nn1,nn2),ll)&
1587 707280 : & *pawtab%dltij(klmn1)*pawrhoij%rhoijp(jrhoij,ispden)*pawrhoij%rhoijp(jrhoij1,ispden)
1588 : end do
1589 : end if
1590 315980 : jrhoij1=jrhoij1+pawrhoij%cplex_rhoij
1591 : end do !irhoij1
1592 : end if
1593 6748 : jrhoij=jrhoij+pawrhoij%cplex_rhoij
1594 : end do !irhoij
1595 : end do ! ispden
1596 42 : eexextemp=eexextemp/two
1597 42 : eexex=eexex+eexextemp*pawtab%exchmix
1598 :
1599 42 : if (abs(pawprtvol)>=2) then
1600 0 : write(message, '(a)' )" Contributions to the direct expression of energy:"
1601 0 : call wrtout(std_out,message,'COLL')
1602 0 : write(message,fmt='(a,e20.10,a)') " HF exchange energy =",eexextemp,ch10
1603 0 : call wrtout(std_out,message,'COLL')
1604 : end if
1605 :
1606 : DBG_EXIT("COLL")
1607 :
1608 42 : end subroutine pawxenergy
1609 : !!***
1610 :
1611 : !----------------------------------------------------------------------
1612 :
1613 : !!****f* m_paw_correlations/setnoccmmp
1614 : !! NAME
1615 : !! setnoccmmp
1616 : !!
1617 : !! FUNCTION
1618 : !! PAW+U only:
1619 : !! Compute density matrix nocc_{m,m_prime}
1620 : !! or
1621 : !! Impose value of density matrix using dmatpawu input array, then symetrize it.
1622 : !!
1623 : !! noccmmp^{\sigma}_{m,m'}=\sum_{ni,nj}[\rho^{\sigma}_{ni,nj}*phiphjint_{ni,nj}]
1624 : !!
1625 : !! INPUTS
1626 : !! compute_dmat= flag: if 1, nocc_{m,mp} is computed
1627 : !! dimdmat=first dimension of dmatpawu array
1628 : !! dmatpawu(dimdmat,dimdmat,nsppol*nspinor,natpawu)=input density matrix to be copied into noccmpp
1629 : !! dmatudiag= flag controlling the use of diagonalization:
1630 : !! 0: no diagonalization of nocc_{m,mp}
1631 : !! 1: diagonalized nocc_{m,mp} matrix is printed
1632 : !! 2: dmatpawu matrix is expressed in the basis where nocc_(m,mp} is diagonal
1633 : !! impose_dmat= flag: if 1, nocc_{m,mp} is replaced by dmatpawu
1634 : !! indsym(4,nsym,natom)=indirect indexing array for atom labels
1635 : !! mpi_atmtab(:)=--optional-- indexes of the atoms treated by current proc
1636 : !! comm_atom=--optional-- MPI communicator over atoms
1637 : !! my_natom=number of atoms treated by current processor
1638 : !! natom=number of atoms in cell
1639 : !! natpawu=number of atoms on which PAW+U is applied
1640 : !! nspinor=number of spinorial components of the wavefunctions
1641 : !! nsppol=number of independant spin components
1642 : !! nsym=number of symmetry elements in space group
1643 : !! ntypat=number of atom types
1644 : !! paw_ij(my_natom) <type(paw_ij_type)>=paw arrays given on (i,j) channels
1645 : !! pawang <type(pawang_type)>=paw angular mesh and related data
1646 : !! pawrhoij(my_natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data
1647 : !! pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data
1648 : !! spinat(3,matom)=initial spin of each atom, in unit of hbar/2
1649 : !! symafm(nsym)=(anti)ferromagnetic part of symmetry operations
1650 : !! typat(natom)=type for each atom
1651 : !! useexexch=1 if local-exact-exchange is activated
1652 : !! usepawu= /=0 if PAW+U is activated
1653 : !!
1654 : !! OUTPUT
1655 : !! paw_ij(natom)%noccmmp(cplex_dij,2*lpawu+1,2*lpawu+1,nsppol or ndij)=density matrix
1656 : !!
1657 : !! NOTES
1658 : !! For non-collinear magnetism,
1659 : !! - nocc_{m,mp} is computed as:noccmmp(:,:,:,1)= n{m,mp}
1660 : !! noccmmp(:,:,:,2)= m_x{m,mp}
1661 : !! noccmmp(:,:,:,3)= m_y{m,mp}
1662 : !! noccmmp(:,:,:,4)= m_z{m,mp}
1663 : !! - but nocc_{m,mp} is stored as: noccmmp(:,:,:,1)= n^{up,up}_{m,mp}
1664 : !! noccmmp(:,:,:,2)= n^{dn,dn}_{m,mp}
1665 : !! noccmmp(:,:,:,3)= n^{up,dn}_{m,mp}
1666 : !! noccmmp(:,:,:,4)= n^{dn,up}_{m,mp}
1667 : !! We choose to have noccmmp complex when ndij=4 (ie nspinor=2)
1668 : !! If ndij=4 and pawspnorb=0, one could keep noccmmp real
1669 : !! with the n11, n22, Re(n12), Im(n21) representation, but it would
1670 : !! less clear to change the representation when pawspnorb is activated.
1671 : !! If ndij=4, nocc_{m,mp} is transformed to the Ylm basis
1672 : !! and then to the J, M_J basis (if cplex_dij==2)
1673 : !!
1674 : !! Note that n_{m,mp}=<mp|hat(n)|m> because rhoij=<p_j|...|p_i>
1675 : !!
1676 : !! SOURCE
1677 :
1678 5458 : subroutine setnoccmmp(compute_dmat,dimdmat,dmatpawu,dmatudiag,impose_dmat,indsym,my_natom,natom,&
1679 5458 : & natpawu,nspinor,nsppol,nsym,ntypat,paw_ij,pawang,pawprtvol,pawrhoij,pawtab,&
1680 2729 : & spinat,symafm,typat,useexexch,usepawu, &
1681 2729 : & mpi_atmtab,comm_atom,l_orbmom,atom_orbmom,my_l_occmat) ! optional arguments (parallelism) and printing lorb mag
1682 :
1683 : !Arguments ---------------------------------------------
1684 : !scalars
1685 : integer,intent(in) :: compute_dmat,dimdmat,dmatudiag,impose_dmat,my_natom,natom,natpawu
1686 : integer,intent(in) :: nspinor,nsppol,nsym,ntypat,useexexch,usepawu
1687 : integer,optional,intent(in) :: comm_atom
1688 : type(pawang_type),intent(in) :: pawang
1689 : integer,intent(in) :: pawprtvol
1690 : !arrays
1691 : integer,intent(in) :: indsym(4,nsym,natom),symafm(nsym),typat(natom)
1692 : integer,optional,target,intent(in) :: mpi_atmtab(:)
1693 : integer,optional,intent(in) :: l_orbmom,atom_orbmom
1694 : real(dp),intent(in) :: dmatpawu(dimdmat,dimdmat,nspinor*nsppol,natpawu*impose_dmat)
1695 : real(dp),intent(in) :: spinat(3,natom)
1696 : type(paw_ij_type),intent(inout) :: paw_ij(my_natom)
1697 : type(pawrhoij_type),intent(in) :: pawrhoij(my_natom)
1698 : type(pawtab_type),intent(in) :: pawtab(ntypat)
1699 :
1700 : !Local variables ---------------------------------------
1701 : !scalars
1702 : integer,parameter :: limp=0 ! could become an input variable
1703 : integer :: at_indx,cplex_dij,cplex_rhoij,dmatudiag_loc,iafm,iatom,iatom_tot,iatpawu,icount
1704 : integer :: ilm,im1,im2,in1,in2,info,iplex,irot,ispden, irhoij,itypat,jlm,jrhoij
1705 : integer :: jspden,klmn,kspden,lcur,ldim,lmax,lmin,lpawu,lwork,my_comm_atom,ndij,nmat,nspden,nsploop
1706 : logical,parameter :: afm_noncoll=.true. ! TRUE if antiferro symmetries are used with non-collinear magnetism
1707 : logical :: antiferro,my_atmtab_allocated,noccsym_error,paral_atom,use_afm
1708 : ! real(dp),parameter :: invsqrt2=one/sqrt2
1709 : real(dp) :: factafm,mnorm,mx,my,mz,ntot,nup,ndn,snorm,sx,sy,szm,szp
1710 : character(len=4) :: wrt_mode
1711 : character(len=500) :: message
1712 : !arrays
1713 : integer :: nsym_used(2)
1714 2729 : integer,pointer :: my_atmtab(:)
1715 : real(dp) :: ro(2),sumocc(2)
1716 2729 : real(dp),allocatable :: eig(:),hdp(:,:,:),hdp2(:,:),noccmmptemp(:,:,:,:),noccmmp_tmp(:,:,:,:)
1717 2729 : real(dp),allocatable :: rwork(:),noccmmp2(:,:,:,:),nocctot2(:)
1718 2729 : complex(dp),allocatable :: noccmmp_ylm(:,:,:),noccmmp_jmj(:,:),noccmmp_slm(:,:,:)
1719 2729 : complex(dp),allocatable :: zhdp(:,:),zhdp2(:,:),znoccmmp_tmp(:,:),zwork(:)
1720 : character(len=9),parameter :: dspin(6)= (/"up ","down ","up-up ","down-down","Re[up-dn]","Im[up-dn]"/)
1721 : character(len=9),parameter :: dspinc(6)= (/"up ","down ","up-up ","down-down","up-dn ","dn-up "/)
1722 : ! character(len=9),parameter :: dspinc2(6)=(/"up ","down ","dn-dn ","up-up ","dn-up ","up-dn "/)
1723 : character(len=9),parameter :: dspinm(6)= (/"dn ","up i ","n ","mx ","my ","mz "/)
1724 2729 : type(coeff4_type),allocatable :: tmp_noccmmp(:)
1725 :
1726 2729 : real(dp),allocatable :: l_noccmmp_tmp(:,:,:,:)
1727 : real(dp),optional,allocatable :: my_l_occmat(:,:,:,:)
1728 : logical :: cal_lmom
1729 : integer :: atom_min,atom_max
1730 : !*********************************************************************
1731 :
1732 : DBG_ENTER("COLL")
1733 : !in case of calculating orbital magnetic moments, only the occupation matrix for atoms atom_orbmom and orbital l_orbmom
1734 : !is calculated and returned in my_l_occmat.
1735 2729 : if (present(l_orbmom) .and. present(atom_orbmom)) then
1736 72 : cal_lmom= .true.
1737 72 : atom_min=atom_orbmom
1738 72 : atom_max=atom_orbmom
1739 : else
1740 : cal_lmom=.false.
1741 : atom_min=1
1742 : atom_max=my_natom
1743 : end if
1744 :
1745 : !Tests
1746 2729 : if (my_natom>0) then
1747 2729 : if (nsppol/=paw_ij(1)%nsppol) then
1748 0 : message='inconsistent values for nsppol!'
1749 0 : ABI_BUG(message)
1750 : end if
1751 2729 : if (compute_dmat>0) then
1752 : if (pawrhoij(1)%nspden/=paw_ij(1)%nspden.and.&
1753 2623 : & pawrhoij(1)%nspden/=4.and.paw_ij(1)%nspden/=1) then
1754 0 : message=' inconsistent values for nspden!'
1755 0 : ABI_BUG(message)
1756 : end if
1757 : end if
1758 2729 : if (pawrhoij(1)%qphase==2) then
1759 0 : message='setnoccmmp not compatible with qphase=2!'
1760 0 : ABI_BUG(message)
1761 : end if
1762 : end if
1763 2729 : if (usepawu/=0.and.useexexch/=0) then
1764 0 : message='usepawu/=0 and useexexch>0 not allowed!'
1765 0 : ABI_BUG(message)
1766 : end if
1767 2729 : if (impose_dmat/=0.and.dimdmat==0) then
1768 0 : message='dmatpawu must be allocated when impose_dmat/=0!'
1769 0 : ABI_BUG(message)
1770 : end if
1771 2729 : if (usepawu>0.and.compute_dmat/=0.and.impose_dmat/=0.and.pawang%nsym==0) then
1772 0 : message='pawang%zarot must be allocated!'
1773 0 : ABI_BUG(message)
1774 : end if
1775 :
1776 : !Some inits
1777 2729 : if (usepawu==0.and.useexexch==0) return
1778 2729 : nspden=1;ndij=1;cplex_dij=1
1779 2729 : if (my_natom>0) then
1780 2729 : nspden=paw_ij(1)%nspden
1781 2729 : ndij=paw_ij(1)%ndij
1782 2729 : cplex_dij=paw_ij(1)%cplex_dij
1783 : end if
1784 2729 : antiferro=(nspden==2.and.nsppol==1)
1785 2729 : use_afm=((antiferro).or.((nspden==4).and.afm_noncoll))
1786 2729 : dmatudiag_loc=dmatudiag
1787 2729 : if (dmatudiag==2.and.(dimdmat==0.or.impose_dmat==0)) dmatudiag_loc=1
1788 :
1789 : !Set up parallelism over atoms
1790 2729 : paral_atom=(present(comm_atom).and.(my_natom/=natom))
1791 2729 : nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
1792 2729 : my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
1793 2729 : call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,my_natom_ref=my_natom) !vz_d
1794 2729 : wrt_mode='COLL';if (paral_atom) wrt_mode='PERS'
1795 :
1796 : !If needed, store dmatpu in suitable format in tmp_noccmmp
1797 2729 : if (usepawu/=0.and.impose_dmat/=0) then
1798 79 : iatpawu=0
1799 474 : ABI_MALLOC(tmp_noccmmp,(natom))
1800 316 : do iatom_tot=1,natom
1801 237 : itypat=typat(iatom_tot)
1802 237 : lpawu=pawtab(itypat)%lpawu
1803 316 : if (lpawu/=-1) then
1804 130 : iatpawu=iatpawu+1
1805 130 : if (ndij/=4) then
1806 636 : ABI_MALLOC(tmp_noccmmp(iatom_tot)%value,(cplex_dij,2*lpawu+1,2*lpawu+1,nsppol))
1807 : tmp_noccmmp(iatom_tot)%value(1,1:2*lpawu+1,1:2*lpawu+1,1:nsppol)=&
1808 5226 : & dmatpawu(1:2*lpawu+1,1:2*lpawu+1,1:nsppol,iatpawu)
1809 : else
1810 144 : ABI_MALLOC(tmp_noccmmp(iatom_tot)%value,(cplex_dij,2*lpawu+1,2*lpawu+1,ndij))
1811 7800 : tmp_noccmmp(iatom_tot)%value=zero
1812 : if (limp==0) then ! default reading
1813 24 : snorm=sqrt(spinat(1,iatom_tot)**2+spinat(1,iatom_tot)**2+spinat(3,iatom_tot)**2)
1814 24 : if (snorm>tol12.and.nspden/=1) then
1815 18 : sx=half*spinat(1,iatom_tot)/snorm
1816 18 : sy=half*spinat(2,iatom_tot)/snorm
1817 18 : szp=half*(one+spinat(3,iatom_tot)/snorm)
1818 18 : szm=half*(one-spinat(3,iatom_tot)/snorm)
1819 : else
1820 : sx=zero;sy=zero
1821 : szp=half;szm=half
1822 : end if
1823 144 : do im2=1,2*lpawu+1
1824 744 : do im1=1,2*lpawu+1
1825 600 : nup=dmatpawu(im1,im2,1,iatpawu);ndn=dmatpawu(im1,im2,2,iatpawu)
1826 : ! if (nspden==1) tmp_noccmmp(iatom_tot)%value(1,im1,im2,1:2)=half*(nup+ndn)
1827 600 : tmp_noccmmp(iatom_tot)%value(1,im1,im2,1)=nup*szp+ndn*szm
1828 600 : tmp_noccmmp(iatom_tot)%value(1,im1,im2,2)=nup*szm+ndn*szp
1829 600 : tmp_noccmmp(iatom_tot)%value(1,im1,im2,3)=(nup-ndn)*sx
1830 720 : tmp_noccmmp(iatom_tot)%value(1,im1,im2,4)=(ndn-nup)*sy
1831 : end do
1832 : end do
1833 :
1834 : else if (limp>=1) then
1835 : ABI_MALLOC(noccmmp_ylm,(2*lpawu+1,2*lpawu+1,ndij))
1836 : noccmmp_ylm=czero
1837 : ABI_MALLOC(noccmmp_slm,(2*lpawu+1,2*lpawu+1,ndij))
1838 : noccmmp_slm=czero
1839 : ABI_MALLOC(noccmmp_jmj,(2*(2*lpawu+1),2*(2*lpawu+1)))
1840 : noccmmp_jmj=czero
1841 : if(limp==1) then ! read input matrix in J,M_J basis (l-1/2, then l+1/2)
1842 : noccmmp_jmj=czero
1843 : do im1=1,2*lpawu+1
1844 : noccmmp_jmj(im1,im1)=cmplx(dmatpawu(im1,im1,1,iatpawu),zero,kind=dp)
1845 : noccmmp_jmj(im1+lpawu,im1+lpawu)=cmplx(dmatpawu(im1+lpawu,im1+lpawu,2,iatpawu),zero,kind=dp)
1846 : end do
1847 : write(message,'(2a,i3,a)') ch10,'== Atom ',iatom_tot,&
1848 : & ' == Imposed occupation matrix (in the J M_J basis: L-1/2 and L+1/2 states)'
1849 : call wrtout(std_out,message,wrt_mode)
1850 : call mat_mlms2jmj(lpawu,noccmmp_ylm,noccmmp_jmj,ndij,&
1851 : & 2,2,pawprtvol,std_out,wrt_mode) ! optspin=1: up spin are first
1852 : end if
1853 : if (limp==2) then ! read input matrix in Ylm basis
1854 : noccmmp_ylm=czero
1855 : do im1=1,2*lpawu+1
1856 : noccmmp_ylm(im1,im1,1)=cmplx(dmatpawu(im1,im1,1,iatpawu),zero,kind=dp)
1857 : noccmmp_ylm(im1,im1,2)=cmplx(dmatpawu(im1,im1,2,iatpawu),zero,kind=dp)
1858 : end do
1859 : write(message,'(2a,i3,a)') ch10,'== Atom ',iatom_tot,&
1860 : & ' == Imposed occupation matrix (in the Ylm basis), for dn and up spin'
1861 : call wrtout(std_out,message,wrt_mode)
1862 : end if
1863 : call mat_slm2ylm(lpawu,noccmmp_ylm,noccmmp_slm,ndij,&
1864 : & 2,2,pawprtvol,std_out,wrt_mode) ! optspin=1 because up spin are first
1865 : ! interchange upup and dndn
1866 : if (limp>=1) then
1867 : tmp_noccmmp(iatom_tot)%value(1,:,:,1)=real(noccmmp_slm(:,:,2))
1868 : tmp_noccmmp(iatom_tot)%value(2,:,:,1)=aimag(noccmmp_slm(:,:,2))
1869 : tmp_noccmmp(iatom_tot)%value(1,:,:,2)=real(noccmmp_slm(:,:,1))
1870 : tmp_noccmmp(iatom_tot)%value(2,:,:,2)=aimag(noccmmp_slm(:,:,1))
1871 : tmp_noccmmp(iatom_tot)%value(1,:,:,3)=real(noccmmp_slm(:,:,4))
1872 : tmp_noccmmp(iatom_tot)%value(2,:,:,3)=aimag(noccmmp_slm(:,:,4))
1873 : tmp_noccmmp(iatom_tot)%value(1,:,:,4)=real(noccmmp_slm(:,:,3))
1874 : tmp_noccmmp(iatom_tot)%value(2,:,:,4)=aimag(noccmmp_slm(:,:,3))
1875 : end if
1876 : if(abs(pawprtvol)>2) then
1877 : write(message, '(2a)' ) ch10,&
1878 : & " Check Imposed density matrix in different basis"
1879 : call wrtout(std_out,message,wrt_mode)
1880 : call mat_slm2ylm(lpawu,noccmmp_slm,noccmmp_ylm,ndij,&
1881 : & 1,2,pawprtvol,std_out,wrt_mode) ! optspin=1 because up spin are first
1882 : call mat_mlms2jmj(lpawu,noccmmp_ylm,noccmmp_jmj,ndij,1,2,&
1883 : & pawprtvol,std_out,wrt_mode) ! optspin=1: up spin are first
1884 : end if
1885 : ABI_FREE(noccmmp_ylm)
1886 : ABI_FREE(noccmmp_jmj)
1887 : ABI_FREE(noccmmp_slm)
1888 : end if
1889 : end if
1890 : end if
1891 : end do
1892 : end if ! impose_dmat/=0
1893 :
1894 : !Print message
1895 2729 : if (usepawu/=0.and.impose_dmat/=0) then
1896 79 : if (dmatudiag_loc/=2) then
1897 79 : write(message,'(6a)') ch10,'Occupation matrix for correlated orbitals is kept constant',ch10,&
1898 79 : & 'and equal to dmatpawu from input file !',ch10,&
1899 158 : & '----------------------------------------------------------'
1900 : else
1901 0 : write(message,'(6a)') ch10,'Occupation matrix for correlated orbitals is imposed',ch10,&
1902 0 : & 'and equal to dmatpawu in the diagonal basis !',ch10,&
1903 0 : & '----------------------------------------------------------'
1904 : end if
1905 79 : call wrtout(std_out,message,'COLL')
1906 : end if
1907 :
1908 2729 : if (usepawu/=0.and.dmatudiag_loc/=0.and.my_natom>0) then
1909 32 : write(message,'(4a)') ch10,'Diagonalized occupation matrix "noccmmp" is printed !',ch10,&
1910 64 : & '-------------------------------------------------------------'
1911 32 : call wrtout(std_out,message,wrt_mode)
1912 : end if
1913 :
1914 : !Loops over atoms
1915 8049 : do iatom=atom_min,atom_max
1916 5320 : iatom_tot=iatom;if (paral_atom) iatom_tot=my_atmtab(iatom)
1917 5320 : itypat=pawrhoij(iatom)%itypat
1918 5320 : cplex_rhoij=pawrhoij(iatom)%cplex_rhoij
1919 :
1920 5320 : if (.not. cal_lmom) then
1921 5248 : if (useexexch/=0) then
1922 4 : lcur=pawtab(itypat)%lexexch
1923 5244 : else if (usepawu/=0) then
1924 5244 : lcur=pawtab(itypat)%lpawu
1925 : end if
1926 : end if
1927 :
1928 72 : if (cal_lmom) lcur=l_orbmom
1929 8049 : if (lcur/=-1) then
1930 :
1931 : ! ########################################################################################
1932 : ! # Compute nocc_mmp
1933 : ! ########################################################################################
1934 3613 : if ((usepawu/=0.and.compute_dmat/=0).or.useexexch/=0) then
1935 :
1936 :
1937 20736 : ABI_MALLOC(l_noccmmp_tmp,(cplex_dij,2*lcur+1,2*lcur+1,ndij))
1938 553306 : l_noccmmp_tmp(:,:,:,:)=zero
1939 :
1940 : ! Loop over spin components
1941 17280 : ABI_MALLOC(noccmmptemp,(cplex_dij,2*lcur+1,2*lcur+1,ndij))
1942 553306 : noccmmptemp(:,:,:,:)=zero
1943 3456 : if(ndij==4) then
1944 4365 : ABI_MALLOC(noccmmp2,(cplex_dij,2*lcur+1,2*lcur+1,ndij))
1945 : end if
1946 3456 : if(ndij==4) then
1947 873 : if(allocated(nocctot2)) then
1948 0 : ABI_FREE(nocctot2)
1949 : end if
1950 2619 : ABI_MALLOC(nocctot2,(ndij))
1951 : end if
1952 3456 : nsploop=ndij
1953 11602 : do ispden=1,nsploop
1954 8146 : jrhoij=1
1955 758981 : do irhoij=1,pawrhoij(iatom)%nrhoijsel
1956 750835 : klmn=pawrhoij(iatom)%rhoijselect(irhoij)
1957 750835 : im1=pawtab(itypat)%klmntomn(1,klmn)
1958 750835 : im2=pawtab(itypat)%klmntomn(2,klmn)
1959 750835 : in1=pawtab(itypat)%klmntomn(3,klmn)
1960 750835 : in2=pawtab(itypat)%klmntomn(4,klmn)
1961 750835 : lmin=pawtab(itypat)%indklmn(3,klmn)
1962 750835 : lmax=pawtab(itypat)%indklmn(4,klmn)
1963 :
1964 750835 : ro(1:2)=zero
1965 2028274 : ro(1:cplex_rhoij)=pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+cplex_rhoij-1,ispden)
1966 798429 : if (ndij==1) ro(1:2)=half*ro(1:2)
1967 : ! Non-collinear magnetism: keep n, m storage because
1968 : ! it is easier for the computation of noccmmp from rhoij)
1969 :
1970 750835 : if(lmin==0.and.lmax==2*lcur) then
1971 257514 : icount=in1+(in2*(in2-1))/2
1972 257514 : if(pawtab(itypat)%ij_proj<icount .and. (.not. cal_lmom) ) then
1973 0 : message='PAW+U: Problem in the loop calculating noccmmp!'
1974 0 : ABI_BUG(message)
1975 : end if
1976 257514 : if(in1/=in2) then
1977 100206 : if(im2<=im1) then
1978 : noccmmptemp(1:cplex_dij,im1,im2,ispden)=noccmmptemp(1:cplex_dij,im1,im2,ispden) &
1979 168844 : & +ro(1:cplex_dij)*pawtab(itypat)%phiphjint(icount)
1980 : end if
1981 : end if
1982 257514 : if(im2>=im1) then
1983 : l_noccmmp_tmp(1:cplex_dij,im1,im2,ispden)=l_noccmmp_tmp(1:cplex_dij,im1,im2,ispden) &
1984 574284 : & +ro(1:cplex_dij)*pawtab(itypat)%phiphjint(icount)
1985 : end if
1986 : end if
1987 758981 : jrhoij=jrhoij+cplex_rhoij
1988 : end do ! irhoij
1989 48706 : do im2=1,2*lcur+1
1990 172255 : do im1=1,im2
1991 : l_noccmmp_tmp(1,im1,im2,ispden)=l_noccmmp_tmp(1,im1,im2,ispden) &
1992 123549 : & +noccmmptemp(1,im2,im1,ispden)
1993 123549 : if(cplex_dij==2) l_noccmmp_tmp(2,im1,im2,ispden)=l_noccmmp_tmp(2,im1,im2,ispden) &
1994 93308 : & -noccmmptemp(2,im2,im1,ispden)
1995 : end do
1996 : end do
1997 52162 : do im1=1,2*lcur+1
1998 172255 : do im2=1,im1
1999 123549 : l_noccmmp_tmp(1,im1,im2,ispden)=l_noccmmp_tmp(1,im2,im1,ispden)
2000 164109 : if(cplex_dij==2) l_noccmmp_tmp(2,im1,im2,ispden)=-l_noccmmp_tmp(2,im2,im1,ispden)
2001 : end do
2002 : end do
2003 : end do ! ispden
2004 3456 : ABI_FREE(noccmmptemp)
2005 : ! Compute noccmmp2, occupation matrix in the spin basis (upup, dndn, updn, dnup)
2006 3456 : if(ndij==4) then
2007 285997 : noccmmp2(:,:,:,:)=zero
2008 5230 : do im1=1,2*lcur+1
2009 27247 : do im2=1,2*lcur+1
2010 22017 : noccmmp2(1,im1,im2,1)=half*(l_noccmmp_tmp(1,im1,im2,1)+l_noccmmp_tmp(1,im1,im2,4))
2011 22017 : noccmmp2(2,im1,im2,1)=half*(l_noccmmp_tmp(2,im1,im2,1)+l_noccmmp_tmp(2,im1,im2,4))
2012 22017 : noccmmp2(1,im1,im2,2)=half*(l_noccmmp_tmp(1,im1,im2,1)-l_noccmmp_tmp(1,im1,im2,4))
2013 22017 : noccmmp2(2,im1,im2,2)=half*(l_noccmmp_tmp(2,im1,im2,1)-l_noccmmp_tmp(2,im1,im2,4))
2014 22017 : noccmmp2(1,im1,im2,3)=half*(l_noccmmp_tmp(1,im1,im2,2)+l_noccmmp_tmp(2,im1,im2,3))
2015 22017 : noccmmp2(2,im1,im2,3)=half*(l_noccmmp_tmp(2,im1,im2,2)-l_noccmmp_tmp(1,im1,im2,3))
2016 22017 : noccmmp2(1,im1,im2,4)=half*(l_noccmmp_tmp(1,im1,im2,2)-l_noccmmp_tmp(2,im1,im2,3))
2017 26374 : noccmmp2(2,im1,im2,4)=half*(l_noccmmp_tmp(2,im1,im2,2)+l_noccmmp_tmp(1,im1,im2,3))
2018 : end do
2019 : end do
2020 873 : if(abs(pawprtvol)>=1 .and. (.not. cal_lmom)) then
2021 783 : write(message,'(2a)') ch10,"== Calculated occupation matrix for correlated orbitals in the n, m basis :"
2022 783 : call wrtout(std_out,message,wrt_mode)
2023 3915 : do ispden=1,ndij
2024 3132 : write(message,'(3a)') ch10,"Calculated occupation matrix for component ",trim(dspinm(ispden+2*(ndij/4)))
2025 3132 : call wrtout(std_out,message,wrt_mode)
2026 19831 : 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 )
2027 15916 : if(cplex_dij==1)&
2028 : & write(message,'(12(1x,9(1x,f10.5)))')&
2029 0 : & (l_noccmmp_tmp(1,im2,im1,ispden),im2=1,lcur*2+1)
2030 15916 : if(cplex_dij==2)&
2031 : ! & write(message,'(12(1x,9(1x,"(",f7.3,",",f7.3,")")))')&
2032 : & write(message,'(12(1x,9(1x,"(",f10.5,",",f10.5,")")))')&
2033 97288 : & (l_noccmmp_tmp(:,im2,im1,ispden),im2=1,lcur*2+1)
2034 19048 : call wrtout(std_out,message,wrt_mode)
2035 : end do
2036 : end do
2037 : end if ! pawprtvol >=1
2038 : end if
2039 :
2040 : ! Compute total number of electrons per spin
2041 3456 : if (.not. cal_lmom) then
2042 11242 : paw_ij(iatom)%nocctot(:)=zero ! contains nmmp in the n m representation
2043 6588 : if(ndij==4) nocctot2(:)=zero ! contains nmmp in the upup dndn updn dnup representation
2044 11242 : do ispden=1,ndij
2045 50650 : do im1=1,2*lcur+1
2046 47266 : if(ndij==4) then
2047 16276 : paw_ij(iatom)%nocctot(ispden)=paw_ij(iatom)%nocctot(ispden)+l_noccmmp_tmp(1,im1,im1,ispden)
2048 16276 : nocctot2(ispden)=nocctot2(ispden)+noccmmp2(1,im1,im1,ispden)
2049 : else
2050 23132 : paw_ij(iatom)%nocctot(ispden)=paw_ij(iatom)%nocctot(ispden)+l_noccmmp_tmp(1,im1,im1,ispden)
2051 : end if
2052 : end do
2053 : end do
2054 : end if
2055 : ! noccmmp will now be in the up up , dn dn... representation and now n_mmp=<m|n|mp> instead of <mp|n|m> !
2056 3456 : if(ndij==4) then
2057 4365 : do ispden=1,ndij
2058 11349 : do iplex=1,cplex_dij
2059 45332 : do im1=1,2*lcur+1
2060 217976 : do im2=1,2*lcur+1
2061 210992 : l_noccmmp_tmp(iplex,im1,im2,ispden)=noccmmp2(iplex,im2,im1,ispden) ! now, noccmmp is in the upup dndn updn dnup representation
2062 : end do
2063 : end do
2064 : end do
2065 : end do
2066 873 : ABI_FREE(noccmmp2)
2067 : end if
2068 : ! Printing of new nocc_mmp
2069 3456 : if (.not. cal_lmom) then
2070 3384 : if ((usepawu/=0.and.abs(usepawu)<10).or.(usepawu>=10.and.pawprtvol>=3)) then
2071 3205 : write(message, '(2a)' ) ch10, &
2072 6410 : & '========== DFT+U DATA =================================================== '
2073 : end if
2074 3384 : if (useexexch/=0) then
2075 2 : write(message, '(2a)' ) ch10, &
2076 4 : & '======= Local ex-exchange (PBE0) DATA =================================== '
2077 : end if
2078 3384 : if (((usepawu/=0.and.abs(usepawu)<10).or.(usepawu>=10.and.pawprtvol>=3)).or.useexexch/=0) then
2079 3207 : call wrtout(std_out,message,wrt_mode)
2080 : end if
2081 3384 : if (usepawu>=10.and.pawprtvol>=3) then
2082 169 : write(message, '(6a)' ) ch10,' ( A DFT+DMFT calculation is carried out ',&
2083 169 : ch10,' Thus, the following DFT+U occupation matrices are not physical ',&
2084 338 : ch10,' and just informative )'
2085 169 : call wrtout(std_out,message,wrt_mode)
2086 : end if
2087 3384 : if(abs(usepawu)<10.or.pawprtvol>=3) then ! Always write except if DMFT and pawprtvol low
2088 3207 : write(message,'(2a,i5,a,i4,a)') ch10,"====== For Atom", iatom_tot,&
2089 6414 : & ", occupations for correlated orbitals. l =",lcur,ch10
2090 3207 : call wrtout(std_out,message,wrt_mode)
2091 3207 : if(ndij==2) then
2092 5682 : do ispden=1,2
2093 3788 : write(message,'(a,i4,3a,f10.5)') "Atom", iatom_tot,". Occupations for spin ",&
2094 7576 : & trim(dspin(ispden))," =",paw_ij(iatom)%nocctot(ispden)
2095 5682 : call wrtout(std_out,message,wrt_mode)
2096 : end do
2097 1894 : write(message,'(a,i4,a,2x,e16.8)') "=> On atom",iatom_tot,", local Mag. is ",&
2098 3788 : & paw_ij(iatom)%nocctot(2)-paw_ij(iatom)%nocctot(1)
2099 1894 : call wrtout(std_out,message,wrt_mode)
2100 : end if
2101 3207 : if(ndij==4) then
2102 801 : ntot=paw_ij(iatom)%nocctot(1)
2103 801 : mx=paw_ij(iatom)%nocctot(2)
2104 801 : my=paw_ij(iatom)%nocctot(3)
2105 801 : mz=paw_ij(iatom)%nocctot(4)
2106 801 : mnorm=sqrt(mx*mx+my*my+mz*mz)
2107 801 : nup=nocctot2(1)
2108 801 : ndn=nocctot2(2)
2109 801 : write(message,'(a,i4,a,2x,e16.8)') "=> On atom",iatom_tot,", local Mag. x is ",mx
2110 801 : call wrtout(std_out,message,wrt_mode)
2111 801 : write(message,'(14x,a,2x,e16.8)') " local Mag. y is ",my
2112 801 : call wrtout(std_out,message,wrt_mode)
2113 801 : write(message,'(14x,a,2x,e16.8)') " local Mag. z is ",mz
2114 801 : call wrtout(std_out,message,wrt_mode)
2115 801 : write(message,'(14x,a,2x,e16.8)') " norm of Mag. is ",mnorm
2116 801 : call wrtout(std_out,message,wrt_mode)
2117 801 : write(message,'(14x,a,2x,f10.5)') " occ. of majority spin is ",half*(ntot+mnorm) ! to be checked versus direct calc from noccmmp
2118 801 : call wrtout(std_out,message,wrt_mode)
2119 801 : if(abs(pawprtvol)>=1) write(message,'(14x,a,2x,f10.5)') " occ. for spin up (along z) ",nup
2120 801 : if(abs(pawprtvol)>=1) then
2121 783 : call wrtout(std_out,message,wrt_mode)
2122 : end if
2123 801 : write(message,'(14x,a,2x,f10.5)') " occ. of minority spin is ",half*(ntot-mnorm)
2124 801 : call wrtout(std_out,message,wrt_mode)
2125 801 : if(abs(pawprtvol)>=1) write(message,'(14x,a,2x,f10.5)') " occ. for spin dn (along z) ",ndn
2126 801 : if(abs(pawprtvol)>=1) then
2127 783 : call wrtout(std_out,message,wrt_mode)
2128 : end if
2129 : if(ndij==4) then
2130 801 : ABI_FREE(nocctot2)
2131 : end if
2132 : end if
2133 3207 : write(message,'(2a)') ch10,"== Calculated occupation matrix for correlated orbitals:"
2134 3207 : call wrtout(std_out,message,wrt_mode)
2135 10711 : do ispden=1,ndij
2136 7504 : write(message,'(3a)') ch10,"Calculated occupation matrix for component ",trim(dspinc(ispden+2*(ndij/4)))
2137 7504 : call wrtout(std_out,message,wrt_mode)
2138 48506 : do im1=1,lcur*2+1
2139 37618 : if(cplex_dij==1)&
2140 : & write(message,'(12(1x,9(1x,f10.5)))')&
2141 21342 : & (l_noccmmp_tmp(1,im1,im2,ispden),im2=1,lcur*2+1)
2142 37618 : if(cplex_dij==2)&
2143 : & write(message,'(12(1x,9(1x,"(",f7.3,",",f7.3,")")))')&
2144 99448 : & (l_noccmmp_tmp(:,im1,im2,ispden),im2=1,lcur*2+1)
2145 45122 : call wrtout(std_out,message,wrt_mode)
2146 : end do
2147 : end do
2148 : end if
2149 :
2150 : ! Transformation matrices: real->complex spherical harmonics (for test)
2151 3384 : if(ndij==4.and.abs(pawprtvol)>=0) then
2152 4005 : ABI_MALLOC(noccmmp_ylm,(2*lcur+1,2*lcur+1,ndij))
2153 103453 : noccmmp_ylm=czero
2154 3204 : ABI_MALLOC(noccmmp_slm,(2*lcur+1,2*lcur+1,ndij))
2155 103453 : noccmmp_slm=czero
2156 3204 : ABI_MALLOC(noccmmp_jmj,(2*(2*lcur+1),2*(2*lcur+1)))
2157 92111 : noccmmp_jmj=czero
2158 : ! go from real notation for complex noccmmp to complex notation in noccmmp_slm
2159 : noccmmp_slm(:,:,:)=cmplx(l_noccmmp_tmp(1,:,:,:)&
2160 103453 : & ,l_noccmmp_tmp(2,:,:,:),kind=dp)
2161 801 : call mat_slm2ylm(lcur,noccmmp_slm,noccmmp_ylm,ndij,1,1,pawprtvol,std_out,wrt_mode) ! optspin=1: up spin are first
2162 :
2163 4005 : do ispden=1,ndij
2164 3204 : write(message,'(3a)') ch10,"Calculated Ylm occupation matrix for component ",trim(dspinc(ispden+2*(ndij/4)))
2165 3204 : call wrtout(std_out,message,wrt_mode)
2166 20281 : do im1=1,lcur*2+1
2167 16276 : write(message,'(12(1x,9(1x,"(",f9.5,",",f9.5,")")))') (noccmmp_ylm(im1,im2,ispden),im2=1,lcur*2+1)
2168 19480 : call wrtout(std_out,message,wrt_mode)
2169 : end do
2170 : end do
2171 801 : call mat_mlms2jmj(lcur,noccmmp_ylm,noccmmp_jmj,ndij,1,1,pawprtvol,std_out,wrt_mode) ! optspin=1: up spin are first
2172 801 : ABI_FREE(noccmmp_ylm)
2173 801 : ABI_FREE(noccmmp_jmj)
2174 801 : ABI_FREE(noccmmp_slm)
2175 : end if !ndij==4
2176 537106 : paw_ij(iatom)%noccmmp(:,:,:,:)=zero
2177 537106 : paw_ij(iatom)%noccmmp(:,:,:,:)=l_noccmmp_tmp(:,:,:,:)
2178 3384 : ABI_FREE(l_noccmmp_tmp)
2179 : else
2180 :
2181 72 : if(allocated(my_l_occmat)) then
2182 72 : ABI_FREE(my_l_occmat)
2183 : end if
2184 :
2185 72 : if(allocated(nocctot2)) then
2186 72 : ABI_FREE(nocctot2)
2187 : end if
2188 360 : ABI_MALLOC(my_l_occmat,(cplex_dij,2*lcur+1,2*lcur+1,ndij))
2189 16200 : my_l_occmat(:,:,:,:)=zero
2190 16272 : my_l_occmat=l_noccmmp_tmp(:,:,:,:)
2191 72 : ABI_FREE(l_noccmmp_tmp)
2192 :
2193 : end if ! not cal_lmom
2194 :
2195 : end if ! impose_dmat==0
2196 :
2197 : ! ########################################################################################
2198 : ! # Diagonalize nocc_mmp
2199 : ! ########################################################################################
2200 3613 : if(usepawu/=0.and.dmatudiag_loc>0.and.(.not. cal_lmom)) then
2201 :
2202 32 : lpawu=lcur;ldim=2*lpawu+1
2203 160 : ABI_MALLOC(noccmmp_tmp,(1,ldim,ldim,ndij))
2204 32 : if (ndij==4) then
2205 0 : ABI_MALLOC(znoccmmp_tmp,(2*ldim,2*ldim))
2206 : end if
2207 :
2208 : ! Select noccmmp for this atom
2209 96 : do ispden=1,ndij
2210 3680 : noccmmp_tmp(1,:,:,ispden)=paw_ij(iatom)%noccmmp(1,:,:,ispden)
2211 : end do
2212 32 : if (ndij==4) then
2213 0 : do im2=1,ldim
2214 0 : do im1=1,ldim
2215 : znoccmmp_tmp(im1 , im2)=cmplx(paw_ij(iatom)%noccmmp(1,im1,im2,1)&
2216 0 : & ,paw_ij(iatom)%noccmmp(2,im1,im2,1),kind=dp)
2217 : znoccmmp_tmp(ldim+im1,ldim+im2)=cmplx(paw_ij(iatom)%noccmmp(1,im1,im2,2)&
2218 0 : & ,paw_ij(iatom)%noccmmp(2,im1,im2,2),kind=dp)
2219 : znoccmmp_tmp( im1,ldim+im2)=cmplx(paw_ij(iatom)%noccmmp(1,im1,im2,3)&
2220 0 : & ,paw_ij(iatom)%noccmmp(2,im1,im2,3),kind=dp)
2221 : znoccmmp_tmp(ldim+im1, im2)=cmplx(paw_ij(iatom)%noccmmp(1,im1,im2,4)&
2222 0 : & ,paw_ij(iatom)%noccmmp(2,im1,im2,4),kind=dp)
2223 : end do
2224 : end do
2225 : end if
2226 :
2227 : ! Diagonalize nocc_mmp
2228 32 : if (ndij/=4) then
2229 128 : ABI_MALLOC(hdp,(ldim,ldim,ndij))
2230 3680 : hdp=zero
2231 32 : lwork=3*ldim-1
2232 96 : ABI_MALLOC(rwork,(lwork))
2233 96 : ABI_MALLOC(eig,(ldim))
2234 96 : do ispden=1,ndij
2235 64 : call dsyev('v','u',ldim,noccmmp_tmp(1,:,:,ispden),ldim,eig,rwork,lwork,info)
2236 64 : if(info/=0) then
2237 0 : message=' Error in diagonalization of noccmmp (DSYEV)!'
2238 0 : ABI_ERROR(message)
2239 : end if
2240 544 : do ilm=1,ldim
2241 512 : hdp(ilm,ilm,ispden)=eig(ilm)
2242 : end do
2243 : end do ! ispden
2244 32 : ABI_FREE(rwork)
2245 32 : ABI_FREE(eig)
2246 : else
2247 0 : ABI_MALLOC(hdp,(2*ldim,2*ldim,1))
2248 0 : hdp=zero
2249 0 : lwork=4*ldim-1
2250 0 : ABI_MALLOC(rwork,(6*ldim-2))
2251 0 : ABI_MALLOC(zwork,(lwork))
2252 0 : ABI_MALLOC(eig,(2*ldim))
2253 0 : call zheev('v','u',2*ldim,znoccmmp_tmp,2*ldim,eig,zwork,lwork,rwork,info)
2254 0 : if(info/=0) then
2255 0 : message=' Error in diagonalization of znoccmmp_tmp (zheev) !'
2256 0 : ABI_ERROR(message)
2257 : end if
2258 0 : do ilm=1,2*ldim
2259 0 : hdp(ilm,ilm,1)=eig(ilm)
2260 : end do
2261 0 : ABI_FREE(rwork)
2262 0 : ABI_FREE(zwork)
2263 0 : ABI_FREE(eig)
2264 : end if
2265 :
2266 : ! Print diagonalized matrix and eigenvectors
2267 96 : do ispden=1,size(hdp,3)
2268 64 : write(message,'(2a,i3,a)') ch10,'== Atom ',iatom_tot,' == Diagonalized Occupation matrix'
2269 64 : if (ndij==1) write(message,fmt='(2a)') trim(message)," for spin up =="
2270 64 : if (ndij==2) write(message,fmt='(2a,i3,a)')trim(message)," for spin ",ispden," =="
2271 64 : if (ndij==4) write(message,fmt='(2a,i3,a)')trim(message)," =="
2272 64 : call wrtout(std_out,message,wrt_mode)
2273 544 : do ilm=1,size(hdp,1)
2274 448 : write(message,'(12(1x,9(1x,f10.5)))') (hdp(ilm,jlm,ispden),jlm=1,size(hdp,2))
2275 512 : call wrtout(std_out,message,wrt_mode)
2276 : end do
2277 : end do ! ispden
2278 32 : if(abs(pawprtvol)>=1) then
2279 32 : if (ndij/=4) then
2280 96 : do ispden=1,ndij
2281 64 : write(message,'(2a,i3,a)') ch10,'== Atom ',iatom_tot,' == Eigenvectors'
2282 64 : if (ndij==1) write(message,fmt='(2a)') trim(message),' for spin up =='
2283 64 : if (ndij==2) write(message,fmt='(2a,i3,a)')trim(message),' for spin ',ispden,' =='
2284 64 : call wrtout(std_out,message,wrt_mode)
2285 544 : do ilm=1,ldim
2286 448 : write(message,'(12(1x,9(1x,f10.5)))') (noccmmp_tmp(1,ilm,jlm,ispden),jlm=1,ldim)
2287 512 : call wrtout(std_out,message,wrt_mode)
2288 : end do
2289 : end do
2290 : else
2291 0 : write(message,'(2a,i3,a)') ch10,'== Atom ',iatom_tot,' == Eigenvectors (spinors) in the real harmonics basis =='
2292 0 : call wrtout(std_out,message,wrt_mode)
2293 0 : do ilm=1,2*ldim
2294 0 : write(message,'(12(1x,9(1x,"(",f7.3,",",f7.3,")")))') (znoccmmp_tmp(ilm,jlm),jlm=1,2*ldim)
2295 0 : call wrtout(std_out,message,wrt_mode)
2296 : end do
2297 : end if
2298 : end if
2299 :
2300 : ! Back rotation of diagonalized matrix and printing
2301 32 : if(abs(pawprtvol)>=1) then
2302 32 : if (ndij/=4) then
2303 128 : ABI_MALLOC(hdp2,(ldim,ldim))
2304 96 : do ispden=1,ndij
2305 64 : call dgemm('n','t',ldim,ldim,ldim,one,hdp(:,:,ispden),ldim,noccmmp_tmp(1,:,:,ispden),ldim,zero,hdp2,ldim)
2306 64 : call dgemm('n','n',ldim,ldim,ldim,one,noccmmp_tmp(1,:,:,ispden),ldim,hdp2,ldim,zero,hdp(:,:,ispden),ldim)
2307 3680 : noccmmp_tmp(1,:,:,ispden)=hdp(:,:,ispden)
2308 : end do ! ispden
2309 32 : ABI_FREE(hdp2)
2310 : else
2311 0 : ABI_MALLOC(zhdp,(2*ldim,2*ldim))
2312 0 : ABI_MALLOC(zhdp2,(2*ldim,2*ldim))
2313 0 : zhdp(:,:)=cmplx(hdp(:,:,1),zero,kind=dp)
2314 0 : zhdp2(:,:)=cmplx(zero,zero,kind=dp)
2315 0 : call zgemm('n','c',2*ldim,2*ldim,2*ldim,cone,zhdp,2*ldim,znoccmmp_tmp,2*ldim,czero,zhdp2,2*ldim)
2316 0 : zhdp(:,:)=cmplx(zero,zero,kind=dp)
2317 0 : call zgemm('n','n',2*ldim,2*ldim,2*ldim,cone,znoccmmp_tmp,2*ldim,zhdp2,2*ldim,czero,zhdp,2*ldim)
2318 0 : znoccmmp_tmp=zhdp
2319 0 : ABI_FREE(zhdp)
2320 0 : ABI_FREE(zhdp2)
2321 : end if
2322 32 : nmat=ndij ; if(ndij==4.and.cplex_dij==2) nmat=1
2323 96 : do ispden=1,nmat
2324 64 : write(message,'(2a,i3,a)') ch10,'== Atom ',iatom_tot,&
2325 128 : & ' == Rotated back diagonalized matrix'
2326 64 : if (ndij==1) write(message,fmt='(2a)') trim(message)," for spin up =="
2327 64 : if (ndij==2) write(message,fmt='(2a,i3,a)')trim(message)," for spin ",ispden," =="
2328 64 : if (ndij==4.and.cplex_dij==2) write(message,fmt='(4a)') trim(message)," for all component "
2329 64 : call wrtout(std_out,message,wrt_mode)
2330 544 : do ilm=1,ldim*cplex_dij
2331 448 : if(ndij==1.or.ndij==2)&
2332 : & write(message,'(12(1x,9(1x,f10.5)))')&
2333 448 : & (noccmmp_tmp(1,ilm,jlm,ispden),jlm=1,ldim)
2334 448 : if(ndij==4.and.cplex_dij==2)&
2335 : & write(message,'(12(1x,18(1x,"(",f7.3,",",f7.3,")")))')&
2336 0 : & (znoccmmp_tmp(ilm,jlm),jlm=1,ldim*cplex_dij)
2337 512 : call wrtout(std_out,message,wrt_mode)
2338 : end do
2339 : end do ! ispden
2340 : end if
2341 32 : ABI_FREE(hdp)
2342 :
2343 : end if ! dmatudiag_loc
2344 :
2345 : ! ########################################################################################
2346 : ! # Impose value of nocc_mmp from dmatpu; symetrize it
2347 : ! ########################################################################################
2348 3613 : if (usepawu/=0.and.impose_dmat/=0.and.(.not.cal_lmom)) then
2349 :
2350 130 : lpawu=lcur
2351 130 : nsploop=nsppol;if (ndij==4) nsploop=4
2352 130 : noccsym_error=.false.
2353 :
2354 : ! Loop over spin components
2355 366 : do ispden=1,nsploop
2356 236 : if (ndij/=4) then
2357 140 : jspden=min(3-ispden,paw_ij(iatom)%nsppol)
2358 96 : else if (ispden<=2) then
2359 48 : jspden=3-ispden
2360 : else
2361 : jspden=ispden
2362 : end if
2363 :
2364 : ! Loops over components of nocc_mmp
2365 1606 : do jlm=1,2*lpawu+1
2366 8096 : do ilm=1,2*lpawu+1
2367 :
2368 7860 : if(nsym>1.and.ndij<4) then
2369 :
2370 4220 : nsym_used(1:2)=0
2371 4220 : sumocc(1:2)=zero
2372 :
2373 : ! Accumulate values of nocc_mmp over symmetries
2374 156380 : do irot=1,nsym
2375 152160 : if ((symafm(irot)/=1).and.(.not.use_afm)) cycle
2376 148560 : kspden=ispden;if (symafm(irot)==-1) kspden=jspden
2377 148560 : factafm=one;if (ispden>3) factafm=dble(symafm(irot))
2378 148560 : iafm=1;if ((antiferro).and.(symafm(irot)==-1)) iafm=2
2379 148560 : nsym_used(iafm)=nsym_used(iafm)+1
2380 148560 : at_indx=indsym(4,irot,iatom_tot)
2381 1036700 : do im2=1,2*lpawu+1
2382 6443520 : do im1=1,2*lpawu+1
2383 : ! Be careful: use here R_rel^-1 in term of spherical harmonics
2384 : ! which is tR_rec in term of spherical harmonics
2385 : ! so, use transpose[zarot]
2386 : sumocc(iafm)=sumocc(iafm)+factafm*tmp_noccmmp(at_indx)%value(1,im1,im2,kspden) &
2387 : & *pawang%zarot(im1,ilm,lpawu+1,irot)&
2388 6291360 : & *pawang%zarot(im2,jlm,lpawu+1,irot)
2389 : ! sumocc(iafm)=sumocc(iafm)+factafm*tmp_noccmmp(at_indx)%value(im1,im2,kspden) &
2390 : ! & *pawang%zarot(ilm,im1,lpawu+1,irot)&
2391 : ! & *pawang%zarot(jlm,im2,lpawu+1,irot)
2392 : end do
2393 : end do
2394 : end do ! End loop over symmetries
2395 :
2396 : ! Store new values of nocc_mmp
2397 4220 : paw_ij(iatom)%noccmmp(1,ilm,jlm,ispden)=sumocc(1)/nsym_used(1)
2398 4220 : if (.not.noccsym_error)&
2399 : & noccsym_error=(abs(paw_ij(iatom)%noccmmp(1,ilm,jlm,ispden) &
2400 3926 : & -tmp_noccmmp(iatom_tot)%value(1,ilm,jlm,ispden))>tol5)
2401 :
2402 : ! Antiferromagnetic case: has to fill up "down" component of nocc_mmp
2403 4220 : if (antiferro.and.nsym_used(2)>0) paw_ij(iatom)%noccmmp(1,ilm,jlm,2)=sumocc(2)/nsym_used(2)
2404 :
2405 : else ! nsym=1
2406 :
2407 : ! Case without symetries
2408 7200 : paw_ij(iatom)%noccmmp(:,ilm,jlm,ispden)= tmp_noccmmp(iatom_tot)%value(:,ilm,jlm,ispden)
2409 : end if
2410 :
2411 : end do !ilm
2412 : end do !jlm
2413 : end do ! ispden
2414 366 : do ispden=1,nsploop
2415 236 : paw_ij(iatom)%nocctot(ispden)=zero ! contains nmmp in the n m representation
2416 1606 : do im1=1,2*lcur+1
2417 1476 : if(ndij==4.and.ispden==1) then
2418 : ! in this case, on computes total number or electron for double counting correction
2419 : paw_ij(iatom)%nocctot(ispden)=paw_ij(iatom)%nocctot(ispden)+&
2420 120 : & paw_ij(iatom)%noccmmp(1,im1,im1,1)+paw_ij(iatom)%noccmmp(1,im1,im1,2)
2421 1120 : else if(ndij==4.and.ispden==2) then
2422 : paw_ij(iatom)%nocctot(ispden)=paw_ij(iatom)%nocctot(ispden)+&
2423 120 : & paw_ij(iatom)%noccmmp(1,im1,im1,3)+paw_ij(iatom)%noccmmp(1,im1,im1,4)
2424 1000 : else if(ndij==4.and.ispden==3) then
2425 : paw_ij(iatom)%nocctot(ispden)=paw_ij(iatom)%nocctot(ispden)-&
2426 120 : & paw_ij(iatom)%noccmmp(2,im1,im1,3)+paw_ij(iatom)%noccmmp(2,im1,im1,4)
2427 880 : else if(ndij==4.and.ispden==4) then
2428 : paw_ij(iatom)%nocctot(ispden)=paw_ij(iatom)%nocctot(ispden)+&
2429 120 : & paw_ij(iatom)%noccmmp(2,im1,im1,1)-paw_ij(iatom)%noccmmp(2,im1,im1,2)
2430 : else
2431 : paw_ij(iatom)%nocctot(ispden)=paw_ij(iatom)%nocctot(ispden)+&
2432 760 : & paw_ij(iatom)%noccmmp(1,im1,im1,ispden)
2433 : end if
2434 : end do
2435 : end do ! ispden
2436 :
2437 : ! Printing of new nocc_mmp
2438 432 : do ispden=1,ndij
2439 302 : if(dmatudiag_loc==2) then
2440 0 : write(message,'(2a,i3,a)') ch10,'== Atom ',iatom_tot,&
2441 0 : & ' == Imposed occupation matrix (in the basis of diagonalization!!)'
2442 : else
2443 302 : write(message,'(2a,i3,a)') ch10,'== Atom ',iatom_tot,&
2444 604 : & ' == Imposed occupation matrix'
2445 : end if
2446 302 : if (ndij==1) write(message,fmt='(2a)') trim(message)," for spin up =="
2447 302 : if (ndij==2) write(message,fmt='(2a,i3,a)')trim(message)," for spin ",ispden," =="
2448 398 : if (ndij==4) write(message,fmt='(4a)') trim(message)," for component ", &
2449 192 : & trim(dspinc(ispden+2*(ndij/4)))," =="
2450 302 : call wrtout(std_out,message,wrt_mode)
2451 2002 : do ilm=1,2*lpawu+1
2452 1570 : if(cplex_dij==1)&
2453 : & write(message,'(12(1x,9(1x,f10.5)))')&
2454 6960 : & (paw_ij(iatom)%noccmmp(1,ilm,jlm,ispden),jlm=1,2*lpawu+1)
2455 1570 : if(cplex_dij==2)&
2456 : & write(message,'(12(1x,9(1x,"(",f7.3,",",f7.3,")")))')&
2457 7680 : & (paw_ij(iatom)%noccmmp(:,ilm,jlm,ispden),jlm=1,2*lpawu+1)
2458 1872 : call wrtout(std_out,message,wrt_mode)
2459 : end do
2460 : end do
2461 :
2462 : ! WARNING if symmetrization changes the matrix
2463 130 : if (noccsym_error) then
2464 : write(message, '(a,i4,6a)' ) &
2465 12 : ' After symmetrization, imposed occupation matrix for atom ',iatom_tot,ch10,&
2466 12 : & ' is different from dmatpawu value set in input file !',ch10,&
2467 12 : & ' It is likely that dmatpawu does not match the symmetry operations of the system.',ch10,&
2468 24 : & ' Action: change dmatpawu in input file or increase precision until 0.00001'
2469 12 : ABI_WARNING(message)
2470 : end if
2471 :
2472 : end if ! impose_dmat/=0
2473 :
2474 : ! ########################################################################################
2475 : ! # Rotate imposed occupation matrix in the non-diagonal basis
2476 : ! ########################################################################################
2477 3613 : if (usepawu/=0.and.impose_dmat/=0.and.dmatudiag_loc==2.and.(.not. cal_lmom)) then
2478 :
2479 0 : lpawu=lcur;ldim=2*lpawu+1
2480 :
2481 : ! Rotation of imposed nocc_mmp
2482 0 : if (ndij/=4) then
2483 0 : ABI_MALLOC(hdp2,(ldim,ldim))
2484 0 : do ispden=1,ndij
2485 : call dgemm('n','t',ldim,ldim,ldim,one,&
2486 0 : & paw_ij(iatom)%noccmmp(1,:,:,ispden),ldim,noccmmp_tmp(1,:,:,ispden),ldim,zero,hdp2,ldim)
2487 : call dgemm('n','n',ldim,ldim,ldim,one,&
2488 0 : & noccmmp_tmp(1,:,:,ispden),ldim,hdp2,ldim,zero,paw_ij(iatom)%noccmmp(1,:,:,ispden),ldim)
2489 : end do ! ispden
2490 0 : ABI_FREE(hdp2)
2491 : else
2492 0 : ABI_MALLOC(zhdp,(2*ldim,2*ldim))
2493 0 : ABI_MALLOC(zhdp2,(2*ldim,2*ldim))
2494 0 : do im2=1,ldim
2495 0 : do im1=1,ldim
2496 0 : zhdp( im1, im2)=cmplx(paw_ij(iatom)%noccmmp(1,im1,im2,1),zero,kind=dp) ! to be checked
2497 0 : zhdp(ldim+im1,ldim+im2)=cmplx(paw_ij(iatom)%noccmmp(1,im1,im2,2),zero,kind=dp) ! to be checked
2498 : zhdp( im1,ldim+im2)=&
2499 0 : & cmplx(paw_ij(iatom)%noccmmp(1,im1,im2,3),+paw_ij(iatom)%noccmmp(1,im1,im2,4),kind=dp) ! to be checked
2500 : zhdp(ldim+im1, im2)=&
2501 0 : & cmplx(paw_ij(iatom)%noccmmp(1,im2,im1,3),-paw_ij(iatom)%noccmmp(1,im2,im1,4),kind=dp) ! to be checked
2502 : end do
2503 : end do
2504 0 : call zgemm('n','c',2*ldim,2*ldim,2*ldim,cone,zhdp,2*ldim,znoccmmp_tmp,2*ldim,czero,zhdp2,2*ldim)
2505 0 : call zgemm('n','n',2*ldim,2*ldim,2*ldim,cone,znoccmmp_tmp,2*ldim,zhdp2,2*ldim,czero,zhdp,2*ldim)
2506 0 : do jlm=1,ldim
2507 0 : do ilm=1,ldim
2508 0 : paw_ij(iatom)%noccmmp(1,ilm,jlm,1)= real(znoccmmp_tmp( ilm, jlm)) ! to be checked
2509 0 : paw_ij(iatom)%noccmmp(1,ilm,jlm,2)= real(znoccmmp_tmp(ldim+ilm,ldim+jlm)) ! to be checked
2510 0 : paw_ij(iatom)%noccmmp(1,ilm,jlm,3)= real(znoccmmp_tmp( ilm,ldim+jlm)) ! to be checked
2511 0 : paw_ij(iatom)%noccmmp(1,ilm,jlm,4)=aimag(znoccmmp_tmp( ilm,ldim+jlm)) ! to be checked
2512 : end do
2513 : end do
2514 0 : ABI_FREE(zhdp)
2515 0 : ABI_FREE(zhdp2)
2516 : end if
2517 :
2518 : ! Printing of rotated imposed matrix
2519 0 : do ispden=1,ndij
2520 0 : write(message,'(2a,i3,a)') ch10,'== Atom ',iatom_tot,&
2521 0 : & ' == Imposed density matrix in original basis'
2522 0 : if (ndij==1) write(message,fmt='(2a)') trim(message)," for spin up =="
2523 0 : if (ndij==2) write(message,fmt='(2a,i3,a)')trim(message)," for spin ",ispden," =="
2524 0 : if (ndij==4) write(message,fmt='(4a)') trim(message)," for component ", &
2525 0 : & trim(dspin(ispden+2*(ndij/4)))," =="
2526 0 : call wrtout(std_out,message,wrt_mode)
2527 0 : do ilm=1,2*lpawu+1
2528 0 : write(message,'(12(1x,9(1x,f10.5)))') (paw_ij(iatom)%noccmmp(1,ilm,jlm,ispden),jlm=1,2*lpawu+1) ! to be checked
2529 0 : call wrtout(std_out,message,wrt_mode)
2530 : end do
2531 : end do ! ispden
2532 :
2533 : end if ! dmatudiag_loc==2
2534 :
2535 3613 : if (usepawu/=0.and.dmatudiag_loc>0) then
2536 32 : ABI_FREE(noccmmp_tmp)
2537 32 : if (ndij==4) then
2538 0 : ABI_FREE(znoccmmp_tmp)
2539 : end if
2540 : end if
2541 :
2542 3613 : paw_ij(iatom)%has_pawu_occ=2
2543 :
2544 : end if ! lcur
2545 : end do ! iatom
2546 :
2547 : !Memory deallocation
2548 2729 : if (usepawu/=0.and.impose_dmat/=0) then
2549 316 : do iatom_tot=1,natom
2550 237 : lpawu=pawtab(typat(iatom_tot))%lpawu
2551 316 : if (lpawu/=-1) then
2552 130 : ABI_FREE(tmp_noccmmp(iatom_tot)%value)
2553 : end if
2554 : end do
2555 316 : ABI_FREE(tmp_noccmmp)
2556 : end if
2557 :
2558 : !Destroy atom table used for parallelism
2559 2729 : call free_my_atmtab(my_atmtab,my_atmtab_allocated)
2560 :
2561 : DBG_EXIT("COLL")
2562 :
2563 5458 : end subroutine setnoccmmp
2564 : !!***
2565 :
2566 : !----------------------------------------------------------------------
2567 :
2568 : !----------------------------------------------------------------------
2569 :
2570 : !!****f* m_paw_correlations/setrhoijpbe0
2571 : !! NAME
2572 : !! setrhoijpbe0
2573 : !!
2574 : !! FUNCTION
2575 : !! PAW local exact exchange only:
2576 : !! Impose value of rhoij for f electrons using an auxiliairy file
2577 : !!
2578 : !! INPUTS
2579 : !! dtset <type(dataset_type)>=all input variables for this dataset
2580 : !! initialized= if 0, the initialization of the gstate run is not yet finished
2581 : !! istep=index of the number of steps in the routine scfcv
2582 : !! istep_mix=index of the number of steps for the SCF mixing (can be <istep)
2583 : !! mpi_atmtab(:)=--optional-- indexes of the atoms treated by current proc
2584 : !! comm_atom=--optional-- MPI communicator over atoms
2585 : !! mpi_comm_read=MPI communicator containing all the processes reading the PBE0 file
2586 : !! my_natom=number of atoms treated by current processor
2587 : !! natom=number of atoms in cell
2588 : !! ntypat=number of types of atoms in unit cell
2589 : !! pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
2590 : !! typat(natom)=type integer for each atom in cell
2591 : !!
2592 : !! SIDE EFFECTS
2593 : !! istep_mix=index of the number of steps for the SCF mixing (can be <istep)
2594 : !! pawrhoij(natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data
2595 : !!
2596 : !! NOTES
2597 : !! Only valid for f electrons !!!
2598 : !!
2599 : !! SOURCE
2600 :
2601 21 : subroutine setrhoijpbe0(dtset,initialized,istep,istep_mix,&
2602 21 : & mpi_comm_read,my_natom,natom,ntypat,pawrhoij,pawtab,typat,&
2603 21 : & mpi_atmtab,comm_atom) ! optional arguments (parallelism)
2604 :
2605 : !Arguments ---------------------------------------------
2606 : !scalars
2607 : integer,intent(in) :: initialized,istep,mpi_comm_read,my_natom,natom,ntypat
2608 : integer,intent(inout) :: istep_mix
2609 : integer,optional,intent(in) :: comm_atom
2610 : type(dataset_type),intent(in) :: dtset
2611 : !arrays
2612 : integer,intent(in) :: typat(natom)
2613 : integer,optional,target,intent(in) :: mpi_atmtab(:)
2614 : type(pawrhoij_type),intent(inout) :: pawrhoij(my_natom)
2615 : type(pawtab_type),intent(in) :: pawtab(ntypat)
2616 :
2617 : !Local variables ---------------------------------------
2618 : !scalars
2619 : integer,parameter :: ll=3
2620 : integer :: cplex_rhoij,iatom,iatom_tot,ierr,ii,ios,iread,irhoij,ispden,itypat,jj
2621 : integer :: klmn,my_comm_atom,my_rank,nselect,nstep1,nstep1_abs,rhoijshft,rhoijsz
2622 : logical :: my_atmtab_allocated,paral_atom,test0
2623 : character(len=9),parameter :: filnam='rhoijpbe0'
2624 : character(len=9),parameter :: dspin(6)=(/"up ","down ","up-up ","down-down","Re[up-dn]","Im[up-dn]"/)
2625 : character(len=500) :: strg, message
2626 : !arrays
2627 21 : integer, allocatable :: nspden_tmp(:)
2628 21 : integer,pointer :: my_atmtab(:)
2629 21 : real(dp),allocatable :: rhoijtmp(:,:),rhoijtmp1(:,:),rhoijtmp2(:,:,:,:)
2630 :
2631 : ! *********************************************************************
2632 :
2633 : DBG_ENTER("COLL")
2634 :
2635 : !Some limitation
2636 21 : if (my_natom>0) then
2637 21 : if (pawrhoij(1)%qphase==2) then
2638 0 : message='setrhoijpbe0 not compatible with qphase=2!'
2639 0 : ABI_BUG(message)
2640 : end if
2641 : end if
2642 :
2643 : !Test existence of file and open it
2644 21 : inquire(file=filnam,iostat=ios,exist=test0)
2645 21 : if(.not.test0) return
2646 :
2647 : !Look for parallelisation over atomic sites
2648 0 : paral_atom=(present(comm_atom).and.(my_natom/=natom))
2649 :
2650 : !Test if exact-exch. is on f electrons
2651 0 : test0=.false.
2652 0 : do itypat=1,ntypat
2653 0 : if (pawtab(itypat)%useexexch/=0.and.pawtab(itypat)%lexexch/=ll) test0=.true.
2654 : end do
2655 0 : if (test0) then
2656 : write(message, '(3a,i1,a)' ) &
2657 0 : & ' Local exact exchange: occ. matrix can only be imposed for l=',ll,' !'
2658 0 : ABI_ERROR(message)
2659 : end if
2660 :
2661 : !============================================================
2662 : !===== First case: no parallelisation over atomic sites =====
2663 : !============================================================
2664 :
2665 0 : if (.not.paral_atom) then
2666 :
2667 : ! Open file
2668 0 : if (open_file(filnam,message,unit=77,form='formatted') /= 0) then
2669 0 : ABI_ERROR(message)
2670 : end if
2671 :
2672 : ! Read step number and eventually exit
2673 0 : nstep1=0;test0=.false.
2674 0 : do while (.not.test0)
2675 0 : read(77,'(A)') strg
2676 0 : test0=(strg(1:1)/="#")
2677 0 : if (test0) read(unit=strg,fmt=*) nstep1
2678 : end do
2679 0 : nstep1_abs=abs(nstep1)
2680 0 : if (nstep1_abs==0.or.istep>nstep1_abs.or.(nstep1>0.and.initialized/=0)) then
2681 0 : close(77)
2682 : ! Reinitalize mixing when rhoij is allowed to change; for experimental purpose...
2683 0 : if (dtset%userib==1234.and.istep==1+nstep1_abs.and.(nstep1<0.or.initialized==0)) istep_mix=1
2684 0 : return
2685 : end if
2686 :
2687 : ! Loop on atoms
2688 0 : do iatom=1,natom
2689 0 : itypat=typat(iatom)
2690 0 : cplex_rhoij=pawrhoij(iatom)%cplex_rhoij
2691 :
2692 0 : if (pawtab(itypat)%useexexch/=0) then
2693 :
2694 : ! Set sizes depending on ll
2695 0 : rhoijsz=4*ll+2
2696 0 : rhoijshft=2*ll*ll
2697 :
2698 : ! Uncompress rhoij
2699 0 : ABI_MALLOC(rhoijtmp,(pawrhoij(iatom)%lmn2_size,pawrhoij(iatom)%nspden))
2700 0 : do ispden=1,pawrhoij(iatom)%nspden
2701 0 : rhoijtmp=zero
2702 0 : do irhoij=1,pawrhoij(iatom)%nrhoijsel
2703 0 : klmn=pawrhoij(iatom)%rhoijselect(irhoij)
2704 0 : rhoijtmp(klmn,ispden)=pawrhoij(iatom)%rhoijp(irhoij,ispden)
2705 : end do
2706 : end do
2707 : ! Read rhoij from file
2708 0 : ABI_MALLOC(rhoijtmp1,(rhoijsz,rhoijsz))
2709 0 : do ispden=1,pawrhoij(iatom)%nspden
2710 0 : do ii=1,rhoijsz
2711 0 : test0=.false.
2712 0 : do while (.not.test0)
2713 0 : read(77,'(A)') strg
2714 0 : test0=(strg(1:1)/="#")
2715 0 : if (test0) read(unit=strg,fmt=*) (rhoijtmp1(ii,jj), jj=1,rhoijsz)
2716 : end do
2717 : end do
2718 :
2719 : ! Impose rhoij
2720 0 : do jj=1,rhoijsz
2721 0 : do ii=1,jj
2722 0 : rhoijtmp((jj+rhoijshft)*((jj+rhoijshft)-1)/2+ii+rhoijshft,ispden)=rhoijtmp1(ii,jj)
2723 : end do
2724 : end do
2725 :
2726 : end do
2727 0 : ABI_FREE(rhoijtmp1)
2728 :
2729 : ! Compress rhoij
2730 0 : nselect=0 ; pawrhoij(iatom)%rhoijselect=0
2731 0 : pawrhoij(iatom)%rhoijp=zero
2732 0 : do klmn=1,pawrhoij(iatom)%lmn2_size
2733 0 : if (any(abs(rhoijtmp(klmn,:))>tol10)) then
2734 0 : nselect=nselect+1 ; ii=cplex_rhoij*(nselect-1)+1
2735 0 : do ispden=1,pawrhoij(iatom)%nspden
2736 0 : pawrhoij(iatom)%rhoijp(ii,ispden)=rhoijtmp(klmn,ispden)
2737 : end do
2738 0 : pawrhoij(iatom)%rhoijselect(nselect)=klmn
2739 : end if
2740 : end do
2741 0 : pawrhoij(iatom)%nrhoijsel=nselect
2742 0 : ABI_FREE(rhoijtmp)
2743 :
2744 : ! Print new rhoij
2745 0 : do ispden=1,pawrhoij(iatom)%nspden
2746 0 : write(message,'(2a,i3,a)') ch10,'== Atom ',iatom,&
2747 0 : & ' == Imposed occupation matrix'
2748 0 : if (pawrhoij(iatom)%nspden==1) write(message,fmt='(2a)') trim(message)," for spin up =="
2749 0 : if (pawrhoij(iatom)%nspden==2) write(message,fmt='(2a,i3,a)')trim(message)," for spin ",ispden," =="
2750 0 : if (pawrhoij(iatom)%nspden==4) write(message,fmt='(4a)') trim(message)," for component ", &
2751 0 : & trim(dspin(ispden+2*(pawrhoij(iatom)%nspden/4)))," =="
2752 0 : call wrtout(std_out,message,'COLL')
2753 : call pawio_print_ij(std_out,pawrhoij(iatom)%rhoijp(:,ispden),pawrhoij(iatom)%nrhoijsel,&
2754 : & pawrhoij(iatom)%cplex_rhoij,pawrhoij(iatom)%lmn_size,ll,&
2755 : & pawtab(itypat)%indlmn(1,1:pawtab(itypat)%lmn_size),&
2756 0 : & 1,-1,pawrhoij(iatom)%rhoijselect(:),-1.d0,1,mode_paral='COLL')
2757 : end do
2758 :
2759 : ! End loop on atoms
2760 : end if
2761 : end do
2762 :
2763 : ! Close file
2764 0 : close (77)
2765 :
2766 : else
2767 :
2768 : ! ============================================================
2769 : ! ====== 2nd case: no parallelisation over atomic sites =====
2770 : ! ============================================================
2771 :
2772 0 : my_rank=xmpi_comm_rank(mpi_comm_read)
2773 :
2774 : ! Read step number and eventually exit
2775 0 : iread=0
2776 0 : if (my_rank==0) then
2777 0 : if (open_file(filnam,message,unit=77,form='formatted') /=0 ) then
2778 0 : ABI_ERROR(message)
2779 : end if
2780 0 : nstep1=0;test0=.false.
2781 0 : do while (.not.test0)
2782 0 : read(77,'(A)') strg
2783 0 : test0=(strg(1:1)/="#")
2784 0 : if (test0) read(unit=strg,fmt=*) nstep1
2785 : end do
2786 0 : nstep1_abs=abs(nstep1)
2787 0 : if (nstep1_abs==0.or.istep>nstep1_abs.or.(nstep1>0.and.initialized/=0)) then
2788 0 : close(77)
2789 : ! Reinitalize mixing when rhoij is allowed to change; for experimental purpose...
2790 0 : if (dtset%userib==1234.and.istep==1+nstep1_abs.and.(nstep1<0.or.initialized==0)) istep_mix=1
2791 0 : iread=1
2792 : end if
2793 : end if
2794 0 : call xmpi_sum(iread,mpi_comm_read,ierr)
2795 0 : if (iread/=0) return
2796 :
2797 : ! Set up parallelism over atoms
2798 0 : nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
2799 0 : my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
2800 0 : call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,my_natom_ref=my_natom)
2801 :
2802 : ! Store number of component for rhoij
2803 0 : ABI_MALLOC(nspden_tmp,(natom))
2804 0 : nspden_tmp(:)=zero
2805 0 : do iatom=1,my_natom
2806 0 : iatom_tot=my_atmtab(iatom)
2807 0 : nspden_tmp(iatom_tot)=pawrhoij(iatom)%nspden
2808 : end do
2809 0 : call xmpi_sum(nspden_tmp,mpi_comm_read,ierr)
2810 :
2811 : ! To be improve if too much memory
2812 0 : ABI_MALLOC(rhoijtmp2,(natom,4,rhoijsz,rhoijsz))
2813 0 : rhoijtmp2=zero
2814 :
2815 : ! Read rhoij from file
2816 0 : if (my_rank==0) then
2817 0 : do iatom=1,natom
2818 0 : itypat=typat(iatom)
2819 0 : if (pawtab(itypat)%useexexch/=0) then
2820 0 : rhoijsz=4*ll+2
2821 0 : do ispden=1,nspden_tmp(iatom)
2822 0 : do ii=1,rhoijsz
2823 0 : test0=.false.
2824 0 : do while (.not.test0)
2825 0 : read(77,'(A)') strg
2826 0 : test0=(strg(1:1)/="#")
2827 0 : if (test0) read(unit=strg,fmt=*) (rhoijtmp2(iatom,ispden,ii,jj),jj=1,rhoijsz)
2828 : end do
2829 : end do
2830 : end do
2831 : end if
2832 : end do
2833 : end if
2834 0 : call xmpi_sum(rhoijtmp2,mpi_comm_read,ierr)
2835 :
2836 : ! Now, distribute rhoij
2837 0 : do iatom=1,my_natom
2838 0 : iatom_tot=my_atmtab(iatom)
2839 0 : itypat=pawrhoij(iatom)%itypat
2840 0 : cplex_rhoij=pawrhoij(iatom)%cplex_rhoij
2841 :
2842 0 : if (pawtab(itypat)%useexexch/=0) then
2843 :
2844 : ! Set sizes depending on ll
2845 0 : rhoijsz=4*ll+2
2846 0 : rhoijshft=2*ll*ll
2847 :
2848 : ! Uncompress rhoij
2849 0 : ABI_MALLOC(rhoijtmp,(pawrhoij(iatom)%lmn2_size,pawrhoij(iatom)%nspden))
2850 0 : do ispden=1,pawrhoij(iatom)%nspden
2851 0 : rhoijtmp=zero
2852 0 : do irhoij=1,pawrhoij(iatom)%nrhoijsel
2853 0 : klmn=pawrhoij(iatom)%rhoijselect(irhoij)
2854 0 : rhoijtmp(klmn,ispden)=pawrhoij(iatom)%rhoijp(irhoij,ispden)
2855 : end do
2856 :
2857 : ! Impose rhoij
2858 0 : do jj=1,rhoijsz
2859 0 : do ii=1,jj
2860 0 : rhoijtmp((jj+rhoijshft)*((jj+rhoijshft)-1)/2+ii+rhoijshft,ispden)=rhoijtmp2(iatom_tot,ispden,ii,jj)
2861 : end do
2862 : end do
2863 :
2864 : end do
2865 :
2866 : ! Compress rhoij
2867 0 : nselect=0 ; pawrhoij(iatom)%rhoijselect=0
2868 0 : pawrhoij(iatom)%rhoijp=zero
2869 0 : do klmn=1,pawrhoij(iatom)%lmn2_size
2870 0 : if (any(abs(rhoijtmp(klmn,:))>tol10)) then
2871 0 : nselect=nselect+1 ; ii=cplex_rhoij*(nselect-1)+1
2872 0 : do ispden=1,pawrhoij(iatom)%nspden
2873 0 : pawrhoij(iatom)%rhoijp(ii,ispden)=rhoijtmp(klmn,ispden)
2874 : end do
2875 0 : pawrhoij(iatom)%rhoijselect(nselect)=klmn
2876 : end if
2877 : end do
2878 0 : pawrhoij(iatom)%nrhoijsel=nselect
2879 0 : ABI_FREE(rhoijtmp)
2880 :
2881 : end if ! useexexch/=0
2882 :
2883 : ! Print new rhoij
2884 0 : do ispden=1,pawrhoij(iatom)%nspden
2885 0 : write(message,'(2a,i3,a)') ch10,'== Atom ',iatom,' == Imposed occupation matrix'
2886 0 : if (pawrhoij(iatom)%nspden==1) write(message,fmt='(2a)') trim(message)," for spin up =="
2887 0 : if (pawrhoij(iatom)%nspden==2) write(message,fmt='(2a,i3,a)')trim(message)," for spin ",ispden," =="
2888 0 : if (pawrhoij(iatom)%nspden==4) write(message,fmt='(4a)') trim(message)," for component ", &
2889 0 : & trim(dspin(ispden+2*(pawrhoij(iatom)%nspden/4)))," =="
2890 0 : call wrtout(std_out,message,'PERS')
2891 : call pawio_print_ij(std_out,pawrhoij(iatom)%rhoijp(:,ispden),pawrhoij(iatom)%nrhoijsel,&
2892 : & pawrhoij(iatom)%cplex_rhoij,pawrhoij(iatom)%lmn_size,ll,&
2893 : & pawtab(itypat)%indlmn(1,1:pawtab(itypat)%lmn_size),&
2894 0 : & 1,-1,pawrhoij(iatom)%rhoijselect(:),-1.d0,1,mode_paral='PERS')
2895 : end do
2896 :
2897 : ! end loop on atoms
2898 : end do
2899 :
2900 0 : ABI_FREE(nspden_tmp)
2901 0 : ABI_FREE(rhoijtmp2)
2902 :
2903 : ! Destroy atom table used for parallelism
2904 0 : call free_my_atmtab(my_atmtab,my_atmtab_allocated)
2905 :
2906 : ! ============================================================
2907 : end if ! paral_atom
2908 :
2909 : DBG_EXIT("COLL")
2910 :
2911 21 : end subroutine setrhoijpbe0
2912 : !!***
2913 :
2914 : !----------------------------------------------------------------------
2915 :
2916 : !!****f* m_paw_correlations/calc_ubare
2917 : !! NAME
2918 : !! calc_ubare
2919 : !!
2920 : !! FUNCTION
2921 : !! Calculate the bare interaction on atomic orbitals
2922 : !!
2923 : !! INPUTS
2924 : !! itypatcor = value of itypat for correlated species
2925 : !! lpawu = angular momentum for correlated species
2926 : !! pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data:
2927 : !! pawang
2928 : !! %lmax=Maximum value of angular momentum l+1
2929 : !! %gntselect((2*l_max-1)**2,l_max**2,l_max**2)=
2930 : !! selection rules for Gaunt coefficients
2931 : !! pawrad(ntypat) <type(pawrad_type)>=paw radial mesh and related data:
2932 : !! %mesh_size=Dimension of radial mesh
2933 : !! %rad(mesh_size)=The coordinates of all the points of the radial mesh
2934 : !! %radfact(mesh_size)=Factor used to compute radial integrals
2935 : !!
2936 : !! OUTPUT
2937 : !!
2938 : !! NOTES
2939 : !!
2940 : !! SOURCE
2941 :
2942 0 : subroutine calc_ubare(itypatcor,lpawu,pawang,pawrad,pawtab,rmax)
2943 :
2944 : !Arguments ------------------------------------
2945 : integer, intent(in) :: itypatcor,lpawu
2946 : type(pawang_type),intent(in) :: pawang
2947 : type(pawrad_type),intent(in) :: pawrad
2948 : type(pawtab_type),target,intent(in) :: pawtab
2949 : real(dp), optional, intent(in) :: rmax
2950 :
2951 : !Local variables ------------------------------
2952 : !scalars
2953 : integer :: ilmn,ilmn1,iln,iln1,isel,isel1,itypat,jlmn,jlmn1,jln,jln1
2954 : integer :: klm,klm1,klmn,klmn1,ll,lm0
2955 : integer :: lmin,lmax,lmn2_size,mesh_size,meshsz,mm
2956 : real(dp) :: norm,r_for_intg,rg,rg1,ubare,uint,uint_tmp
2957 : character(len=800) :: message
2958 : !arrays
2959 0 : real(dp),allocatable :: ff(:),gg(:),phiphj(:),phiphj1(:)
2960 :
2961 : !************************************************************************
2962 :
2963 0 : itypat=itypatcor
2964 : write(message,'(11a,f12.4,2a,i7,2a,f12.4,2a,i7,2a,f12.4)') &
2965 0 : & ch10," =======================================================================",ch10, &
2966 0 : & " == Calculation of diagonal bare Coulomb interaction on ATOMIC orbitals ",ch10, &
2967 0 : & " (it is assumed that the wavefunction for the first reference ",ch10, &
2968 0 : & " energy in PAW atomic data is an atomic eigenvalue)",ch10,ch10, &
2969 0 : & " Max value of the radius in atomic data file =", pawrad%rmax ,ch10, &
2970 0 : & " Max value of the mesh in atomic data file =", pawrad%mesh_size,ch10, &
2971 0 : & " PAW radius is =", pawtab%rpaw,ch10, &
2972 0 : & " PAW value of the mesh for integration is =", pawrad%int_meshsz,ch10, &
2973 0 : & " Integral of atomic wavefunction until rpaw =", pawtab%ph0phiint(1)
2974 0 : if(.not.present(rmax)) then
2975 0 : call wrtout(ab_out,message,'COLL')
2976 0 : call wrtout(std_out,message,'COLL')
2977 : end if
2978 :
2979 0 : mesh_size=pawrad%mesh_size
2980 :
2981 : ! Definition of the mesh used for integration.
2982 : if(present(rmax)) then
2983 0 : if(rmax>pawrad%rmax) then
2984 0 : write(message, '(a)' ) 'calc_ubare: the radius cannot be larger than the maximum radius of the mesh'
2985 0 : ABI_ERROR(message)
2986 : end if
2987 0 : meshsz=pawrad_ifromr(pawrad,rmax)+5
2988 0 : r_for_intg=rmax
2989 : else
2990 0 : meshsz=pawtab%partialwave_mesh_size
2991 0 : r_for_intg=pawrad%rad(meshsz) ! (we could use r_for_intg=-1)
2992 : end if
2993 :
2994 0 : lmn2_size=pawtab%lmn2_size
2995 0 : ABI_MALLOC(ff,(mesh_size))
2996 0 : ABI_MALLOC(gg,(mesh_size))
2997 0 : ABI_MALLOC(phiphj,(mesh_size))
2998 0 : ABI_MALLOC(phiphj1,(mesh_size))
2999 0 : do klmn=1,lmn2_size
3000 0 : ilmn=pawtab%indklmn(7,klmn);jlmn=pawtab%indklmn(8,klmn)
3001 : ! Select lpawu and first projectors il=jl=lpawu and first proj only
3002 : if (( pawtab%indklmn(3,klmn)+pawtab%indklmn(4,klmn)==2*lpawu).and. &
3003 : & (-pawtab%indklmn(3,klmn)+pawtab%indklmn(4,klmn)==2*lpawu).and. &
3004 0 : & (pawtab%indlmn(3,ilmn)==1).and.(pawtab%indlmn(3,jlmn)==1) ) then
3005 0 : klm=pawtab%indklmn(1,klmn);iln=pawtab%indlmn(5,ilmn);jln=pawtab%indlmn(5,jlmn)
3006 0 : lmin=pawtab%indklmn(3,klmn);lmax=pawtab%indklmn(4,klmn)
3007 0 : phiphj(1:meshsz)=pawtab%phi(1:meshsz,iln)*pawtab%phi(1:meshsz,jln)
3008 : !write(6,*) "A",klmn,pawtab%klmntomn(1,klmn),pawtab%klmntomn(2,klmn),&
3009 : !&pawtab%indklmn(7,klmn),pawtab%indklmn(8,klmn),pawtab%klmntomn(3,klmn),pawtab%klmntomn(4,klmn)
3010 0 : do ll=lmin,lmin,2
3011 0 : lm0=ll*ll+ll+1
3012 0 : ff(1:meshsz)=phiphj(1:meshsz)
3013 0 : call simp_gen(norm,ff,pawrad,r_for_intg=r_for_intg)
3014 0 : call poisson(ff,ll,pawrad,gg)
3015 0 : do klmn1=klmn,lmn2_size
3016 0 : ilmn1=pawtab%indklmn(7,klmn);jlmn1=pawtab%indklmn(8,klmn)
3017 : ! Select lpawu and first projectors il=jl=lpawu and first proj only
3018 : if (( pawtab%indklmn(3,klmn1)+pawtab%indklmn(4,klmn1)==2*lpawu).and. &
3019 : & (-pawtab%indklmn(3,klmn1)+pawtab%indklmn(4,klmn1)==2*lpawu).and. &
3020 0 : & (pawtab%indlmn(3,ilmn1)==1).and.(pawtab%indlmn(3,jlmn1)==1) ) then
3021 : !write(6,*) "A1",klmn1,pawtab%klmntomn(1,klmn1),pawtab%klmntomn(2,klmn1),&
3022 : !&pawtab%indklmn(7,klmn1),pawtab%indklmn(8,klmn1),pawtab%klmntomn(3,klmn1),pawtab%klmntomn(4,klmn1)
3023 0 : klm1=pawtab%indklmn(1,klmn1);iln1=pawtab%indlmn(5,ilmn1);jln1=pawtab%indlmn(5,jlmn1)
3024 0 : phiphj1(1:meshsz)=pawtab%phi(1:meshsz,iln1)*pawtab%phi(1:meshsz,jln1)
3025 : uint_tmp=zero
3026 0 : if ((ll==lmin)) then
3027 0 : ff(1)=zero
3028 0 : ff(2:meshsz)=phiphj1(2:meshsz)*gg(2:meshsz)*two/pawrad%rad(2:meshsz)
3029 0 : call simp_gen(uint_tmp,ff,pawrad,r_for_intg=r_for_intg)
3030 : end if
3031 0 : uint=zero
3032 0 : do mm=-ll,ll
3033 0 : isel =pawang%gntselect(lm0+mm,klm)
3034 0 : isel1=pawang%gntselect(lm0+mm,klm1)
3035 0 : if (isel>0.and.isel1>0) then
3036 0 : rg =pawang%realgnt(isel)
3037 0 : rg1=pawang%realgnt(isel1)
3038 0 : uint=uint+uint_tmp*rg*rg1*two_pi
3039 : end if
3040 : end do
3041 : if((pawtab%indklmn(5,klmn)==pawtab%indklmn(6,klmn)).and.&
3042 0 : & (pawtab%indklmn(5,klmn1)==pawtab%indklmn(6,klmn1)).and.&
3043 : & (pawtab%indklmn(5,klmn)==pawtab%indklmn(5,klmn1))) then
3044 0 : ubare=uint*Ha_eV
3045 : end if
3046 : end if
3047 : end do
3048 : end do
3049 : end if
3050 : end do
3051 0 : ABI_FREE(gg)
3052 0 : ABI_FREE(ff)
3053 0 : ABI_FREE(phiphj)
3054 0 : ABI_FREE(phiphj1)
3055 :
3056 0 : write(message,'(a,3(a,f12.4,a),2a,f12.4,a)') ch10," For an atomic wfn truncated at rmax =",r_for_intg,ch10,&
3057 0 : & " The norm of the wfn is =",norm,ch10,&
3058 0 : & " The bare interaction (no renormalization) =",ubare," eV",ch10,&
3059 0 : & " The bare interaction (for a renorm. wfn ) =",ubare/norm/norm," eV"
3060 0 : call wrtout(ab_out,message,'COLL')
3061 0 : call wrtout(std_out,message,'COLL')
3062 0 : if(r_for_intg < 10_dp .and. .not.present(rmax)) then
3063 0 : write(message,'(a,f6.2,4a)') ' ( WARNING: The radial mesh in the atomic data file is cut at',r_for_intg,ch10,&
3064 0 : & ' Use XML atomic data files to compute the bare Coulomb interaction',ch10,&
3065 0 : & ' on a true normalized atomic wavefunction )'
3066 0 : call wrtout(ab_out,message,'COLL')
3067 0 : call wrtout(std_out,message,'COLL')
3068 : end if
3069 0 : if(present(rmax)) then
3070 0 : write(message,'(2a)') " =======================================================================",ch10
3071 0 : call wrtout(ab_out,message,'COLL')
3072 0 : call wrtout(std_out,message,'COLL')
3073 : end if
3074 :
3075 0 : end subroutine calc_ubare
3076 : !!***
3077 : !!****f* m_paw_correlations/loc_orbmom_cal
3078 : !! NAME
3079 : !! loc_orbmom_cal
3080 : !!
3081 : !! FUNCTION
3082 : !! Calculate the orbital magnetic moments in PAW spheres
3083 : !!
3084 : !! INPUTS
3085 : !! INPUTS
3086 : !! compute_dmat= flag: if 1, nocc_{m,mp} is computed
3087 : !! dimdmat=first dimension of dmatpawu array
3088 : !! dmatpawu(dimdmat,dimdmat,nsppol*nspinor,natpawu)=input density matrix to be copied into noccmpp
3089 : !! dmatudiag= flag controlling the use of diagonalization:
3090 : !! 0: no diagonalization of nocc_{m,mp}
3091 : !! 1: diagonalized nocc_{m,mp} matrix is printed
3092 : !! 2: dmatpawu matrix is expressed in the basis where nocc_(m,mp} is diagonal
3093 : !! impose_dmat= flag: if 1, nocc_{m,mp} is replaced by dmatpawu
3094 : !! indsym(4,nsym,natom)=indirect indexing array for atom labels
3095 : !! mpi_atmtab(:)=--optional-- indexes of the atoms treated by current proc
3096 : !! comm_atom=--optional-- MPI communicator over atoms
3097 : !! my_natom=number of atoms treated by current processor
3098 : !! natom=number of atoms in cell
3099 : !! natpawu=number of atoms on which PAW+U is applied
3100 : !! nspinor=number of spinorial components of the wavefunctions
3101 : !! nsppol=number of independant spin components
3102 : !! nsym=number of symmetry elements in space group
3103 : !! ntypat=number of atom types
3104 : !! paw_ij(my_natom) <type(paw_ij_type)>=paw arrays given on (i,j) channels
3105 : !! pawang <type(pawang_type)>=paw angular mesh and related data
3106 : !! pawrhoij(my_natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data
3107 : !! pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data
3108 : !! spinat(3,matom)=initial spin of each atom, in unit of hbar/2
3109 : !! symafm(nsym)=(anti)ferromagnetic part of symmetry operations
3110 : !! typat(natom)=type for each atom
3111 : !! useexexch=1 if local-exact-exchange is activated
3112 : !! usepawu= /=0 if PAW+U is activated
3113 : !!
3114 : !! OUTPUT
3115 : !! printing the values of orbital magnetic moments for atoms in the output file
3116 : !!
3117 : !! SOURCE
3118 :
3119 18 : subroutine loc_orbmom_cal(compute_dmat,dimdmat,dmatpawu,dmatudiag,impose_dmat,indsym,my_natom,natom,&
3120 9 : & natpawu,nspinor,nsppol,nsym,ntypat,paw_ij,pawang,pawrad,pawprtvol,pawrhoij,pawtab,&
3121 9 : & spinat,symafm,typat,useexexch,usepawu,znucl, &
3122 9 : & mpi_atmtab,comm_atom,orb_mom_atom,maxmag,difmag) ! optional arguments (parallelism)
3123 :
3124 : !Arguments ---------------------------------------------
3125 : !scalars
3126 : integer,intent(in) :: compute_dmat,dimdmat,dmatudiag,impose_dmat,my_natom,natom,natpawu
3127 : integer,intent(in) :: nspinor,nsppol,nsym,ntypat,useexexch,usepawu
3128 : integer,optional,intent(in) :: comm_atom
3129 : type(pawang_type),intent(in) :: pawang
3130 : integer,intent(in) :: pawprtvol
3131 : !arrays
3132 : integer,intent(in) :: indsym(4,nsym,natom),symafm(nsym),typat(natom)
3133 : integer,optional,target,intent(in) :: mpi_atmtab(:)
3134 : real(dp),intent(in) :: dmatpawu(dimdmat,dimdmat,nspinor*nsppol,natpawu*impose_dmat)
3135 : real(dp),intent(in) :: spinat(3,natom)
3136 : type(paw_ij_type),intent(in) :: paw_ij(my_natom)
3137 : type(pawrhoij_type),intent(in) :: pawrhoij(my_natom)
3138 : type(pawtab_type),intent(in) :: pawtab(ntypat)
3139 : real(dp),intent(in) :: znucl(ntypat)
3140 9 : integer,pointer :: my_atmtab(:)
3141 : real(dp),intent(inout), optional :: orb_mom_atom(10,3,natom),maxmag,difmag
3142 18 : real(dp):: orb_mom_atom0(10,3,natom)
3143 : !Local variables ---------------------------------------
3144 : !scalars
3145 : logical :: paral_atom,my_atmtab_allocated
3146 : character(len=5) :: orb_char
3147 : integer :: cplex_dij,im1,im2,ndij,itypat,my_comm_atom
3148 : integer :: my_lcur,my_iatom,coor,isp,lmin,lmax,me_atom,mu
3149 9 : real(dp),allocatable :: my_l_occmat(:,:,:,:)
3150 9 : complex(dp),allocatable :: op_l(:,:,:),cmfoccmat(:,:,:)
3151 : real(dp) :: orb_mom(3)
3152 : real(dp) :: sum_orb_mom(3)
3153 : real(dp) :: orbmag_r, orbmag_theta, orbmag_phi
3154 : complex(dp) :: my_sls_val
3155 : character(len=500) :: message
3156 9 : type(paw_ij_type), ABI_CONTIGUOUS pointer :: paw_ij_all(:)
3157 9 : type(pawrhoij_type),ABI_CONTIGUOUS pointer :: pawrhoij_all(:)
3158 : type(atomdata_t) :: atom
3159 : !!****************************************************************
3160 9 : type(pawtab_type),allocatable :: pawtab_tmp(:)
3161 : type(pawrad_type) :: pawrad(ntypat)
3162 27 : real(dp) :: f6of2_sla(ntypat),f4of2_sla(ntypat),upawu(ntypat),jpawu(ntypat)
3163 18 : integer :: llexexch(ntypat),llpawu(ntypat),nn,ii
3164 : !*********************************************************************
3165 9 : orb_char='pdfgh'
3166 9 : if (.not. present(orb_mom_atom)) then
3167 1 : write(message,*) ' '
3168 3 : call wrtout([std_out, ab_out], message)
3169 1 : write(message,*) ' '
3170 3 : call wrtout([std_out, ab_out], message)
3171 1 : write(message,*) 'Integrated orbital magnetic moments inside the PAW spheres:'
3172 3 : call wrtout([std_out, ab_out], message)
3173 1 : write(message,*) '--------------------------------------------------------------------------------------------'
3174 3 : call wrtout([std_out, ab_out], message)
3175 1 : write(message,*) ' Atom orbital |orbmag| orbmag(x) orbmag(y) orbmag(z) orbmag(theta) orbmag(phi) '
3176 3 : call wrtout([std_out, ab_out], message)
3177 1 : write(message,*) '--------------------------------------------------------------------------------------------'
3178 3 : call wrtout([std_out, ab_out], message)
3179 : else
3180 1096 : orb_mom_atom0=orb_mom_atom
3181 : endif
3182 :
3183 :
3184 : !Set up parallelism over atoms
3185 :
3186 9 : paral_atom=(present(comm_atom).and.(my_natom/=natom))
3187 9 : nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
3188 9 : my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
3189 9 : call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,my_natom_ref=my_natom) !vz_d
3190 :
3191 :
3192 9 : if (paral_atom) then
3193 0 : me_atom=xmpi_comm_rank(my_comm_atom)
3194 : else
3195 : me_atom=0
3196 : end if
3197 :
3198 :
3199 :
3200 : !If atomic data are distributed, retrieve all paw_ij on master proc
3201 9 : if (paral_atom) then
3202 0 : if (me_atom==0) then
3203 0 : ABI_MALLOC(paw_ij_all,(natom))
3204 0 : call paw_ij_nullify(paw_ij_all)
3205 : else
3206 0 : ABI_MALLOC(paw_ij_all,(0))
3207 : end if
3208 0 : call paw_ij_gather(paw_ij,paw_ij_all,0,my_comm_atom)
3209 : else
3210 :
3211 63 : ABI_MALLOC(paw_ij_all,(natom))
3212 9 : call paw_ij_nullify(paw_ij_all)
3213 45 : paw_ij_all = paw_ij
3214 : end if
3215 :
3216 :
3217 : !If atomic data are distributed, retrieve all Rhoij on master proc
3218 9 : if (paral_atom) then
3219 0 : if (me_atom==0) then
3220 0 : ABI_MALLOC(pawrhoij_all,(natom))
3221 0 : call pawrhoij_nullify(pawrhoij_all)
3222 : else
3223 0 : ABI_MALLOC(pawrhoij_all,(0))
3224 : end if
3225 : call pawrhoij_gather(pawrhoij,pawrhoij_all,0,my_comm_atom,&
3226 : & with_grhoij=.false.,with_lmnmix=.false.,&
3227 0 : & with_rhoij_=.false.,with_rhoijres=.false.)
3228 :
3229 : else
3230 63 : ABI_MALLOC(pawrhoij_all,(natom))
3231 9 : call pawrhoij_nullify(pawrhoij_all)
3232 45 : pawrhoij_all = pawrhoij
3233 : end if
3234 :
3235 :
3236 9 : if (me_atom==0) then !!!!!!!!
3237 9 : nn=size(pawtab)
3238 45 : ABI_MALLOC(pawtab_tmp,(nn))
3239 9 : call pawtab_nullify(pawtab_tmp)
3240 :
3241 : ! ABI_MALLOC(pawtab_tmp,(nn))
3242 :
3243 9 : if (nn.gt.0) then
3244 27 : do ii=1,nn
3245 18 : pawtab_tmp(ii)%basis_size = pawtab(ii)%basis_size
3246 18 : pawtab_tmp(ii)%lmn_size = pawtab(ii)%lmn_size
3247 18 : pawtab_tmp(ii)%lmn2_size = pawtab(ii)%lmn2_size
3248 18 : pawtab_tmp(ii)%rpaw = pawtab(ii)%rpaw
3249 18 : pawtab_tmp(ii)%rshp = pawtab(ii)%rshp
3250 18 : pawtab_tmp(ii)%shape_type = pawtab(ii)%shape_type
3251 18 : if (allocated(pawtab(ii)%dij0)) then
3252 : ! ABI_MALLOC(pawtab_tmp(ii)%dij0,())
3253 18 : call alloc_copy(pawtab(ii)%dij0, pawtab_tmp(ii)%dij0)
3254 : end if
3255 :
3256 18 : if (allocated(Pawtab(ii)%indlmn)) then
3257 18 : call alloc_copy(pawtab(ii)%indlmn, pawtab_tmp(ii)%indlmn)
3258 : end if
3259 :
3260 18 : if (allocated(Pawtab(ii)%phi)) then
3261 18 : call alloc_copy(pawtab(ii)%phi, pawtab_tmp(ii)%phi)
3262 : end if
3263 :
3264 27 : if (allocated(Pawtab(ii)%indklmn)) then
3265 18 : call alloc_copy(pawtab(ii)%indklmn, pawtab_tmp(ii)%indklmn)
3266 : end if
3267 : end do
3268 :
3269 : ! pawtab_tmp%has_fock = Pawtab%has_fock !
3270 : ! pawtab_tmp%has_kij = Pawtab%has_kij !
3271 : ! pawtab_tmp%has_tproj = Pawtab%has_tproj !
3272 : ! pawtab_tmp%has_tvale = Pawtab%has_tvale
3273 : ! pawtab_tmp%has_coretau = Pawtab%has_coretau
3274 : ! pawtab_tmp%has_vhtnzc = Pawtab%has_vhtnzc
3275 : ! pawtab_tmp%has_vhnzc = Pawtab%has_vhnzc
3276 : ! pawtab_tmp%has_vminushalf = Pawtab%has_vminushalf
3277 : ! pawtab_tmp%has_nabla = Pawtab%has_nabla
3278 : ! pawtab_tmp%has_nablaphi =Pawtab%has_nablaphi
3279 : ! pawtab_tmp%has_shapefncg =Pawtab%has_shapefncg
3280 : ! pawtab_tmp%has_wvl =Pawtab%has_wvl
3281 :
3282 : ! pawtab_tmp%usetcore =Pawtab%usetcore
3283 : ! pawtab_tmp%usexcnhat =Pawtab%usexcnhat
3284 : ! pawtab_tmp%useexexch =Pawtab%useexexch
3285 36 : pawtab_tmp%usepawu = Pawtab%usepawu
3286 : ! pawtab_tmp%usepotzero =Pawtab%usepotzero
3287 36 : pawtab_tmp%usespnorb =Pawtab%usespnorb
3288 : ! pawtab_tmp%mqgrid =Pawtab%mqgrid
3289 : ! pawtab_tmp%mqgrid_shp =Pawtab%mqgrid_shp
3290 :
3291 36 : pawtab_tmp%basis_size =Pawtab%basis_size
3292 36 : pawtab_tmp%ij_proj =Pawtab%ij_proj
3293 36 : pawtab_tmp%ij_size =Pawtab%ij_size
3294 36 : pawtab_tmp%lcut_size =Pawtab%lcut_size
3295 36 : pawtab_tmp%l_size =Pawtab%l_size
3296 36 : pawtab_tmp%lexexch =Pawtab%lexexch
3297 36 : pawtab_tmp%lmn_size =Pawtab%lmn_size
3298 36 : pawtab_tmp%lmn2_size =Pawtab%lmn2_size
3299 36 : pawtab_tmp%lmnmix_sz =Pawtab%lmnmix_sz
3300 36 : pawtab_tmp%lpawu =Pawtab%lpawu
3301 36 : pawtab_tmp%nproju =Pawtab%nproju
3302 36 : pawtab_tmp%option_interaction_pawu =Pawtab%option_interaction_pawu
3303 36 : pawtab_tmp%mesh_size =Pawtab%mesh_size
3304 36 : pawtab_tmp%partialwave_mesh_size =Pawtab%partialwave_mesh_size
3305 36 : pawtab_tmp%core_mesh_size =Pawtab%core_mesh_size
3306 36 : pawtab_tmp%coretau_mesh_size =Pawtab%coretau_mesh_size
3307 36 : pawtab_tmp%vminus_mesh_size =Pawtab%vminus_mesh_size
3308 36 : pawtab_tmp%tnvale_mesh_size =Pawtab%tnvale_mesh_size
3309 36 : pawtab_tmp%shape_type =Pawtab%shape_type
3310 :
3311 : end if
3312 :
3313 9 : sum_orb_mom=zero
3314 :
3315 27 : f4of2_sla = pawtab(:)%f4of2_sla
3316 27 : f6of2_sla = pawtab(:)%f6of2_sla
3317 27 : llexexch= -1*one !(/-1,-1/)
3318 :
3319 9 : lmin=1
3320 45 : do my_iatom=1,natom
3321 36 : itypat=pawrhoij_all(my_iatom)%itypat
3322 36 : lmax=((pawtab(itypat)%l_size)-1)/2
3323 :
3324 108 : do my_lcur=lmin,lmax
3325 :
3326 72 : if (my_lcur == pawtab(itypat)%lpawu) then
3327 54 : upawu= pawtab(:)%upawu
3328 54 : jpawu= pawtab(:)%jpawu
3329 54 : llpawu= pawtab(:)%lpawu
3330 : else
3331 162 : upawu= zero
3332 162 : jpawu= zero
3333 162 : llpawu= one
3334 : end if
3335 :
3336 : call pawpuxinit(2,pawtab(1)%exchmix ,f4of2_sla,f6of2_sla,.False.,jpawu ,llexexch,llpawu,&
3337 : & nspinor,ntypat,pawtab_tmp(itypat)%option_interaction_pawu,pawang,pawprtvol,pawrad,pawtab_tmp,upawu,0,&
3338 : & useexexch,usepawu,&
3339 72 : & 0,lmagCalc=.True.)
3340 72 : cplex_dij=paw_ij_all(my_iatom)%cplex_dij
3341 72 : ndij=paw_ij_all(my_iatom)%ndij
3342 :
3343 :
3344 360 : ABI_MALLOC(op_l,(2*my_lcur+1,2*my_lcur+1,3))
3345 432 : ABI_MALLOC(my_l_occmat,(cplex_dij,2*my_lcur+1,2*my_lcur+1,ndij))
3346 360 : ABI_MALLOC(cmfoccmat,(2*my_lcur+1,2*my_lcur+1,ndij))
3347 :
3348 : call setnoccmmp(compute_dmat,dimdmat,dmatpawu,dmatudiag,impose_dmat,indsym,natom,natom,&
3349 : & natpawu,nspinor,nsppol,nsym,ntypat,paw_ij_all,pawang,pawprtvol,pawrhoij_all,pawtab_tmp,&
3350 : & spinat,symafm,typat,useexexch,usepawu, &
3351 144 : & mpi_atmtab,comm_atom,l_orbmom=my_lcur,atom_orbmom=my_iatom,my_l_occmat=my_l_occmat)
3352 :
3353 6408 : cmfoccmat(:,:,:)=cmplx(my_l_occmat(1,:,:,:),my_l_occmat(2,:,:,:))
3354 :
3355 72 : my_sls_val=zero
3356 72 : orb_mom=zero
3357 4824 : op_l=czero
3358 :
3359 288 : do coor=1,3
3360 1152 : do im1=1,2*my_lcur+1
3361 4752 : do im2=1,2*my_lcur+1
3362 3672 : call slxyzs(my_lcur,im2-(my_lcur+1),coor,my_lcur,im1-(my_lcur+1),my_sls_val)
3363 4536 : op_l(im1,im2,coor)=my_sls_val
3364 : end do
3365 : end do
3366 : end do
3367 :
3368 :
3369 : !!! calculating for spin upup and dndn
3370 :
3371 216 : do isp=1,2 ! upup and dndn spin components are considered
3372 648 : do coor=1,3
3373 2304 : do im1=1,2*my_lcur+1
3374 9504 : do im2=1,2*my_lcur+1
3375 9072 : orb_mom(coor)=orb_mom(coor)+op_l(im2,im1,coor)*conjg(cmfoccmat(im1,im2,isp))
3376 : end do
3377 : end do
3378 : end do
3379 : end do
3380 :
3381 72 : if (allocated(op_l)) then
3382 72 : ABI_FREE(op_l)
3383 : end if
3384 72 : if (allocated(my_l_occmat)) then
3385 72 : ABI_FREE(my_l_occmat)
3386 : end if
3387 72 : if (allocated(cmfoccmat)) then
3388 72 : ABI_FREE(cmfoccmat)
3389 : end if
3390 :
3391 288 : sum_orb_mom=sum_orb_mom+orb_mom
3392 :
3393 72 : call vcart2ylm(orb_mom, orbmag_r, orbmag_theta, orbmag_phi)
3394 72 : call atomdata_from_znucl(atom, znucl(typat(my_iatom)))
3395 108 : if (.not. present(orb_mom_atom)) then
3396 8 : if (my_lcur==1) then
3397 4 : write(message,'(i5,a3,a8,4f12.6,2f13.6)') my_iatom,atom%symbol, orb_char(my_lcur:my_lcur),orbmag_r, orb_mom(1),orb_mom(2),orb_mom(3),orbmag_theta,orbmag_phi
3398 : else
3399 4 : write(message,'(a5,a11,4f12.6,2f13.6)') '', orb_char(my_lcur:my_lcur), orbmag_r,orb_mom(1),orb_mom(2),orb_mom(3) ,orbmag_theta,orbmag_phi
3400 : end if
3401 :
3402 24 : call wrtout([std_out, ab_out], message)
3403 : else
3404 256 : orb_mom_atom(my_lcur,:,my_iatom)=orb_mom
3405 256 : do mu=1,3
3406 192 : maxmag=max(maxmag,abs(orb_mom_atom(my_lcur,mu,my_iatom)))
3407 256 : difmag=max(difmag,abs(orb_mom_atom(my_lcur,mu,my_iatom)-orb_mom_atom0(my_lcur,mu,my_iatom)))
3408 : enddo
3409 : endif
3410 :
3411 :
3412 : end do !!!!!!!!! END DO lcur
3413 45 : if (.not. present(orb_mom_atom)) then
3414 4 : write(message,*) '--------------------------------------------------------------------------------------------'
3415 12 : call wrtout([std_out, ab_out], message)
3416 : endif
3417 : end do !!!!!!!!! END DO natoms
3418 :
3419 :
3420 9 : call vcart2ylm(sum_orb_mom(1:3), orbmag_r, orbmag_theta, orbmag_phi)
3421 9 : if (.not. present(orb_mom_atom)) then
3422 1 : write(message,'(a,3f12.6)') ' Sum (cart.coord.) ', sum_orb_mom(1),sum_orb_mom(2),sum_orb_mom(3)
3423 3 : call wrtout([std_out, ab_out], message)
3424 1 : write(message,'(a,1f11.6,a,2f13.6)') ' (sph.coord.)', orbmag_r,' ',orbmag_theta, orbmag_phi
3425 3 : call wrtout([std_out, ab_out], message)
3426 1 : write(message,*) '--------------------------------------------------------------------------------------------'
3427 3 : call wrtout([std_out, ab_out], message)
3428 1 : write(message,*) ' '
3429 3 : call wrtout([std_out, ab_out], message)
3430 : endif
3431 :
3432 : end if !!!!!!!!!!
3433 :
3434 9 : if (allocated(pawtab_tmp)) then
3435 9 : call pawtab_nullify(pawtab_tmp)
3436 9 : call pawtab_free(pawtab_tmp)
3437 27 : ABI_FREE(pawtab_tmp)
3438 : end if
3439 :
3440 9 : if (paral_atom) then
3441 0 : if (me_atom==0) then
3442 :
3443 0 : do my_iatom=1,natom
3444 0 : if (allocated(paw_ij_all(my_iatom)%noccmmp)) then
3445 0 : ABI_FREE(paw_ij_all(my_iatom)%noccmmp)
3446 : end if
3447 0 : if (allocated(paw_ij_all(my_iatom)%nocctot)) then
3448 0 : ABI_FREE(paw_ij_all(my_iatom)%nocctot)
3449 : end if
3450 : end do
3451 : end if
3452 :
3453 0 : call pawrhoij_nullify(pawrhoij_all)
3454 0 : call pawrhoij_free(pawrhoij_all)
3455 0 : ABI_FREE(pawrhoij_all)
3456 :
3457 0 : call paw_ij_nullify(paw_ij_all)
3458 0 : call paw_ij_free(paw_ij_all)
3459 0 : ABI_FREE(paw_ij_all)
3460 : else
3461 :
3462 9 : call pawrhoij_nullify(pawrhoij_all)
3463 9 : ABI_FREE(pawrhoij_all)
3464 :
3465 9 : call paw_ij_nullify(paw_ij_all)
3466 9 : ABI_FREE(paw_ij_all)
3467 : end if
3468 :
3469 27 : end subroutine loc_orbmom_cal
3470 : !----------------------------------------------------------------------
3471 : END MODULE m_paw_correlations
3472 : !!***
|
