Line data Source code
1 : !!****m* ABINIT/m_newvtr
2 : !! NAME
3 : !! m_newvtr
4 : !!
5 : !! FUNCTION
6 : !!
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 1998-2026 ABINIT group (DCA, XG, 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_newvtr
23 :
24 : use defs_basis
25 : use defs_wvltypes
26 : use m_abicore
27 : use m_errors
28 : use m_abi2big
29 : use m_abi_mixing
30 : use m_cgtools
31 : use m_dtset
32 : use m_precon
33 :
34 : use defs_datatypes, only : pseudopotential_type
35 : use defs_abitypes, only : MPI_type
36 : use m_time, only : timab
37 : use m_geometry, only : metric
38 : use m_pawtab, only : pawtab_type
39 : use m_pawrhoij, only : pawrhoij_type,pawrhoij_filter
40 : use m_prcref, only : prcref_PMA
41 : use m_wvl_rho, only : wvl_prcref
42 : use m_fft, only : fourdp
43 : use m_xctk, only : xcpot
44 :
45 : implicit none
46 :
47 : private
48 : !!***
49 :
50 : public :: newvtr
51 : !!***
52 :
53 : contains
54 : !!***
55 :
56 : !!****f* ABINIT/newvtr
57 : !! NAME
58 : !! newvtr
59 : !!
60 : !! FUNCTION
61 : !! Compute new trial potential by mixing new and old values.
62 : !! Call prcref to compute preconditioned residual potential and forces,
63 : !! Then, call one of the self-consistency drivers,
64 : !! then update vtrial.
65 : !!
66 : !! INPUTS
67 : !! atindx(natom)=index table for atoms (see gstate.f)
68 : !! dielar(7)=input parameters for dielectric matrix:
69 : !! diecut,dielng,diemac,diemix,diegap,dielam,diemixmag.
70 : !! dielinv(2,npwdiel,nspden,npwdiel,nspden)=
71 : !! inverse of the dielectric matrix in rec. space
72 : !! dielstrt=number of the step at which the dielectric preconditioning begins.
73 : !! dtset <type(dataset_type)>=all input variables in this dataset
74 : !! | spinmagntarget=input variable that governs fixed moment calculation
75 : !! | intxc=control xc quadrature
76 : !! | densfor_pred= governs the preconditioning of the atomic charges
77 : !! | iprcel= governs the preconditioning of the potential residual
78 : !! | iprcfc=governs the preconditioning of the forces
79 : !! | iscf=( <= 0 =>non-SCF), >0 => SCF)
80 : !! | iscf =1 => determination of the largest eigenvalue of the SCF cycle
81 : !! | iscf =2 => SCF cycle, simple mixing
82 : !! | iscf =3 => SCF cycle, Anderson mixing
83 : !! | iscf =4 => SCF cycle, Anderson mixing (order 2)
84 : !! | iscf =5 => SCF cycle, CG based on the minimization of the energy
85 : !! | iscf =7 => SCF cycle, Pulay mixing
86 : !! | isecur=level of security of the computation
87 : !! | ixc=exchange-correlation choice parameter.
88 : !! | mffmem=governs the number of FFT arrays which are fit in core memory
89 : !! | it is either 1, in which case the array f_fftgr is used,
90 : !! | or 0, in which case the array f_fftgr_disk is used
91 : !! | natom=number of atoms
92 : !! | nspden=number of spin-density components
93 : !! | occopt=option for occupancies
94 : !! | paral_kgb=option for (kpt,g vectors,bands) parallelism
95 : !! | pawoptmix= - PAW only - 1 if the computed residuals include the PAW (rhoij) part
96 : !! | prtvol=control print volume and debugging
97 : !! | typat(natom)=integer type for each atom in cell
98 : !! etotal=the total energy obtained from the input vtrial
99 : !! fcart(3,natom)=cartesian forces (hartree/bohr)
100 : !! ffttomix(nfft*(1-nfftmix/nfft))=Index of the points of the FFT (fine) grid on the grid used for mixing (coarse)
101 : !! gmet(3,3)=metrix tensor in G space in Bohr**-2.
102 : !! grhf(3,natom)=Hellman-Feynman derivatives of the total energy
103 : !! gsqcut=cutoff on (k+G)^2 (bohr^-2)
104 : !! initialized= if 0, the initialization of the gstate run is not yet finished
105 : !! ispmix=1 if mixing is done in real space, 2 if mixing is done in reciprocal space
106 : !! istep= number of the step in the SCF cycle
107 : !! kg_diel(3,npwdiel)=reduced planewave coordinates for the dielectric matrix.
108 : !! kxc(nfft,nkxc)=exchange-correlation kernel, needed only for electronic!
109 : ! dielectric matrix
110 : !! mgfft=maximum size of 1D FFTs
111 : !! mixtofft(nfftmix*(1-nfftmix/nfft))=Index of the points of the FFT grid used for mixing (coarse) on the FFT (fine) grid
112 : !! moved_atm_inside= if 1, then the preconditioned forces
113 : !! as well as the preconditioned potential residual must be computed;
114 : !! otherwise, compute only the preconditioned potential residual.
115 : !! mpi_enreg=information about MPI parallelization
116 : !! my_natom=number of atoms treated by current processor
117 : !! nattyp(ntypat)=number of atoms of each type in cell.
118 : !! nfft=(effective) number of FFT grid points (for this processor)
119 : !! nfftmix=dimension of FFT grid used to mix the densities (used in PAW only)
120 : !! ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
121 : !! ngfftmix(18)=contain all needed information about 3D FFT, for the grid corresponding to nfftmix
122 : !! nkxc=second dimension of the array kxc, see rhotoxc.f for a description
123 : !! npawmix=-PAW only- number of spherical part elements to be mixed
124 : !! npwdiel=number of planewaves for dielectric matrix
125 : !! nstep=number of steps expected in iterations.
126 : !! ntypat=number of types of atoms in cell.
127 : !! n1xccc=dimension of xccc1d ; 0 if no XC core correction is used
128 : !! pawrhoij(my_natom*usepaw) <type(pawrhoij_type)>= paw rhoij occupancies and related data
129 : !! Use here rhoij residuals (and gradients)
130 : !! psps <type(pseudopotential_type)>=variables related to pseudopotentials
131 : !! rhor(nfft,nspden)=array for electron density in electrons/bohr**3.
132 : !! rprimd(3,3)=dimensional primitive translations in real space (bohr)
133 : !! susmat(2,npwdiel,nspden,npwdiel,nspden)=
134 : !! the susceptibility (or density-density response) matrix in reciprocal space
135 : !! [vtauresid(nfft,nspden*usevtau)]=array for vtau residue (see vtau below))
136 : !! usepaw= 0 for non paw calculation; =1 for paw calculation
137 : !! vhartr(nfft)=array for holding Hartree potential
138 : !! vnew_mean(nspden)=constrained mean value of the future trial potential
139 : !! (might be spin-polarized)
140 : !! vres_mean(nspden)=mean value of the potential residual
141 : !! vpsp(nfft)=array for holding local psp
142 : !! vresid(nfft,nspden)=array for the residual of the potential
143 : !! vxc(nfft,nspden)=exchange-correlation potential (hartree)
144 : !! [vtau(nfftf,dtset%nspden,4*usevtau)]=derivative of XC energy density
145 : !! with respect to kinetic energy density (metaGGA cases) (optional)
146 : !! xred(3,natom)=reduced dimensionless atomic coordinates
147 : !!
148 : !! OUTPUT
149 : !! dbl_nnsclo=1 if nnsclo has to be doubled to secure the convergence.
150 : !!
151 : !! SIDE EFFECTS
152 : !! dtn_pc(3,natom)=preconditioned change of atomic position,
153 : !! in reduced coordinates
154 : !! vtrial(nfft,nspden)= at input, it is the "in" trial potential that gave vresid=(v_out-v_in)
155 : !! at output, it is an updated "mixed" trial potential
156 : !! ===== if size(vtau)>0, i.e. usevtau==1 =====
157 : !! [mix_mgga<type(abi_mixing_object)>]=all data defining the mixing algorithm for
158 : !! the kinetic energy potential
159 : !! ===== if densfor_pred==3 .and. moved_atm_inside==1 =====
160 : !! ph1d(2,3*(2*mgfft+1)*natom)=1-dim structure factor phases
161 : !! ==== if usepaw==1
162 : !! pawrhoij(natom)%nrhoijsel,rhoijselect,rhoijp= several arrays
163 : !! containing new values of rhoij (augmentation occupancies)
164 : !!
165 : !! WARNINGS
166 : !! depending on the value of densfor_pred and moved_atm_inside,
167 : !! the xc potential or the Hxc potential may have been subtracted from vtrial !
168 : !!
169 : !! NOTES
170 : !! In case of PAW calculations:
171 : !! Computations are done either on the fine FFT grid or the coarse grid (depending on dtset%pawmixdg)
172 : !! All variables (nfft,ngfft,mgfft) refer to the fine FFT grid.
173 : !! All arrays (densities/potentials...) are computed on this fine FFT grid.
174 : !! Developpers have to be careful when introducing others arrays:
175 : !! they have to be stored on the fine FFT grid.
176 : !! In case of norm-conserving calculations the FFT grid is the usual FFT grid.
177 : !!
178 : !! Subtility in PAW and non-collinear magnetism:
179 : !! Potentials are stored in (up-up,dn-dn,Re[up-dn],Im[up-dn]) format
180 : !! On-site occupancies (rhoij) are stored in (n,mx,my,mz)
181 : !! This is compatible provided that the mixing factors for n and m are identical
182 : !! and that the residual is not a combination of V_res and rhoij_res (pawoptmix=0).
183 : !!
184 : !! SOURCE
185 :
186 56560 : subroutine newvtr(atindx,dbl_nnsclo,dielar,dielinv,dielstrt,&
187 28280 : & dtn_pc,dtset,etotal,fcart,ffttomix,&
188 28280 : & gmet,grhf,gsqcut,&
189 : & initialized,ispmix,&
190 : & istep,&
191 28280 : & kg_diel,kxc,mgfft,mix,mixtofft,&
192 28280 : & moved_atm_inside,mpi_enreg,my_natom,nattyp,nfft,nfftmix,&
193 : & ngfft,ngfftmix,nkxc,npawmix,npwdiel,&
194 : & nstep,ntypat,n1xccc,&
195 28280 : & pawrhoij,&
196 28280 : & ph1d,&
197 : & precon,&
198 28280 : & psps,rhor,rprimd,susmat,usepaw,&
199 28280 : & vhartr,vnew_mean,vpsp,vresid,vres_mean,&
200 28280 : & vtrial,vxc,xred,&
201 : & nfftf,&
202 28280 : & pawtab,&
203 28280 : & rhog,&
204 : & wvl,&
205 28280 : & mix_mgga,vtau,vtauresid) ! Optional arguments
206 :
207 : !Arguments-------------------------------
208 : ! WARNING
209 : ! BEWARE THERE IS TWO DIFFERENT SIZE DECLARED FOR ARRAY NHAT IN RHOTOV AND RHOHXC
210 : ! THIS MIGHT RESULT IN A BUG
211 : !scalars
212 : integer,intent(in) :: dielstrt,initialized,ispmix,istep,mgfft
213 : integer,intent(in) :: moved_atm_inside,my_natom,n1xccc,nfft
214 : integer,intent(in) :: nfftf,nfftmix,nkxc,npawmix,npwdiel,nstep
215 : integer,intent(in) :: ntypat,usepaw
216 : integer,intent(inout) :: dbl_nnsclo
217 : real(dp),intent(in) :: etotal,gsqcut
218 : type(MPI_type),intent(in) :: mpi_enreg
219 : type(dataset_type),intent(in) :: dtset
220 : type(abi_mixing_object),intent(inout) :: mix
221 : type(abi_mixing_object),intent(inout),optional :: mix_mgga
222 : type(pseudopotential_type),intent(in) :: psps
223 : type(wvl_data), intent(inout) :: wvl
224 : !arrays
225 : integer,intent(in) :: atindx(dtset%natom)
226 : integer,intent(in) :: ffttomix(nfft*(1-nfftmix/nfft))
227 : integer,intent(in) :: kg_diel(3,npwdiel)
228 : integer,intent(in) :: mixtofft(nfftmix*(1-nfftmix/nfft)),nattyp(ntypat)
229 : integer,intent(in) :: ngfft(18),ngfftmix(18)
230 : real(dp),intent(in) :: dielar(7)
231 : real(dp),intent(in) :: fcart(3,dtset%natom),grhf(3,dtset%natom)
232 : real(dp),intent(inout) :: rprimd(3,3)
233 : real(dp),intent(in) :: susmat(2,npwdiel,dtset%nspden,npwdiel,dtset%nspden)
234 : real(dp),intent(in) :: vhartr(nfft),vnew_mean(dtset%nspden),vres_mean(dtset%nspden)
235 : real(dp),intent(in) :: vxc(nfft,dtset%nspden)
236 : real(dp),intent(inout) :: dielinv(2,npwdiel,dtset%nspden,npwdiel,dtset%nspden)
237 : real(dp),intent(inout), target :: dtn_pc(3,dtset%natom)
238 : real(dp),intent(inout) :: gmet(3,3)
239 : real(dp),intent(inout) :: kxc(nfft,nkxc),ph1d(2,3*(2*mgfft+1)*dtset%natom)
240 : real(dp),intent(inout) :: rhog(2,nfftf),vpsp(nfft)
241 : real(dp),intent(inout), target :: rhor(nfft,dtset%nspden)
242 : real(dp),intent(inout) :: vresid(nfft,dtset%nspden),vtrial(nfft,dtset%nspden)
243 : real(dp),intent(inout), target :: xred(3,dtset%natom)
244 : real(dp),intent(inout),optional :: vtau(:,:,:) !vtau(nfft,dtset%nspden,4*usevxctau)
245 : real(dp),intent(inout),optional :: vtauresid(:,:) !vtauresid(nfft,dtset%nspden*usevxctau)
246 : type(pawrhoij_type),intent(inout) :: pawrhoij(my_natom*usepaw)
247 : type(pawtab_type),intent(in) :: pawtab(ntypat*usepaw)
248 : type(precon_object), intent(inout) :: precon
249 :
250 : !Local variables-------------------------------
251 : !scalars
252 : integer :: cplex,dplex,i_vresid1,i_vrespc1
253 : ! integer :: i1,i2,i3,ifft2,ifft3,ifft4,ifft5,ii1,ii2,ii3,ii4,ii5
254 : integer :: errid,iatom,ifft,indx,iq,iq0,irhoij,ispden,jfft,jrhoij,klmn,kklmn,kmix
255 : integer :: mpicomm,mpi_comm_sphgrid,n1,n2,n3,nfftot,qphase,tim_fourdp,usevtau
256 : logical :: mpi_summarize,reset,with_vtau
257 : real(dp) :: dielng,diemix,fact,ucvol,ucvol_local,vme
258 : character(len=500) :: message
259 : !arrays
260 : real(dp),parameter :: identity(4)=(/one,one,zero,zero/)
261 56560 : real(dp) :: gprimd(3,3),rmet(3,3),tsec(2),vmean(dtset%nspden)
262 28280 : real(dp),allocatable :: rhoijrespc(:)
263 28280 : real(dp),allocatable :: rhoijtmp(:,:)
264 28280 : real(dp),allocatable :: vresid0(:,:),vrespc(:,:),vreswk(:,:),vtrialg(:,:,:)
265 28280 : real(dp),allocatable :: vtauresid0(:,:),vtaurespc(:,:),vtaug(:,:,:),vtau0(:,:)
266 28280 : real(dp),pointer :: vtrial0(:,:),vpaw(:)
267 :
268 : ! *************************************************************************
269 :
270 : !DEBUG
271 : !write(std_out,*)' newvtr : enter '
272 : !write(std_out,*)' newvtr : ispmix,nfft,nfftmix=',ispmix,nfft,nfftmix
273 : !ENDDEBUG
274 :
275 28280 : call timab(93,1,tsec)
276 28280 : call timab(901,1,tsec)
277 28280 : tim_fourdp=8
278 :
279 : !mpicomm over spherical grid:
280 28280 : mpi_comm_sphgrid=mpi_enreg%comm_fft
281 28280 : if(dtset%usewvl==1) mpi_comm_sphgrid=mpi_enreg%comm_wvl
282 :
283 : !Compatibility tests
284 28280 : if(nfftmix>nfft) then
285 0 : ABI_BUG(' nfftmix>nfft not allowed !')
286 : end if
287 :
288 28280 : if(ispmix/=2.and.nfftmix/=nfft) then
289 0 : message = ' nfftmix/=nfft allowed only when ispmix=2 !'
290 0 : ABI_BUG(message)
291 : end if
292 :
293 28280 : if(dtset%usewvl==1) then
294 0 : if(dtset%wvl_bigdft_comp==1) then
295 0 : message = 'newvtr: usewvl == 1 and wvl_bigdft_comp==1 not allowed (use wvl_newtr() instead)!'
296 0 : ABI_BUG(message)
297 : end if
298 0 : if(ispmix/=1 .or. nfftmix/=nfft) then
299 0 : ABI_BUG('newvtr: nfftmix/=nfft, ispmix/=1 not allowed for wavelets')
300 : end if
301 : end if
302 :
303 : !Test size of kinetic energy potential Vtau
304 28280 : with_vtau = (present(vtau).and.present(vtauresid))
305 113120 : with_vtau = (with_vtau.and.size(vtau)>0)
306 : if (with_vtau) then
307 834 : if (size(vtau)/=nfft*dtset%nspden*4.or.size(vtauresid)/=nfft*dtset%nspden) then
308 0 : ABI_BUG("Wrong size for vxctau/vtauresid!")
309 : end if
310 139 : if (.not.present(mix_mgga)) then
311 0 : message='mix_gga array missing!'
312 0 : ABI_BUG(message)
313 : end if
314 139 : if (mix_mgga%iscf==ABI_MIXING_CG_ENERGY.or.mix_mgga%iscf==ABI_MIXING_CG_ENERGY_2.or.&
315 : & mix_mgga%iscf==ABI_MIXING_EIG) then
316 0 : message='kinetic energy potential cannot be mixed with the selected mixing algorithm!'
317 0 : ABI_ERROR(message)
318 : end if
319 : end if
320 28141 : usevtau=0;if (with_vtau) usevtau=1
321 :
322 28280 : if(usepaw==1.and.dtset%nspden==4.and.dtset%pawoptmix==1) then
323 0 : message = ' pawoptmix=1 is not compatible with nspden=4 !'
324 0 : ABI_ERROR(message)
325 : end if
326 :
327 28280 : dielng=dielar(2)
328 28280 : diemix=dielar(4)
329 28280 : n1=ngfft(1)
330 28280 : n2=ngfft(2)
331 28280 : n3=ngfft(3)
332 28280 : if (usepaw==1.and.my_natom>0) then
333 782 : cplex=pawrhoij(1)%cplex_rhoij;dplex=cplex-1
334 782 : qphase=pawrhoij(1)%qphase
335 : else
336 : cplex=0;dplex=0 ; qphase=0
337 : end if
338 :
339 : !Get size of FFT grid
340 113120 : nfftot=PRODUCT(ngfft(1:3))
341 :
342 : !Compute different geometric tensor, as well as ucvol, from rprimd
343 28280 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
344 :
345 28280 : if(dtset%usewvl==0) then
346 28280 : ucvol_local=ucvol
347 : #if defined HAVE_BIGDFT
348 : else
349 : ucvol_local = product(wvl%den%denspot%dpbox%hgrids) * real(nfftot, dp)
350 : #endif
351 : end if
352 :
353 : !------Treat the mean of potentiel residual
354 :
355 : !Special care must be taken with components of the
356 : !potential that are associated with NO density change.
357 : !In general, only the global mean of the potential has
358 : !such an anomalous feature. However, in the spin
359 : !polarized case with fixed occupancies, also the
360 : !mean of each spin-potential (independently of the other)
361 : !has such a behaviour. The trick is to remove these
362 : !variables before going in the predictive routines,
363 : !then to put them back
364 :
365 : !Compute the mean of the old vtrial
366 28280 : call mean_fftr(vtrial,vmean,nfft,nfftot,dtset%nspden,mpi_comm_sphgrid)
367 :
368 : !When (collinear) spin-polarized and fixed occupation numbers,
369 : !treat separately spin up and spin down.
370 : !Otherwise, use only global mean
371 63576 : do ispden=1,dtset%nspden
372 35296 : if (dtset%nspden==2.and.dtset%occopt>=3.and. &
373 : & abs(dtset%spinmagntarget+99.99_dp)<1.0d-10)then
374 1340 : vme=(vmean(1)+vmean(2))*half
375 : else
376 33956 : vme=vmean(ispden)
377 : end if
378 577586908 : vtrial(:,ispden)=vtrial(:,ispden)-vme
379 : end do
380 :
381 28280 : call timab(901,2,tsec)
382 :
383 28280 : call timab(902,1,tsec)
384 :
385 : !Select components of potential to be mixed
386 113120 : ABI_MALLOC(vtrial0,(ispmix*nfftmix,dtset%nspden))
387 84840 : ABI_MALLOC(vresid0,(ispmix*nfftmix,dtset%nspden))
388 113120 : ABI_MALLOC(vtau0,(ispmix*nfftmix,dtset%nspden*usevtau))
389 84840 : ABI_MALLOC(vtauresid0,(ispmix*nfftmix,dtset%nspden*usevtau))
390 28280 : if (ispmix==1.and.nfft==nfftmix) then
391 1117817968 : vtrial0=vtrial;vresid0=vresid
392 27714 : if (with_vtau) then
393 15253282 : vtau0(:,:)=vtau(:,:,1);vtauresid0=vtauresid
394 : end if
395 566 : else if (nfft==nfftmix) then
396 0 : do ispden=1,dtset%nspden
397 0 : call fourdp(1,vtrial0(:,ispden),vtrial(:,ispden),-1,mpi_enreg,nfft,1,ngfft,tim_fourdp)
398 0 : call fourdp(1,vresid0(:,ispden),vresid(:,ispden),-1,mpi_enreg,nfft,1,ngfft,tim_fourdp)
399 : end do
400 0 : if (with_vtau) then
401 0 : do ispden=1,dtset%nspden
402 0 : call fourdp(1,vtau0(:,ispden),vtau(:,ispden,1),-1,mpi_enreg,nfft,1,ngfft,tim_fourdp)
403 0 : call fourdp(1,vtauresid0(:,ispden),vtauresid(:,ispden),-1,mpi_enreg,nfft,1,ngfft,tim_fourdp)
404 : end do
405 : end if
406 : else
407 2264 : ABI_MALLOC(vtrialg,(2,nfft,dtset%nspden))
408 1698 : ABI_MALLOC(vreswk,(2,nfft))
409 1484 : do ispden=1,dtset%nspden
410 918 : fact=dielar(4);if (ispden>1) fact=dielar(7)
411 918 : call fourdp(1,vtrialg(:,:,ispden),vtrial(:,ispden),-1,mpi_enreg,nfft,1,ngfft,tim_fourdp)
412 918 : call fourdp(1,vreswk,vresid(:,ispden),-1,mpi_enreg,nfft,1,ngfft,tim_fourdp)
413 18677924 : do ifft=1,nfft
414 18677358 : if (ffttomix(ifft)>0) then
415 3578860 : jfft=2*ffttomix(ifft)
416 3578860 : vtrial0(jfft-1,ispden)=vtrialg(1,ifft,ispden)
417 3578860 : vtrial0(jfft ,ispden)=vtrialg(2,ifft,ispden)
418 3578860 : vresid0(jfft-1,ispden)=vreswk(1,ifft)
419 3578860 : vresid0(jfft ,ispden)=vreswk(2,ifft)
420 : else
421 45292740 : vtrialg(:,ifft,ispden)=vtrialg(:,ifft,ispden)+fact*vreswk(:,ifft)
422 : end if
423 : end do
424 : end do
425 566 : if (with_vtau) then
426 75 : ABI_MALLOC(vtaug,(2,nfft,dtset%nspden))
427 50 : do ispden=1,dtset%nspden
428 25 : fact=dielar(4);if (ispden>1) fact=dielar(7)
429 25 : call fourdp(1,vtaug(:,:,ispden),vtau(:,ispden,1),-1,mpi_enreg,nfft,1,ngfft,tim_fourdp)
430 25 : call fourdp(1,vreswk,vtauresid(:,ispden),-1,mpi_enreg,nfft,1,ngfft,tim_fourdp)
431 345650 : do ifft=1,nfft
432 345625 : if (ffttomix(ifft)>0) then
433 102400 : jfft=2*ffttomix(ifft)
434 102400 : vtau0(jfft-1,ispden)=vtaug(1,ifft,ispden)
435 102400 : vtau0(jfft ,ispden)=vtaug(2,ifft,ispden)
436 102400 : vtauresid0(jfft-1,ispden)=vreswk(1,ifft)
437 102400 : vtauresid0(jfft ,ispden)=vreswk(2,ifft)
438 : else
439 729600 : vtaug(:,ifft,ispden)=vtaug(:,ifft,ispden)+fact*vreswk(:,ifft)
440 : end if
441 : end do
442 : end do
443 : end if
444 566 : ABI_FREE(vreswk)
445 : end if
446 :
447 : !Retrieve "input" Vtau from "output" one and potential residual
448 680 : if (with_vtau) then
449 7831491 : vtau0(:,1:dtset%nspden)=vtau0(:,1:dtset%nspden)-vtauresid0(:,1:dtset%nspden)
450 : end if
451 :
452 : !Choice of preconditioner governed by iprcel, densfor_pred and iprcfc
453 84840 : ABI_MALLOC(vrespc,(ispmix*nfftmix,dtset%nspden))
454 84840 : ABI_MALLOC(vtaurespc,(ispmix*nfftmix,dtset%nspden*usevtau))
455 84840 : ABI_MALLOC(vpaw,(npawmix*usepaw))
456 28280 : if (usepaw==1) then
457 2346 : ABI_MALLOC(rhoijrespc,(npawmix))
458 : else
459 27498 : ABI_MALLOC(rhoijrespc,(0))
460 : end if
461 :
462 28280 : call timab(902,2,tsec)
463 28280 : call timab(903,1,tsec)
464 :
465 28280 : if(dtset%usewvl==0) then
466 : call prcref_PMA(atindx,dielar,dielinv,dielstrt,dtn_pc,dtset,fcart,ffttomix,gmet,gsqcut,&
467 : & istep,kg_diel,kxc,mgfft,moved_atm_inside,mpi_enreg,my_natom,&
468 : & nattyp,nfft,nfftmix,ngfft,ngfftmix,nkxc,npawmix,npwdiel,ntypat,n1xccc,&
469 : & ispmix,0,pawrhoij,ph1d, precon,psps,rhog,rhoijrespc,rhor,rprimd,susmat,&
470 : & vhartr,vpsp,vresid0,vrespc,vxc,xred,&
471 28280 : & etotal,pawtab,wvl)
472 : else
473 : call wvl_prcref(dielar,dtset%iprcel,my_natom,nfftmix,npawmix,dtset%nspden,pawrhoij,&
474 0 : & rhoijrespc,psps%usepaw,vresid0,vrespc)
475 : end if
476 : !At present, only a simple precoditionning for vtau
477 : ! (is Kerker mixing valid for vtau?)
478 28280 : if (with_vtau) then
479 299 : do ispden=1,dtset%nspden
480 160 : fact=dielar(4);if (ispden>1) fact=abs(dielar(7))
481 7831491 : vtaurespc(1:ispmix*nfftmix,ispden)=fact*vtauresid0(1:ispmix*nfftmix,ispden)
482 : end do
483 : end if
484 :
485 28280 : call timab(903,2,tsec)
486 28280 : call timab(904,1,tsec)
487 :
488 : !------Compute new vtrial and eventual new atomic positions
489 :
490 28280 : if (mix%n_fftgr>0) then
491 28280 : i_vresid1=mix%i_vresid(1)
492 28280 : i_vrespc1=mix%i_vrespc(1)
493 : end if
494 :
495 : !Initialise working arrays for the mixing object.
496 28280 : if (moved_atm_inside == 1) then
497 110 : call abi_mixing_use_moving_atoms(mix, dtset%natom, xred, dtn_pc)
498 : end if
499 28280 : call abi_mixing_eval_allocate(mix, istep)
500 : !Copy current step arrays.
501 28280 : if (moved_atm_inside == 1) then
502 : call abi_mixing_copy_current_step(mix, vresid0, errid, message, &
503 110 : & arr_respc = vrespc, arr_atm = grhf)
504 : else
505 : call abi_mixing_copy_current_step(mix, vresid0, errid, message, &
506 28170 : & arr_respc = vrespc)
507 : end if
508 28280 : if (errid /= AB7_NO_ERROR) then
509 0 : ABI_ERROR(message)
510 : end if
511 28280 : if (with_vtau) then
512 139 : call abi_mixing_eval_allocate(mix_mgga, istep)
513 : call abi_mixing_copy_current_step(mix_mgga, vtauresid0, errid, message, &
514 139 : & arr_respc = vtaurespc)
515 139 : if (errid /= AB7_NO_ERROR) then
516 0 : ABI_ERROR(message)
517 : end if
518 : end if
519 28280 : ABI_FREE(vresid0)
520 28280 : ABI_FREE(vrespc)
521 28280 : ABI_FREE(vtauresid0)
522 28280 : ABI_FREE(vtaurespc)
523 :
524 : !PAW: either use the array f_paw or the array f_paw_disk
525 28280 : if (usepaw==1) then
526 782 : indx=-dplex
527 1921 : do iatom=1,my_natom
528 3417 : ABI_MALLOC(rhoijtmp,(cplex*pawrhoij(iatom)%lmn2_size,1))
529 2278 : do iq=1,qphase
530 1139 : iq0=merge(0,cplex*pawrhoij(iatom)%lmn2_size,iq==1)
531 4032 : do ispden=1,pawrhoij(iatom)%nspden
532 243672 : rhoijtmp=zero ; jrhoij=iq0+1
533 69504 : do irhoij=1,pawrhoij(iatom)%nrhoijsel
534 67750 : klmn=cplex*pawrhoij(iatom)%rhoijselect(irhoij)-dplex
535 135500 : rhoijtmp(klmn:klmn+dplex,1)=pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+dplex,ispden)
536 69504 : jrhoij=jrhoij+cplex
537 : end do
538 243057 : do kmix=1,pawrhoij(iatom)%lmnmix_sz
539 240164 : indx=indx+cplex;klmn=cplex*pawrhoij(iatom)%kpawmix(kmix)-dplex; kklmn=iq0+klmn
540 480328 : vpaw(indx:indx+dplex)=rhoijtmp(klmn:klmn+dplex,1)-pawrhoij(iatom)%rhoijres(kklmn:kklmn+dplex,ispden)
541 480328 : mix%f_paw(indx:indx+dplex,i_vresid1)=pawrhoij(iatom)%rhoijres(kklmn:kklmn+dplex,ispden)
542 482082 : mix%f_paw(indx:indx+dplex,i_vrespc1)=rhoijrespc(indx:indx+dplex)
543 : end do
544 : end do
545 : end do
546 1921 : ABI_FREE(rhoijtmp)
547 : end do
548 : end if
549 :
550 : !------Prediction of the components of the potential associated with a density change
551 :
552 : !Init mpicomm
553 28280 : if(mpi_enreg%paral_kgb==1)then
554 2372 : mpicomm=mpi_enreg%comm_fft
555 2372 : mpi_summarize=.true.
556 : else
557 25908 : mpicomm=0
558 25908 : mpi_summarize=.false.
559 : end if
560 :
561 28280 : reset = .false.
562 28280 : if (initialized == 0) reset = .true.
563 : call abi_mixing_eval(mix, vtrial0, istep, nfftot, ucvol_local, &
564 : & mpicomm, mpi_summarize, errid, message, &
565 : & reset = reset, isecur = dtset%isecur, &
566 : & pawopt = dtset%pawoptmix, pawarr = vpaw, etotal = etotal, potden = rhor, &
567 28280 : & comm_atom=mpi_enreg%comm_atom)
568 28280 : if (errid == AB7_ERROR_MIXING_INC_NNSLOOP) then
569 0 : dbl_nnsclo = 1
570 28280 : else if (errid /= AB7_NO_ERROR) then
571 0 : ABI_ERROR(message)
572 : end if
573 28280 : if (with_vtau) then
574 : call abi_mixing_eval(mix_mgga, vtau0, istep, nfftot, ucvol_local, &
575 139 : & mpicomm, mpi_summarize, errid, message, reset = reset)
576 139 : if (errid /= AB7_NO_ERROR) then
577 0 : ABI_ERROR(message)
578 : end if
579 : end if
580 :
581 : !PAW: apply a simple mixing to rhoij (this is temporary)
582 28280 : if(dtset%iscf==5 .or. dtset%iscf==6)then
583 174 : if (usepaw==1) then
584 0 : indx=-dplex
585 0 : do iatom=1,my_natom
586 0 : ABI_MALLOC(rhoijtmp,(cplex*qphase*pawrhoij(iatom)%lmn2_size,pawrhoij(iatom)%nspden))
587 0 : rhoijtmp=zero
588 0 : do iq=1,qphase
589 0 : iq0=merge(0,cplex*pawrhoij(iatom)%lmn2_size,iq==1)
590 0 : if (pawrhoij(iatom)%lmnmix_sz<pawrhoij(iatom)%lmn2_size) then
591 0 : do ispden=1,pawrhoij(iatom)%nspden
592 0 : do kmix=1,pawrhoij(iatom)%lmnmix_sz
593 0 : indx=indx+cplex;klmn=iq0+cplex*pawrhoij(iatom)%kpawmix(kmix)-dplex
594 : rhoijtmp(klmn:klmn+dplex,ispden)=rhoijrespc(indx:indx+dplex) &
595 0 : & -pawrhoij(iatom)%rhoijres(klmn:klmn+dplex,ispden)
596 : end do
597 : end do
598 : end if
599 0 : do ispden=1,pawrhoij(iatom)%nspden
600 0 : jrhoij=iq0+1
601 0 : do irhoij=1,pawrhoij(iatom)%nrhoijsel
602 0 : klmn=iq0+cplex*pawrhoij(iatom)%rhoijselect(irhoij)-dplex
603 : rhoijtmp(klmn:klmn+dplex,ispden)=rhoijtmp(klmn:klmn+dplex,ispden) &
604 0 : & +pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+dplex,ispden)
605 0 : jrhoij=jrhoij+cplex
606 : end do
607 : end do
608 : end do
609 : call pawrhoij_filter(pawrhoij(iatom)%rhoijp,pawrhoij(iatom)%rhoijselect,&
610 : & pawrhoij(iatom)%nrhoijsel,pawrhoij(iatom)%cplex_rhoij,pawrhoij(iatom)%qphase,&
611 0 : & pawrhoij(iatom)%lmn2_size,pawrhoij(iatom)%nspden,rhoij_input=rhoijtmp)
612 0 : ABI_FREE(rhoijtmp)
613 : end do
614 : end if
615 : end if
616 :
617 28280 : ABI_FREE(rhoijrespc)
618 :
619 : !PAW: restore rhoij from compact storage
620 28280 : if (usepaw==1.and.dtset%iscf/=5.and.dtset%iscf/=6) then
621 782 : indx=-dplex
622 1921 : do iatom=1,my_natom
623 4556 : ABI_MALLOC(rhoijtmp,(cplex*qphase*pawrhoij(iatom)%lmn2_size,pawrhoij(iatom)%nspden))
624 243057 : rhoijtmp=zero
625 2278 : do iq=1,qphase
626 1139 : iq0=merge(0,cplex*pawrhoij(iatom)%lmn2_size,iq==1)
627 1139 : if (pawrhoij(iatom)%lmnmix_sz<pawrhoij(iatom)%lmn2_size) then
628 0 : do ispden=1,pawrhoij(iatom)%nspden
629 0 : jrhoij=iq0+1
630 0 : do irhoij=1,pawrhoij(iatom)%nrhoijsel
631 0 : klmn=iq0+cplex*pawrhoij(iatom)%rhoijselect(irhoij)-dplex
632 0 : rhoijtmp(klmn:klmn+dplex,ispden)=pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+dplex,ispden)
633 0 : jrhoij=jrhoij+cplex
634 : end do
635 : end do
636 : end if
637 4032 : do ispden=1,pawrhoij(iatom)%nspden
638 243057 : do kmix=1,pawrhoij(iatom)%lmnmix_sz
639 240164 : indx=indx+cplex;klmn=iq0+cplex*pawrhoij(iatom)%kpawmix(kmix)-dplex
640 482082 : rhoijtmp(klmn:klmn+dplex,ispden)=vpaw(indx:indx+dplex)
641 : end do
642 : end do
643 : end do
644 : call pawrhoij_filter(pawrhoij(iatom)%rhoijp,pawrhoij(iatom)%rhoijselect,&
645 : & pawrhoij(iatom)%nrhoijsel,pawrhoij(iatom)%cplex_rhoij,pawrhoij(iatom)%qphase,&
646 1139 : & pawrhoij(iatom)%lmn2_size,pawrhoij(iatom)%nspden,rhoij_input=rhoijtmp)
647 1921 : ABI_FREE(rhoijtmp)
648 : end do
649 : end if
650 28280 : ABI_FREE(vpaw)
651 :
652 : !Eventually write the data on disk and deallocate f_fftgr_disk
653 28280 : call abi_mixing_eval_deallocate(mix)
654 28280 : if (with_vtau) call abi_mixing_eval_deallocate(mix_mgga)
655 :
656 28280 : call timab(904,2,tsec)
657 :
658 : !Restore potential
659 28280 : if (ispmix==1.and.nfft==nfftmix) then
660 558908984 : vtrial=vtrial0
661 7654241 : if (with_vtau) vtau(:,:,1)=vtau0(:,:)
662 566 : else if (nfft==nfftmix) then
663 0 : do ispden=1,dtset%nspden
664 0 : call fourdp(1,vtrial0(:,ispden),vtrial(:,ispden),+1,mpi_enreg,nfft,1,ngfft,tim_fourdp)
665 : end do
666 0 : if (with_vtau) then
667 0 : do ispden=1,dtset%nspden
668 0 : call fourdp(1,vtau0(:,ispden),vtau(:,ispden,1),+1,mpi_enreg,nfft,1,ngfft,tim_fourdp)
669 : end do
670 : end if
671 : else
672 1484 : do ispden=1,dtset%nspden
673 3579778 : do ifft=1,nfftmix
674 3578860 : jfft=mixtofft(ifft)
675 3578860 : vtrialg(1,jfft,ispden)=vtrial0(2*ifft-1,ispden)
676 3579778 : vtrialg(2,jfft,ispden)=vtrial0(2*ifft ,ispden)
677 : end do
678 1484 : call fourdp(1,vtrialg(:,:,ispden),vtrial(:,ispden),+1,mpi_enreg,nfft,1,ngfft,tim_fourdp)
679 : end do
680 566 : ABI_FREE(vtrialg)
681 566 : if (with_vtau) then
682 50 : do ispden=1,dtset%nspden
683 102425 : do ifft=1,nfftmix
684 102400 : jfft=mixtofft(ifft)
685 102400 : vtaug(1,jfft,ispden)=vtau0(2*ifft-1,ispden)
686 102425 : vtaug(2,jfft,ispden)=vtau0(2*ifft ,ispden)
687 : end do
688 50 : call fourdp(1,vtaug(:,:,ispden),vtau(:,ispden,1),+1,mpi_enreg,nfft,1,ngfft,tim_fourdp)
689 : end do
690 25 : ABI_FREE(vtaug)
691 : end if
692 : end if
693 28280 : ABI_FREE(vtrial0)
694 28280 : ABI_FREE(vtau0)
695 :
696 : !In case of metaGGA, re-compute vtau gradient
697 28280 : if (with_vtau.and.mix_mgga%iscf/=ABI_MIXING_NONE) then
698 0 : call xcpot(1,gprimd,0,0,mpi_enreg,nfft,ngfft,2,dtset%nspden,0,[zero,zero,zero],vxctau=vtau)
699 : end if
700 :
701 28280 : call timab(905,1,tsec)
702 :
703 : !------Treat the mean of the potential
704 :
705 : !Compute the mean of the new vtrial
706 28280 : call mean_fftr(vtrial,vmean,nfft,nfftot,dtset%nspden,mpi_comm_sphgrid)
707 :
708 : !Reset the mean of the new vtrial, to the value vnew_mean
709 : !When spin-polarized and fixed occupation numbers,
710 : !treat separately spin up and spin down.
711 : !Otherwise, use only global mean
712 28280 : if (mix%iscf/=ABI_MIXING_NONE) then
713 63576 : do ispden=1,dtset%nspden
714 35296 : if (dtset%nspden==2.and.dtset%occopt>=3.and. &
715 : & abs(dtset%spinmagntarget+99.99_dp)<1.0d-10)then
716 1340 : vme=(vnew_mean(1)+vnew_mean(2)-vmean(1)-vmean(2))*half
717 : else
718 33956 : vme=vnew_mean(ispden)-vmean(ispden)
719 : end if
720 577586908 : vtrial(:,ispden)=vtrial(:,ispden)+vme
721 : end do
722 : else
723 : ! If no mixing, just re-add the residual mean
724 0 : do ispden=1,dtset%nspden
725 0 : vtrial(:,ispden)=vtrial(:,ispden)+vres_mean(ispden)
726 : end do
727 : end if
728 :
729 28280 : if(moved_atm_inside==1 .and. istep/=nstep )then
730 106 : if(abs(dtset%densfor_pred)==1.or.abs(dtset%densfor_pred)==4)then
731 : ! Subtract current local psp, but also vxc (for core charges)
732 86 : do ispden=1,dtset%nspden
733 461910 : vtrial(:,ispden)=vtrial(:,ispden)-vpsp(:)*identity(ispden)-vxc(:,ispden)
734 : end do
735 : else if(abs(dtset%densfor_pred)==2.or.abs(dtset%densfor_pred)==5.or.abs(dtset%densfor_pred)==6)then
736 : ! Subtract current vpsp+Hxc from vtrial. This should be rationalized later
737 84 : do ispden=1,dtset%nspden
738 1573584 : vtrial(:,ispden)=vtrial(:,ispden)-(vpsp(:)+vhartr(:))*identity(ispden)-vxc(:,ispden)
739 : end do
740 : end if
741 : end if
742 :
743 : !In WVL: copy vtrial to BigDFT object:
744 28280 : if(dtset%usewvl==1) then
745 0 : call wvl_vtrial_abi2big(1,vtrial,wvl%den)
746 : end if
747 :
748 28280 : call timab(905,2,tsec)
749 28280 : call timab(93,2,tsec)
750 :
751 56560 : end subroutine newvtr
752 : !!***
753 :
754 : end module m_newvtr
755 : !!***
|
