TABLE OF CONTENTS


ABINIT/dfpt_vlocal [ Functions ]

[ Top ] [ Functions ]

NAME

 dfpt_vlocal

FUNCTION

 Compute local part of 1st-order potential from the appropriate
 atomic pseudopotential with structure and derivative factor.
 In case of derivative with respect to k or
 electric (magnetic Zeeman) field perturbation, the 1st-order local potential vanishes.

INPUTS

  atindx(natom)=index table for atoms (see gstate.f)
  cplex: if 1, real space 1-order functions on FFT grid
    are REAL, if 2, COMPLEX
  gmet(3,3)=reciprocal space metric (Bohr**-2)
  gsqcut=cutoff G**2 for included G s in fft box.
  idir=direction of atomic displacement (=1,2 or 3 : displacement of
    atom ipert along the 1st, 2nd or 3rd axis).
  ipert=number of the atom being displaced in the frozen-phonon
  mpi_enreg=information about MPI parallelization
  mqgrid=dimension of q grid for pseudopotentials
  natom=number of atoms in cell.
  nattyp(ntypat)=number of atoms of each type in cell.
  nfft=(effective) number of FFT grid points (for this processor)
  ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
  ntypat=number of types of atoms in cell.
  n1,n2,n3=fft grid.
  ph1d(2,3*(2*mgfft+1)*natom)=1-dim structure factor phase information.
  qgrid(mqgrid)=grid of q points from 0 to qmax.
  qphon(3)=wavevector of the phonon
  ucvol=unit cell volume (Bohr**3).
  vlspl(mqgrid,2,ntypat)=spline fit of q^2 V(q) for each type of atom.
  xred(3,natom)=reduced atomic coordinates

OUTPUT

  vpsp1(cplex*nfft)=first-order local crystal pseudopotential in real space
    (including the minus sign, forgotten in the paper non-linear..

PARENTS

      dfpt_looppert,dfpt_nstdy,dfpt_nstpaw,dfptnl_loop

CHILDREN

      fourdp,ptabs_fourdp

SOURCE

 815 subroutine dfpt_vlocal(atindx,cplex,gmet,gsqcut,idir,ipert,&
 816 & mpi_enreg,mqgrid,natom,nattyp,nfft,ngfft,&
 817 & ntypat,n1,n2,n3,paral_kgb,ph1d,qgrid,qphon,ucvol,vlspl,vpsp1,xred)
 818 
 819 
 820 !This section has been created automatically by the script Abilint (TD).
 821 !Do not modify the following lines by hand.
 822 #undef ABI_FUNC
 823 #define ABI_FUNC 'dfpt_vlocal'
 824 !End of the abilint section
 825 
 826  implicit none
 827 
 828 !Arguments -------------------------------
 829 !scalars
 830  integer,intent(in) :: cplex,idir,ipert,mqgrid,n1,n2,n3,natom,nfft,ntypat
 831  integer,intent(in) :: paral_kgb
 832  real(dp),intent(in) :: gsqcut,ucvol
 833  type(MPI_type),intent(in) :: mpi_enreg
 834 !arrays
 835  integer,intent(in) :: atindx(natom),nattyp(ntypat),ngfft(18)
 836  real(dp),intent(in) :: gmet(3,3),ph1d(2,(2*n1+1+2*n2+1+2*n3+1)*natom)
 837  real(dp),intent(in) :: qgrid(mqgrid),qphon(3),vlspl(mqgrid,2,ntypat)
 838  real(dp),intent(in) :: xred(3,natom)
 839  real(dp),intent(out) :: vpsp1(cplex*nfft)
 840 
 841 !Local variables -------------------------
 842 !scalars
 843  integer :: i1,i2,i3,ia1,iatom,id1,id2,id3,ig1,ig2,ig3,ii,ii1,im=2
 844  integer :: itypat,jj,re=1
 845  real(dp),parameter :: tolfix=1.000000001_dp
 846  real(dp) :: aa,bb,cc,cutoff,dd,diff,dq,dq2div6,dqdiv6,dqm1,gmag,gq1
 847  real(dp) :: gq2,gq3,gsquar,phqim,phqre
 848  real(dp) :: qxred2pi,sfi,sfr,vion1,xnorm
 849  logical :: qeq0
 850 !arrays
 851  integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:)
 852  integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:)
 853  real(dp) :: gq(3)
 854  real(dp),allocatable :: work1(:,:)
 855 
 856 ! *********************************************************************
 857 
 858  iatom=ipert
 859 
 860  if(iatom==natom+1 .or. iatom==natom+2 .or. iatom==natom+10  .or. iatom==natom+11 .or. iatom==natom+5)then
 861 
 862 !  (In case of d/dk or an electric field, or magnetic (Zeeman) field->[natom+5] SPr deb )
 863    vpsp1(1:cplex*nfft)=zero
 864 
 865  else
 866 
 867 !  (In case of a phonon perturbation)
 868    ABI_ALLOCATE(work1,(2,nfft))
 869    work1(1:2,1:nfft)=0.0_dp
 870 
 871    dq=(qgrid(mqgrid)-qgrid(1))/dble(mqgrid-1)
 872    dqm1=1.0_dp/dq
 873    dqdiv6=dq/6.0_dp
 874    dq2div6=dq**2/6.0_dp
 875    cutoff=gsqcut*tolfix
 876    id1=n1/2+2
 877    id2=n2/2+2
 878    id3=n3/2+2
 879 
 880    ! Get the distrib associated with this fft_grid
 881    call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
 882 
 883 !  This is to allow q=0
 884    qeq0=.false.
 885    if(qphon(1)**2+qphon(2)**2+qphon(3)**2<1.d-15)qeq0=.true.
 886 
 887 !  Determination of the atom type
 888    ia1=0
 889    itypat=0
 890    do ii=1,ntypat
 891      ia1=ia1+nattyp(ii)
 892      if(atindx(iatom)<=ia1.and.itypat==0)itypat=ii
 893    end do
 894 
 895 !  Determination of phase qxred*
 896    qxred2pi=2.0_dp*pi*(qphon(1)*xred(1,iatom)+ &
 897 &   qphon(2)*xred(2,iatom)+ &
 898 &   qphon(3)*xred(3,iatom) )
 899    phqre=cos(qxred2pi)
 900    phqim=sin(qxred2pi)
 901    ii=0
 902 
 903    do i3=1,n3
 904      ig3=i3-(i3/id3)*n3-1
 905      gq3=dble(ig3)+qphon(3)
 906      gq(3)=gq3
 907      do i2=1,n2
 908        if (fftn2_distrib(i2)==mpi_enreg%me_fft) then
 909          ig2=i2-(i2/id2)*n2-1
 910          gq2=dble(ig2)+qphon(2)
 911          gq(2)=gq2
 912 
 913 !        Note the lower limit of the next loop
 914          ii1=1
 915          if(i3==1 .and. i2==1 .and. qeq0 .and. ig2==0 .and. ig3==0)then
 916            ii1=2
 917            ii=ii+1
 918          end if
 919          do i1=ii1,n1
 920            ig1=i1-(i1/id1)*n1-1
 921            gq1=dble(ig1)+qphon(1)
 922            gq(1)=gq1
 923            ii=ii+1
 924            gsquar=gsq_vl3(gq1,gq2,gq3)
 925 !          Skip G**2 outside cutoff:
 926            if (gsquar<=cutoff) then
 927              gmag=sqrt(gsquar)
 928 
 929 !            Compute vion(G) for given type of atom
 930              jj=1+int(gmag*dqm1)
 931              diff=gmag-qgrid(jj)
 932 
 933 !            Evaluate spline fit from q^2 V(q) to get V(q):
 934 !            (p. 86 Numerical Recipes, Press et al; NOTE error in book for sign
 935 !            of "aa" term in derivative; also see splfit routine.
 936 !            This bug fixed here 27 Jan 1992.)
 937 
 938              bb = diff*dqm1
 939              aa = 1.0_dp-bb
 940              cc = aa*(aa**2-1.0_dp)*dq2div6
 941              dd = bb*(bb**2-1.0_dp)*dq2div6
 942              vion1 = (aa*vlspl(jj,1,itypat)+bb*vlspl(jj+1,1,itypat) + &
 943 &             cc*vlspl(jj,2,itypat)+dd*vlspl(jj+1,2,itypat) ) &
 944 &             / gsquar
 945 
 946 !            Phase   G*xred  (complex conjugate) * -i *2pi*(g+q)*vion
 947              sfr=-phimag_vl3(ig1,ig2,ig3,iatom)*2.0_dp*pi*gq(idir)*vion1
 948              sfi=-phre_vl3(ig1,ig2,ig3,iatom)*2.0_dp*pi*gq(idir)*vion1
 949 
 950 !            Phase   q*xred  (complex conjugate)
 951              work1(re,ii)=sfr*phqre+sfi*phqim
 952              work1(im,ii)=-sfr*phqim+sfi*phqre
 953            end if
 954 
 955          end do
 956        end if
 957      end do
 958    end do
 959 
 960 !  Transform back to real space
 961    call fourdp(cplex,work1,vpsp1,1,mpi_enreg,nfft,ngfft,paral_kgb,0)
 962 
 963    xnorm=1.0_dp/ucvol
 964    vpsp1(1:cplex*nfft)=vpsp1(1:cplex*nfft)*xnorm
 965 
 966    ABI_DEALLOCATE(work1)
 967 
 968 !  End the condition of non-electric-field
 969  end if
 970 
 971  contains
 972 
 973 !Real and imaginary parts of phase.
 974    function phr_vl3(x1,y1,x2,y2,x3,y3)
 975 
 976 
 977 !This section has been created automatically by the script Abilint (TD).
 978 !Do not modify the following lines by hand.
 979 #undef ABI_FUNC
 980 #define ABI_FUNC 'phr_vl3'
 981 !End of the abilint section
 982 
 983    real(dp) :: phr_vl3
 984    real(dp),intent(in) :: x1,x2,x3,y1,y2,y3
 985    phr_vl3=(x1*x2-y1*y2)*x3-(y1*x2+x1*y2)*y3
 986  end function phr_vl3
 987 
 988    function phi_vl3(x1,y1,x2,y2,x3,y3)
 989 
 990 
 991 !This section has been created automatically by the script Abilint (TD).
 992 !Do not modify the following lines by hand.
 993 #undef ABI_FUNC
 994 #define ABI_FUNC 'phi_vl3'
 995 !End of the abilint section
 996 
 997    real(dp) :: phi_vl3
 998    real(dp),intent(in) :: x1,x2,x3,y1,y2,y3
 999    phi_vl3=(x1*x2-y1*y2)*y3+(y1*x2+x1*y2)*x3
