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