Branch data 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 : 34220 : subroutine pawmkrhoij(atindx,atindx1,cprj,dimcprj,istwfk,kptopt,mband,mband_cprj,mcprj,mkmem,mpi_enreg,&
103 [ - + ]: 17110 : & 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 : 17110 : integer, allocatable :: typat(:),req_correl(:,:,:)
133 : : real(dp) :: occup(2)
134 : 17110 : 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 : 17110 : type(pawcprj_type),allocatable :: cprj_tmp(:,:),cwaveprj(:,:),cwaveprjb(:,:)
137 : 17110 : type(pawcprj_type),pointer :: cprj_ptr(:,:)
138 : 17110 : 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 : 17110 : spaceComm=mpi_enreg%comm_kpt
150 : 17110 : me=mpi_enreg%me_kpt
151 : :
152 : : !Check size of cprj
153 : 17110 : my_nspinor=max(1,nspinor/mpi_enreg%nproc_spinor)
154 [ - + ]: 17110 : 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 : 17110 : nprocband=(mband/mband_cprj)
161 [ + + - + ]: 17110 : 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 [ - + - - ]: 17110 : 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 [ + + ]: 17110 : if(paw_dmft%use_sc_dmft/=0) then
174 : 117 : nbandc1=paw_dmft%mbandc
175 : : else
176 : : nbandc1=1
177 : : end if
178 : :
179 : : !Size of pawrhoij datastructure
180 : 17110 : nrhoij=size(pawrhoij)
181 : :
182 : : !Check if pawrhoij is distributed over atomic sites
183 [ - + - - ]: 17110 : 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 [ - + + - : 138672 : ABI_MALLOC(cwaveprj,(natom,nspinor))
- + + - +
- + - - +
- + + + +
+ ]
191 : 17110 : call pawcprj_alloc(cwaveprj,0,dimcprj)
192 [ + + ]: 17110 : if(paw_dmft%use_sc_dmft/=0) then
193 [ - + + - : 907 : ABI_MALLOC(cwaveprjb,(natom,nspinor))
- + + - -
+ - + + +
+ + ]
194 : 117 : call pawcprj_alloc(cwaveprjb,0,dimcprj)
195 : : end if
196 : :
197 [ + + + + ]: 17110 : 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 : 17110 : iorder_cprj=0
208 : :
209 : : !Build and initialize unpacked rhoij (to be computed here)
210 : 17110 : call pawrhoij_init_unpacked(pawrhoij)
211 : :
212 : : !If pawrhoij is MPI-distributed over atomic sites, we have to gather it
213 [ - + ]: 17110 : 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 : 17110 : pawrhoij_all => pawrhoij
231 : : end if
232 : :
233 : : !LOOP OVER SPINS
234 : 17110 : option=1
235 : 17110 : use_timerev=(kptopt>0.and.kptopt<3)
236 [ + + + + ]: 17110 : use_zeromag=(pawrhoij_all(1)%nspden==4.and.nspden==1)
237 : 17110 : bdtot_index=0;ibg=0;jbg=0
238 [ + + ]: 36176 : do isppol=1,nsppol
239 : :
240 : : ! LOOP OVER k POINTS
241 [ + + ]: 212516 : do ikpt=1,nkpt
242 : :
243 : 176340 : nband_k=nband(ikpt+(isppol-1)*nkpt)
244 : 176340 : nband_k_cprj=nband_k/nprocband
245 : 176340 : wtk_k=wtk(ikpt)
246 : :
247 [ + + ]: 176340 : if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me)) then
248 : 39425 : bdtot_index=bdtot_index+nband_k
249 : 39425 : cycle
250 : : end if
251 : :
252 [ + + ]: 136915 : cplex=2;if (istwfk(ikpt)>1) cplex=1
253 : :
254 : : ! In case of spinors parallelism, need some extra storage
255 [ + + ]: 136915 : if (mpi_enreg%paral_spinor==1) then
256 : 496 : nband_k_cprj_used=min(max_nband_cprj,nband_k_cprj)
257 [ - + + - : 30688 : ABI_MALLOC(cprj_tmp,(natom,my_nspinor*nband_k_cprj_used))
- + + - +
- + - - +
- + - + +
+ + + ]
258 [ - + + - : 59392 : ABI_MALLOC(cprj_ptr,(natom, nspinor*nband_k_cprj_used))
- + + - +
- + - - +
- + + + +
+ ]
259 : 496 : call pawcprj_alloc(cprj_tmp,0,dimcprj)
260 : 496 : call pawcprj_alloc(cprj_ptr,0,dimcprj)
261 : : else
262 : 136419 : 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 [ + + + + ]: 136915 : 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 : 136915 : ierr = 0
336 : : ! LOOP OVER BANDS
337 : 136915 : ib_this_proc=0;jb_this_proc=0
338 [ + + ]: 1702346 : do ib=1,nband_k
339 : 1565431 : iband=bdtot_index+ib
340 : :
341 : : ! Parallelization: treat only some bands
342 : : if(xmpi_paral==1)then
343 [ + + ]: 1565431 : if (paral_kgb==1) then
344 [ + + ]: 404484 : if (mod((ib-1)/mpi_enreg%bandpp,mpi_enreg%nproc_band)/=mpi_enreg%me_band) cycle
345 : : else
346 [ - + ]: 1160947 : if (mpi_enreg%proc_distrb(ikpt,ib,isppol)/=me) cycle
347 : : end if
348 : : end if
349 : 1324479 : 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 [ + + ]: 1324479 : if (mpi_enreg%paral_spinor==1) then
354 : 9568 : jb_this_proc=jb_this_proc+1
355 [ + + ]: 9568 : if (mod(jb_this_proc,nband_k_cprj_used)==1) then
356 : 496 : ib_this_proc=1
357 : 496 : nband_k_cprj_read=nband_k_cprj_used
358 [ - + ]: 496 : 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 : 496 : & 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 : 496 : & mpi_enreg%comm_spinor,ierr)
364 : : end if
365 : : end if
366 : :
367 : : ! DMFT: LOOP ON ADDITIONAL BANDS
368 [ + + ]: 3134397 : 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 [ + + ]: 1673003 : if(ibc1 /= 1) then
374 [ + - ]: 348524 : if (paw_dmft%use_sc_dmft == 0) then
375 : : cycle
376 : : else
377 [ + + ]: 348524 : 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 [ + + ]: 1581831 : if(paw_dmft%use_sc_dmft /= 0) then
383 [ + + ]: 273492 : 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 : 266248 : ib1 = paw_dmft%include_bands(ibc1) ! indice reel de la bande
388 : :
389 : 266248 : use_nondiag_occup_dmft = 1
390 : 266248 : locc_test = abs(paw_dmft%occnd(1,ib,ib1,ikpt,isppol))+abs(paw_dmft%occnd(2,ib,ib1,ikpt,isppol))>tol8
391 : :
392 : 266248 : occup(1) = paw_dmft%occnd(1,ib,ib1,ikpt,isppol)
393 : 266248 : 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 [ + + - + ]: 266248 : if (locc_test .or. mkmem == 0) then
398 : :
399 [ + + ]: 122869 : 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 [ + - ]: 121361 : & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
411 : : end if
412 : : end if
413 : : else ! nbandc1=1
414 : 7244 : use_nondiag_occup_dmft=0
415 : 7244 : locc_test = (abs(occ(iband))>tol8)
416 : 7244 : occup(1) = occ(iband)
417 : : end if
418 : : else ! nbandc1=1
419 : 1308339 : use_nondiag_occup_dmft=0
420 : 1308339 : locc_test = (abs(occ(iband))>tol8)
421 : 1308339 : 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 [ + + - + ]: 1581831 : 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 [ - + ]: 966936 : & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
430 : : end if
431 : :
432 : : ! Accumulate contribution from (occupied) current band
433 [ + + ]: 3147262 : if (locc_test) then
434 [ + + ]: 966936 : 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 [ - + ]: 122869 : & occ_k_2=occup(2))
438 : : else
439 : : call pawaccrhoij(atindx,cplex,cwaveprj,cwaveprj ,0,isppol,nrhoij,natom,&
440 [ - + ]: 844067 : & 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 [ + + ]: 136915 : if (mpi_enreg%paral_spinor==1) then
447 : 496 : call pawcprj_free(cprj_tmp)
448 : 496 : call pawcprj_free(cprj_ptr)
449 [ + - + + : 19632 : ABI_FREE(cprj_tmp)
- + - + -
+ ]
450 [ - + ]: 496 : ABI_FREE(cprj_ptr)
451 : : else
452 : 136419 : nullify(cprj_ptr)
453 : : end if
454 : :
455 : 136915 : bdtot_index=bdtot_index+nband_k
456 [ - + ]: 155981 : if (mkmem/=0) then
457 [ + + ]: 136915 : if (mpi_enreg%paral_spinor==0) then
458 : 136419 : ibg=ibg+ nspinor*nband_k_cprj
459 : : else
460 : 496 : 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 : 17110 : call pawcprj_free(cwaveprj)
471 [ + - + + : 67847 : ABI_FREE(cwaveprj)
- + - + -
+ ]
472 : :
473 [ + + ]: 17110 : if(paw_dmft%use_sc_dmft/=0) then
474 : 117 : call pawcprj_free(cwaveprjb)
475 [ + - + + : 551 : ABI_FREE(cwaveprjb)
- + - + -
+ ]
476 : : end if
477 : :
478 [ + + ]: 17110 : 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 [ + + + + ]: 17110 : if (paral_kgb==1.and.nprocband>1) then
485 : 5232 : call pawrhoij_mpisum_unpacked(pawrhoij_all,spaceComm,comm2=mpi_enreg%comm_band)
486 : : else
487 : 11878 : call pawrhoij_mpisum_unpacked(pawrhoij_all,spaceComm)
488 : : end if
489 : :
490 : : !In case of distribution over atomic sites, dispatch rhoij
491 [ - + ]: 17110 : 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 [ - + - + : 34220 : 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 : 3236282 : subroutine pawaccrhoij(atindx,cplex,cwaveprj,cwaveprj1,ipert,isppol,my_natom,natom,&
567 : 1618141 : & nspinor,occ_k,option,pawrhoij,use_timerev,use_zeromag,wtk_k,occ_k_2, &
568 [ + + - + : 1618141 : & 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 : 1618141 : integer,pointer :: my_atmtab(:)
593 : 3236282 : real(dp) :: cpi0(2,nspinor),cpi1(2,nspinor),cpj0(2,nspinor),cpj1(2,nspinor)
594 : 1618141 : real(dp) :: dcpi0(2,nspinor,9),dcpj0(2,nspinor,9)
595 : :
596 : : ! ***********************************************************************
597 : :
598 : : DBG_ENTER("COLL")
599 : :
600 : 1618141 : if (my_natom==0) return
601 : :
602 : 1618141 : ncpgr=0
603 [ + + + + : 1618141 : if (option==2.and.(ipert<=natom.or.ipert==natom+3.or.ipert==natom+4)) ncpgr=1
+ + + + ]
604 [ + + ]: 1618141 : if (option==3) ncpgr=cwaveprj(1,1)%ncpgr
605 : :
606 : : !Tests
607 [ + + + - : 1618141 : 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 [ + + - + ]: 1618141 : 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 [ + + - + ]: 1618141 : 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 [ + + - + ]: 1618141 : 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 [ + + + - ]: 1618141 : paral_atom=(present(comm_atom).and.(my_natom/=natom))
626 [ - + ]: 1618141 : nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
627 [ + + ]: 1618141 : 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 : 1618141 : & my_natom_ref=my_natom)
630 : :
631 : 1618141 : weight=wtk_k*occ_k
632 : 1618141 : weight_2=zero
633 [ + + ]: 1618141 : if(present(occ_k_2)) then
634 [ - + ]: 122869 : 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 : 122869 : weight_2=wtk_k*occ_k_2
640 : : end if
641 [ + + + + : 1618141 : if (pawrhoij(1)%nspden==2.and.pawrhoij(1)%nsppol==1.and.nspinor==1) weight=half*weight
+ - ]
642 [ + + + + : 1618141 : 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 [ + + ]: 1618141 : if (option==1) then
645 : :
646 : : ! ==================================================================
647 : : ! === OPTION 1: Accumulate (n,k) contribution to rhoij =============
648 : : ! ==================================================================
649 [ + + + + ]: 967544 : compute_impart=((.not.use_timerev).and.(pawrhoij(1)%cplex_rhoij==2))
650 [ - + ]: 263997 : compute_impart_cplex=((compute_impart).and.(cplex==2))
651 [ + + ]: 967544 : if (nspinor==1) then
652 : 727402 : cplex_rhoij=pawrhoij(1)%cplex_rhoij
653 [ + + ]: 727402 : if (cplex_rhoij==1) then
654 [ + + ]: 2659178 : do iatom=1,my_natom
655 [ - + ]: 1982058 : iatom1=iatom;if (paral_atom) iatom1=my_atmtab(iatom)
656 : 1982058 : iatm=atindx(iatom1)
657 [ + + ]: 26258493 : do jlmn=1,pawrhoij(iatom)%lmn_size
658 : 23599315 : j0lmn=jlmn*(jlmn-1)/2
659 [ + + ]: 66896199 : cpj0(1:cplex,1)=cwaveprj(iatm,1)%cp(1:cplex,jlmn)
660 [ + + ]: 207909559 : do ilmn=1,jlmn
661 : 182328186 : klmn=j0lmn+ilmn
662 [ + + ]: 527187267 : cpi0(1:cplex,1)=cwaveprj1(iatm,1)%cp(1:cplex,ilmn)
663 : : ro11_re=zero
664 [ + + ]: 527187267 : do iplex=1,cplex
665 : 527187267 : ro11_re=ro11_re+cpi0(iplex,1)*cpj0(iplex,1)
666 : : end do
667 : 182328186 : pawrhoij(iatom)%rhoij_(klmn,isppol)=pawrhoij(iatom)%rhoij_(klmn,isppol)+weight*ro11_re
668 [ + + ]: 205927501 : 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 [ + + ]: 138196 : do iatom=1,my_natom
677 [ - + ]: 87914 : iatom1=iatom;if (paral_atom) iatom1=my_atmtab(iatom)
678 : 87914 : iatm=atindx(iatom1)
679 [ + + ]: 841508 : do jlmn=1,pawrhoij(iatom)%lmn_size
680 : 703312 : j0lmn=jlmn*(jlmn-1)/2
681 [ + + ]: 2109936 : cpj0(1:cplex,1)=cwaveprj(iatm,1)%cp(1:cplex,jlmn)
682 [ + + ]: 3956130 : do ilmn=1,jlmn
683 : 3164904 : klmn=j0lmn+ilmn
684 : 3164904 : klmn_re=cplex_rhoij*(klmn-1)+1
685 [ + + ]: 9494712 : cpi0(1:cplex,1)=cwaveprj1(iatm,1)%cp(1:cplex,ilmn)
686 : : ro11_re=zero
687 [ + + ]: 9494712 : do iplex=1,cplex
688 : 9494712 : ro11_re=ro11_re+cpi0(iplex,1)*cpj0(iplex,1)
689 : : end do
690 : 3164904 : pawrhoij(iatom)%rhoij_(klmn_re,isppol)=pawrhoij(iatom)%rhoij_(klmn_re,isppol)+weight*ro11_re
691 [ - + ]: 3164904 : 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 [ + - ]: 3868216 : if (compute_impart_cplex) then
696 : 3164904 : klmn_im=klmn_re+1
697 : 3164904 : ro11_im=cpi0(1,1)*cpj0(2,1)-cpi0(2,1)*cpj0(1,1)
698 : 3164904 : pawrhoij(iatom)%rhoij_(klmn_im,isppol)=pawrhoij(iatom)%rhoij_(klmn_im,isppol)+weight*ro11_im
699 [ - + ]: 3164904 : 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 [ + + ]: 883615 : do iatom=1,my_natom
709 [ - + ]: 643473 : iatom1=iatom;if (paral_atom) iatom1=my_atmtab(iatom)
710 : 643473 : iatm=atindx(iatom1)
711 : 643473 : cplex_rhoij=pawrhoij(iatom)%cplex_rhoij
712 : 643473 : nspden_rhoij=pawrhoij(iatom)%nspden
713 [ + + ]: 15500819 : do jlmn=1,pawrhoij(iatom)%lmn_size
714 : 14617204 : j0lmn=jlmn*(jlmn-1)/2
715 [ + + ]: 43851612 : cpj0(1:cplex,1)=cwaveprj(iatm,1)%cp(1:cplex,jlmn)
716 [ + + ]: 43851612 : cpj0(1:cplex,2)=cwaveprj(iatm,2)%cp(1:cplex,jlmn)
717 [ + + ]: 213997634 : do ilmn=1,jlmn
718 : 198736957 : klmn=j0lmn+ilmn
719 : 198736957 : klmn_re=cplex_rhoij*(klmn-1)+1
720 [ + + ]: 596210871 : cpi0(1:cplex,1)=cwaveprj1(iatm,1)%cp(1:cplex,ilmn)
721 [ + + ]: 596210871 : cpi0(1:cplex,2)=cwaveprj1(iatm,2)%cp(1:cplex,ilmn)
722 : : ro11_re=zero;ro22_re=zero
723 : 596210871 : ro12_re=zero;ro21_re=zero
724 : 596210871 : ro12_im=zero;ro21_im=zero
725 [ + + ]: 596210871 : do iplex=1,cplex
726 : 397473914 : ro11_re=ro11_re+cpi0(iplex,1)*cpj0(iplex,1)
727 : 596210871 : ro22_re=ro22_re+cpi0(iplex,2)*cpj0(iplex,2)
728 : : end do
729 : 198736957 : pawrhoij(iatom)%rhoij_(klmn_re,1)=pawrhoij(iatom)%rhoij_(klmn_re,1)+weight*(ro11_re+ro22_re)
730 [ + + ]: 198736957 : if (nspden_rhoij>1) then
731 [ + + ]: 594856551 : do iplex=1,cplex
732 : 396571034 : ro12_re=ro12_re+cpi0(iplex,2)*cpj0(iplex,1)
733 : 594856551 : ro21_re=ro21_re+cpi0(iplex,1)*cpj0(iplex,2)
734 : : end do
735 [ + + ]: 198285517 : if (.not.use_zeromag) then
736 : 190245473 : pawrhoij(iatom)%rhoij_(klmn_re,4)=pawrhoij(iatom)%rhoij_(klmn_re,4)+weight*(ro11_re-ro22_re)
737 : 190245473 : pawrhoij(iatom)%rhoij_(klmn_re,2)=pawrhoij(iatom)%rhoij_(klmn_re,2)+weight*(ro12_re+ro21_re)
738 : : end if
739 [ + - ]: 198285517 : 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 : 198285517 : ro12_im=cpi0(1,2)*cpj0(2,1)-cpi0(2,2)*cpj0(1,1)
743 : 198285517 : ro21_im=cpi0(1,1)*cpj0(2,2)-cpi0(2,1)*cpj0(1,2)
744 [ + + ]: 198285517 : if (.not.use_zeromag) then
745 : 190245473 : 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 [ + + ]: 198736957 : if (present(occ_k_2)) then
750 : 165972750 : ro11_im=cpi0(1,1)*cpj0(2,1)-cpi0(2,1)*cpj0(1,1)
751 : 165972750 : ro22_im=cpi0(1,2)*cpj0(2,2)-cpi0(2,2)*cpj0(1,2)
752 : 165972750 : pawrhoij(iatom)%rhoij_(klmn_re,1)=pawrhoij(iatom)%rhoij_(klmn_re,1)+weight_2*(-ro11_im-ro22_im)
753 : 165972750 : pawrhoij(iatom)%rhoij_(klmn_re,2)=pawrhoij(iatom)%rhoij_(klmn_re,2)+weight_2*(-ro21_im-ro12_im)
754 : 165972750 : pawrhoij(iatom)%rhoij_(klmn_re,3)=pawrhoij(iatom)%rhoij_(klmn_re,3)+weight_2*(-ro12_re+ro21_re)
755 : 165972750 : pawrhoij(iatom)%rhoij_(klmn_re,4)=pawrhoij(iatom)%rhoij_(klmn_re,4)+weight_2*(-ro11_im+ro22_im)
756 : : end if
757 [ + + ]: 213354161 : if (compute_impart) then
758 : 192167650 : klmn_im=klmn_re+1
759 [ + - ]: 192167650 : if (nspden_rhoij>1) pawrhoij(iatom)%rhoij_(klmn_im,3)=pawrhoij(iatom)%rhoij_(klmn_im,3)+weight*(ro12_re-ro21_re)
760 [ + - ]: 192167650 : if (cplex==2) then
761 : 192167650 : ro11_im=cpi0(1,1)*cpj0(2,1)-cpi0(2,1)*cpj0(1,1)
762 : 192167650 : ro22_im=cpi0(1,2)*cpj0(2,2)-cpi0(2,2)*cpj0(1,2)
763 : 192167650 : pawrhoij(iatom)%rhoij_(klmn_im,1)=pawrhoij(iatom)%rhoij_(klmn_im,1)+weight*(ro11_im+ro22_im)
764 [ + - ]: 192167650 : if (nspden_rhoij>1) then
765 : 192167650 : pawrhoij(iatom)%rhoij_(klmn_im,4)=pawrhoij(iatom)%rhoij_(klmn_im,4)+weight*(ro11_im-ro22_im)
766 : 192167650 : pawrhoij(iatom)%rhoij_(klmn_im,2)=pawrhoij(iatom)%rhoij_(klmn_im,2)+weight*(ro12_im+ro21_im)
767 : : end if
768 [ + + ]: 192167650 : if (present(occ_k_2)) then
769 : 165972750 : pawrhoij(iatom)%rhoij_(klmn_im,1)=pawrhoij(iatom)%rhoij_(klmn_im,1)+weight_2*( ro11_re+ro22_re)
770 : 165972750 : pawrhoij(iatom)%rhoij_(klmn_im,2)=pawrhoij(iatom)%rhoij_(klmn_im,2)+weight_2*( ro21_re+ro12_re)
771 : 165972750 : pawrhoij(iatom)%rhoij_(klmn_im,3)=pawrhoij(iatom)%rhoij_(klmn_im,3)+weight_2*(-ro12_im+ro21_im)
772 : 165972750 : 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 [ + + ]: 650597 : 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 : 635829 : compute_impart=(pawrhoij(1)%qphase==2)
790 [ + + - + ]: 635829 : compute_impart_cplex=(compute_impart.and.(cplex==2))
791 [ + - ]: 635829 : if (nspinor==1) then
792 [ + + ]: 2124324 : do iatom=1,my_natom
793 [ - + ]: 1488495 : iatom1=iatom;if (paral_atom) iatom1=my_atmtab(iatom)
794 : 1488495 : iatm=atindx(iatom1)
795 : 1488495 : cplex_rhoij=pawrhoij(iatom)%cplex_rhoij
796 : 1488495 : iq0=cplex_rhoij*pawrhoij(iatom)%lmn2_size
797 [ + + ]: 14557324 : do jlmn=1,pawrhoij(iatom)%lmn_size
798 : 12433000 : j0lmn=jlmn*(jlmn-1)/2
799 [ + + ]: 37299000 : cpj0(1:2,1)=cwaveprj (iatm,1)%cp(1:2,jlmn)
800 [ + + ]: 37299000 : cpj1(1:2,1)=cwaveprj1(iatm,1)%cp(1:2,jlmn)
801 [ + + ]: 74314731 : do ilmn=1,jlmn
802 : 60393236 : klmn=j0lmn+ilmn
803 : 60393236 : klmn_re=cplex_rhoij*(klmn-1)+1
804 [ + + ]: 181179708 : cpi0(1:2,1)=cwaveprj (iatm,1)%cp(1:2,ilmn)
805 [ + + ]: 181179708 : cpi1(1:2,1)=cwaveprj1(iatm,1)%cp(1:2,ilmn)
806 : 60393236 : ro11_re=zero
807 [ + + ]: 181179708 : do iplex=1,cplex
808 : 181179708 : ro11_re=ro11_re+cpi0(iplex,1)*cpj1(iplex,1)+cpj0(iplex,1)*cpi1(iplex,1)
809 : : end do
810 : 60393236 : pawrhoij(iatom)%rhoij_(klmn_re,isppol)=pawrhoij(iatom)%rhoij_(klmn_re,isppol)+weight*ro11_re
811 [ + + ]: 72826236 : 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 [ + + ]: 635829 : if (ipert/=natom+2) then
887 : 598165 : compute_impart=(pawrhoij(1)%cplex_rhoij==2)
888 [ - + - - ]: 598165 : compute_impart_cplex=(compute_impart.and.(cplex==2))
889 : 598165 : substract_diagonal=(ipert==natom+3)
890 [ + - ]: 598165 : if (nspinor==1) then
891 [ + + ]: 1999556 : do iatom=1,my_natom
892 [ - + ]: 1401391 : iatom1=iatom;if (paral_atom) iatom1=my_atmtab(iatom)
893 : 1401391 : iatm=atindx(iatom1)
894 [ + + + + ]: 1401391 : if (ipert<=natom.and.iatom/=ipert) cycle
895 : 707649 : cplex_rhoij=pawrhoij(iatom)%cplex_rhoij
896 [ + + ]: 7293246 : do jlmn=1,pawrhoij(iatom)%lmn_size
897 : 5987432 : j0lmn=jlmn*(jlmn-1)/2
898 [ + + ]: 17962296 : cpj0 (1:2,1) =cwaveprj(iatm,1)%cp (1:2 ,jlmn)
899 [ + + ]: 17962296 : dcpj0(1:2,1,1)=cwaveprj(iatm,1)%dcp(1:2,1,jlmn)
900 [ + + ]: 37243315 : do ilmn=1,jlmn
901 : 29854492 : klmn=j0lmn+ilmn
902 : 29854492 : klmn_re=cplex_rhoij*(klmn-1)+1
903 [ + + ]: 89563476 : cpi0 (1:2,1) =cwaveprj(iatm,1)%cp (1:2 ,ilmn)
904 [ + + ]: 89563476 : dcpi0(1:2,1,1)=cwaveprj(iatm,1)%dcp(1:2,1,ilmn)
905 : 29854492 : ro11_re=zero
906 [ + + ]: 89563476 : do iplex=1,cplex
907 : 89563476 : ro11_re=ro11_re+dcpi0(iplex,1,1)*cpj0(iplex,1)+cpi0(iplex,1)*dcpj0(iplex,1,1)
908 : : end do
909 [ + + ]: 29854492 : 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 : 29854492 : 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 [ - + ]: 35841924 : 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 [ + - ]: 14768 : else if (option==3) then
1015 : :
1016 : : ! ==================================================================
1017 : : ! === OPTION 3: Accumulate (n,k) contribution to drhoij/dr =========
1018 : : ! ==================================================================
1019 : :
1020 [ + + - + ]: 14768 : compute_impart=((.not.use_timerev).and.(pawrhoij(1)%cplex_rhoij==2))
1021 [ # # ]: 0 : compute_impart_cplex=((compute_impart).and.(cplex==2))
1022 [ + - ]: 14768 : if (nspinor==1) then
1023 [ + + ]: 42936 : do iatom=1,my_natom
1024 [ - + ]: 28168 : iatom1=iatom;if (paral_atom) iatom1=my_atmtab(iatom)
1025 : 28168 : iatm=atindx(iatom1)
1026 : 28168 : cplex_rhoij=pawrhoij(iatom)%cplex_rhoij
1027 [ + + ]: 313720 : do jlmn=1,pawrhoij(iatom)%lmn_size
1028 : 270784 : j0lmn=jlmn*(jlmn-1)/2
1029 [ + + ]: 812352 : cpj0(1:cplex,1) =cwaveprj(iatm,1)%cp (1:cplex,jlmn)
1030 [ + + + + ]: 5673088 : dcpj0(1:cplex,1,1:ncpgr)=cwaveprj(iatm,1)%dcp(1:cplex,1:ncpgr,jlmn)
1031 [ + + ]: 1925416 : do ilmn=1,jlmn
1032 : 1626464 : klmn=j0lmn+ilmn
1033 : 1626464 : klmn_re=cplex_rhoij*(klmn-1)+1
1034 [ + + ]: 4879392 : cpi0(1:cplex,1) =cwaveprj(iatm,1)%cp (1:cplex,ilmn)
1035 [ + + + + ]: 33684032 : dcpi0(1:cplex,1,1:ncpgr)=cwaveprj(iatm,1)%dcp(1:cplex,1:ncpgr,ilmn)
1036 [ + + ]: 12312320 : do mu=1,ncpgr
1037 : : ro11_re=zero
1038 [ + + ]: 32057568 : do iplex=1,cplex
1039 : 32057568 : ro11_re=ro11_re+dcpi0(iplex,1,mu)*cpj0(iplex,1)+cpi0(iplex,1)*dcpj0(iplex,1,mu)
1040 : : end do
1041 : 12312320 : pawrhoij(iatom)%grhoij(mu,klmn_re,isppol)=pawrhoij(iatom)%grhoij(mu,klmn_re,isppol)+weight*ro11_re
1042 : : end do
1043 [ - + ]: 1897248 : 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 : 1618141 : call free_my_atmtab(my_atmtab,my_atmtab_allocated)
1126 : :
1127 : : DBG_EXIT("COLL")
1128 : :
1129 [ - + ]: 1618141 : 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 : 1251 : subroutine initrhoij(cpxocc,lexexch,lpawu,my_natom,natom,nspden,nspinor,nsppol,&
1174 : 1251 : & ntypat,pawrhoij,pawspnorb,pawtab,qphase,spinat,typat,&
1175 : : & ngrhoij,nlmnmix,use_rhoij_,use_rhoijres,& ! optional arguments
1176 [ + - + + : 1251 : & 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 : 1251 : integer,pointer :: my_atmtab(:),lnspinat(:)
1202 : 1251 : real(dp),allocatable :: occ(:)
1203 : :
1204 : : !************************************************************************
1205 : :
1206 : : DBG_ENTER("COLL")
1207 : :
1208 : : !PAW+U and local exact-exchange restriction
1209 [ + + ]: 3136 : do itypat=1,ntypat
1210 [ + + + + : 3136 : 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 [ + - + + ]: 1251 : paral_atom=(present(comm_atom).and.(my_natom/=natom))
1218 [ + + ]: 1251 : nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
1219 [ + - ]: 1251 : my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
1220 : 1251 : 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 : 1251 : & nspden=nspden,spnorb=pawspnorb,cpxocc=cpxocc,cplex=qphase)
1224 : :
1225 : 1251 : ratio=one;if (nspden_rhoij==2) ratio=half
1226 [ + + + + : 13482 : 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 [ + + ]: 1251 : if (my_natom>0) then
1236 [ + + ]: 1109 : ngrhoij0=0;if (present(ngrhoij)) ngrhoij0=ngrhoij
1237 [ + + ]: 1109 : nlmnmix0=0;if (present(nlmnmix)) nlmnmix0=nlmnmix
1238 [ + + ]: 1109 : use_rhoij_0=0;if (present(use_rhoij_)) use_rhoij_0=use_rhoij_
1239 [ + + ]: 1109 : use_rhoijres0=0;if (present(use_rhoijres)) use_rhoijres0=use_rhoijres
1240 [ + + ]: 1109 : 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 : 789 : & pawtab=pawtab,ngrhoij=ngrhoij0,nlmnmix=nlmnmix0,use_rhoij_=use_rhoij_0,use_rhoijres=use_rhoijres0)
1247 : : end if
1248 : : end if
1249 : :
1250 [ + + ]: 3583 : do iatom_rhoij=1,my_natom
1251 [ + + ]: 2332 : iatom=iatom_rhoij;if (paral_atom) iatom=my_atmtab(iatom_rhoij)
1252 : 2332 : itypat=typat(iatom)
1253 : 2332 : nselect=0
1254 [ - + + - : 6996 : ABI_MALLOC(lnspinat,(pawtab(itypat)%basis_size))
+ - + - -
+ - + ]
1255 [ + + ]: 12510 : lnspinat=-1
1256 : : ! Determine occupancies of each orbital
1257 [ + + ]: 2332 : if (nspden_rhoij==2) then
1258 [ - + + - : 1257 : ABI_MALLOC(occ,(pawtab(itypat)%basis_size))
+ - + - -
+ + - ]
1259 [ + + ]: 2620 : occ=zero
1260 [ + + ]: 6366 : do jlmn=1,pawtab(itypat)%lmn_size
1261 : 5947 : ln=pawtab(itypat)%indlmn(5,jlmn)
1262 : 5947 : klmn=jlmn*(jlmn+1)/2
1263 : 6366 : occ(ln)=occ(ln)+pawtab(itypat)%rhoij0(klmn)
1264 : : end do
1265 [ + + ]: 2620 : do ln=1,pawtab(itypat)%basis_size
1266 [ + + + + ]: 2201 : if(pawtab(itypat)%orbitals(ln)==0.and.occ(ln)==1) lnspinat(ln)=ln
1267 [ + + + + : 2201 : if(pawtab(itypat)%orbitals(ln)==1.and.(occ(ln)>=1.and.occ(ln)<=5)) lnspinat(ln)=ln
+ + ]
1268 [ + + + + : 2201 : if(pawtab(itypat)%orbitals(ln)==2.and.(occ(ln)>=1.and.occ(ln)<=9)) lnspinat(ln)=ln
+ + ]
1269 [ + + + + : 2620 : if(pawtab(itypat)%orbitals(ln)==3.and.(occ(ln)>=1.and.occ(ln)<=13)) lnspinat(ln)=ln
+ - ]
1270 : : end do
1271 [ - + ]: 419 : ABI_FREE(occ)
1272 : : end if
1273 : : lnspinat0=maxval(lnspinat)
1274 : 2332 : lnspinat0=-1
1275 : :
1276 : : ! Determine Z (trace of rhoij0 or part of it)
1277 : 2332 : zz=zero
1278 [ + + ]: 26892 : do jlmn=1,pawtab(itypat)%lmn_size
1279 : 24560 : jl=pawtab(itypat)%indlmn(1,jlmn)
1280 : 24560 : ln=pawtab(itypat)%indlmn(5,jlmn)
1281 : 24560 : j0lmn=jlmn*(jlmn-1)/2
1282 : 24560 : test_lnspinat=(lnspinat0==-1.or.lnspinat(ln)==ln)
1283 [ + + + + ]: 24560 : test_pawu=(lpawu(itypat)==-1.or.lpawu(itypat)==jl)
1284 [ + + + + ]: 24560 : test_exexch=(lexexch(itypat)==-1.or.lexexch(itypat)==jl)
1285 [ + + ]: 212015 : do ilmn=1,jlmn
1286 : 185123 : klmn=j0lmn+ilmn
1287 [ + + + + ]: 185123 : if ((ilmn==jlmn).and.test_pawu.and.test_exexch.and.test_lnspinat) &
1288 : 47008 : & zz=zz+pawtab(itypat)%rhoij0(klmn)
1289 : : end do
1290 : : end do
1291 : :
1292 : : ! Compute rhoij from tabulated value and magnetization
1293 [ + + ]: 5734 : do ispden=1,nspden_rhoij
1294 : :
1295 : 3402 : zratio=zero
1296 : 3402 : roshift=one
1297 : 3402 : ratio=one
1298 [ + + ]: 3402 : if (nspden_rhoij==2) then
1299 : 838 : ratio=half
1300 [ + + + + : 838 : if ((spinat(3,iatom)>tol12.and.ispden==1).or.&
+ + + + ]
1301 : : & (spinat(3,iatom)<tol12.and.ispden==2)) then
1302 [ + + ]: 419 : if(abs(zz)>tol12)then
1303 : 418 : zratio=two*abs(spinat(3,iatom))/zz
1304 : : else
1305 : : zratio=zero
1306 : : end if
1307 : : end if
1308 [ + + + + ]: 2564 : else if (nspden_rhoij==4.and.ispden>=2) then
1309 : 651 : roshift=zero
1310 [ + - + + ]: 651 : 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 : 3402 : nselect=0;nselect1=1-cpxocc
1318 [ + + ]: 45625 : do jlmn=1,pawtab(itypat)%lmn_size
1319 : 39891 : jl=pawtab(itypat)%indlmn(1,jlmn)
1320 : 39891 : ln=pawtab(itypat)%indlmn(5,jlmn)
1321 : 39891 : j0lmn=jlmn*(jlmn-1)/2
1322 : 39891 : test_lnspinat=(lnspinat0==-1.or.lnspinat(ln)==ln)
1323 [ + + + + ]: 39891 : test_pawu=(lpawu(itypat)==-1.or.lpawu(itypat)==jl)
1324 [ + + + + ]: 39891 : test_exexch=(lexexch(itypat)==-1.or.lexexch(itypat)==jl)
1325 [ + + ]: 366052 : do ilmn=1,jlmn
1326 : 322759 : klmn=j0lmn+ilmn
1327 : 322759 : ro=pawtab(itypat)%rhoij0(klmn)
1328 [ + + + + ]: 322759 : if ((ilmn==jlmn).and.test_pawu.and.test_exexch.and.test_lnspinat) then
1329 : 34989 : ro=ro*ratio*(roshift+zratio)
1330 : : else
1331 : 287770 : ro=ro*ratio*roshift
1332 : : end if
1333 : :
1334 : 322759 : klmn1=cpxocc*(klmn-1)+1
1335 [ + + ]: 322759 : if (abs(ro)>tol10) then
1336 : 16825 : pawrhoij(iatom_rhoij)%rhoijp(klmn1,ispden)=ro
1337 : : else
1338 : 305934 : pawrhoij(iatom_rhoij)%rhoijp(klmn1,ispden)=zero
1339 : : end if
1340 : :
1341 [ + + ]: 362650 : if (ispden==nspden_rhoij) then
1342 [ + + + + : 486813 : if (any(abs(pawrhoij(iatom_rhoij)%rhoijp(klmn1,:))>tol10)) then
+ + ]
1343 : 13018 : nselect=nselect+1;nselect1=nselect1+cpxocc
1344 : 13018 : pawrhoij(iatom_rhoij)%rhoijselect(nselect)=klmn
1345 [ + + ]: 34087 : do jspden=1,nspden_rhoij
1346 : 34087 : 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 : 2332 : pawrhoij(iatom_rhoij)%nrhoijsel=nselect
1356 [ + + ]: 2332 : if (nselect<pawrhoij(iatom_rhoij)%lmn2_size) &
1357 [ + + ]: 174225 : & 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 [ - + ]: 3583 : ABI_FREE(lnspinat)
1365 : : end do ! iatom_rhoij
1366 : :
1367 : : !Destroy atom table used for parallelism
1368 : 1251 : call free_my_atmtab(my_atmtab,my_atmtab_allocated)
1369 : :
1370 : : DBG_EXIT("COLL")
1371 : :
1372 : 1251 : end subroutine initrhoij
1373 : : !!***
1374 : :
1375 : : !----------------------------------------------------------------------
1376 : :
1377 : : END MODULE m_paw_occupancies
1378 : : !!***
|