1000  end function phi_vl3
1001 
1002 !  Warning : this function differ from similar ones for ground-state calculations : note the atindx !!
1003    function ph1_vl3(nri,ig1,ia)
1004 
1005 
1006 !This section has been created automatically by the script Abilint (TD).
1007 !Do not modify the following lines by hand.
1008 #undef ABI_FUNC
1009 #define ABI_FUNC 'ph1_vl3'
1010 !End of the abilint section
1011 
1012    real(dp) :: ph1_vl3
1013    integer,intent(in) :: nri,ig1,ia
1014    ph1_vl3=ph1d(nri,ig1+1+n1+(atindx(ia)-1)*(2*n1+1))
1015  end function ph1_vl3
1016 
1017 !  Warning : this function differ from similar ones for ground-state calculations : note the atindx !!
1018    function ph2_vl3(nri,ig2,ia)
1019 
1020 
1021 !This section has been created automatically by the script Abilint (TD).
1022 !Do not modify the following lines by hand.
1023 #undef ABI_FUNC
1024 #define ABI_FUNC 'ph2_vl3'
1025 !End of the abilint section
1026 
1027    real(dp) :: ph2_vl3
1028    integer,intent(in) :: nri,ig2,ia
1029    ph2_vl3=ph1d(nri,ig2+1+n2+(atindx(ia)-1)*(2*n2+1)+natom*(2*n1+1))
1030  end function ph2_vl3
1031 
1032 !  Warning : this function differ from similar ones for ground-state calculations : note the atindx !!
1033    function ph3_vl3(nri,ig3,ia)
1034 
1035 
1036 !This section has been created automatically by the script Abilint (TD).
1037 !Do not modify the following lines by hand.
1038 #undef ABI_FUNC
1039 #define ABI_FUNC 'ph3_vl3'
1040 !End of the abilint section
1041 
1042    real(dp) :: ph3_vl3
1043    integer,intent(in) :: nri,ig3,ia
1044    ph3_vl3=ph1d(nri,ig3+1+n3+(atindx(ia)-1)*(2*n3+1)+natom*(2*n1+1+2*n2+1))
1045  end function ph3_vl3
1046 
1047    function phre_vl3(ig1,ig2,ig3,ia)
1048 
1049 
1050 !This section has been created automatically by the script Abilint (TD).
1051 !Do not modify the following lines by hand.
1052 #undef ABI_FUNC
1053 #define ABI_FUNC 'phre_vl3'
1054 !End of the abilint section
1055 
1056    real(dp) :: phre_vl3
1057    integer,intent(in) :: ig1,ig2,ig3,ia
1058    phre_vl3=phr_vl3(ph1_vl3(re,ig1,ia),ph1_vl3(im,ig1,ia),&
1059 &   ph2_vl3(re,ig2,ia),ph2_vl3(im,ig2,ia),ph3_vl3(re,ig3,ia),ph3_vl3(im,ig3,ia))
1060  end function phre_vl3
1061 
1062    function phimag_vl3(ig1,ig2,ig3,ia)
1063 
1064 
1065 !This section has been created automatically by the script Abilint (TD).
1066 !Do not modify the following lines by hand.
1067 #undef ABI_FUNC
1068 #define ABI_FUNC 'phimag_vl3'
1069 !End of the abilint section
1070 
1071    real(dp) :: phimag_vl3
1072    integer,intent(in) :: ig1,ig2,ig3,ia
1073    phimag_vl3=phi_vl3(ph1_vl3(re,ig1,ia),ph1_vl3(im,ig1,ia),&
1074 &   ph2_vl3(re,ig2,ia),ph2_vl3(im,ig2,ia),ph3_vl3(re,ig3,ia),ph3_vl3(im,ig3,ia))
1075  end function phimag_vl3
1076 
1077    function gsq_vl3(g1,g2,g3)
1078 
1079 
1080 !This section has been created automatically by the script Abilint (TD).
1081 !Do not modify the following lines by hand.
1082 #undef ABI_FUNC
1083 #define ABI_FUNC 'gsq_vl3'
1084 !End of the abilint section
1085 
1086    real(dp) :: gsq_vl3
1087    real(dp),intent(in) :: g1,g2,g3 ! Note that they are real, unlike in other similar function definitions
1088 !Define G^2 based on G space metric gmet.
1089    gsq_vl3=g1*g1*gmet(1,1)+g2*g2*gmet(2,2)+&
1090 &   g3*g3*gmet(3,3)+2.0_dp*g1*g2*gmet(1,2)+&
1091 &   2.0_dp*g2*g3*gmet(2,3)+2.0_dp*g3*g1*gmet(3,1)
1092  end function gsq_vl3
1093 
1094 end subroutine dfpt_vlocal

ABINIT/m_mklocl [ Modules ]

[ Top ] [ Modules ]

NAME

  m_mklocl

FUNCTION

   Routines related to the local part of the pseudopotentials.

COPYRIGHT

  Copyright (C) 1998-2018 ABINIT group (DCA, XG, GMR, MM, DRH)
  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_mklocl
28 
29  use defs_basis
30  use defs_datatypes
31  use defs_abitypes
32  use defs_wvltypes
33  use m_abicore
34  use m_errors
35  use m_xmpi
36 
37  use m_time,     only : timab
38  use m_geometry, only : xred2xcart
39  use m_mpinfo,   only : ptabs_fourdp
40  use m_pawtab,   only : pawtab_type
41  use m_mklocl_realspace, only : mklocl_realspace, mklocl_wavelets
42  use m_fft,      only : fourdp
43 
44 #if defined HAVE_BIGDFT
45  use BigDFT_API, only : ELECTRONIC_DENSITY
46  use m_abi2big, only : wvl_rho_abi2big
47 #endif
48 
49  implicit none
50 
51  private

ABINIT/mklocl [ Functions ]

[ Top ] [ Functions ]

NAME

 mklocl

FUNCTION

 This method is a wrapper for mklocl_recipspace and mklocl_realspace.
 It does some consistency checks before calling one of the two methods.

 Optionally compute :
  option=1 : local ionic potential throughout unit cell
  option=2 : contribution of local ionic potential to E gradient wrt xred
  option=3 : contribution of local ionic potential to
                stress tensor (only with reciprocal space computations)
  option=4 : contribution of local ionic potential to
                second derivative of E wrt xred  (only with reciprocal space computations)

