TABLE OF CONTENTS


ABINIT/m_kg [ Modules ]

[ Top ] [ Modules ]

NAME

 m_kg

FUNCTION

  Low-level functions to operate of G-vectors.

COPYRIGHT

  Copyright (C) 2008-2018 ABINIT group (DCA, XG, GMR, MT, DRH, AR)
  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_kg
28 
29  use defs_abitypes, only : dataset_type
30  use defs_basis
31  use m_errors
32  use m_abicore
33  use m_errors
34  use m_xmpi
35 
36  use m_fftcore,     only : kpgsph, bound
37  use defs_abitypes, only : MPI_type
38  use m_mpinfo,      only : proc_distrb_cycle
39 
40 
41  implicit none
42 
43  private

ABINIT/mknucdipmom_k [ Functions ]

[ Top ] [ Functions ]

NAME

 mknucdipmom_k

FUNCTION

 compute Hamiltonian in reciprocal space due to array of nuclear
 dipole moments, at a given k point

COPYRIGHT

 Copyright (C) 1998-2018 ABINIT group
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .

INPUTS

 gmet(3,3)=metric for reciprocal space vectors
 kg(3,npw)=reduced planewave coordinates at current k point
 kpt(3)=current k point, reduced coordinates
 natom=number of atoms in cell
 nucdipmom(3,natom)=nuclear dipole moment vectors, at each atom
 npw=number of planewaves
 rprimd(3,3)=real space translation vectors
 ucvol=unit cell volume
 xred(3,natom)=location of atoms in unit cell, in reduced coordinates

OUTPUT

  nucdipmom_k(npw*(npw+1)/2) = nuclear dipole moment Hamiltonian matrix, in
                                 lower diagonal Hermitian packed storage, at current k point

SIDE EFFECTS

NOTES

PARENTS

CHILDREN

SOURCE

1313 subroutine mknucdipmom_k(gmet,kg,kpt,natom,nucdipmom,nucdipmom_k,npw,rprimd,ucvol,xred)
1314 
1315 
1316 !This section has been created automatically by the script Abilint (TD).
1317 !Do not modify the following lines by hand.
1318 #undef ABI_FUNC
1319 #define ABI_FUNC 'mknucdipmom_k'
1320 !End of the abilint section
1321 
1322   implicit none
1323 
1324   !Arguments ------------------------------------
1325   !scalars
1326   integer,intent(in) :: natom,npw
1327   real(dp),intent(in) :: ucvol
1328 
1329   !arrays
1330   integer,intent(in) :: kg(3,npw)
1331   real(dp),intent(in) :: gmet(3,3),kpt(3),nucdipmom(3,natom),rprimd(3,3),xred(3,natom)
1332   complex(dpc),intent(out) :: nucdipmom_k(npw*(npw+1)/2)
1333 
1334   !Local variables-------------------------------
1335   !scalars
1336   integer :: atom_nd_tot,col,iatom,ndp_index,row
1337   real(dp) :: crossfac,dg2,permeability,phasefac
1338   complex(dpc) :: cpermfac,cphasefac
1339   !arrays
1340   integer :: atom_nd(natom)
1341   real(dp) :: cprod(3),cprod_cart(3),dgp_red(3), gpk_red(3)
1342 
1343   ! *************************************************************************
1344   !
1345 
1346   ! magnetic permeability mu_0/four_pi in atomic units
1347   ! this constant is also used in m_pawdij.F90/pawdijnd, if you change it here,
1348   ! change it there also for consistency
1349   permeability=5.325135453D-5
1350   ! will need 4*pi*i*(\mu_0/four\pi)
1351   cpermfac = CMPLX(zero,four_pi*permeability)
1352 
1353   ! make list of atoms with non-zero nuclear magnetic dipoles
1354   atom_nd_tot = 0
1355   do iatom = 1, natom
1356      if(any(abs(nucdipmom(:,iatom))>tol12)) then
1357         atom_nd_tot = atom_nd_tot + 1
1358         atom_nd(atom_nd_tot) = iatom
1359      end if
1360   end do
1361 
1362   ndp_index = 0
1363   do col=1,npw ! enumerate plane waves G
1364      ! form k + G at this k point for current plane wave (this is the ket |k+G> )
1365      ! in reduced coordinates
1366      gpk_red(1)=dble(kg(1,col))+kpt(1)
1367      gpk_red(2)=dble(kg(2,col))+kpt(2)
1368      gpk_red(3)=dble(kg(3,col))+kpt(3)
1369 
1370      do row=col,npw ! enumerate lower diagonal from 1 to G
1371         ! index of the current matrix element, in lower triangular packed storage
1372         ! "packed sequentially, column by column"
1373         ndp_index = ndp_index + 1
1374         nucdipmom_k(ndp_index) = czero
1375 
1376         ! form G-G' = \Delta G at this k pt (this is the bra <k+G'| )
1377         ! in reduced coordinates
1378         dgp_red(1)=dble(kg(1,col)-kg(1,row))
1379         dgp_red(2)=dble(kg(2,col)-kg(2,row))
1380         dgp_red(3)=dble(kg(3,col)-kg(3,row))
1381 
1382         ! compute |\Delta G|^2
1383         ! must use gmet metric because G's are in reduced coords in reciprocal space
1384         dg2 = DOT_PRODUCT(dgp_red,MATMUL(gmet,dgp_red))
1385         ! if \Delta G = 0, Hamiltonian term is zero and move on to next one
1386         if (abs(dg2)<tol8) then
1387            nucdipmom_k(ndp_index)=czero
1388            cycle
1389         end if
1390 
1391         ! compute cross product \Delta G \times (k + G)
1392         ! notice that \Delta G and (k + G) are in reduced coords in reciprocal space
1393         cprod(1) = dgp_red(2)*gpk_red(3) - dgp_red(3)*gpk_red(2)
1394         cprod(2) = dgp_red(3)*gpk_red(1) - dgp_red(1)*gpk_red(3)
1395         cprod(3) = dgp_red(1)*gpk_red(2) - dgp_red(2)*gpk_red(1)
1396 
1397         ! proper cross product must account for reduced coords as follows:
1398         ! gprimd*dgp \times gprimd*gpk = (det gprimd)*(gprimd^{-1,T})*(dgp \times gpk)
1399         ! = rprimd * (dgp \times gpk)/ucvol
1400         ! final vector also includes the division by |\Delta G|^2
1401         cprod_cart = MATMUL(rprimd,cprod)/(ucvol*dg2)
1402 
1403         ! loop over the atoms with non-zero nuclear dipoles
1404         ! phase factors exp(i*\Delta G*I) where I is ion position,
1405         ! might be retrievable from ph1d, need to check
1406         do iatom = 1, atom_nd_tot
1407            phasefac = two_pi*DOT_PRODUCT(dgp_red,xred(:,atom_nd(iatom)))
1408            cphasefac = CMPLX(cos(phasefac),sin(phasefac))
1409            crossfac = DOT_PRODUCT(nucdipmom(:,iatom),cprod_cart)
1410            nucdipmom_k(ndp_index) = nucdipmom_k(ndp_index) + cpermfac*cphasefac*crossfac
1411         end do ! end loop over atoms with nonzero dipoles
1412 
1413      end do ! end loop over G' = G to npw
1414 
1415   end do ! end loop over G = 1 to npw
1416 
1417 end subroutine mknucdipmom_k

ABINIT/mkpwind_k [ Functions ]

[ Top ] [ Functions ]

NAME

 mkpwind_k

