TABLE OF CONTENTS
ABINIT/m_mkcore_wvl [ Modules ]
NAME
m_mkcore_wvl
FUNCTION
COPYRIGHT
Copyright (C) 2016-2024 ABINIT group (MT,TRangel) 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
16 #if defined HAVE_CONFIG_H 17 #include "config.h" 18 #endif 19 20 #include "abi_common.h" 21 22 module m_mkcore_wvl 23 24 use defs_basis 25 use m_abicore 26 use m_errors 27 use m_splines 28 use m_xmpi 29 use defs_wvltypes 30 31 use m_time, only : timab 32 use m_sort, only : sort_dp 33 use m_geometry, only : xcart2xred, xred2xcart, metric, strconv 34 use m_paw_numeric, only : paw_splint, paw_splint_der 35 use m_pawrad, only : pawrad_type, pawrad_init, pawrad_free 36 use m_pawtab, only : pawtab_type 37 use m_drivexc, only : mkdenpos 38 39 private
ABINIT/mkcore_inner [ Functions ]
NAME
mkcore_inner
FUNCTION
FIXME: add description.
INPUTS
argin(sizein)=description
OUTPUT
argout(sizeout)=description
SIDE EFFECTS
NOTES
SOURCE
900 subroutine mkcore_inner(corfra,core_mesh,dyfrx2,grxc1,grxc2,grxc3,ifftsph,msz,natom,ncmax,nfft,& 901 & nfgd,nfgd_r0,nspden,n3xccc,option,pawtab,rmet,rr,strdia,vxc,xccc3d,& 902 & rred) ! optional argument 903 904 !Arguments ------------------------------------ 905 !scalars 906 integer ,intent(in) :: msz,natom,ncmax,nfft,nfgd,nfgd_r0,nspden,n3xccc,option 907 real(dp),intent(out) :: grxc1,grxc2,grxc3,strdia 908 !arrays 909 integer,intent(in) :: ifftsph(ncmax) 910 real(dp),intent(in) :: rmet(3,3),rr(ncmax),vxc(nfft,nspden) 911 real(dp),intent(in),optional :: rred(:,:) 912 real(dp),intent(inout) :: corfra(3,3),dyfrx2(3,3,natom),xccc3d(n3xccc) 913 type(pawrad_type),intent(in) :: core_mesh 914 type(pawtab_type),intent(in) :: pawtab 915 916 !Local variables------------------------------- 917 !scalars 918 integer :: iatom,ifgd,ii,jj,mu,nu 919 character(len=500) :: message 920 real(dp) :: term,term2 921 !arrays 922 integer :: iindex(nfgd) 923 real(dp) :: tcore(nfgd),dtcore(nfgd),rr_tmp(nfgd) 924 real(dp),allocatable :: d2tcore(:) 925 926 ! ************************************************************************* 927 928 !Checks 929 if(nspden >1) then 930 write(message, '(a)')'mkcore_inner: this is not yet generalized to npsden>1' 931 ABI_ERROR(message) 932 end if 933 if (present(rred)) then 934 if (option>1.and.size(rred)/=3*ncmax) then 935 write(message, '(a)')'mkcore_inner: incorrect size for rred' 936 ABI_BUG(message) 937 end if 938 else if (option>1) then 939 write(message, '(a)')'mkcore_inner: rred is not present and option >1' 940 ABI_BUG(message) 941 end if 942 943 !Retrieve values of pseudo core density (and derivative) 944 rr_tmp=rr(1:nfgd) 945 iindex(1:nfgd)=(/(ii,ii=1,nfgd)/) 946 call sort_dp(nfgd,rr_tmp,iindex(1:nfgd),tol16) 947 if (option==1.or.option==3) then 948 call paw_splint(msz,core_mesh%rad,pawtab%tcoredens(:,1),pawtab%tcoredens(:,3),& 949 & nfgd,rr_tmp,tcore) 950 end if 951 if (option>=2) then 952 call paw_splint_der(msz,core_mesh%rad,pawtab%tcoredens(:,1),pawtab%tcoredens(:,3),& 953 & nfgd,rr_tmp,dtcore) 954 end if 955 956 !Accumulate contributions to core density on the entire cell 957 if (option==1) then 958 do ifgd=1,nfgd 959 ii=iindex(ifgd);jj=ifftsph(ii) 960 xccc3d(jj)=xccc3d(jj)+tcore(ifgd) 961 end do 962 963 !Accumulate contributions to Exc gradients 964 else if (option==2) then 965 do ifgd=1,nfgd 966 ii=iindex(ifgd);jj=ifftsph(ii) 967 term=vxc(jj,1)*dtcore(ifgd)/rr_tmp(ifgd) 968 grxc1=grxc1-rred(1,ii)*term 969 grxc2=grxc2-rred(2,ii)*term 970 grxc3=grxc3-rred(3,ii)*term 971 end do 972 973 !Accumulate contributions to stress tensor 974 else if (option==3) then 975 ! Write out the 6 symmetric components 976 do ifgd=1,nfgd 977 ii=iindex(ifgd);jj=ifftsph(ii) 978 term=vxc(jj,1)*dtcore(ifgd)/rr_tmp(ifgd) 979 corfra(1,1)=corfra(1,1)+term*rred(1,ifgd)**2 980 corfra(2,2)=corfra(2,2)+term*rred(2,ifgd)**2 981 corfra(3,3)=corfra(3,3)+term*rred(3,ifgd)**2 982 corfra(3,2)=corfra(3,2)+term*rred(3,ifgd)*rred(2,ifgd) 983 corfra(3,1)=corfra(3,1)+term*rred(3,ifgd)*rred(1,ifgd) 984 corfra(2,1)=corfra(2,1)+term*rred(2,ifgd)*rred(1,ifgd) 985 end do 986 ! Also compute diagonal term 987 do ifgd=1,nfgd 988 ii=iindex(ifgd);jj=ifftsph(ii) 989 strdia=strdia+vxc(jj,1)*tcore(ii) 990 end do 991 992 !Compute frozen-wf contribution to Dynamical matrix 993 else if (option==4) then 994 ABI_MALLOC(d2tcore,(nfgd)) 995 if(nfgd>0) then 996 ! Evaluate spline fit of 2nd der of pseudo core density 997 call paw_splint(msz,core_mesh%rad,pawtab%tcoredens(:,3),pawtab%tcoredens(:,5),& 998 & nfgd,rr_tmp,d2tcore) 999 do ifgd=1,nfgd 1000 ii=iindex(ifgd);jj=ifftsph(ii) 1001 term=vxc(jj,1)*dtcore(ifgd)/rr_tmp(ifgd) 1002 term2=term*d2tcore(ifgd)/rr_tmp(ifgd) 1003 do mu=1,3 1004 do nu=1,3 1005 dyfrx2(mu,nu,iatom)=dyfrx2(mu,nu,iatom)& 1006 & +(term2-term/rr_tmp(iindex(ifgd))**2)& 1007 & *rred(mu,iatom)*rred(nu,iatom)+term*rmet(mu,nu) 1008 end do 1009 end do 1010 end do 1011 end if 1012 ! Contributions from |r-R|=0 1013 if (nfgd_r0>0) then 1014 rr_tmp(1)=tol10 1015 call paw_splint(msz,core_mesh%rad,pawtab%tcoredens(:,3),pawtab%tcoredens(:,5),& 1016 & 1,rr_tmp,d2tcore(1)) 1017 ifgd=1 1018 ii=iindex(ifgd);jj=ifftsph(ii) 1019 term=vxc(jj,1)*d2tcore(ifgd) 1020 do mu=1,3 1021 do nu=1,3 1022 dyfrx2(mu,nu,iatom)=dyfrx2(mu,nu,iatom)+term*rmet(mu,nu) 1023 end do 1024 end do 1025 end if 1026 ABI_FREE(d2tcore) 1027 1028 end if !option 1029 1030 end subroutine mkcore_inner
ABINIT/mkcore_wvl [ Functions ]
NAME
mkcore_wvl
FUNCTION
Optionally compute (in a WVL representation): (1) pseudo core electron density throughout unit cell (2) pseudo-core contribution to forces (3) pseudo-core contribution to stress tensor
INPUTS
atindx1(natom)=index table for atoms, inverse of atindx
[mpi_comm_wvl]=MPI communicator (optional)
mpi_enreg=information about MPI parallelization
natom=number of atoms in cell
nattyp(ntypat)=number of atoms of each type
nfft=dimension of vxc (XC potential)
nspden=number of spin-density components
ntypat=number of types of atoms
n1xccc=dimension of xccc1d ; 0 if no XC core correction is used
n3xccc=dimension of xccc3d (pseudo core charge)
option: 1 for computing core charge density
2 for computing core charge contribution to forces
3 for computing core charge contribution to stress tensor
4 for computing contribution to frozen-wf part of dynamical matrix
pawrad(ntypat*usepaw) <type(pawrad_type)>=paw radial mesh and related data
pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
rprimd(3,3)=dimensional primitive translation vectors (bohr)
vxc(nfft,nspden)=exchange-correlation potential (hartree)
xcccrc(ntypat)=XC core correction cutoff radius (bohr) for each atom type
xccc1d(n1xccc,6,ntypat)=1D core charge function and 5 derivatives for each atom type
xred(3,natom)=reduced coordinates for atoms in unit cell
wvl_den=density-potential BigDFT object
wvl_descr=wavelet BigDFT object
OUTPUT
=== if option==1 === xccc3d(n3xccc)=3D core electron density for XC core correction (bohr^-3) === if option==2 === grxc(3,natom)=core charge contribution to forces === if option==3 === corstr(6)=core charge contribution to stress tensor
SIDE EFFECTS
xccc3d(n3xccc)=3D core electron density for XC core correction (bohr^-3) (computed and returned when option=1, needed as input when option=3)
NOTES
Based on mkcore.F90. Adapted to WVL case.
SOURCE
102 subroutine mkcore_wvl(atindx1,corstr,grxc,natom,nattyp,nfft,nspden,ntypat,n1xccc,n3xccc,option,& 103 & pawrad,pawtab,rprimd,vxc,xccc1d,xccc3d,xcccrc,xred,wvl_den,wvl_descr,& 104 & mpi_comm_wvl) ! optional argument 105 106 #if defined HAVE_BIGDFT 107 use BigDFT_API, only : PSPCODE_PAW,ind_positions 108 #endif 109 110 !Arguments ------------------------------------ 111 !scalars 112 integer,intent(in) :: natom,nfft,nspden,ntypat,n1xccc,n3xccc,option 113 integer,intent(in),optional :: mpi_comm_wvl 114 type(wvl_denspot_type), intent(inout) :: wvl_den 115 type(wvl_internal_type), intent(in) :: wvl_descr 116 !arrays 117 integer,intent(in) :: atindx1(natom),nattyp(ntypat) 118 real(dp),intent(in) :: rprimd(3,3),xccc1d(n1xccc,6,ntypat),xcccrc(ntypat),xred(3,natom) 119 real(dp),intent(in),target :: vxc(nfft,nspden) 120 real(dp),intent(out) :: corstr(6),grxc(3,natom) 121 real(dp),intent(inout) :: xccc3d(n3xccc) 122 type(pawrad_type),intent(in) :: pawrad(:) 123 type(pawtab_type),intent(in) :: pawtab(:) 124 125 !Local variables------------------------------- 126 #if defined HAVE_BIGDFT 127 !scalars 128 integer :: iat,iatm,iatom,iex,iey,iez,ind,ioffset,ipts,isx,isy,isz,itypat 129 integer :: iwarn=0,i1,i2,i3,i3s,j1,j2,j3,jj,jpts,me_wvl,nproc_wvl 130 integer :: nbl1,nbr1,nbl2,nbr2,nbl3,nbr3,npts,npts12 131 integer :: ntot,n1,n1i,n2,n2i,n3,n3i,n3pi 132 logical :: perx,pery,perz,gox,goy,goz,USE_PAW 133 real(dp) :: aa,arg,bb,cc,cutoff,dd,delta=0,deltam1=0,delta2div6=0,diff 134 real(dp) :: hxh,hyh,hzh,range,range2,rangem1,rr,rx,ry,rz,r2,strdia 135 real(dp) :: term,ucvol,xx,yy,zz 136 character(len=500) :: msg 137 type(pawrad_type) :: core_mesh 138 !arrays 139 integer,allocatable :: indx(:),ivec(:) 140 real(dp) :: corfra(3,3),corgr(3),gmet(3,3),gprimd(3,3),rmet(3,3),tsec(2),tt(3),xcart(3) 141 real(dp),allocatable :: dtcore(:),d2tcore(:),rnorm(:),tcore(:),vecx(:),vecy(:) 142 real(dp),pointer :: vxc_eff(:) 143 #endif 144 145 !************************************************************************ 146 147 #if defined HAVE_BIGDFT 148 149 call timab(12,1,tsec) 150 151 !Make sure option is acceptable 152 if (option<0.or.option>3) then 153 write(msg,'(a,i2,a)') 'Option= ',option,' is not allowed!' 154 ABI_BUG(MSG) 155 end if 156 if(nfft/=n3xccc)then 157 write(msg,'(a)') 'nfft and n3xccc should be equal!' 158 ABI_BUG(msg) 159 end if 160 161 !MPI 162 nproc_wvl=1;if (present(mpi_comm_wvl)) nproc_wvl=xmpi_comm_size(mpi_comm_wvl) 163 me_wvl=0;if (present(mpi_comm_wvl)) me_wvl=xmpi_comm_rank(mpi_comm_wvl) 164 165 !Zero out only the appropriate array according to option 166 if (option==1) then 167 xccc3d(:)=zero 168 else if (option==2) then 169 grxc(:,:)=zero 170 else if (option==3) then 171 corfra(:,:)=zero 172 corstr(:)=zero 173 strdia=zero 174 end if 175 176 !Nothing to do if no xy plane to handle 177 n3pi=wvl_den%denspot%dpbox%n3pi 178 if (n3pi==0) return 179 180 !Show how calculation runs 181 if (me_wvl==0) then 182 if (option==1) write(std_out,'(a)',advance='no') ' Compute pseudo core density...' 183 if (option==2) write(std_out,'(a)',advance='no') ' Compute forces due to core density...' 184 if (option==3) write(std_out,'(a)',advance='no') ' Compute stresses due to core density...' 185 if (option==4) write(std_out,'(a)',advance='no') ' Compute dyn. matrix due to core density...' 186 end if 187 188 !PAW or NCPP ? 189 USE_PAW=any(wvl_descr%atoms%npspcode==PSPCODE_PAW) 190 191 !Conditions for periodicity in the three directions 192 perx=(wvl_descr%atoms%astruct%geocode /= 'F') 193 pery=(wvl_descr%atoms%astruct%geocode == 'P') 194 perz=(wvl_descr%atoms%astruct%geocode /= 'F') 195 call ext_buffers(perx,nbl1,nbr1) 196 call ext_buffers(pery,nbl2,nbr2) 197 call ext_buffers(perz,nbl3,nbr3) 198 199 if (option>=2) then 200 ! For spin-polarization, replace vxc by (1/2)*(vxc(up)+vxc(down)) 201 ! For non-collinear magnetism, replace vxc by (1/2)*(vxc^{11}+vxc^{22}) 202 if (nspden>=2) then 203 ABI_MALLOC(vxc_eff,(nfft)) 204 do jj=1,nfft 205 vxc_eff(jj)=half*(vxc(jj,1)+vxc(jj,2)) 206 end do 207 else 208 vxc_eff => vxc(1:nfft,1) 209 end if 210 end if 211 212 !Some constants 213 n1=wvl_descr%Glr%d%n1 214 n2=wvl_descr%Glr%d%n2 215 n3=wvl_descr%Glr%d%n3 216 n1i=wvl_den%denspot%dpbox%ndims(1) 217 n2i=wvl_den%denspot%dpbox%ndims(2) 218 n3i=wvl_den%denspot%dpbox%ndims(3) 219 ntot=n1i*n2i*n3i 220 ioffset=n1i*n2i*wvl_den%denspot%dpbox%i3xcsh 221 i3s=1+wvl_den%denspot%dpbox%nscatterarr(me_wvl,3)-wvl_den%denspot%dpbox%i3xcsh 222 hxh=wvl_den%denspot%dpbox%hgrids(1) 223 hyh=wvl_den%denspot%dpbox%hgrids(2) 224 hzh=wvl_den%denspot%dpbox%hgrids(3) 225 call metric(gmet,gprimd,-1,rmet,rprimd,arg) 226 ucvol=real(ntot,dp)*(hxh*hyh*hzh) 227 if (.not.USE_PAW) then 228 delta=one/(n1xccc-1) 229 deltam1=n1xccc-1 230 delta2div6=delta**2/6.0_dp 231 end if 232 233 !Loop over atom types 234 iatm=0 235 do itypat=1,ntypat 236 237 ! Set search range (density cuts off perfectly beyond range) 238 range=xcccrc(itypat);if (USE_PAW) range=pawtab(itypat)%rcore 239 range2=range**2 ; rangem1=one/range 240 241 ! Skip loop if this type has no core charge 242 if (abs(range)<1.d-16) cycle 243 244 ! PAW: create mesh for core density 245 if (USE_PAW) then 246 call pawrad_init(core_mesh,mesh_size=pawtab(itypat)%core_mesh_size,& 247 & mesh_type=pawrad(itypat)%mesh_type,& 248 & rstep=pawrad(itypat)%rstep,lstep=pawrad(itypat)%lstep) 249 end if 250 251 ! Loop over atoms of the type 252 do iat=1,nattyp(itypat) 253 iatm=iatm+1;iatom=atindx1(iatm) 254 255 if (option==2) corgr(:)=zero 256 257 ! Coordinates of the center 258 call xred2xcart(1,rprimd,xcart,xred(:,iatom)) 259 rx=xcart(1) ; ry=xcart(2) ; rz=xcart(3) 260 261 ! Range of points to explore 262 cutoff=1.1_dp*range 263 isx=floor((rx-cutoff)/hxh) 264 isy=floor((ry-cutoff)/hyh) 265 isz=floor((rz-cutoff)/hzh) 266 iex=ceiling((rx+cutoff)/hxh) 267 iey=ceiling((ry+cutoff)/hyh) 268 iez=ceiling((rz+cutoff)/hzh) 269 270 ! Allocate temporary memory 271 !npts12=1+int(((range/hh(1))*(range/hh(2))*pi)) 272 npts12=(iex-isx+1)*(iey-isy+1) 273 ABI_MALLOC(rnorm,(npts12)) 274 ABI_MALLOC(vecx,(npts12)) 275 ABI_MALLOC(vecy,(npts12)) 276 ABI_MALLOC(ivec,(npts12)) 277 ABI_MALLOC(indx,(npts12)) 278 if (option==1.or.option==3) then 279 ABI_MALLOC(tcore,(npts12)) 280 end if 281 if (option>=2) then 282 ABI_MALLOC(dtcore,(npts12)) 283 end if 284 if (option==4) then 285 ABI_MALLOC(d2tcore,(npts12)) 286 end if 287 288 ! Explore range of vectors 289 do i3=isz,iez 290 zz=real(i3,kind=dp)*hzh-rz 291 call ind_positions(perz,i3,n3,j3,goz) 292 j3=j3+nbl3+1 293 294 ! Select the vectors located around the current atom 295 ! TR: all of the following could be done inside or 296 ! outside the loops (i2,i1,i3). 297 ! Outside: the memory consumption increases. 298 ! Inside: the time of calculation increases. 299 ! Here, I choose to do it here, somewhere in the middle. 300 npts=0 301 do i2=isy,iey 302 yy=real(i2,kind=dp)*hyh-ry 303 call ind_positions(pery,i2,n2,j2,goy) 304 do i1=isx,iex 305 xx=real(i1,kind=dp)*hxh-rx 306 call ind_positions(perx,i1,n1,j1,gox) 307 r2=xx**2+yy**2+zz**2 308 if ((j3>=i3s.and.j3<=i3s+n3pi-1) .and. (gox.and.goy)) then 309 ind=j1+1+nbl1+(j2+nbl2)*n1i+(j3-i3s+1-1)*n1i*n2i 310 ! Only accept contributions inside defined range 311 if (r2<range2) then 312 npts=npts+1 ; indx(npts)=npts 313 ivec(npts)=ioffset+ind 314 rnorm(npts)=sqrt(r2) 315 vecx(npts)=xx;vecy(npts)=yy 316 end if 317 end if 318 end do 319 end do 320 if (npts==0) cycle 321 if (npts>npts12) then 322 msg='npts>npts12!' 323 ABI_BUG(msg) 324 end if 325 326 ! Evaluate core density (and derivatives) on the set of selected points 327 if (USE_PAW) then 328 ! PAW: use splint routine 329 call sort_dp(npts,rnorm(1:npts),indx(1:npts),tol16) 330 if (option==1.or.option==3) then 331 ! Evaluate fit of core density 332 call paw_splint(core_mesh%mesh_size,core_mesh%rad, & 333 & pawtab(itypat)%tcoredens(:,1), & 334 & pawtab(itypat)%tcoredens(:,3),& 335 & npts,rnorm(1:npts),tcore(1:npts)) 336 end if 337 if (option>=2) then 338 ! Evaluate fit of 1-der of core density 339 call paw_splint(core_mesh%mesh_size,core_mesh%rad, & 340 & pawtab(itypat)%tcoredens(:,2), & 341 & pawtab(itypat)%tcoredens(:,4),& 342 & npts,rnorm(1:npts),dtcore(1:npts)) 343 end if 344 if (option==4) then 345 ! Evaluate fit of 2nd-der of core density 346 call paw_splint(core_mesh%mesh_size,core_mesh%rad, & 347 & pawtab(itypat)%tcoredens(:,3), & 348 & pawtab(itypat)%tcoredens(:,5),& 349 & npts,rnorm(1:npts),d2tcore(1:npts)) 350 end if 351 else 352 ! Norm-conserving PP: 353 ! Evaluate spline fit with method from Numerical Recipes 354 do ipts=1,npts 355 rr=rnorm(ipts)*rangem1 356 jj=1+int(rr*(n1xccc-1)) 357 diff=rr-(jj-1)*delta 358 bb = diff*deltam1 ; aa = one-bb 359 cc = aa*(aa**2-one)*delta2div6 360 dd = bb*(bb**2-one)*delta2div6 361 if (option==1.or.option==3) then 362 tcore(ipts)=aa*xccc1d(jj,1,itypat)+bb*xccc1d(jj+1,1,itypat) +& 363 & cc*xccc1d(jj,3,itypat)+dd*xccc1d(jj+1,3,itypat) 364 end if 365 if (option>=2) then 366 dtcore(ipts)=aa*xccc1d(jj,2,itypat)+bb*xccc1d(jj+1,2,itypat) +& 367 & cc*xccc1d(jj,4,itypat)+dd*xccc1d(jj+1,4,itypat) 368 end if 369 if (option==4) then 370 d2tcore(ipts)=aa*xccc1d(jj,3,itypat)+bb*xccc1d(jj+1,3,itypat) +& 371 & cc*xccc1d(jj,5,itypat)+dd*xccc1d(jj+1,5,itypat) 372 end if 373 end do 374 end if 375 376 ! Now, perform the loop over selected grid points 377 do ipts=1,npts 378 rr=rnorm(ipts) 379 xx=vecx(indx(ipts)) 380 yy=vecy(indx(ipts)) 381 jpts=ivec(indx(ipts)) 382 383 ! === Evaluate charge density 384 if (option==1) then 385 xccc3d(jpts)=xccc3d(jpts)+tcore(ipts) 386 387 ! === Accumulate contributions to forces 388 else if (option==2) then 389 if (rr>tol10) then 390 term=vxc_eff(jpts)*dtcore(ipts)/rr 391 corgr(1)=corgr(1)+xx*term 392 corgr(2)=corgr(2)+yy*term 393 corgr(3)=corgr(3)+yy*term 394 end if 395 396 ! === Accumulate contributions to stress tensor (in red. coordinates) 397 else if (option==3) then 398 if (rr>tol10) then 399 term=vxc_eff(jpts)*dtcore(ipts)*rangem1/rr/real(ntot,dp) 400 ! Write out the 6 symmetric components 401 corfra(1,1)=corfra(1,1)+term*xx*xx 402 corfra(2,2)=corfra(2,2)+term*yy*yy 403 corfra(3,3)=corfra(3,3)+term*zz*zz 404 corfra(3,2)=corfra(3,2)+term*zz*yy 405 corfra(3,1)=corfra(3,1)+term*zz*xx 406 corfra(2,1)=corfra(2,1)+term*yy*xx 407 ! (the above still needs to be transformed to cartesian coords) 408 end if 409 ! Also compute a diagonal term 410 strdia=strdia+vxc_eff(jpts)*tcore(ipts) 411 412 end if ! Choice of option 413 414 end do ! ipts (i1,i2) 415 416 end do ! i3 417 418 ! Release temporary memory 419 ABI_FREE(rnorm) 420 ABI_FREE(vecx) 421 ABI_FREE(vecy) 422 ABI_FREE(ivec) 423 ABI_FREE(indx) 424 if (allocated(tcore)) then 425 ABI_FREE(tcore) 426 end if 427 if (allocated(dtcore)) then 428 ABI_FREE(dtcore) 429 end if 430 if (allocated(d2tcore)) then 431 ABI_FREE(d2tcore) 432 end if 433 434 if (option==2) then 435 arg=-(ucvol/real(ntot,dp)) 436 !arg=-(ucvol/real(ntot,dp))/range !???? 437 grxc(:,iatom)=corgr(:)*arg 438 end if 439 440 ! End loop on atoms 441 end do 442 443 if (USE_PAW) then 444 call pawrad_free(core_mesh) 445 end if 446 447 !End loop over atom types 448 end do 449 450 if(option>=2.and.nspden>=2) then 451 ABI_FREE(vxc_eff) 452 end if 453 454 !Density: make it positive 455 if (option==1) then 456 call mkdenpos(iwarn,n3xccc,nspden,0,xccc3d,tol20) 457 end if 458 459 !Forces: translate into reduced coordinates 460 if (option==2) then 461 do iatom=1,natom 462 tt(1:3)=grxc(1:3,iatom) 463 grxc(:,iatom)= rprimd(1,:)*tt(1)+rprimd(2,:)*tt(2)+rprimd(3,:)*tt(3) 464 !grxc(:,iatom)=rmet(:,1)*tt(1)+rmet(:,2)*tt(2)+rmet(:,3)*tt(3) 465 end do 466 end if 467 468 !Stress tensor: symmetrize, translate into cartesian coord., add diagonal part 469 if (option==3) then 470 corstr(1)=corfra(1,1) ; corstr(2)=corfra(2,2) 471 corstr(3)=corfra(3,3) ; corstr(4)=corfra(3,2) 472 corstr(5)=corfra(3,1) ; corstr(6)=corfra(2,1) 473 call strconv(corstr,rprimd,corstr) 474 corstr(1)=corstr(1)+strdia/real(ntot,dp) 475 corstr(2)=corstr(2)+strdia/real(ntot,dp) 476 corstr(3)=corstr(3)+strdia/real(ntot,dp) 477 end if 478 479 !If needed, sum over MPI processes 480 if(nproc_wvl>1) then 481 call timab(539,1,tsec) 482 if (option==2) then 483 call xmpi_sum(grxc,mpi_comm_wvl,iex) 484 end if 485 if (option==3) then 486 call xmpi_sum(corstr,mpi_comm_wvl,iex) 487 end if 488 call timab(539,2,tsec) 489 end if 490 491 if (me_wvl==0) write(std_out,'(a)') 'done.' 492 493 call timab(12,1,tsec) 494 495 #else 496 BIGDFT_NOTENABLED_ERROR() 497 ABI_UNUSED(xcccrc) 498 if (.false.) write(std_out,*) natom,nfft,nspden,ntypat,n1xccc,n3xccc,option,mpi_comm_wvl,& 499 & wvl_den%symObj,wvl_descr%h(1),atindx1(1),nattyp(1),rprimd(1,1),vxc(1,1),& 500 & xred(1,1),xccc1d(1,1,1),corstr(1),grxc(1,1),xccc3d(1),pawrad(1)%mesh_size,pawtab(1)%lmn_size 501 #endif 502 503 end subroutine mkcore_wvl
ABINIT/mkcore_wvl_old [ Functions ]
NAME
mkcore_wvl_old
FUNCTION
Optionally compute (1) core electron density throughout unit cell, or (2) contribution to Exc gradient wrt xred, or (3) contribution to stress tensor, or (response function code) (4) contribution to frozen-wavefunction part of the dynamical matrix (part 2)
INPUTS
argin(sizein)=description
OUTPUT
argout(sizeout)=description
SIDE EFFECTS
NOTES
Based on mkcore.F90
SOURCE
534 subroutine mkcore_wvl_old(atindx1,corstr,dyfrx2,geocode,grxc,h,natom,& 535 & nattyp,nfft,nscatterarr,nspden,ntypat,n1,n1i,n2,n2i,n3,n3pi,& 536 & n3xccc,option,pawrad,pawtab,psppar,rprimd,ucvol,& 537 & vxc,xccc3d,xred,mpi_comm_wvl) 538 539 #if defined HAVE_BIGDFT 540 use BigDFT_API, only: ext_buffers,ind_positions 541 #endif 542 543 !Arguments --------------------------------------------- 544 !scalars 545 integer,intent(in) :: natom,ntypat,nfft,nspden 546 integer,intent(in) ::n1,n2,n3,n1i,n2i,n3pi,n3xccc,option 547 integer,intent(in),optional:: mpi_comm_wvl 548 real(dp),intent(in) :: h(3),ucvol 549 character(1),intent(in)::geocode 550 !arrays 551 integer,intent(in) :: atindx1(natom),nattyp(ntypat) 552 integer,intent(in) :: nscatterarr(4) 553 real(dp),intent(in) :: psppar(0:4,0:6,ntypat),rprimd(3,3) 554 real(dp),intent(in)::xred(3,natom) 555 real(dp),intent(in)::vxc(nfft,nspden) 556 real(dp),intent(out)::xccc3d(n3xccc) 557 real(dp),intent(out) :: corstr(6),dyfrx2(3,3,natom),grxc(3,natom) 558 type(pawtab_type),intent(in) :: pawtab(ntypat) 559 type(pawrad_type),intent(in) :: pawrad(ntypat) 560 561 !Local variables ------------------------------ 562 #if defined HAVE_BIGDFT 563 !scalars 564 !buffer to be added at the end of the last dimension of an array to control bounds_check 565 integer :: i1,i2,i3,iat,iatm,iatom,iatom_tot 566 integer :: itypat,iwarn 567 integer :: nproc_wvl=1 568 integer :: ier,iex,iey,iez,isx,isy,isz,ind,i3s 569 integer :: j1,j2,j3,msz 570 integer :: nbl1,nbr1,nbl2,nbr2,nbl3,nbr3 571 integer :: ncmax,nfgd,nfgd_r0,opt_mkdenpos,shift 572 real(dp) :: cutoff,factor,grxc1,grxc2,grxc3 573 real(dp) :: rloc,rr2,rx,ry,rz 574 real(dp) :: rshp,r2shp,rshpm1,strdia,t1,t2,t3 575 real(dp) :: xx,yy,zz,ucvol_ 576 character(len=500) :: message 577 logical :: perx,pery,perz,gox,goy,goz 578 type(pawrad_type)::core_mesh 579 !arrays 580 real(dp) :: hh(3) !fine grid spacing for wavelets 581 real(dp) :: gmet(3,3),gprimd(3,3),rcart(3),rmet(3,3),tsec(2),xcart(3,natom) 582 real(dp) :: corfra(3,3) 583 !allocatable arrays 584 integer,allocatable :: ifftsph_tmp(:) 585 real(dp),allocatable:: rr(:),rred(:,:) 586 #endif 587 588 ! ************************************************************************* 589 590 DBG_ENTER("COLL") 591 592 #if defined HAVE_BIGDFT 593 594 if(nspden >1) then 595 write(message, '(a)')'mkcore_wvl_old: this is not yet generalized to npsden>1' 596 ABI_ERROR(message) 597 end if 598 if(option>4 .or. option<1 )then 599 write(message,'(a,a,a,a,a,a,i6)') ch10,& 600 & ' mkcore_wvl_old: BUG -',ch10,& 601 & ' The argument option should be between 1 and 4,',ch10,& 602 & ' however, option=',option 603 ABI_BUG(message) 604 end if 605 if(nfft .ne. n3xccc)then 606 write(message,'(a,a,a,a,a,a,2i6)') ch10,& 607 & ' mkcore_wvl_old: BUG -',ch10,& 608 & ' nfft and n3xccc should be equal,',ch10,& 609 & ' however, nfft and n3xccc=',nfft,n3xccc 610 ABI_BUG(message) 611 end if 612 613 !MPI 614 if(present(mpi_comm_wvl)) nproc_wvl=xmpi_comm_size(mpi_comm_wvl) 615 i3s =nscatterarr(3)+1-nscatterarr(4) 616 shift=n1i*n2i*nscatterarr(4) 617 618 if (option==1) then 619 ! Zero out array to permit accumulation over atom types below: 620 xccc3d(:)=zero 621 else if (option==2) then 622 ! Zero out gradient of Exc array 623 grxc(:,:)=zero 624 else if (option==3) then 625 ! Zero out locally defined stress array 626 corfra(:,:)=zero 627 strdia=zero 628 else if (option==4) then 629 ! Zero out fr-wf part of the dynamical matrix 630 dyfrx2(:,:,:)=zero 631 else 632 write(message, '(a,a,a,a)' ) ch10,& 633 & ' mkcore_wvl_old: BUG -',ch10,& 634 & ' Can''t be here ! (bad option).' 635 ABI_BUG(message) 636 end if 637 638 write(message,'(a,a)') ch10,& 639 & ' mkcore_wvl_old: Compute core density' 640 call wrtout(std_out,message,'COLL') 641 642 !Fine grid 643 hh(:)=0.5d0*h(:) 644 645 !Compute xcart from xred 646 call xred2xcart(natom,rprimd,xcart,xred) 647 648 !Compute metric tensors and ucvol from rprimd 649 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol_) 650 !correct ucvol for wavelets case is given as an input: 651 !ucvol_local = product(wvl%den%denspot%dpbox%hgrids) * real(product(wvl%den%denspot%dpbox%ndims), dp) 652 !ucvol_local = (half * dtset%wvl_hgrid) ** 3 * ngfft(1)*ngfft(2)*ngfft(3) 653 654 !Conditions for periodicity in the three directions 655 perx=(geocode /= 'F') 656 pery=(geocode == 'P') 657 perz=(geocode /= 'F') 658 659 !Compute values of external buffers 660 call ext_buffers(perx,nbl1,nbr1) 661 call ext_buffers(pery,nbl2,nbr2) 662 call ext_buffers(perz,nbl3,nbr3) 663 664 665 iatm=0 666 !Big loop on atom types 667 do itypat=1,ntypat 668 669 rloc=psppar(0,0,itypat) 670 cutoff=10.d0*rloc 671 672 ! Set radius size: 673 rshp=pawtab(itypat)%rcore 674 r2shp=1.0000001_dp*rshp**2 675 rshpm1=one/rshp 676 677 ! allocate arrays 678 if (n3pi > 0) then 679 ! ncmax=1+int(1.1_dp*nfft*four_pi/(three*ucvol)*rshp**3) 680 ! ncmax=1+int(1.1_dp*nfft*four_pi/(three*ucvol)*rshp**3) 681 ! 1+int(1.1* factors are included just for cautioness 682 ncmax=1+int(1.1d0*((rshp/hh(1))*(rshp/hh(2))*pi)) 683 else 684 ncmax=1 685 end if 686 687 ABI_MALLOC(ifftsph_tmp,(ncmax)) 688 ABI_MALLOC(rr,(ncmax)) 689 if(option>1) then 690 ABI_MALLOC(rred,(3,ncmax)) 691 else 692 ABI_MALLOC(rred,(0,0)) 693 end if 694 695 ! Create mesh_core object 696 ! since core_mesh_size can be bigger than pawrad%mesh_size, 697 msz=pawtab(itypat)%core_mesh_size 698 call pawrad_init(core_mesh,mesh_size=msz,mesh_type=pawrad(itypat)%mesh_type,& 699 & rstep=pawrad(itypat)%rstep,lstep=pawrad(itypat)%lstep) 700 701 ! Big loop on atoms 702 do iat=1,nattyp(itypat) 703 iatm=iatm+1;iatom=atindx1(iatm) 704 iatom_tot=iatom 705 706 if(option==2) then 707 grxc1=zero 708 grxc2=zero 709 grxc3=zero 710 end if 711 712 ! Define a "box" around each atom 713 rx=xcart(1,iatom_tot) 714 ry=xcart(2,iatom_tot) 715 rz=xcart(3,iatom_tot) 716 717 isx=floor((rx-cutoff)/hh(1)) 718 isy=floor((ry-cutoff)/hh(2)) 719 isz=floor((rz-cutoff)/hh(3)) 720 721 iex=ceiling((rx+cutoff)/hh(1)) 722 iey=ceiling((ry+cutoff)/hh(2)) 723 iez=ceiling((rz+cutoff)/hh(3)) 724 725 do i3=isz,iez 726 zz=real(i3,kind=8)*hh(3)-rz 727 call ind_positions(perz,i3,n3,j3,goz) 728 j3=j3+nbl3+1 729 730 ! Initialize counters 731 nfgd=0 732 nfgd_r0=0 733 734 do i2=isy,iey 735 yy=real(i2,kind=8)*hh(2)-ry 736 call ind_positions(pery,i2,n2,j2,goy) 737 738 do i1=isx,iex 739 xx=real(i1,kind=8)*hh(1)-rx 740 call ind_positions(perx,i1,n1,j1,gox) 741 rr2=xx**2+yy**2+zz**2 742 if (j3 >= i3s .and. j3 <= i3s+n3pi-1 .and. goy .and. gox ) then 743 744 if(rr2<=r2shp) then 745 if(rr2>tol5) then 746 ind=j1+1+nbl1+(j2+nbl2)*n1i+(j3-i3s)*n1i*n2i 747 nfgd=nfgd+1 748 rcart(1)=xx; rcart(2)=yy; rcart(3)=zz 749 rr(nfgd)=(rr2)**0.5 750 ifftsph_tmp(nfgd)=shift+ind 751 if(option>1) then 752 call xcart2xred(1,rprimd,rcart,rred(:,nfgd)) 753 end if 754 else if (option==4) then 755 ! We save r=0 vectors only for option==4: 756 ! for other options this is ignored 757 ind=j1+1+nbl1+(j2+nbl2)*n1i+(j3-i3s)*n1i*n2i 758 ! We reuse the same variable "ifftshp_tmp", 759 ! but we start from the higher index 760 nfgd_r0=nfgd_r0+1 761 ifftsph_tmp(ncmax-nfgd_r0+1)=shift+ind 762 763 end if !rr2>tol5 764 end if !rr2<r2shp 765 end if !j3.. 766 end do !i1 767 end do !i2 768 769 ! All of the following could be done inside or outside the loops (i2,i1,i3) 770 ! Outside the loops: the memory consuption increases. 771 ! Inside the inner loop: the time of calculation increases. 772 ! Here, I chose to do it here, somewhere in the middle. 773 774 if(option .ne.4 ) then 775 if(nfgd==0) cycle 776 else 777 if(nfgd==0 .and. nfgd_r0==0) cycle 778 end if 779 call mkcore_inner(corfra,core_mesh,dyfrx2,& 780 & grxc1,grxc2,grxc3,ifftsph_tmp,msz,& 781 & natom,ncmax,nfft,nfgd,nfgd_r0,nspden,n3xccc,option,pawtab(itypat),& 782 & rmet,rr,strdia,vxc,xccc3d,rred=rred) 783 784 end do !i3 785 786 if(option==2) then 787 factor=(ucvol/real(nfft,dp))/rshp 788 grxc(1,iatom)=grxc1*factor 789 grxc(2,iatom)=grxc2*factor 790 grxc(3,iatom)=grxc3*factor 791 792 if( nproc_wvl>1 ) then 793 call timab(539,1,tsec) 794 call xmpi_sum(grxc1,mpi_comm_wvl,ier) 795 call xmpi_sum(grxc2,mpi_comm_wvl,ier) 796 call xmpi_sum(grxc3,mpi_comm_wvl,ier) 797 call timab(539,2,tsec) 798 end if 799 end if 800 801 end do !iat 802 803 ! Deallocate 804 call pawrad_free(core_mesh) 805 ABI_FREE(ifftsph_tmp) 806 ABI_FREE(rr) 807 ABI_FREE(rred) 808 809 end do !itypat 810 811 if (option==2) then 812 ! Apply rmet as needed to get reduced coordinate gradients 813 do iatom=1,natom 814 t1=grxc(1,iatom) 815 t2=grxc(2,iatom) 816 t3=grxc(3,iatom) 817 grxc(:,iatom)=rmet(:,1)*t1+rmet(:,2)*t2+rmet(:,3)*t3 818 end do 819 820 elseif (option==3) then 821 822 ! Transform stress tensor from full storage mode to symmetric storage mode 823 corstr(1)=corfra(1,1) 824 corstr(2)=corfra(2,2) 825 corstr(3)=corfra(3,3) 826 corstr(4)=corfra(3,2) 827 corstr(5)=corfra(3,1) 828 corstr(6)=corfra(2,1) 829 830 ! Transform stress tensor from reduced coordinates to cartesian coordinates 831 call strconv(corstr,rprimd,corstr) 832 833 ! Compute diagonal contribution to stress tensor (need input xccc3d) 834 ! strdia = (1/N) Sum(r) [mu_xc_avg(r) * rho_core(r)] 835 strdia=zero 836 do i3=1,n3pi 837 do i2=1,n2i 838 do i1=1,n1i 839 ind=i1+(i2-1)*n1i+(i3-1)*n1i*n2i 840 strdia=strdia+vxc(ind,1)*xccc3d(ind) 841 ! write(17,'(3(i6),i12,3(1x,1pe24.17))')i1,i2,i3,ind,potion_corr(ind),pot_ion(ind),maxdiff 842 end do 843 end do 844 end do 845 strdia=strdia/real(nfft,dp) 846 847 ! Add diagonal term to stress tensor 848 corstr(1)=corstr(1)+strdia 849 corstr(2)=corstr(2)+strdia 850 corstr(3)=corstr(3)+strdia 851 end if 852 853 if(nproc_wvl > 1) then 854 call timab(539,1,tsec) 855 if(option==3) then 856 call xmpi_sum(corstr,mpi_comm_wvl,ier) 857 end if 858 if(option==2) then 859 call xmpi_sum(grxc,mpi_comm_wvl,ier) 860 end if 861 call timab(539,2,tsec) 862 end if 863 864 !Make xccc3d positive to avoid numerical instabilities in V_xc 865 iwarn=0 ; opt_mkdenpos=0 866 call mkdenpos(iwarn, n3xccc, nspden, opt_mkdenpos, xccc3d, tol20 ) 867 868 #else 869 BIGDFT_NOTENABLED_ERROR() 870 if (.false.) write(std_out,*) natom,ntypat,nfft,nspden,n1,n2,n3,n1i,n2i,n3pi,n3xccc,option,& 871 & mpi_comm_wvl,h(1),ucvol,geocode,atindx1(1),nattyp(1),nscatterarr(1),psppar(1,1,1),rprimd(1,1),& 872 & xred(1,1),vxc(1,1),xccc3d(1),corstr(1),dyfrx2(1,1,1),grxc(1,1),pawtab(1)%mesh_size,pawrad(1)%mesh_size 873 #endif 874 875 DBG_EXIT("COLL") 876 877 end subroutine mkcore_wvl_old