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