Branch data Line data Source code
1 : : !!****m* ABINIT/m_newrho
2 : : !! NAME
3 : : !! m_newrho
4 : : !!
5 : : !! FUNCTION
6 : : !!
7 : : !!
8 : : !! COPYRIGHT
9 : : !! Copyright (C) 2005-2026 ABINIT group (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_newrho
23 : :
24 : : use defs_basis
25 : : use defs_wvltypes
26 : : use m_errors
27 : : use m_abicore
28 : : use m_abi_mixing
29 : : use m_abi2big
30 : : use m_dtset
31 : : use m_precon
32 : :
33 : : use m_rcpaw, only : rcpaw_type
34 : : use m_extfpmd, only : extfpmd_type
35 : : use defs_datatypes, only : pseudopotential_type
36 : : use defs_abitypes, only : MPI_type
37 : : use m_time, only : timab
38 : : use m_geometry, only : metric
39 : : use m_pawtab, only : pawtab_type
40 : : use m_pawrhoij, only : pawrhoij_type,pawrhoij_filter
41 : : use m_prcref, only : prcref
42 : : use m_wvl_rho, only : wvl_prcref
43 : : use m_fft, only : fourdp
44 : :
45 : : implicit none
46 : :
47 : : private
48 : : !!***
49 : :
50 : : public :: newrho
51 : : !!***
52 : :
53 : : contains
54 : : !!***
55 : :
56 : : !!****f* ABINIT/newrho
57 : : !! NAME
58 : : !! newrho
59 : : !!
60 : : !! FUNCTION
61 : : !! Compute new trial density by mixing new and old values.
62 : : !! Call prcref to compute preconditioned residual density and forces,
63 : : !! Then, call one of the self-consistency drivers, then update density.
64 : : !!
65 : : !! INPUTS
66 : : !! atindx(natom)=index table for atoms (see gstate.f)
67 : : !! dielar(7)=input parameters for dielectric matrix:
68 : : !! diecut,dielng,diemac,diemix,diegap,dielam,diemixmag.
69 : : !! dielinv(2,npwdiel,nspden,npwdiel,nspden)=
70 : : !! inverse of the dielectric matrix in rec. space
71 : : !! dielstrt=number of the step at which the dielectric preconditioning begins.
72 : : !! dtset <type(dataset_type)>=all input variables in this dataset
73 : : !! | densfor_pred= governs the preconditioning of the atomic charges
74 : : !! | iprcel= governs the preconditioning of the density residual
75 : : !! | iprcfc= governs the preconditioning of the forces
76 : : !! | iscf=( <= 0 =>non-SCF), >0 => SCF)
77 : : !! | iscf =11 => determination of the largest eigenvalue of the SCF cycle
78 : : !! | iscf =12 => SCF cycle, simple mixing
79 : : !! | iscf =13 => SCF cycle, Anderson mixing
80 : : !! | iscf =14 => SCF cycle, Anderson mixing (order 2)
81 : : !! | iscf =15 => SCF cycle, CG based on the minimization of the energy
82 : : !! | iscf =17 => SCF cycle, Pulay mixing
83 : : !! | isecur=level of security of the computation
84 : : !! | mffmem=governs the number of FFT arrays which are fit in core memory
85 : : !! | it is either 1, in which case the array f_fftgr is used,
86 : : !! | or 0, in which case the array f_fftgr_disk is used
87 : : !! | natom=number of atoms
88 : : !! | nspden=number of spin-density components
89 : : !! | pawoptmix=-PAW- 1 if the computed residuals include the PAW (rhoij) part
90 : : !! | prtvol=control print volume and debugging
91 : : !! etotal=the total energy obtained from the input density
92 : : !! extfpmd <type(extfpmd_type)>=extended first-principles molecular dynamics type
93 : : !! fnametmp_fft=name of _FFT file
94 : : !! fcart(3,natom)=cartesian forces (hartree/bohr)
95 : : !! ffttomix(nfft*(1-nfftmix/nfft))=Index of the points of the FFT (fine) grid on the grid used for mixing (coarse)
96 : : !! gmet(3,3)=metrix tensor in G space in Bohr**-2.
97 : : !! grhf(3,natom)=Hellman-Feynman derivatives of the total energy
98 : : !! gsqcut=cutoff on (k+G)^2 (bohr^-2)
99 : : !! initialized= if 0, the initialization of the gstate run is not yet finished
100 : : !! ispmix=1 if mixing is done in real space, 2 if mixing is done in reciprocal space
101 : : !! istep= number of the step in the SCF cycle
102 : : !! kg_diel(3,npwdiel)=reduced planewave coordinates for the dielectric matrix.
103 : : !! kxc(nfft,nkxc)=exchange-correlation kernel, needed only for electronic
104 : : !! dielectric matrix
105 : : !! mgfft=maximum size of 1D FFTs
106 : : !! mixtofft(nfftmix*(1-nfftmix/nfft))=Index of the points of the FFT grid used for mixing (coarse) on the FFT (fine) grid
107 : : !! moved_atm_inside= if 1, then the preconditioned forces
108 : : !! as well as the preconditioned density residual must be computed;
109 : : !! otherwise, compute only the preconditioned density residual.
110 : : !! mpi_enreg=information about MPI parallelization
111 : : !! my_natom=number of atoms treated by current processor
112 : : !! nattyp(ntypat)=number of atoms of each type in cell.
113 : : !! nfft=(effective) number of FFT grid points (for this processor)
114 : : !! nfftmix=dimension of FFT grid used to mix the densities (used in PAW only)
115 : : !! ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
116 : : !! ngfftmix(18)=contain all needed information about 3D FFT, for the grid corresponding to nfftmix
117 : : !! nkxc=second dimension of the array kxc, see rhotoxc.f for a description
118 : : !! npawmix=-PAW only- number of spherical part elements to be mixed
119 : : !! npwdiel=number of planewaves for dielectric matrix
120 : : !! nresid(nfft,nspden)=array for the residual of the density
121 : : !! ntypat=number of types of atoms in cell.
122 : : !! n1xccc=dimension of xccc1d ; 0 if no XC core correction is used
123 : : !! pawrhoij(my_natom*usepaw) <type(pawrhoij_type)>= paw rhoij occupancies and related data
124 : : !! Use here rhoij residuals (and gradients)
125 : : !! pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
126 : : !! psps <type(pseudopotential_type)>=variables related to pseudopotentials
127 : : !! rprimd(3,3)=dimensional primitive translations in real space (bohr)
128 : : !! susmat(2,npwdiel,nspden,npwdiel,nspden)=
129 : : !! the susceptibility (or density-density response) matrix in reciprocal space
130 : : !! usepaw= 0 for non paw calculation; =1 for paw calculation
131 : : !! vtrial(nfft,nspden)=the trial potential that gave vresid.
132 : : !! xred(3,natom)=reduced dimensionless atomic coordinates
133 : : !! tauresid(nfft,nspden*dtset%usekden)=array for kinetic energy density residue (out - in)
134 : : !!
135 : : !! OUTPUT
136 : : !! dbl_nnsclo=1 if nnsclo has to be doubled to secure the convergence.
137 : : !!
138 : : !! SIDE EFFECTS
139 : : !! dtn_pc(3,natom)=preconditioned change of atomic position,
140 : : !! in reduced coordinates
141 : : !! mix<type(abi_mixing_object)>=all data defining the mixing algorithm for the density
142 : : !! rhor(nfft,nspden)= at input, it is the "out" trial density that gave nresid=(rho_out-rho_in)
143 : : !! at output, it is an updated "mixed" trial density
144 : : !! rhog(2,nfft)= Fourier transform of the new trial density
145 : : !! ===== if usekden==1 =====
146 : : !! [mix_mgga<type(abi_mixing_object)>]=all data defining the mixing algorithm
147 : : !! for the kinetic energy density
148 : : !! ===== if densfor_pred==3 .and. moved_atm_inside==1 =====
149 : : !! ph1d(2,3*(2*mgfft+1)*natom)=1-dim structure factor phases
150 : : !! ==== if usepaw==1
151 : : !! pawrhoij(natom)%nrhoijsel=number of non-zero values of rhoij
152 : : !! pawrhoij(natom)%rhoijp(cplex_rhoij*lmn2_size,nspden)= new (mixed) value of rhoij quantities in PACKED STORAGE
153 : : !! pawrhoij(natom)%rhoijselect(lmn2_size)=select the non-zero values of rhoij
154 : : !! taug(2,nfft*dtset%usekden)=array for Fourier transform of kinetic
155 : : !! energy density
156 : : !! taur(nfft,nspden*dtset%usekden)=array for kinetic energy density
157 : : !!
158 : : !! NOTES
159 : : !! In case of PAW calculations:
160 : : !! Computations are done either on the fine FFT grid or the coarse grid (depending on dtset%pawmixdg)
161 : : !! All variables (nfft,ngfft,mgfft) refer to the fine FFT grid.
162 : : !! All arrays (densities/potentials...) are computed on this fine FFT grid.
163 : : !! ! Developpers have to be careful when introducing others arrays:
164 : : !! they have to be stored on the fine FFT grid (except f_fftgr).
165 : : !! In case of norm-conserving calculations the FFT grid is the usual FFT grid.
166 : : !!
167 : : !! SOURCE
168 : :
169 : 26658 : subroutine newrho(atindx,dbl_nnsclo,dielar,dielinv,dielstrt,dtn_pc,dtset,etotal,fcart,ffttomix,&
170 : 13329 : & gmet,grhf,gsqcut,initialized,ispmix,istep,kg_diel,kxc,mgfft,mix,mixtofft,&
171 : 13329 : & moved_atm_inside,mpi_enreg,my_natom,nattyp,nfft,&
172 : : & nfftmix,nfftmix_per_nfft,ngfft,ngfftmix,nkxc,npawmix,npwdiel,&
173 : 13329 : & nresid,ntypat,n1xccc,pawrhoij,pawtab,&
174 : 13329 : & ph1d,precon,psps,rhog,rhor,rprimd,susmat,usepaw,vtrial,wvl,wvl_den,xred,rcpaw,extfpmd,&
175 : 39987 : & mix_mgga,taug,taur,tauresid)
176 : :
177 : : !Arguments-------------------------------
178 : : !scalars
179 : : integer,intent(in) :: dielstrt,initialized,ispmix,istep,my_natom,mgfft
180 : : integer,intent(in) :: moved_atm_inside,n1xccc,nfft
181 : : integer,intent(in) :: nfftmix,nfftmix_per_nfft
182 : : integer,intent(in) :: nkxc,npawmix,npwdiel,ntypat,usepaw
183 : : integer,intent(inout) :: dbl_nnsclo
184 : : real(dp),intent(in) :: etotal,gsqcut
185 : : type(MPI_type),intent(in) :: mpi_enreg
186 : : type(abi_mixing_object), intent(inout) :: mix
187 : : type(abi_mixing_object), intent(inout),optional :: mix_mgga
188 : : type(dataset_type),intent(in) :: dtset
189 : : type(pseudopotential_type),intent(in) :: psps
190 : : type(wvl_internal_type), intent(in) :: wvl
191 : : type(wvl_denspot_type), intent(inout) :: wvl_den
192 : : type(extfpmd_type), intent(inout), pointer :: extfpmd
193 : : type(rcpaw_type), pointer, intent(inout) :: rcpaw
194 : : type(precon_object), intent(inout) :: precon
195 : : !arrays
196 : : integer,intent(in) :: atindx(dtset%natom)
197 : : integer,intent(in) :: ffttomix(nfft*(nfftmix_per_nfft))
198 : : integer,intent(in) :: kg_diel(3,npwdiel)
199 : : integer,intent(in) :: mixtofft(nfftmix*nfftmix_per_nfft)
200 : : integer,intent(in) :: nattyp(ntypat),ngfft(18),ngfftmix(18)
201 : : real(dp),intent(in) :: dielar(7),fcart(3,dtset%natom),grhf(3,dtset%natom)
202 : : real(dp),intent(inout) :: rprimd(3,3)
203 : : real(dp),intent(in) :: susmat(2,npwdiel,dtset%nspden,npwdiel,dtset%nspden)
204 : : real(dp),intent(in), target :: vtrial(nfft,dtset%nspden)
205 : : real(dp),intent(inout) :: dielinv(2,npwdiel,dtset%nspden,npwdiel,dtset%nspden)
206 : : real(dp),intent(inout), target :: dtn_pc(3,dtset%natom)
207 : : real(dp),intent(inout) :: gmet(3,3)
208 : : !TODO: nresid appears to be only intent in here.
209 : : real(dp),intent(inout) :: kxc(nfft,nkxc),nresid(nfft,dtset%nspden)
210 : : real(dp),intent(inout) :: ph1d(2,3*(2*mgfft+1)*dtset%natom)
211 : : real(dp),intent(inout) :: rhor(nfft,dtset%nspden)
212 : : real(dp),intent(inout), target :: xred(3,dtset%natom)
213 : : real(dp),intent(inout) :: rhog(2,nfft)
214 : : real(dp),intent(inout), optional :: taug(2,nfft*dtset%usekden)
215 : : real(dp),intent(inout), optional :: taur(nfft,dtset%nspden*dtset%usekden)
216 : : real(dp),intent(inout), optional :: tauresid(nfft,dtset%nspden*dtset%usekden)
217 : :
218 : : type(pawrhoij_type),intent(inout) :: pawrhoij(my_natom*psps%usepaw)
219 : : type(pawtab_type),intent(in) :: pawtab(ntypat*psps%usepaw)
220 : :
221 : : !Local variables-------------------------------
222 : : !scalars
223 : : integer,parameter :: tim_fourdp9=9
224 : : integer :: cplex,dplex,errid,i_vresid1,i_vrespc1,iatom,ifft,indx,iq,iq0,irhoij,ispden,jfft
225 : : integer :: jrhoij,klmn,kklmn,kmix,mpicomm,nfftot,qphase,itypat,iln,isppol
226 : : logical :: mpi_summarize,reset
227 : : real(dp) :: fact,ucvol,ucvol_local
228 : : character(len=500) :: message
229 : : !arrays
230 : : real(dp) :: gprimd(3,3),rmet(3,3),ro(2),tsec(2),vhartr_dum(1),vpsp_dum(1)
231 : : real(dp) :: vxc_dum(1,1)
232 : : real(dp) :: nelect_extfpmd_=zero
233 : : real(dp),target :: dum0(0)
234 : 13329 : real(dp),allocatable :: magng(:,:,:),magntaug(:,:,:)
235 : 13329 : real(dp),allocatable :: nresid0(:,:),nrespc(:,:),nreswk(:,:,:)
236 : 13329 : real(dp),allocatable :: rhoijrespc(:),rhoijtmp(:,:)
237 : : ! TODO : these should be allocatables not pointers: is there some reason to
238 : : ! keep them this way, eg an interface somewhere?
239 : 13329 : real(dp), pointer :: rhomag(:,:), npaw(:)
240 : 13329 : real(dp),allocatable :: tauresid0(:,:),taurespc(:,:)
241 : 13329 : real(dp),allocatable :: taumag(:,:)
242 : 13329 : real(dp), pointer :: rcpaw_arr_(:)
243 : 13329 : real(dp),allocatable, target :: rcpaw_arr(:)
244 : :
245 : : ! *************************************************************************
246 : :
247 : : DBG_ENTER("COLL")
248 : 13329 : call timab(94,1,tsec)
249 : :
250 [ + + ]: 53316 : nfftot=PRODUCT(ngfft(1:3))
251 : :
252 : : !Compatibility tests
253 [ - + ]: 13329 : if(nfftmix>nfft) then
254 : 0 : message='nfftmix>nfft not allowed!'
255 : 0 : ABI_BUG(message)
256 : : end if
257 : :
258 [ - + ]: 13329 : if(dtset%usewvl==1) then
259 [ # # # # ]: 0 : if( (ispmix/=1 .or. nfftmix/=nfft)) then
260 : 0 : message='nfftmix/=nfft, ispmix/=1 not allowed for wavelets!'
261 : 0 : ABI_BUG(message)
262 : : end if
263 [ # # ]: 0 : if(dtset%wvl_bigdft_comp==1) then
264 : 0 : message='usewvl == 1 and wvl_bigdft_comp==1 not allowed!'
265 : 0 : ABI_BUG(message)
266 : : end if
267 : : end if
268 : :
269 [ + + - + ]: 13329 : if(ispmix/=2.and.nfftmix/=nfft) then
270 : 0 : message='nfftmix/=nfft allowed only when ispmix=2!'
271 : 0 : ABI_BUG(message)
272 : : end if
273 : :
274 [ + + ]: 13329 : if (dtset%usekden==1) then
275 : : if ((.not.present(tauresid)).or.(.not.present(taug)).or. &
276 [ + - - + ]: 464 : & (.not.present(taur)).or.(.not.present(mix_mgga))) then
277 : 0 : message='Several arrays are missing!'
278 : 0 : ABI_BUG(message)
279 : : end if
280 [ - + ]: 464 : if (mix_mgga%iscf==ABI_MIXING_CG_ENERGY.or.mix_mgga%iscf==ABI_MIXING_CG_ENERGY_2.or.&
281 : : & mix_mgga%iscf==ABI_MIXING_EIG) then
282 : 0 : message='kinetic energy density cannot be mixed with the selected mixing algorithm!'
283 : 0 : ABI_ERROR(message)
284 : : end if
285 : : end if
286 : :
287 [ + + + + ]: 13329 : if (usepaw==1.and.my_natom>0) then
288 : 11668 : cplex=pawrhoij(1)%cplex_rhoij;dplex=cplex-1
289 : 11668 : qphase=pawrhoij(1)%qphase
290 : : else
291 : : cplex = 0;dplex = 0 ; qphase=0
292 : : end if
293 : :
294 : : !Compute different geometric tensor, as well as ucvol, from rprimd
295 : 13329 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
296 : :
297 [ + - ]: 13329 : if(dtset%usewvl==0) then
298 : 13329 : ucvol_local=ucvol
299 : : #if defined HAVE_BIGDFT
300 : : else
301 : : ucvol_local = product(wvl_den%denspot%dpbox%hgrids) * real(nfftot, dp)
302 : : #endif
303 : : end if
304 : :
305 : : !Select components of density to be mixed
306 [ - + + - : 53316 : ABI_MALLOC(rhomag,(ispmix*nfftmix,dtset%nspden))
- + + - +
- + - - +
- + ]
307 [ - + + - : 39987 : ABI_MALLOC(nresid0,(ispmix*nfftmix,dtset%nspden))
- + + - -
+ - + ]
308 [ - + + - : 53316 : ABI_MALLOC(taumag,(ispmix*nfftmix,dtset%nspden*dtset%usekden))
+ + + - +
- + + - +
- + ]
309 [ - + + - : 39987 : ABI_MALLOC(tauresid0,(ispmix*nfftmix,dtset%nspden*dtset%usekden))
+ + + - -
+ - + ]
310 : : ! real space and all fft points are here
311 [ + + + - ]: 13329 : if (ispmix==1.and.nfft==nfftmix) then
312 [ + + + + ]: 93998175 : rhomag(:,1:dtset%nspden)=rhor(:,1:dtset%nspden)
313 [ + + + + ]: 93998175 : nresid0(:,1:dtset%nspden)=nresid(:,1:dtset%nspden)
314 [ + + ]: 1500 : if (dtset%usekden==1) then
315 [ + + + + ]: 7768146 : taumag(:,1:dtset%nspden)=taur(:,1:dtset%nspden)
316 [ + + + + ]: 7768146 : tauresid0(:,1:dtset%nspden)=tauresid(:,1:dtset%nspden)
317 : : end if
318 : : ! recip space and all fft points are here
319 [ + + ]: 11829 : else if (nfft==nfftmix) then
320 [ + + ]: 1484 : do ispden=1,dtset%nspden
321 : 1484 : call fourdp(1,nresid0(:,ispden),nresid(:,ispden),-1,mpi_enreg,nfft,1,ngfft,tim_fourdp9)
322 : : end do
323 [ + + ]: 1368 : rhomag(:,1)=reshape(rhog,(/2*nfft/))
324 [ + + ]: 684 : if (dtset%nspden>1) then
325 [ + + ]: 232 : do ispden=2,dtset%nspden
326 [ - + ]: 232 : call fourdp(1,rhomag(:,ispden),rhor(:,ispden),-1,mpi_enreg,nfft,1,ngfft,tim_fourdp9)
327 : : end do
328 : : end if
329 [ + - ]: 684 : if (dtset%usekden==1) then
330 [ # # ]: 0 : do ispden=1,dtset%nspden
331 [ # # ]: 0 : call fourdp(1,tauresid0(:,ispden),tauresid(:,ispden),-1,mpi_enreg,nfft,1,ngfft,tim_fourdp9)
332 : : end do
333 [ # # # # ]: 0 : taumag(:,1)=reshape(taug,(/2*nfft/))
334 [ # # ]: 0 : if (dtset%nspden>1) then
335 [ # # ]: 0 : do ispden=2,dtset%nspden
336 [ # # ]: 0 : call fourdp(1,taumag(:,ispden),taur(:,ispden),-1,mpi_enreg,nfft,1,ngfft,tim_fourdp9)
337 : : end do
338 : : end if
339 : : end if
340 : : ! not all fft points are here - presumes recip space
341 : : else
342 : 11145 : fact=dielar(4)-1._dp
343 [ - + + - : 44580 : ABI_MALLOC(nreswk,(2,nfft,dtset%nspden))
- + + - +
- + - - +
+ - ]
344 [ + + ]: 26537 : do ispden=1,dtset%nspden
345 : 26537 : call fourdp(1,nreswk(:,:,ispden),nresid(:,ispden),-1,mpi_enreg,nfft,1,ngfft,tim_fourdp9)
346 : : end do
347 [ + + ]: 540211923 : do ifft=1,nfft
348 [ + + ]: 540211923 : if (ffttomix(ifft)>0) then
349 : 185287147 : jfft=2*ffttomix(ifft)
350 [ + + ]: 555861441 : rhomag (jfft-1:jfft,1)=rhog(1:2,ifft)
351 [ + + ]: 555861441 : nresid0(jfft-1:jfft,1)=nreswk(1:2,ifft,1)
352 : : else
353 [ + + ]: 1064740893 : rhog(:,ifft)=rhog(:,ifft)+fact*nreswk(:,ifft,1)
354 : : end if
355 : : end do
356 [ + + ]: 11145 : if (dtset%nspden>1) then
357 [ - + + - : 9644 : ABI_MALLOC(magng,(2,nfft,dtset%nspden-1))
- + + - +
- + - - +
+ - ]
358 [ + + ]: 6658 : do ispden=2,dtset%nspden
359 : 4247 : call fourdp(1,magng(:,:,ispden-1),rhor(:,ispden),-1,mpi_enreg,nfft,1,ngfft,tim_fourdp9)
360 [ + + ]: 102557167 : do ifft=1,nfft
361 [ + + ]: 102554756 : if (ffttomix(ifft)>0) then
362 : 29844700 : jfft=2*ffttomix(ifft)
363 [ + + ]: 89534100 : rhomag (jfft-1:jfft,ispden)=magng (1:2,ifft,ispden-1)
364 [ + + ]: 89534100 : nresid0(jfft-1:jfft,ispden)=nreswk(1:2,ifft,ispden)
365 : : else
366 [ + + ]: 218117427 : magng(:,ifft,ispden-1)=magng(:,ifft,ispden-1)+fact*nreswk(:,ifft,ispden)
367 [ + + + + ]: 167734527 : if (dtset%nspden==2) magng(:,ifft,1)=two*magng(:,ifft,1)-rhog(:,ifft)
368 : : end if
369 : : end do
370 : : end do
371 : : end if
372 [ + + ]: 11145 : if (dtset%usekden==1) then
373 [ + + ]: 846 : do ispden=1,dtset%nspden
374 [ + - ]: 846 : call fourdp(1,nreswk(:,:,ispden),tauresid(:,ispden),-1,mpi_enreg,nfft,1,ngfft,tim_fourdp9)
375 : : end do
376 [ + + ]: 19013607 : do ifft=1,nfft
377 [ + + ]: 19013607 : if (ffttomix(ifft)>0) then
378 : 7982995 : jfft=2*ffttomix(ifft)
379 [ + + ]: 23948985 : taumag (jfft-1:jfft,1)=taug(1:2,ifft)
380 [ + + ]: 23948985 : tauresid0(jfft-1:jfft,1)=nreswk(1:2,ifft,1)
381 : : else
382 [ + + ]: 33090567 : taug(:,ifft)=taug(:,ifft)+fact*nreswk(:,ifft,1)
383 : : end if
384 : : end do
385 [ - + ]: 423 : if (dtset%nspden>1) then
386 [ # # # # : 0 : ABI_MALLOC(magntaug,(2,nfft,dtset%nspden-1))
# # # # #
# # # # #
# # ]
387 [ # # ]: 0 : do ispden=2,dtset%nspden
388 [ # # ]: 0 : call fourdp(1,magntaug(:,:,ispden-1),taur(:,ispden),-1,mpi_enreg,nfft,1,ngfft,tim_fourdp9)
389 [ # # ]: 0 : do ifft=1,nfft
390 [ # # ]: 0 : if (ffttomix(ifft)>0) then
391 : 0 : jfft=2*ffttomix(ifft)
392 [ # # ]: 0 : taumag (jfft-1:jfft,ispden)=magntaug(1:2,ifft,ispden-1)
393 [ # # ]: 0 : tauresid0(jfft-1:jfft,ispden)=nreswk(1:2,ifft,ispden)
394 : : else
395 [ # # ]: 0 : magntaug(:,ifft,ispden-1)=magntaug(:,ifft,ispden-1)+fact*nreswk(:,ifft,ispden)
396 [ # # # # ]: 0 : if (dtset%nspden==2) magntaug(:,ifft,1)=two*magntaug(:,ifft,1)-taug(:,ifft)
397 : : end if
398 : : end do
399 : : end do
400 : : end if
401 : : end if
402 [ - + ]: 11145 : ABI_FREE(nreswk)
403 : : end if
404 : :
405 : : !Retrieve "input" density from "output" density and density residual
406 [ + + + + ]: 544325106 : rhomag(:,1:dtset%nspden)=rhomag(:,1:dtset%nspden)-nresid0(:,1:dtset%nspden)
407 [ + + ]: 13329 : if (dtset%usekden==1) then
408 [ + + + + ]: 23734982 : taumag(:,1:dtset%nspden)=taumag(:,1:dtset%nspden)-tauresid0(:,1:dtset%nspden)
409 : : end if
410 : :
411 : : !If nspden==2, separate density and magnetization
412 [ + + ]: 13329 : if (dtset%nspden==2) then
413 [ + + ]: 41040403 : rhomag (:,2)=two*rhomag (:,2)-rhomag (:,1)
414 [ + + ]: 41040403 : nresid0(:,2)=two*nresid0(:,2)-nresid0(:,1)
415 [ - + ]: 1693 : if (dtset%usekden==1) then
416 [ # # ]: 0 : taumag (:,2)=two*taumag (:,2)-taumag (:,1)
417 [ # # ]: 0 : tauresid0(:,2)=two*tauresid0(:,2)-tauresid0(:,1)
418 : : end if
419 : : end if
420 : :
421 : : !If PAW, handle occupancy matrix
422 [ + + + + ]: 13329 : if (usepaw==1.and.my_natom>0) then
423 [ + + ]: 11668 : if (pawrhoij(1)%nspden==2) then
424 [ + + ]: 5669 : do iatom=1,my_natom
425 [ + + ]: 9713 : do iq=1,qphase
426 [ - + ]: 4044 : iq0=merge(0,cplex*pawrhoij(iatom)%lmn2_size,iq==1)
427 : 4044 : jrhoij=1+iq0
428 [ + + ]: 236442 : do irhoij=1,pawrhoij(iatom)%nrhoijsel
429 [ + + ]: 465393 : ro(1:1+dplex)=pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+dplex,1)
430 [ + + ]: 465393 : pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+dplex,1)=ro(1:1+dplex)+pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+dplex,2)
431 [ + + ]: 465393 : pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+dplex,2)=ro(1:1+dplex)-pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+dplex,2)
432 : 236442 : jrhoij=jrhoij+cplex
433 : : end do
434 [ + + ]: 496029 : do kmix=1,pawrhoij(iatom)%lmnmix_sz
435 : 487941 : klmn=iq0+cplex*pawrhoij(iatom)%kpawmix(kmix)-dplex
436 [ + + ]: 976530 : ro(1:1+dplex)=pawrhoij(iatom)%rhoijres(klmn:klmn+dplex,1)
437 [ + + ]: 976530 : pawrhoij(iatom)%rhoijres(klmn:klmn+dplex,1)=ro(1:1+dplex)+pawrhoij(iatom)%rhoijres(klmn:klmn+dplex,2)
438 [ + + ]: 980574 : pawrhoij(iatom)%rhoijres(klmn:klmn+dplex,2)=ro(1:1+dplex)-pawrhoij(iatom)%rhoijres(klmn:klmn+dplex,2)
439 : : end do
440 : : end do
441 : : end do
442 : : end if
443 : : end if
444 : :
445 : : !Choice of preconditioner governed by iprcel, densfor_pred and iprcfc
446 [ - + + - : 39987 : ABI_MALLOC(nrespc,(ispmix*nfftmix,dtset%nspden))
- + + - -
+ - + ]
447 [ - + + - : 39987 : ABI_MALLOC(taurespc,(ispmix*nfftmix,dtset%nspden*dtset%usekden))
+ + + - -
+ - + ]
448 [ + + + - : 39987 : ABI_MALLOC(npaw,(npawmix*usepaw))
+ - + + -
+ - + ]
449 [ + + ]: 13329 : if (usepaw==1) then
450 [ + + + - : 39762 : ABI_MALLOC(rhoijrespc,(npawmix))
+ - + + -
+ - + ]
451 : : else
452 [ - + ]: 75 : ABI_MALLOC(rhoijrespc,(0))
453 : : end if
454 [ + - ]: 13329 : if(dtset%usewvl==0) then
455 : : call prcref(atindx,dielar,dielinv,&
456 : : & dielstrt,dtn_pc,dtset,etotal,fcart,ffttomix,gmet,gsqcut,&
457 : : & istep,kg_diel,kxc,&
458 : : & mgfft,moved_atm_inside,mpi_enreg,my_natom,&
459 : : & nattyp,nfft,nfftmix,ngfft,ngfftmix,nkxc,npawmix,npwdiel,ntypat,n1xccc,&
460 : : & ispmix,1,pawrhoij,pawtab,ph1d,precon,psps,rhog,rhoijrespc,rhor,rprimd,&
461 : 13329 : & susmat,vhartr_dum,vpsp_dum,nresid0,nrespc,vxc_dum,wvl,wvl_den,xred,rcpaw,extfpmd)
462 : : else
463 : : call wvl_prcref(dielar,dtset%iprcel,my_natom,nfftmix,npawmix,dtset%nspden,pawrhoij,&
464 : 0 : & rhoijrespc,psps%usepaw,nresid0,nrespc)
465 : : end if
466 : : !At present, only a simple precoditionning for the kinetic energy density
467 : : ! (is Kerker mixing valid for tau?)
468 [ + + ]: 13329 : if (dtset%usekden==1) then
469 [ + + ]: 928 : do ispden=1,dtset%nspden
470 [ - + ]: 464 : fact=dielar(4);if (ispden>1) fact=abs(dielar(7))
471 [ + + ]: 23734982 : taurespc(1:ispmix*nfftmix,ispden)=fact*tauresid0(1:ispmix*nfftmix,ispden)
472 : : end do
473 : : end if
474 : :
475 : : !------Compute new trial density and eventual new atomic positions
476 : :
477 [ + - ]: 13329 : if (mix%n_fftgr>0) then
478 : 13329 : i_vresid1=mix%i_vresid(1)
479 : 13329 : i_vrespc1=mix%i_vrespc(1)
480 : : end if
481 : :
482 [ + + ]: 13329 : if(associated(extfpmd)) then
483 : 83 : mix%useextfpmd=1
484 : : endif
485 [ + + ]: 13329 : if(associated(rcpaw)) then
486 : 18 : mix%use_rcpaw=1
487 : 18 : mix%n_rcpawmix=0
488 [ + + ]: 36 : do itypat=1,size(rcpaw%atm)
489 : 36 : mix%n_rcpawmix=mix%n_rcpawmix+rcpaw%atm(itypat)%ln_size*rcpaw%atm(itypat)%nsppol
490 : : enddo
491 [ - + + - : 54 : ABI_MALLOC(rcpaw_arr,(mix%n_rcpawmix))
+ - + - -
+ - + ]
492 : : endif
493 : :
494 : : !Initialise working arrays for the mixing object.
495 [ - + ]: 13329 : if (moved_atm_inside == 1) then
496 : 0 : call abi_mixing_use_moving_atoms(mix, dtset%natom, xred, dtn_pc)
497 : : end if
498 : 13329 : call abi_mixing_eval_allocate(mix, istep)
499 : :
500 : : !Copy current step arrays.
501 [ - + ]: 13329 : if (moved_atm_inside == 1) then
502 : 0 : call abi_mixing_copy_current_step(mix, nresid0, errid, message, arr_respc = nrespc, arr_atm = grhf)
503 : : else
504 : 13329 : call abi_mixing_copy_current_step(mix, nresid0, errid, message, arr_respc = nrespc)
505 : : end if
506 [ - + ]: 13329 : if (errid /= AB7_NO_ERROR) then
507 : 0 : ABI_ERROR(message)
508 : : end if
509 : :
510 : : !Same treatment for the kinetic energy density
511 [ + + ]: 13329 : if (dtset%usekden==1) then
512 : 464 : call abi_mixing_eval_allocate(mix_mgga, istep)
513 : 464 : call abi_mixing_copy_current_step(mix_mgga, tauresid0, errid, message, arr_respc = taurespc)
514 [ - + ]: 464 : if (errid /= AB7_NO_ERROR) then
515 : 0 : ABI_ERROR(message)
516 : : end if
517 : : end if
518 : :
519 [ - + ]: 13329 : ABI_FREE(nresid0)
520 [ - + ]: 13329 : ABI_FREE(nrespc)
521 [ - + ]: 13329 : ABI_FREE(tauresid0)
522 [ - + ]: 13329 : ABI_FREE(taurespc)
523 : :
524 : : !PAW: either use the array f_paw or the array f_paw_disk
525 [ + + ]: 13329 : if (usepaw==1) then
526 : 13254 : indx=-dplex
527 [ + + ]: 43946 : do iatom=1,my_natom
528 [ - + + - : 92076 : ABI_MALLOC(rhoijtmp,(cplex*pawrhoij(iatom)%lmn2_size,1))
+ - + - -
+ - + ]
529 [ + + ]: 61384 : do iq=1,qphase
530 [ - + ]: 30692 : iq0=merge(0,cplex*pawrhoij(iatom)%lmn2_size,iq==1)
531 [ + + ]: 103773 : do ispden=1,pawrhoij(iatom)%nspden
532 [ + + + + ]: 4026226 : rhoijtmp=zero ; jrhoij=1+iq0
533 [ + + ]: 2176360 : do irhoij=1,pawrhoij(iatom)%nrhoijsel
534 : 2133971 : klmn=cplex*pawrhoij(iatom)%rhoijselect(irhoij)-dplex
535 [ + + ]: 5127403 : rhoijtmp(klmn:klmn+dplex,1)=pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+dplex,ispden)
536 : 2176360 : jrhoij=jrhoij+cplex
537 : : end do
538 [ + + ]: 3098105 : do kmix=1,pawrhoij(iatom)%lmnmix_sz
539 : 3025024 : indx=indx+cplex;klmn=cplex*pawrhoij(iatom)%kpawmix(kmix)-dplex ; kklmn=klmn+iq0
540 [ + + ]: 6964312 : npaw(indx:indx+dplex)=rhoijtmp(klmn:klmn+dplex,1)-pawrhoij(iatom)%rhoijres(kklmn:kklmn+dplex,ispden)
541 [ + + ]: 6964312 : mix%f_paw(indx:indx+dplex,i_vresid1)=pawrhoij(iatom)%rhoijres(kklmn:kklmn+dplex,ispden)
542 [ + + ]: 7006701 : mix%f_paw(indx:indx+dplex,i_vrespc1)=rhoijrespc(indx:indx+dplex)
543 : : end do
544 : : end do
545 : : end do
546 [ - + ]: 43946 : ABI_FREE(rhoijtmp)
547 : : end do
548 : : end if
549 : :
550 : 13329 : rcpaw_arr_=>dum0
551 [ + + ]: 13329 : if(associated(rcpaw)) then
552 : 18 : indx=0
553 [ + + ]: 36 : do itypat=1,size(rcpaw%atm)
554 [ + + ]: 54 : do isppol=1,rcpaw%atm(itypat)%nsppol
555 [ + + ]: 90 : do iln=1,rcpaw%atm(itypat)%ln_size
556 : 54 : indx=indx+1
557 : 54 : mix%f_rcpaw(indx,i_vresid1)=rcpaw%atm(itypat)%occ_res(iln,isppol)
558 : 54 : mix%f_rcpaw(indx,i_vrespc1)=rcpaw%atm(itypat)%occ_respc(iln,isppol)
559 : 72 : rcpaw_arr(indx)=rcpaw%atm(itypat)%occ(iln,isppol)-rcpaw%atm(itypat)%occ_res(iln,isppol)
560 : : enddo
561 : : enddo
562 : : enddo
563 : 18 : rcpaw_arr_=>rcpaw_arr
564 : : endif
565 : :
566 [ + + ]: 13329 : if(associated(extfpmd)) then
567 : 83 : mix%f_extfpmd(i_vresid1)=extfpmd%nelect_res
568 : 83 : mix%f_extfpmd(i_vrespc1)=extfpmd%nelect_respc
569 : 83 : nelect_extfpmd_=extfpmd%nelect-extfpmd%nelect_res
570 : : endif
571 : :
572 : :
573 : :
574 : : !------Prediction of the components of the density
575 : :
576 : : !Init mpicomm
577 [ + + ]: 13329 : if(mpi_enreg%paral_kgb==1)then
578 : 4504 : mpicomm=mpi_enreg%comm_fft
579 : 4504 : mpi_summarize=.true.
580 : : else
581 : 8825 : mpicomm=0
582 : 8825 : mpi_summarize=.false.
583 : : end if
584 [ - + ]: 13329 : if(dtset%usewvl==1) then
585 : 0 : mpicomm=mpi_enreg%comm_wvl
586 : 0 : mpi_summarize=(mpi_enreg%nproc_wvl > 1)
587 : : end if
588 : :
589 : 13329 : reset = .false.
590 [ + + ]: 13329 : if (initialized == 0) reset = .true.
591 : :
592 : : !Electronic density mixing
593 : : call abi_mixing_eval(mix, rhomag, istep, nfftot, ucvol_local, &
594 : : & mpicomm, mpi_summarize, errid, message, &
595 : : & reset = reset, isecur = dtset%isecur,&
596 : : & pawopt = dtset%pawoptmix, pawarr = npaw, &
597 : : & etotal = etotal, potden = vtrial, &
598 : : & nelect_extfpmd = nelect_extfpmd_,&
599 : : & rcpaw_arr = rcpaw_arr_ ,&
600 [ - + - + : 13329 : & comm_atom=mpi_enreg%comm_atom)
- + ]
601 [ - + ]: 13329 : if (errid == AB7_ERROR_MIXING_INC_NNSLOOP) then
602 : 0 : dbl_nnsclo = 1
603 [ - + ]: 13329 : else if (errid /= AB7_NO_ERROR) then
604 : 0 : ABI_ERROR(message)
605 : : end if
606 : : !Kinetic energy density mixing (if any)
607 [ + + ]: 13329 : if (dtset%usekden==1) then
608 : : call abi_mixing_eval(mix_mgga, taumag, istep, nfftot, ucvol_local, &
609 : 464 : & mpicomm, mpi_summarize, errid, message, reset = reset)
610 [ - + ]: 464 : if (errid /= AB7_NO_ERROR) then
611 : 0 : ABI_ERROR(message)
612 : : end if
613 : : end if
614 : :
615 [ + + ]: 13329 : if(associated(rcpaw)) then
616 : 18 : indx=0
617 : 18 : rcpaw%nelect_core=zero
618 [ + + ]: 36 : do itypat=1,size(rcpaw%atm)
619 : 18 : rcpaw%atm(itypat)%zcore=zero
620 [ + + ]: 36 : do isppol=1,rcpaw%atm(itypat)%nsppol
621 [ + + ]: 90 : do iln=1,rcpaw%atm(itypat)%ln_size
622 : 54 : indx=indx+1
623 : 54 : rcpaw%atm(itypat)%occ(iln,isppol)=rcpaw_arr_(indx)
624 : 72 : rcpaw%atm(itypat)%zcore=rcpaw%atm(itypat)%zcore+rcpaw%atm(itypat)%occ(iln,isppol)
625 : : enddo
626 : : enddo
627 : 36 : rcpaw%nelect_core=rcpaw%nelect_core+rcpaw%atm(itypat)%zcore*rcpaw%atm(itypat)%mult
628 : : enddo
629 [ - + ]: 18 : ABI_FREE(rcpaw_arr)
630 : : endif
631 : 13329 : nullify(rcpaw_arr_)
632 : :
633 [ + + ]: 13329 : if(associated(extfpmd)) then
634 : 83 : extfpmd%nelect=nelect_extfpmd_
635 : : endif
636 : :
637 : : !PAW: apply a simple mixing to rhoij (this is temporary)
638 [ + + ]: 13329 : if(dtset%iscf==15 .or. dtset%iscf==16)then
639 [ - + ]: 12 : if (usepaw==1) then
640 : 0 : indx=-dplex
641 [ # # ]: 0 : do iatom=1,my_natom
642 [ # # # # : 0 : ABI_MALLOC(rhoijtmp,(cplex*qphase*pawrhoij(iatom)%lmn2_size,pawrhoij(iatom)%nspden))
# # # # #
# # # # #
# # ]
643 [ # # # # ]: 0 : rhoijtmp=zero
644 [ # # ]: 0 : do iq=1,qphase
645 [ # # ]: 0 : iq0=merge(0,cplex*pawrhoij(iatom)%lmn2_size,iq==1)
646 [ # # ]: 0 : if (pawrhoij(iatom)%lmnmix_sz<pawrhoij(iatom)%lmn2_size) then
647 [ # # ]: 0 : do ispden=1,pawrhoij(iatom)%nspden
648 [ # # ]: 0 : do kmix=1,pawrhoij(iatom)%lmnmix_sz
649 : 0 : indx=indx+cplex;klmn=iq0+cplex*pawrhoij(iatom)%kpawmix(kmix)-dplex
650 : : rhoijtmp(klmn:klmn+dplex,ispden)=rhoijrespc(indx:indx+dplex) &
651 [ # # ]: 0 : & -pawrhoij(iatom)%rhoijres(klmn:klmn+dplex,ispden)
652 : : end do
653 : : end do
654 : : end if
655 [ # # ]: 0 : if (pawrhoij(iatom)%nspden/=2) then
656 [ # # ]: 0 : do ispden=1,pawrhoij(iatom)%nspden
657 : 0 : jrhoij=iq0+1
658 [ # # ]: 0 : do irhoij=1,pawrhoij(iatom)%nrhoijsel
659 : 0 : klmn=iq0+cplex*pawrhoij(iatom)%rhoijselect(irhoij)-dplex
660 : : rhoijtmp(klmn:klmn+dplex,ispden)=rhoijtmp(klmn:klmn+dplex,ispden) &
661 [ # # ]: 0 : & +pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+dplex,ispden)
662 : 0 : jrhoij=jrhoij+cplex
663 : : end do
664 : : end do
665 : : else
666 : 0 : jrhoij=iq0+1
667 [ # # ]: 0 : do irhoij=1,pawrhoij(iatom)%nrhoijsel
668 : 0 : klmn=iq0+cplex*pawrhoij(iatom)%rhoijselect(irhoij)-dplex
669 [ # # ]: 0 : ro(1:1+dplex)=rhoijtmp(klmn:klmn+dplex,1)
670 : : rhoijtmp(klmn:klmn+dplex,1)=half*(ro(1:1+dplex)+rhoijtmp(klmn:klmn+dplex,2)) &
671 [ # # ]: 0 : & +pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+dplex,1)
672 : : rhoijtmp(klmn:klmn+dplex,2)=half*(ro(1:1+dplex)-rhoijtmp(klmn:klmn+dplex,2)) &
673 [ # # ]: 0 : & +pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+dplex,2)
674 : 0 : jrhoij=jrhoij+cplex
675 : : end do
676 : : end if
677 : : end do
678 : : call pawrhoij_filter(pawrhoij(iatom)%rhoijp,pawrhoij(iatom)%rhoijselect,&
679 : : & pawrhoij(iatom)%nrhoijsel,pawrhoij(iatom)%cplex_rhoij,pawrhoij(iatom)%qphase,&
680 : 0 : & pawrhoij(iatom)%lmn2_size,pawrhoij(iatom)%nspden,rhoij_input=rhoijtmp)
681 [ # # ]: 0 : ABI_FREE(rhoijtmp)
682 : : end do
683 : : end if
684 : : end if
685 : :
686 : : !if (usepaw==1) then
687 [ - + ]: 13329 : ABI_FREE(rhoijrespc)
688 : : !end if
689 : :
690 : : !PAW: restore rhoij from compact storage
691 [ + + + - : 13329 : if (usepaw==1.and.dtset%iscf/=15.and.dtset%iscf/=16) then
+ - ]
692 : 13254 : indx=-dplex
693 [ + + ]: 43946 : do iatom=1,my_natom
694 [ - + + - : 122768 : ABI_MALLOC(rhoijtmp,(cplex*qphase*pawrhoij(iatom)%lmn2_size,pawrhoij(iatom)%nspden))
- + + - +
- + - - +
- + ]
695 [ + + + + ]: 4014529 : rhoijtmp=zero
696 [ + + ]: 61384 : do iq=1,qphase
697 [ - + ]: 30692 : iq0=merge(0,cplex*pawrhoij(iatom)%lmn2_size,iq==1)
698 [ + + ]: 30692 : if (pawrhoij(iatom)%lmnmix_sz<pawrhoij(iatom)%lmn2_size) then
699 [ + + ]: 24 : do ispden=1,pawrhoij(iatom)%nspden
700 : 16 : jrhoij=iq0+1
701 [ + + ]: 456 : do irhoij=1,pawrhoij(iatom)%nrhoijsel
702 : 432 : klmn=iq0+cplex*pawrhoij(iatom)%rhoijselect(irhoij)-dplex
703 [ + + ]: 864 : rhoijtmp(klmn:klmn+dplex,ispden)=pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+dplex,ispden)
704 : 448 : jrhoij=jrhoij+cplex
705 : : end do
706 : : end do
707 : : end if
708 [ + + ]: 73081 : do ispden=1,pawrhoij(iatom)%nspden
709 [ + + ]: 3098105 : do kmix=1,pawrhoij(iatom)%lmnmix_sz
710 : 3025024 : indx=indx+cplex;klmn=iq0+cplex*pawrhoij(iatom)%kpawmix(kmix)-dplex
711 [ + + ]: 7006701 : rhoijtmp(klmn:klmn+dplex,ispden)=npaw(indx:indx+dplex)
712 : : end do
713 : : end do
714 [ + + ]: 61384 : if (pawrhoij(iatom)%nspden==2) then
715 : 4044 : jrhoij=iq0+1
716 [ + + ]: 236442 : do irhoij=1,pawrhoij(iatom)%nrhoijsel
717 : 232398 : klmn=iq0+cplex*pawrhoij(iatom)%rhoijselect(irhoij)-dplex
718 [ + + ]: 465393 : ro(1:1+dplex)=rhoijtmp(klmn:klmn+dplex,1)
719 [ + + ]: 465393 : rhoijtmp(klmn:klmn+dplex,1)=half*(ro(1:1+dplex)+rhoijtmp(klmn:klmn+dplex,2))
720 [ + + ]: 465393 : rhoijtmp(klmn:klmn+dplex,2)=half*(ro(1:1+dplex)-rhoijtmp(klmn:klmn+dplex,2))
721 : 236442 : jrhoij=jrhoij+cplex
722 : : end do
723 : : end if
724 : : end do
725 : : call pawrhoij_filter(pawrhoij(iatom)%rhoijp,pawrhoij(iatom)%rhoijselect,&
726 : : & pawrhoij(iatom)%nrhoijsel,pawrhoij(iatom)%cplex_rhoij,pawrhoij(iatom)%qphase,&
727 : 30692 : & pawrhoij(iatom)%lmn2_size,pawrhoij(iatom)%nspden,rhoij_input=rhoijtmp)
728 [ - + ]: 43946 : ABI_FREE(rhoijtmp)
729 : : end do
730 : : end if ! usepaw==1.and.dtset%iscf/=15.and.dtset%iscf/=16
731 [ - + ]: 13329 : ABI_FREE(npaw)
732 : :
733 : : !Eventually write the data on disk and deallocate f_fftgr_disk
734 : 13329 : call abi_mixing_eval_deallocate(mix)
735 [ + + ]: 13329 : if (dtset%usekden==1) call abi_mixing_eval_deallocate(mix_mgga)
736 : :
737 : : !Fourier transform the density
738 [ + + + - ]: 13329 : if (ispmix==1.and.nfft==nfftmix) then
739 : : !Real space mixing, no need to transform rhomag
740 [ + + + + ]: 93998175 : rhor(:,1:dtset%nspden)=rhomag(:,1:dtset%nspden)
741 [ + - ]: 1500 : if(dtset%usewvl==0) then
742 : : !Get rhog from rhor(:,1)
743 : 1500 : call fourdp(1,rhog,rhor(:,1),-1,mpi_enreg,nfft,1,ngfft,tim_fourdp9)
744 : : end if
745 [ + + ]: 1500 : if (dtset%usekden==1) then
746 [ + + + + ]: 7768146 : taur(:,1:dtset%nspden)=taumag(:,1:dtset%nspden)
747 [ + - ]: 41 : if(dtset%usewvl==0) then
748 [ + - ]: 41 : call fourdp(1,taug,taur(:,1),-1,mpi_enreg,nfft,1,ngfft,tim_fourdp9)
749 : : end if
750 : : end if
751 [ + + ]: 11829 : else if (nfft==nfftmix) then
752 : : !Reciprocal mixing space mixing, need to generate rhor in real space from rhomag in reciprocal space
753 [ + + ]: 1484 : do ispden=1,dtset%nspden
754 [ - + ]: 1484 : call fourdp(1,rhomag(:,ispden),rhor(:,ispden),+1,mpi_enreg,nfft,1,ngfft,tim_fourdp9)
755 : : end do
756 [ + + ]: 2052 : rhog(:,:)=reshape(rhomag(:,1),(/2,nfft/))
757 [ + - ]: 684 : if (dtset%usekden==1) then
758 [ # # ]: 0 : do ispden=1,dtset%nspden
759 [ # # ]: 0 : call fourdp(1,taumag(:,ispden),taur(:,ispden),+1,mpi_enreg,nfft,1,ngfft,tim_fourdp9)
760 : : end do
761 [ # # # # ]: 0 : taug(:,:)=reshape(taumag(:,1),(/2,nfft/))
762 : : end if
763 : : else
764 [ + + ]: 185298292 : do ifft=1,nfftmix
765 : 185287147 : jfft=mixtofft(ifft)
766 [ + + ]: 555872586 : rhog(1:2,jfft)=rhomag(2*ifft-1:2*ifft,1)
767 : : end do
768 : 11145 : call fourdp(1,rhog,rhor(:,1),+1,mpi_enreg,nfft,1,ngfft,tim_fourdp9)
769 [ + + ]: 11145 : if (dtset%nspden>1) then
770 [ + + ]: 6658 : do ispden=2,dtset%nspden
771 [ + + ]: 29848947 : do ifft=1,nfftmix
772 : 29844700 : jfft=mixtofft(ifft)
773 [ + + ]: 89538347 : magng(1:2,jfft,ispden-1)=rhomag(2*ifft-1:2*ifft,ispden)
774 : : end do
775 : 6658 : call fourdp(1,magng(:,:,ispden-1),rhor(:,ispden),+1,mpi_enreg,nfft,1,ngfft,tim_fourdp9)
776 : : end do
777 [ - + ]: 2411 : ABI_FREE(magng)
778 : : end if
779 [ + + ]: 11145 : if (dtset%usekden==1) then
780 [ + + ]: 7983418 : do ifft=1,nfftmix
781 : 7982995 : jfft=mixtofft(ifft)
782 [ + + ]: 23949408 : taug(1:2,jfft)=taumag(2*ifft-1:2*ifft,1)
783 : : end do
784 [ + - ]: 423 : call fourdp(1,taug,taur(:,1),+1,mpi_enreg,nfft,1,ngfft,tim_fourdp9)
785 [ - + ]: 423 : if (dtset%nspden>1) then
786 [ # # ]: 0 : do ispden=2,dtset%nspden
787 [ # # ]: 0 : do ifft=1,nfftmix
788 : 0 : jfft=mixtofft(ifft)
789 [ # # ]: 0 : magntaug(1:2,jfft,ispden-1)=taumag(2*ifft-1:2*ifft,ispden)
790 : : end do
791 [ # # ]: 0 : call fourdp(1,magntaug(:,:,ispden-1),taur(:,ispden),+1,mpi_enreg,nfft,1,ngfft,tim_fourdp9)
792 : : end do
793 [ # # ]: 0 : ABI_FREE(magntaug)
794 : : end if
795 : : end if
796 : : end if
797 [ - + ]: 13329 : ABI_FREE(rhomag)
798 [ - + ]: 13329 : ABI_FREE(taumag)
799 : :
800 : : !Set back rho in (up+dn,up) form if nspden=2
801 [ + + ]: 13329 : if (dtset%nspden==2) then
802 [ + + ]: 69282998 : rhor(:,2)=half*(rhor(:,1)+rhor(:,2))
803 [ - + - - ]: 1693 : if (dtset%usekden==1) taur(:,2)=half*(taur(:,1)+taur(:,2))
804 : : end if
805 : :
806 : : !In WVL: copy density to BigDFT object:
807 [ - + ]: 13329 : if(dtset%usewvl==1) then
808 : 0 : call wvl_rho_abi2big(1,rhor,wvl_den)
809 : : end if
810 : :
811 : 13329 : call timab(94,2,tsec)
812 : :
813 : : DBG_EXIT("COLL")
814 : :
815 [ + - + - : 66645 : end subroutine newrho
+ - - + -
+ - + ]
816 : : !!***
817 : :
818 : : end module m_newrho
819 : : !!***
|
