Branch data 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 : 56742 : subroutine newvtr(atindx,dbl_nnsclo,dielar,dielinv,dielstrt,&
187 : 28371 : & dtn_pc,dtset,etotal,fcart,ffttomix,&
188 : 28371 : & gmet,grhf,gsqcut,&
189 : : & initialized,ispmix,&
190 : : & istep,&
191 : 28371 : & kg_diel,kxc,mgfft,mix,mixtofft,&
192 : 28371 : & moved_atm_inside,mpi_enreg,my_natom,nattyp,nfft,nfftmix,&
193 : : & ngfft,ngfftmix,nkxc,npawmix,npwdiel,&
194 : : & nstep,ntypat,n1xccc,&
195 : 28371 : & pawrhoij,&
196 : 28371 : & ph1d,&
197 : : & precon,&
198 : 28371 : & psps,rhor,rprimd,susmat,usepaw,&
199 : 28371 : & vhartr,vnew_mean,vpsp,vresid,vres_mean,&
200 : 28371 : & vtrial,vxc,xred,&
201 : : & nfftf,&
202 : 28371 : & pawtab,&
203 : 28371 : & rhog,&
204 : : & wvl,&
205 [ + - + - : 28371 : & 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 : 56742 : real(dp) :: gprimd(3,3),rmet(3,3),tsec(2),vmean(dtset%nspden)
262 : 28371 : real(dp),allocatable :: rhoijrespc(:)
263 : 28371 : real(dp),allocatable :: rhoijtmp(:,:)
264 : 28371 : real(dp),allocatable :: vresid0(:,:),vrespc(:,:),vreswk(:,:),vtrialg(:,:,:)
265 : 28371 : real(dp),allocatable :: vtauresid0(:,:),vtaurespc(:,:),vtaug(:,:,:),vtau0(:,:)
266 : 28371 : 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 : 28371 : call timab(93,1,tsec)
276 : 28371 : call timab(901,1,tsec)
277 : 28371 : tim_fourdp=8
278 : :
279 : : !mpicomm over spherical grid:
280 : 28371 : mpi_comm_sphgrid=mpi_enreg%comm_fft
281 [ - + ]: 28371 : if(dtset%usewvl==1) mpi_comm_sphgrid=mpi_enreg%comm_wvl
282 : :
283 : : !Compatibility tests
284 [ - + ]: 28371 : if(nfftmix>nfft) then
285 : 0 : ABI_BUG(' nfftmix>nfft not allowed !')
286 : : end if
287 : :
288 [ + + - + ]: 28371 : 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 [ - + ]: 28371 : 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 : 28371 : with_vtau = (present(vtau).and.present(vtauresid))
305 [ + + + - : 113484 : 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 : 28232 : usevtau=0;if (with_vtau) usevtau=1
321 : :
322 [ + + + + : 28371 : 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 : 28371 : dielng=dielar(2)
328 : 28371 : diemix=dielar(4)
329 : 28371 : n1=ngfft(1)
330 : 28371 : n2=ngfft(2)
331 : 28371 : n3=ngfft(3)
332 [ + + + - ]: 28371 : 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 [ + + ]: 113484 : nfftot=PRODUCT(ngfft(1:3))
341 : :
342 : : !Compute different geometric tensor, as well as ucvol, from rprimd
343 : 28371 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
344 : :
345 [ + - ]: 28371 : if(dtset%usewvl==0) then
346 : 28371 : 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 : 28371 : 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 [ + + ]: 63938 : do ispden=1,dtset%nspden
372 [ + + + + : 35567 : 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 : 34227 : vme=vmean(ispden)
377 : : end if
378 [ + + ]: 582522190 : vtrial(:,ispden)=vtrial(:,ispden)-vme
379 : : end do
380 : :
381 : 28371 : call timab(901,2,tsec)
382 : :
383 : 28371 : call timab(902,1,tsec)
384 : :
385 : : !Select components of potential to be mixed
386 [ - + + - : 113484 : ABI_MALLOC(vtrial0,(ispmix*nfftmix,dtset%nspden))
- + + - +
- + - - +
- + ]
387 [ - + + - : 85113 : ABI_MALLOC(vresid0,(ispmix*nfftmix,dtset%nspden))
- + + - -
+ - + ]
388 [ - + + - : 113484 : ABI_MALLOC(vtau0,(ispmix*nfftmix,dtset%nspden*usevtau))
+ + + - +
- + + - +
- + ]
389 [ - + + - : 85113 : ABI_MALLOC(vtauresid0,(ispmix*nfftmix,dtset%nspden*usevtau))
+ + + - -
+ - + ]
390 [ + + + - ]: 28371 : if (ispmix==1.and.nfft==nfftmix) then
391 [ + + + + : 1127688532 : vtrial0=vtrial;vresid0=vresid
+ - - + -
- + + +
+ ]
392 [ + + ]: 27805 : 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 [ - + + - : 85113 : ABI_MALLOC(vrespc,(ispmix*nfftmix,dtset%nspden))
- + + - -
+ - + ]
454 [ - + + - : 85113 : ABI_MALLOC(vtaurespc,(ispmix*nfftmix,dtset%nspden*usevtau))
+ + + - -
+ - + ]
455 [ + + + - : 85113 : ABI_MALLOC(vpaw,(npawmix*usepaw))
+ - + + -
+ - + ]
456 [ + + ]: 28371 : if (usepaw==1) then
457 [ - + + - : 2346 : ABI_MALLOC(rhoijrespc,(npawmix))
+ - + - -
+ - + ]
458 : : else
459 [ - + ]: 27589 : ABI_MALLOC(rhoijrespc,(0))
460 : : end if
461 : :
462 : 28371 : call timab(902,2,tsec)
463 : 28371 : call timab(903,1,tsec)
464 : :
465 [ + - ]: 28371 : 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 : 28371 : & 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 [ + + ]: 28371 : 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 : 28371 : call timab(903,2,tsec)
486 : 28371 : call timab(904,1,tsec)
487 : :
488 : : !------Compute new vtrial and eventual new atomic positions
489 : :
490 [ + - ]: 28371 : if (mix%n_fftgr>0) then
491 : 28371 : i_vresid1=mix%i_vresid(1)
492 : 28371 : i_vrespc1=mix%i_vrespc(1)
493 : : end if
494 : :
495 : : !Initialise working arrays for the mixing object.
496 [ + + ]: 28371 : if (moved_atm_inside == 1) then
497 : 110 : call abi_mixing_use_moving_atoms(mix, dtset%natom, xred, dtn_pc)
498 : : end if
499 : 28371 : call abi_mixing_eval_allocate(mix, istep)
500 : : !Copy current step arrays.
501 [ + + ]: 28371 : 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 : 28261 : & arr_respc = vrespc)
507 : : end if
508 [ - + ]: 28371 : if (errid /= AB7_NO_ERROR) then
509 : 0 : ABI_ERROR(message)
510 : : end if
511 [ + + ]: 28371 : 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 [ - + ]: 28371 : ABI_FREE(vresid0)
520 [ - + ]: 28371 : ABI_FREE(vrespc)
521 [ - + ]: 28371 : ABI_FREE(vtauresid0)
522 [ - + ]: 28371 : ABI_FREE(vtaurespc)
523 : :
524 : : !PAW: either use the array f_paw or the array f_paw_disk
525 [ + + ]: 28371 : 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 [ + + ]: 28371 : if(mpi_enreg%paral_kgb==1)then
554 : 2372 : mpicomm=mpi_enreg%comm_fft
555 : 2372 : mpi_summarize=.true.
556 : : else
557 : 25999 : mpicomm=0
558 : 25999 : mpi_summarize=.false.
559 : : end if
560 : :
561 : 28371 : reset = .false.
562 [ + + ]: 28371 : 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 [ - + - + ]: 28371 : & comm_atom=mpi_enreg%comm_atom)
568 [ - + ]: 28371 : if (errid == AB7_ERROR_MIXING_INC_NNSLOOP) then
569 : 0 : dbl_nnsclo = 1
570 [ - + ]: 28371 : else if (errid /= AB7_NO_ERROR) then
571 : 0 : ABI_ERROR(message)
572 : : end if
573 [ + + ]: 28371 : 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 [ + + ]: 28371 : 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 [ - + ]: 28371 : ABI_FREE(rhoijrespc)
618 : :
619 : : !PAW: restore rhoij from compact storage
620 [ + + + - : 28371 : 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 [ - + ]: 28371 : ABI_FREE(vpaw)
651 : :
652 : : !Eventually write the data on disk and deallocate f_fftgr_disk
653 : 28371 : call abi_mixing_eval_deallocate(mix)
654 [ + + ]: 28371 : if (with_vtau) call abi_mixing_eval_deallocate(mix_mgga)
655 : :
656 : 28371 : call timab(904,2,tsec)
657 : :
658 : : !Restore potential
659 [ + + + - ]: 28371 : if (ispmix==1.and.nfft==nfftmix) then
660 [ + + + + ]: 563844266 : vtrial=vtrial0
661 [ + + + + : 7654332 : 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 [ - + ]: 28371 : ABI_FREE(vtrial0)
694 [ - + ]: 28371 : ABI_FREE(vtau0)
695 : :
696 : : !In case of metaGGA, re-compute vtau gradient
697 [ + + - + ]: 28371 : 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 : 28371 : call timab(905,1,tsec)
702 : :
703 : : !------Treat the mean of the potential
704 : :
705 : : !Compute the mean of the new vtrial
706 : 28371 : 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 [ - + ]: 28371 : if (mix%iscf/=ABI_MIXING_NONE) then
713 [ + + ]: 63938 : do ispden=1,dtset%nspden
714 [ + + + + : 35567 : 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 : 34227 : vme=vnew_mean(ispden)-vmean(ispden)
719 : : end if
720 [ + + ]: 582522190 : 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 [ + + + + ]: 28371 : 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 [ - + ]: 28371 : if(dtset%usewvl==1) then
745 : 0 : call wvl_vtrial_abi2big(1,vtrial,wvl%den)
746 : : end if
747 : :
748 : 28371 : call timab(905,2,tsec)
749 : 28371 : call timab(93,2,tsec)
750 : :
751 [ - + - + ]: 56742 : end subroutine newvtr
752 : : !!***
753 : :
754 : : end module m_newvtr
755 : : !!***
|
