TABLE OF CONTENTS
ABINIT/m_newrho [ Modules ]
NAME
m_newrho
FUNCTION
COPYRIGHT
Copyright (C) 2005-2024 ABINIT group (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_newrho 23 24 use defs_basis 25 use defs_wvltypes 26 use m_errors 27 use m_abicore 28 use m_ab7_mixing 29 use m_abi2big 30 use m_dtset 31 32 use defs_datatypes, only : pseudopotential_type 33 use defs_abitypes, only : MPI_type 34 use m_time, only : timab 35 use m_geometry, only : metric 36 use m_pawtab, only : pawtab_type 37 use m_pawrhoij, only : pawrhoij_type,pawrhoij_filter 38 use m_prcref, only : prcref 39 use m_wvl_rho, only : wvl_prcref 40 use m_fft, only : fourdp 41 42 implicit none 43 44 private
ABINIT/newrho [ Functions ]
NAME
newrho
FUNCTION
Compute new trial density by mixing new and old values. Call prcref to compute preconditioned residual density and forces, Then, call one of the self-consistency drivers, then update density.
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 | densfor_pred= governs the preconditioning of the atomic charges | iprcel= governs the preconditioning of the density residual | iprcfc= governs the preconditioning of the forces | iscf=( <= 0 =>non-SCF), >0 => SCF) | iscf =11 => determination of the largest eigenvalue of the SCF cycle | iscf =12 => SCF cycle, simple mixing | iscf =13 => SCF cycle, Anderson mixing | iscf =14 => SCF cycle, Anderson mixing (order 2) | iscf =15 => SCF cycle, CG based on the minimization of the energy | iscf =17 => SCF cycle, Pulay mixing | 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 | pawoptmix=-PAW- 1 if the computed residuals include the PAW (rhoij) part | prtvol=control print volume and debugging etotal=the total energy obtained from the input density fnametmp_fft=name of _FFT file 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 dielectric matrix 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 density residual must be computed; otherwise, compute only the preconditioned density 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 nresid(nfft,nspden)=array for the residual of the density 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) pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data psps <type(pseudopotential_type)>=variables related to pseudopotentials 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 vtrial(nfft,nspden)=the trial potential that gave vresid. xred(3,natom)=reduced dimensionless atomic coordinates tauresid(nfft,nspden*dtset%usekden)=array for kinetic energy density residue (out - in)
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 mix<type(ab7_mixing_object)>=all data defining the mixing algorithm for the density rhor(nfft,nspden)= at input, it is the "out" trial density that gave nresid=(rho_out-rho_in) at output, it is an updated "mixed" trial density rhog(2,nfft)= Fourier transform of the new trial density ===== if usekden==1 ===== [mix_mgga<type(ab7_mixing_object)>]=all data defining the mixing algorithm for the kinetic energy density ===== 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=number of non-zero values of rhoij pawrhoij(natom)%rhoijp(cplex_rhoij*lmn2_size,nspden)= new (mixed) value of rhoij quantities in PACKED STORAGE pawrhoij(natom)%rhoijselect(lmn2_size)=select the non-zero values of rhoij taug(2,nfft*dtset%usekden)=array for Fourier transform of kinetic energy density taur(nfft,nspden*dtset%usekden)=array for kinetic energy density
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 (except f_fftgr). In case of norm-conserving calculations the FFT grid is the usual FFT grid.
SOURCE
165 subroutine newrho(atindx,dbl_nnsclo,dielar,dielinv,dielstrt,dtn_pc,dtset,etotal,fcart,ffttomix,& 166 & gmet,grhf,gsqcut,initialized,ispmix,istep,kg_diel,kxc,mgfft,mix,mixtofft,& 167 & moved_atm_inside,mpi_enreg,my_natom,nattyp,nfft,& 168 & nfftmix,nfftmix_per_nfft,ngfft,ngfftmix,nkxc,npawmix,npwdiel,& 169 & nresid,ntypat,n1xccc,pawrhoij,pawtab,& 170 & ph1d,psps,rhog,rhor,rprimd,susmat,usepaw,vtrial,wvl,wvl_den,xred,& 171 & mix_mgga,taug,taur,tauresid) 172 173 !Arguments------------------------------- 174 !scalars 175 integer,intent(in) :: dielstrt,initialized,ispmix,istep,my_natom,mgfft 176 integer,intent(in) :: moved_atm_inside,n1xccc,nfft 177 integer,intent(in) :: nfftmix,nfftmix_per_nfft 178 integer,intent(in) :: nkxc,npawmix,npwdiel,ntypat,usepaw 179 integer,intent(inout) :: dbl_nnsclo 180 real(dp),intent(in) :: etotal,gsqcut 181 type(MPI_type),intent(in) :: mpi_enreg 182 type(ab7_mixing_object), intent(inout) :: mix 183 type(ab7_mixing_object), intent(inout),optional :: mix_mgga 184 type(dataset_type),intent(in) :: dtset 185 type(pseudopotential_type),intent(in) :: psps 186 type(wvl_internal_type), intent(in) :: wvl 187 type(wvl_denspot_type), intent(inout) :: wvl_den 188 !arrays 189 integer,intent(in) :: atindx(dtset%natom) 190 integer,intent(in) :: ffttomix(nfft*(nfftmix_per_nfft)) 191 integer,intent(in) :: kg_diel(3,npwdiel) 192 integer,intent(in) :: mixtofft(nfftmix*nfftmix_per_nfft) 193 integer,intent(in) :: nattyp(ntypat),ngfft(18),ngfftmix(18) 194 real(dp),intent(in) :: dielar(7),fcart(3,dtset%natom),grhf(3,dtset%natom) 195 real(dp),intent(inout) :: rprimd(3,3) 196 real(dp),intent(in) :: susmat(2,npwdiel,dtset%nspden,npwdiel,dtset%nspden) 197 real(dp),intent(in), target :: vtrial(nfft,dtset%nspden) 198 real(dp),intent(inout) :: dielinv(2,npwdiel,dtset%nspden,npwdiel,dtset%nspden) 199 real(dp),intent(inout), target :: dtn_pc(3,dtset%natom) 200 real(dp),intent(inout) :: gmet(3,3) 201 !TODO: nresid appears to be only intent in here. 202 real(dp),intent(inout) :: kxc(nfft,nkxc),nresid(nfft,dtset%nspden) 203 real(dp),intent(inout) :: ph1d(2,3*(2*mgfft+1)*dtset%natom) 204 real(dp),intent(inout) :: rhor(nfft,dtset%nspden) 205 real(dp),intent(inout), target :: xred(3,dtset%natom) 206 real(dp),intent(inout) :: rhog(2,nfft) 207 real(dp),intent(inout), optional :: taug(2,nfft*dtset%usekden) 208 real(dp),intent(inout), optional :: taur(nfft,dtset%nspden*dtset%usekden) 209 real(dp),intent(inout), optional :: tauresid(nfft,dtset%nspden*dtset%usekden) 210 211 type(pawrhoij_type),intent(inout) :: pawrhoij(my_natom*psps%usepaw) 212 type(pawtab_type),intent(in) :: pawtab(ntypat*psps%usepaw) 213 214 !Local variables------------------------------- 215 !scalars 216 integer,parameter :: tim_fourdp9=9 217 integer :: cplex,dplex,errid,i_vresid1,i_vrespc1,iatom,ifft,indx,iq,iq0,irhoij,ispden,jfft 218 integer :: jrhoij,klmn,kklmn,kmix,mpicomm,nfftot,qphase 219 logical :: mpi_summarize,reset 220 real(dp) :: fact,ucvol,ucvol_local 221 character(len=500) :: message 222 !arrays 223 real(dp) :: gprimd(3,3),rmet(3,3),ro(2),tsec(2),vhartr_dum(1),vpsp_dum(1) 224 real(dp) :: vxc_dum(1,1) 225 real(dp),allocatable :: magng(:,:,:),magntaug(:,:,:) 226 real(dp),allocatable :: nresid0(:,:),nrespc(:,:),nreswk(:,:,:) 227 real(dp),allocatable :: rhoijrespc(:),rhoijtmp(:,:) 228 ! TODO : these should be allocatables not pointers: is there some reason to 229 ! keep them this way, eg an interface somewhere? 230 real(dp), pointer :: rhomag(:,:), npaw(:) 231 real(dp),allocatable :: tauresid0(:,:),taurespc(:,:) 232 real(dp),allocatable :: taumag(:,:) 233 234 ! ************************************************************************* 235 236 DBG_ENTER("COLL") 237 call timab(94,1,tsec) 238 239 nfftot=PRODUCT(ngfft(1:3)) 240 241 !Compatibility tests 242 if(nfftmix>nfft) then 243 message='nfftmix>nfft not allowed!' 244 ABI_BUG(message) 245 end if 246 247 if(dtset%usewvl==1) then 248 if( (ispmix/=1 .or. nfftmix/=nfft)) then 249 message='nfftmix/=nfft, ispmix/=1 not allowed for wavelets!' 250 ABI_BUG(message) 251 end if 252 if(dtset%wvl_bigdft_comp==1) then 253 message='usewvl == 1 and wvl_bigdft_comp==1 not allowed!' 254 ABI_BUG(message) 255 end if 256 end if 257 258 if(ispmix/=2.and.nfftmix/=nfft) then 259 message='nfftmix/=nfft allowed only when ispmix=2!' 260 ABI_BUG(message) 261 end if 262 263 if (dtset%usekden==1) then 264 if ((.not.present(tauresid)).or.(.not.present(taug)).or. & 265 & (.not.present(taur)).or.(.not.present(mix_mgga))) then 266 message='Several arrays are missing!' 267 ABI_BUG(message) 268 end if 269 if (mix_mgga%iscf==AB7_MIXING_CG_ENERGY.or.mix_mgga%iscf==AB7_MIXING_CG_ENERGY_2.or.& 270 & mix_mgga%iscf==AB7_MIXING_EIG) then 271 message='kinetic energy density cannot be mixed with the selected mixing algorithm!' 272 ABI_ERROR(message) 273 end if 274 end if 275 276 if (usepaw==1.and.my_natom>0) then 277 cplex=pawrhoij(1)%cplex_rhoij;dplex=cplex-1 278 qphase=pawrhoij(1)%qphase 279 else 280 cplex = 0;dplex = 0 ; qphase=0 281 end if 282 283 !Compute different geometric tensor, as well as ucvol, from rprimd 284 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 285 286 if(dtset%usewvl==0) then 287 ucvol_local=ucvol 288 #if defined HAVE_BIGDFT 289 else 290 ucvol_local = product(wvl_den%denspot%dpbox%hgrids) * real(nfftot, dp) 291 #endif 292 end if 293 294 !Select components of density to be mixed 295 ABI_MALLOC(rhomag,(ispmix*nfftmix,dtset%nspden)) 296 ABI_MALLOC(nresid0,(ispmix*nfftmix,dtset%nspden)) 297 ABI_MALLOC(taumag,(ispmix*nfftmix,dtset%nspden*dtset%usekden)) 298 ABI_MALLOC(tauresid0,(ispmix*nfftmix,dtset%nspden*dtset%usekden)) 299 ! real space and all fft points are here 300 if (ispmix==1.and.nfft==nfftmix) then 301 rhomag(:,1:dtset%nspden)=rhor(:,1:dtset%nspden) 302 nresid0(:,1:dtset%nspden)=nresid(:,1:dtset%nspden) 303 if (dtset%usekden==1) then 304 taumag(:,1:dtset%nspden)=taur(:,1:dtset%nspden) 305 tauresid0(:,1:dtset%nspden)=tauresid(:,1:dtset%nspden) 306 end if 307 ! recip space and all fft points are here 308 else if (nfft==nfftmix) then 309 do ispden=1,dtset%nspden 310 call fourdp(1,nresid0(:,ispden),nresid(:,ispden),-1,mpi_enreg,nfft,1,ngfft,tim_fourdp9) 311 end do 312 rhomag(:,1)=reshape(rhog,(/2*nfft/)) 313 if (dtset%nspden>1) then 314 do ispden=2,dtset%nspden 315 call fourdp(1,rhomag(:,ispden),rhor(:,ispden),-1,mpi_enreg,nfft,1,ngfft,tim_fourdp9) 316 end do 317 end if 318 if (dtset%usekden==1) then 319 do ispden=1,dtset%nspden 320 call fourdp(1,tauresid0(:,ispden),tauresid(:,ispden),-1,mpi_enreg,nfft,1,ngfft,tim_fourdp9) 321 end do 322 taumag(:,1)=reshape(taug,(/2*nfft/)) 323 if (dtset%nspden>1) then 324 do ispden=2,dtset%nspden 325 call fourdp(1,taumag(:,ispden),taur(:,ispden),-1,mpi_enreg,nfft,1,ngfft,tim_fourdp9) 326 end do 327 end if 328 end if 329 ! not all fft points are here - presumes recip space 330 else 331 fact=dielar(4)-1._dp 332 ABI_MALLOC(nreswk,(2,nfft,dtset%nspden)) 333 do ispden=1,dtset%nspden 334 call fourdp(1,nreswk(:,:,ispden),nresid(:,ispden),-1,mpi_enreg,nfft,1,ngfft,tim_fourdp9) 335 end do 336 do ifft=1,nfft 337 if (ffttomix(ifft)>0) then 338 jfft=2*ffttomix(ifft) 339 rhomag (jfft-1:jfft,1)=rhog(1:2,ifft) 340 nresid0(jfft-1:jfft,1)=nreswk(1:2,ifft,1) 341 else 342 rhog(:,ifft)=rhog(:,ifft)+fact*nreswk(:,ifft,1) 343 end if 344 end do 345 if (dtset%nspden>1) then 346 ABI_MALLOC(magng,(2,nfft,dtset%nspden-1)) 347 do ispden=2,dtset%nspden 348 call fourdp(1,magng(:,:,ispden-1),rhor(:,ispden),-1,mpi_enreg,nfft,1,ngfft,tim_fourdp9) 349 do ifft=1,nfft 350 if (ffttomix(ifft)>0) then 351 jfft=2*ffttomix(ifft) 352 rhomag (jfft-1:jfft,ispden)=magng (1:2,ifft,ispden-1) 353 nresid0(jfft-1:jfft,ispden)=nreswk(1:2,ifft,ispden) 354 else 355 magng(:,ifft,ispden-1)=magng(:,ifft,ispden-1)+fact*nreswk(:,ifft,ispden) 356 if (dtset%nspden==2) magng(:,ifft,1)=two*magng(:,ifft,1)-rhog(:,ifft) 357 end if 358 end do 359 end do 360 end if 361 if (dtset%usekden==1) then 362 do ispden=1,dtset%nspden 363 call fourdp(1,nreswk(:,:,ispden),tauresid(:,ispden),-1,mpi_enreg,nfft,1,ngfft,tim_fourdp9) 364 end do 365 do ifft=1,nfft 366 if (ffttomix(ifft)>0) then 367 jfft=2*ffttomix(ifft) 368 taumag (jfft-1:jfft,1)=taug(1:2,ifft) 369 tauresid0(jfft-1:jfft,1)=nreswk(1:2,ifft,1) 370 else 371 taug(:,ifft)=taug(:,ifft)+fact*nreswk(:,ifft,1) 372 end if 373 end do 374 if (dtset%nspden>1) then 375 ABI_MALLOC(magntaug,(2,nfft,dtset%nspden-1)) 376 do ispden=2,dtset%nspden 377 call fourdp(1,magntaug(:,:,ispden-1),taur(:,ispden),-1,mpi_enreg,nfft,1,ngfft,tim_fourdp9) 378 do ifft=1,nfft 379 if (ffttomix(ifft)>0) then 380 jfft=2*ffttomix(ifft) 381 taumag (jfft-1:jfft,ispden)=magntaug(1:2,ifft,ispden-1) 382 tauresid0(jfft-1:jfft,ispden)=nreswk(1:2,ifft,ispden) 383 else 384 magntaug(:,ifft,ispden-1)=magntaug(:,ifft,ispden-1)+fact*nreswk(:,ifft,ispden) 385 if (dtset%nspden==2) magntaug(:,ifft,1)=two*magntaug(:,ifft,1)-taug(:,ifft) 386 end if 387 end do 388 end do 389 end if 390 end if 391 ABI_FREE(nreswk) 392 end if 393 394 !Retrieve "input" density from "output" density and density residual 395 rhomag(:,1:dtset%nspden)=rhomag(:,1:dtset%nspden)-nresid0(:,1:dtset%nspden) 396 if (dtset%usekden==1) then 397 taumag(:,1:dtset%nspden)=taumag(:,1:dtset%nspden)-tauresid0(:,1:dtset%nspden) 398 end if 399 400 !If nspden==2, separate density and magnetization 401 if (dtset%nspden==2) then 402 rhomag (:,2)=two*rhomag (:,2)-rhomag (:,1) 403 nresid0(:,2)=two*nresid0(:,2)-nresid0(:,1) 404 if (dtset%usekden==1) then 405 taumag (:,2)=two*taumag (:,2)-taumag (:,1) 406 tauresid0(:,2)=two*tauresid0(:,2)-tauresid0(:,1) 407 end if 408 end if 409 410 !If PAW, handle occupancy matrix 411 if (usepaw==1.and.my_natom>0) then 412 if (pawrhoij(1)%nspden==2) then 413 do iatom=1,my_natom 414 do iq=1,qphase 415 iq0=merge(0,cplex*pawrhoij(iatom)%lmn2_size,iq==1) 416 jrhoij=1+iq0 417 do irhoij=1,pawrhoij(iatom)%nrhoijsel 418 ro(1:1+dplex)=pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+dplex,1) 419 pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+dplex,1)=ro(1:1+dplex)+pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+dplex,2) 420 pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+dplex,2)=ro(1:1+dplex)-pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+dplex,2) 421 jrhoij=jrhoij+cplex 422 end do 423 do kmix=1,pawrhoij(iatom)%lmnmix_sz 424 klmn=iq0+cplex*pawrhoij(iatom)%kpawmix(kmix)-dplex 425 ro(1:1+dplex)=pawrhoij(iatom)%rhoijres(klmn:klmn+dplex,1) 426 pawrhoij(iatom)%rhoijres(klmn:klmn+dplex,1)=ro(1:1+dplex)+pawrhoij(iatom)%rhoijres(klmn:klmn+dplex,2) 427 pawrhoij(iatom)%rhoijres(klmn:klmn+dplex,2)=ro(1:1+dplex)-pawrhoij(iatom)%rhoijres(klmn:klmn+dplex,2) 428 end do 429 end do 430 end do 431 end if 432 end if 433 434 !Choice of preconditioner governed by iprcel, densfor_pred and iprcfc 435 ABI_MALLOC(nrespc,(ispmix*nfftmix,dtset%nspden)) 436 ABI_MALLOC(taurespc,(ispmix*nfftmix,dtset%nspden*dtset%usekden)) 437 ABI_MALLOC(npaw,(npawmix*usepaw)) 438 if (usepaw==1) then 439 ABI_MALLOC(rhoijrespc,(npawmix)) 440 else 441 ABI_MALLOC(rhoijrespc,(0)) 442 end if 443 if(dtset%usewvl==0) then 444 call prcref(atindx,dielar,dielinv,& 445 & dielstrt,dtn_pc,dtset,etotal,fcart,ffttomix,gmet,gsqcut,& 446 & istep,kg_diel,kxc,& 447 & mgfft,moved_atm_inside,mpi_enreg,my_natom,& 448 & nattyp,nfft,nfftmix,ngfft,ngfftmix,nkxc,npawmix,npwdiel,ntypat,n1xccc,& 449 & ispmix,1,pawrhoij,pawtab,ph1d,psps,rhog,rhoijrespc,rhor,rprimd,& 450 & susmat,vhartr_dum,vpsp_dum,nresid0,nrespc,vxc_dum,wvl,wvl_den,xred) 451 else 452 call wvl_prcref(dielar,dtset%iprcel,my_natom,nfftmix,npawmix,dtset%nspden,pawrhoij,& 453 & rhoijrespc,psps%usepaw,nresid0,nrespc) 454 end if 455 !At present, only a simple precoditionning for the kinetic energy density 456 ! (is Kerker mixing valid for tau?) 457 if (dtset%usekden==1) then 458 do ispden=1,dtset%nspden 459 fact=dielar(4);if (ispden>1) fact=abs(dielar(7)) 460 taurespc(1:ispmix*nfftmix,ispden)=fact*tauresid0(1:ispmix*nfftmix,ispden) 461 end do 462 end if 463 464 !------Compute new trial density and eventual new atomic positions 465 466 if (mix%n_fftgr>0) then 467 i_vresid1=mix%i_vresid(1) 468 i_vrespc1=mix%i_vrespc(1) 469 end if 470 471 !Initialise working arrays for the mixing object. 472 if (moved_atm_inside == 1) then 473 call ab7_mixing_use_moving_atoms(mix, dtset%natom, xred, dtn_pc) 474 end if 475 call ab7_mixing_eval_allocate(mix, istep) 476 477 !Copy current step arrays. 478 if (moved_atm_inside == 1) then 479 call ab7_mixing_copy_current_step(mix, nresid0, errid, message, arr_respc = nrespc, arr_atm = grhf) 480 else 481 call ab7_mixing_copy_current_step(mix, nresid0, errid, message, arr_respc = nrespc) 482 end if 483 if (errid /= AB7_NO_ERROR) then 484 ABI_ERROR(message) 485 end if 486 487 !Same treatment for the kinetic energy density 488 if (dtset%usekden==1) then 489 call ab7_mixing_eval_allocate(mix_mgga, istep) 490 call ab7_mixing_copy_current_step(mix_mgga, tauresid0, errid, message, arr_respc = taurespc) 491 if (errid /= AB7_NO_ERROR) then 492 ABI_ERROR(message) 493 end if 494 end if 495 496 ABI_FREE(nresid0) 497 ABI_FREE(nrespc) 498 ABI_FREE(tauresid0) 499 ABI_FREE(taurespc) 500 501 !PAW: either use the array f_paw or the array f_paw_disk 502 if (usepaw==1) then 503 indx=-dplex 504 do iatom=1,my_natom 505 ABI_MALLOC(rhoijtmp,(cplex*pawrhoij(iatom)%lmn2_size,1)) 506 do iq=1,qphase 507 iq0=merge(0,cplex*pawrhoij(iatom)%lmn2_size,iq==1) 508 do ispden=1,pawrhoij(iatom)%nspden 509 rhoijtmp=zero ; jrhoij=1+iq0 510 do irhoij=1,pawrhoij(iatom)%nrhoijsel 511 klmn=cplex*pawrhoij(iatom)%rhoijselect(irhoij)-dplex 512 rhoijtmp(klmn:klmn+dplex,1)=pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+dplex,ispden) 513 jrhoij=jrhoij+cplex 514 end do 515 do kmix=1,pawrhoij(iatom)%lmnmix_sz 516 indx=indx+cplex;klmn=cplex*pawrhoij(iatom)%kpawmix(kmix)-dplex ; kklmn=klmn+iq0 517 npaw(indx:indx+dplex)=rhoijtmp(klmn:klmn+dplex,1)-pawrhoij(iatom)%rhoijres(kklmn:kklmn+dplex,ispden) 518 mix%f_paw(indx:indx+dplex,i_vresid1)=pawrhoij(iatom)%rhoijres(kklmn:kklmn+dplex,ispden) 519 mix%f_paw(indx:indx+dplex,i_vrespc1)=rhoijrespc(indx:indx+dplex) 520 end do 521 end do 522 end do 523 ABI_FREE(rhoijtmp) 524 end do 525 end if 526 527 !------Prediction of the components of the density 528 529 !Init mpicomm 530 if(mpi_enreg%paral_kgb==1)then 531 mpicomm=mpi_enreg%comm_fft 532 mpi_summarize=.true. 533 else 534 mpicomm=0 535 mpi_summarize=.false. 536 end if 537 if(dtset%usewvl==1) then 538 mpicomm=mpi_enreg%comm_wvl 539 mpi_summarize=(mpi_enreg%nproc_wvl > 1) 540 end if 541 542 reset = .false. 543 if (initialized == 0) reset = .true. 544 545 !Electronic density mixing 546 call ab7_mixing_eval(mix, rhomag, istep, nfftot, ucvol_local, & 547 & mpicomm, mpi_summarize, errid, message, & 548 & reset = reset, isecur = dtset%isecur,& 549 & pawopt = dtset%pawoptmix, pawarr = npaw, & 550 & etotal = etotal, potden = vtrial, & 551 & comm_atom=mpi_enreg%comm_atom) 552 if (errid == AB7_ERROR_MIXING_INC_NNSLOOP) then 553 dbl_nnsclo = 1 554 else if (errid /= AB7_NO_ERROR) then 555 ABI_ERROR(message) 556 end if 557 !Kinetic energy density mixing (if any) 558 if (dtset%usekden==1) then 559 call ab7_mixing_eval(mix_mgga, taumag, istep, nfftot, ucvol_local, & 560 & mpicomm, mpi_summarize, errid, message, reset = reset) 561 if (errid /= AB7_NO_ERROR) then 562 ABI_ERROR(message) 563 end if 564 end if 565 566 !PAW: apply a simple mixing to rhoij (this is temporary) 567 if(dtset%iscf==15 .or. dtset%iscf==16)then 568 if (usepaw==1) then 569 indx=-dplex 570 do iatom=1,my_natom 571 ABI_MALLOC(rhoijtmp,(cplex*qphase*pawrhoij(iatom)%lmn2_size,pawrhoij(iatom)%nspden)) 572 rhoijtmp=zero 573 do iq=1,qphase 574 iq0=merge(0,cplex*pawrhoij(iatom)%lmn2_size,iq==1) 575 if (pawrhoij(iatom)%lmnmix_sz<pawrhoij(iatom)%lmn2_size) then 576 do ispden=1,pawrhoij(iatom)%nspden 577 do kmix=1,pawrhoij(iatom)%lmnmix_sz 578 indx=indx+cplex;klmn=iq0+cplex*pawrhoij(iatom)%kpawmix(kmix)-dplex 579 rhoijtmp(klmn:klmn+dplex,ispden)=rhoijrespc(indx:indx+dplex) & 580 & -pawrhoij(iatom)%rhoijres(klmn:klmn+dplex,ispden) 581 end do 582 end do 583 end if 584 if (pawrhoij(iatom)%nspden/=2) then 585 do ispden=1,pawrhoij(iatom)%nspden 586 jrhoij=iq0+1 587 do irhoij=1,pawrhoij(iatom)%nrhoijsel 588 klmn=iq0+cplex*pawrhoij(iatom)%rhoijselect(irhoij)-dplex 589 rhoijtmp(klmn:klmn+dplex,ispden)=rhoijtmp(klmn:klmn+dplex,ispden) & 590 & +pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+dplex,ispden) 591 jrhoij=jrhoij+cplex 592 end do 593 end do 594 else 595 jrhoij=iq0+1 596 do irhoij=1,pawrhoij(iatom)%nrhoijsel 597 klmn=iq0+cplex*pawrhoij(iatom)%rhoijselect(irhoij)-dplex 598 ro(1:1+dplex)=rhoijtmp(klmn:klmn+dplex,1) 599 rhoijtmp(klmn:klmn+dplex,1)=half*(ro(1:1+dplex)+rhoijtmp(klmn:klmn+dplex,2)) & 600 & +pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+dplex,1) 601 rhoijtmp(klmn:klmn+dplex,2)=half*(ro(1:1+dplex)-rhoijtmp(klmn:klmn+dplex,2)) & 602 & +pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+dplex,2) 603 jrhoij=jrhoij+cplex 604 end do 605 end if 606 end do 607 call pawrhoij_filter(pawrhoij(iatom)%rhoijp,pawrhoij(iatom)%rhoijselect,& 608 & pawrhoij(iatom)%nrhoijsel,pawrhoij(iatom)%cplex_rhoij,pawrhoij(iatom)%qphase,& 609 & pawrhoij(iatom)%lmn2_size,pawrhoij(iatom)%nspden,rhoij_input=rhoijtmp) 610 ABI_FREE(rhoijtmp) 611 end do 612 end if 613 end if 614 615 !if (usepaw==1) then 616 ABI_FREE(rhoijrespc) 617 !end if 618 619 !PAW: restore rhoij from compact storage 620 if (usepaw==1.and.dtset%iscf/=15.and.dtset%iscf/=16) then 621 indx=-dplex 622 do iatom=1,my_natom 623 ABI_MALLOC(rhoijtmp,(cplex*qphase*pawrhoij(iatom)%lmn2_size,pawrhoij(iatom)%nspden)) 624 rhoijtmp=zero 625 do iq=1,qphase 626 iq0=merge(0,cplex*pawrhoij(iatom)%lmn2_size,iq==1) 627 if (pawrhoij(iatom)%lmnmix_sz<pawrhoij(iatom)%lmn2_size) then 628 do ispden=1,pawrhoij(iatom)%nspden 629 jrhoij=iq0+1 630 do irhoij=1,pawrhoij(iatom)%nrhoijsel 631 klmn=iq0+cplex*pawrhoij(iatom)%rhoijselect(irhoij)-dplex 632 rhoijtmp(klmn:klmn+dplex,ispden)=pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+dplex,ispden) 633 jrhoij=jrhoij+cplex 634 end do 635 end do 636 end if 637 do ispden=1,pawrhoij(iatom)%nspden 638 do kmix=1,pawrhoij(iatom)%lmnmix_sz 639 indx=indx+cplex;klmn=iq0+cplex*pawrhoij(iatom)%kpawmix(kmix)-dplex 640 rhoijtmp(klmn:klmn+dplex,ispden)=npaw(indx:indx+dplex) 641 end do 642 end do 643 if (pawrhoij(iatom)%nspden==2) then 644 jrhoij=iq0+1 645 do irhoij=1,pawrhoij(iatom)%nrhoijsel 646 klmn=iq0+cplex*pawrhoij(iatom)%rhoijselect(irhoij)-dplex 647 ro(1:1+dplex)=rhoijtmp(klmn:klmn+dplex,1) 648 rhoijtmp(klmn:klmn+dplex,1)=half*(ro(1:1+dplex)+rhoijtmp(klmn:klmn+dplex,2)) 649 rhoijtmp(klmn:klmn+dplex,2)=half*(ro(1:1+dplex)-rhoijtmp(klmn:klmn+dplex,2)) 650 jrhoij=jrhoij+cplex 651 end do 652 end if 653 end do 654 call pawrhoij_filter(pawrhoij(iatom)%rhoijp,pawrhoij(iatom)%rhoijselect,& 655 & pawrhoij(iatom)%nrhoijsel,pawrhoij(iatom)%cplex_rhoij,pawrhoij(iatom)%qphase,& 656 & pawrhoij(iatom)%lmn2_size,pawrhoij(iatom)%nspden,rhoij_input=rhoijtmp) 657 ABI_FREE(rhoijtmp) 658 end do 659 end if ! usepaw==1.and.dtset%iscf/=15.and.dtset%iscf/=16 660 ABI_FREE(npaw) 661 662 !Eventually write the data on disk and deallocate f_fftgr_disk 663 call ab7_mixing_eval_deallocate(mix) 664 if (dtset%usekden==1) call ab7_mixing_eval_deallocate(mix_mgga) 665 666 !Fourier transform the density 667 if (ispmix==1.and.nfft==nfftmix) then 668 !Real space mixing, no need to transform rhomag 669 rhor(:,1:dtset%nspden)=rhomag(:,1:dtset%nspden) 670 if(dtset%usewvl==0) then 671 !Get rhog from rhor(:,1) 672 call fourdp(1,rhog,rhor(:,1),-1,mpi_enreg,nfft,1,ngfft,tim_fourdp9) 673 end if 674 if (dtset%usekden==1) then 675 taur(:,1:dtset%nspden)=taumag(:,1:dtset%nspden) 676 if(dtset%usewvl==0) then 677 call fourdp(1,taug,taur(:,1),-1,mpi_enreg,nfft,1,ngfft,tim_fourdp9) 678 end if 679 end if 680 else if (nfft==nfftmix) then 681 !Reciprocal mixing space mixing, need to generate rhor in real space from rhomag in reciprocal space 682 do ispden=1,dtset%nspden 683 call fourdp(1,rhomag(:,ispden),rhor(:,ispden),+1,mpi_enreg,nfft,1,ngfft,tim_fourdp9) 684 end do 685 rhog(:,:)=reshape(rhomag(:,1),(/2,nfft/)) 686 if (dtset%usekden==1) then 687 do ispden=1,dtset%nspden 688 call fourdp(1,taumag(:,ispden),taur(:,ispden),+1,mpi_enreg,nfft,1,ngfft,tim_fourdp9) 689 end do 690 taug(:,:)=reshape(taumag(:,1),(/2,nfft/)) 691 end if 692 else 693 do ifft=1,nfftmix 694 jfft=mixtofft(ifft) 695 rhog(1:2,jfft)=rhomag(2*ifft-1:2*ifft,1) 696 end do 697 call fourdp(1,rhog,rhor(:,1),+1,mpi_enreg,nfft,1,ngfft,tim_fourdp9) 698 if (dtset%nspden>1) then 699 do ispden=2,dtset%nspden 700 do ifft=1,nfftmix 701 jfft=mixtofft(ifft) 702 magng(1:2,jfft,ispden-1)=rhomag(2*ifft-1:2*ifft,ispden) 703 end do 704 call fourdp(1,magng(:,:,ispden-1),rhor(:,ispden),+1,mpi_enreg,nfft,1,ngfft,tim_fourdp9) 705 end do 706 ABI_FREE(magng) 707 end if 708 if (dtset%usekden==1) then 709 do ifft=1,nfftmix 710 jfft=mixtofft(ifft) 711 taug(1:2,jfft)=taumag(2*ifft-1:2*ifft,1) 712 end do 713 call fourdp(1,taug,taur(:,1),+1,mpi_enreg,nfft,1,ngfft,tim_fourdp9) 714 if (dtset%nspden>1) then 715 do ispden=2,dtset%nspden 716 do ifft=1,nfftmix 717 jfft=mixtofft(ifft) 718 magntaug(1:2,jfft,ispden-1)=taumag(2*ifft-1:2*ifft,ispden) 719 end do 720 call fourdp(1,magntaug(:,:,ispden-1),taur(:,ispden),+1,mpi_enreg,nfft,1,ngfft,tim_fourdp9) 721 end do 722 ABI_FREE(magntaug) 723 end if 724 end if 725 end if 726 ABI_FREE(rhomag) 727 ABI_FREE(taumag) 728 729 !Set back rho in (up+dn,up) form if nspden=2 730 if (dtset%nspden==2) then 731 rhor(:,2)=half*(rhor(:,1)+rhor(:,2)) 732 if (dtset%usekden==1) taur(:,2)=half*(taur(:,1)+taur(:,2)) 733 end if 734 735 !In WVL: copy density to BigDFT object: 736 if(dtset%usewvl==1) then 737 call wvl_rho_abi2big(1,rhor,wvl_den) 738 end if 739 740 call timab(94,2,tsec) 741 742 DBG_EXIT("COLL") 743 744 end subroutine newrho