FUNCTION

 Make plane wave index at k point for basis at second k point,
 needed to compute overlaps $\langle u_{k,n}|u_{k+b,n}\rangle$
 as appear in Berry phase derived quantities

COPYRIGHT

 Copyright (C) 2003-2017 ABINIT  group
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .
 For the initials of contributors, see ~abinit/doc/developers/contributors.txt.

INPUTS

 dk(3)=real vector difference of ket kpt - bra kpt
 dtset <type(dataset_type)>=all input variables in this dataset
 fnkpt=number of kpts in full BZ
 fkptns=kpts in full BZ
 gmet(3,3)=metric in reciprocal space
 indkk_f2ibz(fnkpt,6)=information on folding from FBZ to IBZ (see initberry or initorbmag)
 ikpt=index of bra k pt
 ikpt1=index of neighbour ket k pt
 kg(3,dtset%mpw*dtset%mkmem)=planewave basis data
 kgindex(dtset%nkpt)= index of kg per kpt
 mpi_enreg=information about MPI parallelization
 npw_k=number of planewaves at k
 symrec(3,3,nsym) = symmetries in reciprocal space in terms of
   reciprocal space primitive translations

OUTPUT

 pwind_k1(dtset%mpw)=output index of ikpt1 basis states refered to ikpt

SIDE EFFECTS

TODO

NOTES

PARENTS

      chern_number

CHILDREN

      kpgsph

SOURCE

1161 subroutine mkpwind_k(dk,dtset,fnkpt,fkptns,gmet,indkk_f2ibz,ikpt,ikpt1,&
1162 & kg,kgindex,mpi_enreg,npw_k,pwind_k1,symrec)
1163 
1164 
1165 !This section has been created automatically by the script Abilint (TD).
1166 !Do not modify the following lines by hand.
1167 #undef ABI_FUNC
1168 #define ABI_FUNC 'mkpwind_k'
1169 !End of the abilint section
1170 
1171   implicit none
1172 
1173   !Arguments ------------------------------------
1174   !scalars
1175   integer,intent(in) :: fnkpt,ikpt,ikpt1,npw_k
1176   type(dataset_type),intent(in) :: dtset
1177   type(MPI_type), intent(inout) :: mpi_enreg
1178 
1179   !arrays
1180   integer,intent(in) :: indkk_f2ibz(fnkpt,6),kg(3,dtset%mpw*dtset%mkmem),kgindex(dtset%nkpt)
1181   integer,intent(in) :: symrec(3,3,dtset%nsym)
1182   integer,intent(out) :: pwind_k1(dtset%mpw)
1183   real(dp),intent(in) :: dk(3),fkptns(3,fnkpt),gmet(3,3)
1184 
1185   !Local variables -------------------------
1186   !scalars
1187   integer :: exchn2n3d,idum1,ikg1,ipw,istwf_k,isym,isym1,jpw,npw_k1
1188   real(dp) :: ecut_eff
1189 
1190   !arrays
1191   integer,allocatable :: kg1_k(:,:)
1192   real(dp) :: dg(3),dum33(3,3),kpt1(3),iadum(3),iadum1(3)
1193 
1194   ! ***********************************************************************
1195 
1196   ABI_ALLOCATE(kg1_k,(3,dtset%mpw))
1197 
1198   ecut_eff = dtset%ecut*(dtset%dilatmx)**2
1199   exchn2n3d = 0 ; istwf_k = 1 ; ikg1 = 0
1200 
1201   ! Build basis sphere of plane waves for the nearest neighbour of the k-point
1202 
1203   kg1_k(:,:) = 0
1204   kpt1(:) = dtset%kptns(:,ikpt1)
1205   call kpgsph(ecut_eff,exchn2n3d,gmet,ikg1,ikpt1,istwf_k,kg1_k,kpt1,1,mpi_enreg,dtset%mpw,npw_k1)
1206 
1207   !
1208   !        Deal with symmetry transformations
1209   !
1210 
1211   !        bra k-point k(b) and IBZ k-point kIBZ(b) related by
1212   !        k(b) = alpha(b) S(b)^t kIBZ(b) + G(b)
1213   !        where alpha(b), S(b) and G(b) are given by indkk_f2ibz
1214   !
1215   !        For the ket k-point:
1216   !        k(k) = alpha(k) S(k)^t kIBZ(k) + G(k) - GBZ(k)
1217   !        where GBZ(k) takes k(k) to the BZ
1218   !
1219 
1220   isym  = indkk_f2ibz(ikpt,2)
1221   isym1 = indkk_f2ibz(ikpt1,2)
1222 
1223   !        Construct transformed G vector that enters the matching condition:
1224   !        alpha(k) S(k)^{t,-1} ( -G(b) - GBZ(k) + G(k) )
1225 
1226   dg(:) = -indkk_f2ibz(ikpt,3:5) &
1227        & - nint(-fkptns(:,ikpt) - dk(:) - tol10 + fkptns(:,ikpt1)) &
1228        & + indkk_f2ibz(ikpt1,3:5)
1229 
1230   iadum(:) = MATMUL(TRANSPOSE(dtset%symrel(:,:,isym1)),dg(:))
1231 
1232   dg(:) = iadum(:)
1233 
1234   !        Construct S(k)^{t,-1} S(b)^{t}
1235 
1236   dum33(:,:) = MATMUL(TRANSPOSE(dtset%symrel(:,:,isym1)),symrec(:,:,isym))
1237 
1238   !        Construct alpha(k) alpha(b)
1239 
1240   pwind_k1(:) = 0
1241   do ipw = 1, npw_k
1242 
1243      !          NOTE: the bra G vector is taken for the sym-related IBZ k point,
1244      !          not for the FBZ k point
1245      iadum(:) = kg(:,kgindex(ikpt) + ipw)
1246 
1247      !          to determine r.l.v. matchings, we transformed the bra vector
1248      !          Rotation
1249      iadum1(:)=0
1250      do idum1=1,3
1251         iadum1(:)=iadum1(:)+dum33(:,idum1)*iadum(idum1)
1252      end do
1253      iadum(:)=iadum1(:)
1254      iadum(:) = iadum(:) + dg(:)
1255 
1256      do jpw = 1, npw_k1
1257         iadum1(1:3) = kg1_k(1:3,jpw)
1258         if ( (iadum(1) == iadum1(1)).and. &
1259              &     (iadum(2) == iadum1(2)).and. &
1260              &     (iadum(3) == iadum1(3)) ) then
1261            pwind_k1(ipw) = jpw
1262            ! write(std_out,'(a,2i4)')'JWZ debug : bg ipw == jpw ',ipw,jpw
1263            exit
1264         end if
1265      end do
1266   end do
1267 
1268   ABI_DEALLOCATE(kg1_k)
1269 
1270 end subroutine mkpwind_k

m_kg/getcut [ Functions ]

[ Top ] [ m_kg ] [ Functions ]

NAME

 getcut

FUNCTION

 For input kpt, fft box dim ngfft(1:3), recip space metric gmet,
 and kinetic energy cutoff ecut (hartree), COMPUTES:
 if iboxcut==0:
   gsqcut: cut-off on G^2 for "large sphere" of radius double that
            of the basis sphere corresponding to ecut
   boxcut: where boxcut == gcut(box)/gcut(sphere).
                 boxcut >=2 for no aliasing.
                 boxcut < 1 is wrong and halts subroutine.
 if iboxcut==1:
   gsqcut: cut-off on G^2 for "large sphere"
            containing the whole fft box
   boxcut: no meaning (zero)