INPUTS

  if(option==3) eei=local pseudopotential part of total energy (hartree)
  gmet(3,3)=reciprocal space metric ($\textrm{Bohr}^{-2}$).
  gprimd(3,3)=reciprocal space dimensional primitive translations
  gsqcut=cutoff on $|G|^2$: see setup1 for definition (doubled sphere).
  mgfft=maximum size of 1D FFTs
  mpi_enreg=information about MPI parallelization
  natom=number of atoms in unit cell.
  nattyp(ntypat)=number of atoms of each type in cell.
  nfft=(effective) number of FFT grid points (for this processor)
  ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
  nspden=number of spin-density components
  ntypat=number of types of atoms.
  option= (see above)
  pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data
  ph1d(2,3*(2*mgfft+1)*natom)=1-dim structure factor phase information.
  psps <type(pseudopotential_type)>=variables related to pseudopotentials
  qprtrb(3)= integer wavevector of possible perturbing potential
   in basis of reciprocal lattice translations
  rhog(2,nfft)=electron density rho(G) (electrons/$\textrm{Bohr}^3$)
    (needed if option==2 or if option==4)
  rhor(nfft,nspden)=electron density in electrons/bohr**3.
    (needed if option==2 or if option==4)
  rprimd(3,3)=dimensional primitive translations in real space (bohr)
  ucvol=unit cell volume ($\textrm{Bohr}^3$).
  vprtrb(2)=complex amplitude of possible perturbing potential; if nonzero,
   perturbing potential is added of the form
   $V(G)=(vprtrb(1)+I*vprtrb(2))/2$ at the values G=qprtrb and
   $(vprtrb(1)-I*vprtrb(2))/2$ at $G=-qprtrb$ (integers)
  xred(3,natom)=reduced dimensionless atomic coordinates

OUTPUT

  (if option==1) vpsp(nfft)=local crystal pseudopotential in real space.
  (if option==2) grtn(3,natom)=grads of Etot wrt tn.
  (if option==3) lpsstr(6)=components of local psp part of stress tensor
   (Cartesian coordinates, symmetric tensor) in hartree/$\textrm{bohr}^3$
   Store 6 unique components in order 11, 22, 33, 32, 31, 21
  (if option==4) dyfrlo(3,3,natom)=d2 Eei/d tn(i)/d tn(j).  (Hartrees)

NOTES

 Note that the present routine is tightly connected to the dfpt_vlocal.f routine,
 that compute the derivative of the local ionic potential
 with respect to one atomic displacement. The argument list
 and the internal loops to be considered were sufficiently different
 as to make the two routine different.

PARENTS

      forces,prcref,prcref_PMA,respfn,setvtr

CHILDREN

      mklocl_realspace,mklocl_recipspace,mklocl_wavelets,wvl_rho_abi2big
      xred2xcart

SOURCE

134 subroutine mklocl(dtset, dyfrlo,eei,gmet,gprimd,grtn,gsqcut,lpsstr,mgfft,&
135 &  mpi_enreg,natom,nattyp,nfft,ngfft,nspden,ntypat,option,pawtab,ph1d,psps,qprtrb,&
136 &  rhog,rhor,rprimd,ucvol,vprtrb,vpsp,wvl,wvl_den,xred)
137 
138 
139 !This section has been created automatically by the script Abilint (TD).
140 !Do not modify the following lines by hand.
141 #undef ABI_FUNC
142 #define ABI_FUNC 'mklocl'
143 !End of the abilint section
144 
145  implicit none
146 
147 !Arguments ------------------------------------
148 !scalars
149  integer,intent(in) :: mgfft,natom,nfft,nspden,ntypat,option
150  real(dp),intent(in) :: eei,gsqcut,ucvol
151  type(MPI_type),intent(in) :: mpi_enreg
152  type(dataset_type),intent(in) :: dtset
153  type(pseudopotential_type),intent(in) :: psps
154  type(wvl_internal_type), intent(in) :: wvl
155  type(wvl_denspot_type), intent(inout) :: wvl_den
156  type(pawtab_type),intent(in)  :: pawtab(ntypat*psps%usepaw)
157 !arrays
158  integer,intent(in) :: nattyp(ntypat),ngfft(18),qprtrb(3)
159  real(dp),intent(in) :: gmet(3,3),gprimd(3,3),ph1d(2,3*(2*mgfft+1)*natom)
160  real(dp),intent(in) :: rhog(2,nfft),rprimd(3,3)
161  real(dp),intent(in) :: vprtrb(2),xred(3,natom)
162  real(dp),intent(in),target :: rhor(nfft,nspden)
163  real(dp),intent(out) :: dyfrlo(3,3,natom),grtn(3,natom),lpsstr(6)
164  real(dp),intent(inout) :: vpsp(nfft)
165 
166 !Local variables-------------------------------
167 !scalars
168  character(len=500) :: message
169 !arrays
170  real(dp),allocatable :: xcart(:,:)
171 #if defined HAVE_BIGDFT
172  real(dp),pointer :: rhor_ptr(:,:)
173 #endif
174 
175 ! *************************************************************************
176 
177  if (option < 1 .or. option > 4) then
178    write(message,'(a,i0,a,a)')&
179 &   'From the calling routine, option=',option,ch10,&
180 &   'The only allowed values are between 1 and 4.'
181    MSG_ERROR(message)
182  end if
183  if (option > 2 .and. .not.psps%vlspl_recipSpace) then
184    write(message,'(a,i0,a,a,a,a)')&
185 &   'From the calling routine, option=',option,ch10,&
186 &   'but the local part of the pseudo-potential is in real space.',ch10,&
187 &   'Action: set icoulomb = 0 to turn-off real space computations.'
188    MSG_ERROR(message)
189  end if
190  if (option > 2 .and. dtset%usewvl == 1) then
191    write(message,'(a,i0,a,a)')&
192 &   'From the calling routine, option=',option,ch10,&
193 &   'but this is not implemented yet from wavelets.'
194    MSG_ERROR(message)
195  end if
196 
197  if (dtset%usewvl == 0) then
198 !  Plane wave case
199    if (psps%vlspl_recipSpace) then
200      call mklocl_recipspace(dyfrlo,eei,gmet,gprimd,grtn,gsqcut,lpsstr,mgfft, &
201 &     mpi_enreg,psps%mqgrid_vl,natom,nattyp,nfft,ngfft, &
202 &     ntypat,option,dtset%paral_kgb,ph1d,psps%qgrid_vl,qprtrb,rhog,ucvol, &
203 &     psps%vlspl,vprtrb,vpsp)
204    else
205      call mklocl_realspace(grtn,dtset%icoulomb,mpi_enreg,natom,nattyp,nfft, &
206 &     ngfft,dtset%nscforder,nspden,ntypat,option,pawtab,psps,rhog,rhor, &
207 &     rprimd,dtset%typat,ucvol,dtset%usewvl,vpsp,xred)
208    end if
209  else
210 !  Store xcart for each atom
211    ABI_ALLOCATE(xcart,(3, dtset%natom))
212    call xred2xcart(dtset%natom, rprimd, xcart, xred)
213 !  Eventually retrieve density
214 #if defined HAVE_BIGDFT
215    if (option>1.and.wvl_den%denspot%rhov_is/=ELECTRONIC_DENSITY) then
216      rhor_ptr => rhor ! Just to bypass intent(inout)
217      call wvl_rho_abi2big(1,rhor_ptr,wvl_den)
218    end if
219 #endif
220 !  Wavelets case
221    call mklocl_wavelets(dtset%efield, grtn, mpi_enreg, dtset%natom, &
222 &   nfft, nspden, option, rprimd, vpsp, &
223 &   wvl_den, wvl, xcart)
224    ABI_DEALLOCATE(xcart)
225  end if
226 
227 end subroutine mklocl

ABINIT/mklocl_recipspace [ Functions ]

[ Top ] [ Functions ]

NAME

 mklocl_recipspace

FUNCTION

 Optionally compute :
  option=1 : local ionic potential throughout unit cell
  option=2 : contribution of local ionic potential to E gradient wrt xred
  option=3 : contribution of local ionic potential to stress tensor
  option=4 : contribution of local ionic potential to
                second derivative of E wrt xred

INPUTS

  if(option==3) eei=local pseudopotential part of total energy (hartree)
  gmet(3,3)=reciprocal space metric ($\textrm{Bohr}^{-2}$).
  gprimd(3,3)=reciprocal space dimensional primitive translations
  gsqcut=cutoff on $|G|^2$: see setup1 for definition (doubled sphere).
  mgfft=maximum size of 1D FFTs
  mpi_enreg=information about MPI parallelization
  mqgrid=number of grid pts in q array for f(q) spline.
  natom=number of atoms in unit cell.
  nattyp(ntypat)=number of atoms of each type in cell.
  nfft=(effective) number of FFT grid points (for this processor)
  ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
  ntypat=number of types of atoms.
  option= (see above)
  ph1d(2,3*(2*mgfft+1)*natom)=1-dim structure factor phase information.
  qgrid(mqgrid)=q grid for spline from 0 to qmax.
  qprtrb(3)= integer wavevector of possible perturbing potential
   in basis of reciprocal lattice translations
  rhog(2,nfft)=electron density rho(G) (electrons/$\textrm{Bohr}^3$)
    (needed if option==2 or if option==4)
  ucvol=unit cell volume ($\textrm{Bohr}^3$).
  vlspl(mqgrid,2,ntypat)=q^2 v(q) spline for each type of atom.
  vprtrb(2)=complex amplitude of possible perturbing potential; if nonzero,
   perturbing potential is added of the form
   $V(G)=(vprtrb(1)+I*vprtrb(2))/2$ at the values G=qprtrb and
   $(vprtrb(1)-I*vprtrb(2))/2$ at $G=-qprtrb$ (integers)

