TABLE OF CONTENTS
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.
COPYRIGHT
Copyright (C) 2005-2018 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 . 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 | 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
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
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 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*lmn2_size,nspden)= new (mixed) value of rhoij quantities in PACKED STORAGE
pawrhoij(natom)%rhoijselect(lmn2_size)=select the non-zero values of rhoij
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.
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,metric prcref,timab,wvl_prcref,wvl_rho_abi2big
SOURCE
122 #if defined HAVE_CONFIG_H 123 #include "config.h" 124 #endif 125 126 #include "abi_common.h" 127 128 subroutine newrho(atindx,dbl_nnsclo,dielar,dielinv,dielstrt,dtn_pc,dtset,etotal,fcart,ffttomix,& 129 & gmet,grhf,gsqcut,initialized,ispmix,istep,kg_diel,kxc,mgfft,mix,mixtofft,& 130 & moved_atm_inside,mpi_enreg,my_natom,nattyp,nfft,& 131 & nfftmix,nfftmix_per_nfft,ngfft,ngfftmix,nkxc,npawmix,npwdiel,& 132 & nresid,ntypat,n1xccc,pawrhoij,pawtab,& 133 & ph1d,psps,rhog,rhor,rprimd,susmat,usepaw,vtrial,wvl,wvl_den,xred) 134 135 use defs_basis 136 use defs_datatypes 137 use defs_abitypes 138 use defs_wvltypes 139 use m_errors 140 use m_profiling_abi 141 use m_ab7_mixing 142 use m_abi2big 143 144 use m_geometry, only : metric 145 use m_pawtab, only : pawtab_type 146 use m_pawrhoij, only : pawrhoij_type 147 148 !This section has been created automatically by the script Abilint (TD). 149 !Do not modify the following lines by hand. 150 #undef ABI_FUNC 151 #define ABI_FUNC 'newrho' 152 use interfaces_18_timing 153 use interfaces_53_ffts 154 use interfaces_68_rsprc, except_this_one => newrho 155 !End of the abilint section 156 157 implicit none 158 159 !Arguments------------------------------- 160 !scalars 161 integer,intent(in) :: dielstrt,initialized,ispmix,istep,my_natom,mgfft 162 integer,intent(in) :: moved_atm_inside,n1xccc,nfft 163 integer,intent(in) :: nfftmix,nfftmix_per_nfft 164 integer,intent(in) :: nkxc,npawmix,npwdiel,ntypat,usepaw 165 integer,intent(inout) :: dbl_nnsclo 166 real(dp),intent(in) :: etotal,gsqcut 167 type(MPI_type),intent(in) :: mpi_enreg 168 type(ab7_mixing_object), intent(inout) :: mix 169 type(dataset_type),intent(in) :: dtset 170 type(pseudopotential_type),intent(in) :: psps 171 type(wvl_internal_type), intent(in) :: wvl 172 type(wvl_denspot_type), intent(inout) :: wvl_den 173 !arrays 174 integer,intent(in) :: atindx(dtset%natom) 175 integer,intent(in) :: ffttomix(nfft*(nfftmix_per_nfft)) 176 integer,intent(in) :: kg_diel(3,npwdiel) 177 integer,intent(in) :: mixtofft(nfftmix*nfftmix_per_nfft) 178 integer,intent(in) :: nattyp(ntypat),ngfft(18),ngfftmix(18) 179 real(dp),intent(in) :: dielar(7),fcart(3,dtset%natom),grhf(3,dtset%natom) 180 real(dp),intent(in) :: rprimd(3,3) 181 real(dp),intent(in) :: susmat(2,npwdiel,dtset%nspden,npwdiel,dtset%nspden) 182 real(dp),intent(in), target :: vtrial(nfft,dtset%nspden) 183 real(dp),intent(inout) :: dielinv(2,npwdiel,dtset%nspden,npwdiel,dtset%nspden) 184 real(dp),intent(inout), target :: dtn_pc(3,dtset%natom) 185 real(dp), intent(inout) :: gmet(3,3) 186 real(dp),intent(inout) :: kxc(nfft,nkxc),nresid(nfft,dtset%nspden) 187 real(dp),intent(inout) :: ph1d(2,3*(2*mgfft+1)*dtset%natom) 188 real(dp),intent(inout) :: rhor(nfft,dtset%nspden) 189 real(dp), intent(inout), target :: xred(3,dtset%natom) 190 real(dp),intent(inout) :: rhog(2,nfft) !vz_i 191 type(pawrhoij_type),intent(inout) :: pawrhoij(my_natom*psps%usepaw) 192 type(pawtab_type),intent(in) :: pawtab(ntypat*psps%usepaw) 193 194 !Local variables------------------------------- 195 !scalars 196 integer,parameter :: tim_fourdp9=9 197 integer :: cplex,dplex,errid,i_vresid1,i_vrespc1,iatom,ifft,indx,irhoij,ispden,jfft 198 integer :: jrhoij,klmn,kmix,mpicomm,nfftot,nselect 199 logical :: mpi_summarize,reset 200 real(dp) :: fact,ucvol,ucvol_local 201 character(len=500) :: message 202 !arrays 203 real(dp) :: gprimd(3,3),rmet(3,3),ro(2),tsec(2),vhartr_dum(1),vpsp_dum(1) 204 real(dp) :: vxc_dum(1,1) 205 real(dp),allocatable :: magng(:,:,:) 206 real(dp),allocatable :: nresid0(:,:),nrespc(:,:),nreswk(:,:,:) 207 real(dp),allocatable :: rhoijrespc(:),rhoijtmp(:,:) 208 real(dp), pointer :: rhomag(:,:), npaw(:) 209 210 ! ************************************************************************* 211 212 DBG_ENTER("COLL") 213 call timab(94,1,tsec) 214 215 nfftot=PRODUCT(ngfft(1:3)) 216 217 !Compatibility tests 218 if(nfftmix>nfft) then 219 MSG_BUG('nfftmix>nfft not allowed !') 220 end if 221 222 if(dtset%usewvl==1) then 223 if( (ispmix/=1 .or. nfftmix/=nfft)) then 224 MSG_BUG('newrho: nfftmix/=nfft, ispmix/=1 not allowed for wavelets') 225 end if 226 if(dtset%wvl_bigdft_comp==1) then 227 message = 'newrho: usewvl == 1 and wvl_bigdft_comp==1 not allowed!' 228 MSG_BUG(message) 229 end if 230 end if 231 232 if(ispmix/=2.and.nfftmix/=nfft) then 233 MSG_BUG('nfftmix/=nfft allowed only when ispmix=2 !') 234 end if 235 236 if (usepaw==1.and.my_natom>0) then 237 cplex=pawrhoij(1)%cplex;dplex=cplex-1 238 else 239 cplex = 0;dplex = 0 240 end if 241 242 !Compute different geometric tensor, as well as ucvol, from rprimd 243 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 244 245 if(dtset%usewvl==0) then 246 ucvol_local=ucvol 247 #if defined HAVE_BIGDFT 248 else 249 ucvol_local = product(wvl_den%denspot%dpbox%hgrids) * real(nfftot, dp) 250 #endif 251 end if 252 253 !Select components of density to be mixed 254 ABI_ALLOCATE(rhomag,(ispmix*nfftmix,dtset%nspden)) 255 ABI_ALLOCATE(nresid0,(ispmix*nfftmix,dtset%nspden)) 256 if (ispmix==1.and.nfft==nfftmix) then 257 rhomag(:,1:dtset%nspden)=rhor(:,1:dtset%nspden) 258 nresid0(:,1:dtset%nspden)=nresid(:,1:dtset%nspden) 259 else if (nfft==nfftmix) then 260 do ispden=1,dtset%nspden 261 call fourdp(1,nresid0(:,ispden),nresid(:,ispden),-1,mpi_enreg,nfft,ngfft,dtset%paral_kgb,tim_fourdp9) 262 end do 263 rhomag(:,1)=reshape(rhog,(/2*nfft/)) 264 if (dtset%nspden>1) then 265 do ispden=2,dtset%nspden 266 call fourdp(1,rhomag(:,ispden),rhor(:,ispden),-1,mpi_enreg,nfft,ngfft,dtset%paral_kgb,tim_fourdp9) 267 end do 268 end if 269 else 270 fact=dielar(4)-1._dp 271 ABI_ALLOCATE(nreswk,(2,nfft,dtset%nspden)) 272 do ispden=1,dtset%nspden 273 call fourdp(1,nreswk(:,:,ispden),nresid(:,ispden),-1,mpi_enreg,nfft,ngfft,dtset%paral_kgb,tim_fourdp9) 274 end do 275 do ifft=1,nfft 276 if (ffttomix(ifft)>0) then 277 jfft=2*ffttomix(ifft) 278 rhomag (jfft-1:jfft,1)=rhog(1:2,ifft) 279 nresid0(jfft-1:jfft,1)=nreswk(1:2,ifft,1) 280 else 281 rhog(:,ifft)=rhog(:,ifft)+fact*nreswk(:,ifft,1) 282 end if 283 end do 284 if (dtset%nspden>1) then 285 ABI_ALLOCATE(magng,(2,nfft,dtset%nspden-1)) 286 do ispden=2,dtset%nspden 287 call fourdp(1,magng(:,:,ispden-1),rhor(:,ispden),-1,mpi_enreg,nfft,ngfft,dtset%paral_kgb,tim_fourdp9) 288 do ifft=1,nfft 289 if (ffttomix(ifft)>0) then 290 jfft=2*ffttomix(ifft) 291 rhomag (jfft-1:jfft,ispden)=magng (1:2,ifft,ispden-1) 292 nresid0(jfft-1:jfft,ispden)=nreswk(1:2,ifft,ispden) 293 else 294 magng(:,ifft,ispden-1)=magng(:,ifft,ispden-1)+fact*nreswk(:,ifft,ispden) 295 if (dtset%nspden==2) magng(:,ifft,1)=two*magng(:,ifft,1)-rhog(:,ifft) 296 end if 297 end do 298 end do 299 end if 300 ABI_DEALLOCATE(nreswk) 301 end if 302 303 !Retrieve "input" density from "output" density and density residual 304 rhomag(:,1:dtset%nspden)=rhomag(:,1:dtset%nspden)-nresid0(:,1:dtset%nspden) 305 306 !If nspden==2, separate density and magnetization 307 if (dtset%nspden==2) then 308 rhomag (:,2)=two*rhomag (:,2)-rhomag (:,1) 309 nresid0(:,2)=two*nresid0(:,2)-nresid0(:,1) 310 end if 311 if (usepaw==1.and.my_natom>0) then 312 if (pawrhoij(1)%nspden==2) then 313 do iatom=1,my_natom 314 jrhoij=1 315 do irhoij=1,pawrhoij(iatom)%nrhoijsel 316 ro(1:1+dplex)=pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+dplex,1) 317 pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+dplex,1)=ro(1:1+dplex)+pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+dplex,2) 318 pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+dplex,2)=ro(1:1+dplex)-pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+dplex,2) 319 jrhoij=jrhoij+cplex 320 end do 321 do kmix=1,pawrhoij(iatom)%lmnmix_sz 322 klmn=cplex*pawrhoij(iatom)%kpawmix(kmix)-dplex 323 ro(1:1+dplex)=pawrhoij(iatom)%rhoijres(klmn:klmn+dplex,1) 324 pawrhoij(iatom)%rhoijres(klmn:klmn+dplex,1)=ro(1:1+dplex)+pawrhoij(iatom)%rhoijres(klmn:klmn+dplex,2) 325 pawrhoij(iatom)%rhoijres(klmn:klmn+dplex,2)=ro(1:1+dplex)-pawrhoij(iatom)%rhoijres(klmn:klmn+dplex,2) 326 end do 327 end do 328 end if 329 end if 330 331 !Choice of preconditioner governed by iprcel, densfor_pred and iprcfc 332 ABI_ALLOCATE(nrespc,(ispmix*nfftmix,dtset%nspden)) 333 ABI_ALLOCATE(npaw,(npawmix*usepaw)) 334 if (usepaw==1) then 335 ABI_ALLOCATE(rhoijrespc,(npawmix)) 336 else 337 ABI_ALLOCATE(rhoijrespc,(0)) 338 end if 339 if(dtset%usewvl==0) then 340 call prcref(atindx,dielar,dielinv,& 341 & dielstrt,dtn_pc,dtset,etotal,fcart,ffttomix,gmet,gsqcut,& 342 & istep,kg_diel,kxc,& 343 & mgfft,moved_atm_inside,mpi_enreg,my_natom,& 344 & nattyp,nfft,nfftmix,ngfft,ngfftmix,nkxc,npawmix,npwdiel,ntypat,n1xccc,& 345 & ispmix,1,pawrhoij,pawtab,ph1d,psps,rhog,rhoijrespc,rhor,rprimd,& 346 & susmat,vhartr_dum,vpsp_dum,nresid0,nrespc,vxc_dum,wvl,wvl_den,xred) 347 else 348 call wvl_prcref(dielar,dtset%iprcel,my_natom,nfftmix,npawmix,dtset%nspden,pawrhoij,& 349 & rhoijrespc,psps%usepaw,nresid0,nrespc) 350 end if 351 352 !------Compute new trial density and eventual new atomic positions 353 354 i_vresid1=mix%i_vresid(1) 355 i_vrespc1=mix%i_vrespc(1) 356 357 !Initialise working arrays for the mixing object. 358 if (moved_atm_inside == 1) then 359 call ab7_mixing_use_moving_atoms(mix, dtset%natom, xred, dtn_pc) 360 end if 361 call ab7_mixing_eval_allocate(mix, istep) 362 !Copy current step arrays. 363 if (moved_atm_inside == 1) then 364 call ab7_mixing_copy_current_step(mix, nresid0, errid, message, & 365 & arr_respc = nrespc, arr_atm = grhf) 366 else 367 call ab7_mixing_copy_current_step(mix, nresid0, errid, message, & 368 & arr_respc = nrespc) 369 end if 370 if (errid /= AB7_NO_ERROR) then 371 MSG_ERROR(message) 372 end if 373 ABI_DEALLOCATE(nresid0) 374 ABI_DEALLOCATE(nrespc) 375 376 !PAW: either use the array f_paw or the array f_paw_disk 377 if (usepaw==1) then 378 indx=-dplex 379 do iatom=1,my_natom 380 do ispden=1,pawrhoij(iatom)%nspden 381 ABI_ALLOCATE(rhoijtmp,(cplex*pawrhoij(iatom)%lmn2_size,1)) 382 rhoijtmp=zero 383 jrhoij=1 384 do irhoij=1,pawrhoij(iatom)%nrhoijsel 385 klmn=cplex*pawrhoij(iatom)%rhoijselect(irhoij)-dplex 386 rhoijtmp(klmn:klmn+dplex,1)=pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+dplex,ispden) 387 jrhoij=jrhoij+cplex 388 end do 389 do kmix=1,pawrhoij(iatom)%lmnmix_sz 390 indx=indx+cplex;klmn=cplex*pawrhoij(iatom)%kpawmix(kmix)-dplex 391 npaw(indx:indx+dplex)=rhoijtmp(klmn:klmn+dplex,1)-pawrhoij(iatom)%rhoijres(klmn:klmn+dplex,ispden) 392 mix%f_paw(indx:indx+dplex,i_vresid1)=pawrhoij(iatom)%rhoijres(klmn:klmn+dplex,ispden) 393 mix%f_paw(indx:indx+dplex,i_vrespc1)=rhoijrespc(indx:indx+dplex) 394 end do 395 ABI_DEALLOCATE(rhoijtmp) 396 end do 397 end do 398 end if 399 400 !------Prediction of the components of the density 401 402 !Init mpicomm 403 if(mpi_enreg%paral_kgb==1)then 404 mpicomm=mpi_enreg%comm_fft 405 mpi_summarize=.true. 406 else 407 mpicomm=0 408 mpi_summarize=.false. 409 end if 410 if(dtset%usewvl==1) then 411 mpicomm=mpi_enreg%comm_wvl 412 mpi_summarize=(mpi_enreg%nproc_wvl > 1) 413 end if 414 415 reset = .false. 416 if (initialized == 0) reset = .true. 417 call ab7_mixing_eval(mix, rhomag, istep, nfftot, ucvol_local, & 418 & mpicomm, mpi_summarize, errid, message, & 419 & reset = reset, isecur = dtset%isecur,& 420 & pawopt = dtset%pawoptmix, pawarr = npaw, & 421 & etotal = etotal, potden = vtrial, & 422 & comm_atom=mpi_enreg%comm_atom) 423 424 if (errid == AB7_ERROR_MIXING_INC_NNSLOOP) then 425 dbl_nnsclo = 1 426 else if (errid /= AB7_NO_ERROR) then 427 MSG_ERROR(message) 428 end if 429 430 !PAW: apply a simple mixing to rhoij (this is temporary) 431 if(dtset%iscf==15 .or. dtset%iscf==16)then 432 if (usepaw==1) then 433 indx=-dplex 434 do iatom=1,my_natom 435 ABI_ALLOCATE(rhoijtmp,(cplex*pawrhoij(iatom)%lmn2_size,pawrhoij(iatom)%nspden)) 436 rhoijtmp=zero 437 if (pawrhoij(iatom)%lmnmix_sz<pawrhoij(iatom)%lmn2_size) then 438 do ispden=1,pawrhoij(iatom)%nspden 439 do kmix=1,pawrhoij(iatom)%lmnmix_sz 440 indx=indx+cplex;klmn=cplex*pawrhoij(iatom)%kpawmix(kmix)-dplex 441 rhoijtmp(klmn:klmn+dplex,ispden)=rhoijrespc(indx:indx+dplex) & 442 & -pawrhoij(iatom)%rhoijres(klmn:klmn+dplex,ispden) 443 end do 444 end do 445 end if 446 if (pawrhoij(iatom)%nspden/=2) then 447 do ispden=1,pawrhoij(iatom)%nspden 448 jrhoij=1 449 do irhoij=1,pawrhoij(iatom)%nrhoijsel 450 klmn=cplex*pawrhoij(iatom)%rhoijselect(irhoij)-dplex 451 rhoijtmp(klmn:klmn+dplex,ispden)=rhoijtmp(klmn:klmn+dplex,ispden) & 452 & +pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+dplex,ispden) 453 jrhoij=jrhoij+cplex 454 end do 455 end do 456 else 457 jrhoij=1 458 do irhoij=1,pawrhoij(iatom)%nrhoijsel 459 klmn=cplex*pawrhoij(iatom)%rhoijselect(irhoij)-dplex 460 ro(1:1+dplex)=rhoijtmp(klmn:klmn+dplex,1) 461 rhoijtmp(klmn:klmn+dplex,1)=half*(ro(1:1+dplex)+rhoijtmp(klmn:klmn+dplex,2)) & 462 & +pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+dplex,1) 463 rhoijtmp(klmn:klmn+dplex,2)=half*(ro(1:1+dplex)-rhoijtmp(klmn:klmn+dplex,2)) & 464 & +pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+dplex,2) 465 jrhoij=jrhoij+cplex 466 end do 467 end if 468 nselect=0 469 do klmn=1,pawrhoij(iatom)%lmn2_size 470 if (any(abs(rhoijtmp(cplex*klmn-dplex:cplex*klmn,:))>tol10)) then 471 nselect=nselect+1 472 pawrhoij(iatom)%rhoijselect(nselect)=klmn 473 do ispden=1,pawrhoij(iatom)%nspden 474 pawrhoij(iatom)%rhoijp(cplex*nselect-dplex:cplex*nselect,ispden)=& 475 & rhoijtmp(cplex*klmn-dplex:cplex*klmn,ispden) 476 end do 477 end if 478 end do 479 pawrhoij(iatom)%nrhoijsel=nselect 480 ABI_DEALLOCATE(rhoijtmp) 481 end do 482 end if 483 end if 484 485 !if (usepaw==1) then 486 ABI_DEALLOCATE(rhoijrespc) 487 !end if 488 489 !PAW: restore rhoij from compact storage 490 if (usepaw==1.and.dtset%iscf/=15.and.dtset%iscf/=16) then 491 indx=-dplex 492 do iatom=1,my_natom 493 ABI_ALLOCATE(rhoijtmp,(cplex*pawrhoij(iatom)%lmn2_size,pawrhoij(iatom)%nspden)) 494 rhoijtmp=zero 495 if (pawrhoij(iatom)%lmnmix_sz<pawrhoij(iatom)%lmn2_size) then 496 do ispden=1,pawrhoij(iatom)%nspden 497 jrhoij=1 498 do irhoij=1,pawrhoij(iatom)%nrhoijsel 499 klmn=cplex*pawrhoij(iatom)%rhoijselect(irhoij)-dplex 500 rhoijtmp(klmn:klmn+dplex,ispden)=pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+dplex,ispden) 501 jrhoij=jrhoij+cplex 502 end do 503 end do 504 end if 505 do ispden=1,pawrhoij(iatom)%nspden 506 do kmix=1,pawrhoij(iatom)%lmnmix_sz 507 indx=indx+cplex;klmn=cplex*pawrhoij(iatom)%kpawmix(kmix)-dplex 508 rhoijtmp(klmn:klmn+dplex,ispden)=npaw(indx:indx+dplex) 509 end do 510 end do 511 if (pawrhoij(iatom)%nspden==2) then 512 jrhoij=1 513 do irhoij=1,pawrhoij(iatom)%nrhoijsel 514 klmn=cplex*pawrhoij(iatom)%rhoijselect(irhoij)-dplex 515 ro(1:1+dplex)=rhoijtmp(klmn:klmn+dplex,1) 516 rhoijtmp(klmn:klmn+dplex,1)=half*(ro(1:1+dplex)+rhoijtmp(klmn:klmn+dplex,2)) 517 rhoijtmp(klmn:klmn+dplex,2)=half*(ro(1:1+dplex)-rhoijtmp(klmn:klmn+dplex,2)) 518 jrhoij=jrhoij+cplex 519 end do 520 end if 521 nselect=0 522 if (cplex==1) then 523 do klmn=1,pawrhoij(iatom)%lmn2_size 524 if (any(abs(rhoijtmp(klmn,:))>tol10)) then 525 nselect=nselect+1 526 pawrhoij(iatom)%rhoijselect(nselect)=klmn 527 do ispden=1,pawrhoij(iatom)%nspden 528 pawrhoij(iatom)%rhoijp(nselect,ispden)=rhoijtmp(klmn,ispden) 529 end do 530 end if 531 end do 532 else 533 do klmn=1,pawrhoij(iatom)%lmn2_size 534 if (any(abs(rhoijtmp(2*klmn-1:2*klmn,:))>tol10)) then 535 nselect=nselect+1 536 pawrhoij(iatom)%rhoijselect(nselect)=klmn 537 do ispden=1,pawrhoij(iatom)%nspden 538 pawrhoij(iatom)%rhoijp(2*nselect-1:2*nselect,ispden)=rhoijtmp(2*klmn-1:2*klmn,ispden) 539 end do 540 end if 541 end do 542 end if 543 pawrhoij(iatom)%nrhoijsel=nselect 544 ABI_DEALLOCATE(rhoijtmp) 545 end do 546 end if 547 ABI_DEALLOCATE(npaw) 548 549 !Eventually write the data on disk and deallocate f_fftgr_disk 550 call ab7_mixing_eval_deallocate(mix) 551 552 !Fourier transform the density 553 if (ispmix==1.and.nfft==nfftmix) then 554 rhor(:,1:dtset%nspden)=rhomag(:,1:dtset%nspden) 555 if(dtset%usewvl==0) then 556 call fourdp(1,rhog,rhor(:,1),-1,mpi_enreg,nfft,ngfft,dtset%paral_kgb,tim_fourdp9) 557 end if 558 else if (nfft==nfftmix) then 559 do ispden=1,dtset%nspden 560 call fourdp(1,rhomag(:,ispden),rhor(:,ispden),+1,mpi_enreg,nfft,ngfft,dtset%paral_kgb,tim_fourdp9) 561 end do 562 rhog(:,:)=reshape(rhomag(:,1),(/2,nfft/)) 563 else 564 do ifft=1,nfftmix 565 jfft=mixtofft(ifft) 566 rhog(1:2,jfft)=rhomag(2*ifft-1:2*ifft,1) 567 end do 568 call fourdp(1,rhog,rhor(:,1),+1,mpi_enreg,nfft,ngfft,dtset%paral_kgb,tim_fourdp9) 569 if (dtset%nspden>1) then 570 do ispden=2,dtset%nspden 571 do ifft=1,nfftmix 572 jfft=mixtofft(ifft) 573 magng(1:2,jfft,ispden-1)=rhomag(2*ifft-1:2*ifft,ispden) 574 end do 575 call fourdp(1,magng(:,:,ispden-1),rhor(:,ispden),+1,mpi_enreg,nfft,ngfft,dtset%paral_kgb,tim_fourdp9) 576 end do 577 ABI_DEALLOCATE(magng) 578 end if 579 end if 580 ABI_DEALLOCATE(rhomag) 581 582 !Set back rho in (up+dn,up) form if nspden=2 583 if (dtset%nspden==2) rhor(:,2)=half*(rhor(:,1)+rhor(:,2)) 584 585 !In WVL: copy density to BigDFT object: 586 if(dtset%usewvl==1) then 587 call wvl_rho_abi2big(1,rhor,wvl_den) 588 end if 589 590 call timab(94,2,tsec) 591 592 DBG_EXIT("COLL") 593 594 end subroutine newrho