INPUTS

 ecut=kinetic energy cutoff for planewave sphere (hartree)
 gmet(3,3)=reciprocal space metric (bohr^-2)
 iboxcut=0: compute gsqcut and boxcut with boxcut>=1
         1: compute gsqcut for boxcut=1 (sphere_cutoff=box_cutoff)
 iout=unit number for output file
 kpt(3)=input k vector (reduced coordinates--in terms of reciprocal lattice primitive translations)
 ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft

OUTPUT

 boxcut=defined above (dimensionless), ratio of basis sphere
  diameter to fft box length (smallest value)
 gsqcut=Fourier cutoff on G^2 for "large sphere" of radius double
  that of the basis sphere--appropriate for charge density rho(G),
  Hartree potential, and pseudopotentials

NOTES

 2*gcut arises from rho(g)=sum g prime (psi(g primt)*psi(g prime+g))
               where psi(g) is only nonzero for |g| <= gcut).
 ecut (currently in hartree) is proportional to gcut(sphere)**2.

PARENTS

      dfpt_looppert,dfpt_scfcv,m_pawfgr,nonlinear,respfn,scfcv,setup1
      setup_positron

CHILDREN

      bound,wrtout

SOURCE

109 subroutine getcut(boxcut,ecut,gmet,gsqcut,iboxcut,iout,kpt,ngfft)
110 
111 
112 !This section has been created automatically by the script Abilint (TD).
113 !Do not modify the following lines by hand.
114 #undef ABI_FUNC
115 #define ABI_FUNC 'getcut'
116 !End of the abilint section
117 
118  implicit none
119 
120 !Arguments ------------------------------------
121 !scalars
122  integer,intent(in) :: iboxcut,iout
123  real(dp),intent(in) :: ecut
124  real(dp),intent(out) :: boxcut,gsqcut
125 !arrays
126  integer,intent(in) :: ngfft(18)
127  real(dp),intent(in) :: gmet(3,3),kpt(3)
128 
129 !Local variables-------------------------------
130 !scalars
131  integer :: plane
132  real(dp) :: boxsq,cutrad,ecut_pw,effcut,largesq,sphsq
133  character(len=500) :: message
134 !arrays
135  integer :: gbound(3)
136 
137 ! *************************************************************************
138 
139 !This is to treat the case where ecut has not been initialized,
140 !for wavelet computations. The default for ecut is -1.0 , allowed
141 !only for wavelets calculations
142  ecut_pw=ecut
143  if(ecut<-tol8)ecut_pw=ten
144 
145 !gcut(box)**2=boxsq; gcut(sphere)**2=sphsq
146 !get min. d**2 to boundary of fft box:
147 !(gmet sets dimensions: bohr**-2)
148 !ecut(sphere)=0.5*(2 pi)**2 * sphsq:
149  call bound(largesq,boxsq,gbound,gmet,kpt,ngfft,plane)
150  effcut=0.5_dp * (two_pi)**2 * boxsq
151  sphsq=2._dp*ecut_pw/two_pi**2
152 
153  if (iboxcut/=0) then
154    boxcut=10._dp
155    gsqcut=(largesq/sphsq)*(2.0_dp*ecut)/two_pi**2
156 
157    write(message, '(a,a,3f8.4,a,3i4,a,a,f11.3,a,a)' ) ch10,&
158 &   ' getcut: wavevector=',kpt,'  ngfft=',ngfft(1:3),ch10,&
159 &   '         ecut(hartree)=',ecut_pw+tol8,ch10,'=> whole FFT box selected'
160    if(iout/=std_out) then
161      call wrtout(iout,message,'COLL')
162    end if
163    call wrtout(std_out,message,'COLL')
164 
165  else
166 
167 !  Get G^2 cutoff for sphere of double radius of basis sphere
168 !  for selecting G s for rho(G), V_Hartree(G), and V_psp(G)--
169 !  cut off at fft box boundary or double basis sphere radius, whichever
170 !  is smaller.  If boxcut were 2, then relation would be
171 !$ecut_eff = (1/2) * (2 Pi Gsmall)^2 and gsqcut=4*Gsmall^2$.
172    boxcut = sqrt(boxsq/sphsq)
173    cutrad = min(2.0_dp,boxcut)
174    gsqcut = (cutrad**2)*(2.0_dp*ecut_pw)/two_pi**2
175 
176    if(ecut>-tol8)then
177 
178      write(message, '(a,a,3f8.4,a,3i4,a,a,f11.3,3x,a,f10.5)' ) ch10,&
179 &     ' getcut: wavevector=',kpt,'  ngfft=',ngfft(1:3),ch10,&
180 &     '         ecut(hartree)=',ecut+tol8,'=> boxcut(ratio)=',boxcut+tol8
181      if(iout/=std_out) then
182        call wrtout(iout,message,'COLL')
183      end if
184      call wrtout(std_out,message,'COLL')
185 
186      if (boxcut<1.0_dp) then
187        write(message, '(a,a,a,a,a,a,a,a,a,f12.6)' )&
188 &       '  Choice of acell, ngfft, and ecut',ch10,&
189 &       '  ===> basis sphere extends BEYOND fft box !',ch10,&
190 &       '  Recall that boxcut=Gcut(box)/Gcut(sphere)  must be > 1.',ch10,&
191 &       '  Actio: try larger ngfft or smaller ecut.',ch10,&
192 &       '  Note that ecut=effcut/boxcut**2 and effcut=',effcut+tol8
193        if(iout/=std_out) then
194          call wrtout(iout,message,'COLL')
195        end if
196        MSG_ERROR(message)
197      end if
198 
199      if (boxcut>2.2_dp) then
200        write(message, '(a,a,a,a,a,a,a,a,a,a,a,f12.6,a,a)' ) ch10,&
201 &       ' getcut : COMMENT -',ch10,&
202 &       '  Note that boxcut > 2.2 ; recall that',' boxcut=Gcut(box)/Gcut(sphere) = 2',ch10,&
203 &       '  is sufficient for exact treatment of convolution.',ch10,&
204 &       '  Such a large boxcut is a waste : you could raise ecut',ch10,&
205 &       '  e.g. ecut=',effcut*0.25_dp+tol8,' Hartrees makes boxcut=2',ch10
206        if(iout/=std_out) then
207          call wrtout(iout,message,'COLL')
208        end if
209        call wrtout(std_out,message,'COLL')
210      end if
211 
212      if (boxcut<1.5_dp) then
213        write(message, '(a,a,a,a,a,a,a,a,a,a)' ) ch10,&
214 &       ' getcut : WARNING -',ch10,&
215 &       '  Note that boxcut < 1.5; this usually means',ch10,&
216 &       '  that the forces are being fairly strongly affected by','  the smallness of the fft box.',ch10,&
217 &       '  Be sure to test with larger ngfft(1:3) values.',ch10
218        if(iout/=std_out) then
219          call wrtout(iout,message,'COLL')
220        end if
221        call wrtout(std_out,message,'COLL')
222      end if
223 
224    end if
225 
226  end if  ! iboxcut
227 
228 end subroutine getcut

m_kg/getmpw [ Functions ]

[ Top ] [ m_kg ] [ Functions ]

NAME

 getmpw

FUNCTION

 From input ecut, combined with ucvol and gmet, compute recommended mpw
 mpw is the maximum number of plane-waves in the wave-function basis
  for one processor of the WF group

INPUTS

 ecut=plane-wave cutoff energy in Hartrees
 exchn2n3d=if n1, n2 and n3 are exchanged
 gmet(3,3)=reciprocal space metric (bohr**-2).
 istwfk(nkpt)=input option parameter that describes the storage of wfs
 kptns(3,nkpt)=real(dp) array for k points (normalisation is already taken into account)
 mpi_enreg=information about MPI parallelization
 nkpt=integer number of k points in the calculation

