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 25844 : subroutine newrho(atindx,dbl_nnsclo,dielar,dielinv,dielstrt,dtn_pc,dtset,etotal,fcart,ffttomix,&
170 12922 : & gmet,grhf,gsqcut,initialized,ispmix,istep,kg_diel,kxc,mgfft,mix,mixtofft,&
171 12922 : & moved_atm_inside,mpi_enreg,my_natom,nattyp,nfft,&
172 : & nfftmix,nfftmix_per_nfft,ngfft,ngfftmix,nkxc,npawmix,npwdiel,&
173 12922 : & nresid,ntypat,n1xccc,pawrhoij,pawtab,&
174 12922 : & ph1d,precon,psps,rhog,rhor,rprimd,susmat,usepaw,vtrial,wvl,wvl_den,xred,rcpaw,extfpmd,&
175 38766 : & 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 12922 : real(dp),allocatable :: magng(:,:,:),magntaug(:,:,:)
235 12922 : real(dp),allocatable :: nresid0(:,:),nrespc(:,:),nreswk(:,:,:)
236 12922 : 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 12922 : real(dp), pointer :: rhomag(:,:), npaw(:)
240 12922 : real(dp),allocatable :: tauresid0(:,:),taurespc(:,:)
241 12922 : real(dp),allocatable :: taumag(:,:)
242 12922 : real(dp), pointer :: rcpaw_arr_(:)
243 12922 : real(dp),allocatable, target :: rcpaw_arr(:)
244 :
245 : ! *************************************************************************
246 :
247 : DBG_ENTER("COLL")
248 12922 : call timab(94,1,tsec)
249 :
250 51688 : nfftot=PRODUCT(ngfft(1:3))
251 :
252 : !Compatibility tests
253 12922 : if(nfftmix>nfft) then
254 0 : message='nfftmix>nfft not allowed!'
255 0 : ABI_BUG(message)
256 : end if
257 :
258 12922 : 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 12922 : 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 12922 : 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 12922 : if (usepaw==1.and.my_natom>0) then
288 11261 : cplex=pawrhoij(1)%cplex_rhoij;dplex=cplex-1
289 11261 : 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 12922 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
296 :
297 12922 : if(dtset%usewvl==0) then
298 12922 : 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 51688 : ABI_MALLOC(rhomag,(ispmix*nfftmix,dtset%nspden))
307 38766 : ABI_MALLOC(nresid0,(ispmix*nfftmix,dtset%nspden))
308 51688 : ABI_MALLOC(taumag,(ispmix*nfftmix,dtset%nspden*dtset%usekden))
309 38766 : ABI_MALLOC(tauresid0,(ispmix*nfftmix,dtset%nspden*dtset%usekden))
310 : ! real space and all fft points are here
311 12922 : if (ispmix==1.and.nfft==nfftmix) then
312 45997983 : rhomag(:,1:dtset%nspden)=rhor(:,1:dtset%nspden)
313 45997983 : nresid0(:,1:dtset%nspden)=nresid(:,1:dtset%nspden)
314 1404 : 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 11518 : 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 10834 : fact=dielar(4)-1._dp
343 43336 : ABI_MALLOC(nreswk,(2,nfft,dtset%nspden))
344 25876 : do ispden=1,dtset%nspden
345 25876 : call fourdp(1,nreswk(:,:,ispden),nresid(:,ispden),-1,mpi_enreg,nfft,1,ngfft,tim_fourdp9)
346 : end do
347 532220152 : do ifft=1,nfft
348 532220152 : if (ffttomix(ifft)>0) then
349 182919307 : jfft=2*ffttomix(ifft)
350 548757921 : rhomag (jfft-1:jfft,1)=rhog(1:2,ifft)
351 548757921 : nresid0(jfft-1:jfft,1)=nreswk(1:2,ifft,1)
352 : else
353 1047870033 : rhog(:,ifft)=rhog(:,ifft)+fact*nreswk(:,ifft,1)
354 : end if
355 : end do
356 10834 : if (dtset%nspden>1) then
357 9488 : ABI_MALLOC(magng,(2,nfft,dtset%nspden-1))
358 6580 : do ispden=2,dtset%nspden
359 4208 : call fourdp(1,magng(:,:,ispden-1),rhor(:,ispden),-1,mpi_enreg,nfft,1,ngfft,tim_fourdp9)
360 98244001 : do ifft=1,nfft
361 98241629 : if (ffttomix(ifft)>0) then
362 28566748 : jfft=2*ffttomix(ifft)
363 85700244 : rhomag (jfft-1:jfft,ispden)=magng (1:2,ifft,ispden-1)
364 85700244 : nresid0(jfft-1:jfft,ispden)=nreswk(1:2,ifft,ispden)
365 : else
366 209012019 : magng(:,ifft,ispden-1)=magng(:,ifft,ispden-1)+fact*nreswk(:,ifft,ispden)
367 158629119 : 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 10834 : 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 10834 : ABI_FREE(nreswk)
403 : end if
404 :
405 : !Retrieve "input" density from "output" density and density residual
406 489032669 : rhomag(:,1:dtset%nspden)=rhomag(:,1:dtset%nspden)-nresid0(:,1:dtset%nspden)
407 12922 : 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 12922 : if (dtset%nspden==2) then
413 38484460 : rhomag (:,2)=two*rhomag (:,2)-rhomag (:,1)
414 38484460 : nresid0(:,2)=two*nresid0(:,2)-nresid0(:,1)
415 1654 : 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 12922 : if (usepaw==1.and.my_natom>0) then
423 11261 : if (pawrhoij(1)%nspden==2) then
424 5474 : do iatom=1,my_natom
425 9362 : do iq=1,qphase
426 3888 : iq0=merge(0,cplex*pawrhoij(iatom)%lmn2_size,iq==1)
427 3888 : jrhoij=1+iq0
428 230124 : do irhoij=1,pawrhoij(iatom)%nrhoijsel
429 453069 : ro(1:1+dplex)=pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+dplex,1)
430 453069 : pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+dplex,1)=ro(1:1+dplex)+pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+dplex,2)
431 453069 : pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+dplex,2)=ro(1:1+dplex)-pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+dplex,2)
432 230124 : jrhoij=jrhoij+cplex
433 : end do
434 485811 : do kmix=1,pawrhoij(iatom)%lmnmix_sz
435 478035 : klmn=iq0+cplex*pawrhoij(iatom)%kpawmix(kmix)-dplex
436 956718 : ro(1:1+dplex)=pawrhoij(iatom)%rhoijres(klmn:klmn+dplex,1)
437 956718 : pawrhoij(iatom)%rhoijres(klmn:klmn+dplex,1)=ro(1:1+dplex)+pawrhoij(iatom)%rhoijres(klmn:klmn+dplex,2)
438 960606 : 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 38766 : ABI_MALLOC(nrespc,(ispmix*nfftmix,dtset%nspden))
447 38766 : ABI_MALLOC(taurespc,(ispmix*nfftmix,dtset%nspden*dtset%usekden))
448 38766 : ABI_MALLOC(npaw,(npawmix*usepaw))
449 12922 : if (usepaw==1) then
450 38541 : ABI_MALLOC(rhoijrespc,(npawmix))
451 : else
452 75 : ABI_MALLOC(rhoijrespc,(0))
453 : end if
454 12922 : 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 12922 : & 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 12922 : 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 12922 : if (mix%n_fftgr>0) then
478 12922 : i_vresid1=mix%i_vresid(1)
479 12922 : i_vrespc1=mix%i_vrespc(1)
480 : end if
481 :
482 12922 : if(associated(extfpmd)) then
483 83 : mix%useextfpmd=1
484 : endif
485 12922 : 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 12922 : if (moved_atm_inside == 1) then
496 0 : call abi_mixing_use_moving_atoms(mix, dtset%natom, xred, dtn_pc)
497 : end if
498 12922 : call abi_mixing_eval_allocate(mix, istep)
499 :
500 : !Copy current step arrays.
501 12922 : 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 12922 : call abi_mixing_copy_current_step(mix, nresid0, errid, message, arr_respc = nrespc)
505 : end if
506 12922 : if (errid /= AB7_NO_ERROR) then
507 0 : ABI_ERROR(message)
508 : end if
509 :
510 : !Same treatment for the kinetic energy density
511 12922 : 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 12922 : ABI_FREE(nresid0)
520 12922 : ABI_FREE(nrespc)
521 12922 : ABI_FREE(tauresid0)
522 12922 : ABI_FREE(taurespc)
523 :
524 : !PAW: either use the array f_paw or the array f_paw_disk
525 12922 : if (usepaw==1) then
526 12847 : indx=-dplex
527 42743 : do iatom=1,my_natom
528 89688 : ABI_MALLOC(rhoijtmp,(cplex*pawrhoij(iatom)%lmn2_size,1))
529 59792 : do iq=1,qphase
530 29896 : iq0=merge(0,cplex*pawrhoij(iatom)%lmn2_size,iq==1)
531 101145 : do ispden=1,pawrhoij(iatom)%nspden
532 3968630 : rhoijtmp=zero ; jrhoij=1+iq0
533 2138064 : do irhoij=1,pawrhoij(iatom)%nrhoijsel
534 2096711 : klmn=cplex*pawrhoij(iatom)%rhoijselect(irhoij)-dplex
535 5046187 : rhoijtmp(klmn:klmn+dplex,1)=pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+dplex,ispden)
536 2138064 : jrhoij=jrhoij+cplex
537 : end do
538 3048237 : do kmix=1,pawrhoij(iatom)%lmnmix_sz
539 2976988 : indx=indx+cplex;klmn=cplex*pawrhoij(iatom)%kpawmix(kmix)-dplex ; kklmn=klmn+iq0
540 6860752 : npaw(indx:indx+dplex)=rhoijtmp(klmn:klmn+dplex,1)-pawrhoij(iatom)%rhoijres(kklmn:kklmn+dplex,ispden)
541 6860752 : mix%f_paw(indx:indx+dplex,i_vresid1)=pawrhoij(iatom)%rhoijres(kklmn:kklmn+dplex,ispden)
542 6902105 : mix%f_paw(indx:indx+dplex,i_vrespc1)=rhoijrespc(indx:indx+dplex)
543 : end do
544 : end do
545 : end do
546 42743 : ABI_FREE(rhoijtmp)
547 : end do
548 : end if
549 :
550 12922 : rcpaw_arr_=>dum0
551 12922 : 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 12922 : 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 12922 : if(mpi_enreg%paral_kgb==1)then
578 4408 : mpicomm=mpi_enreg%comm_fft
579 4408 : mpi_summarize=.true.
580 : else
581 8514 : mpicomm=0
582 8514 : mpi_summarize=.false.
583 : end if
584 12922 : 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 12922 : reset = .false.
590 12922 : 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 12922 : & comm_atom=mpi_enreg%comm_atom)
601 12922 : if (errid == AB7_ERROR_MIXING_INC_NNSLOOP) then
602 0 : dbl_nnsclo = 1
603 12922 : else if (errid /= AB7_NO_ERROR) then
604 0 : ABI_ERROR(message)
605 : end if
606 : !Kinetic energy density mixing (if any)
607 12922 : 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 12922 : 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 12922 : nullify(rcpaw_arr_)
632 :
633 12922 : 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 12922 : 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 12922 : ABI_FREE(rhoijrespc)
688 : !end if
689 :
690 : !PAW: restore rhoij from compact storage
691 12922 : if (usepaw==1.and.dtset%iscf/=15.and.dtset%iscf/=16) then
692 12847 : indx=-dplex
693 42743 : do iatom=1,my_natom
694 119584 : ABI_MALLOC(rhoijtmp,(cplex*qphase*pawrhoij(iatom)%lmn2_size,pawrhoij(iatom)%nspden))
695 3957173 : rhoijtmp=zero
696 59792 : do iq=1,qphase
697 29896 : iq0=merge(0,cplex*pawrhoij(iatom)%lmn2_size,iq==1)
698 29896 : 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 71249 : do ispden=1,pawrhoij(iatom)%nspden
709 3048237 : do kmix=1,pawrhoij(iatom)%lmnmix_sz
710 2976988 : indx=indx+cplex;klmn=iq0+cplex*pawrhoij(iatom)%kpawmix(kmix)-dplex
711 6902105 : rhoijtmp(klmn:klmn+dplex,ispden)=npaw(indx:indx+dplex)
712 : end do
713 : end do
714 59792 : if (pawrhoij(iatom)%nspden==2) then
715 3888 : jrhoij=iq0+1
716 230124 : do irhoij=1,pawrhoij(iatom)%nrhoijsel
717 226236 : klmn=iq0+cplex*pawrhoij(iatom)%rhoijselect(irhoij)-dplex
718 453069 : ro(1:1+dplex)=rhoijtmp(klmn:klmn+dplex,1)
719 453069 : rhoijtmp(klmn:klmn+dplex,1)=half*(ro(1:1+dplex)+rhoijtmp(klmn:klmn+dplex,2))
720 453069 : rhoijtmp(klmn:klmn+dplex,2)=half*(ro(1:1+dplex)-rhoijtmp(klmn:klmn+dplex,2))
721 230124 : 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 29896 : & pawrhoij(iatom)%lmn2_size,pawrhoij(iatom)%nspden,rhoij_input=rhoijtmp)
728 42743 : ABI_FREE(rhoijtmp)
729 : end do
730 : end if ! usepaw==1.and.dtset%iscf/=15.and.dtset%iscf/=16
731 12922 : ABI_FREE(npaw)
732 :
733 : !Eventually write the data on disk and deallocate f_fftgr_disk
734 12922 : call abi_mixing_eval_deallocate(mix)
735 12922 : if (dtset%usekden==1) call abi_mixing_eval_deallocate(mix_mgga)
736 :
737 : !Fourier transform the density
738 12922 : if (ispmix==1.and.nfft==nfftmix) then
739 : !Real space mixing, no need to transform rhomag
740 45997983 : rhor(:,1:dtset%nspden)=rhomag(:,1:dtset%nspden)
741 1404 : if(dtset%usewvl==0) then
742 : !Get rhog from rhor(:,1)
743 1404 : call fourdp(1,rhog,rhor(:,1),-1,mpi_enreg,nfft,1,ngfft,tim_fourdp9)
744 : end if
745 1404 : 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 11518 : 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 182930141 : do ifft=1,nfftmix
765 182919307 : jfft=mixtofft(ifft)
766 548768755 : rhog(1:2,jfft)=rhomag(2*ifft-1:2*ifft,1)
767 : end do
768 10834 : call fourdp(1,rhog,rhor(:,1),+1,mpi_enreg,nfft,1,ngfft,tim_fourdp9)
769 10834 : if (dtset%nspden>1) then
770 6580 : do ispden=2,dtset%nspden
771 28570956 : do ifft=1,nfftmix
772 28566748 : jfft=mixtofft(ifft)
773 85704452 : magng(1:2,jfft,ispden-1)=rhomag(2*ifft-1:2*ifft,ispden)
774 : end do
775 6580 : call fourdp(1,magng(:,:,ispden-1),rhor(:,ispden),+1,mpi_enreg,nfft,1,ngfft,tim_fourdp9)
776 : end do
777 2372 : ABI_FREE(magng)
778 : end if
779 10834 : 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 12922 : ABI_FREE(rhomag)
798 12922 : ABI_FREE(taumag)
799 :
800 : !Set back rho in (up+dn,up) form if nspden=2
801 12922 : if (dtset%nspden==2) then
802 64969871 : rhor(:,2)=half*(rhor(:,1)+rhor(:,2))
803 1654 : if (dtset%usekden==1) taur(:,2)=half*(taur(:,1)+taur(:,2))
804 : end if
805 :
806 : !In WVL: copy density to BigDFT object:
807 12922 : if(dtset%usewvl==1) then
808 0 : call wvl_rho_abi2big(1,rhor,wvl_den)
809 : end if
810 :
811 12922 : call timab(94,2,tsec)
812 :
813 : DBG_EXIT("COLL")
814 :
815 64610 : end subroutine newrho
816 : !!***
817 :
818 : end module m_newrho
819 : !!***
|
