Branch data 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 : 1525 : subroutine pawpuxinit(dmatpuopt,exchmix,f4of2_sla,f6of2_sla,is_dfpt,jpawu,llexexch,llpawu,&
123 : 1525 : & 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 : 1525 : integer,contiguous, pointer :: indlmn(:,:)
167 : : real(dp) :: euijkl_temp(3),euijkl_temp2(3),euijkl_dc(3)
168 : 1525 : real(dp),allocatable :: ff(:),fk(:),gg(:)
169 : 1525 : type(pawrad_type) :: pawrad_tmp
170 : : ! *************************************************************************
171 : :
172 : : DBG_ENTER("COLL")
173 : 1525 : Loc_prtvol = 3
174 : 1525 : lmagCalc_ = .False.
175 : 1525 : 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 [ + + + + : 1525 : if(useexexch==0.and.usepawu==0.and.use_dmft==0) then
+ - ]
182 [ + + ]: 3138 : do itypat=1,ntypat
183 : 3138 : 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 [ + + - + ]: 237 : 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 [ + + ]: 238 : & write(message, '(3a)' ) ch10,ch10," ******************************************"
205 : : if(usepawu==1) then
206 : 145 : 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 [ + + ]: 237 : if(useexexch/=0) write(message, '(3a)' ) trim(message),ch10," PAW Local Exact exchange: PBE0"
222 [ + + + + : 237 : if((abs(usepawu)>=1.and.abs(usepawu)<=4).or.useexexch/=0 .and.(.not.lmagCalc_)) then
+ - ]
223 [ + + ]: 185 : if (nspinor==2) then
224 : 99 : write(message, '(3a,i1)' ) trim(message),ch10," Magnetic DC : option_interaction = ",option_interaction
225 : : end if
226 : 185 : write(message, '(3a)' ) trim(message),ch10," ******************************************"
227 : : end if
228 [ + + + - : 237 : if(use_dmft==0 .and. abs(usepawu)<=4 .and.(.not.lmagCalc_)) then
+ + ]
229 : 113 : call wrtout(ab_out,message,'COLL')
230 : 113 : 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 : 237 : option_interaction_ = option_interaction
237 [ + + + + : 237 : 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 [ + + + + : 237 : 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 [ + + ]: 641 : do itypat=1,ntypat
249 : 404 : indlmn => pawtab(itypat)%indlmn
250 : 404 : lmn_size=pawtab(itypat)%lmn_size
251 : 404 : lmn2_size=pawtab(itypat)%lmn2_size
252 : 404 : mesh_size=pawtab(itypat)%mesh_size
253 : 404 : int_meshsz=pawrad(itypat)%int_meshsz
254 : 404 : lcur=-1
255 : :
256 [ + + ]: 404 : 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 [ + + - + ]: 404 : if (usepawu/=0.or.use_dmft>0) then
266 : 402 : lcur=llpawu(itypat)
267 : 402 : pawtab(itypat)%lpawu=lcur
268 [ + + ]: 402 : if(lcur/=-1) then
269 : 302 : pawtab(itypat)%usepawu=usepawu
270 : 302 : pawtab(itypat)%upawu=upawu(itypat)
271 : 302 : pawtab(itypat)%jpawu=jpawu(itypat)
272 : 302 : pawtab(itypat)%f4of2_sla=f4of2_sla(itypat)
273 : 302 : pawtab(itypat)%f6of2_sla=f6of2_sla(itypat)
274 : 302 : pawtab(itypat)%option_interaction_pawu=option_interaction_
275 : : else
276 : 100 : pawtab(itypat)%usepawu=0
277 : 100 : pawtab(itypat)%upawu=zero
278 : 100 : pawtab(itypat)%jpawu=zero
279 : 100 : pawtab(itypat)%f4of2_sla=zero
280 : 100 : pawtab(itypat)%f6of2_sla=zero
281 : 100 : pawtab(itypat)%option_interaction_pawu=option_interaction_
282 : : end if
283 : : end if
284 : :
285 : : ! Local exact-echange data
286 [ + + ]: 404 : 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 [ + + ]: 641 : if(lcur/=-1) then
296 : :
297 : : ! Compute number of projectors for DFT+U/local exact-exchange/DFT+DMFT
298 [ + + + + ]: 5544 : icount=count(indlmn(1,1:lmn_size)==lcur)
299 : 303 : pawtab(itypat)%nproju=icount/(2*lcur+1)
300 [ + + - + ]: 303 : 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 [ - + ]: 303 : 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 [ + + ]: 303 : if ((.not.lmagCalc_)) then
311 : 177 : write(message, '(a,a,i4,a,a,i4)' ) ch10,&
312 : 177 : & ' pawpuxinit : for species ',itypat,ch10,&
313 : 354 : & ' number of projectors is',pawtab(itypat)%nproju
314 : 177 : call wrtout(std_out,message,'COLL')
315 : :
316 : : end if
317 : 303 : pawtab(itypat)%ij_proj=pawtab(itypat)%nproju*(pawtab(itypat)%nproju+1)/2
318 : :
319 : : ! ==================================================
320 : : ! A-define useful indexes
321 : : ! --------------------------------------------------
322 [ + + ]: 303 : if (allocated(pawtab(itypat)%lnproju)) then
323 : 204 : ABI_FREE(pawtab(itypat)%lnproju)
324 : : end if
325 [ - + + - : 909 : ABI_MALLOC(pawtab(itypat)%lnproju,(pawtab(itypat)%nproju))
+ - + - -
+ - + -
+ ]
326 : 303 : icount=0
327 [ + + ]: 5544 : do ilmn=1,lmn_size
328 [ + + ]: 5544 : if(indlmn(1,ilmn)==lcur) then
329 : 2361 : itemp=icount/(2*lcur+1)
330 [ + + ]: 2361 : if (itemp*(2*lcur+1)==icount) then
331 : 545 : pawtab(itypat)%lnproju(itemp+1)=indlmn(5,ilmn)
332 : : end if
333 : 2361 : icount=icount+1
334 : : end if
335 : : end do
336 : :
337 [ + + ]: 303 : if (allocated(pawtab(itypat)%klmntomn)) then
338 : 204 : ABI_FREE(pawtab(itypat)%klmntomn)
339 : : end if
340 [ - + + - : 909 : ABI_MALLOC(pawtab(itypat)%klmntomn,(4,lmn2_size))
+ - + - -
+ - + -
+ ]
341 [ + + ]: 5544 : do jlmn=1,lmn_size
342 : 5241 : jl= indlmn(1,jlmn)
343 : 5241 : j0lmn=jlmn*(jlmn-1)/2
344 [ + + ]: 57968 : do ilmn=1,jlmn
345 : 52424 : il= indlmn(1,ilmn)
346 : 52424 : klmn=j0lmn+ilmn
347 : 52424 : pawtab(itypat)%klmntomn(1,klmn)=indlmn(2,ilmn)+il+1
348 : 52424 : pawtab(itypat)%klmntomn(2,klmn)=indlmn(2,jlmn)+jl+1
349 : 52424 : pawtab(itypat)%klmntomn(3,klmn)=indlmn(3,ilmn)
350 : 57665 : 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 [ + + + - ]: 303 : 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 [ + + + + ]: 235 : else if(dmatpuopt==2 .and.(.not.lmagCalc_)) then
369 : 85 : write(message, '(6a)' ) ch10,&
370 : 85 : & ' pawpuxinit : dmatpuopt=2 ',ch10,&
371 : 85 : & ' PAW+U: dens. mat. constructed by selecting contribution',ch10,&
372 : 170 : & ' for each angular momentum to the density (inside PAW augm. region(s))'
373 : 85 : 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 [ - + + - : 909 : ABI_MALLOC(ff,(mesh_size))
+ - + - -
+ - + -
+ ]
388 [ + + ]: 191896 : ff(:)=zero
389 : :
390 [ + + ]: 303 : if (allocated(pawtab(itypat)%ph0phiint)) then
391 : 204 : ABI_FREE(pawtab(itypat)%ph0phiint)
392 : : end if
393 [ + + ]: 303 : if (allocated(pawtab(itypat)%zioneff)) then
394 : 204 : ABI_FREE(pawtab(itypat)%zioneff)
395 : : end if
396 [ - + + - : 909 : ABI_MALLOC(pawtab(itypat)%ph0phiint,(pawtab(itypat)%nproju))
+ - + - -
+ - + -
+ ]
397 [ - + + - : 606 : ABI_MALLOC(pawtab(itypat)%zioneff,(pawtab(itypat)%nproju))
- + - + -
+ ]
398 : :
399 : 303 : icount=0
400 [ + + ]: 848 : 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 : 545 : ju=pawtab(itypat)%mesh_size-1
413 : 545 : ak=pawtab(itypat)%phi(ju,pawtab(itypat)%lnproju(iu))/pawtab(itypat)%phi(ju+1,pawtab(itypat)%lnproju(iu))
414 : 545 : ak=ak*(pawrad(itypat)%rad(ju+1)/pawrad(itypat)%rad(ju))**(pawtab(itypat)%lpawu-1)
415 : 545 : 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 [ + + ]: 362381 : & *pawtab(itypat)%phi(1:mesh_size,pawtab(itypat)%lnproju(iu))
419 : 545 : call simp_gen(int1,ff,pawrad(itypat))
420 : 848 : pawtab(itypat)%ph0phiint(iu)=int1
421 : : end do
422 [ + + ]: 303 : 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 [ + + ]: 303 : if (allocated(pawtab(itypat)%phiphjint)) then
436 : 204 : ABI_FREE(pawtab(itypat)%phiphjint)
437 : : end if
438 [ - + + - : 909 : ABI_MALLOC(pawtab(itypat)%phiphjint,(pawtab(itypat)%ij_proj))
+ - + - -
+ - + -
+ ]
439 : :
440 : 303 : icount=0
441 [ + + ]: 848 : do ju=1,pawtab(itypat)%nproju
442 [ + + ]: 1635 : do iu=1,ju
443 : 787 : icount=icount+1
444 [ + + + - ]: 1332 : if ((dmatpuopt==1).and.(useexexch==0)) then
445 : : pawtab(itypat)%phiphjint(icount)=pawtab(itypat)%ph0phiint(iu)*&
446 : 144 : & pawtab(itypat)%ph0phiint(ju)
447 [ + + - + ]: 643 : else if((dmatpuopt==2).or.(useexexch/=0)) then
448 : : ff(1:mesh_size)=pawtab(itypat)%phi(1:mesh_size,pawtab(itypat)%lnproju(iu))&
449 [ + + ]: 334866 : & *pawtab(itypat)%phi(1:mesh_size,pawtab(itypat)%lnproju(ju))
450 : 571 : call simp_gen(int1,ff,pawrad(itypat))
451 : 571 : 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 [ - + ]: 303 : 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 [ - + ]: 303 : ABI_FREE(ff)
468 [ + + ]: 303 : 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 [ + + ]: 303 : if (usepawu/=0) then
482 : 302 : lpawu=lcur
483 : :
484 [ + + ]: 302 : if (allocated(pawtab(itypat)%vee)) then
485 : 204 : ABI_FREE(pawtab(itypat)%vee)
486 : : end if
487 : 302 : sz1=2*lpawu+1
488 [ - + + - : 1812 : 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 : 302 : & 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 [ + + ]: 302 : if (.not.lmagCalc_) then
502 : 176 : write(message,'(a)') ch10
503 : 176 : call wrtout(std_out,message,'COLL')
504 : 176 : write(message,'(a)') " Matrix of interaction vee(m1,m2,m1,m2)"
505 : 176 : call wrtout(std_out,message,'COLL')
506 [ + + ]: 1090 : do m1=1,2*lpawu+1
507 [ + - + + ]: 5770 : write(message,'(2x,14(f20.14,2x))') (pawtab(itypat)%vee(m1,m2,m1,m2)*Ha_eV,m2=1,2*lpawu+1)
508 : 1090 : 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 : 176 : write(message,'(a)') ch10
514 : 176 : 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 [ + + + + ]: 302 : compute_euijkl=(is_dfpt.or.usepawu<0)
539 : : if (compute_euijkl) then
540 [ + + ]: 12 : if (allocated(pawtab(itypat)%euijkl)) then
541 : 7 : ABI_FREE(pawtab(itypat)%euijkl)
542 : : end if
543 [ - + + - : 72 : ABI_MALLOC(pawtab(itypat)%euijkl,(3,lmn_size,lmn_size,lmn_size,lmn_size))
- + + - -
+ + - - +
+ - + - +
- - + - +
- + ]
544 [ + + + + : 4494038 : pawtab(itypat)%euijkl = zero
+ + + + +
+ ]
545 : 12 : compute_euij_fll = .false.
546 : 12 : euijkl_temp2=zero
547 [ + + ]: 12 : if (abs(usepawu)==1.or.abs(usepawu)==4) then ! Only for FLL
548 [ + + ]: 10 : if (allocated(pawtab(itypat)%euij_fll)) then ! allocate euij_fll for FLL
549 : 5 : ABI_FREE(pawtab(itypat)%euij_fll)
550 : : end if
551 [ - + + - : 30 : ABI_MALLOC(pawtab(itypat)%euij_fll,(lmn2_size))
+ - + - -
+ - + -
+ ]
552 [ + + ]: 1560 : pawtab(itypat)%euij_fll = zero
553 : : compute_euij_fll = .true.
554 : : end if
555 : :
556 : : ! loop on i,j
557 [ + + ]: 1904 : do klmna=1,lmn2_size
558 : 1892 : ilmn=pawtab(itypat)%indklmn(7,klmna) ! i
559 : 1892 : jlmn=pawtab(itypat)%indklmn(8,klmna) ! j
560 [ + + + - ]: 1904 : if (pawtab(itypat)%indlmn(1,ilmn)==lpawu.and.pawtab(itypat)%indlmn(1,jlmn)==lpawu) then ! only correlated orbitals
561 : 580 : iu = pawtab(itypat)%indlmn(3,ilmn) ! ni
562 : 580 : ju = pawtab(itypat)%indlmn(3,jlmn) ! nj
563 : 580 : phiint_ij = pawtab(itypat)%phiphjint(iu+(ju*(ju-1))/2) ! iu <= ju by construction (ilmn<=jlmn)
564 : 580 : m2 = pawtab(itypat)%indlmn(2,ilmn) ! mi
565 : 580 : m21=m2+lpawu+1
566 : 580 : m1 = pawtab(itypat)%indlmn(2,jlmn) ! mj
567 : 580 : m11=m1+lpawu+1
568 : :
569 [ + + ]: 580 : if (compute_euij_fll.and.m1==m2) then ! FLL
570 : 130 : 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 [ + + ]: 97360 : do klmnb=1,lmn2_size
575 : 96780 : ilmnp=pawtab(itypat)%indklmn(7,klmnb) ! ip (=k)
576 : 96780 : jlmnp=pawtab(itypat)%indklmn(8,klmnb) ! jp (=l)
577 [ + + + - ]: 97360 : if (pawtab(itypat)%indlmn(1,ilmnp)==lpawu.and.pawtab(itypat)%indlmn(1,jlmnp)==lpawu) then ! correlated orbitals
578 : 30700 : iup = pawtab(itypat)%indlmn(3,ilmnp) ! nip
579 : 30700 : jup = pawtab(itypat)%indlmn(3,jlmnp) ! njp
580 : 30700 : phiint_ipjp = pawtab(itypat)%phiphjint(iup+(jup*(jup-1))/2) ! iup <= jup by construction (ilmnp<=jlmnp)
581 : 30700 : m4 = pawtab(itypat)%indlmn(2,ilmnp) ! mip
582 : 30700 : m41=m4+lpawu+1
583 : 30700 : m3 = pawtab(itypat)%indlmn(2,jlmnp) ! mjp
584 : 30700 : m31=m3+lpawu+1
585 : :
586 : 30700 : 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 [ + + ]: 30700 : if (m1==m2.and.m3==m4) then ! In that case, we have to add the double-counting term
590 : :
591 [ + + + + ]: 2300 : if (abs(usepawu)==1.and.nspinor==1) then ! FLL
592 : :
593 : : euijkl_dc(1) = &
594 : 475 : & phiint_ij * phiint_ipjp * ( pawtab(itypat)%upawu - pawtab(itypat)%jpawu )
595 : : euijkl_dc(2) = &
596 : 475 : & 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 : 30700 : euijkl_temp(:) = zero
628 : 30700 : euijkl_temp2(:) = zero
629 : :
630 : 30700 : 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 [ + + ]: 92100 : euijkl_temp(1:2) = phiint_ij * phiint_ipjp * vee1
638 : :
639 : 30700 : 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 : 30700 : euijkl_temp2(1) = phiint_ij * phiint_ipjp * vee2
650 : 30700 : euijkl_temp2(3) = phiint_ij * phiint_ipjp * vee2
651 [ + + ]: 122800 : pawtab(itypat)%euijkl(:,ilmn,jlmn,ilmnp,jlmnp) = euijkl_temp(:) - euijkl_temp2(:) - euijkl_dc(:)
652 [ + + ]: 122800 : 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 : 30700 : vee2 = pawtab(itypat)%vee(m21,m31,m41,m11)
656 : 30700 : euijkl_temp2(1) = phiint_ij * phiint_ipjp * vee2
657 : 30700 : euijkl_temp2(3) = phiint_ij * phiint_ipjp * vee2
658 [ + + ]: 122800 : 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 : 30700 : vee2 = pawtab(itypat)%vee(m11,m41,m31,m21)
662 : 30700 : euijkl_temp2(1) = phiint_ij * phiint_ipjp * vee2
663 : 30700 : euijkl_temp2(3) = phiint_ij * phiint_ipjp * vee2
664 [ + + ]: 122800 : 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 [ + + ]: 303 : 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 [ + + ]: 303 : if (present(ucrpa)) then
819 [ - + ]: 293 : 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 [ + + ]: 303 : 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 [ + + - + : 1762 : 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 : 302 : 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 : 302 : real(dp),allocatable :: fk(:)
1026 : :
1027 : : ! *************************************************************************
1028 : :
1029 : : DBG_ENTER("COLL")
1030 : :
1031 : :
1032 : 302 : prtvol_ = 3
1033 [ + - ]: 302 : if (present(prtvol)) then
1034 : 302 : prtvol_ = prtvol
1035 : : end if
1036 : : ! Select only atoms with +U
1037 [ + - ]: 302 : 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 [ + - - + : 906 : ABI_MALLOC(fk,(lpawu+1))
+ - + - +
- - + -
+ ]
1046 : 302 : 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 [ + - ]: 302 : if(lpawu==0) then
1050 : : fk(1)=fk(1)
1051 [ + + ]: 302 : else if(lpawu==1) then
1052 : 114 : fk(2)=jpawu*5._dp
1053 [ + + ]: 188 : else if(lpawu==2) then
1054 : : ! f4of2=0._dp
1055 [ + + ]: 165 : if(f4of2_sla<-0.1_dp) then
1056 : 147 : f4of2=0.625_dp
1057 : 147 : f4of2_sla=f4of2
1058 : : else
1059 : : f4of2=f4of2_sla
1060 : : end if
1061 : 165 : fk(2)=jpawu*14._dp/(One+f4of2)
1062 : 165 : fk(3)=fk(2)*f4of2
1063 [ + + ]: 165 : if(abs(prtvol_)>=2) then
1064 : 147 : write(message,'(a,3x,a,f9.4,f9.4,f9.4,f9.4)') ch10,&
1065 : 294 : & "Slater parameters F^0, F^2, F^4 are",fk(1),fk(2),fk(3)
1066 : 147 : 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 [ + + + + : 207082 : vee=zero
+ + + + ]
1105 : 302 : lmpawu=(lpawu-1)**2+2*(lpawu-1)+1 ! number of m value below correlated orbitals
1106 : 302 : klm0u=lmpawu*(lmpawu+1)/2 ! value of klmn just below correlated orbitals
1107 : : ! --------- 4 loops for interaction matrix
1108 [ + + ]: 1630 : do m1=-lpawu,lpawu
1109 : 1328 : m11=m1+lpawu+1
1110 [ + + ]: 5433 : do m2=-lpawu,m1
1111 : 3803 : 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 : 3803 : klma=klm0u+m11*lmpawu+(m11-1)*m11/2+m21
1117 [ + + ]: 24066 : do m3=-lpawu,lpawu
1118 : 18935 : m31=m3+lpawu+1
1119 [ + + ]: 81999 : do m4=-lpawu,m3
1120 : 59261 : m41=m4+lpawu+1
1121 : 59261 : klmb=klm0u+m31*lmpawu+(m31-1)*m31/2+m41
1122 : : ! --------- loop on k=1,2,3 (4 if f orbitals)
1123 [ + - ]: 59261 : do kyc=1,2*lpawu+1,2
1124 : 191711 : lkyc=kyc-1
1125 : 191711 : lmkyc=(lkyc+1)*(lkyc)+1
1126 : 191711 : ak=zero
1127 [ + + ]: 1278106 : do mkyc=-lkyc,lkyc,1
1128 : 1086395 : isela=pawang%gntselect(lmkyc+mkyc,klma)
1129 : 1086395 : iselb=pawang%gntselect(lmkyc+mkyc,klmb)
1130 [ + + ]: 1278106 : 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 : 191711 : ak=ak/(two*dble(lkyc)+one)
1134 [ + + ]: 191711 : vee(m11,m31,m21,m41)=ak*fk(lkyc/2+1)+vee(m11,m31,m21,m41)
1135 : : end do !kyc
1136 : 59261 : vee(m11,m31,m21,m41)=vee(m11,m31,m21,m41)*four_pi
1137 : 59261 : vee(m21,m31,m11,m41)=vee(m11,m31,m21,m41)
1138 : 59261 : vee(m11,m41,m21,m31)=vee(m11,m31,m21,m41)
1139 : 78196 : vee(m21,m41,m11,m31)=vee(m11,m31,m21,m41)
1140 : : end do
1141 : : end do
1142 : : end do
1143 : : end do
1144 : 302 : ABI_FREE(fk)
1145 : : endif
1146 : :
1147 : 302 : 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 [ - + - + ]: 2787 : 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 : 2787 : real(dp),allocatable :: n12_sig(:),n34_msig(:),n34_sig(:)
1198 : : character(len=4) :: tag
1199 : : character(len=500) :: message
1200 : : ! *****************************************************
1201 : :
1202 : 2787 : nspden=size(nocctot)
1203 : 2787 : cplex_occ=size(noccmmp,1)
1204 : :
1205 [ - + ]: 2787 : 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 [ - + ]: 2787 : if (pawtab%usepawu<0) then
1210 : 0 : message='not allowed for usepawu<0!'
1211 : 0 : ABI_BUG(message)
1212 : : end if
1213 [ + + ]: 2787 : 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 : 2787 : lpawu=pawtab%lpawu
1226 [ + + ]: 2787 : upawu=pawtab%upawu;if(present(u_dmft)) upawu=u_dmft
1227 [ + + ]: 2787 : jpawu=pawtab%jpawu;if(present(j_dmft)) jpawu=j_dmft
1228 : :
1229 : : !======================================================
1230 : : !Compute DFT+U Energy
1231 : : !-----------------------------------------------------
1232 : :
1233 : 2787 : edftutemp=zero
1234 : 2787 : edcdc_opt3=zero
1235 : 2787 : eks_opt3=zero
1236 : :
1237 [ - + + - : 8361 : ABI_MALLOC(n12_sig,(cplex_occ))
+ - + - -
+ - + ]
1238 [ - + + - : 5574 : ABI_MALLOC(n34_msig,(cplex_occ))
- + - + ]
1239 [ - + + - : 5574 : ABI_MALLOC(n34_sig,(cplex_occ))
- + - + ]
1240 [ + + ]: 8011 : do ispden=1,min(nspden,2)
1241 : 5224 : jspden=min(nspden,2)-ispden+1
1242 : :
1243 : : ! Compute n_sigs and n_msigs for pawtab%usepawu=3
1244 [ + + ]: 5224 : if (nspden<=2) then
1245 : 3954 : n_sig =nocctot(ispden)
1246 : 3954 : n_msig=nocctot(jspden)
1247 : 3954 : 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 : 5224 : n_sigs =n_sig/(float(2*lpawu+1))
1267 : 5224 : 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 [ + + ]: 34117 : do m1=-lpawu,lpawu
1286 : 26106 : m11=m1+lpawu+1
1287 [ + + ]: 165346 : do m2=-lpawu,lpawu
1288 : 134016 : m21=m2+lpawu+1
1289 [ + + ]: 311481 : n12_sig(:)=noccmmp(:,m11,m21,ispden)
1290 [ + + + - : 134016 : if(m21==m11.and.(pawtab%usepawu==3.or.dmftdc==3)) n12_sig(1)=n12_sig(1)-n_sigs
- + ]
1291 [ + + ]: 865356 : do m3=-lpawu,lpawu
1292 : 705234 : m31=m3+lpawu+1
1293 [ + + ]: 4643754 : do m4=-lpawu,lpawu
1294 : 3804504 : m41=m4+lpawu+1
1295 [ + + ]: 8814009 : n34_sig(:) =noccmmp(:,m31,m41,ispden)
1296 [ + + ]: 8814009 : n34_msig(:)=noccmmp(:,m31,m41,jspden)
1297 [ + + + - : 3804504 : 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 : 3804504 : & + n12_sig(1)*n34_sig(1) *(pawtab%vee(m11,m31,m21,m41)-pawtab%vee(m11,m31,m41,m21))
1304 [ + + ]: 3804504 : 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 [ + - - + ]: 4509738 : 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 [ + + ]: 2787 : if (nspden==1) edftutemp=two*edftutemp ! Non-magn. system: sum up and dn energies
1329 [ - + ]: 2787 : ABI_FREE(n12_sig)
1330 [ - + ]: 2787 : ABI_FREE(n34_msig)
1331 [ - + ]: 2787 : ABI_FREE(n34_sig)
1332 : :
1333 : : !Non-collinear magnetism: add non-diagonal term; see (Eq 3) in PRB 72, 024458 (2005) [[cite:Shurikov2005]]
1334 [ + + ]: 2787 : 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 : 2787 : 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 : 2787 : n_upup=n_sig
1386 : 2787 : n_dndn=n_msig
1387 : :
1388 : 2787 : edcdctemp=zero;edctemp=zero
1389 : :
1390 : : !Full localized limit
1391 [ + + + + : 2787 : if((pawtab%usepawu==1.or.pawtab%usepawu==4).or.(dmftdc==1.or.dmftdc==4.or.dmftdc==5)) then
+ + ]
1392 : 2691 : jpawu_dc=jpawu
1393 [ + - ]: 2691 : if(dmftdc==4) then
1394 : 2717 : jpawu_dc=zero
1395 : : end if
1396 : 2717 : edcdctemp=edcdctemp-half*upawu*n_tot**2
1397 : 2717 : edctemp =edctemp +half*upawu*(n_tot*(n_tot-one))
1398 [ + + - + ]: 2717 : if (nspden/=4.or.pawtab%option_interaction_pawu==2) then
1399 [ + + + + ]: 2082 : if(dmftdc/=5.and.pawtab%usepawu/=4) then
1400 : 1886 : edcdctemp=edcdctemp+half*jpawu_dc*(n_upup**2+n_dndn**2)
1401 : 1886 : 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 : 2787 : edftumdc =edftumdc +edftutemp-edctemp
1453 : 2787 : edftumdcdc=edftumdcdc-edftutemp-edcdctemp
1454 : :
1455 : : !if(pawtab%usepawu/=10.or.pawprtvol>=3) then
1456 [ + + ]: 2787 : 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 [ + + ]: 2787 : if(present(e_ee)) e_ee=e_ee+edftutemp
1499 [ + + ]: 2787 : if(present(e_dc)) e_dc=e_dc+edctemp
1500 [ + + ]: 2787 : if(present(e_dcdc)) e_dcdc=e_dcdc+edcdctemp
1501 : :
1502 : : DBG_EXIT("COLL")
1503 : :
1504 : 2787 : 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 : 5542 : subroutine setnoccmmp(compute_dmat,dimdmat,dmatpawu,dmatudiag,impose_dmat,indsym,my_natom,natom,&
1679 : 5542 : & natpawu,nspinor,nsppol,nsym,ntypat,paw_ij,pawang,pawprtvol,pawrhoij,pawtab,&
1680 : 2771 : & spinat,symafm,typat,useexexch,usepawu, &
1681 [ + + + + : 2771 : & 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 : 2771 : integer,pointer :: my_atmtab(:)
1715 : : real(dp) :: ro(2),sumocc(2)
1716 : 2771 : real(dp),allocatable :: eig(:),hdp(:,:,:),hdp2(:,:),noccmmptemp(:,:,:,:),noccmmp_tmp(:,:,:,:)
1717 : 2771 : real(dp),allocatable :: rwork(:),noccmmp2(:,:,:,:),nocctot2(:)
1718 : 2771 : complex(dp),allocatable :: noccmmp_ylm(:,:,:),noccmmp_jmj(:,:),noccmmp_slm(:,:,:)
1719 : 2771 : 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 : 2771 : type(coeff4_type),allocatable :: tmp_noccmmp(:)
1725 : :
1726 : 2771 : 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 : 2771 : 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 [ + - ]: 2771 : if (my_natom>0) then
1747 [ - + ]: 2771 : if (nsppol/=paw_ij(1)%nsppol) then
1748 : 0 : message='inconsistent values for nsppol!'
1749 : 0 : ABI_BUG(message)
1750 : : end if
1751 [ + + ]: 2771 : if (compute_dmat>0) then
1752 : : if (pawrhoij(1)%nspden/=paw_ij(1)%nspden.and.&
1753 [ + + - + : 2665 : & 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 [ - + ]: 2771 : 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 [ + + - + ]: 2771 : 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 [ + + - + ]: 2771 : 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 [ + + + + : 2771 : 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 [ + + + - ]: 2771 : if (usepawu==0.and.useexexch==0) return
1778 : 2771 : nspden=1;ndij=1;cplex_dij=1
1779 [ + - ]: 2771 : if (my_natom>0) then
1780 : 2771 : nspden=paw_ij(1)%nspden
1781 : 2771 : ndij=paw_ij(1)%ndij
1782 : 2771 : cplex_dij=paw_ij(1)%cplex_dij
1783 : : end if
1784 [ + + + + ]: 2771 : antiferro=(nspden==2.and.nsppol==1)
1785 : 2771 : use_afm=((antiferro).or.((nspden==4).and.afm_noncoll))
1786 : 2771 : dmatudiag_loc=dmatudiag
1787 [ - + - - : 2771 : if (dmatudiag==2.and.(dimdmat==0.or.impose_dmat==0)) dmatudiag_loc=1
- - ]
1788 : :
1789 : : !Set up parallelism over atoms
1790 [ + + + + ]: 2771 : paral_atom=(present(comm_atom).and.(my_natom/=natom))
1791 [ + + ]: 2771 : nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
1792 [ + + ]: 2771 : my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
1793 : 2771 : call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,my_natom_ref=my_natom) !vz_d
1794 [ + + ]: 2771 : wrt_mode='COLL';if (paral_atom) wrt_mode='PERS'
1795 : :
1796 : : !If needed, store dmatpu in suitable format in tmp_noccmmp
1797 [ + + + + ]: 2771 : 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 [ + + + + ]: 2771 : 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 [ + + + + : 2771 : 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 [ + + ]: 8259 : do iatom=atom_min,atom_max
1916 [ + + ]: 5488 : iatom_tot=iatom;if (paral_atom) iatom_tot=my_atmtab(iatom)
1917 : 5488 : itypat=pawrhoij(iatom)%itypat
1918 : 5488 : cplex_rhoij=pawrhoij(iatom)%cplex_rhoij
1919 : :
1920 [ + + ]: 5488 : if (.not. cal_lmom) then
1921 [ + + ]: 5416 : if (useexexch/=0) then
1922 : 4 : lcur=pawtab(itypat)%lexexch
1923 [ + - ]: 5412 : else if (usepawu/=0) then
1924 : 5412 : lcur=pawtab(itypat)%lpawu
1925 : : end if
1926 : : end if
1927 : :
1928 : 72 : if (cal_lmom) lcur=l_orbmom
1929 [ + + ]: 8259 : if (lcur/=-1) then
1930 : :
1931 : : ! ########################################################################################
1932 : : ! # Compute nocc_mmp
1933 : : ! ########################################################################################
1934 [ + + + + : 3697 : if ((usepawu/=0.and.compute_dmat/=0).or.useexexch/=0) then
+ + ]
1935 : :
1936 : :
1937 [ - + + - : 21240 : ABI_MALLOC(l_noccmmp_tmp,(cplex_dij,2*lcur+1,2*lcur+1,ndij))
+ - - + +
- - + + -
- + + - +
- + - - +
- + ]
1938 [ + + + + : 562798 : l_noccmmp_tmp(:,:,:,:)=zero
+ + + + ]
1939 : :
1940 : : ! Loop over spin components
1941 [ - + + - : 17700 : ABI_MALLOC(noccmmptemp,(cplex_dij,2*lcur+1,2*lcur+1,ndij))
- + + - -
+ + - - +
+ - - + -
+ ]
1942 [ + + + + : 562798 : noccmmptemp(:,:,:,:)=zero
+ + + + ]
1943 [ + + ]: 3540 : if(ndij==4) then
1944 [ - + + - : 4365 : ABI_MALLOC(noccmmp2,(cplex_dij,2*lcur+1,2*lcur+1,ndij))
- + + - -
+ + - - +
+ - - + -
+ - + ]
1945 : : end if
1946 [ + + ]: 3540 : 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 : 3540 : nsploop=ndij
1953 [ + + ]: 11854 : do ispden=1,nsploop
1954 : 8314 : jrhoij=1
1955 [ + + ]: 766241 : do irhoij=1,pawrhoij(iatom)%nrhoijsel
1956 : 757927 : klmn=pawrhoij(iatom)%rhoijselect(irhoij)
1957 : 757927 : im1=pawtab(itypat)%klmntomn(1,klmn)
1958 : 757927 : im2=pawtab(itypat)%klmntomn(2,klmn)
1959 : 757927 : in1=pawtab(itypat)%klmntomn(3,klmn)
1960 : 757927 : in2=pawtab(itypat)%klmntomn(4,klmn)
1961 : 757927 : lmin=pawtab(itypat)%indklmn(3,klmn)
1962 : 757927 : lmax=pawtab(itypat)%indklmn(4,klmn)
1963 : :
1964 : 757927 : ro(1:2)=zero
1965 [ + + ]: 2042458 : ro(1:cplex_rhoij)=pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+cplex_rhoij-1,ispden)
1966 [ + + + + ]: 805521 : 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 [ + + + + ]: 757927 : if(lmin==0.and.lmax==2*lcur) then
1971 : 259666 : icount=in1+(in2*(in2-1))/2
1972 [ - + - - ]: 259666 : 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 [ + + ]: 259666 : 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 [ + + ]: 259666 : if(im2>=im1) then
1983 : : l_noccmmp_tmp(1:cplex_dij,im1,im2,ispden)=l_noccmmp_tmp(1:cplex_dij,im1,im2,ispden) &
1984 [ + + ]: 578588 : & +ro(1:cplex_dij)*pawtab(itypat)%phiphjint(icount)
1985 : : end if
1986 : : end if
1987 : 766241 : jrhoij=jrhoij+cplex_rhoij
1988 : : end do ! irhoij
1989 [ + + ]: 49714 : do im2=1,2*lcur+1
1990 [ + + ]: 175783 : do im1=1,im2
1991 : : l_noccmmp_tmp(1,im1,im2,ispden)=l_noccmmp_tmp(1,im1,im2,ispden) &
1992 : 126069 : & +noccmmptemp(1,im2,im1,ispden)
1993 [ + + ]: 126069 : if(cplex_dij==2) l_noccmmp_tmp(2,im1,im2,ispden)=l_noccmmp_tmp(2,im1,im2,ispden) &
1994 : 94148 : & -noccmmptemp(2,im2,im1,ispden)
1995 : : end do
1996 : : end do
1997 [ + + ]: 53254 : do im1=1,2*lcur+1
1998 [ + + ]: 175783 : do im2=1,im1
1999 : 126069 : l_noccmmp_tmp(1,im1,im2,ispden)=l_noccmmp_tmp(1,im2,im1,ispden)
2000 [ + + ]: 167469 : 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 [ - + ]: 3540 : ABI_FREE(noccmmptemp)
2005 : : ! Compute noccmmp2, occupation matrix in the spin basis (upup, dndn, updn, dnup)
2006 [ + + ]: 3540 : 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 [ + + ]: 3540 : if (.not. cal_lmom) then
2042 [ + + ]: 11494 : paw_ij(iatom)%nocctot(:)=zero ! contains nmmp in the n m representation
2043 [ + + + + ]: 6672 : if(ndij==4) nocctot2(:)=zero ! contains nmmp in the upup dndn updn dnup representation
2044 [ + + ]: 11494 : do ispden=1,ndij
2045 [ + + ]: 51742 : do im1=1,2*lcur+1
2046 [ + + ]: 48274 : 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 : 23972 : 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 [ + + ]: 3540 : 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 [ + + ]: 3540 : if (.not. cal_lmom) then
2070 [ + + + + : 3468 : if ((usepawu/=0.and.abs(usepawu)<10).or.(usepawu>=10.and.pawprtvol>=3)) then
+ + + + ]
2071 : 3289 : write(message, '(2a)' ) ch10, &
2072 : 6578 : & '========== DFT+U DATA =================================================== '
2073 : : end if
2074 [ + + ]: 3468 : if (useexexch/=0) then
2075 : 2 : write(message, '(2a)' ) ch10, &
2076 : 4 : & '======= Local ex-exchange (PBE0) DATA =================================== '
2077 : : end if
2078 [ + + + + : 3468 : if (((usepawu/=0.and.abs(usepawu)<10).or.(usepawu>=10.and.pawprtvol>=3)).or.useexexch/=0) then
+ + + + +
+ ]
2079 : 3291 : call wrtout(std_out,message,wrt_mode)
2080 : : end if
2081 [ + + + + ]: 3468 : 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 [ + + + + ]: 3468 : if(abs(usepawu)<10.or.pawprtvol>=3) then ! Always write except if DMFT and pawprtvol low
2088 : 3291 : write(message,'(2a,i5,a,i4,a)') ch10,"====== For Atom", iatom_tot,&
2089 : 6582 : & ", occupations for correlated orbitals. l =",lcur,ch10
2090 : 3291 : call wrtout(std_out,message,wrt_mode)
2091 [ + + ]: 3291 : if(ndij==2) then
2092 [ + + ]: 5934 : do ispden=1,2
2093 : 3956 : write(message,'(a,i4,3a,f10.5)') "Atom", iatom_tot,". Occupations for spin ",&
2094 [ + - ]: 7912 : & trim(dspin(ispden))," =",paw_ij(iatom)%nocctot(ispden)
2095 : 5934 : call wrtout(std_out,message,wrt_mode)
2096 : : end do
2097 : 1978 : write(message,'(a,i4,a,2x,e16.8)') "=> On atom",iatom_tot,", local Mag. is ",&
2098 : 3956 : & paw_ij(iatom)%nocctot(2)-paw_ij(iatom)%nocctot(1)
2099 : 1978 : call wrtout(std_out,message,wrt_mode)
2100 : : end if
2101 [ + + ]: 3291 : 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 : 3291 : write(message,'(2a)') ch10,"== Calculated occupation matrix for correlated orbitals:"
2134 : 3291 : call wrtout(std_out,message,wrt_mode)
2135 [ + + ]: 10963 : do ispden=1,ndij
2136 [ + - ]: 7672 : write(message,'(3a)') ch10,"Calculated occupation matrix for component ",trim(dspinc(ispden+2*(ndij/4)))
2137 : 7672 : call wrtout(std_out,message,wrt_mode)
2138 [ + + ]: 49598 : do im1=1,lcur*2+1
2139 [ + + ]: 38458 : if(cplex_dij==1)&
2140 : : & write(message,'(12(1x,9(1x,f10.5)))')&
2141 : 22182 : & (l_noccmmp_tmp(1,im1,im2,ispden),im2=1,lcur*2+1)
2142 [ + + ]: 38458 : 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 : 46130 : 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 [ + + ]: 3468 : 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 [ + + + + : 546598 : paw_ij(iatom)%noccmmp(:,:,:,:)=zero
+ + + + ]
2177 [ + + + + : 546598 : paw_ij(iatom)%noccmmp(:,:,:,:)=l_noccmmp_tmp(:,:,:,:)
+ + + + ]
2178 [ - + ]: 3468 : 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 [ + + + + ]: 3697 : 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 [ + + + + : 3697 : 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 [ + + + + : 3697 : 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 [ + + + + ]: 3697 : 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 : 3697 : paw_ij(iatom)%has_pawu_occ=2
2543 : :
2544 : : end if ! lcur
2545 : : end do ! iatom
2546 : :
2547 : : !Memory deallocation
2548 [ + + + + ]: 2771 : 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 : 2771 : call free_my_atmtab(my_atmtab,my_atmtab_allocated)
2560 : :
2561 : : DBG_EXIT("COLL")
2562 : :
2563 [ + + - + : 5542 : 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), contiguous, pointer :: paw_ij_all(:)
3157 : 9 : type(pawrhoij_type),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 : : !!***
|