OUTPUT

 mpw=maximal number of plane waves over all k points of the processor
  (for one processor of the WF group)

PARENTS

      dfpt_looppert,memory_eval,mpi_setup,scfcv

CHILDREN

      kpgsph,wrtout

SOURCE

261 subroutine getmpw(ecut,exchn2n3d,gmet,istwfk,kptns,mpi_enreg,mpw,nkpt)
262 
263 
264 !This section has been created automatically by the script Abilint (TD).
265 !Do not modify the following lines by hand.
266 #undef ABI_FUNC
267 #define ABI_FUNC 'getmpw'
268 !End of the abilint section
269 
270  implicit none
271 
272 !Arguments ------------------------------------
273 !scalars
274  integer,intent(in) :: exchn2n3d,nkpt
275  integer,intent(out) :: mpw
276  real(dp),intent(in) :: ecut
277  type(MPI_type),intent(inout) :: mpi_enreg
278 !arrays
279  integer,intent(in) :: istwfk(nkpt)
280  real(dp),intent(in) :: gmet(3,3),kptns(3,nkpt)
281 
282 !Local variables-------------------------------
283 !scalars
284  integer :: ikpt,istwf_k,npw
285 ! integer :: npwwrk,pad=50
286 ! real(dp) :: scale=1.3_dp
287  character(len=500) :: message
288 !arrays
289  integer,allocatable :: kg(:,:)
290  real(dp) :: kpoint(3)
291 
292 ! *************************************************************************
293 
294 !An upper bound for mpw, might be obtained as follows
295 !the average number of plane-waves in the cutoff sphere is
296 !npwave = (2*ecut)**(3/2) * ucvol / (6*pi**2)
297 !the upper bound is calculated as
298 !npwwrk = int(scale * npwave) + pad
299 !rescale so an upper bound
300 !npwave=nint(ucvol*(2.0_dp*ecut)**1.5_dp/(6.0_dp*pi**2))
301 !npwwrk=nint(dble(npwave)*scale)+pad
302 
303  ABI_ALLOCATE(kg,(3,100))
304 
305 !set mpw to zero, as needed for only counting in kpgsph
306  mpw = 0
307 
308 !Might be parallelized over k points ? !
309  do ikpt = 1,nkpt
310 !  Do computation of G sphere, returning npw
311    kpoint(:)=kptns(:,ikpt)
312    istwf_k=istwfk(ikpt)
313    call kpgsph(ecut,exchn2n3d,gmet,0,ikpt,istwf_k,kg,kpoint,0,mpi_enreg,0,npw)
314    mpw = max(npw,mpw)
315  end do
316 
317  write(message,'(a,i0)') ' getmpw: optimal value of mpw= ',mpw
318  call wrtout(std_out,message,'COLL')
319 
320  ABI_DEALLOCATE(kg)
321 
322 end subroutine getmpw

m_kg/getph [ Functions ]

[ Top ] [ m_kg ] [ Functions ]

NAME

 getph

FUNCTION

 Compute three factors of one-dimensional structure factor phase
 for input atomic coordinates, for all planewaves which fit in fft box.
 The storage of these atomic factors is made according to the
 values provided by the index table atindx. This will save time in nonlop.

INPUTS

  atindx(natom)=index table for atoms (see gstate.f)
  natom=number of atoms in cell.
  n1,n2,n3=dimensions of fft box (ngfft(3)).
  xred(3,natom)=reduced atomic coordinates.

OUTPUT

  ph1d(2,(2*n1+1)*natom+(2*n2+1)*natom+(2*n3+1)*natom)=exp(2Pi i G.xred) for
   integer vector G with components ranging from -nj <= G <= nj.
   Real and imag given in usual Fortran convention.

PARENTS

      afterscfloop,bethe_salpeter,cg_rotate,dfpt_looppert,dfptnl_loop
      extrapwf,gstate,m_cut3d,m_fock,m_gkk,m_hamiltonian,m_phgamma,m_phpi
      m_sigmaph,m_wfd,partial_dos_fractions,prcref,prcref_PMA,respfn,scfcv
      screening,sigma,update_e_field_vars,wfconv

CHILDREN

SOURCE

821 subroutine getph(atindx,natom,n1,n2,n3,ph1d,xred)
822 
823 
824 !This section has been created automatically by the script Abilint (TD).
825 !Do not modify the following lines by hand.
826 #undef ABI_FUNC
827 #define ABI_FUNC 'getph'
828 !End of the abilint section
829 
830  implicit none
831 
832 !Arguments ------------------------------------
833 !scalars
834  integer,intent(in) :: n1,n2,n3,natom
835 !arrays
836  integer,intent(in) :: atindx(natom)
837  real(dp),intent(in) :: xred(3,natom)
838  real(dp),intent(out) :: ph1d(:,:)
839 
840 !Local variables-------------------------------
841 !scalars
842  integer,parameter :: im=2,re=1
843  integer :: i1,i2,i3,ia,ii,ph1d_size1,ph1d_size2,ph1d_sizemin
844  !character(len=500) :: msg
845  real(dp) :: arg
846 
847 ! *************************************************************************
848 
849  ph1d_size1=size(ph1d,1);ph1d_size2=size(ph1d,2)
850  ph1d_sizemin=(2*n1+1+2*n2+1+2*n3+1)*natom
851  if (ph1d_size1/=2.or.ph1d_size2<ph1d_sizemin) then
852    MSG_BUG('Wrong ph1d sizes!')
853  end if
854 
855  do ia=1,natom
856 
857 !  Store the phase factor of atom number ia in place atindx(ia)
858    i1=(atindx(ia)-1)*(2*n1+1)
859    i2=(atindx(ia)-1)*(2*n2+1)+natom*(2*n1+1)
860    i3=(atindx(ia)-1)*(2*n3+1)+natom*(2*n1+1+2*n2+1)
861 
862    do ii=1,2*n1+1
863      arg=two_pi*dble(ii-1-n1)*xred(1,ia)
864      ph1d(re,ii+i1)=dcos(arg)
865      ph1d(im,ii+i1)=dsin(arg)
866    end do
867 
868    do ii=1,2*n2+1
869      arg=two_pi*dble(ii-1-n2)*xred(2,ia)
870      ph1d(re,ii+i2)=dcos(arg)
871      ph1d(im,ii+i2)=dsin(arg)
872    end do
873 
874    do ii=1,2*n3+1
875      arg=two_pi*dble(ii-1-n3)*xred(3,ia)
876      ph1d(re,ii+i3)=dcos(arg)
877      ph1d(im,ii+i3)=dsin(arg)
878    end do
879 
880  end do
881 
882 !This is to avoid uninitialized ph1d values
883  if (ph1d_sizemin<ph1d_size2) then
884    ph1d(:,ph1d_sizemin+1:ph1d_size2)=zero
885  end if
886 
887 end subroutine getph

m_kg/kpgio [ Functions ]

[ Top ] [ m_kg ] [ Functions ]

NAME

 kpgio

FUNCTION

 Do initialization of kg data.

INPUTS

  ecut=kinetic energy planewave cutoff (hartree)
  exchn2n3d=if 1, n2 and n3 are exchanged
  gmet(3,3)=reciprocal space metric (bohr^-2)
  istwfk(nkpt)=input option parameter that describes the storage of wfs
  kptns(3,nkpt)=reduced coords of k points
  mkmem =number of k points treated by this node.
  character(len=4) : mode_paral=either 'COLL' or 'PERS', tells whether
   the loop over k points must be done by all processors or not,
   in case of parallel execution.
  mpi_enreg=information about MPI parallelization
  mpw=maximum number of planewaves as dimensioned in calling routine
  nband(nkpt*nsppol)=number of bands at each k point
  nkpt=number of k points
  nsppol=1 for unpolarized, 2 for polarized

