TABLE OF CONTENTS
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.
COPYRIGHT
Copyright (C) 1998-2018 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 . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
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 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 vpsp(nfft)=array for holding local psp vresid(nfft,nspden)=array for the residual of the potential vxc(nfft,nspden)=exchange-correlation potential (hartree) 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 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).
PARENTS
scfcv
CHILDREN
ab7_mixing_copy_current_step,ab7_mixing_eval,ab7_mixing_eval_allocate ab7_mixing_eval_deallocate,ab7_mixing_use_moving_atoms,fourdp,mean_fftr metric,prcref_pma,timab,wvl_prcref,wvl_vtrial_abi2big
SOURCE
140 #if defined HAVE_CONFIG_H 141 #include "config.h" 142 #endif 143 144 #include "abi_common.h" 145 146 subroutine newvtr(atindx,dbl_nnsclo,dielar,dielinv,dielstrt,& 147 & dtn_pc,dtset,etotal,fcart,ffttomix,& 148 & gmet,grhf,gsqcut,& 149 & initialized,ispmix,& 150 & istep,& 151 & kg_diel,kxc,mgfft,mix,mixtofft,& 152 & moved_atm_inside,mpi_enreg,my_natom,nattyp,nfft,nfftmix,& 153 & ngfft,ngfftmix,nkxc,npawmix,npwdiel,& 154 & nstep,ntypat,n1xccc,& 155 & pawrhoij,& 156 & ph1d,& 157 & psps,rhor,rprimd,susmat,usepaw,& 158 & vhartr,vnew_mean,vpsp,vresid,& 159 & vtrial,vxc,xred,& 160 & nfftf,& 161 & pawtab,& 162 & rhog,& 163 & wvl) 164 165 use defs_basis 166 use defs_datatypes 167 use defs_abitypes 168 use defs_wvltypes 169 use m_profiling_abi 170 use m_errors 171 use m_abi2big 172 use m_ab7_mixing 173 use m_cgtools 174 175 use m_geometry, only : metric 176 use m_pawtab, only : pawtab_type 177 use m_pawrhoij,only : pawrhoij_type 178 179 !This section has been created automatically by the script Abilint (TD). 180 !Do not modify the following lines by hand. 181 #undef ABI_FUNC 182 #define ABI_FUNC 'newvtr' 183 use interfaces_18_timing 184 use interfaces_53_ffts 185 use interfaces_68_rsprc, except_this_one => newvtr 186 !End of the abilint section 187 188 implicit none 189 190 !Arguments------------------------------- 191 ! WARNING 192 ! BEWARE THERE IS TWO DIFFERENT SIZE DECLARED FOR ARRAY NHAT IN RHOTOV AND RHOHXC 193 ! THIS MIGHT RESULT IN A BUG 194 !scalars 195 integer,intent(in) :: dielstrt,initialized,ispmix,istep,mgfft 196 integer,intent(in) :: moved_atm_inside,my_natom,n1xccc,nfft 197 integer,intent(in) :: nfftf,nfftmix,nkxc,npawmix,npwdiel,nstep 198 integer,intent(in) :: ntypat,usepaw 199 integer,intent(inout) :: dbl_nnsclo 200 real(dp),intent(in) :: etotal,gsqcut 201 type(MPI_type),intent(in) :: mpi_enreg 202 type(dataset_type),intent(in) :: dtset 203 type(ab7_mixing_object),intent(inout) :: mix 204 type(pseudopotential_type),intent(in) :: psps 205 type(wvl_data), intent(inout) :: wvl 206 !arrays 207 integer,intent(in) :: atindx(dtset%natom) 208 integer,intent(in) :: ffttomix(nfft*(1-nfftmix/nfft)) 209 integer,intent(in) :: kg_diel(3,npwdiel) 210 integer,intent(in) :: mixtofft(nfftmix*(1-nfftmix/nfft)),nattyp(ntypat) 211 integer,intent(in) :: ngfft(18),ngfftmix(18) 212 real(dp),intent(in) :: dielar(7) 213 real(dp),intent(in) :: fcart(3,dtset%natom),grhf(3,dtset%natom) 214 real(dp),intent(in) :: rprimd(3,3) 215 real(dp),intent(in) :: susmat(2,npwdiel,dtset%nspden,npwdiel,dtset%nspden) 216 real(dp),intent(in) :: vhartr(nfft),vnew_mean(dtset%nspden) 217 real(dp),intent(in) :: vxc(nfft,dtset%nspden) 218 real(dp),intent(inout) :: dielinv(2,npwdiel,dtset%nspden,npwdiel,dtset%nspden) 219 real(dp),intent(inout), target :: dtn_pc(3,dtset%natom) 220 real(dp),intent(inout) :: gmet(3,3) 221 real(dp),intent(inout) :: kxc(nfft,nkxc),ph1d(2,3*(2*mgfft+1)*dtset%natom) 222 real(dp),intent(inout) :: rhog(2,nfftf),vpsp(nfft) 223 real(dp),intent(inout), target :: rhor(nfft,dtset%nspden) 224 real(dp),intent(inout) :: vresid(nfft,dtset%nspden),vtrial(nfft,dtset%nspden) 225 real(dp),intent(inout), target :: xred(3,dtset%natom) 226 type(pawrhoij_type),intent(inout) :: pawrhoij(my_natom*usepaw) 227 type(pawtab_type),intent(in) :: pawtab(ntypat*usepaw) 228 229 !Local variables------------------------------- 230 !scalars 231 integer :: cplex,dplex,i_vresid1,i_vrespc1 232 ! integer :: i1,i2,i3,ifft2,ifft3,ifft4,ifft5,ii1,ii2,ii3,ii4,ii5 233 integer :: iatom,ifft,indx,irhoij,mpicomm,mpi_comm_sphgrid 234 integer :: ispden,jfft,jrhoij,klmn,kmix,n1,n2,n3,nfftot,nselect 235 integer :: errid,tim_fourdp 236 logical :: mpi_summarize,reset 237 real(dp) :: dielng,diemix,fact,ucvol,ucvol_local,vme 238 character(len=500) :: message 239 !arrays 240 real(dp),parameter :: identity(4)=(/one,one,zero,zero/) 241 real(dp) :: gprimd(3,3),rmet(3,3),tsec(2),vmean(dtset%nspden) 242 real(dp),allocatable :: rhoijrespc(:) 243 real(dp),allocatable :: rhoijtmp(:,:) 244 real(dp),allocatable :: vresid0(:,:),vrespc(:,:),vreswk(:,:),vtrialg(:,:,:) 245 real(dp),pointer :: vtrial0(:,:),vpaw(:) 246 247 ! ************************************************************************* 248 249 !DEBUG 250 !write(std_out,*)' newvtr : enter ' 251 !write(std_out,*)' newvtr : ispmix,nfft,nfftmix=',ispmix,nfft,nfftmix 252 !ENDDEBUG 253 254 call timab(93,1,tsec) 255 call timab(901,1,tsec) 256 tim_fourdp=8 257 258 !mpicomm over spherical grid: 259 mpi_comm_sphgrid=mpi_enreg%comm_fft 260 if(dtset%usewvl==1) mpi_comm_sphgrid=mpi_enreg%comm_wvl 261 262 !Compatibility tests 263 if(nfftmix>nfft) then 264 MSG_BUG(' nfftmix>nfft not allowed !') 265 end if 266 267 if(ispmix/=2.and.nfftmix/=nfft) then 268 message = ' nfftmix/=nfft allowed only when ispmix=2 !' 269 MSG_BUG(message) 270 end if 271 272 if(dtset%usewvl==1) then 273 if(dtset%wvl_bigdft_comp==1) then 274 message = 'newvtr: usewvl == 1 and wvl_bigdft_comp==1 not allowed (use wvl_newtr() instead)!' 275 MSG_BUG(message) 276 end if 277 if(ispmix/=1 .or. nfftmix/=nfft) then 278 MSG_BUG('newvtr: nfftmix/=nfft, ispmix/=1 not allowed for wavelets') 279 end if 280 end if 281 282 if(usepaw==1.and.dtset%nspden==4.and.dtset%pawoptmix==1) then 283 message = ' pawoptmix=1 is not compatible with nspden=4 !' 284 MSG_ERROR(message) 285 end if 286 287 dielng=dielar(2) 288 diemix=dielar(4) 289 n1=ngfft(1) 290 n2=ngfft(2) 291 n3=ngfft(3) 292 if (usepaw==1.and.my_natom>0) then 293 cplex=pawrhoij(1)%cplex;dplex=cplex-1 294 else 295 cplex=0;dplex=0 296 end if 297 298 !Get size of FFT grid 299 nfftot=PRODUCT(ngfft(1:3)) 300 301 !Compute different geometric tensor, as well as ucvol, from rprimd 302 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 303 304 if(dtset%usewvl==0) then 305 ucvol_local=ucvol 306 #if defined HAVE_BIGDFT 307 else 308 ucvol_local = product(wvl%den%denspot%dpbox%hgrids) * real(nfftot, dp) 309 #endif 310 end if 311 312 !------Treat the mean of potentiel residual 313 314 !Special care must be taken with components of the 315 !potential that are associated with NO density change. 316 !In general, only the global mean of the potential has 317 !such an anomalous feature. However, in the spin 318 !polarized cas with fixed occupancies, also the 319 !mean of each spin-potential (independently of the other) 320 !has such a behaviour. The trick is to remove these 321 !variables before going in the predictive routines, 322 !then to put them back 323 324 !Compute the mean of the old vtrial 325 call mean_fftr(vtrial,vmean,nfft,nfftot,dtset%nspden,mpi_comm_sphgrid) 326 327 !When (collinear) spin-polarized and fixed occupation numbers, 328 !treat separately spin up and spin down. 329 !Otherwise, use only global mean 330 do ispden=1,dtset%nspden 331 if (dtset%nspden==2.and.dtset%occopt>=3.and. & 332 & abs(dtset%spinmagntarget+99.99_dp)<1.0d-10)then 333 vme=(vmean(1)+vmean(2))*half 334 else 335 vme=vmean(ispden) 336 end if 337 vtrial(:,ispden)=vtrial(:,ispden)-vme 338 end do 339 340 call timab(901,2,tsec) 341 342 !Select components of potential to be mixed 343 ABI_ALLOCATE(vtrial0,(ispmix*nfftmix,dtset%nspden)) 344 ABI_ALLOCATE(vresid0,(ispmix*nfftmix,dtset%nspden)) 345 if (ispmix==1.and.nfft==nfftmix) then 346 vtrial0=vtrial;vresid0=vresid 347 else if (nfft==nfftmix) then 348 do ispden=1,dtset%nspden 349 call fourdp(1,vtrial0(:,ispden),vtrial(:,ispden),-1,mpi_enreg,nfft,ngfft,dtset%paral_kgb,tim_fourdp) 350 call fourdp(1,vresid0(:,ispden),vresid(:,ispden),-1,mpi_enreg,nfft,ngfft,dtset%paral_kgb,tim_fourdp) 351 end do 352 else 353 ABI_ALLOCATE(vtrialg,(2,nfft,dtset%nspden)) 354 ABI_ALLOCATE(vreswk,(2,nfft)) 355 do ispden=1,dtset%nspden 356 fact=dielar(4);if (ispden>1) fact=dielar(7) 357 call fourdp(1,vtrialg(:,:,ispden),vtrial(:,ispden),-1,mpi_enreg,nfft,ngfft,dtset%paral_kgb,tim_fourdp) 358 call fourdp(1,vreswk,vresid(:,ispden),-1,mpi_enreg,nfft,ngfft,dtset%paral_kgb,tim_fourdp) 359 do ifft=1,nfft 360 if (ffttomix(ifft)>0) then 361 jfft=2*ffttomix(ifft) 362 vtrial0(jfft-1,ispden)=vtrialg(1,ifft,ispden) 363 vtrial0(jfft ,ispden)=vtrialg(2,ifft,ispden) 364 vresid0(jfft-1,ispden)=vreswk(1,ifft) 365 vresid0(jfft ,ispden)=vreswk(2,ifft) 366 else 367 vtrialg(:,ifft,ispden)=vtrialg(:,ifft,ispden)+fact*vreswk(:,ifft) 368 end if 369 end do 370 end do 371 ABI_DEALLOCATE(vreswk) 372 end if 373 374 call timab(902,1,tsec) 375 376 !Choice of preconditioner governed by iprcel, densfor_pred and iprcfc 377 ABI_ALLOCATE(vrespc,(ispmix*nfftmix,dtset%nspden)) 378 ABI_ALLOCATE(vpaw,(npawmix*usepaw)) 379 if (usepaw==1) then 380 ABI_ALLOCATE(rhoijrespc,(npawmix)) 381 else 382 ABI_ALLOCATE(rhoijrespc,(0)) 383 end if 384 385 call timab(902,2,tsec) 386 call timab(903,1,tsec) 387 388 if(dtset%usewvl==0) then 389 call prcref_PMA(atindx,dielar,dielinv,dielstrt,dtn_pc,dtset,fcart,ffttomix,gmet,gsqcut,& 390 & istep,kg_diel,kxc,mgfft,moved_atm_inside,mpi_enreg,my_natom,& 391 & nattyp,nfft,nfftmix,ngfft,ngfftmix,nkxc,npawmix,npwdiel,ntypat,n1xccc,& 392 & ispmix,0,pawrhoij,ph1d,psps,rhog,rhoijrespc,rhor,rprimd,susmat,& 393 & vhartr,vpsp,vresid0,vrespc,vxc,xred,& 394 & etotal,pawtab,wvl) 395 else 396 call wvl_prcref(dielar,dtset%iprcel,my_natom,nfftmix,npawmix,dtset%nspden,pawrhoij,& 397 & rhoijrespc,psps%usepaw,vresid0,vrespc) 398 end if 399 400 call timab(903,2,tsec) 401 call timab(904,1,tsec) 402 403 !------Compute new vtrial and eventual new atomic positions 404 405 i_vresid1=mix%i_vresid(1) 406 i_vrespc1=mix%i_vrespc(1) 407 408 !Initialise working arrays for the mixing object. 409 if (moved_atm_inside == 1) then 410 call ab7_mixing_use_moving_atoms(mix, dtset%natom, xred, dtn_pc) 411 end if 412 call ab7_mixing_eval_allocate(mix, istep) 413 !Copy current step arrays. 414 if (moved_atm_inside == 1) then 415 call ab7_mixing_copy_current_step(mix, vresid0, errid, message, & 416 & arr_respc = vrespc, arr_atm = grhf) 417 else 418 call ab7_mixing_copy_current_step(mix, vresid0, errid, message, & 419 & arr_respc = vrespc) 420 end if 421 if (errid /= AB7_NO_ERROR) then 422 MSG_ERROR(message) 423 end if 424 ABI_DEALLOCATE(vresid0) 425 426 !PAW: either use the array f_paw or the array f_paw_disk 427 if (usepaw==1) then 428 indx=-dplex 429 do iatom=1,my_natom 430 do ispden=1,pawrhoij(iatom)%nspden 431 ABI_ALLOCATE(rhoijtmp,(cplex*pawrhoij(iatom)%lmn2_size,1)) 432 rhoijtmp=zero 433 jrhoij=1 434 do irhoij=1,pawrhoij(iatom)%nrhoijsel 435 klmn=cplex*pawrhoij(iatom)%rhoijselect(irhoij)-dplex 436 rhoijtmp(klmn:klmn+dplex,1)=pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+dplex,ispden) 437 jrhoij=jrhoij+cplex 438 end do 439 do kmix=1,pawrhoij(iatom)%lmnmix_sz 440 indx=indx+cplex;klmn=cplex*pawrhoij(iatom)%kpawmix(kmix)-dplex 441 vpaw(indx:indx+dplex)=rhoijtmp(klmn:klmn+dplex,1)-pawrhoij(iatom)%rhoijres(klmn:klmn+dplex,ispden) 442 mix%f_paw(indx:indx+dplex,i_vresid1)=pawrhoij(iatom)%rhoijres(klmn:klmn+dplex,ispden) 443 mix%f_paw(indx:indx+dplex,i_vrespc1)=rhoijrespc(indx:indx+dplex) 444 end do 445 ABI_DEALLOCATE(rhoijtmp) 446 end do 447 end do 448 end if 449 450 !------Prediction of the components of the potential associated with a density change 451 452 !Init mpicomm 453 if(mpi_enreg%paral_kgb==1)then 454 mpicomm=mpi_enreg%comm_fft 455 mpi_summarize=.true. 456 else 457 mpicomm=0 458 mpi_summarize=.false. 459 end if 460 461 reset = .false. 462 if (initialized == 0) reset = .true. 463 call ab7_mixing_eval(mix, vtrial0, istep, nfftot, ucvol_local, & 464 & mpicomm, mpi_summarize, errid, message, & 465 & reset = reset, isecur = dtset%isecur, & 466 & pawopt = dtset%pawoptmix, pawarr = vpaw, etotal = etotal, potden = rhor, & 467 & comm_atom=mpi_enreg%comm_atom) 468 469 if (errid == AB7_ERROR_MIXING_INC_NNSLOOP) then 470 dbl_nnsclo = 1 471 else if (errid /= AB7_NO_ERROR) then 472 MSG_ERROR(message) 473 end if 474 475 !PAW: apply a simple mixing to rhoij (this is temporary) 476 if(dtset%iscf==5 .or. dtset%iscf==6)then 477 if (usepaw==1) then 478 indx=-dplex 479 do iatom=1,my_natom 480 ABI_ALLOCATE(rhoijtmp,(cplex*pawrhoij(iatom)%lmn2_size,pawrhoij(iatom)%nspden)) 481 rhoijtmp=zero 482 if (pawrhoij(iatom)%lmnmix_sz<pawrhoij(iatom)%lmn2_size) then 483 do ispden=1,pawrhoij(iatom)%nspden 484 do kmix=1,pawrhoij(iatom)%lmnmix_sz 485 indx=indx+cplex;klmn=cplex*pawrhoij(iatom)%kpawmix(kmix)-dplex 486 rhoijtmp(klmn:klmn+dplex,ispden)=rhoijrespc(indx:indx+dplex) & 487 & -pawrhoij(iatom)%rhoijres(klmn:klmn+dplex,ispden) 488 end do 489 end do 490 end if 491 do ispden=1,pawrhoij(iatom)%nspden 492 jrhoij=1 493 do irhoij=1,pawrhoij(iatom)%nrhoijsel 494 klmn=cplex*pawrhoij(iatom)%rhoijselect(irhoij)-dplex 495 rhoijtmp(klmn:klmn+dplex,ispden)=rhoijtmp(klmn:klmn+dplex,ispden) & 496 & +pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+dplex,ispden) 497 jrhoij=jrhoij+cplex 498 end do 499 end do 500 nselect=0 501 do klmn=1,pawrhoij(iatom)%lmn2_size 502 if (any(abs(rhoijtmp(cplex*klmn-dplex:cplex*klmn,:))>tol10)) then 503 nselect=nselect+1 504 pawrhoij(iatom)%rhoijselect(nselect)=klmn 505 do ispden=1,pawrhoij(iatom)%nspden 506 pawrhoij(iatom)%rhoijp(cplex*nselect-dplex:cplex*nselect,ispden)=& 507 & rhoijtmp(cplex*klmn-dplex:cplex*klmn,ispden) 508 end do 509 end if 510 end do 511 pawrhoij(iatom)%nrhoijsel=nselect 512 ABI_DEALLOCATE(rhoijtmp) 513 end do 514 end if 515 end if 516 517 ABI_DEALLOCATE(rhoijrespc) 518 519 !PAW: restore rhoij from compact storage 520 if (usepaw==1.and.dtset%iscf/=5.and.dtset%iscf/=6) then 521 indx=-dplex 522 do iatom=1,my_natom 523 ABI_ALLOCATE(rhoijtmp,(cplex*pawrhoij(iatom)%lmn2_size,pawrhoij(iatom)%nspden)) 524 rhoijtmp=zero 525 if (pawrhoij(iatom)%lmnmix_sz<pawrhoij(iatom)%lmn2_size) then 526 do ispden=1,pawrhoij(iatom)%nspden 527 jrhoij=1 528 do irhoij=1,pawrhoij(iatom)%nrhoijsel 529 klmn=cplex*pawrhoij(iatom)%rhoijselect(irhoij)-dplex 530 rhoijtmp(klmn:klmn+dplex,ispden)=pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+dplex,ispden) 531 jrhoij=jrhoij+cplex 532 end do 533 end do 534 end if 535 do ispden=1,pawrhoij(iatom)%nspden 536 do kmix=1,pawrhoij(iatom)%lmnmix_sz 537 indx=indx+cplex;klmn=cplex*pawrhoij(iatom)%kpawmix(kmix)-dplex 538 rhoijtmp(klmn:klmn+dplex,ispden)=vpaw(indx:indx+dplex) 539 end do 540 end do 541 nselect=0 542 if (cplex==1) then 543 do klmn=1,pawrhoij(iatom)%lmn2_size 544 if (any(abs(rhoijtmp(klmn,:))>tol10)) then 545 nselect=nselect+1 546 pawrhoij(iatom)%rhoijselect(nselect)=klmn 547 do ispden=1,pawrhoij(iatom)%nspden 548 pawrhoij(iatom)%rhoijp(nselect,ispden)=rhoijtmp(klmn,ispden) 549 end do 550 end if 551 end do 552 else 553 do klmn=1,pawrhoij(iatom)%lmn2_size 554 if (any(abs(rhoijtmp(2*klmn-1:2*klmn,:))>tol10)) then 555 nselect=nselect+1 556 pawrhoij(iatom)%rhoijselect(nselect)=klmn 557 do ispden=1,pawrhoij(iatom)%nspden 558 pawrhoij(iatom)%rhoijp(2*nselect-1:2*nselect,ispden)=rhoijtmp(2*klmn-1:2*klmn,ispden) 559 end do 560 end if 561 end do 562 end if 563 pawrhoij(iatom)%nrhoijsel=nselect 564 ABI_DEALLOCATE(rhoijtmp) 565 end do 566 end if 567 ABI_DEALLOCATE(vpaw) 568 ABI_DEALLOCATE(vrespc) 569 570 !Eventually write the data on disk and deallocate f_fftgr_disk 571 call ab7_mixing_eval_deallocate(mix) 572 573 call timab(904,2,tsec) 574 575 !Restore potential 576 if (ispmix==1.and.nfft==nfftmix) then 577 vtrial=vtrial0 578 else if (nfft==nfftmix) then 579 do ispden=1,dtset%nspden 580 call fourdp(1,vtrial0(:,ispden),vtrial(:,ispden),+1,mpi_enreg,nfft,ngfft,dtset%paral_kgb,tim_fourdp) 581 end do 582 else 583 do ispden=1,dtset%nspden 584 do ifft=1,nfftmix 585 jfft=mixtofft(ifft) 586 vtrialg(1,jfft,ispden)=vtrial0(2*ifft-1,ispden) 587 vtrialg(2,jfft,ispden)=vtrial0(2*ifft ,ispden) 588 end do 589 call fourdp(1,vtrialg(:,:,ispden),vtrial(:,ispden),+1,mpi_enreg,nfft,ngfft,dtset%paral_kgb,tim_fourdp) 590 end do 591 ABI_DEALLOCATE(vtrialg) 592 end if 593 ABI_DEALLOCATE(vtrial0) 594 595 call timab(905,1,tsec) 596 597 !------Treat the mean of the potential 598 599 !Compute the mean of the new vtrial 600 call mean_fftr(vtrial,vmean,nfft,nfftot,dtset%nspden,mpi_comm_sphgrid) 601 602 !Reset the mean of the new vtrial, to the value vnew_mean 603 !When spin-polarized and fixed occupation numbers, 604 !treat separately spin up and spin down. 605 !Otherwise, use only global mean 606 do ispden=1,dtset%nspden 607 if (dtset%nspden==2.and.dtset%occopt>=3.and. & 608 & abs(dtset%spinmagntarget+99.99_dp)<1.0d-10)then 609 vme=(vnew_mean(1)+vnew_mean(2)-vmean(1)-vmean(2))*half 610 else 611 vme=vnew_mean(ispden)-vmean(ispden) 612 end if 613 vtrial(:,ispden)=vtrial(:,ispden)+vme 614 end do 615 616 if(moved_atm_inside==1 .and. istep/=nstep )then 617 if(abs(dtset%densfor_pred)==1.or.abs(dtset%densfor_pred)==4)then 618 ! Subtract current local psp, but also vxc (for core charges) 619 do ispden=1,dtset%nspden 620 vtrial(:,ispden)=vtrial(:,ispden)-vpsp(:)*identity(ispden)-vxc(:,ispden) 621 end do 622 else if(abs(dtset%densfor_pred)==2.or.abs(dtset%densfor_pred)==5.or.abs(dtset%densfor_pred)==6)then 623 ! Subtract current vpsp+Hxc from vtrial. This should be rationalized later 624 do ispden=1,dtset%nspden 625 vtrial(:,ispden)=vtrial(:,ispden)-(vpsp(:)+vhartr(:))*identity(ispden)-vxc(:,ispden) 626 end do 627 end if 628 end if 629 630 !In WVL: copy vtrial to BigDFT object: 631 if(dtset%usewvl==1) then 632 call wvl_vtrial_abi2big(1,vtrial,wvl%den) 633 end if 634 635 call timab(905,2,tsec) 636 call timab(93,2,tsec) 637 638 end subroutine newvtr