Line data Source code
1 : !!****m* m_paw_occupancies/m_paw_occupancies
2 : !! NAME
3 : !! m_paw_occupancies
4 : !!
5 : !! FUNCTION
6 : !! This module contains routines related to the computation of PAW on-site occupancies (rhoij).
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 2018-2026 ABINIT group (FJ, MT)
10 : !! This file is distributed under the terms of the
11 : !! GNU General Public License, see ~abinit/COPYING
12 : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : !!
14 : !! SOURCE
15 :
16 : #if defined HAVE_CONFIG_H
17 : #include "config.h"
18 : #endif
19 :
20 : #include "abi_common.h"
21 :
22 : MODULE m_paw_occupancies
23 :
24 : use defs_basis
25 : use m_abicore
26 : use m_errors
27 : use m_xmpi
28 :
29 : use defs_abitypes, only : MPI_type
30 : use m_pawtab, only : pawtab_type
31 : use m_pawrhoij, only : pawrhoij_type,pawrhoij_init_unpacked,pawrhoij_mpisum_unpacked, &
32 : & pawrhoij_alloc,pawrhoij_free,pawrhoij_inquire_dim
33 : use m_pawcprj, only : pawcprj_type,pawcprj_alloc,pawcprj_get, &
34 : & pawcprj_gather_spin, pawcprj_free, pawcprj_mpi_send, &
35 : & pawcprj_mpi_recv, pawcprj_copy, pawcprj_unpack, pawcprj_pack
36 : use m_paw_io, only : pawio_print_ij
37 : use m_paral_atom, only : get_my_atmtab,free_my_atmtab
38 : use m_paw_dmft, only : paw_dmft_type
39 : use m_mpinfo, only : proc_distrb_cycle
40 :
41 : implicit none
42 :
43 : private
44 :
45 : !public procedures.
46 : public :: pawmkrhoij ! Compute the PAW occupancies rhoij
47 : public :: pawaccrhoij ! Accumulate the contribution of one band to the PAW occupancies rhoij
48 : public :: initrhoij ! Initialize the PAW occupancies rhoij from atomic data
49 :
50 : CONTAINS !========================================================================================
51 : !!***
52 :
53 : !----------------------------------------------------------------------
54 :
55 : !!****f* m_paw_occupancies/pawmkrhoij
56 : !!
57 : !! NAME
58 : !! pawmkrhoij
59 : !!
60 : !! FUNCTION
61 : !! Calculate the PAW quantities rhoij (augmentation occupancies)
62 : !! Remember:for each atom, rho_ij=Sum_{n,k} {occ(n,k)*<Cnk|p_i><p_j|Cnk>}
63 : !!
64 : !! INPUTS
65 : !! atindx1(natom)=index table for atoms, inverse of atindx
66 : !! cprj(natom,mcprj)= wave functions projected with non-local projectors:
67 : !! cprj_nk(i)=<p_i|Cnk> where p_i is a non-local projector.
68 : !! dimcprj(natom)=array of dimensions of array cprj (ordered by atom-type)
69 : !! istwfk(nkpt)=parameter that describes the storage of wfs
70 : !! kptopt=option for the generation of k points
71 : !! mband=maximum number of bands
72 : !! mband_cprj=maximum number of bands used in the dimensioning of cprj array (usually mband/nproc_band)
73 : !! mcprj=size of projected wave-functions array (cprj) =nspinor*mband*mkmem*nsppol
74 : !! mkmem =number of k points treated by this node.
75 : !! mpi_enreg=information about MPI parallelization
76 : !! natom=number of atoms in cell
77 : !! nband=number of bands for all k points
78 : !! nkpt=number of k points
79 : !! nspden=number of spin components for the density
80 : !! nspinor=number of spinorial components of the wavefunctions
81 : !! nsppol=1 for unpolarized, 2 for spin-polarized
82 : !! occ(mband*nkpt*nsppol)=occupation number for each band for each k
83 : !! paral_kgb=Flag related to the kpoint-band-fft parallelism
84 : !! paw_dmft <type(paw_dmft_type)>= paw+dmft related data
85 : !! unpaw=unit number for cprj PAW data (if used)
86 : !! wtk(nkpt)=weight assigned to each k point
87 : !!
88 : !! SIDE EFFECTS
89 : !! pawrhoij(natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data
90 : !! On input: arrays dimensions
91 : !! On output:
92 : !! pawrhoij(:)%rhoij_(lmn2_size,nspden)=
93 : !! Sum_{n,k} {occ(n,k)*conjugate[cprj_nk(ii)].cprj_nk(jj)} (non symetrized)
94 : !!
95 : !! NOTES
96 : !! The cprj are distributed over band processors.
97 : !! Only the mod((iband-1)/mpi_enreg%bandpp,mpi_enreg%nproc_band) projectors
98 : !! are stored on each proc.
99 : !!
100 : !! SOURCE
101 :
102 33362 : subroutine pawmkrhoij(atindx,atindx1,cprj,dimcprj,istwfk,kptopt,mband,mband_cprj,mcprj,mkmem,mpi_enreg,&
103 16681 : & natom,nband,nkpt,nspden,nspinor,nsppol,occ,paral_kgb,paw_dmft,pawrhoij,unpaw,usewvl,wtk)
104 : !Arguments ---------------------------------------------
105 : !scalars
106 : integer,intent(in) :: kptopt,mband,mband_cprj,mcprj,mkmem,natom,nkpt,nspden,nspinor,nsppol
107 : integer,intent(in) :: paral_kgb,unpaw,usewvl
108 : type(MPI_type),intent(in) :: mpi_enreg
109 : !arrays
110 : integer,intent(in) :: atindx(natom),atindx1(natom),dimcprj(natom),istwfk(nkpt)
111 : integer,intent(in) :: nband(nkpt*nsppol)
112 : real(dp),intent(in) :: occ(mband*nkpt*nsppol),wtk(nkpt)
113 : type(pawcprj_type),target,intent(in) :: cprj(natom,mcprj)
114 : type(paw_dmft_type),intent(in) :: paw_dmft
115 : type(pawrhoij_type),intent(inout),target:: pawrhoij(:)
116 :
117 : !Local variables ---------------------------------------
118 : !scalars
119 : integer,parameter :: max_nband_cprj=100
120 : integer :: bdtot_index,cplex
121 : integer :: cplex_rhoij,iatom,iatom_tot,ib,ib1,iband,ibc1,ibg,ib_this_proc,ierr
122 : integer :: ikpt,iorder_cprj,isppol,jb_this_proc,jbg,me,my_nspinor,nband_k,nband_k_cprj
123 : integer :: nbandc1,nband_k_cprj_read,nband_k_cprj_used,nprocband,nrhoij,nspden_rhoij
124 : integer :: option,spaceComm,use_nondiag_occup_dmft
125 : logical :: locc_test,paral_atom,use_timerev,use_zeromag
126 : integer :: ib1_this_proc, ib_loop, proc_sender, proc_recver
127 : real(dp) :: wtk_k
128 : character(len=500) :: msg
129 :
130 : !arrays
131 : integer :: n2buff
132 16681 : integer, allocatable :: typat(:),req_correl(:,:,:)
133 : real(dp) :: occup(2)
134 16681 : real(dp) ABI_ASYNC, allocatable :: buffer_cprj_correl(:,:,:)
135 : character(len=8),parameter :: dspin(6)=(/"up ","down ","dens (n)","magn (x)","magn (y)","magn (z)"/)
136 16681 : type(pawcprj_type),allocatable :: cprj_tmp(:,:),cwaveprj(:,:),cwaveprjb(:,:)
137 16681 : type(pawcprj_type),pointer :: cprj_ptr(:,:)
138 16681 : type(pawrhoij_type),pointer :: pawrhoij_all(:)
139 :
140 : !************************************************************************
141 :
142 : DBG_ENTER("COLL")
143 :
144 0 : ABI_CHECK(mkmem/=0,"mkmem==0 not supported anymore!")
145 :
146 : !Init MPI data
147 : ! spaceComm=mpi_enreg%comm_cell
148 : ! if(mpi_enreg%paral_kgb==1) spaceComm=mpi_enreg%comm_kpt
149 16681 : spaceComm=mpi_enreg%comm_kpt
150 16681 : me=mpi_enreg%me_kpt
151 :
152 : !Check size of cprj
153 16681 : my_nspinor=max(1,nspinor/mpi_enreg%nproc_spinor)
154 16681 : if (mcprj/=my_nspinor*mband_cprj*mkmem*nsppol) then
155 0 : msg=' wrong size for cprj !'
156 0 : ABI_BUG(msg)
157 : end if
158 :
159 : !Check if cprj is distributed over bands
160 16681 : nprocband=(mband/mband_cprj)
161 16681 : if (paral_kgb==1.and.nprocband/=mpi_enreg%nproc_band) then
162 0 : msg='mband/mband_cprj must be equal to nproc_band!'
163 0 : ABI_BUG(msg)
164 : end if
165 :
166 16681 : if( usewvl==1 .and. (nprocband/=1)) then
167 0 : write(msg,'(2a)') ch10,&
168 0 : & 'Parallelization over bands is not compatible with WAVELETS!'
169 0 : ABI_ERROR(msg)
170 : end if
171 :
172 : !Initialise and check dmft variables
173 16681 : if(paw_dmft%use_sc_dmft/=0) then
174 113 : nbandc1=paw_dmft%mbandc
175 : else
176 : nbandc1=1
177 : end if
178 :
179 : !Size of pawrhoij datastructure
180 16681 : nrhoij=size(pawrhoij)
181 :
182 : !Check if pawrhoij is distributed over atomic sites
183 16681 : paral_atom=(nrhoij/=natom.and.mpi_enreg%nproc_atom>1)
184 0 : if (paral_atom.and.nrhoij/=mpi_enreg%my_natom) then
185 0 : msg='Size of pawrhoij should be natom or my_natom!'
186 0 : ABI_BUG(msg)
187 : end if
188 :
189 : !Allocate temporary cwaveprj storage
190 135627 : ABI_MALLOC(cwaveprj,(natom,nspinor))
191 16681 : call pawcprj_alloc(cwaveprj,0,dimcprj)
192 16681 : if(paw_dmft%use_sc_dmft/=0) then
193 855 : ABI_MALLOC(cwaveprjb,(natom,nspinor))
194 113 : call pawcprj_alloc(cwaveprjb,0,dimcprj)
195 : end if
196 :
197 16681 : if (paw_dmft%use_sc_dmft /= 0 .and. mpi_enreg%paral_kgb /= 0) then
198 40 : if(paw_dmft%use_bandc(mpi_enreg%me_band+1)) then
199 120 : n2buff = nspinor*sum(dimcprj)
200 80 : ABI_MALLOC(buffer_cprj_correl,(2,n2buff,nbandc1))
201 100 : ABI_MALLOC(req_correl,(nbandc1, nkpt, nsppol))
202 1340 : req_correl(:,:,:) = 0
203 : end if
204 : end if
205 :
206 : !Initialize temporary file (if used)
207 16681 : iorder_cprj=0
208 :
209 : !Build and initialize unpacked rhoij (to be computed here)
210 16681 : call pawrhoij_init_unpacked(pawrhoij)
211 :
212 : !If pawrhoij is MPI-distributed over atomic sites, we have to gather it
213 16681 : if (paral_atom) then
214 0 : ABI_MALLOC(pawrhoij_all,(natom))
215 0 : ABI_MALLOC(typat,(natom))
216 0 : typat(:)=0;cplex_rhoij=0;nspden_rhoij=0
217 0 : do iatom=1,nrhoij
218 0 : iatom_tot=mpi_enreg%my_atmtab(iatom)
219 0 : typat(iatom_tot)=pawrhoij(iatom)%itypat
220 0 : cplex_rhoij=max(cplex_rhoij,pawrhoij(iatom)%cplex_rhoij)
221 0 : nspden_rhoij=max(nspden_rhoij,pawrhoij(iatom)%nspden)
222 : end do
223 0 : call xmpi_sum(typat,mpi_enreg%comm_atom,ierr)
224 0 : call xmpi_max_ip(cplex_rhoij,mpi_enreg%comm_atom,ierr)
225 0 : call xmpi_max_ip(nspden_rhoij,mpi_enreg%comm_atom,ierr)
226 : call pawrhoij_alloc(pawrhoij_all,cplex_rhoij,nspden_rhoij,nspinor,&
227 0 : & nsppol,typat,lmnsize=dimcprj,use_rhoijp=0,use_rhoij_=1)
228 0 : ABI_FREE(typat)
229 : else
230 16681 : pawrhoij_all => pawrhoij
231 : end if
232 :
233 : !LOOP OVER SPINS
234 16681 : option=1
235 16681 : use_timerev=(kptopt>0.and.kptopt<3)
236 16681 : use_zeromag=(pawrhoij_all(1)%nspden==4.and.nspden==1)
237 16681 : bdtot_index=0;ibg=0;jbg=0
238 35278 : do isppol=1,nsppol
239 :
240 : ! LOOP OVER k POINTS
241 209515 : do ikpt=1,nkpt
242 :
243 174237 : nband_k=nband(ikpt+(isppol-1)*nkpt)
244 174237 : nband_k_cprj=nband_k/nprocband
245 174237 : wtk_k=wtk(ikpt)
246 :
247 174237 : if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me)) then
248 39377 : bdtot_index=bdtot_index+nband_k
249 39377 : cycle
250 : end if
251 :
252 134860 : cplex=2;if (istwfk(ikpt)>1) cplex=1
253 :
254 : ! In case of spinors parallelism, need some extra storage
255 134860 : if (mpi_enreg%paral_spinor==1) then
256 464 : nband_k_cprj_used=min(max_nband_cprj,nband_k_cprj)
257 25760 : ABI_MALLOC(cprj_tmp,(natom,my_nspinor*nband_k_cprj_used))
258 49664 : ABI_MALLOC(cprj_ptr,(natom, nspinor*nband_k_cprj_used))
259 464 : call pawcprj_alloc(cprj_tmp,0,dimcprj)
260 464 : call pawcprj_alloc(cprj_ptr,0,dimcprj)
261 : else
262 134396 : cprj_ptr => cprj
263 : end if
264 :
265 : ! In case of band parallelism combined with self consistent DMFT, need to
266 : ! exchange bands cprj
267 134860 : if (paw_dmft%use_sc_dmft /= 0 .and. mpi_enreg%paral_kgb /= 0) then
268 200 : if (paw_dmft%use_bandc(mpi_enreg%me_band+1)) then
269 : ! only proc using correlated band have to do this
270 1300 : do ibc1=1,nbandc1
271 1200 : proc_sender = paw_dmft%bandc_proc(ibc1)
272 1300 : if(proc_sender == mpi_enreg%me_band) then
273 :
274 : ! get the index of band local to this proc
275 600 : ib1 = paw_dmft%include_bands(ibc1)
276 600 : ib1_this_proc = 0
277 10500 : do ib_loop=1,ib1-1
278 10500 : if (mod((ib_loop-1)/mpi_enreg%bandpp,mpi_enreg%nproc_band) == mpi_enreg%me_band) then
279 2300 : ib1_this_proc = ib1_this_proc+1
280 : end if
281 : end do
282 600 : ib1_this_proc = ib1_this_proc+1
283 :
284 : ! extract the band
285 : call pawcprj_get(atindx1,cwaveprjb,cprj_ptr,natom,ib1_this_proc,ibg,ikpt,&
286 : & iorder_cprj,isppol,mband_cprj,mkmem,natom,1,nband_k_cprj,nspinor,nsppol,&
287 600 : & unpaw,mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
288 :
289 600 : call pawcprj_pack(dimcprj,cwaveprjb,buffer_cprj_correl(:,:,ibc1))
290 3000 : do proc_recver=0,mpi_enreg%nproc_band-1
291 3000 : if (proc_sender /= proc_recver .and. paw_dmft%use_bandc(proc_recver+1)) then
292 : ! locc_test = At least one of the bands used by proc_recver have a non neglectable occnd
293 : locc_test = .false.
294 7800 : do ib_loop=1,nbandc1
295 7800 : if(proc_recver == paw_dmft%bandc_proc(ib_loop)) then
296 3500 : ib = paw_dmft%include_bands(ib_loop)
297 : locc_test = locc_test .or. (abs(paw_dmft%occnd(1,ib,ib1,ikpt,isppol))+&
298 3500 : & abs(paw_dmft%occnd(2,ib,ib1,ikpt,isppol))>tol8)
299 : end if
300 : end do
301 600 : if(locc_test) then
302 : ! send to proc_recver
303 : ierr = 0
304 : call xmpi_isend(buffer_cprj_correl(:,:,ibc1),proc_recver,&
305 : & 10000+ibc1+nbandc1*(ikpt+nsppol*isppol),mpi_enreg%comm_band,&
306 259 : & req_correl(ibc1,ikpt,isppol),ierr)
307 : ! force sending or buffering
308 259 : call xmpi_wait(req_correl(ibc1,ikpt,isppol), ierr)
309 : end if
310 : end if
311 : end do
312 : else
313 : ! locc_test = At least one of the bands used by this proc have a non neglectable occnd
314 : locc_test = .false.
315 7800 : do ib_loop=1,nbandc1
316 7800 : if(mpi_enreg%me_band == paw_dmft%bandc_proc(ib_loop)) then
317 3500 : ib = paw_dmft%include_bands(ib_loop)
318 3500 : ib1 = paw_dmft%include_bands(ibc1)
319 : locc_test = locc_test .or. (abs(paw_dmft%occnd(1,ib,ib1,ikpt,isppol))+&
320 3500 : & abs(paw_dmft%occnd(2,ib,ib1,ikpt,isppol))>tol8)
321 : end if
322 : end do
323 600 : if(locc_test) then
324 : ! recv from proc_sender
325 : ierr = 0
326 : call xmpi_irecv(buffer_cprj_correl(:,:,ibc1),proc_sender,&
327 : & 10000+ibc1+nbandc1*(ikpt+nsppol*isppol),mpi_enreg%comm_band,&
328 259 : & req_correl(ibc1,ikpt,isppol),ierr)
329 : end if
330 : end if
331 : end do
332 : end if
333 : end if
334 :
335 134860 : ierr = 0
336 : ! LOOP OVER BANDS
337 134860 : ib_this_proc=0;jb_this_proc=0
338 1684663 : do ib=1,nband_k
339 1549803 : iband=bdtot_index+ib
340 :
341 : ! Parallelization: treat only some bands
342 : if(xmpi_paral==1)then
343 1549803 : if (paral_kgb==1) then
344 402484 : if (mod((ib-1)/mpi_enreg%bandpp,mpi_enreg%nproc_band)/=mpi_enreg%me_band) cycle
345 : else
346 1147319 : if (mpi_enreg%proc_distrb(ikpt,ib,isppol)/=me) cycle
347 : end if
348 : end if
349 1309051 : ib_this_proc=ib_this_proc+1
350 :
351 : ! In case of spinors parallelism, gather cprj because we need both components together
352 : ! We do that nband_k_cprj_used by nband_k_cprj_used bands
353 1309051 : if (mpi_enreg%paral_spinor==1) then
354 7968 : jb_this_proc=jb_this_proc+1
355 7968 : if (mod(jb_this_proc,nband_k_cprj_used)==1) then
356 464 : ib_this_proc=1
357 464 : nband_k_cprj_read=nband_k_cprj_used
358 464 : if (nband_k_cprj<jb_this_proc+nband_k_cprj_used-1) nband_k_cprj_read=nband_k_cprj-jb_this_proc+1
359 : call pawcprj_get(atindx1,cprj_tmp,cprj,natom,jb_this_proc,jbg,ikpt,iorder_cprj,isppol,&
360 : & mband_cprj,mkmem,natom,nband_k_cprj_read,nband_k_cprj,my_nspinor,nsppol,unpaw,&
361 464 : & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
362 : call pawcprj_gather_spin(cprj_tmp,cprj_ptr,natom,nband_k_cprj_read,my_nspinor,nspinor,&
363 464 : & mpi_enreg%comm_spinor,ierr)
364 : end if
365 : end if
366 :
367 : ! DMFT: LOOP ON ADDITIONAL BANDS
368 2934830 : do ibc1=1,nbandc1
369 : ! check if dmft and occupations
370 : ! write(std_out,*) 'ib,ibc1 ',ib,ibc1
371 :
372 : ! if ib is not part a band correlated in dmft do not repeat the following
373 1490919 : if(ibc1 /= 1) then
374 181868 : if (paw_dmft%use_sc_dmft == 0) then
375 : cycle
376 : else
377 181868 : if (.not.(paw_dmft%band_in(ib))) cycle
378 : end if
379 : end if
380 :
381 : ! DMFT stuff: extract cprj and occupations for additional band
382 1452163 : if(paw_dmft%use_sc_dmft /= 0) then
383 157268 : if(paw_dmft%band_in(ib)) then
384 : ! write(std_out,*) 'use_sc_dmft=1 ib,ib1',ib,ib1
385 : ! write(std_out,*) 'ib, ib1 ',paw_dmft%band_in(ib),paw_dmft%band_in(ib1)
386 :
387 150648 : ib1 = paw_dmft%include_bands(ibc1) ! indice reel de la bande
388 :
389 150648 : use_nondiag_occup_dmft = 1
390 150648 : locc_test = abs(paw_dmft%occnd(1,ib,ib1,ikpt,isppol))+abs(paw_dmft%occnd(2,ib,ib1,ikpt,isppol))>tol8
391 :
392 150648 : occup(1) = paw_dmft%occnd(1,ib,ib1,ikpt,isppol)
393 150648 : occup(2) = paw_dmft%occnd(2,ib,ib1,ikpt,isppol)
394 :
395 : ! write(std_out,*) 'use_sc_dmft=1,band_in(ib)=1, ib,ibc1',ib,ib1,locc_test
396 : !
397 150648 : if (locc_test .or. mkmem == 0) then
398 :
399 28779 : if (paral_kgb==1) then
400 : ! cprj have already been extracted
401 1508 : if (paw_dmft%bandc_proc(ibc1) /= mpi_enreg%me_band) then
402 : ! if the band is not on this proc, wait for the recv to complete
403 : ierr = 0
404 616 : call xmpi_wait(req_correl(ibc1,ikpt,isppol), ierr)
405 : end if
406 1508 : call pawcprj_unpack(dimcprj,cwaveprjb,buffer_cprj_correl(:,:,ibc1))
407 : else ! paral_kgb /= 0
408 : call pawcprj_get(atindx1,cwaveprjb,cprj_ptr,natom,ib1,ibg,ikpt,iorder_cprj,isppol,&
409 : & mband_cprj,mkmem,natom,1,nband_k_cprj,nspinor,nsppol,unpaw,&
410 27271 : & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
411 : end if
412 : end if
413 : else ! nbandc1=1
414 6620 : use_nondiag_occup_dmft=0
415 6620 : locc_test = (abs(occ(iband))>tol8)
416 6620 : occup(1) = occ(iband)
417 : end if
418 : else ! nbandc1=1
419 1294895 : use_nondiag_occup_dmft=0
420 1294895 : locc_test = (abs(occ(iband))>tol8)
421 1294895 : occup(1) = occ(iband)
422 : end if
423 :
424 : ! Extract cprj for current band
425 : ! Must read cprj when mkmem=0 (even if unused) to have right pointer inside _PAW file
426 1452163 : if (locc_test.or.mkmem==0) then
427 : call pawcprj_get(atindx1,cwaveprj,cprj_ptr,natom,ib_this_proc,ibg,ikpt,iorder_cprj,isppol,&
428 : & mband_cprj,mkmem,natom,1,nband_k_cprj,nspinor,nsppol,unpaw,&
429 866539 : & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
430 : end if
431 :
432 : ! Accumulate contribution from (occupied) current band
433 3001966 : if (locc_test) then
434 866539 : if(use_nondiag_occup_dmft == 1) then
435 : call pawaccrhoij(atindx,cplex,cwaveprj,cwaveprjb,0,isppol,nrhoij,natom,&
436 : & nspinor,occup(1),option,pawrhoij_all,use_timerev,use_zeromag,wtk_k,&
437 28779 : & occ_k_2=occup(2))
438 : else
439 : call pawaccrhoij(atindx,cplex,cwaveprj,cwaveprj ,0,isppol,nrhoij,natom,&
440 837760 : & nspinor,occup(1),option,pawrhoij_all,use_timerev,use_zeromag,wtk_k)
441 : end if
442 : end if
443 : end do ! ib1c
444 : end do ! ib
445 :
446 134860 : if (mpi_enreg%paral_spinor==1) then
447 464 : call pawcprj_free(cprj_tmp)
448 464 : call pawcprj_free(cprj_ptr)
449 16400 : ABI_FREE(cprj_tmp)
450 464 : ABI_FREE(cprj_ptr)
451 : else
452 134396 : nullify(cprj_ptr)
453 : end if
454 :
455 134860 : bdtot_index=bdtot_index+nband_k
456 153457 : if (mkmem/=0) then
457 134860 : if (mpi_enreg%paral_spinor==0) then
458 134396 : ibg=ibg+ nspinor*nband_k_cprj
459 : else
460 464 : jbg=jbg+my_nspinor*nband_k_cprj
461 : end if
462 : end if
463 :
464 : end do ! ikpt
465 : end do ! isppol
466 :
467 : !call xmpi_barrier(mpi_enreg%comm_band)
468 :
469 : !deallocate temporary cwaveprj/cprj storage
470 16681 : call pawcprj_free(cwaveprj)
471 66530 : ABI_FREE(cwaveprj)
472 :
473 16681 : if(paw_dmft%use_sc_dmft/=0) then
474 113 : call pawcprj_free(cwaveprjb)
475 515 : ABI_FREE(cwaveprjb)
476 : end if
477 :
478 16681 : if (allocated(buffer_cprj_correl)) then
479 20 : ABI_FREE(buffer_cprj_correl)
480 20 : ABI_FREE(req_correl)
481 : end if
482 :
483 : !MPI: need to exchange rhoij_ between procs
484 16681 : if (paral_kgb==1.and.nprocband>1) then
485 5132 : call pawrhoij_mpisum_unpacked(pawrhoij_all,spaceComm,comm2=mpi_enreg%comm_band)
486 : else
487 11549 : call pawrhoij_mpisum_unpacked(pawrhoij_all,spaceComm)
488 : end if
489 :
490 : !In case of distribution over atomic sites, dispatch rhoij
491 16681 : if (paral_atom) then
492 0 : do iatom=1,nrhoij
493 0 : iatom_tot=mpi_enreg%my_atmtab(iatom)
494 0 : pawrhoij(iatom)%rhoij_(:,:)=pawrhoij_all(iatom_tot)%rhoij_(:,:)
495 : end do
496 0 : call pawrhoij_free(pawrhoij_all)
497 0 : ABI_FREE(pawrhoij_all)
498 : end if
499 :
500 : DBG_EXIT("COLL")
501 :
502 33362 : end subroutine pawmkrhoij
503 : !!***
504 :
505 : !----------------------------------------------------------------------
506 :
507 : !!****f* m_paw_occupancies/pawaccrhoij
508 : !!
509 : !! NAME
510 : !! pawaccrhoij
511 : !!
512 : !! FUNCTION
513 : !! Accumulate the PAW quantities rhoij (augmentation occupancies)
514 : !! or their 1-st order change or their gradient vs r
515 : !! Add the contribution of a given k-point and band
516 : !! Remember: for each atom, rho_ij=Sum_{n,k} {occ(n,k)*<Cnk|p_i><p_j|Cnk>}
517 : !!
518 : !! INPUTS
519 : !! atindx(natom)=index table for atoms (sorted-->random), inverse of atindx.
520 : !! cplex: if 1, WFs (or 1st-order WFs) are REAL, if 2, COMPLEX
521 : !! cwaveprj(natom,nspinor)= LEFT wave function at given n,k
522 : !! projected with non-local projectors: cwaveprj=<p_i|Cnk>
523 : !! cwaveprj1(natom,nspinor)= RIGHT wave function at n,k,q
524 : !! projected with non-local projectors: cwaveprj1=<p_i|C1nk,q>
525 : !! * USED for RF : C1nk is the first-order wave function
526 : !! * USED for DMFT: C1nk is the RIGHT wave function
527 : !! * NOT USED in usual GS case; can be set to cwaveprj in that case
528 : !! ipert=index of perturbation (RF only, i.e. option=2)
529 : !! isppol=index of current spin component
530 : !! mpi_atmtab(:)=--optional-- indexes of the atoms treated by current proc
531 : !! comm_atom=--optional-- MPI communicator over atoms
532 : !! my_natom=number of atoms treated by current processor
533 : !! natom=number of atoms in cell
534 : !! nspinor=number of spinorial components (on current proc)
535 : !! occ_k=occupation number for current band n,k
536 : !! option: choice of calculation:
537 : !! 1: update rhoij (Ground-State)
538 : !! 2: update 1-st order rhoij (Response Function) according to ipert
539 : !! 3: update gradients of rhoij with respect to r,strain of both
540 : !! use_timerev=.TRUE. if time-reversal symmetry is used (WF(-k)=Conjg[WF(k)])
541 : !! use_zeromag=.TRUE. if rhoij "magnetization" is enforced to be zero
542 : !! Applies only when nspden_rhoij=4 (note: only the real part is set to zero)
543 : !! wtk_k=weight assigned to current k-point
544 : !! [occ_k_2]=??
545 : !!
546 : !! SIDE EFFECTS
547 : !! pawrhoij(natom) <type(pawrhoij_type)>= GS: paw rhoij occupancies and related data
548 : !! RF: 1-st order paw rhoij occupancies and related data
549 : !! On output, has been updated with the contribution of current n,k
550 : !! === option=1:
551 : !! pawrhoij(:)%rhoij_(lmn2_size,nspden)= (non symetrized)
552 : !! Sum_{n,k} {occ(n,k)*conjugate[cprj_nk(ii)].cprj_nk(jj)}
553 : !! === option=2:
554 : !! pawrhoij(:)%rhoij_(lmn2_size,nspden)= (non symetrized)
555 : !! Sum_{n,k} {occ(n,k)*(conjugate[cprj_nk(ii)].cprj1_nk,q(jj)
556 : !! conjugate[cprj_nk(jj)].cprj1_nk,q(ii)}
557 : !! + Sum_{n,k} {occ(n,k)*(conjugate[dcprj_nk(ii)/dlambda].cprj_nk(jj)
558 : !! +conjugate[cprj_nk(ii)].dcprj_nk(jj)/dlambda)}
559 : !! === option=3:
560 : !! pawrhoij(:)%grhoij(lmn2_size,mu,nspden)= (non symetrized)
561 : !! Sum_{n,k} {occ(n,k)*(conjugate[dcprj_nk(ii)/dr_mu].cprj_nk(jj)
562 : !! +conjugate[cprj_nk(ii)].dcprj_nk(jj)/dr_mu)}
563 : !!
564 : !! SOURCE
565 :
566 2997664 : subroutine pawaccrhoij(atindx,cplex,cwaveprj,cwaveprj1,ipert,isppol,my_natom,natom,&
567 1498832 : & nspinor,occ_k,option,pawrhoij,use_timerev,use_zeromag,wtk_k,occ_k_2, &
568 1498832 : & comm_atom,mpi_atmtab ) ! optional (parallelism)
569 :
570 : !Arguments ---------------------------------------------
571 : !scalars
572 : integer,intent(in) :: cplex,ipert,isppol,my_natom,natom,nspinor,option
573 : integer,optional,intent(in) :: comm_atom
574 : logical,intent(in) :: use_timerev,use_zeromag
575 : real(dp),intent(in) :: occ_k,wtk_k
576 : real(dp),optional,intent(in) :: occ_k_2
577 : !arrays
578 : integer,intent(in) :: atindx(natom)
579 : integer,optional,target,intent(in) :: mpi_atmtab(:)
580 : type(pawcprj_type),intent(in) :: cwaveprj(natom,nspinor),cwaveprj1(natom,nspinor)
581 : type(pawrhoij_type),intent(inout) :: pawrhoij(my_natom)
582 :
583 : !Local variables ---------------------------------------
584 : !scalars
585 : integer :: cplex_rhoij,iatm,iatom,iatom1,ilmn,iplex,iq0,j0lmn,jlmn,klmn,klmn_im,klmn_re
586 : integer :: mu,my_comm_atom,ncpgr,nspden_rhoij
587 : logical :: compute_impart,compute_impart_cplex,substract_diagonal
588 : logical :: my_atmtab_allocated,paral_atom
589 : real(dp) :: ro11_im,ro11_re,ro12_im,ro12_re,ro21_im,ro21_re,ro22_im,ro22_re,weight,weight_2
590 : character(len=500) :: message
591 : !arrays
592 1498832 : integer,pointer :: my_atmtab(:)
593 2997664 : real(dp) :: cpi0(2,nspinor),cpi1(2,nspinor),cpj0(2,nspinor),cpj1(2,nspinor)
594 1498832 : real(dp) :: dcpi0(2,nspinor,9),dcpj0(2,nspinor,9)
595 :
596 : ! ***********************************************************************
597 :
598 : DBG_ENTER("COLL")
599 :
600 1498832 : if (my_natom==0) return
601 :
602 1498832 : ncpgr=0
603 1498832 : if (option==2.and.(ipert<=natom.or.ipert==natom+3.or.ipert==natom+4)) ncpgr=1
604 1498832 : if (option==3) ncpgr=cwaveprj(1,1)%ncpgr
605 :
606 : !Tests
607 1498832 : if(option==2.and.(ipert==natom+1.or.ipert==natom+10.or.ipert==natom+11)) then
608 0 : message = 'Not relevant for ipert=natom+1 or ipert=natom+10 or ipert=natom+11!'
609 0 : ABI_BUG(message)
610 : end if
611 1498832 : if(option==2.and.cwaveprj(1,1)%ncpgr<ncpgr) then
612 0 : message = 'Error on cwaveprj1 factors derivatives!'
613 0 : ABI_BUG(message)
614 : end if
615 1498832 : if(option==3.and.cwaveprj(1,1)%ncpgr/=ncpgr) then
616 0 : message = 'Error on cwaveprj factors derivatives!'
617 0 : ABI_BUG(message)
618 : end if
619 1498832 : if (pawrhoij(1)%qphase==2.and.option/=2) then
620 0 : message = 'pawaccrhoij: qphase=2 only allowed with option=2 (1st-order rhoij)!'
621 0 : ABI_BUG(message)
622 : end if
623 :
624 : !Set up parallelism over atoms
625 1498832 : paral_atom=(present(comm_atom).and.(my_natom/=natom))
626 1498832 : nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
627 1498832 : my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
628 : call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,&
629 1498832 : & my_natom_ref=my_natom)
630 :
631 1498832 : weight=wtk_k*occ_k
632 1498832 : weight_2=zero
633 1498832 : if(present(occ_k_2)) then
634 28779 : if (cplex==1) then
635 : !DMFT need complex cprj
636 0 : message='DMFT computation must be done with complex cprj. Check that istwfk = 1!'
637 0 : ABI_ERROR(message)
638 : end if
639 28779 : weight_2=wtk_k*occ_k_2
640 : end if
641 1498832 : if (pawrhoij(1)%nspden==2.and.pawrhoij(1)%nsppol==1.and.nspinor==1) weight=half*weight
642 1498832 : if (pawrhoij(1)%nspden==2.and.pawrhoij(1)%nsppol==1.and.nspinor==1.and.present(occ_k_2)) weight_2=half*weight_2
643 :
644 1498832 : if (option==1) then
645 :
646 : ! ==================================================================
647 : ! === OPTION 1: Accumulate (n,k) contribution to rhoij =============
648 : ! ==================================================================
649 867147 : compute_impart=((.not.use_timerev).and.(pawrhoij(1)%cplex_rhoij==2))
650 174740 : compute_impart_cplex=((compute_impart).and.(cplex==2))
651 867147 : if (nspinor==1) then
652 711583 : cplex_rhoij=pawrhoij(1)%cplex_rhoij
653 711583 : if (cplex_rhoij==1) then
654 2596961 : do iatom=1,my_natom
655 1935460 : iatom1=iatom;if (paral_atom) iatom1=my_atmtab(iatom)
656 1935460 : iatm=atindx(iatom1)
657 25742102 : do jlmn=1,pawrhoij(iatom)%lmn_size
658 23145141 : j0lmn=jlmn*(jlmn-1)/2
659 65875531 : cpj0(1:cplex,1)=cwaveprj(iatm,1)%cp(1:cplex,jlmn)
660 204824494 : do ilmn=1,jlmn
661 179743893 : klmn=j0lmn+ilmn
662 521513241 : cpi0(1:cplex,1)=cwaveprj1(iatm,1)%cp(1:cplex,ilmn)
663 : ro11_re=zero
664 521513241 : do iplex=1,cplex
665 521513241 : ro11_re=ro11_re+cpi0(iplex,1)*cpj0(iplex,1)
666 : end do
667 179743893 : pawrhoij(iatom)%rhoij_(klmn,isppol)=pawrhoij(iatom)%rhoij_(klmn,isppol)+weight*ro11_re
668 202889034 : if (present(occ_k_2)) then
669 13235786 : ro11_im=cpi0(1,1)*cpj0(2,1)-cpi0(2,1)*cpj0(1,1)
670 13235786 : pawrhoij(iatom)%rhoij_(klmn,isppol)=pawrhoij(iatom)%rhoij_(klmn,isppol)-weight_2*ro11_im
671 : end if
672 : end do
673 : end do
674 : end do
675 : else
676 137796 : do iatom=1,my_natom
677 87714 : iatom1=iatom;if (paral_atom) iatom1=my_atmtab(iatom)
678 87714 : iatm=atindx(iatom1)
679 839508 : do jlmn=1,pawrhoij(iatom)%lmn_size
680 701712 : j0lmn=jlmn*(jlmn-1)/2
681 2105136 : cpj0(1:cplex,1)=cwaveprj(iatm,1)%cp(1:cplex,jlmn)
682 3947130 : do ilmn=1,jlmn
683 3157704 : klmn=j0lmn+ilmn
684 3157704 : klmn_re=cplex_rhoij*(klmn-1)+1
685 9473112 : cpi0(1:cplex,1)=cwaveprj1(iatm,1)%cp(1:cplex,ilmn)
686 : ro11_re=zero
687 9473112 : do iplex=1,cplex
688 9473112 : ro11_re=ro11_re+cpi0(iplex,1)*cpj0(iplex,1)
689 : end do
690 3157704 : pawrhoij(iatom)%rhoij_(klmn_re,isppol)=pawrhoij(iatom)%rhoij_(klmn_re,isppol)+weight*ro11_re
691 3157704 : if (present(occ_k_2)) then
692 0 : ro11_im=cpi0(1,1)*cpj0(2,1)-cpi0(2,1)*cpj0(1,1)
693 0 : pawrhoij(iatom)%rhoij_(klmn_re,isppol)=pawrhoij(iatom)%rhoij_(klmn_re,isppol)-weight_2*ro11_im
694 : end if
695 3859416 : if (compute_impart_cplex) then
696 3157704 : klmn_im=klmn_re+1
697 3157704 : ro11_im=cpi0(1,1)*cpj0(2,1)-cpi0(2,1)*cpj0(1,1)
698 3157704 : pawrhoij(iatom)%rhoij_(klmn_im,isppol)=pawrhoij(iatom)%rhoij_(klmn_im,isppol)+weight*ro11_im
699 3157704 : if (present(occ_k_2)) then
700 0 : pawrhoij(iatom)%rhoij_(klmn_im,isppol)=pawrhoij(iatom)%rhoij_(klmn_im,isppol)+weight_2*ro11_re
701 : end if
702 : end if
703 : end do
704 : end do
705 : end do
706 : end if
707 : else ! nspinor=2
708 446150 : do iatom=1,my_natom
709 290586 : iatom1=iatom;if (paral_atom) iatom1=my_atmtab(iatom)
710 290586 : iatm=atindx(iatom1)
711 290586 : cplex_rhoij=pawrhoij(iatom)%cplex_rhoij
712 290586 : nspden_rhoij=pawrhoij(iatom)%nspden
713 4709836 : do jlmn=1,pawrhoij(iatom)%lmn_size
714 4263686 : j0lmn=jlmn*(jlmn-1)/2
715 12791058 : cpj0(1:cplex,1)=cwaveprj(iatm,1)%cp(1:cplex,jlmn)
716 12791058 : cpj0(1:cplex,2)=cwaveprj(iatm,2)%cp(1:cplex,jlmn)
717 41124266 : do ilmn=1,jlmn
718 36569994 : klmn=j0lmn+ilmn
719 36569994 : klmn_re=cplex_rhoij*(klmn-1)+1
720 109709982 : cpi0(1:cplex,1)=cwaveprj1(iatm,1)%cp(1:cplex,ilmn)
721 109709982 : cpi0(1:cplex,2)=cwaveprj1(iatm,2)%cp(1:cplex,ilmn)
722 : ro11_re=zero;ro22_re=zero
723 109709982 : ro12_re=zero;ro21_re=zero
724 109709982 : ro12_im=zero;ro21_im=zero
725 109709982 : do iplex=1,cplex
726 73139988 : ro11_re=ro11_re+cpi0(iplex,1)*cpj0(iplex,1)
727 109709982 : ro22_re=ro22_re+cpi0(iplex,2)*cpj0(iplex,2)
728 : end do
729 36569994 : pawrhoij(iatom)%rhoij_(klmn_re,1)=pawrhoij(iatom)%rhoij_(klmn_re,1)+weight*(ro11_re+ro22_re)
730 36569994 : if (nspden_rhoij>1) then
731 108355662 : do iplex=1,cplex
732 72237108 : ro12_re=ro12_re+cpi0(iplex,2)*cpj0(iplex,1)
733 108355662 : ro21_re=ro21_re+cpi0(iplex,1)*cpj0(iplex,2)
734 : end do
735 36118554 : if (.not.use_zeromag) then
736 26823990 : pawrhoij(iatom)%rhoij_(klmn_re,4)=pawrhoij(iatom)%rhoij_(klmn_re,4)+weight*(ro11_re-ro22_re)
737 26823990 : pawrhoij(iatom)%rhoij_(klmn_re,2)=pawrhoij(iatom)%rhoij_(klmn_re,2)+weight*(ro12_re+ro21_re)
738 : end if
739 36118554 : if (cplex==2) then
740 : !Important note: the present implementation follows eq(15) in Hobbs et al, PRB 62, 11556(2000)
741 : ! rho^alpha,beta_ij = Sum[<Psi^beta|pi><pj|Psi^alpha] (alpha and beta exponents inverted)
742 36118554 : ro12_im=cpi0(1,2)*cpj0(2,1)-cpi0(2,2)*cpj0(1,1)
743 36118554 : ro21_im=cpi0(1,1)*cpj0(2,2)-cpi0(2,1)*cpj0(1,2)
744 36118554 : if (.not.use_zeromag) then
745 26823990 : pawrhoij(iatom)%rhoij_(klmn_re,3)=pawrhoij(iatom)%rhoij_(klmn_re,3)+weight*(ro21_im-ro12_im)
746 : end if
747 : end if
748 : end if
749 36569994 : if (present(occ_k_2)) then
750 844800 : ro11_im=cpi0(1,1)*cpj0(2,1)-cpi0(2,1)*cpj0(1,1)
751 844800 : ro22_im=cpi0(1,2)*cpj0(2,2)-cpi0(2,2)*cpj0(1,2)
752 844800 : pawrhoij(iatom)%rhoij_(klmn_re,1)=pawrhoij(iatom)%rhoij_(klmn_re,1)+weight_2*(-ro11_im-ro22_im)
753 844800 : pawrhoij(iatom)%rhoij_(klmn_re,2)=pawrhoij(iatom)%rhoij_(klmn_re,2)+weight_2*(-ro21_im-ro12_im)
754 844800 : pawrhoij(iatom)%rhoij_(klmn_re,3)=pawrhoij(iatom)%rhoij_(klmn_re,3)+weight_2*(-ro12_re+ro21_re)
755 844800 : pawrhoij(iatom)%rhoij_(klmn_re,4)=pawrhoij(iatom)%rhoij_(klmn_re,4)+weight_2*(-ro11_im+ro22_im)
756 : end if
757 40833680 : if (compute_impart) then
758 28468869 : klmn_im=klmn_re+1
759 28468869 : if (nspden_rhoij>1) pawrhoij(iatom)%rhoij_(klmn_im,3)=pawrhoij(iatom)%rhoij_(klmn_im,3)+weight*(ro12_re-ro21_re)
760 28468869 : if (cplex==2) then
761 28468869 : ro11_im=cpi0(1,1)*cpj0(2,1)-cpi0(2,1)*cpj0(1,1)
762 28468869 : ro22_im=cpi0(1,2)*cpj0(2,2)-cpi0(2,2)*cpj0(1,2)
763 28468869 : pawrhoij(iatom)%rhoij_(klmn_im,1)=pawrhoij(iatom)%rhoij_(klmn_im,1)+weight*(ro11_im+ro22_im)
764 28468869 : if (nspden_rhoij>1) then
765 28468869 : pawrhoij(iatom)%rhoij_(klmn_im,4)=pawrhoij(iatom)%rhoij_(klmn_im,4)+weight*(ro11_im-ro22_im)
766 28468869 : pawrhoij(iatom)%rhoij_(klmn_im,2)=pawrhoij(iatom)%rhoij_(klmn_im,2)+weight*(ro12_im+ro21_im)
767 : end if
768 28468869 : if (present(occ_k_2)) then
769 844800 : pawrhoij(iatom)%rhoij_(klmn_im,1)=pawrhoij(iatom)%rhoij_(klmn_im,1)+weight_2*( ro11_re+ro22_re)
770 844800 : pawrhoij(iatom)%rhoij_(klmn_im,2)=pawrhoij(iatom)%rhoij_(klmn_im,2)+weight_2*( ro21_re+ro12_re)
771 844800 : pawrhoij(iatom)%rhoij_(klmn_im,3)=pawrhoij(iatom)%rhoij_(klmn_im,3)+weight_2*(-ro12_im+ro21_im)
772 844800 : pawrhoij(iatom)%rhoij_(klmn_im,4)=pawrhoij(iatom)%rhoij_(klmn_im,4)+weight_2*( ro11_re-ro22_re)
773 : end if
774 : end if
775 : end if
776 : end do
777 : end do
778 : end do
779 : end if
780 :
781 631685 : else if (option==2) then
782 :
783 : ! ==================================================================
784 : ! === OPTION 2: Accumulate (n,k) contribution to 1st-order rhoij ===
785 : ! ==================================================================
786 :
787 : ! Accumulate (n,k) contribution to rhoij1
788 : ! due to derivative of wave-function
789 617397 : compute_impart=(pawrhoij(1)%qphase==2)
790 617397 : compute_impart_cplex=(compute_impart.and.(cplex==2))
791 617397 : if (nspinor==1) then
792 2032164 : do iatom=1,my_natom
793 1414767 : iatom1=iatom;if (paral_atom) iatom1=my_atmtab(iatom)
794 1414767 : iatm=atindx(iatom1)
795 1414767 : cplex_rhoij=pawrhoij(iatom)%cplex_rhoij
796 1414767 : iq0=cplex_rhoij*pawrhoij(iatom)%lmn2_size
797 13691020 : do jlmn=1,pawrhoij(iatom)%lmn_size
798 11658856 : j0lmn=jlmn*(jlmn-1)/2
799 34976568 : cpj0(1:2,1)=cwaveprj (iatm,1)%cp(1:2,jlmn)
800 34976568 : cpj1(1:2,1)=cwaveprj1(iatm,1)%cp(1:2,jlmn)
801 68785131 : do ilmn=1,jlmn
802 55711508 : klmn=j0lmn+ilmn
803 55711508 : klmn_re=cplex_rhoij*(klmn-1)+1
804 167134524 : cpi0(1:2,1)=cwaveprj (iatm,1)%cp(1:2,ilmn)
805 167134524 : cpi1(1:2,1)=cwaveprj1(iatm,1)%cp(1:2,ilmn)
806 55711508 : ro11_re=zero
807 167134524 : do iplex=1,cplex
808 167134524 : ro11_re=ro11_re+cpi0(iplex,1)*cpj1(iplex,1)+cpj0(iplex,1)*cpi1(iplex,1)
809 : end do
810 55711508 : pawrhoij(iatom)%rhoij_(klmn_re,isppol)=pawrhoij(iatom)%rhoij_(klmn_re,isppol)+weight*ro11_re
811 67370364 : if (compute_impart_cplex) then
812 31140504 : klmn_im=klmn_re+iq0
813 31140504 : ro11_im=cpi0(1,1)*cpj1(2,1)-cpi0(2,1)*cpj1(1,1)+cpj0(1,1)*cpi1(2,1)-cpj0(2,1)*cpi1(1,1)
814 31140504 : pawrhoij(iatom)%rhoij_(klmn_im,isppol)=pawrhoij(iatom)%rhoij_(klmn_im,isppol)+weight*ro11_im
815 : end if
816 : end do
817 : end do
818 : end do
819 : else ! nspinor=2
820 0 : do iatom=1,my_natom
821 0 : iatom1=iatom;if (paral_atom) iatom1=my_atmtab(iatom)
822 0 : iatm=atindx(iatom1)
823 0 : cplex_rhoij=pawrhoij(iatom)%cplex_rhoij
824 0 : nspden_rhoij=pawrhoij(iatom)%nspden
825 0 : iq0=cplex_rhoij*pawrhoij(iatom)%lmn2_size
826 0 : do jlmn=1,pawrhoij(iatom)%lmn_size
827 0 : j0lmn=jlmn*(jlmn-1)/2
828 0 : cpj0(1:2,1)=cwaveprj (iatm,1)%cp(1:2,jlmn)
829 0 : cpj0(1:2,2)=cwaveprj (iatm,2)%cp(1:2,jlmn)
830 0 : cpj1(1:2,1)=cwaveprj1(iatm,1)%cp(1:2,jlmn)
831 0 : cpj1(1:2,2)=cwaveprj1(iatm,2)%cp(1:2,jlmn)
832 0 : do ilmn=1,jlmn
833 0 : klmn=j0lmn+ilmn
834 0 : klmn_re=cplex_rhoij*(klmn-1)+1
835 0 : cpi0(1:2,1)=cwaveprj (iatm,1)%cp(1:2,ilmn)
836 0 : cpi0(1:2,2)=cwaveprj (iatm,2)%cp(1:2,ilmn)
837 0 : cpi1(1:2,1)=cwaveprj1(iatm,1)%cp(1:2,ilmn)
838 0 : cpi1(1:2,2)=cwaveprj1(iatm,2)%cp(1:2,ilmn)
839 0 : ro11_re=zero;ro22_re=zero
840 0 : ro12_re=zero;ro21_re=zero
841 0 : ro12_im=zero;ro21_im=zero
842 0 : do iplex=1,cplex
843 0 : ro11_re=ro11_re+cpj0(iplex,1)*cpi1(iplex,1)+cpi0(iplex,1)*cpj1(iplex,1)
844 0 : ro22_re=ro22_re+cpj0(iplex,2)*cpi1(iplex,2)+cpi0(iplex,2)*cpj1(iplex,2)
845 : end do
846 0 : pawrhoij(iatom)%rhoij_(klmn_re,1)=pawrhoij(iatom)%rhoij_(klmn_re,1)+weight*(ro11_re+ro22_re)
847 0 : if (nspden_rhoij>1) then
848 0 : do iplex=1,cplex
849 0 : ro12_re=ro12_re+cpj0(iplex,1)*cpi1(iplex,2)+cpi0(iplex,2)*cpj1(iplex,1)
850 0 : ro21_re=ro21_re+cpj0(iplex,2)*cpi1(iplex,1)+cpi0(iplex,1)*cpj1(iplex,2)
851 : end do
852 0 : if (.not.use_zeromag) then
853 0 : pawrhoij(iatom)%rhoij_(klmn_re,4)=pawrhoij(iatom)%rhoij_(klmn_re,4)+weight*(ro11_re-ro22_re)
854 0 : pawrhoij(iatom)%rhoij_(klmn_re,2)=pawrhoij(iatom)%rhoij_(klmn_re,2)+weight*(ro12_re+ro21_re)
855 : end if
856 0 : if (cplex==2) then
857 : !Important note: the present implementation follows eq(15) in Hobbs et al, PRB 62, 11556(2000)
858 : ! rho^alpha,beta_ij = Sum[<Psi^beta|pi><pj|Psi^alpha] (alpha and beta exponents inverted)
859 0 : ro12_im=cpj0(2,1)*cpi1(1,2)-cpi1(2,2)*cpj0(1,1)+cpi0(1,2)*cpj1(2,1)-cpj1(1,1)*cpi0(2,2)
860 0 : ro21_im=cpj0(2,2)*cpi1(1,1)-cpi1(2,1)*cpj0(1,2)+cpi0(1,1)*cpj1(2,2)-cpj1(1,2)*cpi0(2,1)
861 0 : if (.not.use_zeromag) then
862 0 : pawrhoij(iatom)%rhoij_(klmn_re,3)=pawrhoij(iatom)%rhoij_(klmn_re,3)+weight*(ro21_im-ro12_im)
863 : end if
864 : end if
865 : end if
866 0 : if (compute_impart) then
867 0 : klmn_im=klmn_re+iq0
868 0 : if (nspden_rhoij>1) pawrhoij(iatom)%rhoij_(klmn_re,3)=pawrhoij(iatom)%rhoij_(klmn_re,3)+weight*(ro12_re-ro21_re)
869 0 : if (cplex==2) then
870 0 : ro11_im=cpj0(2,1)*cpi1(1,1)-cpi1(2,1)*cpj0(1,1)+cpi0(1,1)*cpj1(2,1)-cpj1(1,1)*cpi0(2,1)
871 0 : ro22_im=cpj0(2,2)*cpi1(1,2)-cpi1(2,2)*cpj0(1,2)+cpi0(1,2)*cpj1(2,2)-cpj1(1,2)*cpi0(2,2)
872 0 : pawrhoij(iatom)%rhoij_(klmn_im,1)=pawrhoij(iatom)%rhoij_(klmn_im,1)+weight*(ro11_im+ro22_im)
873 0 : if (nspden_rhoij>1) then
874 0 : pawrhoij(iatom)%rhoij_(klmn_im,4)=pawrhoij(iatom)%rhoij_(klmn_im,4)+weight*(ro11_im-ro22_im)
875 0 : pawrhoij(iatom)%rhoij_(klmn_im,2)=pawrhoij(iatom)%rhoij_(klmn_re,2)+weight*(ro12_im+ro21_im)
876 : end if
877 : end if
878 : end if
879 : end do
880 : end do
881 : end do
882 : end if
883 :
884 : ! Accumulate (n,k) contribution to rhoij1
885 : ! due to derivative of projectors
886 617397 : if (ipert/=natom+2) then
887 579733 : compute_impart=(pawrhoij(1)%cplex_rhoij==2)
888 579733 : compute_impart_cplex=(compute_impart.and.(cplex==2))
889 579733 : substract_diagonal=(ipert==natom+3)
890 579733 : if (nspinor==1) then
891 1907396 : do iatom=1,my_natom
892 1327663 : iatom1=iatom;if (paral_atom) iatom1=my_atmtab(iatom)
893 1327663 : iatm=atindx(iatom1)
894 1327663 : if (ipert<=natom.and.iatom/=ipert) cycle
895 689217 : cplex_rhoij=pawrhoij(iatom)%cplex_rhoij
896 7062846 : do jlmn=1,pawrhoij(iatom)%lmn_size
897 5793896 : j0lmn=jlmn*(jlmn-1)/2
898 17381688 : cpj0 (1:2,1) =cwaveprj(iatm,1)%cp (1:2 ,jlmn)
899 17381688 : dcpj0(1:2,1,1)=cwaveprj(iatm,1)%dcp(1:2,1,jlmn)
900 35805619 : do ilmn=1,jlmn
901 28684060 : klmn=j0lmn+ilmn
902 28684060 : klmn_re=cplex_rhoij*(klmn-1)+1
903 86052180 : cpi0 (1:2,1) =cwaveprj(iatm,1)%cp (1:2 ,ilmn)
904 86052180 : dcpi0(1:2,1,1)=cwaveprj(iatm,1)%dcp(1:2,1,ilmn)
905 28684060 : ro11_re=zero
906 86052180 : do iplex=1,cplex
907 86052180 : ro11_re=ro11_re+dcpi0(iplex,1,1)*cpj0(iplex,1)+cpi0(iplex,1)*dcpj0(iplex,1,1)
908 : end do
909 28684060 : if (substract_diagonal) then
910 8953848 : do iplex=1,cplex
911 8953848 : ro11_re=ro11_re-cpi0(iplex,1)*cpj0(iplex,1)
912 : end do
913 : end if
914 28684060 : pawrhoij(iatom)%rhoij_(klmn_re,isppol)=pawrhoij(iatom)%rhoij_(klmn_re,isppol)+weight*ro11_re
915 : ! This imaginary part does not have to be computed
916 : ! It is cancelled because rho_ij+rho_ji is stored in rho_ij
917 34477956 : if (compute_impart_cplex) then
918 0 : klmn_im=klmn_re+1
919 0 : ro11_im=dcpi0(1,1,1)*cpj0(2,1)-dcpi0(2,1,1)*cpj0(1,1)+cpi0(1,1)*dcpj0(2,1,1)-cpi0(2,1)*dcpj0(1,1,1)
920 0 : if (substract_diagonal) then
921 0 : ro11_im=ro11_im-cpi0(1,1)*cpj0(2,1)+cpi0(2,1)*cpj0(1,1)
922 : end if
923 0 : pawrhoij(iatom)%rhoij_(klmn_im,isppol)=pawrhoij(iatom)%rhoij_(klmn_im,isppol)+weight*ro11_im
924 : end if
925 : end do
926 : end do
927 : end do
928 : else ! nspinor=2
929 0 : do iatom=1,my_natom
930 0 : iatom1=iatom;if (paral_atom) iatom1=my_atmtab(iatom)
931 0 : iatm=atindx(iatom1)
932 0 : if (ipert<=natom.and.iatom/=ipert) cycle
933 0 : cplex_rhoij=pawrhoij(iatom)%cplex_rhoij
934 0 : nspden_rhoij=pawrhoij(iatom)%nspden
935 0 : do jlmn=1,pawrhoij(iatom)%lmn_size
936 0 : j0lmn=jlmn*(jlmn-1)/2
937 0 : cpj0 (1:2,1) =cwaveprj(iatm,1)%cp (1:2 ,jlmn)
938 0 : dcpj0(1:2,1,1)=cwaveprj(iatm,1)%dcp(1:2,1,jlmn)
939 0 : cpj0 (1:2,2) =cwaveprj(iatm,2)%cp (1:2 ,jlmn)
940 0 : dcpj0(1:2,2,1)=cwaveprj(iatm,2)%dcp(1:2,1,jlmn)
941 0 : do ilmn=1,jlmn
942 0 : klmn=j0lmn+ilmn
943 0 : klmn_re=cplex_rhoij*(klmn-1)+1
944 0 : cpi0 (1:2,1) =cwaveprj(iatm,1)%cp (1:2 ,ilmn)
945 0 : dcpi0(1:2,1,1)=cwaveprj(iatm,1)%dcp(1:2,1,ilmn)
946 0 : cpi0 (1:2,2) =cwaveprj(iatm,2)%cp (1:2 ,ilmn)
947 0 : dcpi0(1:2,2,1)=cwaveprj(iatm,2)%dcp(1:2,1,ilmn)
948 0 : ro11_re=zero;ro22_re=zero
949 0 : ro12_re=zero;ro21_re=zero
950 0 : ro12_im=zero;ro21_im=zero
951 0 : do iplex=1,cplex
952 0 : ro11_re=ro11_re+dcpi0(iplex,1,1)*cpj0(iplex,1)+cpi0(iplex,1)*dcpj0(iplex,1,1)
953 0 : ro22_re=ro22_re+dcpi0(iplex,2,1)*cpj0(iplex,2)+cpi0(iplex,2)*dcpj0(iplex,2,1)
954 : end do
955 0 : if (substract_diagonal) then
956 0 : do iplex=1,cplex
957 0 : ro11_re=ro11_re-cpi0(iplex,1)*cpj0(iplex,1)
958 0 : ro22_re=ro22_re-cpi0(iplex,2)*cpj0(iplex,2)
959 : end do
960 : end if
961 0 : pawrhoij(iatom)%rhoij_(klmn_re,1)=pawrhoij(iatom)%rhoij_(klmn_re,1)+weight*(ro11_re+ro22_re)
962 0 : if (nspden_rhoij>1) then
963 0 : do iplex=1,cplex
964 0 : ro12_re=ro12_re+dcpi0(iplex,2,1)*cpj0(iplex,1)+cpi0(iplex,2)*dcpj0(iplex,1,1)
965 0 : ro21_re=ro21_re+dcpi0(iplex,1,1)*cpj0(iplex,2)+cpi0(iplex,1)*dcpj0(iplex,2,1)
966 : end do
967 0 : if (substract_diagonal) then
968 0 : do iplex=1,cplex
969 0 : ro12_re=ro12_re-cpi0(iplex,2)*cpj0(iplex,1)
970 0 : ro21_re=ro21_re-cpi0(iplex,1)*cpj0(iplex,2)
971 : end do
972 : end if
973 0 : if (.not.use_zeromag) then
974 0 : pawrhoij(iatom)%rhoij_(klmn_re,4)=pawrhoij(iatom)%rhoij_(klmn_re,4)+weight*(ro11_re-ro22_re)
975 0 : pawrhoij(iatom)%rhoij_(klmn_re,2)=pawrhoij(iatom)%rhoij_(klmn_re,2)+weight*(ro12_re+ro21_re)
976 : end if
977 0 : if (cplex==2) then
978 : !Important note: the present implementation follows eq(15) in Hobbs et al, PRB 62, 11556(2000)
979 : ! rho^alpha,beta_ij = Sum[<Psi^beta|pi><pj|Psi^alpha] (alpha and beta exponents inverted)
980 0 : ro12_im=dcpi0(1,2,1)*cpj0(2,1)-dcpi0(2,2,1)*cpj0(1,1)+cpi0(1,2)*dcpj0(2,1,1)-cpi0(2,2)*dcpj0(1,1,1)
981 0 : ro21_im=dcpi0(1,1,1)*cpj0(2,2)-dcpi0(2,1,1)*cpj0(1,2)+cpi0(1,1)*dcpj0(2,2,1)-cpi0(2,1)*dcpj0(1,2,1)
982 0 : if (substract_diagonal) then
983 0 : ro12_im=ro12_im-cpi0(1,2)*cpj0(2,1)+cpi0(2,2)*cpj0(1,1)
984 0 : ro21_im=ro21_im-cpi0(1,1)*cpj0(2,2)+cpi0(2,1)*cpj0(1,2)
985 : end if
986 0 : if (.not.use_zeromag) then
987 0 : pawrhoij(iatom)%rhoij_(klmn_re,3)=pawrhoij(iatom)%rhoij_(klmn_re,3)+weight*(ro21_im-ro12_im)
988 : end if
989 : end if
990 : end if
991 0 : if (compute_impart) then
992 0 : klmn_im=klmn_re+1
993 0 : if (nspden_rhoij>1) pawrhoij(iatom)%rhoij_(klmn_im,3)=pawrhoij(iatom)%rhoij_(klmn_im,3)+weight*(ro12_re-ro21_re)
994 0 : if (cplex==2) then
995 0 : ro11_im=dcpi0(1,1,1)*cpj0(2,1)-dcpi0(2,1,1)*cpj0(1,1)+cpi0(1,1)*dcpj0(2,1,1)-cpi0(2,1)*dcpj0(1,1,1)
996 0 : ro22_im=dcpi0(1,2,1)*cpj0(2,2)-dcpi0(2,2,1)*cpj0(1,2)+cpi0(1,2)*dcpj0(2,2,1)-cpi0(2,2)*dcpj0(1,2,1)
997 0 : if (substract_diagonal) then
998 0 : ro11_im=ro11_im-cpi0(1,1)*cpj0(2,1)+cpi0(2,1)*cpj0(1,1)
999 0 : ro22_im=ro22_im-cpi0(1,2)*cpj0(2,2)+cpi0(2,2)*cpj0(1,2)
1000 : end if
1001 0 : pawrhoij(iatom)%rhoij_(klmn_im,1)=pawrhoij(iatom)%rhoij_(klmn_im,1)+weight*(ro11_im+ro22_im)
1002 0 : if (nspden_rhoij>1) then
1003 0 : pawrhoij(iatom)%rhoij_(klmn_im,4)=pawrhoij(iatom)%rhoij_(klmn_im,4)+weight*(ro11_im-ro22_im)
1004 0 : pawrhoij(iatom)%rhoij_(klmn_im,2)=pawrhoij(iatom)%rhoij_(klmn_im,2)+weight*(ro12_im+ro21_im)
1005 : end if
1006 : end if
1007 : end if
1008 : end do
1009 : end do
1010 : end do
1011 : end if
1012 : end if
1013 :
1014 14288 : else if (option==3) then
1015 :
1016 : ! ==================================================================
1017 : ! === OPTION 3: Accumulate (n,k) contribution to drhoij/dr =========
1018 : ! ==================================================================
1019 :
1020 14288 : compute_impart=((.not.use_timerev).and.(pawrhoij(1)%cplex_rhoij==2))
1021 0 : compute_impart_cplex=((compute_impart).and.(cplex==2))
1022 14288 : if (nspinor==1) then
1023 40536 : do iatom=1,my_natom
1024 26248 : iatom1=iatom;if (paral_atom) iatom1=my_atmtab(iatom)
1025 26248 : iatm=atindx(iatom1)
1026 26248 : cplex_rhoij=pawrhoij(iatom)%cplex_rhoij
1027 291160 : do jlmn=1,pawrhoij(iatom)%lmn_size
1028 250624 : j0lmn=jlmn*(jlmn-1)/2
1029 751872 : cpj0(1:cplex,1) =cwaveprj(iatm,1)%cp (1:cplex,jlmn)
1030 5290048 : dcpj0(1:cplex,1,1:ncpgr)=cwaveprj(iatm,1)%dcp(1:cplex,1:ncpgr,jlmn)
1031 1781416 : do ilmn=1,jlmn
1032 1504544 : klmn=j0lmn+ilmn
1033 1504544 : klmn_re=cplex_rhoij*(klmn-1)+1
1034 4513632 : cpi0(1:cplex,1) =cwaveprj(iatm,1)%cp (1:cplex,ilmn)
1035 31367552 : dcpi0(1:cplex,1,1:ncpgr)=cwaveprj(iatm,1)%dcp(1:cplex,1:ncpgr,ilmn)
1036 11458880 : do mu=1,ncpgr
1037 : ro11_re=zero
1038 29863008 : do iplex=1,cplex
1039 29863008 : ro11_re=ro11_re+dcpi0(iplex,1,mu)*cpj0(iplex,1)+cpi0(iplex,1)*dcpj0(iplex,1,mu)
1040 : end do
1041 11458880 : pawrhoij(iatom)%grhoij(mu,klmn_re,isppol)=pawrhoij(iatom)%grhoij(mu,klmn_re,isppol)+weight*ro11_re
1042 : end do
1043 1755168 : if (compute_impart_cplex) then
1044 0 : klmn_im=klmn_re+1
1045 0 : do mu=1,ncpgr
1046 0 : ro11_im=dcpi0(1,1,mu)*cpj0(2,1)+cpi0(1,1)*dcpj0(2,1,mu)-dcpi0(2,1,mu)*cpj0(1,1)-cpi0(2,1)*dcpj0(1,1,mu)
1047 0 : pawrhoij(iatom)%grhoij(mu,klmn_im,isppol)=pawrhoij(iatom)%grhoij(mu,klmn_im,isppol)+weight*ro11_im
1048 : end do
1049 : end if
1050 : end do
1051 : end do
1052 : end do
1053 : else ! nspinor=2
1054 0 : do iatom=1,my_natom
1055 0 : iatom1=iatom;if (paral_atom) iatom1=my_atmtab(iatom)
1056 0 : iatm=atindx(iatom1)
1057 0 : cplex_rhoij=pawrhoij(iatom)%cplex_rhoij
1058 0 : nspden_rhoij=pawrhoij(iatom)%nspden
1059 0 : do jlmn=1,pawrhoij(iatom)%lmn_size
1060 0 : j0lmn=jlmn*(jlmn-1)/2
1061 0 : cpj0(1:cplex,1) =cwaveprj(iatm,1)%cp (1:cplex,jlmn)
1062 0 : cpj0(1:cplex,2) =cwaveprj(iatm,2)%cp (1:cplex,jlmn)
1063 0 : dcpj0(1:cplex,1,1:ncpgr)=cwaveprj(iatm,1)%dcp(1:cplex,1:ncpgr,jlmn)
1064 0 : dcpj0(1:cplex,2,1:ncpgr)=cwaveprj(iatm,2)%dcp(1:cplex,1:ncpgr,jlmn)
1065 0 : do ilmn=1,jlmn
1066 0 : klmn=j0lmn+ilmn
1067 0 : klmn_re=cplex_rhoij*(klmn-1)+1
1068 0 : klmn_im=klmn_re+1
1069 0 : cpi0(1:cplex,1) =cwaveprj(iatm,1)%cp (1:cplex,ilmn)
1070 0 : cpi0(1:cplex,2) =cwaveprj(iatm,2)%cp (1:cplex,ilmn)
1071 0 : dcpi0(1:cplex,1,1:ncpgr)=cwaveprj(iatm,1)%dcp(1:cplex,1:ncpgr,ilmn)
1072 0 : dcpi0(1:cplex,2,1:ncpgr)=cwaveprj(iatm,2)%dcp(1:cplex,1:ncpgr,ilmn)
1073 0 : do mu=1,ncpgr
1074 : ro11_re=zero;ro22_re=zero
1075 0 : ro12_re=zero;ro21_re=zero
1076 0 : ro12_im=zero;ro21_im=zero
1077 0 : do iplex=1,cplex
1078 0 : ro11_re=ro11_re+dcpi0(iplex,1,mu)*cpj0(iplex,1)+cpi0(iplex,1)*dcpj0(iplex,1,mu)
1079 0 : ro22_re=ro22_re+dcpi0(iplex,2,mu)*cpj0(iplex,2)+cpi0(iplex,2)*dcpj0(iplex,2,mu)
1080 : end do
1081 0 : pawrhoij(iatom)%grhoij(mu,klmn_re,1)=pawrhoij(iatom)%grhoij(mu,klmn_re,1)+weight*(ro11_re+ro22_re)
1082 0 : if (nspden_rhoij>1) then
1083 0 : do iplex=1,cplex
1084 0 : ro12_re=ro12_re+dcpi0(iplex,2,mu)*cpj0(iplex,1)+cpi0(iplex,2)*dcpj0(iplex,1,mu)
1085 0 : ro21_re=ro21_re+dcpi0(iplex,1,mu)*cpj0(iplex,2)+cpi0(iplex,1)*dcpj0(iplex,2,mu)
1086 : end do
1087 0 : if (.not.use_zeromag) then
1088 0 : pawrhoij(iatom)%grhoij(mu,klmn_re,4)=pawrhoij(iatom)%grhoij(mu,klmn_re,4)+weight*(ro11_re-ro22_re)
1089 0 : pawrhoij(iatom)%grhoij(mu,klmn_re,2)=pawrhoij(iatom)%grhoij(mu,klmn_re,2)+weight*(ro12_re+ro21_re)
1090 : end if
1091 0 : if (cplex==2) then
1092 : !Important note: the present implementation follows eq(15) in Hobbs et al, PRB 62, 11556(2000)
1093 : ! rho^alpha,beta_ij = Sum[<Psi^beta|pi><pj|Psi^alpha] (alpha and beta exponents inverted)
1094 0 : ro12_im=dcpi0(1,2,mu)*cpj0(2,1)+cpi0(1,2)*dcpj0(2,1,mu)-dcpi0(2,2,mu)*cpj0(1,1)-cpi0(2,2)*dcpj0(1,1,mu)
1095 0 : ro21_im=dcpi0(1,1,mu)*cpj0(2,2)+cpi0(1,1)*dcpj0(2,2,mu)-dcpi0(2,1,mu)*cpj0(1,2)-cpi0(2,1)*dcpj0(1,2,mu)
1096 0 : if (.not.use_zeromag) then
1097 0 : pawrhoij(iatom)%grhoij(mu,klmn_re,3)=pawrhoij(iatom)%grhoij(mu,klmn_re,3)+weight*(ro21_im-ro12_im)
1098 : end if
1099 : end if
1100 : end if
1101 0 : if (compute_impart) then
1102 0 : if (nspden_rhoij>1) then
1103 0 : pawrhoij(iatom)%grhoij(mu,klmn_im,3)=pawrhoij(iatom)%grhoij(mu,klmn_im,3)+weight*(ro12_re-ro21_re)
1104 : end if
1105 0 : if (cplex==2) then
1106 0 : ro11_im=dcpi0(1,1,mu)*cpj0(2,1)+cpi0(1,1)*dcpj0(2,1,mu)-dcpi0(2,1,mu)*cpj0(1,1)-cpi0(2,1)*dcpj0(1,1,mu)
1107 0 : ro22_im=dcpi0(1,2,mu)*cpj0(2,2)+cpi0(1,2)*dcpj0(2,2,mu)-dcpi0(2,2,mu)*cpj0(1,2)-cpi0(2,2)*dcpj0(1,2,mu)
1108 0 : pawrhoij(iatom)%grhoij(mu,klmn_im,1)=pawrhoij(iatom)%grhoij(mu,klmn_im,1)+weight*(ro11_im+ro22_im)
1109 0 : if (nspden_rhoij>1) then
1110 0 : pawrhoij(iatom)%grhoij(mu,klmn_im,4)=pawrhoij(iatom)%grhoij(mu,klmn_im,4)+weight*(ro11_im-ro22_im)
1111 0 : pawrhoij(iatom)%grhoij(mu,klmn_im,2)=pawrhoij(iatom)%grhoij(mu,klmn_im,2)+weight*(ro12_im+ro21_im)
1112 : end if
1113 : end if
1114 : end if
1115 : end do
1116 : end do
1117 : end do
1118 : end do
1119 : end if
1120 :
1121 : ! End
1122 : end if ! option
1123 :
1124 : !Destroy atom table used for parallelism
1125 1498832 : call free_my_atmtab(my_atmtab,my_atmtab_allocated)
1126 :
1127 : DBG_EXIT("COLL")
1128 :
1129 1498832 : end subroutine pawaccrhoij
1130 : !!***
1131 :
1132 : !----------------------------------------------------------------------
1133 :
1134 : !!****f* m_paw_occupancies/initrhoij
1135 : !! NAME
1136 : !! initrhoij
1137 : !!
1138 : !! FUNCTION
1139 : !! Initialize PAW rhoij occupancies (in packed storage)
1140 : !! from atomic ones
1141 : !!
1142 : !! INPUTS
1143 : !! cpxocc=1 if rhoij are real, 2 if they are complex
1144 : !! lexexch(ntypat)=l on which local exact-exchange is applied for a given type of atom
1145 : !! lpawu(ntypat)=l on which U is applied for a given type of atom (PAW+U)
1146 : !! mpi_atmtab(:)=--optional-- indexes of the atoms treated by current proc
1147 : !! comm_atom=--optional-- MPI communicator over atoms
1148 : !! my_natom=number of atoms treated by current processor
1149 : !! natom=number of atoms
1150 : !! nspden=number of spin-density components
1151 : !! nspinor=number of spinorial components
1152 : !! nsppol=1 for unpolarized, 2 for spin-polarized
1153 : !! ntypat=number of atom types
1154 : !! pawspnorb=flag: 1 if spin-orbit coupling is activated in PAW augmentation regions
1155 : !! pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data
1156 : !! (containing initial rhoij)
1157 : !! qphase=2 if rhoij have a exp(iqR) phase, 1 if not (typical use: response function at q<>0)
1158 : !! spinat(3,natom)=initial spin of each atom, in unit of hbar/2.
1159 : !! typat(natom)=type of each atom
1160 : !! === Optional arguments
1161 : !! ngrhoij=number of gradients to be allocated (OPTIONAL, default=0)
1162 : !! nlmnmix=number of rhoij elements to be mixed during SCF cycle (OPTIONAL, default=0)
1163 : !! use_rhoij_=1 if pawrhoij(:)%rhoij_ has to be allocated (OPTIONAL, default=0)
1164 : !! use_rhoijres=1 if pawrhoij(:)%rhoijres has to be allocated (OPTIONAL, default=0)
1165 :
1166 : !!
1167 : !! OUTPUT
1168 : !! pawrhoij(natom) <type(pawrhoij_type)>=rhoij quantities for each atom
1169 : !! in packed storage
1170 : !!
1171 : !! SOURCE
1172 :
1173 1217 : subroutine initrhoij(cpxocc,lexexch,lpawu,my_natom,natom,nspden,nspinor,nsppol,&
1174 1217 : & ntypat,pawrhoij,pawspnorb,pawtab,qphase,spinat,typat,&
1175 : & ngrhoij,nlmnmix,use_rhoij_,use_rhoijres,& ! optional arguments
1176 1217 : & mpi_atmtab,comm_atom) ! optional arguments (parallelism)
1177 :
1178 : !Arguments ---------------------------------------------
1179 : !scalars
1180 : integer,intent(in) :: cpxocc,my_natom,natom,nspden,nspinor,nsppol,ntypat,pawspnorb,qphase
1181 : integer,intent(in),optional :: comm_atom,ngrhoij,nlmnmix,use_rhoij_,use_rhoijres
1182 : character(len=500) :: message
1183 : !arrays
1184 : integer,intent(in) :: lexexch(ntypat),lpawu(ntypat)
1185 : integer,intent(in) :: typat(natom)
1186 : integer,optional,target,intent(in) :: mpi_atmtab(:)
1187 : real(dp),intent(in) :: spinat(3,natom)
1188 : type(pawrhoij_type),intent(inout) :: pawrhoij(my_natom)
1189 : type(pawtab_type),intent(in) :: pawtab(ntypat)
1190 :
1191 : !Local variables ---------------------------------------
1192 : !Arrays
1193 : !scalars
1194 : integer :: cplex_rhoij,iatom,iatom_rhoij,ilmn,ispden,itypat,j0lmn,jl,jlmn,jspden
1195 : integer :: klmn,klmn1,ln,lnspinat0,my_comm_atom
1196 : integer :: ngrhoij0,nlmnmix0,nselect,nselect1,nspden_rhoij,qphase_rhoij
1197 : integer :: use_rhoij_0,use_rhoijres0
1198 : real(dp) :: ratio,ro,roshift,zratio,zz
1199 : logical :: my_atmtab_allocated,paral_atom,spinat_zero,test_exexch,test_pawu,test_lnspinat
1200 : !arrays
1201 1217 : integer,pointer :: my_atmtab(:),lnspinat(:)
1202 1217 : real(dp),allocatable :: occ(:)
1203 :
1204 : !************************************************************************
1205 :
1206 : DBG_ENTER("COLL")
1207 :
1208 : !PAW+U and local exact-exchange restriction
1209 3042 : do itypat=1,ntypat
1210 3042 : if (lpawu(itypat)/=lexexch(itypat).and. lpawu(itypat)/=-1.and.lexexch(itypat)/=-1) then
1211 0 : message = ' lpawu must be equal to lexexch !'
1212 0 : ABI_ERROR(message)
1213 : end if
1214 : end do
1215 :
1216 : !Set up parallelism over atoms
1217 1217 : paral_atom=(present(comm_atom).and.(my_natom/=natom))
1218 1217 : nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
1219 1217 : my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
1220 1217 : call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,my_natom_ref=my_natom)
1221 :
1222 : call pawrhoij_inquire_dim(cplex_rhoij=cplex_rhoij,qphase_rhoij=qphase_rhoij,nspden_rhoij=nspden_rhoij,&
1223 1217 : & nspden=nspden,spnorb=pawspnorb,cpxocc=cpxocc,cplex=qphase)
1224 :
1225 1217 : ratio=one;if (nspden_rhoij==2) ratio=half
1226 13188 : spinat_zero=(all(abs(spinat(:,:))<tol10).or.(nspden_rhoij==4.and.nspden==1))
1227 :
1228 : ! if (nspden_rhoij==4.and.nspden==1.and.(.not.spinat_zero)) then
1229 : ! write(message,'(5a)') 'You are performing a unpolarized calculation (nspden==1)',ch10,&
1230 : !& 'but you start with a magnetization on atom (spinat/=0).',ch10,&
1231 : !& 'This is not expected and my produce unphysical results!'
1232 : ! ABI_WARNING(message)
1233 : ! end if
1234 :
1235 1217 : if (my_natom>0) then
1236 1075 : ngrhoij0=0;if (present(ngrhoij)) ngrhoij0=ngrhoij
1237 1075 : nlmnmix0=0;if (present(nlmnmix)) nlmnmix0=nlmnmix
1238 1075 : use_rhoij_0=0;if (present(use_rhoij_)) use_rhoij_0=use_rhoij_
1239 1075 : use_rhoijres0=0;if (present(use_rhoijres)) use_rhoijres0=use_rhoijres
1240 1075 : if (paral_atom) then
1241 : call pawrhoij_alloc(pawrhoij,cplex_rhoij,nspden_rhoij,nspinor,nsppol,typat,&
1242 : & ngrhoij=ngrhoij0,nlmnmix=nlmnmix0,use_rhoij_=use_rhoij_0,use_rhoijres=use_rhoijres0,&
1243 320 : & qphase=qphase_rhoij,pawtab=pawtab,comm_atom=my_comm_atom,mpi_atmtab=my_atmtab)
1244 : else
1245 : call pawrhoij_alloc(pawrhoij,cplex_rhoij,nspden_rhoij,nspinor,nsppol,typat,qphase=qphase_rhoij,&
1246 755 : & pawtab=pawtab,ngrhoij=ngrhoij0,nlmnmix=nlmnmix0,use_rhoij_=use_rhoij_0,use_rhoijres=use_rhoijres0)
1247 : end if
1248 : end if
1249 :
1250 3477 : do iatom_rhoij=1,my_natom
1251 2260 : iatom=iatom_rhoij;if (paral_atom) iatom=my_atmtab(iatom_rhoij)
1252 2260 : itypat=typat(iatom)
1253 2260 : nselect=0
1254 6780 : ABI_MALLOC(lnspinat,(pawtab(itypat)%basis_size))
1255 12088 : lnspinat=-1
1256 : ! Determine occupancies of each orbital
1257 2260 : if (nspden_rhoij==2) then
1258 1233 : ABI_MALLOC(occ,(pawtab(itypat)%basis_size))
1259 2576 : occ=zero
1260 6274 : do jlmn=1,pawtab(itypat)%lmn_size
1261 5863 : ln=pawtab(itypat)%indlmn(5,jlmn)
1262 5863 : klmn=jlmn*(jlmn+1)/2
1263 6274 : occ(ln)=occ(ln)+pawtab(itypat)%rhoij0(klmn)
1264 : end do
1265 2576 : do ln=1,pawtab(itypat)%basis_size
1266 2165 : if(pawtab(itypat)%orbitals(ln)==0.and.occ(ln)==1) lnspinat(ln)=ln
1267 2165 : if(pawtab(itypat)%orbitals(ln)==1.and.(occ(ln)>=1.and.occ(ln)<=5)) lnspinat(ln)=ln
1268 2165 : if(pawtab(itypat)%orbitals(ln)==2.and.(occ(ln)>=1.and.occ(ln)<=9)) lnspinat(ln)=ln
1269 2576 : if(pawtab(itypat)%orbitals(ln)==3.and.(occ(ln)>=1.and.occ(ln)<=13)) lnspinat(ln)=ln
1270 : end do
1271 411 : ABI_FREE(occ)
1272 : end if
1273 : lnspinat0=maxval(lnspinat)
1274 2260 : lnspinat0=-1
1275 :
1276 : ! Determine Z (trace of rhoij0 or part of it)
1277 2260 : zz=zero
1278 25886 : do jlmn=1,pawtab(itypat)%lmn_size
1279 23626 : jl=pawtab(itypat)%indlmn(1,jlmn)
1280 23626 : ln=pawtab(itypat)%indlmn(5,jlmn)
1281 23626 : j0lmn=jlmn*(jlmn-1)/2
1282 23626 : test_lnspinat=(lnspinat0==-1.or.lnspinat(ln)==ln)
1283 23626 : test_pawu=(lpawu(itypat)==-1.or.lpawu(itypat)==jl)
1284 23626 : test_exexch=(lexexch(itypat)==-1.or.lexexch(itypat)==jl)
1285 201618 : do ilmn=1,jlmn
1286 175732 : klmn=j0lmn+ilmn
1287 175732 : if ((ilmn==jlmn).and.test_pawu.and.test_exexch.and.test_lnspinat) &
1288 45388 : & zz=zz+pawtab(itypat)%rhoij0(klmn)
1289 : end do
1290 : end do
1291 :
1292 : ! Compute rhoij from tabulated value and magnetization
1293 5528 : do ispden=1,nspden_rhoij
1294 :
1295 3268 : zratio=zero
1296 3268 : roshift=one
1297 3268 : ratio=one
1298 3268 : if (nspden_rhoij==2) then
1299 822 : ratio=half
1300 822 : if ((spinat(3,iatom)>tol12.and.ispden==1).or.&
1301 : & (spinat(3,iatom)<tol12.and.ispden==2)) then
1302 411 : if(abs(zz)>tol12)then
1303 410 : zratio=two*abs(spinat(3,iatom))/zz
1304 : else
1305 : zratio=zero
1306 : end if
1307 : end if
1308 2446 : else if (nspden_rhoij==4.and.ispden>=2) then
1309 597 : roshift=zero
1310 597 : if(abs(zz)>tol12.and.(.not.spinat_zero)) then
1311 192 : zratio=spinat(ispden-1,iatom)/zz
1312 : else
1313 : zratio=zero
1314 : end if
1315 : end if
1316 :
1317 3268 : nselect=0;nselect1=1-cpxocc
1318 42985 : do jlmn=1,pawtab(itypat)%lmn_size
1319 37457 : jl=pawtab(itypat)%indlmn(1,jlmn)
1320 37457 : ln=pawtab(itypat)%indlmn(5,jlmn)
1321 37457 : j0lmn=jlmn*(jlmn-1)/2
1322 37457 : test_lnspinat=(lnspinat0==-1.or.lnspinat(ln)==ln)
1323 37457 : test_pawu=(lpawu(itypat)==-1.or.lpawu(itypat)==jl)
1324 37457 : test_exexch=(lexexch(itypat)==-1.or.lexexch(itypat)==jl)
1325 332309 : do ilmn=1,jlmn
1326 291584 : klmn=j0lmn+ilmn
1327 291584 : ro=pawtab(itypat)%rhoij0(klmn)
1328 291584 : if ((ilmn==jlmn).and.test_pawu.and.test_exexch.and.test_lnspinat) then
1329 33483 : ro=ro*ratio*(roshift+zratio)
1330 : else
1331 258101 : ro=ro*ratio*roshift
1332 : end if
1333 :
1334 291584 : klmn1=cpxocc*(klmn-1)+1
1335 291584 : if (abs(ro)>tol10) then
1336 16263 : pawrhoij(iatom_rhoij)%rhoijp(klmn1,ispden)=ro
1337 : else
1338 275321 : pawrhoij(iatom_rhoij)%rhoijp(klmn1,ispden)=zero
1339 : end if
1340 :
1341 329041 : if (ispden==nspden_rhoij) then
1342 447601 : if (any(abs(pawrhoij(iatom_rhoij)%rhoijp(klmn1,:))>tol10)) then
1343 12512 : nselect=nselect+1;nselect1=nselect1+cpxocc
1344 12512 : pawrhoij(iatom_rhoij)%rhoijselect(nselect)=klmn
1345 32227 : do jspden=1,nspden_rhoij
1346 32227 : pawrhoij(iatom_rhoij)%rhoijp(nselect1,jspden)=pawrhoij(iatom_rhoij)%rhoijp(klmn1,jspden)
1347 : end do
1348 : end if
1349 : end if
1350 :
1351 : end do
1352 : end do
1353 :
1354 : end do
1355 2260 : pawrhoij(iatom_rhoij)%nrhoijsel=nselect
1356 2260 : if (nselect<pawrhoij(iatom_rhoij)%lmn2_size) &
1357 165268 : & pawrhoij(iatom_rhoij)%rhoijselect(nselect+1:pawrhoij(iatom_rhoij)%lmn2_size)=0
1358 :
1359 : ! Non-collinear magnetism: avoid zero magnetization, because it produces numerical instabilities
1360 : ! Add a small real to the magnetization ; not yet activated => must be tested.
1361 : ! if (pawrhoij(iatom_rhoij)%nspden==4.and.spinat_zero) then
1362 : ! pawrhoij(iatom_rhoij)%rhoijp(:,4)=pawrhoij(iatom_rhoij)%rhoijp(:,4)+tol10
1363 : ! end if
1364 3477 : ABI_FREE(lnspinat)
1365 : end do ! iatom_rhoij
1366 :
1367 : !Destroy atom table used for parallelism
1368 1217 : call free_my_atmtab(my_atmtab,my_atmtab_allocated)
1369 :
1370 : DBG_EXIT("COLL")
1371 :
1372 1217 : end subroutine initrhoij
1373 : !!***
1374 :
1375 : !----------------------------------------------------------------------
1376 :
1377 : END MODULE m_paw_occupancies
1378 : !!***
|
