TABLE OF CONTENTS
ABINIT/prcref [ Functions ]
NAME
prcref
FUNCTION
Compute preconditioned residual potential (or density) and forces. iprcel, densfor_pred and iprcfc govern the choice of the preconditioner. Three tasks are done : 1) Preconditioning of the forces (residual has already been included) using the approximate force constant matrix. Get proposed change of atomic positions. 2) Precondition the residual, get first part of proposed trial potential change. 3) PAW only: precondition the rhoij residuals (simple preconditionning) 4) Take into account the proposed change of atomic positions to modify the proposed trial potential change. NOTE This routine is almost similar to prcref_PMA.F90 which is employed in case of potential mixing. Yet it has undergone strong changes simultaneously from two different sources at the same time which resulted in a splitting.
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. dielstrt=number of the step at which the dielectric preconditioning begins. dtset <type(dataset_type)>=all input variables in this dataset | intxc=control xc quadrature | densfor_pred= not yet used here | iprcel= governs the preconditioning of the potential residual | 0 => simple model dielectric matrix, described by the | parameters dielng, diemac, diemix and diemixmag contained in dielar. | between 21 and 39 => until istep=dielstart, same as iprcel=0, then uses | the RPA dielectric matrix (routine dielmt) | between 41 and 49 => uses the RPA dielectric matrix (routine dielmt). | between 51 and 59 => uses the RPA dielectric matrix (routine dieltcel). | between 61 and 69 => uses the electronic dielectric matr (routine dieltcel). | between 71 and 79 => uses the real-space preconditioner based on Kerker prc (prcrskerkerN) | between 81 and 99 => reserved for futur version of the real-space preconditioner | between 141 and 169 -> same as between 41 and 69 but with a different periodicity: modulo(iprcel modulo (10)) | iprcfc= governs the preconditioning of the forces | 0 => hessian is the identity matrix | 1 => hessian is 0.5 times the identity matrix | 2 => hessian is 0.25 times the identity matrix | ixc=exchange-correlation choice parameter. | natom=number of atoms | nspden=number of spin-density components | occopt=option for occupancies | prtvol=control print volume and debugging | typat(natom)=integer type for each atom in cell etotal=total ennergy fcart(3,natom)=cartesian forces (hartree/bohr) ffttomix(nfft*(1-nfftprc/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. gsqcut=cutoff on (k+G)^2 (bohr^-2) istep= number of the step in the SCF cycle kg_diel(3,npwdiel)=reduced planewave coordinates for the dielectric matrix. mgfft=maximum size of 1D FFTs 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=number of fft grid points nfftprc=size of FFT grid on which the potential residual will be preconditionned ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft ngfftprc(18)=contain all needed information about 3D FFT for the grid corresponding to nfftprc nkxc=second dimension of the array kxc, see rhotoxc.F90 for a description npawmix=-PAW only- number of spherical part elements to be mixed npwdiel=number of planewaves for dielectric matrix ntypat=number of types of atoms in cell. n1xccc=dimension of xccc1d ; 0 if no XC core correction is used optreal=1 if residual potential is is REAL space, 2 if it is in RECIPROCAL SPACE optres=0: the array vresid contains a potential residual 1: the array vresid contains a density residual 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 rhog(2,nfft)=array for electron density in reciprocal space 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 vresid(optreal*nfftprc,nspden)=residual potential vxc(nfft,nspden)=exchange-correlation potential (hartree) vhartr(nfft)=array for holding Hartree potential vlspl(mqgrid,2,ntypat)=q^2 v(q) spline for each type of atom. vpsp(nfft)=array for holding local psp xred(3,natom)=reduced dimensionless atomic coordinates
OUTPUT
dtn_pc(3,natom)=preconditioned change of atomic position,
in reduced coordinates
vrespc(optreal*nfftprc,nspden)=preconditioned residual of the potential
==== if psps%usepaw==1
rhoijrespc(npawmix)= preconditionned rhoij residuals at output
SIDE EFFECT
dielinv(2,npwdiel,nspden,npwdiel,nspden)=
inverse of the dielectric matrix in rec. space
kxc(nfft,nkxc)=exchange-correlation kernel,
needed if the electronic dielectric matrix is computed
===== if densfor_pred==3 .and. moved_atm_inside==1 =====
ph1d(2,3*(2*mgfft+1)*natom)=1-dim structure factor phases
PARENTS
newrho
CHILDREN
atm2fft,dielmt,dieltcel,fourdp,fresid,getph,hartre indirect_parallel_fourier,kgindex,mean_fftr,metric,mkcore,mklocl moddiel,prcrskerker1,prcrskerker2,rhotoxc,testsusmat,xcart2xred xcdata_init,xmpi_sum,zerosym
SOURCE
129 #if defined HAVE_CONFIG_H 130 #include "config.h" 131 #endif 132 133 #include "abi_common.h" 134 135 subroutine prcref(atindx,dielar,dielinv,& 136 & dielstrt,dtn_pc,dtset,etotal,fcart,ffttomix,gmet,gsqcut,& 137 & istep,kg_diel,kxc,& 138 & mgfft,moved_atm_inside,mpi_enreg,my_natom,& 139 & nattyp,nfft,nfftprc,ngfft,ngfftprc,nkxc,npawmix,npwdiel,ntypat,n1xccc,& 140 & optreal,optres,pawrhoij,pawtab,ph1d,psps,rhog,rhoijrespc,rhor,rprimd,& 141 & susmat,vhartr,vpsp,vresid,vrespc,vxc,wvl,wvl_den,xred) 142 143 use defs_basis 144 use defs_datatypes 145 use defs_abitypes 146 use defs_wvltypes 147 use m_errors 148 use m_profiling_abi 149 use m_xmpi 150 use m_cgtools 151 use m_xcdata 152 153 use m_geometry, only : xcart2xred, metric 154 use m_pawtab, only : pawtab_type 155 use m_pawrhoij, only : pawrhoij_type 156 use m_fft, only : zerosym 157 use m_kg, only : getph 158 use m_dtset, only : testsusmat 159 160 !This section has been created automatically by the script Abilint (TD). 161 !Do not modify the following lines by hand. 162 #undef ABI_FUNC 163 #define ABI_FUNC 'prcref' 164 use interfaces_52_fft_mpi_noabirule 165 use interfaces_53_ffts 166 use interfaces_56_xc 167 use interfaces_64_psp 168 use interfaces_67_common 169 use interfaces_68_rsprc, except_this_one => prcref 170 !End of the abilint section 171 172 implicit none 173 174 !Arguments------------------------------- 175 !scalars 176 integer,intent(in) :: dielstrt,istep,my_natom,mgfft,moved_atm_inside,n1xccc 177 integer,intent(in) :: nfft,nfftprc,nkxc,npawmix,npwdiel,ntypat,optreal,optres 178 real(dp),intent(in) :: etotal,gsqcut 179 type(MPI_type),intent(in) :: mpi_enreg 180 type(dataset_type),intent(in) :: dtset 181 type(pseudopotential_type),intent(in) :: psps 182 type(wvl_internal_type), intent(in) :: wvl 183 type(wvl_denspot_type), intent(inout) :: wvl_den 184 !arrays 185 integer,intent(in) :: atindx(dtset%natom),ffttomix(nfft*(1-nfftprc/nfft)) 186 integer,intent(in) :: kg_diel(3,npwdiel),nattyp(ntypat),ngfft(18),ngfftprc(18) 187 real(dp),intent(in) :: dielar(7),fcart(3,dtset%natom),rhog(2,nfft) 188 real(dp),intent(in) :: rhor(nfft,dtset%nspden),rprimd(3,3) 189 real(dp),intent(in) :: susmat(2,npwdiel,dtset%nspden,npwdiel,dtset%nspden) 190 real(dp),intent(in) :: vhartr(nfft),vresid(nfftprc*optreal,dtset%nspden) 191 real(dp),intent(in) :: vxc(nfft,dtset%nspden) 192 real(dp),intent(inout) :: dielinv(2,npwdiel,dtset%nspden,npwdiel,dtset%nspden) 193 real(dp),intent(inout) :: gmet(3,3),kxc(nfft,nkxc) 194 real(dp),intent(inout) :: ph1d(2,3*(2*mgfft+1)*dtset%natom),vpsp(nfft) 195 real(dp),intent(inout) :: xred(3,dtset%natom) 196 real(dp),intent(out) :: dtn_pc(3,dtset%natom),rhoijrespc(npawmix) 197 real(dp),intent(out) :: vrespc(nfftprc*optreal,dtset%nspden) 198 type(pawrhoij_type),intent(inout) :: pawrhoij(my_natom*psps%usepaw) 199 type(pawtab_type),intent(in) :: pawtab(ntypat*psps%usepaw) 200 201 !Local variables------------------------------- 202 !scalars 203 integer :: coredens_method,cplex,dielop,iatom,ier,ifft,ii,index,ipw1 204 integer :: ipw2,ispden,klmn,kmix,n1,n2,n3,n3xccc,nfftot,nk3xc,optatm 205 integer :: optdyfr,opteltfr,optgr,option,optn,optn2,optstr,optv,vloc_method 206 real(dp) :: ai,ar,diemix,diemixmag,eei,enxc 207 real(dp) :: mixfac 208 real(dp) :: mixfac_eff,mixfacmag,ucvol,vxcavg 209 logical :: computediel 210 logical :: non_magnetic_xc 211 character(len=500) :: message 212 type(xcdata_type) :: xcdata 213 !arrays 214 integer :: qprtrb(3) 215 integer,allocatable :: indpw_prc(:) 216 real(dp) :: dummy6(6),dummy7(6),gprimd(3,3),qphon(3),rmet(3,3),strsxc(6) 217 real(dp) :: vmean(dtset%nspden),vprtrb(2) 218 real(dp),allocatable :: dummy(:),dummy1(:),dummy2(:),dummy3(:),dummy4(:),dummy5(:),dummy8(:),dummy9(:) 219 real(dp),allocatable :: dyfrlo_indx(:,:,:),dyfrx2(:,:,:) 220 real(dp),allocatable :: fcart_pc(:,:),gresid(:,:),grtn_indx(:,:) 221 real(dp),allocatable :: grxc(:,:),grxc_indx(:,:),rhog_wk(:,:),rhor_new(:,:) 222 real(dp),allocatable :: rhor_wk(:,:),rhor_wk0(:,:),vhartr_wk(:),vpsp_wk(:) 223 real(dp),allocatable :: vres_diel(:,:),vxc_wk(:,:),work(:),work1(:,:),work2(:) 224 real(dp),allocatable :: work3(:,:),xccc3d(:),xred_wk(:,:) 225 logical,allocatable :: mask(:) 226 227 ! ************************************************************************* 228 229 !Compute different geometric tensor, as well as ucvol, from rprimd 230 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 231 232 ! Initialise non_magnetic_xc for rhohxc 233 non_magnetic_xc=(dtset%usepawu==4).or.(dtset%usepawu==14) 234 235 !1) Eventually take care of the forces 236 237 if(moved_atm_inside==1)then 238 ABI_ALLOCATE(fcart_pc,(3,dtset%natom)) 239 240 if(dtset%iprcfc==0)then 241 fcart_pc(:,:)=fcart(:,:) 242 else 243 fcart_pc(:,:)= (two**dtset%iprcfc) * fcart(:,:) 244 end if 245 246 ! Compute preconditioned delta xred from preconditioned fcart and rprimd 247 call xcart2xred(dtset%natom,rprimd,fcart_pc,dtn_pc) 248 249 ABI_DEALLOCATE(fcart_pc) 250 end if 251 252 !####################################################################### 253 254 !2) Take care of the potential residual 255 256 !Compute the residuals corresponding to the solution 257 !of an approximate realspace dielectric function according 258 !to X. Gonze PRB vol54 nb7 p4383 (1996) 259 if(dtset%iprcel>=71.and.dtset%iprcel<=79) then 260 if (nfft==nfftprc) then 261 if (dtset%iprcel<=78) then 262 call prcrskerker1(dtset,mpi_enreg,nfft,dtset%nspden,ngfft,dielar,etotal, & 263 & gprimd,vresid,vrespc,rhor(:,1)) 264 else 265 call prcrskerker2(dtset,nfft,dtset%nspden,ngfft,dielar,gprimd,rprimd, & 266 & vresid,vrespc,dtset%natom,xred,mpi_enreg,ucvol) 267 end if 268 else 269 ! If preconditionning has to be done on a coarse grid, 270 ! has to transfer several arrays 271 ABI_ALLOCATE(work1,(nfftprc,dtset%nspden)) 272 ABI_ALLOCATE(work3,(nfftprc,dtset%nspden)) 273 ABI_ALLOCATE(work,(2*nfftprc)) 274 do ispden=1,dtset%nspden 275 work(:)=vresid(:,ispden) 276 call fourdp(1,work,work1(:,ispden),+1,mpi_enreg,nfftprc,ngfftprc,dtset%paral_kgb,0) 277 end do 278 ABI_DEALLOCATE(work) 279 if (dtset%iprcel<=78) then 280 ABI_ALLOCATE(rhog_wk,(2,nfftprc)) 281 rhog_wk(:,:)=zero 282 if (mpi_enreg%nproc_fft>1.and. mpi_enreg%paral_kgb==1) then 283 nfftot=PRODUCT(ngfft(1:3)) 284 call indirect_parallel_Fourier(ffttomix,rhog_wk,mpi_enreg,ngfftprc,& 285 & ngfft,nfftprc,nfft,dtset%paral_kgb,rhog,nfftot) 286 else 287 do ii=1,nfft 288 if (ffttomix(ii)>0) rhog_wk(:,ffttomix(ii))=rhog(:,ii) 289 end do 290 end if 291 call zerosym(rhog_wk,2,ngfftprc(1),ngfftprc(2),ngfftprc(3),& 292 & comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft) 293 ABI_ALLOCATE(work,(nfftprc)) 294 call fourdp(1,rhog_wk,work,+1,mpi_enreg,nfftprc,ngfftprc,dtset%paral_kgb,0) 295 call prcrskerker1(dtset,mpi_enreg,nfftprc,dtset%nspden,ngfftprc,dielar,etotal, & 296 & gprimd,work1,work3,work) 297 ABI_DEALLOCATE(work) 298 else 299 call prcrskerker2(dtset,nfftprc,dtset%nspden,ngfftprc,dielar,gprimd,rprimd, & 300 & work1,work3,dtset%natom,xred,mpi_enreg,ucvol) 301 end if 302 do ispden=1,dtset%nspden 303 call fourdp(1,vrespc(:,ispden),work3(:,ispden),-1,mpi_enreg,nfftprc,ngfftprc,dtset%paral_kgb,0) 304 end do 305 ABI_DEALLOCATE(work1) 306 ABI_DEALLOCATE(work3) 307 end if 308 309 else 310 311 if(dtset%iprcel==0 .or. (dtset%iprcel<40.and.istep<dielstrt) )then 312 cplex=optreal 313 qphon(:)=zero 314 ! Simple scalar multiplication, or model dielectric function 315 call moddiel(cplex,dielar,mpi_enreg,nfftprc,ngfftprc,dtset%nspden,optreal,optres,dtset%paral_kgb,qphon,rprimd,vresid,vrespc) 316 317 ! Use the inverse dielectric matrix in a small G sphere 318 else if( (istep>=dielstrt .and. dtset%iprcel>=21) .or. modulo(dtset%iprcel,100)>=41 )then 319 320 ! With dielop=1, the matrices will be computed when istep=dielstrt 321 ! With dielop=2, the matrices will be computed when istep=dielstrt and 1 322 dielop=1 323 if(modulo(dtset%iprcel,100)>=41)dielop=2 324 call testsusmat(computediel,dielop,dielstrt,dtset,istep) !test if the matrix is to be computed 325 if(computediel) then 326 ! Compute the inverse dielectric matrix from the susceptibility matrix 327 ! There are two routines for the RPA matrix, while for the electronic 328 ! dielectric matrix, only dieltcel will do the work 329 if(modulo(dtset%iprcel,100)<=49)then 330 call dielmt(dielinv,gmet,kg_diel,& 331 & npwdiel,dtset%nspden,dtset%occopt,dtset%prtvol,susmat) 332 else 333 option=1 334 if(modulo(dtset%iprcel,100)>=61)option=2 335 call dieltcel(dielinv,gmet,kg_diel,kxc,& 336 & nfft,ngfft,nkxc,npwdiel,dtset%nspden,dtset%occopt,option,dtset%paral_kgb,dtset%prtvol,susmat) 337 end if 338 end if 339 340 ABI_ALLOCATE(work1,(2,nfftprc)) 341 ABI_ALLOCATE(work2,(optreal*nfftprc)) 342 343 ! Presently, one uses the inverse of the RPA dielectric matrix, 344 ! for which spin must be averaged. 345 346 ! Do fft from real space (work2) to G space (work1) 347 if (optreal==1) then 348 work2(:)=vresid(:,1) 349 ! Must average over spins in the case of a potential residual 350 if(dtset%nspden/=1.and.optres==0)work2(:)=(work2(:)+vresid(:,2))*half 351 call fourdp(1,work1,work2,-1,mpi_enreg,nfftprc,ngfftprc,dtset%paral_kgb,0) 352 else 353 work1(:,:)=reshape(vresid(:,1),(/2,nfftprc/)) 354 if(dtset%nspden/=1.and.optres==0)work1(:,:)=(work1(:,:)+reshape(vresid(:,2),(/2,nfftprc/)))*half 355 end if 356 357 ! Multiply by restricted inverse of dielectric matrix. 358 ! Must first copy relevant elements of work1 to a npwdiel-dimensioned array, 359 ! then zero work1, operate with the dielinv matrix, and store in work1. 360 361 ABI_ALLOCATE(vres_diel,(2,npwdiel)) 362 ABI_ALLOCATE(indpw_prc,(npwdiel)) 363 ABI_ALLOCATE(mask,(npwdiel)) 364 mask(:)=.true. 365 call kgindex(indpw_prc,kg_diel,mask,mpi_enreg,ngfftprc,npwdiel) 366 367 do ipw1=1,npwdiel 368 if(mask(ipw1)) then 369 vres_diel(1,ipw1)=work1(1,indpw_prc(ipw1)) 370 vres_diel(2,ipw1)=work1(2,indpw_prc(ipw1)) 371 end if 372 end do 373 work1(:,:)=zero 374 do ipw1=1,npwdiel 375 ar=zero ; ai=zero 376 377 ! Use inverse of dielectric matrix (potential mixing) 378 if (optres==0) then 379 do ipw2=1,npwdiel 380 if(mask(ipw2))then 381 ar=ar+dielinv(1,ipw1,1,ipw2,1)*vres_diel(1,ipw2) & 382 & -dielinv(2,ipw1,1,ipw2,1)*vres_diel(2,ipw2) 383 ai=ai+dielinv(2,ipw1,1,ipw2,1)*vres_diel(1,ipw2) & 384 & +dielinv(1,ipw1,1,ipw2,1)*vres_diel(2,ipw2) 385 end if 386 end do 387 else 388 ! Use symetric of inverse of dielectric matrix (density mixing) 389 do ipw2=1,npwdiel 390 if(mask(ipw2))then 391 ar=ar+dielinv(1,ipw2,1,ipw1,1)*vres_diel(1,ipw2) & 392 & +dielinv(2,ipw2,1,ipw1,1)*vres_diel(2,ipw2) 393 ai=ai-dielinv(2,ipw2,1,ipw1,1)*vres_diel(1,ipw2) & 394 & +dielinv(1,ipw2,1,ipw1,1)*vres_diel(2,ipw2) 395 end if 396 end do 397 end if 398 ! Must be careful not to count the diagonal 1 twice : it is added later, 399 ! so must be subtracted now. 400 call xmpi_sum(ar,mpi_enreg%comm_fft,ier) 401 call xmpi_sum(ai,mpi_enreg%comm_fft,ier) 402 if(mask(ipw1)) then 403 work1(1,indpw_prc(ipw1))=ar-vres_diel(1,ipw1) 404 work1(2,indpw_prc(ipw1))=ai-vres_diel(2,ipw1) 405 end if !mask(ipw1) 406 end do ! ipw1 407 ABI_DEALLOCATE(vres_diel) 408 ABI_DEALLOCATE(indpw_prc) 409 ABI_DEALLOCATE(mask) 410 ! Fourier transform 411 if (optreal==1) then 412 call fourdp(1,work1,work2,1,mpi_enreg,nfftprc,ngfftprc,dtset%paral_kgb,0) 413 else 414 work2(:)=reshape(work1(:,:),(/nfftprc*2/)) 415 end if 416 417 ! Add to get the preconditioned vresid, must be careful about spins. 418 if(dtset%iprcel>=30)then 419 diemix=dielar(4);diemixmag=abs(dielar(7)) 420 vrespc(:,1)=diemix*(vresid(:,1)+work2(:)) 421 if(dtset%nspden/=1.and.optres==0)vrespc(:,2)=diemixmag*(vresid(:,2)+work2(:)) 422 if(dtset%nspden==4.and.optres==0)vrespc(:,3:4)=diemixmag*vresid(:,3:4) 423 if(dtset%nspden/=1.and.optres==1)vrespc(:,2:dtset%nspden)=diemixmag*vresid(:,2:dtset%nspden) 424 else 425 vrespc(:,1)=vresid(:,1)+work2(:) 426 if(dtset%nspden/=1.and.optres==0)vrespc(:,2)=vresid(:,2)+work2(:) 427 if(dtset%nspden==4.and.optres==0)vrespc(:,3:4)=vresid(:,3:4) 428 if(dtset%nspden/=1.and.optres==1)vrespc(:,2:dtset%nspden)=vresid(:,2:dtset%nspden) 429 end if 430 431 ABI_DEALLOCATE(work1) 432 ABI_DEALLOCATE(work2) 433 434 ! Other choice ? 435 else 436 write(message, '(a,i3,a,a,a,a)' )& 437 & 'From the calling routine, iprcel=',dtset%iprcel,ch10,& 438 & 'The only allowed values are 0 or larger than 20.',ch10,& 439 & 'Action: correct your input file.' 440 MSG_ERROR(message) 441 end if 442 end if 443 !####################################################################### 444 445 !3) PAW only : precondition the rhoij quantities (augmentation 446 !occupancies) residuals. Use a simple preconditionning 447 !with the same mixing factor as the model dielectric function. 448 449 if (psps%usepaw==1.and.my_natom>0) then 450 if (istep>=dielstrt.and.dtset%iprcel>=21.and.dtset%iprcel<30) then 451 mixfac=one;mixfacmag=one 452 else 453 mixfac=dielar(4);mixfacmag=abs(dielar(7)) 454 end if 455 if (pawrhoij(1)%cplex==1) then 456 index=0 457 do iatom=1,my_natom 458 do ispden=1,pawrhoij(iatom)%nspden 459 mixfac_eff=mixfac;if (ispden>1) mixfac_eff=mixfacmag 460 do kmix=1,pawrhoij(iatom)%lmnmix_sz 461 index=index+1;klmn=pawrhoij(iatom)%kpawmix(kmix) 462 rhoijrespc(index)=mixfac_eff*pawrhoij(iatom)%rhoijres(klmn,ispden) 463 end do 464 end do 465 end do 466 else 467 index=-1 468 do iatom=1,my_natom 469 do ispden=1,pawrhoij(iatom)%nspden 470 mixfac_eff=mixfac;if (ispden>1) mixfac_eff=mixfacmag 471 do kmix=1,pawrhoij(iatom)%lmnmix_sz 472 index=index+2;klmn=2*pawrhoij(iatom)%kpawmix(kmix)-1 473 rhoijrespc(index:index+1)=mixfac_eff*pawrhoij(iatom)%rhoijres(klmn:klmn+1,ispden) 474 end do 475 end do 476 end do 477 end if 478 end if 479 480 !####################################################################### 481 482 !4) Take care of the change of atomic positions 483 !Note : this part is very demanding on memory... 484 !however, since this algorithm is still in development, 485 !it was NOT included in the estimation provided by memory.f 486 if(abs(dtset%densfor_pred)==3 .and. moved_atm_inside==1)then 487 488 ! Not yet compatible with resid given in reciprocal space 489 if (optreal/=1) then 490 write(message, '(5a)' )& 491 & 'From the calling routine, densfor_pred=3',ch10,& 492 & 'You cannot use residuals in reciprocal space.',ch10,& 493 & 'Action: correct your input file.' 494 MSG_ERROR(message) 495 end if 496 ! Not compatible with non-collinear magnetism 497 if(dtset%nspden==4)then 498 MSG_ERROR('densfor_pred=3 does not work for nspden=4 !') 499 end if 500 501 n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3) 502 nfftot=PRODUCT(ngfft(1:3)) 503 504 if (optres==0) then ! Array vresid contains a potential residual 505 ! ----------------------------------------------------------------- 506 507 ! First subtract the current local, hartree and exchange correlation potentials 508 do ispden=1,min(dtset%nspden,2) 509 vrespc(:,ispden)=vrespc(:,ispden)-vpsp(:)-vhartr(:)-vxc(:,ispden) 510 end do 511 if (dtset%nspden==4) then 512 do ispden=3,4 513 vrespc(:,ispden)=vrespc(:,ispden)-vxc(:,ispden) 514 end do 515 end if 516 517 ! Compute the modified density, in rhor_wk 518 option=2 519 ABI_ALLOCATE(gresid,(3,dtset%natom)) 520 ABI_ALLOCATE(grxc,(3,dtset%natom)) 521 ABI_ALLOCATE(rhor_wk,(nfft,dtset%nspden)) 522 ABI_ALLOCATE(rhor_wk0,(nfft,dtset%nspden)) 523 ABI_ALLOCATE(xred_wk,(3,dtset%natom)) 524 xred_wk(:,:)=xred(:,:)+dtn_pc(:,:) 525 call fresid(dtset,gresid,mpi_enreg,nfft,ngfft,& 526 & ntypat,option,pawtab,rhor,rprimd,& 527 & ucvol,rhor_wk,xred_wk,xred,psps%znuclpsp) 528 529 ! Compute up+down rhog_wk(G) by fft 530 ABI_ALLOCATE(work,(nfft)) 531 ABI_ALLOCATE(rhog_wk,(2,nfft)) 532 work(:)=rhor_wk(:,1) 533 call fourdp(1,rhog_wk,work,-1,mpi_enreg,nfft,ngfft,dtset%paral_kgb,0) 534 ABI_DEALLOCATE(work) 535 536 ! Compute structure factor phases for new atomic pos: 537 call getph(atindx,dtset%natom,n1,n2,n3,ph1d,xred_wk) 538 539 ! Compute local ionic pseudopotential vpsp: 540 ! and core electron density xccc3d, if needed. 541 n3xccc=0;if (n1xccc/=0) n3xccc=nfft 542 ABI_ALLOCATE(xccc3d,(n3xccc)) 543 ABI_ALLOCATE(vpsp_wk,(nfft)) 544 vprtrb(1:2)=zero 545 546 ! Determine by which method the local ionic potential and/or 547 ! the pseudo core charge density contributions have to be computed 548 ! Local ionic potential: 549 ! Method 1: PAW 550 ! Method 2: Norm-conserving PP, icoulomb>0, wavelets 551 vloc_method=1;if (psps%usepaw==0) vloc_method=2 552 if (dtset%icoulomb>0) vloc_method=2 553 if (psps%usewvl==1) vloc_method=2 554 ! Pseudo core charge density: 555 ! Method 1: PAW, nc_xccc_gspace 556 ! Method 2: Norm-conserving PP, wavelets 557 coredens_method=1;if (psps%usepaw==0) coredens_method=2 558 if (psps%nc_xccc_gspace==1) coredens_method=1 559 if (psps%nc_xccc_gspace==0) coredens_method=2 560 if (psps%usewvl==1) coredens_method=2 561 562 ! Local ionic potential and/or pseudo core charge by method 1 563 if (vloc_method==1.or.coredens_method==1) then 564 optv=0;if (vloc_method==1) optv=1 565 optn=0;if (coredens_method==1) optn=n3xccc/nfft 566 optatm=1;optdyfr=0;optgr=0;optstr=0;optn2=1;opteltfr=0 567 ! Note: atindx1 should be passed to atm2fft (instead of atindx) but it is unused... 568 call atm2fft(atindx,xccc3d,vpsp,dummy,dummy2,dummy9,dummy1,gmet,gprimd,dummy3,dummy4,gsqcut,& 569 & mgfft,psps%mqgrid_vl,dtset%natom,nattyp,nfft,ngfft,ntypat,optatm,optdyfr,opteltfr,optgr,optn,optn2,& 570 & optstr,optv,psps,pawtab,ph1d,psps%qgrid_vl,qprtrb,dummy5,dummy6,dummy7,& 571 & ucvol,psps%usepaw,dummy8,dummy8,dummy8,vprtrb,psps%vlspl,& 572 & comm_fft=mpi_enreg%comm_fft,me_g0=mpi_enreg%me_g0,& 573 & paral_kgb=mpi_enreg%paral_kgb,distribfft=mpi_enreg%distribfft) 574 end if 575 576 ! Local ionic potential by method 2 577 if (vloc_method==2) then 578 option=1 579 ABI_ALLOCATE(dyfrlo_indx,(3,3,dtset%natom)) 580 ABI_ALLOCATE(grtn_indx,(3,dtset%natom)) 581 call mklocl(dtset,dyfrlo_indx,eei,gmet,gprimd,grtn_indx,gsqcut,dummy6,& 582 & mgfft,mpi_enreg,dtset%natom,nattyp,nfft,ngfft,dtset%nspden,& 583 & ntypat,option,pawtab,ph1d,psps,qprtrb,rhog_wk,rhor_wk,rprimd,& 584 & ucvol,vprtrb,vpsp_wk,wvl,wvl_den,xred) 585 ABI_DEALLOCATE(dyfrlo_indx) 586 ABI_DEALLOCATE(grtn_indx) 587 end if 588 589 ! Pseudo core electron density by method 2 590 if (coredens_method==2.and.n1xccc/=0) then 591 option=1 592 ABI_ALLOCATE(dyfrx2,(3,3,dtset%natom)) 593 ABI_ALLOCATE(grxc_indx,(3,dtset%natom)) 594 call mkcore(dummy6,dyfrx2,grxc_indx,mpi_enreg,dtset%natom,nfft,dtset%nspden,ntypat,& 595 & n1,n1xccc,n2,n3,option,rprimd,dtset%typat,ucvol,vxc,psps%xcccrc,& 596 & psps%xccc1d,xccc3d,xred_wk) 597 ABI_DEALLOCATE(dyfrx2) 598 ABI_DEALLOCATE(grxc_indx) 599 end if 600 601 ! Compute Hartree+xc potentials 602 ABI_ALLOCATE(vxc_wk,(nfft,dtset%nspden)) 603 ABI_ALLOCATE(vhartr_wk,(nfft)) 604 option=1 605 606 call hartre(1,gsqcut,psps%usepaw,mpi_enreg,nfft,ngfft,dtset%paral_kgb,rhog_wk,rprimd,vhartr_wk) 607 608 ! Prepare the call to rhotoxc 609 call xcdata_init(xcdata,dtset=dtset) 610 nk3xc=1 611 call rhotoxc(enxc,kxc,mpi_enreg,nfft,ngfft,& 612 & work,0,work,0,nkxc,nk3xc,non_magnetic_xc,n3xccc,option,dtset%paral_kgb,rhor_wk,rprimd,strsxc,1,& 613 & vxc_wk,vxcavg,xccc3d,xcdata,vhartr=vhartr_wk) 614 ABI_DEALLOCATE(xccc3d) 615 616 ! Sum all contributions 617 do ispden=1,min(dtset%nspden,2) 618 do ifft=1,nfft 619 vrespc(ifft,ispden)=vrespc(ifft,ispden)+vpsp_wk(ifft)+vhartr_wk(ifft)+vxc_wk(ifft,ispden) 620 end do 621 end do 622 if (dtset%nspden==4) then 623 do ispden=3,4 624 do ifft=1,nfft 625 vrespc(ifft,ispden)=vrespc(ifft,ispden)+vxc_wk(ifft,ispden) 626 end do 627 end do 628 end if 629 call mean_fftr(vrespc,vmean,nfft,nfftot,dtset%nspden,mpi_comm_sphgrid=mpi_enreg%comm_fft) 630 if(dtset%nspden==2) then 631 vmean(1)=half*(vmean(1)+vmean(2)) 632 vmean(2)=vmean(1) 633 end if 634 do ispden=1,dtset%nspden 635 vrespc(:,ispden)=vrespc(:,ispden)-vmean(ispden) 636 end do 637 ABI_DEALLOCATE(gresid) 638 ABI_DEALLOCATE(grxc) 639 ABI_DEALLOCATE(rhog_wk) 640 ABI_DEALLOCATE(rhor_wk) 641 ABI_DEALLOCATE(rhor_wk0) 642 ABI_DEALLOCATE(xred_wk) 643 ABI_DEALLOCATE(vhartr_wk) 644 ABI_DEALLOCATE(vpsp_wk) 645 ABI_DEALLOCATE(vxc_wk) 646 647 else ! Array vresid contains a density residual 648 ! ----------------------------------------------------------------- 649 650 ! Only have to compute the modified preconditionned density residual 651 option=2 652 ABI_ALLOCATE(gresid,(3,dtset%natom)) 653 ABI_ALLOCATE(grxc,(3,dtset%natom)) 654 ABI_ALLOCATE(rhor_new,(nfft,dtset%nspden)) 655 ABI_ALLOCATE(rhor_wk,(nfft,dtset%nspden)) 656 ABI_ALLOCATE(rhor_wk0,(nfft,dtset%nspden)) 657 ABI_ALLOCATE(xred_wk,(3,dtset%natom)) 658 xred_wk(:,:)=xred(:,:)+dtn_pc(:,:) 659 rhor_new(:,1)=rhor(:,1)+vrespc(:,1) 660 if (dtset%nspden==2) then 661 rhor_new(:,1)=rhor_new(:,1)+vrespc(:,2) 662 rhor_new(:,2)=rhor(:,2)+vrespc(:,1) 663 end if 664 call fresid(dtset,gresid,mpi_enreg,nfft,ngfft,& 665 & ntypat,option,pawtab,rhor,rprimd,& 666 & ucvol,rhor_wk0,xred_wk,xred,psps%znuclpsp) 667 call fresid(dtset,gresid,mpi_enreg,nfft,ngfft,& 668 & ntypat,option,pawtab,rhor_new,rprimd,& 669 & ucvol,rhor_wk,xred_wk,xred,psps%znuclpsp) 670 vrespc(:,1)=rhor_wk(:,dtset%nspden)-rhor_wk0(:,dtset%nspden) 671 if (dtset%nspden==2) vrespc(:,2)=rhor_wk(:,1)-rhor_wk0(:,1)-vrespc(:,1) 672 ABI_DEALLOCATE(gresid) 673 ABI_DEALLOCATE(grxc) 674 ABI_DEALLOCATE(rhor_new) 675 ABI_DEALLOCATE(rhor_wk) 676 ABI_DEALLOCATE(rhor_wk0) 677 ABI_DEALLOCATE(xred_wk) 678 end if 679 680 end if 681 682 end subroutine prcref