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