OUTPUT

  npwarr(nkpt)=array holding npw for each k point, taking into account
   the effect of istwfk, and the spreading over processors
  npwtot(nkpt)=array holding the total number of plane waves for each k point,
  kg(3,mpw*mkmem)=dimensionless coords of G vecs in basis sphere at k point

NOTES

 Note that in case of band parallelism, the number of spin-up
 and spin-down bands must be equal at each k points

PARENTS

      dfpt_looppert,dfptff_initberry,gstate,initmv,inwffil,m_cut3d,nonlinear
      respfn,scfcv

CHILDREN

      kpgsph,wrtout,xmpi_sum

SOURCE

536 subroutine kpgio(ecut,exchn2n3d,gmet,istwfk,kg,kptns,mkmem,nband,nkpt,&
537 & mode_paral,mpi_enreg,mpw,npwarr,npwtot,nsppol)
538 
539 
540 !This section has been created automatically by the script Abilint (TD).
541 !Do not modify the following lines by hand.
542 #undef ABI_FUNC
543 #define ABI_FUNC 'kpgio'
544 !End of the abilint section
545 
546  implicit none
547 
548 !Arguments ------------------------------------
549 !scalars
550  integer,intent(in) :: exchn2n3d,mkmem,mpw,nkpt,nsppol
551  real(dp),intent(in) :: ecut
552  character(len=4),intent(in) :: mode_paral
553  type(MPI_type),intent(inout) :: mpi_enreg
554 !arrays
555  integer,intent(in) :: istwfk(nkpt),nband(nkpt*nsppol)
556  integer,intent(out) :: kg(3,mpw*mkmem),npwarr(nkpt),npwtot(nkpt)
557  real(dp),intent(in) :: gmet(3,3),kptns(3,nkpt)
558 
559 !Local variables-------------------------------
560 !scalars
561  integer :: ierr,ikg,ikpt,istwf_k,me,nband_down,nband_k,npw1
562  logical :: test_npw
563  character(len=500) :: message
564 !arrays
565  real(dp) :: kpoint(3)
566 
567 ! *************************************************************************
568 
569 !DEBUG
570 !write(std_out,*)' kpgio : enter '
571 !ENDDEBUG
572 
573 !Define me
574  me=mpi_enreg%me_kpt
575 
576  if((mpi_enreg%paralbd==1) .and. (mode_paral=='PERS')) then
577    if(nsppol==2)then
578      do ikpt=1,nkpt
579        nband_k=nband(ikpt)
580        nband_down=nband(ikpt+nkpt)
581        if(nband_k/=nband_down)then
582          write(message,'(a,a,a,a,a,a,a,a,i4,a,i4,a,a,a,i4,a,a,a)')ch10,&
583 &         ' kpgio : ERROR -',ch10,&
584 &         '  Band parallel case, one must have same number',ch10,&
585 &         '  of spin up and spin down bands, but input is :',ch10,&
586 &         '  nband(up)=',nband_k,', nband(down)=',nband_down,',',ch10,&
587 &         '  for ikpt=',ikpt,'.',ch10,&
588 &         '  Action: correct nband in your input file.'
589 !        MG: Tests v3(10,11,17) and v6(67) fail if this test is enabled
590 !        call wrtout(std_out,message,mode_paral)
591        end if
592      end do
593    end if
594  end if
595 
596  npwarr(:)=0
597  npwtot(:)=0
598 
599  kg=0
600  ikg=0
601 !Find (k+G) sphere for each k.
602 
603  do ikpt=1,nkpt
604 
605    nband_k = nband(ikpt)
606 
607    if(mode_paral=='PERS')then
608      if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,-1,me)) cycle
609    end if
610 
611    kpoint(:)=kptns(:,ikpt)
612    istwf_k=istwfk(ikpt)
613    call kpgsph(ecut,exchn2n3d,gmet,ikg,ikpt,istwf_k,kg,kpoint,mkmem,mpi_enreg,mpw,npw1)
614 
615    test_npw=.true.
616    if (xmpi_paral==1)then
617      if (mode_paral=='PERS')then
618        test_npw=(minval(mpi_enreg%proc_distrb(ikpt,1:nband_k,1:nsppol))==me)
619      end if
620    end if
621    if (test_npw) npwarr(ikpt)=npw1
622 
623 !  Make sure npw < nband never happens:
624 !  if (npw1<nband(ikpt)) then
625 !  write(message, '(a,a,a,a,i5,a,3f8.4,a,a,i10,a,i10,a,a,a,a)' )ch10,&
626 !  &   ' kpgio : ERROR -',ch10,&
627 !  &   '  At k point number',ikpt,' k=',(kptns(mu,ikpt),mu=1,3),ch10,&
628 !  &   '  npw=',npw1,' < nband=',nband(ikpt),ch10,&
629 !  &   '  Indicates not enough planewaves for desired number of bands.',ch10,&
630 !  &   '  Action: change either ecut or nband in input file.'
631 !  MSG_ERROR(message)
632 !  end if
633 
634 !  Find boundary of G sphere for efficient zero padding,
635 !    Shift to next section of each array kg
636    ikg=ikg+npw1
637  end do !  End of the loop over k points
638 
639  if(mode_paral == 'PERS') then
640    call xmpi_sum(npwarr,mpi_enreg%comm_kpt,ierr)
641  end if
642 
643  if (mpi_enreg%nproc>1) then
644    call wrtout(std_out,' kpgio: loop on k-points done in parallel','COLL')
645  end if
646 
647 !XG030513 MPIWF : now, one should sum npwarr over all processors
648 !of the WF group, to get npwtot (to be spread on all procs of the WF group
649  npwtot(:)=npwarr(:)
650 
651 !Taking into account istwfk
652  do ikpt=1,nkpt
653    if(istwfk(ikpt)>1)then
654      if(istwfk(ikpt)==2)then
655        npwtot(ikpt)=2*npwtot(ikpt)-1
656      else
657        npwtot(ikpt)=2*npwtot(ikpt)
658      end if
659    end if
660  end do
661 
662 !DEBUG
663 !write(std_out,*)' kpgio : exit '
664 !ENDDEBUG
665 
666 end subroutine kpgio

m_kg/kpgstr [ Functions ]

[ Top ] [ m_kg ] [ Functions ]

NAME

 kpgstr

FUNCTION

 Compute elements of the derivative the kinetic energy operator in reciprocal
 space at given k point wrt a single cartesian strain component

INPUTS

  ecut=cut-off energy for plane wave basis sphere (Ha)
  ecutsm=smearing energy for plane wave kinetic energy (Ha)
  effmass_free=effective mass for electrons (1. in common case)
  gmet(3,3) = reciprocal lattice metric tensor (Bohr**-2)
  gprimd(3,3)=reciprocal space dimensional primitive translations
  istr=1,...6 specifies cartesian strain component 11,22,33,32,31,21
  kg(3,npw) = integer coordinates of planewaves in basis sphere.
  kpt(3)    = reduced coordinates of k point
  npw       = number of plane waves at kpt.

OUTPUT

  dkinpw(npw)=d/deps(istr) ( (1/2)*(2 pi)**2 * (k+G)**2 )

NOTES

  Src_6response/kpg3.f

PARENTS

      dfpt_nsteltwf,dfpt_nstpaw,dfpt_rhofermi,getgh1c

CHILDREN

SOURCE

 922 subroutine kpgstr(dkinpw,ecut,ecutsm,effmass_free,gmet,gprimd,istr,kg,kpt,npw)
 923 
 924 
 925 !This section has been created automatically by the script Abilint (TD).
 926 !Do not modify the following lines by hand.
 927 #undef ABI_FUNC
 928 #define ABI_FUNC 'kpgstr'
 929 !End of the abilint section
 930 
 931  implicit none
 932 
 933 !Arguments -------------------------------
 934 !scalars
 935  integer,intent(in) :: istr,npw
 936  real(dp),intent(in) :: ecut,ecutsm,effmass_free
 937 !arrays
 938  integer,intent(in) :: kg(3,npw)
 939  real(dp),intent(in) :: gmet(3,3),gprimd(3,3),kpt(3)
 940  real(dp),intent(out) :: dkinpw(npw)
 941 
 942 !Local variables -------------------------
 943 !scalars
 944  integer :: ig,ii,ka,kb
 945  real(dp) :: dfsm,dkinetic,dkpg2,ecutsm_inv,fsm,gpk1,gpk2,gpk3,htpisq
 946 ! real(dp) :: d2fsm ! used in commented section below
 947  real(dp) :: kpg2,xx
 948  character(len=500) :: message
 949 !arrays
 950  integer,save :: idx(12)=(/1,1,2,2,3,3,3,2,3,1,2,1/)
 951  real(dp) :: dgmetds(3,3)
 952 
 953 ! *********************************************************************
 954 
 955 !htpisq is (1/2) (2 Pi) **2:
 956  htpisq=0.5_dp*(two_pi)**2
 957 
 958  ecutsm_inv=0.0_dp
 959  if(ecutsm>1.0d-20)ecutsm_inv=1/ecutsm
 960 
 961 !Compute derivative of metric tensor wrt strain component istr
 962  if(istr<1 .or. istr>6)then
 963    write(message, '(a,i10,a,a,a)' )&
 964 &   'Input istr=',istr,' not allowed.',ch10,&
 965 &   'Possible values are 1,2,3,4,5,6 only.'
 966    MSG_BUG(message)
 967  end if
 968 
 969  ka=idx(2*istr-1);kb=idx(2*istr)
 970  do ii = 1,3
 971    dgmetds(:,ii)=-(gprimd(ka,:)*gprimd(kb,ii)+gprimd(kb,:)*gprimd(ka,ii))
 972  end do
 973 !For historical reasons:
 974  dgmetds(:,:)=0.5_dp*dgmetds(:,:)
 975 
 976  do ig=1,npw
 977    gpk1=dble(kg(1,ig))+kpt(1)
 978    gpk2=dble(kg(2,ig))+kpt(2)
 979    gpk3=dble(kg(3,ig))+kpt(3)
 980    kpg2=htpisq*&
 981 &   ( gmet(1,1)*gpk1**2+         &
 982 &   gmet(2,2)*gpk2**2+         &
 983 &   gmet(3,3)*gpk3**2          &
 984 &   +2.0_dp*(gpk1*gmet(1,2)*gpk2+  &
 985 &   gpk1*gmet(1,3)*gpk3+  &
 986 &   gpk2*gmet(2,3)*gpk3 )  )
 987    dkpg2=htpisq*2.0_dp*&
 988 &   (gpk1*(dgmetds(1,1)*gpk1+dgmetds(1,2)*gpk2+dgmetds(1,3)*gpk3)+  &
 989 &   gpk2*(dgmetds(2,1)*gpk1+dgmetds(2,2)*gpk2+dgmetds(2,3)*gpk3)+  &
 990 &   gpk3*(dgmetds(3,1)*gpk1+dgmetds(3,2)*gpk2+dgmetds(3,3)*gpk3) )
 991    dkinetic=dkpg2
 992    if(kpg2>ecut-ecutsm)then
 993      if(kpg2>ecut-tol12)then
 994 !      The wavefunction has been filtered : no derivative
 995        dkinetic=0.0_dp
 996      else
 997        xx=(ecut-kpg2)*ecutsm_inv
 998 !      This kinetic cutoff smoothing function and its xx derivatives
 999 !      were produced with Mathematica and the fortran code has been
1000 !      numerically checked against Mathematica.
1001        fsm=1.0_dp/(xx**2*(3+xx*(1+xx*(-6+3*xx))))
1002        dfsm=-3.0_dp*(-1+xx)**2*xx*(2+5*xx)*fsm**2
1003 !      d2fsm=6.0_dp*xx**2*(9+xx*(8+xx*(-52+xx*(-3+xx*(137+xx*&
1004 !      &                        (-144+45*xx))))))*fsm**3
1005        dkinetic=dkpg2*(fsm-ecutsm_inv*kpg2*dfsm)
1006      end if
1007    end if
1008    dkinpw(ig)=dkinetic/effmass_free
1009  end do
1010 
1011 end subroutine kpgstr

m_kg/mkkin [ Functions ]

[ Top ] [ m_kg ] [ Functions ]

NAME

 mkkin

FUNCTION

 compute elements of kinetic energy operator in reciprocal space at a given k point

INPUTS

  ecut=cut-off energy for plane wave basis sphere (Ha)
  ecutsm=smearing energy for plane wave kinetic energy (Ha)
  effmass_free=effective mass for electrons (1. in common case)
  gmet(3,3)=reciprocal lattice metric tensor ($\textrm{Bohr}^{-2}$)
  idir1 = 1st direction of the derivative (if 1 <= idir1 <= 3, not used otherwise)
  idir2 = 2st direction of the derivative (if 1 <= idir1,idir2 <= 3, not used otherwise))
  kg(3,npw)=integer coordinates of planewaves in basis sphere.
  kpt(3)=reduced coordinates of k point
  npw=number of plane waves at kpt.

