TABLE OF CONTENTS


ABINIT/m_mkcore_wvl [ Modules ]

[ Top ] [ 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 ]

[ Top ] [ 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 ]

[ Top ] [ 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 ]

[ Top ] [ 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