TABLE OF CONTENTS
ABINIT/dfpt_newvtr [ Functions ]
NAME
dfpt_newvtr
FUNCTION
Compute new first-order trial potential by mixing new and old values. First, compute preconditioned residual first-order potential. Then, call one of the self-consistency drivers, and 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 .
INPUTS
cplex: if 1, real space 1-order functions on FFT grid are REAL, if 2, COMPLEX dielar(7)=input parameters for dielectric matrix: diecut,dielng,diemac,diemix,diegap,dielam,diemixmag. dtset <type(dataset_type)>=all input variables in this dataset | isecur=level of security of the computation | 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 | paral_kgb=option for (kpt,g vectors,bands) parallelism | pawoptmix= - PAW only - 1 if the computed residuals include the PAW (rhoij) part etotal=the total energy obtained from the input vtrial ffttomix(nfft*(1-nfftmix/nfft))=Index of the points of the FFT (fine) grid on the grid used for mixing (coarse) initialized= if 0 the initialization of the RF run is not yet finished 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 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 mixtofft(nfftmix*(1-nfftmix/nfft))=Index of the points of the FFT grid used for mixing (coarse) on the FFT (fine) grid mpi_enreg=information about MPI parallelization my_natom=number of atoms treated by current processor 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 npawmix=-PAW only- number of spherical part elements to be mixed qphon(3)=reduced coordinates for the phonon wavelength (needed if cplex==2). pawrhoij(natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data Use here rhoij residuals (and gradients) rhor(cplex*nfft,nspden)=array for 1st-order electron density in electrons/bohr**3. rprimd(3,3)=dimensional primitive translations in real space (bohr) usepaw= 0 for non paw calculation; =1 for paw calculation vresid(cplex*nfft,nspden)=array for the residual of the potential xred(3,natom)=reduced dimensionless atomic coordinates
OUTPUT
dbl_nnsclo=1 if nnsclo has to be doubled to secure the convergence.
SIDE EFFECTS
vtrial(cplex*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 usepaw==1
pawrhoij(natom)%nrhoijsel,rhoijselect,rhoijp= several arrays
containing new values of rhoij (augmentation occupancies)
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
dfpt_scfcv
CHILDREN
ab7_mixing_copy_current_step,ab7_mixing_eval,ab7_mixing_eval_allocate ab7_mixing_eval_deallocate,fourdp,metric,moddiel,timab
SOURCE
94 #if defined HAVE_CONFIG_H 95 #include "config.h" 96 #endif 97 98 #include "abi_common.h" 99 100 subroutine dfpt_newvtr(cplex,dbl_nnsclo,dielar,dtset,etotal,ffttomix,& 101 & initialized,iscf,ispmix,istep,mix,mixtofft,& 102 & mpi_enreg,my_natom,nfft,nfftmix,ngfft,ngfftmix,npawmix,pawrhoij,& 103 & qphon,rhor,rprimd,usepaw,vresid,vtrial) 104 105 use defs_basis 106 use defs_abitypes 107 use m_profiling_abi 108 use m_ab7_mixing 109 use m_errors 110 111 use m_geometry, only : metric 112 use m_pawrhoij, only : pawrhoij_type 113 114 !This section has been created automatically by the script Abilint (TD). 115 !Do not modify the following lines by hand. 116 #undef ABI_FUNC 117 #define ABI_FUNC 'dfpt_newvtr' 118 use interfaces_18_timing 119 use interfaces_53_ffts 120 use interfaces_67_common 121 !End of the abilint section 122 123 implicit none 124 125 !Arguments------------------------------- 126 !scalars 127 integer,intent(in) :: cplex,initialized,iscf,ispmix,istep,my_natom,nfft 128 integer,intent(in) :: nfftmix,npawmix,usepaw 129 integer,intent(inout) :: dbl_nnsclo !vz_i 130 real(dp),intent(in) :: etotal 131 type(MPI_type),intent(in) :: mpi_enreg 132 type(ab7_mixing_object), intent(inout) :: mix 133 type(dataset_type),intent(in) :: dtset 134 !arrays 135 integer,intent(in) :: ffttomix(nfft*(1-nfftmix/nfft)) 136 integer,intent(in) :: mixtofft(nfftmix*(1-nfftmix/nfft)),ngfft(18) 137 integer,intent(in) :: ngfftmix(18) 138 real(dp),intent(in) :: dielar(7),qphon(3) 139 real(dp), intent(in), target :: rhor(cplex*nfft,dtset%nspden) 140 real(dp),intent(in) :: rprimd(3,3) 141 real(dp),intent(inout) :: vresid(cplex*nfft,dtset%nspden) 142 real(dp),intent(inout) :: vtrial(cplex*nfft,dtset%nspden) 143 type(pawrhoij_type),intent(inout) :: pawrhoij(my_natom*usepaw) 144 145 !Local variables------------------------------- 146 !scalars 147 integer :: cplex_mix,cplex_rhoij,dplex,i_vresid1,i_vrespc1,iatom,ifft,indx 148 integer :: irhoij,ispden,jfft,jrhoij,klmn,kmix,moved_atm_inside,nfftot,nselect 149 integer :: mpicomm,errid 150 logical :: mpi_summarize,reset 151 real(dp) :: fact,mixfac,mixfac_eff,mixfacmag,ucvol 152 character(len=500) :: msg 153 !arrays 154 real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3),tsec(2) 155 real(dp),allocatable :: rhoijrespc(:),rhoijtmp(:,:) 156 real(dp),allocatable :: vresid0(:,:),vrespc(:,:),vreswk(:,:) 157 real(dp), pointer :: vtrial0(:,:),vpaw(:) 158 real(dp),allocatable :: vtrialg(:,:,:) 159 160 ! ************************************************************************* 161 162 DBG_ENTER("COLL") 163 164 call timab(158,1,tsec) 165 166 !Compatibility tests 167 if(usepaw==1) then 168 if(dtset%nspden==4.and.dtset%pawoptmix==1) then 169 msg='pawoptmix=1 is not compatible with nspden=4 !' 170 MSG_ERROR(msg) 171 end if 172 if (my_natom>0) then 173 if (pawrhoij(1)%cplex<cplex) then 174 msg='pawrhoij()%cplex must be >=cplex !' 175 MSG_ERROR(msg) 176 end if 177 end if 178 end if 179 180 nfftot=ngfft(1)*ngfft(2)*ngfft(3) 181 cplex_mix=max(cplex,ispmix) 182 if (usepaw==1.and.my_natom>0) cplex_rhoij=pawrhoij(1)%cplex 183 184 !Compute different geometric tensor, as well as ucvol, from rprimd 185 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 186 moved_atm_inside=0 187 188 !Select components of potential to be mixed 189 ABI_ALLOCATE(vtrial0,(cplex_mix*nfftmix,dtset%nspden)) 190 ABI_ALLOCATE(vresid0,(cplex_mix*nfftmix,dtset%nspden)) 191 if (ispmix==1.and.nfft==nfftmix) then 192 vtrial0=vtrial;vresid0=vresid 193 else if (nfft==nfftmix) then 194 do ispden=1,dtset%nspden 195 call fourdp(cplex,vtrial0(:,ispden),vtrial(:,ispden),-1,mpi_enreg,nfft,ngfft,dtset%paral_kgb,0) 196 call fourdp(cplex,vresid0(:,ispden),vresid(:,ispden),-1,mpi_enreg,nfft,ngfft,dtset%paral_kgb,0) 197 end do 198 else 199 ABI_ALLOCATE(vtrialg,(2,nfft,dtset%nspden)) 200 ABI_ALLOCATE(vreswk,(2,nfft)) 201 do ispden=1,dtset%nspden 202 fact=dielar(4);if (ispden>1) fact=dielar(7) 203 call fourdp(cplex,vtrialg(:,:,ispden),vtrial(:,ispden),-1,mpi_enreg,nfft,ngfft,dtset%paral_kgb,0) 204 call fourdp(cplex,vreswk,vresid(:,ispden),-1,mpi_enreg,nfft,ngfft,dtset%paral_kgb,0) 205 do ifft=1,nfft 206 if (ffttomix(ifft)>0) then 207 jfft=2*ffttomix(ifft) 208 vtrial0(jfft-1,ispden)=vtrialg(1,ifft,ispden) 209 vtrial0(jfft ,ispden)=vtrialg(2,ifft,ispden) 210 vresid0(jfft-1,ispden)=vreswk(1,ifft) 211 vresid0(jfft ,ispden)=vreswk(2,ifft) 212 else 213 vtrialg(:,ifft,ispden)=vtrialg(:,ifft,ispden)+fact*vreswk(:,ifft) 214 end if 215 end do 216 end do 217 ABI_DEALLOCATE(vreswk) 218 end if 219 220 !Precondition the potential residual: 221 !Use a model dielectric function preconditioning, or simple mixing 222 ABI_ALLOCATE(vrespc,(cplex_mix*nfftmix,dtset%nspden)) 223 call moddiel(cplex_mix,dielar,mpi_enreg,nfftmix,ngfftmix,dtset%nspden,ispmix,0,& 224 & dtset%paral_kgb,qphon,rprimd,vresid0,vrespc) 225 226 !PAW only : precondition the rhoij quantities (augmentation occupancies) residuals. 227 !Use a simple preconditionning with the same mixing factor 228 !as the model dielectric function. 229 if (usepaw==1.and.my_natom>0) then 230 ABI_ALLOCATE(rhoijrespc,(npawmix)) 231 mixfac=dielar(4);mixfacmag=abs(dielar(7)) 232 if (cplex_rhoij==1) then 233 indx=0 234 do iatom=1,my_natom 235 do ispden=1,pawrhoij(iatom)%nspden 236 mixfac_eff=mixfac;if (ispden>1) mixfac_eff=mixfacmag 237 do kmix=1,pawrhoij(iatom)%lmnmix_sz 238 indx=indx+1;klmn=pawrhoij(iatom)%kpawmix(kmix) 239 rhoijrespc(indx)=mixfac_eff*pawrhoij(iatom)%rhoijres(klmn,ispden) 240 end do 241 end do 242 end do 243 else 244 indx=-1 245 do iatom=1,my_natom 246 do ispden=1,pawrhoij(iatom)%nspden 247 mixfac_eff=mixfac;if (ispden>1) mixfac_eff=mixfacmag 248 do kmix=1,pawrhoij(iatom)%lmnmix_sz 249 indx=indx+2;klmn=2*pawrhoij(iatom)%kpawmix(kmix)-1 250 rhoijrespc(indx:indx+1)=mixfac_eff*pawrhoij(iatom)%rhoijres(klmn:klmn+1,ispden) 251 end do 252 end do 253 end do 254 end if 255 end if 256 257 !------Compute new vtrial 258 259 i_vresid1=mix%i_vresid(1) 260 i_vrespc1=mix%i_vrespc(1) 261 262 !Initialise working arrays for the mixing object. 263 call ab7_mixing_eval_allocate(mix, istep) 264 265 !Copy current step arrays. 266 call ab7_mixing_copy_current_step(mix, vresid0, errid, msg, arr_respc = vrespc) 267 268 if (errid /= AB7_NO_ERROR) then 269 MSG_ERROR(msg) 270 end if 271 272 ABI_DEALLOCATE(vrespc) 273 ABI_DEALLOCATE(vresid0) 274 275 !PAW: either use the array f_paw or the array f_paw_disk 276 ABI_ALLOCATE(vpaw,(npawmix*usepaw)) 277 if (usepaw==1.and.my_natom>0) then 278 dplex=cplex_rhoij-1 279 indx=-dplex 280 do iatom=1,my_natom 281 do ispden=1,pawrhoij(iatom)%nspden 282 ABI_ALLOCATE(rhoijtmp,(cplex_rhoij*pawrhoij(iatom)%lmn2_size,1)) 283 rhoijtmp=zero 284 jrhoij=1 285 do irhoij=1,pawrhoij(iatom)%nrhoijsel 286 klmn=cplex_rhoij*pawrhoij(iatom)%rhoijselect(irhoij)-dplex 287 rhoijtmp(klmn:klmn+dplex,1)=pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+dplex,ispden) 288 jrhoij=jrhoij+cplex_rhoij 289 end do 290 do kmix=1,pawrhoij(iatom)%lmnmix_sz 291 indx=indx+cplex_rhoij;klmn=cplex_rhoij*pawrhoij(iatom)%kpawmix(kmix)-dplex 292 vpaw(indx:indx+dplex)=rhoijtmp(klmn:klmn+dplex,1)-pawrhoij(iatom)%rhoijres(klmn:klmn+dplex,ispden) 293 mix%f_paw(indx:indx+dplex,i_vresid1)=pawrhoij(iatom)%rhoijres(klmn:klmn+dplex,ispden) 294 mix%f_paw(indx:indx+dplex,i_vrespc1)=rhoijrespc(indx:indx+dplex) 295 end do 296 ABI_DEALLOCATE(rhoijtmp) 297 end do 298 end do 299 end if 300 301 !Unlike for GS, no need to modify the mean of vtrial 302 303 mpicomm=0;mpi_summarize=.false. 304 reset=.false.;if (initialized==0) reset=.true. 305 call ab7_mixing_eval(mix, vtrial0, istep, nfftot, ucvol, & 306 & mpicomm, mpi_summarize, errid, msg, & 307 & reset = reset, isecur = dtset%isecur, & 308 & pawopt = dtset%pawoptmix, response = 1, pawarr = vpaw, & 309 & etotal = etotal, potden = rhor, comm_atom=mpi_enreg%comm_atom) 310 311 if (errid == AB7_ERROR_MIXING_INC_NNSLOOP) then 312 dbl_nnsclo = 1 313 else if (errid /= AB7_NO_ERROR) then 314 ! MG FIXME, Why this? 315 ! One should propagate the error so that we can handle it 316 ! in the caller! 317 MSG_ERROR(msg) 318 end if 319 320 !Do here the mixing of the potential 321 if(iscf==2 .or. iscf==3 .or. iscf==7)then 322 ! PAW: restore rhoij from compact storage 323 if (usepaw==1.and.my_natom>0) then 324 dplex=cplex_rhoij-1 325 indx=-dplex 326 do iatom=1,my_natom 327 ABI_ALLOCATE(rhoijtmp,(cplex_rhoij*pawrhoij(iatom)%lmn2_size,pawrhoij(iatom)%nspden)) 328 rhoijtmp=zero 329 if (pawrhoij(iatom)%lmnmix_sz<pawrhoij(iatom)%lmn2_size) then 330 do ispden=1,pawrhoij(iatom)%nspden 331 jrhoij=1 332 do irhoij=1,pawrhoij(iatom)%nrhoijsel 333 klmn=cplex_rhoij*pawrhoij(iatom)%rhoijselect(irhoij)-dplex 334 rhoijtmp(klmn:klmn+dplex,ispden)=pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+dplex,ispden) 335 jrhoij=jrhoij+cplex_rhoij 336 end do 337 end do 338 end if 339 do ispden=1,pawrhoij(iatom)%nspden 340 do kmix=1,pawrhoij(iatom)%lmnmix_sz 341 indx=indx+cplex_rhoij;klmn=cplex_rhoij*pawrhoij(iatom)%kpawmix(kmix)-dplex 342 rhoijtmp(klmn:klmn+dplex,ispden)=vpaw(indx:indx+dplex) 343 end do 344 end do 345 nselect=0 346 do klmn=1,pawrhoij(iatom)%lmn2_size 347 if (any(abs(rhoijtmp(cplex_rhoij*klmn-dplex:cplex_rhoij*klmn,:))>tol10)) then 348 nselect=nselect+1 349 pawrhoij(iatom)%rhoijselect(nselect)=klmn 350 do ispden=1,pawrhoij(iatom)%nspden 351 pawrhoij(iatom)%rhoijp(cplex_rhoij*nselect-dplex:cplex_rhoij*nselect,ispden)=& 352 & rhoijtmp(cplex_rhoij*klmn-dplex:cplex_rhoij*klmn,ispden) 353 end do 354 end if 355 end do 356 pawrhoij(iatom)%nrhoijsel=nselect 357 ABI_DEALLOCATE(rhoijtmp) 358 end do 359 end if 360 361 else if(iscf==5 .or. iscf==6)then 362 if(ispmix/=1) then 363 msg = ' Mixing on reciprocal space not allowed with iscf=5 or 6.' 364 MSG_ERROR(msg) 365 end if 366 ! PAW: apply a simple mixing to rhoij (this is temporary) 367 if (usepaw==1.and.my_natom>0) then 368 indx=1-cplex_rhoij 369 do iatom=1,my_natom 370 ABI_ALLOCATE(rhoijtmp,(cplex_rhoij*pawrhoij(iatom)%lmn2_size,pawrhoij(iatom)%nspden)) 371 rhoijtmp=zero 372 if (pawrhoij(iatom)%lmnmix_sz<pawrhoij(iatom)%lmn2_size) then 373 do ispden=1,pawrhoij(iatom)%nspden 374 do kmix=1,pawrhoij(iatom)%lmnmix_sz 375 indx=indx+cplex_rhoij;klmn=cplex_rhoij*pawrhoij(iatom)%kpawmix(kmix)-dplex 376 rhoijtmp(klmn:klmn+dplex,ispden)=rhoijrespc(indx:indx+dplex) & 377 & -pawrhoij(iatom)%rhoijres(klmn:klmn+dplex,ispden) 378 end do 379 end do 380 end if 381 do ispden=1,pawrhoij(iatom)%nspden 382 jrhoij=1 383 do irhoij=1,pawrhoij(iatom)%nrhoijsel 384 klmn=cplex_rhoij*pawrhoij(iatom)%rhoijselect(irhoij)-dplex 385 rhoijtmp(klmn:klmn+dplex,ispden)=rhoijtmp(klmn:klmn+dplex,ispden) & 386 & +pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+dplex,ispden) 387 jrhoij=jrhoij+cplex_rhoij 388 end do 389 end do 390 nselect=0 391 do klmn=1,pawrhoij(iatom)%lmn2_size 392 if (any(abs(rhoijtmp(cplex_rhoij*klmn-dplex:cplex_rhoij*klmn,:))>tol10)) then 393 nselect=nselect+1 394 pawrhoij(iatom)%rhoijselect(nselect)=klmn 395 do ispden=1,pawrhoij(iatom)%nspden 396 pawrhoij(iatom)%rhoijp(cplex_rhoij*nselect-dplex:cplex_rhoij*nselect,ispden)=& 397 & rhoijtmp(cplex_rhoij*klmn-dplex:cplex_rhoij*klmn,ispden) 398 end do 399 end if 400 end do 401 pawrhoij(iatom)%nrhoijsel=nselect 402 ABI_DEALLOCATE(rhoijtmp) 403 end do 404 end if 405 end if 406 407 ABI_DEALLOCATE(vpaw) 408 if (usepaw==1.and.my_natom>0) then 409 ABI_DEALLOCATE(rhoijrespc) 410 end if 411 412 !Eventually write the data on disk and deallocate f_fftgr_disk 413 call ab7_mixing_eval_deallocate(mix) 414 415 !Restore potential 416 if (ispmix==1.and.nfft==nfftmix) then 417 vtrial=vtrial0 418 else if (nfft==nfftmix) then 419 do ispden=1,dtset%nspden 420 call fourdp(cplex,vtrial0(:,ispden),vtrial(:,ispden),+1,mpi_enreg,nfft,ngfft,dtset%paral_kgb,0) 421 end do 422 else 423 do ispden=1,dtset%nspden 424 do ifft=1,nfftmix 425 jfft=mixtofft(ifft) 426 vtrialg(1,jfft,ispden)=vtrial0(2*ifft-1,ispden) 427 vtrialg(2,jfft,ispden)=vtrial0(2*ifft ,ispden) 428 end do 429 call fourdp(cplex,vtrialg(:,:,ispden),vtrial(:,ispden),+1,mpi_enreg,nfft,ngfft,dtset%paral_kgb,0) 430 end do 431 ABI_DEALLOCATE(vtrialg) 432 end if 433 ABI_DEALLOCATE(vtrial0) 434 435 call timab(158,2,tsec) 436 437 DBG_ENTER("COLL") 438 439 end subroutine dfpt_newvtr