OUTPUT

  kinpw(npw)=(modified) kinetic energy (or derivative) for each plane wave (Hartree)

NOTES

 Usually, the kinetic energy expression is $(1/2) (2 \pi)^2 (k+G)^2 $
 However, the present implementation allows for a modification
 of this kinetic energy, in order to obtain smooth total energy
 curves with respect to the cut-off energy or the cell size and shape.
 Thus the usual expression is kept if it is lower then ecut-ecutsm,
 zero is returned beyond ecut, and in between, the kinetic
 energy is DIVIDED by a smearing factor (to make it infinite at the
 cut-off energy). The smearing factor is $x^2 (3-2x)$, where
 x = (ecut- unmodified energy)/ecutsm.
 This smearing factor is also used to derived a modified kinetic
 contribution to stress, in another routine (forstrnps.f)

 Also, in order to break slightly the symmetry between axes, that causes
 sometimes a degeneracy of eigenvalues and do not allow to obtain
 the same results on different machines, there is a modification
 by one part over 1.0e12 of the metric tensor elements (1,1) and (3,3)

PARENTS

      calc_vhxc_me,d2frnl,dfpt_nsteltwf,dfpt_nstpaw,dfpt_rhofermi,energy
      getgh1c,ks_ddiago,m_io_kss,m_vkbr,mkffnl,vtorho

CHILDREN

SOURCE

