TABLE OF CONTENTS
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.
COPYRIGHT
Copyright (C) 1998-2018 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 . For the initials of contributors, see ~abinit/doc/developers/contributors.txt.
INPUTS
[add_tfw]=flag controling the addition of Weiszacker gradient correction to Thomas-Fermi kin energy 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=informations 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 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) [taug(2,nfftf*dtset%usekden)]=array for Fourier transform of kinetic energy density [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) kxc(nfft,nkxc)=exchange-correlation kernel, needed only if optxc==2. 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)
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.
PARENTS
scfcv
CHILDREN
dotprod_vn,hartre,mag_constr,mean_fftr,psolver_rhohxc,rhohxcpositron rhotoxc,sqnorm_v,timab,wvl_psitohpsi,wvl_vtrial_abi2big,xcdata_init xchybrid_ncpp_cc,xred2xcart
SOURCE
110 #if defined HAVE_CONFIG_H 111 #include "config.h" 112 #endif 113 114 #include "abi_common.h" 115 116 subroutine rhotov(dtset,energies,gprimd,gsqcut,istep,kxc,mpi_enreg,nfft,ngfft,& 117 & nhat,nhatgr,nhatgrdim,nkxc,vresidnew,n3xccc,optene,optres,optxc,& 118 & rhog,rhor,rprimd,strsxc,ucvol,usepaw,usexcnhat,& 119 & vhartr,vnew_mean,vpsp,vres_mean,vres2,vtrial,vxcavg,vxc,wvl,xccc3d,xred,& 120 & electronpositron,taug,taur,vxc_hybcomp,vxctau,add_tfw) ! optional arguments 121 122 use defs_basis 123 use defs_abitypes 124 use defs_wvltypes 125 use m_errors 126 use m_profiling_abi 127 use m_ab7_mixing 128 use m_abi2big 129 use m_xmpi 130 use m_cgtools 131 use m_xcdata 132 133 use m_geometry, only : xred2xcart 134 use m_energies, only : energies_type 135 use m_electronpositron, only : electronpositron_type,electronpositron_calctype 136 use libxc_functionals, only : libxc_functionals_is_hybrid 137 138 !This section has been created automatically by the script Abilint (TD). 139 !Do not modify the following lines by hand. 140 #undef ABI_FUNC 141 #define ABI_FUNC 'rhotov' 142 use interfaces_18_timing 143 use interfaces_56_xc 144 use interfaces_62_poisson 145 use interfaces_62_wvl_wfs 146 use interfaces_67_common, except_this_one => rhotov 147 !End of the abilint section 148 149 implicit none 150 151 !Arguments ------------------------------------ 152 !scalars 153 integer,intent(in) :: n3xccc,nfft,nhatgrdim,nkxc,optene,optres,optxc,usepaw,istep 154 integer,intent(in) :: usexcnhat 155 logical,intent(in),optional :: add_tfw 156 real(dp),intent(in) :: gsqcut,ucvol 157 real(dp),intent(out) :: vres2,vxcavg 158 type(MPI_type),intent(inout) :: mpi_enreg 159 type(dataset_type),intent(in) :: dtset 160 type(electronpositron_type),pointer,optional :: electronpositron 161 type(energies_type),intent(inout) :: energies 162 type(wvl_data), intent(inout) :: wvl 163 !arrays 164 integer,intent(in) :: ngfft(18) 165 real(dp),intent(in) :: gprimd(3,3),nhat(nfft,dtset%nspden*usepaw) 166 real(dp),intent(in) :: nhatgr(nfft,dtset%nspden,3*nhatgrdim),rhog(2,nfft) 167 real(dp),intent(in) :: rprimd(3,3) 168 real(dp),intent(inout) :: rhor(nfft,dtset%nspden),vhartr(nfft),vpsp(nfft) 169 real(dp),intent(inout) :: vtrial(nfft,dtset%nspden),vxc(nfft,dtset%nspden) 170 real(dp),intent(inout) :: xccc3d(n3xccc),xred(3,dtset%natom) 171 real(dp),intent(out) :: kxc(nfft,nkxc),strsxc(6),vnew_mean(dtset%nspden) 172 real(dp),intent(out) :: vres_mean(dtset%nspden),vresidnew(nfft,dtset%nspden) 173 real(dp),intent(in),optional :: taug(2,nfft*dtset%usekden) 174 real(dp),intent(in),optional :: taur(nfft,dtset%nspden*dtset%usekden) 175 real(dp),intent(out),optional :: vxctau(nfft,dtset%nspden,4*dtset%usekden) 176 real(dp),intent(out),optional :: vxc_hybcomp(:,:) ! (nfft,nspden) 177 178 !Local variables------------------------------- 179 !scalars 180 integer :: nk3xc,ifft,ipositron,ispden,nfftot,offset 181 integer :: mpi_comm_sphgrid,ixc_current 182 !integer :: ii,jj,kk,ipt,nx,ny,nz !SPr: debug 183 !real(dp):: rx,ry,rz !SPr: debug 184 real(dp) :: doti,e_xcdc_vxctau 185 logical :: add_tfw_,calc_xcdc,non_magnetic_xc,with_vxctau 186 logical :: is_hybrid_ncpp,wvlbigdft=.false. 187 type(xcdata_type) :: xcdata 188 !arrays 189 real(dp) :: evxc,tsec(2),vmean(dtset%nspden),vzeeman(dtset%nspden) 190 real(dp),allocatable :: rhowk(:,:),Vmagconstr(:,:),vnew(:,:),xcart(:,:) 191 !real(dp),allocatable :: vzeemanHarm(:,:) !SPr: debug Zeeman field q/=0 real space 192 193 ! ********************************************************************* 194 195 DBG_ENTER("COLL") 196 197 call timab(940,1,tsec) 198 199 ! Initialise non_magnetic_xc for rhohxc 200 non_magnetic_xc=(dtset%usepawu==4).or.(dtset%usepawu==14) 201 202 !Check that usekden is not 0 if want to use vxctau 203 with_vxctau = (present(vxctau).and.present(taur).and.(dtset%usekden/=0)) 204 205 !Check if we're in hybrid norm conserving pseudopotential with a core correction 206 is_hybrid_ncpp=(usepaw==0 .and. n3xccc/=0 .and. & 207 & (dtset%ixc==41.or.dtset%ixc==42.or.libxc_functionals_is_hybrid())) 208 209 !If usewvl: wvlbigdft indicates that the BigDFT workflow will be followed 210 wvlbigdft=(dtset%usewvl==1.and.dtset%wvl_bigdft_comp==1) 211 212 !mpi communicator for spherical grid 213 mpi_comm_sphgrid=mpi_enreg%comm_fft 214 if(dtset%usewvl==1) mpi_comm_sphgrid=mpi_enreg%comm_wvl 215 216 !Get size of FFT grid 217 nfftot=PRODUCT(ngfft(1:3)) 218 219 ipositron=0;if (present(electronpositron)) ipositron=electronpositron_calctype(electronpositron) 220 add_tfw_=.false.;if (present(add_tfw)) add_tfw_=add_tfw 221 222 !------Compute Hartree and xc potentials---------------------------------- 223 224 !allocate vnew here. 225 !In wvl: vnew is used at call to wvl_psitohpsi 226 if (optres==0) then 227 ABI_ALLOCATE(vnew,(nfft,dtset%nspden)) 228 vmean(:)=zero ; vnew_mean(:)=zero 229 end if 230 231 if (ipositron/=1) then 232 ! Compute xc potential (separate up and down if spin-polarized) 233 if (dtset%icoulomb == 0 .and. dtset%usewvl == 0) then 234 call hartre(1,gsqcut,usepaw,mpi_enreg,nfft,ngfft,dtset%paral_kgb,rhog,rprimd,vhartr) 235 !Use the proper exchange_correlation energy : either the origin one, or the auxiliary one 236 ixc_current=dtset%ixc 237 if(mod(dtset%fockoptmix,100)==11)ixc_current=dtset%auxc_ixc 238 call xcdata_init(xcdata,dtset=dtset,ixc=ixc_current) 239 240 ! Use the periodic solver to compute Hxc. 241 nk3xc=1 242 !write(80,*) "rhotov" 243 !xccc3d=zero 244 !DEBUG 245 write(std_out,*)' rhotov : is_hybrid_ncpp=',is_hybrid_ncpp 246 write(std_out,*)' rhotov : n3xccc=',n3xccc 247 !ENDDEBUG 248 249 call timab(941,1,tsec) 250 if (ipositron==0) then 251 if(.not.is_hybrid_ncpp .or. mod(dtset%fockoptmix,100)==11)then 252 call rhotoxc(energies%e_xc,kxc,mpi_enreg,nfft,ngfft,& 253 & nhat,usepaw,nhatgr,nhatgrdim,nkxc,nk3xc,non_magnetic_xc,n3xccc,optxc,dtset%paral_kgb,& 254 & rhor,rprimd,strsxc,usexcnhat,vxc,vxcavg,xccc3d,xcdata,& 255 & taug=taug,taur=taur,vhartr=vhartr,vxctau=vxctau,add_tfw=add_tfw_) 256 if(mod(dtset%fockoptmix,100)==11)then 257 energies%e_xc=energies%e_xc*dtset%auxc_scal 258 vxc(:,:)=vxc(:,:)*dtset%auxc_scal 259 end if 260 else 261 call xchybrid_ncpp_cc(dtset,energies%e_xc,mpi_enreg,nfft,ngfft,n3xccc,rhor,rprimd,& 262 & strsxc,vxcavg,xccc3d,vxc=vxc) 263 end if 264 else 265 call rhotoxc(energies%e_xc,kxc,mpi_enreg,nfft,ngfft,& 266 & nhat,usepaw,nhatgr,nhatgrdim,nkxc,nk3xc,non_magnetic_xc,n3xccc,optxc,dtset%paral_kgb,& 267 & rhor,rprimd,strsxc,usexcnhat,vxc,vxcavg,xccc3d,xcdata,& 268 & taug=taug,taur=taur,vhartr=vhartr,vxctau=vxctau,add_tfw=add_tfw_,& 269 & electronpositron=electronpositron) 270 end if 271 call timab(941,2,tsec) 272 273 elseif (.not. wvlbigdft) then 274 ! Use the free boundary solver. 275 call timab(943,1,tsec) 276 call psolver_rhohxc(energies%e_hartree, energies%e_xc, evxc, & 277 & dtset%icoulomb, dtset%ixc, & 278 & mpi_enreg, nfft, & 279 & ngfft, nhat,usepaw,& 280 & dtset%nscforder, dtset%nspden, n3xccc, rhor,rprimd,& 281 & usexcnhat,dtset%usepaw,dtset%usewvl,vhartr, vxc, vxcavg,& 282 & wvl%descr,wvl%den,wvl%e,& 283 & xccc3d,dtset%xclevel,dtset%xc_denpos) 284 call timab(943,2,tsec) 285 end if 286 ! For icoulomb==0 and usewvl Ehartree is calculated in psolver_rhohxc(). 287 ! For PAW we recalculate this since nhat was not taken into account 288 ! in psolver_rhohxc: E_H= int v_H (n+nhat) dr 289 if(.not. wvlbigdft .and. (dtset%icoulomb==0 .or. dtset%usepaw==1 ) ) then 290 291 call timab(942,1,tsec) 292 call dotprod_vn(1,rhor,energies%e_hartree,doti,nfft,nfftot,1,1,vhartr,ucvol,mpi_comm_sphgrid=mpi_comm_sphgrid) 293 energies%e_hartree=half*energies%e_hartree 294 call timab(942,2,tsec) 295 end if 296 else 297 call timab(944,1,tsec) 298 energies%e_hartree=zero;energies%e_xc=zero 299 call rhohxcpositron(electronpositron,gprimd,kxc,mpi_enreg,nfft,ngfft,nhat,nkxc,dtset%nspden,n3xccc,& 300 & dtset%paral_kgb,rhor,strsxc,ucvol,usexcnhat,usepaw,vhartr,vxc,vxcavg,xccc3d,dtset%xc_denpos) 301 call timab(944,2,tsec) 302 end if 303 304 call timab(945,1,tsec) 305 if (ipositron/=0) then 306 call dotprod_vn(1,rhor,electronpositron%e_hartree,doti,& 307 & nfft,nfftot,1,1,electronpositron%vha_ep,ucvol,mpi_comm_sphgrid=mpi_comm_sphgrid) 308 vhartr=vhartr+electronpositron%vha_ep 309 end if 310 311 !------Compute parts of total energy depending on potentials-------- 312 313 if ( (optene==0.or.optene==2 ).and. .not. wvlbigdft) then 314 ! Compute local psp energy energies%e_localpsp 315 call dotprod_vn(1,rhor,energies%e_localpsp,doti,nfft,nfftot,1,1,vpsp,ucvol,& 316 & mpi_comm_sphgrid=mpi_comm_sphgrid) 317 end if 318 319 if(mod(dtset%fockoptmix,100)==11)then 320 if (.not. wvlbigdft) then 321 ! Compute second compensation energy for hybrid functionals 322 call dotprod_vn(1,rhor,energies%e_hybcomp_v,doti,nfft,nfftot,1,1,vxc_hybcomp,ucvol,& 323 & mpi_comm_sphgrid=mpi_comm_sphgrid) 324 end if 325 end if 326 327 calc_xcdc=.false. 328 if (optene==1.or.optene==2) calc_xcdc=.true. 329 if (dtset%usewvl==1.and.dtset%nnsclo>0) calc_xcdc=.true. 330 if (wvlbigdft) calc_xcdc=.false. 331 if (dtset%usefock==1) calc_xcdc=.true. 332 333 if (calc_xcdc) then 334 335 ! Compute double-counting XC energy energies%e_xcdc 336 if (ipositron/=1) then 337 if (usepaw==0.or.usexcnhat/=0) then 338 call dotprod_vn(1,rhor,energies%e_xcdc,doti,nfft,nfftot,dtset%nspden,1,vxc,ucvol,& 339 & mpi_comm_sphgrid=mpi_comm_sphgrid) 340 if (with_vxctau)then 341 call dotprod_vn(1,taur,e_xcdc_vxctau,doti,nfft,nfftot,dtset%nspden,1,vxctau(:,:,1),& 342 & ucvol,mpi_comm_sphgrid=mpi_comm_sphgrid) 343 energies%e_xcdc=energies%e_xcdc+e_xcdc_vxctau 344 end if 345 else 346 ABI_ALLOCATE(rhowk,(nfft,dtset%nspden)) 347 rhowk=rhor-nhat 348 call dotprod_vn(1,rhowk,energies%e_xcdc,doti,nfft,nfftot,dtset%nspden,1,vxc,ucvol,& 349 & mpi_comm_sphgrid=mpi_comm_sphgrid) 350 ABI_DEALLOCATE(rhowk) 351 end if 352 if (ipositron==2) energies%e_xcdc=energies%e_xcdc-electronpositron%e_xcdc 353 else 354 energies%e_xcdc=zero 355 end if 356 357 end if 358 359 !------Produce residual vector and square norm of it------------- 360 !(only if requested ; if optres==0) 361 362 !Set up array for Zeeman field 363 !EB vzeeman(:) = factor*( Hz, Hx+iHy; Hx-iHy, -Hz) 364 !EB factor = -g/2 * mu_B * mu_0 = -1/2*B in a.u. 365 !EB <-- vzeeman might have to be allocated correctly --> to be checked 366 ! vzeeman = 1/2 ( -B_z, -B_x + iB_y ; -B_x - iB_y , B_z) 367 vzeeman(:) = zero 368 ! ABI_ALLOCATE(vzeemanHarm,(nfft,dtset%nspden)) ! SPr: debug stuff 369 ! vzeemanHarm(:,:) = zero ! 370 if (any(abs(dtset%zeemanfield(:))>tol8)) then 371 if(dtset%nspden==2)then 372 ! EB The collinear case has to be checked : 373 ! EB Is it vzeeman(1) or (2) that has to be added here? to be checked in setvtr and energy as well 374 ! SPr: the density components are: rhor(1) => n_upup + n_dwndwn 375 ! rhor(2) => n_upup 376 ! the convention for the potential components is different: 377 ! v(1) => v_upup 378 ! v(2) => v_dndn 379 ! verified by comparing collinear and non-collinear calculations 380 381 vzeeman(1) =-half*dtset%zeemanfield(3) ! v_upup 382 vzeeman(2) = half*dtset%zeemanfield(3) ! v_dndn 383 384 !vzeeman(1) = zero ! v_upup 385 !vzeeman(2) = zero ! v_dndn 386 387 !nx=ngfft(1); ny=ngfft(2); nz=ngfft(3) 388 !do kk=0,nz-1 389 ! do jj=0,ny-1 390 ! do ii=0,nx-1 391 ! ipt=1+ii+nx*(jj+ny*kk) 392 ! !rx=(dble(ii)/nx)*rprimd(1,1)+(dble(jj)/ny)*rprimd(1,2)+(dble(kk)/nz)*rprimd(1,3) 393 ! !ry=(dble(ii)/nx)*rprimd(2,1)+(dble(jj)/ny)*rprimd(2,2)+(dble(kk)/nz)*rprimd(2,3) 394 ! !rz=(dble(ii)/nx)*rprimd(3,1)+(dble(jj)/ny)*rprimd(3,2)+(dble(kk)/nz)*rprimd(3,3) 395 ! vzeemanHarm(ipt,1)= -half*dtset%zeemanfield(3)*cos(2*PI*(dble(ii)/dble(nx))) 396 ! vzeemanHarm(ipt,2)= half*dtset%zeemanfield(3)*cos(2*PI*(dble(ii)/dble(nx))) 397 ! end do 398 ! end do 399 !end do 400 401 else if(dtset%nspden==4)then 402 403 vzeeman(1)=-half*dtset%zeemanfield(3) ! v_upup 404 vzeeman(2)= half*dtset%zeemanfield(3) ! v_dndn 405 vzeeman(3)=-half*dtset%zeemanfield(1) ! Re(v_updn) 406 vzeeman(4)= half*dtset%zeemanfield(2) ! Im(v_updn) 407 408 !vzeeman(1)=0.0 409 !vzeeman(2)=0.0 410 !vzeeman(3)=0.0 411 !vzeeman(4)=0.0 412 413 !nx=ngfft(1); ny=ngfft(2); nz=ngfft(3) 414 !do kk=0,nz-1 415 ! do jj=0,ny-1 416 ! do ii=0,nx-1 417 ! ipt=1+ii+nx*(jj+ny*kk) 418 ! !rx=(dble(ii)/nx)*rprimd(1,1)+(dble(jj)/ny)*rprimd(1,2)+(dble(kk)/nz)*rprimd(1,3) 419 ! !ry=(dble(ii)/nx)*rprimd(2,1)+(dble(jj)/ny)*rprimd(2,2)+(dble(kk)/nz)*rprimd(2,3) 420 ! !rz=(dble(ii)/nx)*rprimd(3,1)+(dble(jj)/ny)*rprimd(3,2)+(dble(kk)/nz)*rprimd(3,3) 421 ! vzeemanHarm(ipt,1)= -half*dtset%zeemanfield(3)*cos(2*PI*(dble(ii)/dble(nx))) 422 ! vzeemanHarm(ipt,2)= half*dtset%zeemanfield(3)*cos(2*PI*(dble(ii)/dble(nx))) 423 ! vzeemanHarm(ipt,3)= -half*dtset%zeemanfield(1)*cos(2*PI*(dble(ii)/dble(nx))) 424 ! vzeemanHarm(ipt,4)= half*dtset%zeemanfield(2)*cos(2*PI*(dble(ii)/dble(nx))) 425 ! end do 426 ! end do 427 !end do 428 429 end if 430 end if 431 432 !Compute the constrained potential for the magnetic moments 433 ABI_ALLOCATE(Vmagconstr, (nfft,dtset%nspden)) 434 if (dtset%magconon==1.or.dtset%magconon==2) then 435 Vmagconstr = zero 436 call mag_constr(dtset%natom,dtset%spinat,dtset%nspden,dtset%magconon,dtset%magcon_lambda,rprimd, & 437 & mpi_enreg,nfft,ngfft,dtset%ntypat,dtset%ratsph,rhor,dtset%typat,Vmagconstr,xred) 438 else 439 ! NOTE: mjv 25 May 2013 added this for ibm6 - otherwise gives NaN in vnew after 440 ! the addition below, in case magconon==0 and for certain libxc 441 ! functionals!! May need to copy this to setvtr.F90 if the same effect appears. 442 ! should check libxc test with a proper memory checker (valgrind). 443 do ispden=1,dtset%nspden 444 do ifft=1,nfft 445 Vmagconstr(ifft,ispden) = zero 446 end do 447 end do 448 end if 449 450 if (optres==0) then 451 452 453 ! ------ Compute potential residual ------------- 454 455 if (.not. wvlbigdft) then 456 !$OMP PARALLEL DO COLLAPSE(2) 457 do ispden=1,min(dtset%nspden,2) 458 do ifft=1,nfft 459 vnew(ifft,ispden)=vhartr(ifft)+vpsp(ifft)+vxc(ifft,ispden)+vzeeman(ispden)+Vmagconstr(ifft,ispden) 460 !vnew(ifft,ispden)=vnew(ifft,ispden)+vzeemanHarm(ifft,ispden) 461 if(mod(dtset%fockoptmix,100)==11)vnew(ifft,ispden)=vnew(ifft,ispden)+vxc_hybcomp(ifft,ispden) 462 vresidnew(ifft,ispden)=vnew(ifft,ispden)-vtrial(ifft,ispden) 463 end do 464 end do 465 if(dtset%nspden==4)then 466 !$OMP PARALLEL DO COLLAPSE(2) 467 do ispden=3,4 468 do ifft=1,nfft 469 vnew(ifft,ispden)=vxc(ifft,ispden)+vzeeman(ispden)+Vmagconstr(ifft,ispden) 470 !vnew(ifft,ispden)=vnew(ifft,ispden)+vzeemanHarm(ifft,ispden) 471 if(mod(dtset%fockoptmix,100)==11)vnew(ifft,ispden)=vnew(ifft,ispden)+vxc_hybcomp(ifft,ispden) 472 vresidnew(ifft,ispden)=vnew(ifft,ispden)-vtrial(ifft,ispden) 473 end do 474 end do 475 end if 476 477 offset = 0 478 479 if (dtset%iscf==0) vtrial=vnew 480 481 ! Pass vtrial to BigDFT object 482 if(dtset%usewvl==1) then 483 call wvl_vtrial_abi2big(1,vnew,wvl%den) 484 ! call wvl_vtrial_abi2big(1,vtrial,wvl%den) 485 end if 486 487 else 488 ! Compute with covering comms the different part of the potential. 489 ! only for wvlbigdft 490 ABI_ALLOCATE(xcart,(3, dtset%natom)) 491 call xred2xcart(dtset%natom, rprimd, xcart, xred) 492 call wvl_psitohpsi(dtset%diemix,energies%e_exactX, energies%e_xc, energies%e_hartree, & 493 & energies%e_kinetic, energies%e_localpsp, energies%e_nonlocalpsp, energies%e_sicdc, & 494 & istep + 1, 1, dtset%iscf, mpi_enreg%me_wvl, dtset%natom, dtset%nfft,& 495 & mpi_enreg%nproc_wvl, dtset%nspden, & 496 & vres2, .true., energies%e_xcdc, wvl,& 497 & wvlbigdft, xcart, strsxc,vtrial=vnew,vxc=vxc) 498 ABI_DEALLOCATE(xcart) 499 500 vresidnew = vnew - vtrial 501 vtrial = vnew 502 503 call mean_fftr(vxc, vmean(1:1), nfft, nfftot, dtset%nspden,& 504 & mpi_comm_sphgrid=mpi_comm_sphgrid) 505 vxcavg = vmean(1) 506 offset = 0 507 end if 508 509 ! Compute mean values of potential and residual 510 call mean_fftr(vnew(1+offset, 1),vnew_mean,nfft,nfftot,dtset%nspden,& 511 & mpi_comm_sphgrid=mpi_comm_sphgrid) 512 call mean_fftr(vresidnew(1+offset, 1),vmean,nfft,nfftot,dtset%nspden,& 513 & mpi_comm_sphgrid=mpi_comm_sphgrid) 514 515 ABI_DEALLOCATE(vnew) 516 517 ! Subtract the mean of the residual 518 ! Must take into account fixed occupation number in case of spin-polarized 519 do ispden=1,dtset%nspden 520 if (dtset%nspden==2.and.dtset%occopt>=3.and. abs(dtset%spinmagntarget+99.99_dp)<1.0d-10)then 521 vres_mean(ispden)=(vmean(1)+vmean(2))*half 522 else 523 vres_mean(ispden)=vmean(ispden) 524 end if 525 end do 526 527 !$OMP PARALLEL DO COLLAPSE(2) 528 do ispden=1,dtset%nspden 529 do ifft=1,nfft 530 vresidnew(ifft,ispden)=vresidnew(ifft,ispden)-vres_mean(ispden) 531 end do 532 end do 533 534 ! Compute square norm vres2 of potential residual vresid 535 call sqnorm_v(1,nfft,vres2,dtset%nspden,optres,vresidnew(1+offset, 1),mpi_comm_sphgrid=mpi_comm_sphgrid) 536 537 else ! optres/=0 538 539 ! ------Produce new value of trial potential------------- 540 541 if (.not. wvlbigdft) then 542 !$OMP PARALLEL DO COLLAPSE(2) 543 do ispden=1,min(dtset%nspden,2) 544 do ifft=1,nfft 545 vtrial(ifft,ispden)=vhartr(ifft)+vpsp(ifft)+vxc(ifft,ispden)+vzeeman(ispden)+Vmagconstr(ifft,ispden) 546 !vtrial(ifft,ispden)=vtrial(ifft,ispden)+vzeemanHarm(ifft,ispden) 547 if(mod(dtset%fockoptmix,100)==11)vtrial(ifft,ispden)=vtrial(ifft,ispden)+vxc_hybcomp(ifft,ispden) 548 end do 549 end do 550 if(dtset%nspden==4) then 551 !$OMP PARALLEL DO 552 do ifft=1,nfft 553 vtrial(ifft,3:4)=vxc(ifft,3:4)+vzeeman(3:4)+Vmagconstr(ifft,3:4) 554 !vtrial(ifft,3:4)=vtrial(ifft,3:4)+vzeemanHarm(ifft,3:4) 555 if(mod(dtset%fockoptmix,100)==11)vtrial(ifft,3:4)=vtrial(ifft,3:4)+vxc_hybcomp(ifft,3:4) 556 end do 557 end if 558 ! Pass vtrial to BigDFT object 559 if(dtset%usewvl==1) then 560 call wvl_vtrial_abi2big(1,vtrial,wvl%den) 561 end if 562 else 563 ! Compute with covering comms the different part of the potential. 564 ABI_ALLOCATE(xcart,(3, dtset%natom)) 565 call xred2xcart(dtset%natom, rprimd, xcart, xred) 566 call wvl_psitohpsi(dtset%diemix,energies%e_exactX, energies%e_xc, energies%e_hartree, & 567 & energies%e_kinetic, energies%e_localpsp, energies%e_nonlocalpsp, energies%e_sicdc, & 568 & istep + 1, 1, dtset%iscf, mpi_enreg%me_wvl, & 569 & dtset%natom, dtset%nfft, mpi_enreg%nproc_wvl,& 570 & dtset%nspden,vres2, .true.,energies%e_xcdc, wvl,& 571 & wvlbigdft, xcart, strsxc, vtrial, vxc) 572 ABI_DEALLOCATE(xcart) 573 ! Compute vxcavg 574 call mean_fftr(vxc, vmean(1:1), nfft, nfftot, dtset%nspden,& 575 & mpi_comm_sphgrid=mpi_comm_sphgrid) 576 vxcavg = vmean(1) 577 end if 578 579 end if 580 581 ABI_DEALLOCATE(Vmagconstr) 582 !ABI_DEALLOCATE(vzeemanHarm) !SPr: debug for q/=0 magnetic field 583 584 call timab(945,2,tsec) 585 call timab(940,2,tsec) 586 587 DBG_EXIT("COLL") 588 589 end subroutine rhotov