TABLE OF CONTENTS
ABINIT/m_rhotov [ Modules ]
NAME
m_rhotov
FUNCTION
COPYRIGHT
Copyright (C) 1998-2024 ABINIT group (XG, GMR, MT, EB) 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
15 #if defined HAVE_CONFIG_H 16 #include "config.h" 17 #endif 18 19 #include "abi_common.h" 20 21 module m_rhotov 22 23 use defs_basis 24 use defs_wvltypes 25 use m_errors 26 use m_abicore 27 use m_ab7_mixing 28 use m_abi2big 29 use m_xmpi 30 use m_cgtools 31 use m_xcdata 32 use m_dtset 33 34 use defs_abitypes, only : MPI_type 35 use m_time, only : timab 36 use m_geometry, only : xred2xcart 37 use m_energies, only : energies_type 38 use m_electronpositron, only : electronpositron_type, electronpositron_calctype, rhohxcpositron 39 use libxc_functionals, only : libxc_functionals_is_hybrid 40 use m_spacepar, only : hartre 41 use m_dens, only : constrained_dft_t,mag_penalty,constrained_residual 42 use m_rhotoxc, only : rhotoxc 43 use m_xchybrid, only : xchybrid_ncpp_cc 44 use m_psolver, only : psolver_rhohxc 45 use m_wvl_psi, only : wvl_psitohpsi 46 use m_pawang, only : pawang_type 47 use m_pawrad, only : pawrad_type 48 use m_pawrhoij, only : pawrhoij_type 49 use m_pawtab, only : pawtab_type 50 use m_xc_tb09, only : xc_tb09_update_c 51 52 implicit none 53 54 private
ABINIT/rhotov [ Functions ]
NAME
rhotov
FUNCTION
This routine is called to compute, from a given total density the trial (local) potential and the residual potential.
INPUTS
[add_tfw]=flag controling the addition of Weiszacker gradient correction to Thomas-Fermi kin energy constrained_dft <type(constrained_dft_t>=data for constrained dft calculations dtset <type(dataset_type)>=all input variables in this dataset | spinmagntarget=input variable that governs fixed moment calculation | natom=number of atoms in cell. | nspden=number of spin-density components | ntypat=number of types of atoms in unit cell. | occopt=option for occupancies | typat(natom)=type (integer) for each atom gprimd(3,3)=dimensional reciprocal space primitive translations gsqcut=cutoff on (k+G)^2 (bohr^-2) mpi_enreg=information about MPI parallelization nfft=(effective) number of FFT grid points (for this processor) ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft nhat(nfft,nspden*usepaw)= -PAW only- compensation density nhatgr(nfft,nspden,3*nhatgrdim)= -PAW only- cartesian gradients of compensation density nhatgrdim= -PAW only- 0 if nhatgr array is not used ; 1 otherwise nkxc=second dimension of the array kxc, see rhotoxc.F90 for a description n3xccc=dimension of the xccc3d array (0 or nfft). optene=option for the computation of additional energies optres=0: the trial potential residual is computed ; the input potential value is kept 1: the new value of the trial potential is computed in place of the input value optxc=option to be used for the call to rhotoxc pawang <type(pawang_type)> =paw angular mesh and related data pawrad(ntypat*usepaw) <type(pawrad_type)>=paw radial mesh and related data pawrhoij <type(pawrhoij_type)>= paw rhoij occupancies and related data (for the current atom) pawtab(ntypat*dtset%usepaw) <type(pawtab_type)>=paw tabulated starting data rhog(2,nfft)=array for Fourier transform of electron density rhor(nfft,nspden)=array for electron density in electrons/bohr**3. | definition for spin components: | case of nspden = 2 | rhor(:,1) => rho_up + rho_dwn | rhor(:,2) => rho_up | case of nspden = 4 | rhor(:,1) => rho_upup + rho_dwndwn | rhor(:,2:4) => {m_x,m_y,m_z} rprimd(3,3)=dimensional primitive translations in real space (bohr) [taur(nfftf,nspden*dtset%usekden)]=array for kinetic energy density ucvol = unit cell volume (Bohr**3) usepaw= 0 for non paw calculation; =1 for paw calculation usexcnhat= -PAW only- flag controling use of compensation density in Vxc vpsp(nfft)=array for holding local psp [vxc_hybcomp(nfft,nspden)= compensation xc potential (Hartree) in case of hybrids] Optional output i.e. difference between the hybrid Vxc at fixed density and the auxiliary Vxc at fixed density xccc3d(n3xccc)=3D core electron density for XC core correction (bohr^-3) ==== if optres==0 vtrial(nfft,nspden)= old value of trial potential
OUTPUT
energies <type(energies_type)>=all part of total energy. | e_hartree=Hartree part of total energy (hartree units) | e_xc=exchange-correlation energy (hartree) | In case of hybrid compensation algorithm: | e_hybcomp_v=self-consistent potential compensation term for the exchange-correlation energy (hartree) ==== if optene==0.or.2 | e_localpsp=local psp energy (hartree) ==== if optene==1.or.2 | e_xcdc=exchange-correlation double-counting energy (hartree) grcondft(3,natom)=d(E_constrained_DFT)/d(xred) (hartree) intgres(nspden,ngrcondft)=integrated residuals from constrained DFT. They are also Lagrange parameters, or gradients with respect to constraints. kxc(nfft,nkxc)=exchange-correlation kernel, needed only if optxc==2. strscondft(6)=constrained DFT contribution to stress tensor (hartree/bohr^3) strsxc(6)=xc contribution to stress tensor (hartree/bohr^3) vxc(nfft,nspden)=Vxc(r) (already computed above; gets recomputed below too) vxcavg=mean of the vxc potential ==== if optres==0 vresidnew(nfft,nspden)=potential residual vnew_mean(nspden)=mean of the potential formed from vpsp, vhartr and vxc, might be spin-dependent vres_mean(nspden)=mean of the potential residual, might be spin-dependent vres2=square of the norm of the residual [vxctau(nfftf,dtset%nspden,4*dtset%usekden)]=derivative of XC energy density with respect to kinetic energy density (metaGGA cases) (optional output) [vtauresid(nfft,nspden*dtset%usekden)]=array for vxctau residue (see vtau below))
SIDE EFFECTS
Input/Output: electronpositron <type(electronpositron_type)>=quantities for the electron-positron annihilation (optional argument) vhartr(nfft)=array for holding Hartree potential ==== if optres==1 vtrial(nfft,nspden)= new value of trial potential
NOTES
In case of PAW calculations:
All computations are done on the fine FFT grid.
All variables (nfft,ngfft,mgfft) refer to this fine FFT grid.
All arrays (densities/potentials...) are computed on this fine FFT grid.
! Developpers have to be careful when introducing others arrays:
they have to be stored on the fine FFT grid.
In case of norm-conserving calculations the FFT grid is the usual FFT grid.
SOURCE
164 subroutine rhotov(constrained_dft,dtset,energies,gprimd,grcondft,gsqcut,intgres,istep,kxc,mpi_enreg,nfft,ngfft,& 165 & nhat,nhatgr,nhatgrdim,nkxc,vresidnew,n3xccc,optene,optres,optxc,& 166 & pawang,pawrad,pawrhoij,pawtab,rhog,rhor,rprimd,strscondft,strsxc,ucvol,usepaw,usexcnhat,& 167 & vhartr,vnew_mean,vpsp,vres_mean,vres2,vtrial,vxcavg,vxc,wvl,xccc3d,xred,& 168 & electronpositron,taur,vxc_hybcomp,vxctau,vtauresid,add_tfw,xcctau3d) ! optional arguments 169 170 !Arguments ------------------------------------ 171 !scalars 172 integer,intent(in) :: n3xccc,nfft,nhatgrdim,nkxc,optene,optres,optxc,usepaw,istep 173 integer,intent(in) :: usexcnhat 174 logical,intent(in),optional :: add_tfw 175 real(dp),intent(in) :: gsqcut,ucvol 176 real(dp),intent(out) :: vres2,vxcavg 177 type(MPI_type),intent(inout) :: mpi_enreg 178 type(constrained_dft_t),intent(inout) :: constrained_dft 179 type(dataset_type),intent(in) :: dtset 180 type(electronpositron_type),pointer,optional :: electronpositron 181 type(energies_type),intent(inout) :: energies 182 type(pawang_type),intent(in) :: pawang 183 type(wvl_data), intent(inout) :: wvl 184 !arrays 185 integer,intent(in) :: ngfft(18) 186 real(dp),intent(in) :: gprimd(3,3),nhat(nfft,dtset%nspden*usepaw) 187 real(dp),intent(in) :: nhatgr(nfft,dtset%nspden,3*nhatgrdim),rhog(2,nfft) 188 real(dp),intent(in) :: rprimd(3,3) 189 real(dp),intent(inout) :: rhor(nfft,dtset%nspden),vhartr(nfft),vpsp(nfft) 190 real(dp),intent(inout) :: vtrial(nfft,dtset%nspden),vxc(nfft,dtset%nspden) 191 real(dp),intent(inout) :: xccc3d(n3xccc),xred(3,dtset%natom) 192 real(dp),intent(out) :: grcondft(:,:) ! (3,ngrcondft) ngrcondft=natom when condft is activated 193 real(dp),intent(out) :: intgres(:,:) ! (nspden,ngrcondft) ngrcondft=natom when condft is activated 194 real(dp),intent(out) :: kxc(nfft,nkxc),strsxc(6),vnew_mean(dtset%nspden) 195 real(dp),intent(out) :: strscondft(6) 196 real(dp),intent(out) :: vres_mean(dtset%nspden),vresidnew(nfft,dtset%nspden) 197 real(dp),intent(in),optional :: taur(nfft,dtset%nspden*dtset%usekden) 198 real(dp),intent(inout),optional :: vtauresid(nfft,dtset%nspden*dtset%usekden) 199 real(dp),intent(out),optional,target :: vxctau(nfft,dtset%nspden,4*dtset%usekden) 200 real(dp),intent(out),optional :: vxc_hybcomp(:,:) ! (nfft,nspden) 201 real(dp),intent(out),optional :: xcctau3d(n3xccc) 202 type(pawrhoij_type),intent(in) :: pawrhoij(:) 203 type(pawrad_type),intent(in) :: pawrad(dtset%ntypat*dtset%usepaw) 204 type(pawtab_type),intent(in) :: pawtab(dtset%ntypat*dtset%usepaw) 205 206 !Local variables------------------------------- 207 !scalars 208 integer :: nk3xc,ifft,ipositron,ispden,nfftot,offset 209 integer :: mpi_comm_sphgrid,ixc_current 210 !integer :: ii,jj,kk,ipt,nx,ny,nz !SPr: debug 211 !real(dp):: rx,ry,rz !SPr: debug 212 real(dp) :: doti,e_xcdc_vxctau 213 logical :: add_tfw_,calc_xcdc,non_magnetic_xc,with_vxctau 214 logical :: is_hybrid_ncpp,wvlbigdft=.false. 215 type(xcdata_type) :: xcdata 216 !arrays 217 real(dp) :: evxc,tsec(2),vmean(dtset%nspden),vzeeman(dtset%nspden) 218 real(dp),target :: vxctau_dum(0,0,0) 219 real(dp),allocatable :: rhowk(:,:),v_constr_dft_r(:,:),vnew(:,:),xcart(:,:) 220 real(dp),pointer :: vxctau_(:,:,:) 221 !real(dp),allocatable :: vzeemanHarm(:,:) !SPr: debug Zeeman field q/=0 real space 222 223 ! ********************************************************************* 224 225 DBG_ENTER("COLL") 226 227 call timab(940,1,tsec) 228 229 !Check that usekden is not 0 if want to use vxctau 230 with_vxctau = (present(vxctau).and.present(taur).and.(dtset%usekden/=0)) 231 vxctau_ => vxctau_dum ; if (with_vxctau) vxctau_ => vxctau 232 if (with_vxctau.and.optres==0.and.(.not.present(vtauresid))) then 233 ABI_BUG('need vtauresid!') 234 end if 235 236 !Check if we're in hybrid norm conserving pseudopotential with a core correction 237 is_hybrid_ncpp=(usepaw==0 .and. n3xccc/=0 .and. & 238 & (dtset%ixc==41.or.dtset%ixc==42.or.libxc_functionals_is_hybrid())) 239 240 !If usewvl: wvlbigdft indicates that the BigDFT workflow will be followed 241 wvlbigdft=(dtset%usewvl==1.and.dtset%wvl_bigdft_comp==1) 242 243 !mpi communicator for spherical grid 244 mpi_comm_sphgrid=mpi_enreg%comm_fft 245 if(dtset%usewvl==1) mpi_comm_sphgrid=mpi_enreg%comm_wvl 246 247 !Get size of FFT grid 248 nfftot=PRODUCT(ngfft(1:3)) 249 250 ipositron=0;if (present(electronpositron)) ipositron=electronpositron_calctype(electronpositron) 251 add_tfw_=.false.;if (present(add_tfw)) add_tfw_=add_tfw 252 253 !------Compute Hartree and xc potentials---------------------------------- 254 255 !allocate vnew here. 256 !In wvl: vnew is used at call to wvl_psitohpsi 257 if (optres==0) then 258 ABI_MALLOC(vnew,(nfft,dtset%nspden)) 259 vmean(:)=zero ; vnew_mean(:)=zero 260 end if 261 262 if (ipositron/=1) then 263 ! if metaGGA, save current value of vxctau potential 264 if (with_vxctau) vtauresid(:,:)=vxctau(:,:,1) 265 ! Compute xc potential (separate up and down if spin-polarized) 266 if (dtset%icoulomb == 0 .and. dtset%usewvl == 0) then 267 268 ! >>>> Hartree potential 269 call hartre(1,gsqcut,dtset%icutcoul,usepaw,mpi_enreg,nfft,ngfft,& 270 &dtset%nkpt,dtset%rcut,rhog,rprimd,dtset%vcutgeo,vhartr) 271 272 ! >>>> Exchange-correlation potential 273 !Use the proper exchange_correlation energy : either the origin one, or the auxiliary one 274 ixc_current=dtset%ixc 275 if(mod(dtset%fockoptmix,100)==11)ixc_current=dtset%auxc_ixc 276 call xcdata_init(xcdata,dtset=dtset,ixc=ixc_current) 277 non_magnetic_xc=(dtset%usepaw==1.and.mod(abs(dtset%usepawu),10)==4) 278 nk3xc=1 279 280 ! If we use the XC Tran-Blaha 2009 (modified BJ) functional, update the c value 281 if (dtset%xc_tb09_c>99._dp) then 282 call xc_tb09_update_c(dtset%intxc,dtset%ixc,mpi_enreg,dtset%natom, & 283 & nfft,ngfft,nhat,usepaw,nhatgr,nhatgrdim,dtset%nspden,dtset%ntypat,n3xccc, & 284 & pawang,pawrad,pawrhoij,pawtab,dtset%pawxcdev,rhor,rprimd,usepaw, & 285 & xccc3d,dtset%xc_denpos,comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab, & 286 & computation_type='all') 287 end if 288 289 ! Use the periodic solver to compute Hxc. 290 call timab(941,1,tsec) 291 if (ipositron==0) then 292 if(.not.is_hybrid_ncpp .or. mod(dtset%fockoptmix,100)==11)then 293 call rhotoxc(energies%e_xc,kxc,mpi_enreg,nfft,ngfft,& 294 & nhat,usepaw,nhatgr,nhatgrdim,nkxc,nk3xc,non_magnetic_xc,n3xccc,optxc,& 295 & rhor,rprimd,strsxc,usexcnhat,vxc,vxcavg,xccc3d,xcdata,& 296 & taur=taur,vhartr=vhartr,vxctau=vxctau_,add_tfw=add_tfw_,xcctau3d=xcctau3d) 297 if(mod(dtset%fockoptmix,100)==11)then 298 energies%e_xc=energies%e_xc*dtset%auxc_scal 299 vxc(:,:)=vxc(:,:)*dtset%auxc_scal 300 end if 301 else 302 call xchybrid_ncpp_cc(dtset,energies%e_xc,mpi_enreg,nfft,ngfft,n3xccc,rhor,rprimd,& 303 & strsxc,vxcavg,xccc3d,vxc=vxc) 304 end if 305 else 306 call rhotoxc(energies%e_xc,kxc,mpi_enreg,nfft,ngfft,& 307 & nhat,usepaw,nhatgr,nhatgrdim,nkxc,nk3xc,non_magnetic_xc,n3xccc,optxc,& 308 & rhor,rprimd,strsxc,usexcnhat,vxc,vxcavg,xccc3d,xcdata,& 309 & taur=taur,vhartr=vhartr,vxctau=vxctau_,add_tfw=add_tfw_,& 310 & electronpositron=electronpositron,xcctau3d=xcctau3d) 311 end if 312 call timab(941,2,tsec) 313 elseif (.not. wvlbigdft) then 314 ! Use the free boundary solver. 315 call timab(943,1,tsec) 316 call psolver_rhohxc(energies%e_hartree, energies%e_xc, evxc, & 317 & dtset%icoulomb, dtset%ixc, & 318 & mpi_enreg, nfft, & 319 & ngfft, nhat,usepaw,& 320 & dtset%nscforder, dtset%nspden, n3xccc, rhor,rprimd,& 321 & usexcnhat,dtset%usepaw,dtset%usewvl,vhartr, vxc, vxcavg,& 322 & wvl%descr,wvl%den,wvl%e,& 323 & xccc3d,dtset%xclevel,dtset%xc_denpos) 324 call timab(943,2,tsec) 325 end if 326 ! For icoulomb==0 and usewvl Ehartree is calculated in psolver_rhohxc(). 327 ! For PAW we recalculate this since nhat was not taken into account 328 ! in psolver_rhohxc: E_H= int v_H (n+nhat) dr 329 if(.not. wvlbigdft .and. (dtset%icoulomb==0 .or. dtset%usepaw==1 ) ) then 330 call timab(942,1,tsec) 331 call dotprod_vn(1,rhor,energies%e_hartree,doti,nfft,nfftot,1,1,vhartr,ucvol,mpi_comm_sphgrid=mpi_comm_sphgrid) 332 energies%e_hartree=half*energies%e_hartree 333 call timab(942,2,tsec) 334 end if 335 else 336 call timab(944,1,tsec) 337 energies%e_hartree=zero;energies%e_xc=zero 338 call rhohxcpositron(electronpositron,gprimd,kxc,mpi_enreg,nfft,ngfft,nhat,nkxc,dtset%nspden,n3xccc,& 339 & dtset%paral_kgb,rhor,strsxc,ucvol,usexcnhat,usepaw,vhartr,vxc,vxcavg,xccc3d,dtset%xc_denpos) 340 call timab(944,2,tsec) 341 end if 342 343 call timab(945,1,tsec) 344 if (ipositron/=0) then 345 call dotprod_vn(1,rhor,electronpositron%e_hartree,doti,& 346 & nfft,nfftot,1,1,electronpositron%vha_ep,ucvol,mpi_comm_sphgrid=mpi_comm_sphgrid) 347 vhartr=vhartr+electronpositron%vha_ep 348 end if 349 350 !------Compute parts of total energy depending on potentials-------- 351 352 if ( (optene==0.or.optene==2 ).and. .not. wvlbigdft) then 353 ! Compute local psp energy energies%e_localpsp 354 call dotprod_vn(1,rhor,energies%e_localpsp,doti,nfft,nfftot,1,1,vpsp,ucvol,& 355 & mpi_comm_sphgrid=mpi_comm_sphgrid) 356 end if 357 358 if(mod(dtset%fockoptmix,100)==11)then 359 if (.not. wvlbigdft) then 360 ! Compute second compensation energy for hybrid functionals 361 call dotprod_vn(1,rhor,energies%e_hybcomp_v,doti,nfft,nfftot,1,1,vxc_hybcomp,ucvol,& 362 & mpi_comm_sphgrid=mpi_comm_sphgrid) 363 end if 364 end if 365 366 calc_xcdc=.false. 367 if (optene==1.or.optene==2) calc_xcdc=.true. 368 if (dtset%usewvl==1.and.dtset%nnsclo>0) calc_xcdc=.true. 369 if (wvlbigdft) calc_xcdc=.false. 370 if (dtset%usefock==1) calc_xcdc=.true. 371 372 if (calc_xcdc) then 373 374 ! Compute double-counting XC energy energies%e_xcdc 375 if (ipositron/=1) then 376 if (usepaw==0.or.usexcnhat/=0) then 377 call dotprod_vn(1,rhor,energies%e_xcdc,doti,nfft,nfftot,dtset%nspden,1,vxc,ucvol,& 378 & mpi_comm_sphgrid=mpi_comm_sphgrid) 379 else 380 ABI_MALLOC(rhowk,(nfft,dtset%nspden)) 381 rhowk=rhor-nhat 382 call dotprod_vn(1,rhowk,energies%e_xcdc,doti,nfft,nfftot,dtset%nspden,1,vxc,ucvol,& 383 & mpi_comm_sphgrid=mpi_comm_sphgrid) 384 ABI_FREE(rhowk) 385 end if 386 if (with_vxctau)then 387 call dotprod_vn(1,taur,e_xcdc_vxctau,doti,nfft,nfftot,dtset%nspden,1,vxctau(:,:,1),& 388 & ucvol,mpi_comm_sphgrid=mpi_comm_sphgrid) 389 energies%e_xcdc=energies%e_xcdc+e_xcdc_vxctau 390 end if 391 if (ipositron==2) energies%e_xcdc=energies%e_xcdc-electronpositron%e_xcdc 392 else 393 energies%e_xcdc=zero 394 end if 395 396 end if 397 398 !------Produce residual vector and square norm of it------------- 399 !(only if requested ; if optres==0) 400 401 !Set up array for Zeeman field 402 !EB vzeeman(:) = factor*( Hz, Hx+iHy; Hx-iHy, -Hz) 403 !EB factor = -g/2 * mu_B * mu_0 = -1/2*B in a.u. 404 !EB <-- vzeeman might have to be allocated correctly --> to be checked 405 ! vzeeman = 1/2 ( -B_z, -B_x + iB_y ; -B_x - iB_y , B_z) 406 vzeeman(:) = zero 407 ! ABI_MALLOC(vzeemanHarm,(nfft,dtset%nspden)) ! SPr: debug stuff 408 ! vzeemanHarm(:,:) = zero ! 409 if (any(abs(dtset%zeemanfield(:))>tol8)) then 410 if(dtset%nspden==2)then 411 ! EB The collinear case has to be checked : 412 ! EB Is it vzeeman(1) or (2) that has to be added here? to be checked in setvtr and energy as well 413 ! SPr: the density components are: rhor(1) => n_upup + n_dwndwn 414 ! rhor(2) => n_upup 415 ! the convention for the potential components is different: 416 ! v(1) => v_upup 417 ! v(2) => v_dndn 418 ! verified by comparing collinear and non-collinear calculations 419 420 vzeeman(1) =-half*dtset%zeemanfield(3) ! v_upup 421 vzeeman(2) = half*dtset%zeemanfield(3) ! v_dndn 422 423 !vzeeman(1) = zero ! v_upup 424 !vzeeman(2) = zero ! v_dndn 425 426 !nx=ngfft(1); ny=ngfft(2); nz=ngfft(3) 427 !do kk=0,nz-1 428 ! do jj=0,ny-1 429 ! do ii=0,nx-1 430 ! ipt=1+ii+nx*(jj+ny*kk) 431 ! !rx=(dble(ii)/nx)*rprimd(1,1)+(dble(jj)/ny)*rprimd(1,2)+(dble(kk)/nz)*rprimd(1,3) 432 ! !ry=(dble(ii)/nx)*rprimd(2,1)+(dble(jj)/ny)*rprimd(2,2)+(dble(kk)/nz)*rprimd(2,3) 433 ! !rz=(dble(ii)/nx)*rprimd(3,1)+(dble(jj)/ny)*rprimd(3,2)+(dble(kk)/nz)*rprimd(3,3) 434 ! vzeemanHarm(ipt,1)= -half*dtset%zeemanfield(3)*cos(2*PI*(dble(ii)/dble(nx))) 435 ! vzeemanHarm(ipt,2)= half*dtset%zeemanfield(3)*cos(2*PI*(dble(ii)/dble(nx))) 436 ! end do 437 ! end do 438 !end do 439 440 else if(dtset%nspden==4)then 441 442 vzeeman(1)=-half*dtset%zeemanfield(3) ! v_upup 443 vzeeman(2)= half*dtset%zeemanfield(3) ! v_dndn 444 vzeeman(3)=-half*dtset%zeemanfield(1) ! Re(v_updn) 445 vzeeman(4)= half*dtset%zeemanfield(2) ! Im(v_updn) 446 447 !vzeeman(1)=0.0 448 !vzeeman(2)=0.0 449 !vzeeman(3)=0.0 450 !vzeeman(4)=0.0 451 452 !nx=ngfft(1); ny=ngfft(2); nz=ngfft(3) 453 !do kk=0,nz-1 454 ! do jj=0,ny-1 455 ! do ii=0,nx-1 456 ! ipt=1+ii+nx*(jj+ny*kk) 457 ! !rx=(dble(ii)/nx)*rprimd(1,1)+(dble(jj)/ny)*rprimd(1,2)+(dble(kk)/nz)*rprimd(1,3) 458 ! !ry=(dble(ii)/nx)*rprimd(2,1)+(dble(jj)/ny)*rprimd(2,2)+(dble(kk)/nz)*rprimd(2,3) 459 ! !rz=(dble(ii)/nx)*rprimd(3,1)+(dble(jj)/ny)*rprimd(3,2)+(dble(kk)/nz)*rprimd(3,3) 460 ! vzeemanHarm(ipt,1)= -half*dtset%zeemanfield(3)*cos(2*PI*(dble(ii)/dble(nx))) 461 ! vzeemanHarm(ipt,2)= half*dtset%zeemanfield(3)*cos(2*PI*(dble(ii)/dble(nx))) 462 ! vzeemanHarm(ipt,3)= -half*dtset%zeemanfield(1)*cos(2*PI*(dble(ii)/dble(nx))) 463 ! vzeemanHarm(ipt,4)= half*dtset%zeemanfield(2)*cos(2*PI*(dble(ii)/dble(nx))) 464 ! end do 465 ! end do 466 !end do 467 468 end if 469 end if 470 471 !Compute the constrained potential for the magnetic moments 472 ABI_MALLOC(v_constr_dft_r, (nfft,dtset%nspden)) 473 v_constr_dft_r = zero 474 if (dtset%magconon==1.or.dtset%magconon==2) then 475 call mag_penalty(constrained_dft,mpi_enreg,rhor,v_constr_dft_r,xred) 476 end if 477 478 if (optres==0) then 479 480 481 ! ------ Compute potential residual ------------- 482 483 if (.not. wvlbigdft) then 484 !$OMP PARALLEL DO COLLAPSE(2) 485 do ispden=1,min(dtset%nspden,2) 486 do ifft=1,nfft 487 vnew(ifft,ispden)=vhartr(ifft)+vpsp(ifft)+vxc(ifft,ispden)+vzeeman(ispden)+v_constr_dft_r(ifft,ispden) 488 !vnew(ifft,ispden)=vnew(ifft,ispden)+vzeemanHarm(ifft,ispden) 489 if(mod(dtset%fockoptmix,100)==11)vnew(ifft,ispden)=vnew(ifft,ispden)+vxc_hybcomp(ifft,ispden) 490 vresidnew(ifft,ispden)=vnew(ifft,ispden)-vtrial(ifft,ispden) 491 end do 492 end do 493 if(dtset%nspden==4)then 494 !$OMP PARALLEL DO COLLAPSE(2) 495 do ispden=3,4 496 do ifft=1,nfft 497 vnew(ifft,ispden)=vxc(ifft,ispden)+vzeeman(ispden)+v_constr_dft_r(ifft,ispden) 498 !vnew(ifft,ispden)=vnew(ifft,ispden)+vzeemanHarm(ifft,ispden) 499 if(mod(dtset%fockoptmix,100)==11)vnew(ifft,ispden)=vnew(ifft,ispden)+vxc_hybcomp(ifft,ispden) 500 vresidnew(ifft,ispden)=vnew(ifft,ispden)-vtrial(ifft,ispden) 501 end do 502 end do 503 end if 504 505 !If constrained_dft, must take into account the constraints, and recompute the residual and the new potential 506 if( any(dtset%constraint_kind(:)/=0))then 507 call constrained_residual(constrained_dft,energies%e_constrained_dft,& 508 & grcondft,intgres,mpi_enreg,rhor,strscondft,vresidnew,xred) 509 vnew(:,1:dtset%nspden)=vtrial(:,1:dtset%nspden)+vresidnew(:,1:dtset%nspden) 510 endif 511 512 offset = 0 513 514 if (dtset%iscf==0) vtrial=vnew 515 516 ! Pass vtrial to BigDFT object 517 if(dtset%usewvl==1) then 518 call wvl_vtrial_abi2big(1,vnew,wvl%den) 519 ! call wvl_vtrial_abi2big(1,vtrial,wvl%den) 520 end if 521 522 else 523 ! Compute with covering comms the different part of the potential. 524 ! only for wvlbigdft 525 ABI_MALLOC(xcart,(3, dtset%natom)) 526 call xred2xcart(dtset%natom, rprimd, xcart, xred) 527 call wvl_psitohpsi(dtset%diemix,energies%e_exactX, energies%e_xc, energies%e_hartree, & 528 & energies%e_kinetic, energies%e_localpsp, energies%e_nlpsp_vfock, energies%e_sicdc, & 529 & istep + 1, 1, dtset%iscf, mpi_enreg%me_wvl, dtset%natom, dtset%nfft,& 530 & mpi_enreg%nproc_wvl, dtset%nspden, & 531 & vres2, .true., energies%e_xcdc, wvl,& 532 & wvlbigdft, xcart, strsxc,vtrial=vnew,vxc=vxc) 533 ABI_FREE(xcart) 534 535 vresidnew = vnew - vtrial 536 vtrial = vnew 537 538 call mean_fftr(vxc, vmean(1:1), nfft, nfftot, dtset%nspden,& 539 & mpi_comm_sphgrid=mpi_comm_sphgrid) 540 vxcavg = vmean(1) 541 offset = 0 542 end if 543 544 ! Compute mean values of potential and residual 545 call mean_fftr(vnew(1+offset, 1),vnew_mean,nfft,nfftot,dtset%nspden,& 546 & mpi_comm_sphgrid=mpi_comm_sphgrid) 547 call mean_fftr(vresidnew(1+offset, 1),vmean,nfft,nfftot,dtset%nspden,& 548 & mpi_comm_sphgrid=mpi_comm_sphgrid) 549 550 ABI_FREE(vnew) 551 552 ! Subtract the mean of the residual 553 ! Must take into account fixed occupation number in case of spin-polarized 554 do ispden=1,dtset%nspden 555 if (dtset%nspden==2.and.dtset%occopt>=3.and. abs(dtset%spinmagntarget+99.99_dp)<1.0d-10)then 556 vres_mean(ispden)=(vmean(1)+vmean(2))*half 557 else 558 vres_mean(ispden)=vmean(ispden) 559 end if 560 end do 561 562 !$OMP PARALLEL DO COLLAPSE(2) 563 do ispden=1,dtset%nspden 564 do ifft=1,nfft 565 vresidnew(ifft,ispden)=vresidnew(ifft,ispden)-vres_mean(ispden) 566 end do 567 end do 568 569 ! Compute square norm vres2 of potential residual vresid 570 call sqnorm_v(1,nfft,vres2,dtset%nspden,optres,vresidnew(1+offset, 1),mpi_comm_sphgrid=mpi_comm_sphgrid) 571 572 ! Now take care of Vxctau residual (metaGGA) 573 if (with_vxctau) then 574 if (ipositron/=1) vtauresid(:,:)=vxctau(:,:,1)-vtauresid(:,:) 575 if (ipositron==1) vtauresid(:,:)=zero 576 end if 577 578 else ! optres/=0 579 580 ! ------Produce new value of trial potential------------- 581 582 if (.not. wvlbigdft) then 583 !$OMP PARALLEL DO COLLAPSE(2) 584 do ispden=1,min(dtset%nspden,2) 585 do ifft=1,nfft 586 vtrial(ifft,ispden)=vhartr(ifft)+vpsp(ifft)+vxc(ifft,ispden)+vzeeman(ispden)+v_constr_dft_r(ifft,ispden) 587 !vtrial(ifft,ispden)=vtrial(ifft,ispden)+vzeemanHarm(ifft,ispden) 588 if(mod(dtset%fockoptmix,100)==11)vtrial(ifft,ispden)=vtrial(ifft,ispden)+vxc_hybcomp(ifft,ispden) 589 end do 590 end do 591 if(dtset%nspden==4) then 592 !$OMP PARALLEL DO 593 do ifft=1,nfft 594 vtrial(ifft,3:4)=vxc(ifft,3:4)+vzeeman(3:4)+v_constr_dft_r(ifft,3:4) 595 !vtrial(ifft,3:4)=vtrial(ifft,3:4)+vzeemanHarm(ifft,3:4) 596 if(mod(dtset%fockoptmix,100)==11)vtrial(ifft,3:4)=vtrial(ifft,3:4)+vxc_hybcomp(ifft,3:4) 597 end do 598 end if 599 ! Pass vtrial to BigDFT object 600 if(dtset%usewvl==1) then 601 call wvl_vtrial_abi2big(1,vtrial,wvl%den) 602 end if 603 else 604 ! Compute with covering comms the different part of the potential. 605 ABI_MALLOC(xcart,(3, dtset%natom)) 606 call xred2xcart(dtset%natom, rprimd, xcart, xred) 607 call wvl_psitohpsi(dtset%diemix,energies%e_exactX, energies%e_xc, energies%e_hartree, & 608 & energies%e_kinetic, energies%e_localpsp, energies%e_nlpsp_vfock, energies%e_sicdc, & 609 & istep + 1, 1, dtset%iscf, mpi_enreg%me_wvl, & 610 & dtset%natom, dtset%nfft, mpi_enreg%nproc_wvl,& 611 & dtset%nspden,vres2, .true.,energies%e_xcdc, wvl,& 612 & wvlbigdft, xcart, strsxc, vtrial, vxc) 613 ABI_FREE(xcart) 614 ! Compute vxcavg 615 call mean_fftr(vxc, vmean(1:1), nfft, nfftot, dtset%nspden,& 616 & mpi_comm_sphgrid=mpi_comm_sphgrid) 617 vxcavg = vmean(1) 618 end if 619 620 end if 621 622 ABI_FREE(v_constr_dft_r) 623 !ABI_FREE(vzeemanHarm) !SPr: debug for q/=0 magnetic field 624 625 call timab(945,2,tsec) 626 call timab(940,2,tsec) 627 628 DBG_EXIT("COLL") 629 630 end subroutine rhotov