372 subroutine mkkin (ecut,ecutsm,effmass_free,gmet,kg,kinpw,kpt,npw,idir1,idir2)
373 
374 
375 !This section has been created automatically by the script Abilint (TD).
376 !Do not modify the following lines by hand.
377 #undef ABI_FUNC
378 #define ABI_FUNC 'mkkin'
379 !End of the abilint section
380 
381  implicit none
382 
383 !Arguments ------------------------------------
384 !scalars
385  integer,intent(in) :: npw
386  integer,intent(in) :: idir1,idir2
387  real(dp),intent(in) :: ecut,ecutsm,effmass_free
388 !arrays
389  integer,intent(in) :: kg(3,npw)
390  real(dp),intent(in) :: gmet(3,3),kpt(3)
391  real(dp),intent(out) :: kinpw(npw)
392 
393 !Local variables-------------------------------
394 !scalars
395  integer :: ig,order
396  real(dp),parameter :: break_symm=1.0d-11
397  real(dp) :: ecutsm_inv,fsm,gpk1,gpk2,gpk3,htpisq,kinetic,kpg2,dkpg2,xx
398  real(dp) :: d1kpg2,d2kpg2,ddfsm, dfsm
399 !arrays
400  real(dp) :: gmet_break(3,3)
401 
402 ! *************************************************************************
403 !
404 !htpisq is (1/2) (2 Pi) **2:
405  htpisq=0.5_dp*(two_pi)**2
406 
407  ecutsm_inv=0.0_dp
408  if(ecutsm>1.0d-20)ecutsm_inv=1/ecutsm
409 
410  gmet_break(:,:)=gmet(:,:)
411  gmet_break(1,1)=(1.0_dp+break_symm)*gmet(1,1)
412  gmet_break(3,3)=(1.0_dp-break_symm)*gmet(3,3)
413 
414  order=0 ! Compute the kinetic operator
415  if (idir1>0.and.idir1<4) then
416    order=1 ! Compute the 1st derivative of the kinetic operator
417    if (idir2>0.and.idir2<4) then
418      order=2 ! Compute the 2nd derivative of the kinetic operator
419    end if
420  end if
421 
422 !$OMP PARALLEL DO PRIVATE(dkpg2,d1kpg2,d2kpg2,gpk1,gpk2,gpk3,ig,kinetic,kpg2,xx,fsm,dfsm,ddfsm) &
423 !$OMP SHARED(kinpw,ecut,ecutsm,ecutsm_inv) &
424 !$OMP SHARED(gmet_break,htpisq,idir1,idir2,kg,kpt,npw)
425  do ig=1,npw
426    gpk1=dble(kg(1,ig))+kpt(1)
427    gpk2=dble(kg(2,ig))+kpt(2)
428    gpk3=dble(kg(3,ig))+kpt(3)
429    kpg2=htpisq*&
430 &   ( gmet_break(1,1)*gpk1**2+         &
431 &   gmet_break(2,2)*gpk2**2+         &
432 &   gmet_break(3,3)*gpk3**2          &
433 &   +2.0_dp*(gpk1*gmet_break(1,2)*gpk2+  &
434 &   gpk1*gmet_break(1,3)*gpk3+  &
435 &   gpk2*gmet_break(2,3)*gpk3 )  )
436    select case (order)
437    case(0)
438      kinetic=kpg2
439    case(1)
440      dkpg2=htpisq*2.0_dp*&
441 &     (gmet_break(idir1,1)*gpk1+gmet_break(idir1,2)*gpk2+gmet_break(idir1,3)*gpk3)
442      kinetic=dkpg2
443    case(2)
444      dkpg2=htpisq*2.0_dp*gmet_break(idir1,idir2)
445      kinetic=dkpg2
446    end select
447    if(kpg2>ecut-ecutsm)then
448      if(kpg2>ecut-tol12)then
449        if(order==0) then
450 !        Will filter the wavefunction, based on this value, in cgwf.f, getghc.f and precon.f
451          kinetic=huge(0.0_dp)*1.d-10
452        else
453 !        The wavefunction has been filtered : no derivative
454          kinetic=0
455        end if
456      else
457        if(order==0) then
458          xx=max( (ecut-kpg2)*ecutsm_inv , 1.0d-20)
459        else
460          xx=(ecut-kpg2)*ecutsm_inv
461        end if
462        if(order==2) then
463          d1kpg2=htpisq*2.0_dp*&
464 &         (gmet_break(idir1,1)*gpk1+gmet_break(idir1,2)*gpk2+gmet_break(idir1,3)*gpk3)
465          d2kpg2=htpisq*2.0_dp*&
466 &         (gmet_break(idir2,1)*gpk1+gmet_break(idir2,2)*gpk2+gmet_break(idir2,3)*gpk3)
467        end if
468 !      This kinetic cutoff smoothing function and its xx derivatives
469 !      were produced with Mathematica and the fortran code has been
470 !      numerically checked against Mathematica.
471        fsm=1.0_dp/(xx**2*(3+xx*(1+xx*(-6+3*xx))))
472        if(order>0) dfsm=-3.0_dp*(-1+xx)**2*xx*(2+5*xx)*fsm**2
473        if(order>1) ddfsm=6.0_dp*xx**2*(9+xx*(8+xx*(-52+xx*(-3+xx*(137+xx*(-144+45*xx))))))*fsm**3
474        select case (order)
475        case(0)
476          kinetic=kpg2*fsm
477        case(1)
478          kinetic=dkpg2*(fsm-ecutsm_inv*kpg2*dfsm)
479        case(2)
480          kinetic=dkpg2*fsm&
481 &         -2.0_dp*d1kpg2*dfsm*ecutsm_inv*d2kpg2&
482 &         +kpg2*ddfsm*(ecutsm_inv**2)*d1kpg2*d2kpg2&
483 &         -kpg2*dfsm*ecutsm_inv*dkpg2
484        end select
485      end if
486    end if
487    kinpw(ig)=kinetic/effmass_free
488  end do
489 !$OMP END PARALLEL DO
490 
491 end subroutine mkkin

m_kg/mkkpg [ Functions ]

[ Top ] [ m_kg ] [ Functions ]

NAME

 mkkpg

FUNCTION

 Compute all (k+G) vectors (dp, in reduced coordinates) for given k point,
 from integer coordinates of G vectors.
 Eventually compute related data.

INPUTS

  kg(3,npw)=integer coords of planewaves in basis sphere
  kpt(3)=k point in terms of recip. translations
  nkpg=second dimension of array kpg
  npw=number of plane waves in reciprocal space

OUTPUT

  kpg(npw,3)= (k+G) components
  === if nkpg==9 ===
    kpg(npw,4:9)= [(k+G)_a].[(k+G)_b] quantities

PARENTS

      ctocprj,d2frnl,debug_tools,dfpt_nstpaw,dfpt_nstwf,dfpt_rhofermi
      dfptnl_resp,fock2ACE,forstrnps,getcprj,getgh1c,ks_ddiago,m_io_kss
      m_shirley,m_wfd,nonlop_test,nonlop_ylm,prep_bandfft_tabs,vtorho
      wfd_vnlpsi

CHILDREN

SOURCE