OUTPUT

  (if option==1) vpsp(nfft)=local crystal pseudopotential in real space.
  (if option==2) grtn(3,natom)=grads of Etot wrt tn.
  (if option==3) lpsstr(6)=components of local psp part of stress tensor
   (Cartesian coordinates, symmetric tensor) in hartree/$\textrm{bohr}^3$
   Store 6 unique components in order 11, 22, 33, 32, 31, 21
  (if option==4) dyfrlo(3,3,natom)=d2 Eei/d tn(i)/d tn(j).  (Hartrees)

NOTES

 Note that the present routine is tightly connected to the dfpt_vlocal.f routine,
 that compute the derivative of the local ionic potential
 with respect to one atomic displacement. The argument list
 and the internal loops to be considered were sufficiently different
 as to make the two routine different.

PARENTS

      dfpt_dyfro,mklocl,stress

CHILDREN

      fourdp,ptabs_fourdp,timab,wrtout,xmpi_sum

SOURCE

292 subroutine mklocl_recipspace(dyfrlo,eei,gmet,gprimd,grtn,gsqcut,lpsstr,mgfft,&
293 &  mpi_enreg,mqgrid,natom,nattyp,nfft,ngfft,ntypat,option,paral_kgb,ph1d,qgrid,qprtrb,&
294 &  rhog,ucvol,vlspl,vprtrb,vpsp)
295 
296 
297 !This section has been created automatically by the script Abilint (TD).
298 !Do not modify the following lines by hand.
299 #undef ABI_FUNC
300 #define ABI_FUNC 'mklocl_recipspace'
301 !End of the abilint section
302 
303  implicit none
304 
305 !Arguments ------------------------------------
306 !scalars
307  integer,intent(in) :: mgfft,mqgrid,natom,nfft,ntypat,option,paral_kgb
308  real(dp),intent(in) :: eei,gsqcut,ucvol
309  type(MPI_type),intent(in) :: mpi_enreg
310 !arrays
311  integer,intent(in) :: nattyp(ntypat),ngfft(18),qprtrb(3)
312  real(dp),intent(in) :: gmet(3,3),gprimd(3,3),ph1d(2,3*(2*mgfft+1)*natom)
313  real(dp),intent(in) :: qgrid(mqgrid),rhog(2,nfft),vlspl(mqgrid,2,ntypat)
314  real(dp),intent(in) :: vprtrb(2)
315  real(dp),intent(out) :: dyfrlo(3,3,natom),grtn(3,natom),lpsstr(6) !vz_i
316  real(dp),intent(inout) :: vpsp(nfft) !vz_i
317 
318 !Local variables-------------------------------
319 !scalars
320  integer,parameter :: im=2,re=1
321  integer :: i1,i2,i3,ia,ia1,ia2,id1,id2,id3,ierr,ig1,ig2,ig3,ii,itypat
322  integer :: jj,me_fft,me_g0,n1,n2,n3,nproc_fft,shift1
323  integer :: shift2,shift3
324  real(dp),parameter :: tolfix=1.0000001_dp
325  real(dp) :: aa,bb,cc,cutoff,dbl_ig1,dbl_ig2,dbl_ig3,dd,diff,dq,dq2div6,dqdiv6
326  real(dp) :: dqm1,ee,ff,gmag,gsquar,ph12i,ph12r,ph1i,ph1r,ph2i,ph2r
327  real(dp) :: ph3i,ph3r,phimag_igia,phre_igia,sfi,sfr
328  real(dp) :: svion,svioni,svionr,term,vion1,vion2,xnorm
329  character(len=500) :: message
330 !arrays
331  integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:)
332  integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:)
333  real(dp) :: gcart(3),tsec(2)
334  real(dp),allocatable :: work1(:,:)
335 
336 ! *************************************************************************
337 
338 !Define G^2 based on G space metric gmet.
339 ! gsq(i1,i2,i3)=dble(i1*i1)*gmet(1,1)+dble(i2*i2)*gmet(2,2)+&
340 !& dble(i3*i3)*gmet(3,3)+dble(2*i1*i2)*gmet(1,2)+&
341 !& dble(2*i2*i3)*gmet(2,3)+dble(2*i3*i1)*gmet(3,1)
342 
343 !Real and imaginary parts of phase--statment functions:
344 ! phr(x1,y1,x2,y2,x3,y3)=(x1*x2-y1*y2)*x3-(y1*x2+x1*y2)*y3
345 ! phi(x1,y1,x2,y2,x3,y3)=(x1*x2-y1*y2)*y3+(y1*x2+x1*y2)*x3
346 ! ph1(nri,i1,ia)=ph1d(nri,i1+1+n1+(ia-1)*(2*n1+1))
347 ! ph2(nri,i2,ia)=ph1d(nri,i2+1+n2+(ia-1)*(2*n2+1)+&
348 !& natom*(2*n1+1))
349 ! ph3(nri,i3,ia)=ph1d(nri,i3+1+n3+(ia-1)*(2*n3+1)+&
350 !& natom*(2*n1+1+2*n2+1))
351 ! phre(i1,i2,i3,ia)=phr(ph1(re,i1,ia),ph1(im,i1,ia),ph2(re,i2,ia),&
352 !& ph2(im,i2,ia),ph3(re,i3,ia),ph3(im,i3,ia))
353 ! phimag(i1,i2,i3,ia)=phi(ph1(re,i1,ia),ph1(im,i1,ia),ph2(re,i2,ia),&
354 !& ph2(im,i2,ia),ph3(re,i3,ia),ph3(im,i3,ia))
355 
356 !-----
357 
358 !Keep track of total time spent in mklocl
359  if(option==2)then
360    call timab(72,1,tsec)
361  end if
362  n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3)
363  me_fft=ngfft(11)
364  nproc_fft=ngfft(10)
365 
366 !Get the distrib associated with this fft_grid
367  call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
368 
369 !Zero out array to permit accumulation over atom types below:
370  if(option==1)then
371    ABI_ALLOCATE(work1,(2,nfft))
372    work1(:,:)=zero
373  end if
374 !
375  dq=(qgrid(mqgrid)-qgrid(1))/dble(mqgrid-1)
376  dqm1=1.0_dp/dq
377  dqdiv6=dq/6.0_dp
378  dq2div6=dq**2/6.0_dp
379  cutoff=gsqcut*tolfix
380  id1=n1/2+2
381  id2=n2/2+2
382  id3=n3/2+2
383  grtn(:,:)=zero
384  lpsstr(:)=zero
385  dyfrlo(:,:,:)=zero
386  me_g0=0
387  ia1=1
388 
389  do itypat=1,ntypat
390 !  ia1,ia2 sets range of loop over atoms:
391    ia2=ia1+nattyp(itypat)-1
392 
393    ii=0
394    do i3=1,n3
395      ig3=i3-(i3/id3)*n3-1
396      do i2=1,n2
397        ig2=i2-(i2/id2)*n2-1
398        if(fftn2_distrib(i2) == me_fft ) then
399          do i1=1,n1
400            ig1=i1-(i1/id1)*n1-1
401 
402            ii=ii+1
403 !          ***     GET RID OF THIS THESE IF STATEMENTS (if they slow code)
404 !          Skip G=0:
405 !          if (ii==1) cycle
406            if (ig1==0 .and. ig2==0 .and. ig3==0) me_g0=1
407            if (ig1==0 .and. ig2==0 .and. ig3==0) cycle
408 
409            gsquar=gsq_mk(ig1,ig2,ig3)
410 !          Skip G**2 outside cutoff:
411            if (gsquar<=cutoff) then
412              gmag=sqrt(gsquar)
413 
414 !            Compute vion(G) for given type of atom
415              jj=1+int(gmag*dqm1)
416              diff=gmag-qgrid(jj)
417 
418 !            Evaluate spline fit from q^2 V(q) to get V(q):
419 !            (p. 86 Numerical Recipes, Press et al;
420 !            NOTE error in book for sign
421 !            of "aa" term in derivative; also see splfit routine).
422 
423              bb = diff*dqm1
424              aa = 1.0_dp-bb
425              cc = aa*(aa**2-1.0_dp)*dq2div6
426              dd = bb*(bb**2-1.0_dp)*dq2div6
427 
428              vion1 = (aa*vlspl(jj,1,itypat)+bb*vlspl(jj+1,1,itypat) +&
429 &             cc*vlspl(jj,2,itypat)+dd*vlspl(jj+1,2,itypat) ) / gsquar
430 
431              if(option==1)then
432 
433 !              Assemble structure factor over all atoms of given type:
434                sfr=zero
435                sfi=zero
436                do ia=ia1,ia2
437                  sfr=sfr+phre_mk(ig1,ig2,ig3,ia)
438                  sfi=sfi-phimag_mk(ig1,ig2,ig3,ia)
439                end do
440 !              Multiply structure factor times vion:
441                work1(re,ii)=work1(re,ii)+sfr*vion1
442                work1(im,ii)=work1(im,ii)+sfi*vion1
443 
444              else if(option==2 .or. option==4)then
445 
446 !              Compute Re and Im part of (2Pi)*Vion(G)*rho(G):
447                svionr=(two_pi*vion1)*rhog(re,ii)
448                svioni=(two_pi*vion1)*rhog(im,ii)
449 
450 !              Loop over all atoms of this type:
451                do ia=ia1,ia2
452                  shift1=1+n1+(ia-1)*(2*n1+1)
453                  shift2=1+n2+(ia-1)*(2*n2+1)+natom*(2*n1+1)
454                  shift3=1+n3+(ia-1)*(2*n3+1)+natom*(2*n1+1+2*n2+1)
455                  ph1r=ph1d(1,ig1+shift1)
456                  ph1i=ph1d(2,ig1+shift1)
457                  ph2r=ph1d(1,ig2+shift2)
458                  ph2i=ph1d(2,ig2+shift2)
459                  ph3r=ph1d(1,ig3+shift3)
460                  ph3i=ph1d(2,ig3+shift3)
461                  ph12r=ph1r*ph2r-ph1i*ph2i
462                  ph12i=ph1r*ph2i+ph1i*ph2r
463                  phre_igia=ph12r*ph3r-ph12i*ph3i
464                  phimag_igia=ph12r*ph3i+ph12i*ph3r
465 
466                  if(option==2)then
467 
468 !                  Compute "Vion" part of gradient
469 !                  svion=svioni*phre(ig1,ig2,ig3,ia)+svionr*phimag(ig1,ig2,ig3,ia)
470                    svion=svioni*phre_igia+svionr*phimag_igia
471 
472 !                  Open loop over 3-index for speed:
473                    grtn(1,ia)=grtn(1,ia)-dble(ig1)*svion
474                    grtn(2,ia)=grtn(2,ia)-dble(ig2)*svion
475                    grtn(3,ia)=grtn(3,ia)-dble(ig3)*svion
476 
477                  else
478 
479 !                  Compute "Vion" part of the second derivative
480 !                  svion=two_pi*
481 !                  (svionr*phre(ig1,ig2,ig3,ia)-svioni*phimag(ig1,ig2,ig3,ia))
482                    svion=two_pi*(svionr*phre_igia-svioni*phimag_igia)
483 
484 !                  Open loop over 3-index for speed
485                    dbl_ig1=dble(ig1) ; dbl_ig2=dble(ig2) ; dbl_ig3=dble(ig3)
486                    dyfrlo(1,1,ia)=dyfrlo(1,1,ia)-dbl_ig1*dbl_ig1*svion
487                    dyfrlo(1,2,ia)=dyfrlo(1,2,ia)-dbl_ig1*dbl_ig2*svion
488                    dyfrlo(1,3,ia)=dyfrlo(1,3,ia)-dbl_ig1*dbl_ig3*svion
489                    dyfrlo(2,2,ia)=dyfrlo(2,2,ia)-dbl_ig2*dbl_ig2*svion
490                    dyfrlo(2,3,ia)=dyfrlo(2,3,ia)-dbl_ig2*dbl_ig3*svion
491                    dyfrlo(3,3,ia)=dyfrlo(3,3,ia)-dbl_ig3*dbl_ig3*svion
492 
493                  end if
494 
495                end do
496 
497              else if(option==3)then
498 
499 !              Also get (dV(q)/dq)/q:
500 !              (note correction of Numerical Recipes sign error
501 !              before (3._dp*aa**2-1._dp)
502 !              ee*dqm1 + ff*dqdiv6 is the best estimate of dV(q)/dq from splines
503                ee= vlspl(jj+1,1,itypat)-vlspl(jj,1,itypat)
504                ff=  (3._dp*bb**2-1._dp)*vlspl(jj+1,2,itypat) &
505 &               - (3._dp*aa**2-1._dp)*vlspl(jj,2,itypat)
506                vion2 = ( ( ee*dqm1 + ff*dqdiv6 )/gmag&
507 &               - 2.0_dp*vion1                 ) / gsquar
508 
509                gcart(1)=gprimd(1,1)*dble(ig1)+gprimd(1,2)*dble(ig2)+&
510 &               gprimd(1,3)*dble(ig3)
511                gcart(2)=gprimd(2,1)*dble(ig1)+gprimd(2,2)*dble(ig2)+&
512 &               gprimd(2,3)*dble(ig3)
513                gcart(3)=gprimd(3,1)*dble(ig1)+gprimd(3,2)*dble(ig2)+&
514 &               gprimd(3,3)*dble(ig3)
515 !              Assemble structure over all atoms of given type
516                sfr=zero
517                sfi=zero
518                do ia=ia1,ia2
519                  sfr=sfr+phre_mk(ig1,ig2,ig3,ia)
520                  sfi=sfi-phimag_mk(ig1,ig2,ig3,ia)
521                end do
522 
523 !              Compute Re( rho^*(G)* sf ) * [(dV(G)/dG)/|G|]
524                term=(rhog(re,ii)*sfr+rhog(im,ii)*sfi)*vion2
525 
526 !              Compute contribution to stress tensor
527                lpsstr(1)=lpsstr(1)-term*gcart(1)*gcart(1)
528                lpsstr(2)=lpsstr(2)-term*gcart(2)*gcart(2)
529                lpsstr(3)=lpsstr(3)-term*gcart(3)*gcart(3)
530                lpsstr(4)=lpsstr(4)-term*gcart(3)*gcart(2)
531                lpsstr(5)=lpsstr(5)-term*gcart(3)*gcart(1)
532                lpsstr(6)=lpsstr(6)-term*gcart(2)*gcart(1)
533 
534              else
535                write(message, '(a,i0,a)' )' mklocl: Option=',option,' not allowed.'
536                MSG_BUG(message)
537              end if ! End option choice
538 
539 !            End skip G**2 outside cutoff:
540            end if
541 
542 !          End loop on n1, n2, n3. There is a "cycle" inside the loop
543          end do
544        end if ! this plane is for me_fft
545      end do
546    end do
547 
548 !  Symmetrize the dynamical matrix with respect to indices
549    do ia=ia1,ia2
550      dyfrlo(2,1,ia)=dyfrlo(1,2,ia)
551      dyfrlo(3,1,ia)=dyfrlo(1,3,ia)
552      dyfrlo(3,2,ia)=dyfrlo(2,3,ia)
553    end do
554 
555    ia1=ia2+1
556 
557 !  End loop on type of atoms
558  end do
559 
560  if(option==1)then
561 !  Dont't change work1 on g=0 if Poisson solver is used since work1
562 !  hold not the potential but the density generated by the pseudo.
563    if(me_g0 == 1) then
564 !    Set Vloc(G=0)=0:
565      work1(re,1)=zero
566      work1(im,1)=zero
567    end if
568 
569 !  DEBUG
570 !  write(std_out,*) ' mklocl_recipspace : will add potential with strength vprtrb(:)=',vprtrb(:)
571 !  ENDDEBUG
572 
573 !  Allow for the addition of a perturbing potential
574    if ((vprtrb(1)**2+vprtrb(2)**2) > 1.d-30) then
575 !    Find the linear indices which correspond with the input
576 !    wavevector qprtrb
577 !    The double modulus handles both i>=n and i<0, mapping into [0,n-1];
578 !    then add 1 to get range [1,n] for each
579      i3=1+mod(n3+mod(qprtrb(3),n3),n3)
580      i2=1+mod(n2+mod(qprtrb(2),n2),n2)
581      i1=1+mod(n1+mod(qprtrb(1),n1),n1)
582 !    Compute the linear index in the 3 dimensional array
583      ii=i1+n1*((ffti2_local(i2)-1)+(n2/nproc_fft)*(i3-1))
584 !    Add in the perturbation at G=qprtrb
585      work1(re,ii)=work1(re,ii)+0.5_dp*vprtrb(1)
586      work1(im,ii)=work1(im,ii)+0.5_dp*vprtrb(2)
587 !    Same thing for G=-qprtrb
588      i3=1+mod(n3+mod(-qprtrb(3),n3),n3)
589      i2=1+mod(n2+mod(-qprtrb(2),n2),n2)
590      i1=1+mod(n1+mod(-qprtrb(1),n1),n1)
591 !    ii=i1+n1*((i2-1)+n2*(i3-1))
592      work1(re,ii)=work1(re,ii)+0.5_dp*vprtrb(1)
593      work1(im,ii)=work1(im,ii)-0.5_dp*vprtrb(2)
594      write(message, '(a,1p,2e12.4,a,0p,3i4,a)' )&
595 &     ' mklocl: perturbation of vprtrb=', vprtrb,&
596 &     ' and q=',qprtrb,' has been added'
597      call wrtout(std_out,message,'COLL')
598    end if
599 
600 !  Transform back to real space
601    call fourdp(1,work1,vpsp,1,mpi_enreg,nfft,ngfft,paral_kgb,0)
602 
603 !  Divide by unit cell volume
604    xnorm=1.0_dp/ucvol
605    vpsp(:)=vpsp(:)*xnorm
606 
607    ABI_DEALLOCATE(work1)
608 
609  end if
610 
611  if(option==2)then
612 !  Init mpi_comm
613    if(mpi_enreg%nproc_fft>1)then
614      call timab(48,1,tsec)
615      call xmpi_sum(grtn,mpi_enreg%comm_fft ,ierr)
616      call timab(48,2,tsec)
617    end if
618    call timab(72,2,tsec)
619  end if
620 
621  if(option==3)then
622 !  Init mpi_comm
623    if(mpi_enreg%nproc_fft>1)then
624      call timab(48,1,tsec)
625      call xmpi_sum(lpsstr,mpi_enreg%comm_fft ,ierr)
626      call timab(48,2,tsec)
627    end if
628 
629 !  Normalize and add term -eei/ucvol on diagonal
630 !  (see page 802 of notes)
631    lpsstr(1)=(lpsstr(1)-eei)/ucvol
632    lpsstr(2)=(lpsstr(2)-eei)/ucvol
633    lpsstr(3)=(lpsstr(3)-eei)/ucvol
634    lpsstr(4)=lpsstr(4)/ucvol
635    lpsstr(5)=lpsstr(5)/ucvol
636    lpsstr(6)=lpsstr(6)/ucvol
637 
638  end if
639 
640  if(option==4)then
641 !  Init mpi_comm
642    if(mpi_enreg%nproc_fft>1)then
643      call timab(48,1,tsec)
644      call xmpi_sum(dyfrlo,mpi_enreg%comm_fft ,ierr)
645      call timab(48,2,tsec)
646    end if
647  end if
648 
649  contains
650 
651 !Real and imaginary parts of phase--statment functions:
652    function phr_mk(x1,y1,x2,y2,x3,y3)
653 
654 
655 !This section has been created automatically by the script Abilint (TD).
656 !Do not modify the following lines by hand.
657 #undef ABI_FUNC
658 #define ABI_FUNC 'phr_mk'
659 !End of the abilint section
660 
661    real(dp) :: phr_mk,x1,x2,x3,y1,y2,y3
662    phr_mk=(x1*x2-y1*y2)*x3-(y1*x2+x1*y2)*y3
663  end function phr_mk
664 
665    function phi_mk(x1,y1,x2,y2,x3,y3)
666 
667 
668 !This section has been created automatically by the script Abilint (TD).
669 !Do not modify the following lines by hand.
670 #undef ABI_FUNC
671 #define ABI_FUNC 'phi_mk'
672 !End of the abilint section
673 
674    real(dp):: phi_mk,x1,x2,x3,y1,y2,y3
675    phi_mk=(x1*x2-y1*y2)*y3+(y1*x2+x1*y2)*x3
676  end function phi_mk
677 
678    function ph1_mk(nri,ig1,ia)
679 
680 
681 !This section has been created automatically by the script Abilint (TD).
682 !Do not modify the following lines by hand.
683 #undef ABI_FUNC
684 #define ABI_FUNC 'ph1_mk'
685 !End of the abilint section
686 
687    real(dp):: ph1_mk
688    integer :: nri,ig1,ia
689    ph1_mk=ph1d(nri,ig1+1+n1+(ia-1)*(2*n1+1))
690  end function ph1_mk
691 
692    function ph2_mk(nri,ig2,ia)
693 
694 
695 !This section has been created automatically by the script Abilint (TD).
696 !Do not modify the following lines by hand.
697 #undef ABI_FUNC
698 #define ABI_FUNC 'ph2_mk'
699 !End of the abilint section
700 
701    real(dp):: ph2_mk
702    integer :: nri,ig2,ia
703    ph2_mk=ph1d(nri,ig2+1+n2+(ia-1)*(2*n2+1)+natom*(2*n1+1))
704  end function ph2_mk
705 
706    function ph3_mk(nri,ig3,ia)
707 
708 
709 !This section has been created automatically by the script Abilint (TD).
710 !Do not modify the following lines by hand.
711 #undef ABI_FUNC
712 #define ABI_FUNC 'ph3_mk'
713 !End of the abilint section
714 
715    real(dp):: ph3_mk
716    integer :: nri,ig3,ia
717    ph3_mk=ph1d(nri,ig3+1+n3+(ia-1)*(2*n3+1)+natom*(2*n1+1+2*n2+1))
718  end function ph3_mk
719 
720    function phre_mk(ig1,ig2,ig3,ia)
721 
722 
723 !This section has been created automatically by the script Abilint (TD).
724 !Do not modify the following lines by hand.
725 #undef ABI_FUNC
726 #define ABI_FUNC 'phre_mk'
727 !End of the abilint section
728 
729    real(dp):: phre_mk
730    integer :: ig1,ig2,ig3,ia
731    phre_mk=phr_mk(ph1_mk(re,ig1,ia),ph1_mk(im,ig1,ia),&
732 &   ph2_mk(re,ig2,ia),ph2_mk(im,ig2,ia),ph3_mk(re,ig3,ia),ph3_mk(im,ig3,ia))
733  end function phre_mk
734 
735    function phimag_mk(ig1,ig2,ig3,ia)
736 
737 
738 !This section has been created automatically by the script Abilint (TD).
739 !Do not modify the following lines by hand.
740 #undef ABI_FUNC
741 #define ABI_FUNC 'phimag_mk'
742 !End of the abilint section
743 
744    real(dp) :: phimag_mk
745    integer :: ig1,ig2,ig3,ia
746    phimag_mk=phi_mk(ph1_mk(re,ig1,ia),ph1_mk(im,ig1,ia),&
747 &   ph2_mk(re,ig2,ia),ph2_mk(im,ig2,ia),ph3_mk(re,ig3,ia),ph3_mk(im,ig3,ia))
748  end function phimag_mk
749 
750    function gsq_mk(i1,i2,i3)
751 
752 
753 !This section has been created automatically by the script Abilint (TD).
754 !Do not modify the following lines by hand.
755 #undef ABI_FUNC
756 #define ABI_FUNC 'gsq_mk'
757 !End of the abilint section
758 
759    real(dp) :: gsq_mk
760    integer :: i1,i2,i3
761    gsq_mk=dble(i1*i1)*gmet(1,1)+dble(i2*i2)*gmet(2,2)+&
762 &   dble(i3*i3)*gmet(3,3)+dble(2*i1*i2)*gmet(1,2)+&
763 &   dble(2*i2*i3)*gmet(2,3)+dble(2*i3*i1)*gmet(3,1)
764  end function gsq_mk
765 
766 end subroutine mklocl_recipspace

ABINIT/vlocalstr [ Functions ]

[ Top ] [ Functions ]

NAME

 vlocalstr

FUNCTION

 Compute strain derivatives of local ionic potential
                second derivative of E wrt xred

INPUTS

  gmet(3,3)=reciprocal space metric ($\textrm{Bohr}^{-2}$).
  gprimd(3,3)=reciprocal space dimensional primitive translations
  gsqcut=cutoff on $|G|^2$: see setup1 for definition (doubled sphere).
  istr=1,...6 specifies cartesian strain component 11,22,33,32,31,21
  mgfft=maximum size of 1D FFTs
  mpi_enreg=information about MPI parallelization
  mqgrid=number of grid pts in q array for f(q) spline.
  natom=number of atoms in unit cell.
  nattyp(ntypat)=number of atoms of each type in cell.
  nfft=(effective) number of FFT grid points (for this processor)
  ngfft(18)=contain all needed information about 3D FFT,
    see ~abinit/doc/variables/vargs.htm#ngfft
  ntypat=number of types of atoms.
  ph1d(2,3*(2*mgfft+1)*natom)=1-dim structure factor phase information.
  qgrid(mqgrid)=q grid for spline from 0 to qmax.
  ucvol=unit cell volume ($\textrm{Bohr}^3$).
  vlspl(mqgrid,2,ntypat)=q^2 v(q) spline for each type of atom.

OUTPUT

  vpsp1(nfft)=first-order local crystal pseudopotential in real space.

NOTES

 * Note that the present routine is tightly connected to the dfpt_vlocal.f routine,
 that compute the derivative of the local ionic potential
 with respect to one atomic displacement. The argument list
 and the internal loops to be considered were sufficiently different
 as to make the two routines different.
 * The routine was adapted from mklocl.F90

PARENTS

      dfpt_looppert,dfpt_nselt,dfpt_nstpaw

CHILDREN

      fourdp,ptabs_fourdp

SOURCE

1143 subroutine vlocalstr(gmet,gprimd,gsqcut,istr,mgfft,mpi_enreg,&
1144 &  mqgrid,natom,nattyp,nfft,ngfft,ntypat,paral_kgb,ph1d,qgrid,&
1145 &  ucvol,vlspl,vpsp1)
1146 
1147 
1148 !This section has been created automatically by the script Abilint (TD).
1149 !Do not modify the following lines by hand.
1150 #undef ABI_FUNC
1151 #define ABI_FUNC 'vlocalstr'
1152 !End of the abilint section
1153 
1154  implicit none
1155 
1156 !Arguments ------------------------------------
1157 !scalars
1158  integer,intent(in) :: istr,mgfft,mqgrid,natom,nfft,ntypat,paral_kgb
1159  real(dp),intent(in) :: gsqcut,ucvol
1160  type(MPI_type),intent(in) :: mpi_enreg
1161 !arrays
1162  integer,intent(in) :: nattyp(ntypat),ngfft(18)
1163  real(dp),intent(in) :: gmet(3,3),gprimd(3,3),ph1d(2,3*(2*mgfft+1)*natom)
1164  real(dp),intent(in) :: qgrid(mqgrid),vlspl(mqgrid,2,ntypat)
1165  real(dp),intent(out) :: vpsp1(nfft)
1166 
1167 !Local variables-------------------------------
1168 !scalars
1169  integer,parameter :: im=2,re=1
1170  integer :: i1,i2,i3,ia,ia1,ia2,id1,id2,id3,ig1,ig2,ig3,ii,itypat,jj
1171  integer :: ka,kb,n1,n2,n3
1172  real(dp),parameter :: tolfix=1.0000001_dp
1173  real(dp) :: aa,bb,cc,cutoff,dd,dgsquards,diff
1174  real(dp) :: dq,dq2div6,dqdiv6,dqm1,ee,ff,gmag,gsquar
1175  real(dp) :: sfi,sfr,term,vion1,vion2
1176  real(dp) :: xnorm
1177  character(len=500) :: message
1178 !arrays
1179  integer,save :: idx(12)=(/1,1,2,2,3,3,3,2,3,1,2,1/)
1180  integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:)
1181  integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:)
1182  real(dp) :: dgmetds(3,3)
1183  real(dp),allocatable :: work1(:,:)
1184 
1185 ! *************************************************************************
1186 
1187 !Define G^2 based on G space metric gmet.
1188 ! gsq_vl(i1,i2,i3)=dble(i1*i1)*gmet(1,1)+dble(i2*i2)*gmet(2,2)+&
1189 !& dble(i3*i3)*gmet(3,3)+dble(2*i1*i2)*gmet(1,2)+&
1190 !& dble(2*i2*i3)*gmet(2,3)+dble(2*i3*i1)*gmet(3,1)
1191 
1192 !Define dG^2/ds based on G space metric derivative dgmetds.
1193 ! dgsqds_vl(i1,i2,i3)=dble(i1*i1)*dgmetds(1,1)+dble(i2*i2)*dgmetds(2,2)+&
1194 !& dble(i3*i3)*dgmetds(3,3)+&
1195 !& dble(i1*i2)*(dgmetds(1,2)+dgmetds(2,1))+&
1196 !& dble(i1*i3)*(dgmetds(1,3)+dgmetds(3,1))+&
1197 !& dble(i2*i3)*(dgmetds(2,3)+dgmetds(3,2))
1198 
1199 !Real and imaginary parts of phase--statment functions:
1200 ! phr_vl(x1,y1,x2,y2,x3,y3)=(x1*x2-y1*y2)*x3-(y1*x2+x1*y2)*y3
1201 ! phi_vl(x1,y1,x2,y2,x3,y3)=(x1*x2-y1*y2)*y3+(y1*x2+x1*y2)*x3
1202 ! ph1_vl(nri,i1,ia)=ph1d(nri,i1+1+n1+(ia-1)*(2*n1+1))
1203 ! ph2_vl(nri,i2,ia)=ph1d(nri,i2+1+n2+(ia-1)*(2*n2+1)+&
1204 !& natom*(2*n1+1))
1205 ! ph3_vl(nri,i3,ia)=ph1d(nri,i3+1+n3+(ia-1)*(2*n3+1)+&
1206 !& natom*(2*n1+1+2*n2+1))
1207 ! phre_vl(i1,i2,i3,ia)=phr_vl(ph1_vl(re,i1,ia),ph1_vl(im,i1,ia),ph2_vl(re,i2,ia),&
1208 !& ph2_vl(im,i2,ia),ph3_vl(re,i3,ia),ph3_vl(im,i3,ia))
1209 ! phimag_vl(i1,i2,i3,ia)=phi_vl(ph1_vl(re,i1,ia),ph1_vl(im,i1,ia),ph2_vl(re,i2,ia),&
1210 !& ph2_vl(im,i2,ia),ph3_vl(re,i3,ia),ph3_vl(im,i3,ia))
1211 
1212 !-----
1213 !Compute derivative of metric tensor wrt strain component istr
1214  if(istr<1 .or. istr>6)then
1215    write(message, '(a,i10,a,a,a)' )&
1216 &   ' Input istr=',istr,' not allowed.',ch10,&
1217 &   ' Possible values are 1,2,3,4,5,6 only.'
1218    MSG_BUG(message)
1219  end if
1220 
1221  ka=idx(2*istr-1);kb=idx(2*istr)
1222  do ii = 1,3
1223    dgmetds(:,ii)=-(gprimd(ka,:)*gprimd(kb,ii)+gprimd(kb,:)*gprimd(ka,ii))
1224  end do
1225 !For historical reasons:
1226  dgmetds(:,:)=0.5_dp*dgmetds(:,:)
1227 
1228  n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3)
1229 
1230 !Get the distrib associated with this fft_grid
1231  call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
1232 
1233 !Zero out array to permit accumulation over atom types below:
1234  ABI_ALLOCATE(work1,(2,nfft))
1235  work1(:,:)=0.0_dp
1236 !
1237  dq=(qgrid(mqgrid)-qgrid(1))/dble(mqgrid-1)
1238  dqm1=1.0_dp/dq
1239  dqdiv6=dq/6.0_dp
1240  dq2div6=dq**2/6.0_dp
1241  cutoff=gsqcut*tolfix
1242  id1=n1/2+2
1243  id2=n2/2+2
1244  id3=n3/2+2
1245 
1246  ia1=1
1247  do itypat=1,ntypat
1248 !  ia1,ia2 sets range of loop over atoms:
1249    ia2=ia1+nattyp(itypat)-1
1250 
1251    ii=0
1252    do i3=1,n3
1253      ig3=i3-(i3/id3)*n3-1
1254      do i2=1,n2
1255        if (fftn2_distrib(i2)==mpi_enreg%me_fft) then
1256          ig2=i2-(i2/id2)*n2-1
1257          do i1=1,n1
1258            ig1=i1-(i1/id1)*n1-1
1259            ii=ii+1
1260 !          ***     GET RID OF THIS THESE IF STATEMENTS (if they slow code)
1261 !          Skip G=0:
1262 !          if (ii==1) cycle
1263            if (ig1==0 .and. ig2==0 .and. ig3==0) cycle
1264            gsquar=gsq_vl(ig1,ig2,ig3)
1265 
1266 !          Skip G**2 outside cutoff:
1267            if (gsquar<=cutoff) then
1268              gmag=sqrt(gsquar)
1269              dgsquards=dgsqds_vl(ig1,ig2,ig3)
1270 !            Compute vion(G) for given type of atom
1271              jj=1+int(gmag*dqm1)
1272              diff=gmag-qgrid(jj)
1273 
1274 !            Evaluate spline fit from q^2 V(q) to get V(q):
1275 !            (p. 86 Numerical Recipes, Press et al;
1276 !            NOTE error in book for sign
1277 !            of "aa" term in derivative; also see splfit routine).
1278 
1279              bb = diff*dqm1
1280              aa = 1.0_dp-bb
1281              cc = aa*(aa**2-1.0_dp)*dq2div6
1282              dd = bb*(bb**2-1.0_dp)*dq2div6
1283 
1284              vion1 = (aa*vlspl(jj,1,itypat)+bb*vlspl(jj+1,1,itypat) +&
1285 &             cc*vlspl(jj,2,itypat)+dd*vlspl(jj+1,2,itypat) ) &
1286 &             / gsquar
1287 
1288 !            Also get (dV(q)/dq)/q:
1289 !            (note correction of Numerical Recipes sign error
1290 !            before (3._dp*aa**2-1._dp)
1291              ee= vlspl(jj+1,1,itypat)-vlspl(jj,1,itypat)
1292              ff=  (3._dp*bb**2-1._dp)*vlspl(jj+1,2,itypat) &
1293 &             - (3._dp*aa**2-1._dp)*vlspl(jj,2,itypat)
1294              vion2 = ( ( ee*dqm1 + ff*dqdiv6 )/gmag&
1295 &             - 2.0_dp*vion1                 ) / gsquar
1296 
1297 
1298 !            Assemble structure factor over all atoms of given type:
1299              sfr=0.0_dp
1300              sfi=0.0_dp
1301              do ia=ia1,ia2
1302                sfr=sfr+phre_vl(ig1,ig2,ig3,ia)
1303                sfi=sfi-phimag_vl(ig1,ig2,ig3,ia)
1304              end do
1305 
1306              term=dgsquards*vion2
1307 !            Add potential for diagonal strain components
1308              if(istr <=3) then
1309                term=term-vion1
1310              end if
1311 
1312 !            Multiply structure factor times vion derivatives:
1313              work1(re,ii)=work1(re,ii)+sfr*term
1314              work1(im,ii)=work1(im,ii)+sfi*term
1315 
1316 !            End skip G**2 outside cutoff:
1317            end if
1318 !          End loop on n1, n2, n3. There is a "cycle" inside the loop
1319          end do
1320        end if
1321      end do
1322    end do
1323 
1324    ia1=ia2+1
1325 
1326 !  End loop on type of atoms
1327  end do
1328 
1329 
1330 !Set Vloc(G=0)=0:
1331  work1(re,1)=0.0_dp
1332  work1(im,1)=0.0_dp
1333 
1334 
1335 !Transform back to real space
1336  call fourdp(1,work1,vpsp1,1,mpi_enreg,nfft,ngfft,paral_kgb,0)
1337 
1338 !Divide by unit cell volume
1339  xnorm=1.0_dp/ucvol
1340  vpsp1(:)=vpsp1(:)*xnorm
1341 
1342  ABI_DEALLOCATE(work1)
1343 
1344  contains
1345 
1346 !Real and imaginary parts of phase.
1347    function phr_vl(x1,y1,x2,y2,x3,y3)
1348 
1349 
1350 !This section has been created automatically by the script Abilint (TD).
1351 !Do not modify the following lines by hand.
1352 #undef ABI_FUNC
1353 #define ABI_FUNC 'phr_vl'
1354 !End of the abilint section
1355 
1356    real(dp) :: phr_vl,x1,x2,x3,y1,y2,y3
1357    phr_vl=(x1*x2-y1*y2)*x3-(y1*x2+x1*y2)*y3
1358  end function phr_vl
1359 
1360    function phi_vl(x1,y1,x2,y2,x3,y3)
1361 
1362 
1363 !This section has been created automatically by the script Abilint (TD).
1364 !Do not modify the following lines by hand.
1365 #undef ABI_FUNC
1366 #define ABI_FUNC 'phi_vl'
1367 !End of the abilint section
1368 
1369    real(dp):: phi_vl,x1,x2,x3,y1,y2,y3
1370    phi_vl=(x1*x2-y1*y2)*y3+(y1*x2+x1*y2)*x3
1371  end function phi_vl
1372 
1373    function ph1_vl(nri,ig1,ia)
1374 
1375 
1376 !This section has been created automatically by the script Abilint (TD).
1377 !Do not modify the following lines by hand.
1378 #undef ABI_FUNC
1379 #define ABI_FUNC 'ph1_vl'
1380 !End of the abilint section
1381 
1382    real(dp):: ph1_vl
1383    integer :: nri,ig1,ia
1384    ph1_vl=ph1d(nri,ig1+1+n1+(ia-1)*(2*n1+1))
1385  end function ph1_vl
1386 
1387    function ph2_vl(nri,ig2,ia)
1388 
1389 
1390 !This section has been created automatically by the script Abilint (TD).
1391 !Do not modify the following lines by hand.
1392 #undef ABI_FUNC
1393 #define ABI_FUNC 'ph2_vl'
1394 !End of the abilint section
1395 
1396    real(dp):: ph2_vl
1397    integer :: nri,ig2,ia
1398    ph2_vl=ph1d(nri,ig2+1+n2+(ia-1)*(2*n2+1)+natom*(2*n1+1))
1399  end function ph2_vl
1400 
1401    function ph3_vl(nri,ig3,ia)
1402 
1403 
1404 !This section has been created automatically by the script Abilint (TD).
1405 !Do not modify the following lines by hand.
1406 #undef ABI_FUNC
1407 #define ABI_FUNC 'ph3_vl'
1408 !End of the abilint section
1409 
1410    real(dp):: ph3_vl
1411    integer :: nri,ig3,ia
1412    ph3_vl=ph1d(nri,ig3+1+n3+(ia-1)*(2*n3+1)+natom*(2*n1+1+2*n2+1))
1413  end function ph3_vl
1414 
1415    function phre_vl(ig1,ig2,ig3,ia)
1416 
1417 
1418 !This section has been created automatically by the script Abilint (TD).
1419 !Do not modify the following lines by hand.
1420 #undef ABI_FUNC
1421 #define ABI_FUNC 'phre_vl'
1422 !End of the abilint section
1423 
1424    real(dp):: phre_vl
1425    integer :: ig1,ig2,ig3,ia
1426    phre_vl=phr_vl(ph1_vl(re,ig1,ia),ph1_vl(im,ig1,ia),&
1427 &   ph2_vl(re,ig2,ia),ph2_vl(im,ig2,ia),ph3_vl(re,ig3,ia),ph3_vl(im,ig3,ia))
1428  end function phre_vl
1429 
1430    function phimag_vl(ig1,ig2,ig3,ia)
1431 
1432 
1433 !This section has been created automatically by the script Abilint (TD).
1434 !Do not modify the following lines by hand.
1435 #undef ABI_FUNC
1436 #define ABI_FUNC 'phimag_vl'
1437 !End of the abilint section
1438 
1439    real(dp) :: phimag_vl
1440    integer :: ig1,ig2,ig3,ia
1441    phimag_vl=phi_vl(ph1_vl(re,ig1,ia),ph1_vl(im,ig1,ia),&
1442 &   ph2_vl(re,ig2,ia),ph2_vl(im,ig2,ia),ph3_vl(re,ig3,ia),ph3_vl(im,ig3,ia))
1443  end function phimag_vl
1444 
1445    function gsq_vl(i1,i2,i3)
1446 
1447 
1448 !This section has been created automatically by the script Abilint (TD).
1449 !Do not modify the following lines by hand.
1450 #undef ABI_FUNC
1451 #define ABI_FUNC 'gsq_vl'
1452 !End of the abilint section
1453 
1454    real(dp) :: gsq_vl
1455    integer :: i1,i2,i3
1456 !Define G^2 based on G space metric gmet.
1457    gsq_vl=dble(i1*i1)*gmet(1,1)+dble(i2*i2)*gmet(2,2)+&
1458 &   dble(i3*i3)*gmet(3,3)+dble(2*i1*i2)*gmet(1,2)+&
1459 &   dble(2*i2*i3)*gmet(2,3)+dble(2*i3*i1)*gmet(3,1)
1460  end function gsq_vl
1461 
1462    function dgsqds_vl(i1,i2,i3)
1463 
1464 
1465 !This section has been created automatically by the script Abilint (TD).
1466 !Do not modify the following lines by hand.
1467 #undef ABI_FUNC
1468 #define ABI_FUNC 'dgsqds_vl'
1469 !End of the abilint section
1470 
1471    real(dp) :: dgsqds_vl
1472    integer :: i1,i2,i3
1473 !Define dG^2/ds based on G space metric derivative dgmetds.
1474    dgsqds_vl=dble(i1*i1)*dgmetds(1,1)+dble(i2*i2)*dgmetds(2,2)+&
1475 &   dble(i3*i3)*dgmetds(3,3)+&
1476 &   dble(i1*i2)*(dgmetds(1,2)+dgmetds(2,1))+&
1477 &   dble(i1*i3)*(dgmetds(1,3)+dgmetds(3,1))+&
1478 &   dble(i2*i3)*(dgmetds(2,3)+dgmetds(3,2))
1479  end function dgsqds_vl
1480 
1481 end subroutine vlocalstr