TABLE OF CONTENTS
ABINIT/nres2vres [ Functions ]
NAME
nres2vres
FUNCTION
Convert a density residual into a potential residual
using a first order formula:
V^res(r)=dV/dn.n^res(r)
=V_hartree(n^res)(r) + Kxc.n^res(r)
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (MT) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
INPUTS
dtset <type(dataset_type)>=all input variables in this dataset | icoulomb=0 periodic treatment of Hartree potential, 1 use of Poisson solver | natom= number of atoms in cell | nspden=number of spin-density components | ntypat=number of atom types | typat(natom)=type (integer) for each atom gsqcut=cutoff value on G**2 for sphere inside fft box izero=if 1, unbalanced components of Vhartree(g) are set to zero kxc(nfft,nkxc)=exchange-correlation kernel, needed only if nkxc>0 mpi_enreg=information about MPI parallelization my_natom=number of atoms treated by current processor nfft=(effective) number of FFT grid points (for this processor) ngfft(18)=contain all needed information about 3D FFT nhat(nfft,nspden*usepaw)= -PAW only- compensation density nkxc=second dimension of the array kxc, see rhotoxc.F90 for a description nresid(nfft,nspden)= the input density residual n3xccc=dimension of the xccc3d array (0 or nfft). optnc=option for non-collinear magnetism (nspden=4): 1: the whole 2x2 Vres matrix is computed 2: only Vres^{11} and Vres^{22} are computed optxc=0 if LDA part of XC kernel has only to be taken into account (even for GGA) 1 if XC kernel has to be fully taken into -1 if XC kernel does not have to be taken into account pawang <type(pawang_type)>=paw angular mesh and related data pawfgrtab(my_natom*usepaw) <type(pawfgrtab_type)>=atomic data given on fine rectangular grid pawrhoij(my_natom*usepaw) <type(pawrhoij_type)>= paw rhoij occupancies and related data pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data rhor(nfft,nspden)=electron density in real space (used only if Kxc was not computed before) rprimd(3,3)=dimensional primitive translation vectors (bohr) usepaw= 0 for non paw calculation; =1 for paw calculation xccc3d(n3xccc)=3D core electron density for XC core correction (bohr^-3) xred(3,natom)=reduced dimensionless atomic coordinates === optional inputs === vxc(cplex*nfft,nspden)=XC GS potential
OUTPUT
vresid(nfft,nspden)= the output potential residual
PARENTS
etotfor,forstr
CHILDREN
dfpt_mkvxc,dfpt_mkvxc_noncoll,fourdp,hartre,metric,pawmknhat psolver_hartree,rhotoxc,xcdata_init
SOURCE
69 #if defined HAVE_CONFIG_H 70 #include "config.h" 71 #endif 72 73 #include "abi_common.h" 74 75 subroutine nres2vres(dtset,gsqcut,izero,kxc,mpi_enreg,my_natom,nfft,ngfft,nhat,& 76 & nkxc,nresid,n3xccc,optnc,optxc,pawang,pawfgrtab,pawrhoij,pawtab,& 77 & rhor,rprimd,usepaw,vresid,xccc3d,xred,vxc) 78 79 use defs_basis 80 use defs_abitypes 81 use m_errors 82 use m_xmpi 83 use m_xcdata 84 85 use m_geometry, only : metric 86 use m_pawang, only : pawang_type 87 use m_pawtab, only : pawtab_type 88 use m_pawfgrtab,only : pawfgrtab_type 89 use m_pawrhoij, only : pawrhoij_type 90 91 !This section has been created automatically by the script Abilint (TD). 92 !Do not modify the following lines by hand. 93 #undef ABI_FUNC 94 #define ABI_FUNC 'nres2vres' 95 use interfaces_53_ffts 96 use interfaces_56_xc 97 use interfaces_62_poisson 98 use interfaces_65_paw 99 !End of the abilint section 100 101 implicit none 102 103 !Arguments ------------------------------------ 104 !scalars 105 integer,intent(in) :: izero,my_natom,n3xccc,nfft,nkxc,optnc,optxc,usepaw 106 real(dp),intent(in) :: gsqcut 107 type(MPI_type),intent(in) :: mpi_enreg 108 type(dataset_type),intent(in) :: dtset 109 type(pawang_type),intent(in) :: pawang 110 !arrays 111 integer,intent(in) :: ngfft(18) 112 real(dp),intent(in) :: kxc(nfft,nkxc),nresid(nfft,dtset%nspden) 113 real(dp),intent(in) :: rhor(nfft,dtset%nspden),rprimd(3,3),xccc3d(n3xccc),xred(3,dtset%natom) 114 real(dp),intent(inout) :: nhat(nfft,dtset%nspden*usepaw) 115 real(dp),intent(out) :: vresid(nfft,dtset%nspden) 116 type(pawfgrtab_type),intent(inout) :: pawfgrtab(my_natom*usepaw) 117 type(pawrhoij_type),intent(in) :: pawrhoij(my_natom*usepaw) 118 type(pawtab_type),intent(in) :: pawtab(dtset%ntypat*usepaw) 119 real(dp),intent(in) :: vxc(nfft,dtset%nspden) !FR TODO:cplex? 120 121 !Local variables------------------------------- 122 !scalars 123 integer :: cplex,ider,idir,ipert,ispden,nhatgrdim,nkxc_cur,option,me,nproc,comm,usexcnhat 124 logical :: has_nkxc_gga 125 logical :: non_magnetic_xc 126 real(dp) :: dum,energy,m_norm_min,ucvol,vxcavg 127 character(len=500) :: message 128 type(xcdata_type) :: xcdata 129 !arrays 130 integer :: nk3xc 131 real(dp) :: dummy6(6),gmet(3,3),gprimd(3,3),qq(3),rmet(3,3) 132 real(dp),allocatable :: dummy(:),kxc_cur(:,:),nhatgr(:,:,:) 133 real(dp),allocatable :: nresg(:,:),rhor0(:,:),vhres(:) 134 135 ! ************************************************************************* 136 137 !Compatibility tests: 138 has_nkxc_gga=(nkxc==7.or.nkxc==19) 139 140 ! Initialise non_magnetic_xc for rhohxc 141 non_magnetic_xc=(dtset%usepawu==4).or.(dtset%usepawu==14) 142 143 if(optxc<-1.or.optxc>1)then 144 write(message,'(a,i0)')' Wrong value for optxc ',optxc 145 MSG_BUG(message) 146 end if 147 148 if((optnc/=1.and.optnc/=2).or.(dtset%nspden/=4.and.optnc/=1))then 149 write(message,'(a,i0)')' Wrong value for optnc ',optnc 150 MSG_BUG(message) 151 end if 152 153 if(dtset%icoulomb==1.and.optxc/=-1)then 154 write(message,'(a)')' This routine is not compatible with icoulomb==1 and optxc/=-1 !' 155 MSG_BUG(message) 156 end if 157 158 if(dtset%nspden==4.and.dtset%xclevel==2.and.optxc==1.and.(.not.has_nkxc_gga))then 159 MSG_ERROR(' Wrong values for optxc and nkxc !') 160 end if 161 162 qq=zero 163 nkxc_cur=0 164 m_norm_min=EPSILON(0.0_dp)**2 165 usexcnhat=0;if (usepaw==1) usexcnhat=maxval(pawtab(1:dtset%ntypat)%usexcnhat) 166 if (dtset%xclevel==1.or.optxc==0) nkxc_cur= 2*min(dtset%nspden,2)-1 ! LDA: nkxc=1,3 167 if (dtset%xclevel==2.and.optxc==1)nkxc_cur=12*min(dtset%nspden,2)-5 ! GGA: nkxc=7,19 168 ABI_ALLOCATE(vhres,(nfft)) 169 170 !Compute different geometric tensor, as well as ucvol, from rprimd 171 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 172 173 !Compute density residual in reciprocal space 174 if (dtset%icoulomb==0) then 175 ABI_ALLOCATE(nresg,(2,nfft)) 176 ABI_ALLOCATE(dummy,(nfft)) 177 dummy(:)=nresid(:,1) 178 call fourdp(1,nresg,dummy,-1,mpi_enreg,nfft,ngfft,dtset%paral_kgb,0) 179 ABI_DEALLOCATE(dummy) 180 end if 181 182 !For GGA, has to recompute gradients of nhat 183 nhatgrdim=0 184 if ((nkxc==nkxc_cur.and.has_nkxc_gga).or.(optxc==-1.and.has_nkxc_gga).or.& 185 & (optxc/=-1.and.nkxc/=nkxc_cur)) then 186 if (usepaw==1.and.dtset%xclevel==2.and.usexcnhat>0.and.dtset%pawnhatxc>0) then 187 nhatgrdim=1 188 ABI_ALLOCATE(nhatgr,(nfft,dtset%nspden,3)) 189 ider=1;cplex=1;ipert=0;idir=0 190 call pawmknhat(dum,cplex,ider,idir,ipert,izero,gprimd,my_natom,dtset%natom,& 191 & nfft,ngfft,nhatgrdim,dtset%nspden,dtset%ntypat,pawang,pawfgrtab,& 192 & nhatgr,nhat,pawrhoij,pawrhoij,pawtab,qq,rprimd,ucvol,dtset%usewvl,xred,& 193 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,& 194 & comm_fft=mpi_enreg%comm_fft,paral_kgb=dtset%paral_kgb,me_g0=mpi_enreg%me_g0,& 195 & distribfft=mpi_enreg%distribfft,mpi_comm_wvl=mpi_enreg%comm_wvl) 196 else 197 ABI_ALLOCATE(nhatgr,(0,0,0)) 198 end if 199 else 200 ABI_ALLOCATE(nhatgr,(0,0,0)) 201 end if 202 203 ABI_ALLOCATE(dummy,(0)) 204 !First case: Kxc has already been computed 205 !----------------------------------------- 206 if (nkxc==nkxc_cur.or.optxc==-1) then 207 208 ! Compute VH(n^res)(r) 209 if (dtset%icoulomb == 0) then 210 call hartre(1,gsqcut,izero,mpi_enreg,nfft,ngfft,dtset%paral_kgb,nresg,rprimd,vhres) 211 else 212 comm=mpi_enreg%comm_cell 213 nproc=xmpi_comm_size(comm) 214 me=xmpi_comm_rank(comm) 215 call psolver_hartree(energy, (/ rprimd(1,1) / dtset%ngfft(1), & 216 & rprimd(2,2) / dtset%ngfft(2), rprimd(3,3) / dtset%ngfft(3) /), dtset%icoulomb, & 217 & me, comm, dtset%nfft, dtset%ngfft(1:3), nproc, dtset%nscforder, dtset%nspden, & 218 & nresid(:,1), vhres, dtset%usewvl) 219 end if 220 221 ! Compute Kxc(r).n^res(r) 222 if (optxc/=-1) then 223 224 ! Collinear magnetism or non-polarized 225 if (dtset%nspden/=4) then 226 !Note: imposing usexcnhat=1 avoid nhat to be substracted 227 call dfpt_mkvxc(1,dtset%ixc,kxc,mpi_enreg,nfft,ngfft,nhat,usepaw,nhatgr,nhatgrdim,& 228 & nkxc,dtset%nspden,0,2,dtset%paral_kgb,qq,nresid,rprimd,1,vresid,dummy) 229 else 230 !FR call routine for Non-collinear magnetism 231 ABI_ALLOCATE(rhor0,(nfft,dtset%nspden)) 232 rhor0(:,:)=rhor(:,:)-nresid(:,:) 233 !Note: imposing usexcnhat=1 avoid nhat to be substracted 234 call dfpt_mkvxc_noncoll(1,dtset%ixc,kxc,mpi_enreg,nfft,ngfft,nhat,usepaw,nhat,usepaw,nhatgr,nhatgrdim,& 235 & nkxc,dtset%nspden,0,2,2,dtset%paral_kgb,qq,rhor0,nresid,rprimd,1,vxc,vresid,xccc3d) 236 ABI_DEALLOCATE(rhor0) 237 end if 238 239 else 240 vresid=zero 241 end if 242 243 end if 244 245 !2nd case: Kxc has to be computed 246 !-------------------------------- 247 if (nkxc/=nkxc_cur.and.optxc/=-1) then 248 249 ! Has to use the "initial" density to compute Kxc 250 ABI_ALLOCATE(rhor0,(nfft,dtset%nspden)) 251 rhor0(:,:)=rhor(:,:)-nresid(:,:) 252 253 ! Compute VH(n^res) and XC kernel (Kxc) together 254 ABI_ALLOCATE(kxc_cur,(nfft,nkxc_cur)) 255 256 option=2;if (dtset%xclevel==2.and.optxc==0) option=12 257 258 call hartre(1,gsqcut,izero,mpi_enreg,nfft,ngfft,dtset%paral_kgb,nresg,rprimd,vhres) 259 call xcdata_init(xcdata,dtset=dtset) 260 261 ! To be adjusted for the call to rhotoxc 262 nk3xc=1 263 call rhotoxc(energy,kxc_cur,mpi_enreg,nfft,ngfft,& 264 & nhat,usepaw,nhatgr,nhatgrdim,nkxc_cur,nk3xc,non_magnetic_xc,n3xccc,option,dtset%paral_kgb,& 265 & rhor0,rprimd,dummy6,usexcnhat,vresid,vxcavg,xccc3d,xcdata,vhartr=vhres) !vresid=work space 266 if (dtset%nspden/=4) then 267 ABI_DEALLOCATE(rhor0) 268 end if 269 270 ! Compute Kxc(r).n^res(r) 271 272 ! Collinear magnetism or non-polarized 273 if (dtset%nspden/=4) then 274 call dfpt_mkvxc(1,dtset%ixc,kxc_cur,mpi_enreg,nfft,ngfft,nhat,usepaw,nhatgr,nhatgrdim,& 275 & nkxc_cur,dtset%nspden,0,2,dtset%paral_kgb,qq,nresid,rprimd,1,vresid,dummy) 276 else 277 !FR call routine for Non-collinear magnetism 278 ABI_ALLOCATE(rhor0,(nfft,dtset%nspden)) 279 rhor0(:,:)=rhor(:,:)-nresid(:,:) 280 call dfpt_mkvxc_noncoll(1,dtset%ixc,kxc_cur,mpi_enreg,nfft,ngfft,nhat,usepaw,nhat,usepaw,nhatgr,nhatgrdim,& 281 & nkxc,dtset%nspden,0,2,2,dtset%paral_kgb,qq,rhor0,nresid,rprimd,1,vxc,vresid,xccc3d) 282 ABI_DEALLOCATE(rhor0) 283 end if 284 285 ABI_DEALLOCATE(kxc_cur) 286 end if 287 288 !if (nhatgrdim>0) then 289 ABI_DEALLOCATE(nhatgr) 290 !end if 291 292 !Assemble potential residual: V^res(r)=VH(n^res)(r) + Kxc(r).n^res(r) 293 !-------------------------------------------------------------------- 294 do ispden=1,dtset%nspden/optnc 295 vresid(:,ispden)=vresid(:,ispden)+vhres(:) 296 end do 297 298 if (dtset%icoulomb==0) then 299 ABI_DEALLOCATE(nresg) 300 end if 301 ABI_DEALLOCATE(vhres) 302 ABI_DEALLOCATE(dummy) 303 304 end subroutine nres2vres