1045 subroutine mkkpg(kg,kpg,kpt,nkpg,npw)
1046 
1047 
1048 !This section has been created automatically by the script Abilint (TD).
1049 !Do not modify the following lines by hand.
1050 #undef ABI_FUNC
1051 #define ABI_FUNC 'mkkpg'
1052 !End of the abilint section
1053 
1054  implicit none
1055 
1056 !Arguments ------------------------------------
1057 !scalars
1058  integer,intent(in) :: nkpg,npw
1059 !arrays
1060  integer,intent(in) :: kg(3,npw)
1061  real(dp),intent(in) :: kpt(3)
1062  real(dp),intent(out) :: kpg(npw,nkpg)
1063 
1064 !Local variables-------------------------------
1065 !scalars
1066  integer :: ipw,mu,mua,mub
1067  character(len=500) :: message
1068 !arrays
1069  integer,parameter :: alpha(6)=(/1,2,3,3,3,2/),beta(6)=(/1,2,3,2,1,1/)
1070 
1071 ! *************************************************************************
1072 
1073  DBG_ENTER("COLL")
1074 
1075  if (nkpg==0) return
1076 
1077 !-- Test nkpg --
1078  if (nkpg/=3.and.nkpg/=9) then
1079    write(message, '(a,i0)' )' Bad value for nkpg !',nkpg
1080    MSG_BUG(message)
1081  end if
1082 
1083 !-- Compute (k+G) --
1084 !$OMP PARALLEL DO COLLAPSE(2) &
1085 !$OMP PRIVATE(mu,ipw)
1086  do ipw=1,npw
1087    do mu=1,3
1088      kpg(ipw,mu)=kpt(mu)+dble(kg(mu,ipw))
1089    end do
1090  end do
1091 !$OMP END PARALLEL DO
1092 
1093 !-- Compute [(k+G)_a].[(k+G)_b] --
1094  if (nkpg==9) then
1095 !$OMP PARALLEL DO COLLAPSE(2) &
1096 !$OMP PRIVATE(ipw,mu,mua,mub)
1097    do ipw=1,npw
1098      do mu=4,9
1099        mua=alpha(mu-3);mub=beta(mu-3)
1100        kpg(ipw,mu)=kpg(ipw,mua)*kpg(ipw,mub)
1101      end do
1102    end do
1103 !$OMP END PARALLEL DO
1104  end if
1105 
1106  DBG_EXIT("COLL")
1107 
1108 end subroutine mkkpg

m_kg/ph1d3d [ Functions ]

[ Top ] [ m_kg ] [ Functions ]

NAME

 ph1d3d

FUNCTION

 Compute the three-dimensional phase factor $e^{i 2 \pi (k+G) cdot xred}$
 from the three one-dimensional factors, the k point coordinates,
 and the atom positions, for all planewaves which fit in the fft box.

INPUTS

  iatom, jatom= bounds of atom indices in ph1d for which ph3d has to be computed
  kg_k(3,npw_k)=reduced planewave coordinates.
  matblk= dimension of ph3d
  natom= dimension of ph1d
  npw=number of plane waves
  n1,n2,n3=dimensions of fft box (ngfft(3)).
  phkxred(2,natom)=phase factors exp(2 pi k.xred)
  ph1d(2,(2*n1+1)*natom+(2*n2+1)*natom+(2*n3+1)*natom)=exp(2Pi i G xred) for
   vectors (Gx,0,0), (0,Gy,0) and (0,0,Gz)
   with components ranging from -nj <= Gj <= nj

OUTPUT

  ph3d(2,npw_k,matblk)=$e^{2 i \pi (k+G) cdot xred}$ for vectors (Gx,Gy,Gz),
   and for atoms in the range iatom to jatom with respect to ph1d

PARENTS

      cg_rotate,ctocprj,getcprj,m_cut3d,m_epjdos,m_fock,m_hamiltonian,m_wfd
      nonlop_pl,nonlop_ylm,partial_dos_fractions,suscep_stat,wfconv

CHILDREN

SOURCE

702 subroutine ph1d3d(iatom,jatom,kg_k,matblk,natom,npw_k,n1,n2,n3,phkxred,ph1d,ph3d)
703 
704 
705 !This section has been created automatically by the script Abilint (TD).
706 !Do not modify the following lines by hand.
707 #undef ABI_FUNC
708 #define ABI_FUNC 'ph1d3d'
709 !End of the abilint section
710 
711  implicit none
712 
713 !Arguments ------------------------------------
714 !scalars
715  integer,intent(in) :: iatom,jatom,matblk,n1,n2,n3,natom,npw_k
716 !arrays
717  integer,intent(in) :: kg_k(3,npw_k)
718  real(dp),intent(in) :: ph1d(2,(2*n1+1+2*n2+1+2*n3+1)*natom)
719  real(dp),intent(in) :: phkxred(2,natom)
720  real(dp),intent(out) :: ph3d(2,npw_k,matblk)
721 
722 !Local variables-------------------------------
723 !scalars
724  integer :: i1,ia,iatblk,ig,shift1,shift2,shift3
725  real(dp) :: ph12i,ph12r,ph1i,ph1r,ph2i,ph2r,ph3i,ph3r,phkxi,phkxr
726  character(len=500) :: message
727 !arrays
728  real(dp),allocatable :: ph1kxred(:,:)
729 
730 ! *************************************************************************
731 
732  if(matblk-1 < jatom-iatom)then
733    write(message,'(a,a,a,a,a,i0,a,a,i0,a,i0,a)')&
734 &   'Input natom-1 must be larger or equal to jatom-iatom,',ch10,&
735 &   'while their value is : ',ch10,&
736 &   'natom-1 = ',natom-1,ch10,&
737 &   'jatom=',jatom,', iatom=',iatom,'.'
738    MSG_BUG(message)
739  end if
740 
741  ABI_ALLOCATE(ph1kxred,(2,-n1:n1))
742 
743 !ia runs from iatom to jatom
744  do ia=iatom,jatom
745 
746 !  iatblk runs from 1 to matblk
747    iatblk=ia-iatom+1
748 !write(87,*) iatblk
749    shift1=1+n1+(ia-1)*(2*n1+1)
750    shift2=1+n2+(ia-1)*(2*n2+1)+natom*(2*n1+1)
751    shift3=1+n3+(ia-1)*(2*n3+1)+natom*(2*n1+1+2*n2+1)
752 !  Compute product of phkxred by phase for the first component of G vector
753    phkxr=phkxred(1,ia)
754    phkxi=phkxred(2,ia)
755 !  DEBUG (needed to compare with version prior to 2.0)
756 !  phkxr=1.0d0
757 !  phkxi=0.0d0
758 !  ENDDEBUG
759    do i1=-n1,n1
760      ph1kxred(1,i1)=ph1d(1,i1+shift1)*phkxr-ph1d(2,i1+shift1)*phkxi
761      ph1kxred(2,i1)=ph1d(2,i1+shift1)*phkxr+ph1d(1,i1+shift1)*phkxi
762    end do
763 
764 !  Compute tri-dimensional phase factor
765 !$OMP PARALLEL DO PRIVATE(ig,ph1r,ph1i,ph2r,ph2i,ph3r,ph3i,ph12r,ph12i)
766    do ig=1,npw_k
767      ph1r=ph1kxred(1,kg_k(1,ig))
768      ph1i=ph1kxred(2,kg_k(1,ig))
769      ph2r=ph1d(1,kg_k(2,ig)+shift2)
770      ph2i=ph1d(2,kg_k(2,ig)+shift2)
771      ph3r=ph1d(1,kg_k(3,ig)+shift3)
772      ph3i=ph1d(2,kg_k(3,ig)+shift3)
773      ph12r=ph1r*ph2r-ph1i*ph2i
774      ph12i=ph1r*ph2i+ph1i*ph2r
775 !if(ig==487) then
776 !write(87,*)iatblk,ph3d(1,ig,iatblk),ph12r,ph3r,ph12i,ph3i
777 !endif
778      ph3d(1,ig,iatblk)=ph12r*ph3r-ph12i*ph3i
779      ph3d(2,ig,iatblk)=ph12r*ph3i+ph12i*ph3r
780    end do
781 !$OMP END PARALLEL DO
782  end do
783 !write(87,*)ph3d(1,487,8)
784  ABI_DEALLOCATE(ph1kxred)
785 
786 end subroutine ph1d3d