TABLE OF CONTENTS


ABINIT/calc_efg [ Functions ]

[ Top ] [ Functions ]

NAME

 calc_efg

FUNCTION

 calculation and output of electric field gradient tensor at each atomic site

COPYRIGHT

 Copyright (C) 2005-2018 ABINIT group (JZ,MT)
 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

  mpi_atmtab(:)=--optional-- indexes of the atoms treated by current proc
  comm_atom=--optional-- MPI communicator over atoms
  mpi_enreg=information about MPI parallelization
  my_natom=number of atoms treated by current processor
  natom=number of atoms in cell.
  nfft=number of points on fft grid
  ngfft(18)=details of fft
  nspden=number of spin densities
  nsym=number of symmetries in space group
  ntypat=number of atom types
  paral_kgb
  ptcharge(ntypat)=user input charges on atoms to make simple point charge calc
  paw_an(my_natom) <type(paw_an_type)>=paw arrays given on angular mesh
  pawang <type(pawang_type)>=paw angular mesh and related data
  pawrad(ntypat) <type(pawrad_type)>=paw radial mesh and related data
  pawrhoij(my_natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data
  pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data
  prtefg=1 to print summary output, 2 for detailed output
  quadmom(ntypat)=quadrupole moments in barns of different atomic nuclei
  rhor(nfft,nspden)=electron density on grid (strictly $\tilde{n}+\hat{n}$)
  rprimd(3,3)=matrix relating cartesian coordinates to crystal coordinates
  symrel(3,3,nsym)=symmetry operators in terms of action on primitive translations
  tnons(3,nsym)=nonsymmorphic translations
  typat(natom)=type (integer) for each atom
  ucvol=unit cell volume in Bohr^3
  usepaw=1 if we are using PAW formalism, 0 else
  xred(3,natom)=vectors locating each atom in the unit cell, in crystal coords
  zion(ntypat)=net core charge on each type of atom

OUTPUT

  (only writing, printing)

SIDE EFFECTS

NOTES

PARENTS

      outscfcv

CHILDREN

      dsyev,free_my_atmtab,get_my_atmtab,make_efg_el,make_efg_ion
      make_efg_onsite,wrtout

SOURCE

124   subroutine calc_efg(mpi_enreg,my_natom,natom,nfft,ngfft,nspden,nsym,ntypat,paral_kgb,&
125        &                    paw_an,pawang,pawrad,pawrhoij,pawtab,&
126        &                    ptcharge,prtefg,quadmom,rhor,rprimd,symrel,tnons,typat,ucvol,usepaw,xred,zion,&
127        &                    mpi_atmtab,comm_atom) ! optional arguments (parallelism)
128 
129 
130 !This section has been created automatically by the script Abilint (TD).
131 !Do not modify the following lines by hand.
132 #undef ABI_FUNC
133 #define ABI_FUNC 'calc_efg'
134 !End of the abilint section
135 
136     implicit none
137 
138     !Arguments ------------------------------------
139     !scalars
140     integer,intent(in) :: my_natom,natom,nfft,nspden,nsym,ntypat,paral_kgb,prtefg,usepaw
141     integer,optional,intent(in) :: comm_atom
142     real(dp),intent(in) :: ucvol
143     type(MPI_type),intent(in) :: mpi_enreg
144     type(pawang_type),intent(in) :: pawang
145     !arrays
146     integer,intent(in) :: ngfft(18),symrel(3,3,nsym),typat(natom)
147     integer,optional,target,intent(in) :: mpi_atmtab(:)
148     real(dp),intent(in) :: ptcharge(ntypat)
149     real(dp),intent(in) :: quadmom(ntypat),rhor(nfft,nspden),rprimd(3,3)
150     real(dp),intent(in) :: tnons(3,nsym),zion(ntypat)
151     real(dp),intent(inout) :: xred(3,natom)
152     type(paw_an_type),intent(in) :: paw_an(my_natom)
153     type(pawrad_type),intent(in) :: pawrad(ntypat)
154     type(pawrhoij_type),intent(in) :: pawrhoij(my_natom)
155     type(pawtab_type),intent(in) :: pawtab(ntypat)
156 
157     !Local variables-------------------------------
158     !scalars
159     integer :: INFO,LDA,LWORK,N,iatom,my_comm_atom
160     logical :: my_atmtab_allocated,paral_atom
161     real(dp) :: cq,eta,vxx,vyy,vzz
162     character(len=500) :: message
163     !arrays
164     integer,pointer :: my_atmtab(:)
165     real(dp) :: eigval(3),matr(3,3),work(8)
166     real(dp),allocatable :: efg(:,:,:),efg_el(:,:,:),efg_ion(:,:,:),efg_paw(:,:,:)
167     real(dp),allocatable :: efg_point_charge(:,:,:)
168 
169     ! ************************************************************************
170 
171     !Compatibility tests
172     if (usepaw /= 1) then
173        message = ' usepaw /= 1 but EFG calculation requires PAW '
174        MSG_ERROR(message)
175     end if
176 
177     !Set up parallelism over atoms
178     paral_atom=(present(comm_atom).and.(my_natom/=natom))
179     nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
180     my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
181     call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,my_natom_ref=my_natom)
182 
183     ABI_ALLOCATE(efg,(3,3,natom))
184     ABI_ALLOCATE(efg_el,(3,3,natom))
185     ABI_ALLOCATE(efg_ion,(3,3,natom))
186     ABI_ALLOCATE(efg_paw,(3,3,natom))
187     ABI_ALLOCATE(efg_point_charge,(3,3,natom))
188     efg_el(:,:,:) = zero
189     efg_ion(:,:,:) = zero
190     efg_paw(:,:,:) = zero
191     efg_point_charge(:,:,:) = zero
192 
193     call make_efg_el(efg_el,mpi_enreg,natom,nfft,ngfft,nspden,nsym,paral_kgb,rhor,rprimd,symrel,tnons,xred)
194 
195     call make_efg_ion(efg_ion,natom,nsym,ntypat,rprimd,symrel,tnons,typat,ucvol,xred,zion)
196 
197     if (paral_atom) then
198        call make_efg_onsite(efg_paw,my_natom,natom,nsym,ntypat,paw_an,pawang,pawrhoij,pawrad,pawtab,&
199             &   rprimd,symrel,tnons,xred,comm_atom=my_comm_atom,mpi_atmtab=my_atmtab)
200     else
201        call make_efg_onsite(efg_paw,my_natom,natom,nsym,ntypat,paw_an,pawang,pawrhoij,pawrad,pawtab,&
202             &   rprimd,symrel,tnons,xred)
203     end if
204 
205     !calculate efg due to pure point charges, as input in variable ptcharge(ntypat)
206     !note here all atoms of the same type will have the same valence; in the future this
207     !could be made more flexible by having ptcharge(natom) but that will require a slightly
208     !different version than the existing make_efg_ion routine
209     if(prtefg > 2) then
210        call make_efg_ion(efg_point_charge,natom,nsym,ntypat,rprimd,symrel,tnons,typat,ucvol,xred,ptcharge)
211     end if
212 
213     efg(:,:,:) = efg_el(:,:,:) + efg_ion(:,:,:) + efg_paw(:,:,:)
214 
215     write(message,'(a,a,a)' ) ch10,' Electric Field Gradient Calculation ',ch10
216     call wrtout(ab_out,message,'COLL')
217 
218     LDA=3; LWORK=8;N=3 ! these parameters are needed for the LAPACK dsyev routine
219     do iatom = 1, natom
220        matr(:,:) = efg(:,:,iatom)
221        call dsyev('V','U',N,matr,LDA,eigval,work,LWORK,INFO) ! get eigenvalues and eigenvectors
222        if (eigval(3) > abs(eigval(1)) ) then ! In NMR, the convention is that whatever component is
223           !    largest in magnitude is called Vzz, next comes Vxx, then Vyy
224           vzz = eigval(3)
225           vxx = eigval(1)
226           vyy = eigval(2)
227        else
228           vzz = eigval(1)
229           vxx = eigval(3)
230           vyy = eigval(2)
231        end if
232        if (abs(quadmom(typat(iatom))) > tol8 ) then ! only relevant when quadmom > 0 for a given atom
233           !    cq = (eQ)*Vzz/h, where Q is the electric quadrupole moment and Vzz is the largest in magnitude
234           !    principal component of the EFG tensor. Q is input in quadmom in barns, and Vzz is computed in atomic
235           !    units. The factor 2349647.81 Ha^{-1}Bohr^2 fm^{-2} sec^-1 converts from atomic units to frequency (see
236           !    http://www.ismar.org/ISMARpedia/index.php/Nuclear_Quadrupole_Resonance for discussion); we divide by
237           !    10^6 to convert to MHz from Hz and multiply by 100 to convert from fm^2 to Barns.
238           cq = vzz*quadmom(typat(iatom))*2349647.81/1.0E4
239           if(abs(cq) > tol6 )then ! if Cq is non-zero, eta is meaningful, otherwise it s numerical noise
240              eta = abs(vxx - vyy)/abs(vzz)
241           else
242              eta=zero
243           end if
244        else
245           cq =zero
246           eta =zero
247        end if
248        !  we always write Cq and eta, these are the NMR observables
249        write(message,'(a,i3,a,i3,a,f13.6,a,f13.6)') ' Atom ',iatom,', typat ',typat(iatom),': Cq = ',cq,' MHz     eta = ',eta
250        call wrtout(ab_out,message,'COLL')
251        if (prtefg > 1) then ! print detailed results on component EFG's
252           write(message,'(a,a,f13.6,a,a,3f13.6)')ch10,'      efg eigval : ',eigval(1),ch10,&
253                &     '-         eigvec : ',matr(1,1),matr(2,1),matr(3,1)
254           call wrtout(ab_out,message,'COLL')
255           write(message,'(a,f13.6,a,a,3f13.6)')'      efg eigval : ',eigval(2),ch10,&
256                &     '-         eigvec : ',matr(1,2),matr(2,2),matr(3,2)
257           call wrtout(ab_out,message,'COLL')
258           write(message,'(a,f13.6,a,a,3f13.6)')'      efg eigval : ',eigval(3),ch10,&
259                &     '-         eigvec : ',matr(1,3),matr(2,3),matr(3,3)
260           call wrtout(ab_out,message,'COLL')
261           write(message,'(a,a,3f13.6)')ch10,'      total efg : ',efg(1,1,iatom),efg(1,2,iatom),efg(1,3,iatom)
262           call wrtout(ab_out,message,'COLL')
263           write(message,'(a,3f13.6)')'      total efg : ',efg(2,1,iatom),efg(2,2,iatom),efg(2,3,iatom)
264           call wrtout(ab_out,message,'COLL')
265           write(message,'(a,3f13.6,a)')'      total efg : ',efg(3,1,iatom),efg(3,2,iatom),efg(3,3,iatom),ch10
266           call wrtout(ab_out,message,'COLL')
267           write(message,'(a,a,3f13.6)')ch10,'      efg_el : ',efg_el(1,1,iatom),efg_el(1,2,iatom),efg_el(1,3,iatom)
268           call wrtout(ab_out,message,'COLL')
269           write(message,'(a,3f13.6)')'      efg_el : ',efg_el(2,1,iatom),efg_el(2,2,iatom),efg_el(2,3,iatom)
270           call wrtout(ab_out,message,'COLL')
271           write(message,'(a,3f13.6,a)')'      efg_el : ',efg_el(3,1,iatom),efg_el(3,2,iatom),efg_el(3,3,iatom),ch10
272           call wrtout(ab_out,message,'COLL')
273           write(message,'(a,3f13.6)')'      efg_ion : ',efg_ion(1,1,iatom),efg_ion(1,2,iatom),efg_ion(1,3,iatom)
274           call wrtout(ab_out,message,'COLL')
275           write(message,'(a,3f13.6)')'      efg_ion : ',efg_ion(2,1,iatom),efg_ion(2,2,iatom),efg_ion(2,3,iatom)
276           call wrtout(ab_out,message,'COLL')
277           write(message,'(a,3f13.6,a)')'      efg_ion : ',efg_ion(3,1,iatom),efg_ion(3,2,iatom),efg_ion(3,3,iatom),ch10
278           call wrtout(ab_out,message,'COLL')
279           write(message,'(a,3f13.6)')'      efg_paw : ',efg_paw(1,1,iatom),efg_paw(1,2,iatom),efg_paw(1,3,iatom)
280           call wrtout(ab_out,message,'COLL')
281           write(message,'(a,3f13.6)')'      efg_paw : ',efg_paw(2,1,iatom),efg_paw(2,2,iatom),efg_paw(2,3,iatom)
282           call wrtout(ab_out,message,'COLL')
283           write(message,'(a,3f13.6,a)')'      efg_paw : ',efg_paw(3,1,iatom),efg_paw(3,2,iatom),efg_paw(3,3,iatom),ch10
284           call wrtout(ab_out,message,'COLL')
285        end if
286        if (prtefg > 2) then ! write output of pure pointcharge calculation
287           matr(:,:) = efg_point_charge(:,:,iatom)
288           call dsyev('V','U',N,matr,LDA,eigval,work,LWORK,INFO) ! get eigenvalues and eigenvectors
289           if (eigval(3) > abs(eigval(1)) ) then ! In NMR, the convention is that whatever component is
290              !      largest in magnitude is called Vzz, next comes Vxx, then Vyy
291              vzz = eigval(3)
292              vxx = eigval(1)
293              vyy = eigval(2)
294           else
295              vzz = eigval(1)
296              vxx = eigval(3)
297              vyy = eigval(2)
298           end if
299           if (abs(quadmom(typat(iatom))) > tol8 ) then ! only relevant when quadmom > 0 for a given atom
300              !      cq = e2Qq/h, where Vzz = eq and quadmom = Q; the other factors convert from atomic units to MHz
301              cq = vzz*quadmom(typat(iatom))*2349647.81/1.0E4
302              if(abs(cq) > tol6 )then ! if Cq is non-zero, eta is meaningful, otherwise it s numerical noise
303                 eta = abs(vxx - vyy)/abs(vzz)
304              else
305                 eta=zero
306              end if
307           else
308              cq =zero
309              eta =zero
310           end if
311           !    we always write Cq and eta, these are the NMR observables
312           write(message,'(a,i3,a,i3,a,f13.6,a,f13.6)') ' Atom ',iatom,', typat ',typat(iatom),&
313                &     ': Point charge Cq = ',cq,' MHz     eta = ',eta
314           call wrtout(ab_out,message,'COLL')
315           write(message,'(a,a,f13.6,a,a,3f13.6)')ch10,'      point charge efg eigval : ',eigval(1),ch10,&
316                &     '-         eigvec : ',matr(1,1),matr(2,1),matr(3,1)
317           call wrtout(ab_out,message,'COLL')
318           write(message,'(a,f13.6,a,a,3f13.6)')'      point charge efg eigval : ',eigval(2),ch10,&
319                &     '-         eigvec : ',matr(1,2),matr(2,2),matr(3,2)
320           call wrtout(ab_out,message,'COLL')
321           write(message,'(a,f13.6,a,a,3f13.6)')'      point charge efg eigval : ',eigval(3),ch10,&
322                &     '-         eigvec : ',matr(1,3),matr(2,3),matr(3,3)
323           call wrtout(ab_out,message,'COLL')
324           write(message,'(a,a,3f13.6)')ch10,'      point charge efg : ',efg_point_charge(1,1,iatom),&
325                &     efg_point_charge(1,2,iatom),efg_point_charge(1,3,iatom)
326           call wrtout(ab_out,message,'COLL')
327           write(message,'(a,3f13.6)')'      point charge efg : ',efg_point_charge(2,1,iatom),&
328                &     efg_point_charge(2,2,iatom),efg_point_charge(2,3,iatom)
329           call wrtout(ab_out,message,'COLL')
330           write(message,'(a,3f13.6,a)')'      point charge efg : ',efg_point_charge(3,1,iatom),&
331                &     efg_point_charge(3,2,iatom),efg_point_charge(3,3,iatom),ch10
332           call wrtout(ab_out,message,'COLL')
333        end if
334     end do
335     write(message,'(3a)')ch10,ch10,ch10
336     call wrtout(ab_out,message,'COLL')
337 
338     ABI_DEALLOCATE(efg)
339     ABI_DEALLOCATE(efg_el)
340     ABI_DEALLOCATE(efg_ion)
341     ABI_DEALLOCATE(efg_paw)
342     ABI_DEALLOCATE(efg_point_charge)
343 
344     !Destroy atom table used for parallelism
345     call free_my_atmtab(my_atmtab,my_atmtab_allocated)
346 
347     !DEBUG
348     !write(std_out,*)' calc_efg : exit '
349     !stop
350     !ENDDEBUG
351 
352   end subroutine calc_efg

ABINIT/calc_fc [ Functions ]

[ Top ] [ Functions ]

NAME

 calc_fc

FUNCTION

 calculation and output of Fermi-contact term at each atomic site

COPYRIGHT

 Copyright (C) 2009-2018 ABINIT group (JWZ,MT)
 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

  mpi_atmtab(:)=--optional-- indexes of the atoms treated by current proc
  comm_atom=--optional-- MPI communicator over atoms
  my_natom=number of atoms treated by current processor
  natom=number of atoms in cell.
  nspden=number of spin density components
  ntypat=number of atom types
  pawrad(ntypat) <type(pawrad_type)>=paw radial mesh and related data
  pawrhoij(natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data
  pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data
  typat(natom)=type (integer) for each atom
  usepaw=1 if PAW is activated

OUTPUT

  (only writing, printing)

SIDE EFFECTS

NOTES

PARENTS

      outscfcv

CHILDREN

      free_my_atmtab,get_my_atmtab,make_fc_paw,wrtout

SOURCE

399   subroutine calc_fc(my_natom,natom,nspden,ntypat,pawrad,pawrhoij,pawtab,typat,usepaw,&
400        &                  mpi_atmtab,comm_atom) ! optional arguments (parallelism)
401 
402 
403 !This section has been created automatically by the script Abilint (TD).
404 !Do not modify the following lines by hand.
405 #undef ABI_FUNC
406 #define ABI_FUNC 'calc_fc'
407 !End of the abilint section
408 
409     implicit none
410 
411     !Arguments ------------------------------------
412     !scalars
413     integer,intent(in) :: my_natom,natom,nspden,ntypat,usepaw
414     integer,optional,intent(in) :: comm_atom
415     !arrays
416     integer,intent(in) :: typat(natom)
417     integer,optional,target,intent(in) :: mpi_atmtab(:)
418     type(pawrad_type),intent(in) :: pawrad(ntypat)
419     type(pawrhoij_type),intent(in) :: pawrhoij(my_natom)
420     type(pawtab_type),intent(in) :: pawtab(ntypat)
421 
422     !Local variables-------------------------------
423     !scalars
424     integer :: iatom,my_comm_atom
425     logical :: my_atmtab_allocated,paral_atom
426     character(len=500) :: message
427     !arrays
428     integer,pointer :: my_atmtab(:)
429     real(dp),allocatable :: fc(:,:)
430 
431 !***********************************************************************
432 
433     !Compatibility tests
434     if (usepaw /= 1) then
435        message = ' usepaw /= 1 but Fermi-contact calculation requires PAW '
436        MSG_ERROR(message)
437     end if
438 
439     !Set up parallelism over atoms
440     paral_atom=(present(comm_atom).and.(my_natom/=natom))
441     nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
442     my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
443     call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,my_natom_ref=my_natom)
444 
445     !Initialization
446     ABI_ALLOCATE(fc,(nspden,natom))
447 
448     !Computation
449     if (paral_atom) then
450        call make_fc_paw(fc,my_natom,natom,nspden,ntypat,pawrhoij,pawrad,pawtab,&
451             &   comm_atom=my_comm_atom,mpi_atmtab=my_atmtab)
452     else
453        call make_fc_paw(fc,my_natom,natom,nspden,ntypat,pawrhoij,pawrad,pawtab)
454     end if
455 
456     !Printing
457     write(message,'(a,a,a)' ) ch10,' Fermi-contact Term Calculation ',ch10
458     call wrtout(ab_out,message,'COLL')
459 
460     do iatom = 1, natom
461        if (nspden == 2) then
462           write(message,'(a,i3,a,i3,a,f12.4)') ' Atom ',iatom,', typat ',typat(iatom),': FC total = ',&
463                &     fc(1,iatom)+fc(2,iatom)
464           call wrtout(ab_out,message,'COLL')
465           write(message,'(a,i3,a,i3,a,f12.4)') ' Atom ',iatom,', typat ',typat(iatom),': FC up - down = ',&
466                &     fc(1,iatom)-fc(2,iatom)
467           call wrtout(ab_out,message,'COLL')
468        else
469           write(message,'(a,i3,a,i3,a,f12.4)') ' Atom ',iatom,', typat ',typat(iatom),': FC = ',&
470                &     fc(1,iatom)
471           call wrtout(ab_out,message,'COLL')
472        end if
473     end do
474 
475     write(message,'(3a)')ch10,ch10,ch10
476     call wrtout(ab_out,message,'COLL')
477 
478     !Memory deallocation
479     ABI_DEALLOCATE(fc)
480 
481     !Destroy atom table used for parallelism
482     call free_my_atmtab(my_atmtab,my_atmtab_allocated)
483 
484   end subroutine calc_fc

ABINIT/m_nucprop [ Modules ]

[ Top ] [ Modules ]

NAME

  m_nucprop

FUNCTION

  routines used to compute properties at the nuclear sites, including
  electric field gradient and Fermi contact

COPYRIGHT

 Copyright (C) 1998-2018 ABINIT group (MT, JWZ)
  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

22 #if defined HAVE_CONFIG_H
23 #include "config.h"
24 #endif
25 
26 #include "abi_common.h"
27 
28 module m_nucprop
29 
30   use defs_basis
31   use m_abicore
32   use m_errors
33 
34   use defs_abitypes, only : MPI_type
35   use m_mpinfo,   only : ptabs_fourdp
36   use m_xmpi, only : xmpi_comm_self, xmpi_sum
37   use m_geometry,       only : xred2xcart
38   use m_linalg_interfaces, only: dsyev
39   use m_paw_an,      only : paw_an_type
40   use m_pawang,      only : pawang_type
41   use m_pawrad,     only : pawrad_type
42   use m_pawtab,     only : pawtab_type
43   use m_pawrhoij,   only : pawrhoij_type
44   use m_paw_nmr,    only : make_efg_onsite,make_fc_paw
45   use m_paral_atom, only : get_my_atmtab,free_my_atmtab
46   use m_special_funcs,  only : abi_derfc
47   use m_symtk,          only : matr3inv, matpointsym
48   use m_fft,           only : fourdp
49 
50   implicit none
51 
52   private

ABINIT/make_efg_el [ Functions ]

[ Top ] [ Functions ]

NAME

 make_efg_el

FUNCTION

 compute the electric field gradient due to electron density

COPYRIGHT

 Copyright (C) 2005-2018 ABINIT group (JWZ)
 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

 mpi_enreg=information about MPI parallelization
 natom, number of atoms in unit cell
 nfft,ngfft(18), number of FFT points and details of FFT
 nspden, number of spin components
 nsym=number of symmetries in space group
 rhor(nfft,nspden), valence electron density, here $\tilde{n} + \hat{n}$
 rprimd(3,3), conversion from crystal coordinates to cartesian coordinates
 symrel(3,3,nsym)=symmetry operators in terms of action on primitive translations
 tnons(3,nsym) = nonsymmorphic translations
 xred(3,natom), location of atoms in crystal coordinates.

OUTPUT

 efg(3,3,natom), the 3x3 efg tensor at each atomic site due to rhor

NOTES

 This routine computes the electric field gradient, specifically the components
 $\partial^2 V/\partial x_\alpha \partial x_\beta$ of the potential generated by the valence
 electrons, at each atomic site in the unit cell. Key references: Kresse and Joubert, ``From
 ultrasoft pseudopotentials to the projector augmented wave method'', Phys. Rev. B. 59, 1758--1775 (1999) [[cite:Kresse1999]],
 and Profeta, Mauri, and Pickard, ``Accurate first principles prediction of $^{17}$O NMR parameters in
 SiO$_2$: Assignment of the zeolite ferrierite spectrum'', J. Am. Chem. Soc. 125, 541--548 (2003) [[cite:Profeta2003]]. This
 routine computes the second derivatives of the potential generated by $\tilde{n}$ (see Kresse and Joubert
 for notation, Fourier-transforming the density, doing the sum in G space, and then transforming back at
 each atomic site. The final formula is
 \begin{displaymath}
 \frac{\partial^2 V}{\partial x_\alpha\partial x_\beta} = -4\pi^2\sum_G (G_\alpha G_\beta - \delta_{\alpha,\beta}G^2/3)
 \left(\frac{\tilde{n}(G)}{\pi G^2}\right)e^{2\pi i G\cdot R}
 \end{displaymath}

PARENTS

      calc_efg

CHILDREN

      fourdp,matpointsym,matr3inv,ptabs_fourdp,xmpi_sum,xred2xcart

SOURCE

782 subroutine make_efg_el(efg,mpi_enreg,natom,nfft,ngfft,nspden,nsym,paral_kgb,rhor,rprimd,symrel,tnons,xred)
783 
784 
785 !This section has been created automatically by the script Abilint (TD).
786 !Do not modify the following lines by hand.
787 #undef ABI_FUNC
788 #define ABI_FUNC 'make_efg_el'
789 !End of the abilint section
790 
791   implicit none
792 
793   !Arguments ------------------------------------
794   !scalars
795   integer,intent(in) :: natom,nfft,nspden,nsym,paral_kgb
796   type(MPI_type),intent(in) :: mpi_enreg
797   !arrays
798   integer,intent(in) :: ngfft(18),symrel(3,3,nsym)
799   real(dp),intent(in) :: rhor(nfft,nspden),rprimd(3,3),tnons(3,nsym),xred(3,natom)
800   real(dp),intent(out) :: efg(3,3,natom)
801 
802   !Local variables-------------------------------
803   !scalars
804   integer :: cplex,fftdir,fofg_index,iatom,i1,i2,i2_local,i23,i3,id1,id2,id3
805   integer :: ierr,ig,ig2,ig3,ii,ii1,ing,jj
806   integer :: me_fft,n1,n2,n3,nproc_fft,tim_fourdp
807   real(dp) :: cph,derivs,phase,sph,trace
808   ! type(MPI_type) :: mpi_enreg_seq
809   !arrays
810   integer :: id(3)
811   integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:)
812   integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:)
813   real(dp) :: gprimd(3,3),gred(3),gvec(3),ratom(3)
814   real(dp),allocatable :: fofg(:,:),fofr(:),gq(:,:),xcart(:,:)
815 
816   ! ************************************************************************
817 
818   !DEBUG
819   !write(std_out,*)' make_efg_el : enter'
820   !ENDDEBUG
821 
822   ABI_ALLOCATE(fofg,(2,nfft))
823   ABI_ALLOCATE(fofr,(nfft))
824   ABI_ALLOCATE(xcart,(3,natom))
825 
826   efg(:,:,:) = zero
827   call xred2xcart(natom,rprimd,xcart,xred) ! get atomic locations in cartesian coords
828   call matr3inv(rprimd,gprimd)
829 
830   tim_fourdp = 0 ! timing code, not using
831   fftdir = -1 ! FT from R to G
832   cplex = 1 ! fofr is real
833   !here we are only interested in the total charge density including nhat, which is rhor(:,1)
834   !regardless of the value of nspden. This may change in the future depending on
835   !developments with noncollinear magnetization and so forth. Such a change will
836   !require an additional loop over nspden.
837   !Multiply by -1 to convert the electron particle density to the charge density
838   fofr(:) = -rhor(:,1)
839 
840   ! Get the distrib associated with this fft_grid  See hartre.F90 for another example where
841   ! this is done
842   n1=ngfft(1); n2=ngfft(2); n3=ngfft(3)
843   nproc_fft = mpi_enreg%nproc_fft; me_fft = mpi_enreg%me_fft
844   call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
845 
846   call fourdp(cplex,fofg,fofr,fftdir,mpi_enreg,nfft,ngfft,paral_kgb,tim_fourdp) ! construct charge density in G space
847 
848   ! the following loops over G vectors has been copied from hartre.F90 in order to be compatible with
849   ! possible FFT parallelism
850 
851   ! In order to speed the routine, precompute the components of g
852   ! Also check if the booked space was large enough...
853   ABI_ALLOCATE(gq,(3,max(n1,n2,n3)))
854   do ii=1,3
855      id(ii)=ngfft(ii)/2+2
856      do ing=1,ngfft(ii)
857         ig=ing-(ing/id(ii))*ngfft(ii)-1
858         gq(ii,ing)=ig
859      end do
860   end do
861   id1=n1/2+2;id2=n2/2+2;id3=n3/2+2
862 
863   ! Triple loop on each dimension
864   do i3=1,n3
865      ig3=i3-(i3/id3)*n3-1
866      gred(3) = gq(3,i3)
867 
868      do i2=1,n2
869         ig2=i2-(i2/id2)*n2-1
870         if (fftn2_distrib(i2) == me_fft) then
871 
872            gred(2) = gq(2,i2)
873            i2_local = ffti2_local(i2)
874            i23=n1*(i2_local-1 +(n2/nproc_fft)*(i3-1))
875            ! Do the test that eliminates the Gamma point outside of the inner loop
876            ii1=1
877            if(i23==0 .and. ig2==0 .and. ig3==0) ii1=2
878 
879            ! Final inner loop on the first dimension (note the lower limit)
880            do i1=ii1,n1
881               !         gs=gs2+ gq(1,i1)*(gq(1,i1)*gmet(1,1)+gqg2p3)
882               gred(1) = gq(1,i1)
883               gvec(1:3) = MATMUL(gprimd,gred)
884               fofg_index=i1+i23
885               trace = dot_product(gvec,gvec)
886               do ii = 1, 3 ! sum over components of efg tensor
887                  do jj = 1, 3 ! sum over components of efg tensor
888                     derivs = gvec(ii)*gvec(jj) ! This term is $G_\alpha G_\beta$
889                     if (ii == jj) derivs = derivs - trace/three
890                     do iatom = 1, natom ! sum over atoms in unit cell
891                        ratom(:) = xcart(:,iatom) ! extract location of atom iatom
892                        phase = two_pi*dot_product(gvec,ratom) ! argument of $e^{2\pi i G\cdot R}$
893                        cph = cos(phase)
894                        sph = sin(phase)
895                        efg(ii,jj,iatom) = efg(ii,jj,iatom) - &
896                             &               four_pi*derivs*(fofg(1,fofg_index)*cph-fofg(2,fofg_index)*sph)/trace ! real part of efg tensor
897                     end do ! end loop over atoms in cell
898                  end do ! end loop over jj in V_ij
899               end do ! end loop over ii in V_ij
900            end do ! End loop on i1
901         end if
902      end do ! End loop on i2
903   end do ! End loop on i3
904 
905   call xmpi_sum(efg,mpi_enreg%comm_fft,ierr)
906 
907   ! symmetrize tensor at each atomic site using point symmetry operations
908   do iatom = 1, natom
909      call matpointsym(iatom,efg(:,:,iatom),natom,nsym,rprimd,symrel,tnons,xred)
910   end do
911 
912   ABI_DEALLOCATE(fofg)
913   ABI_DEALLOCATE(fofr)
914   ABI_DEALLOCATE(xcart)
915   ABI_DEALLOCATE(gq)
916 
917   !DEBUG
918   !write(std_out,*)' make_efg_el : exit '
919   !stop
920   !ENDDEBUG
921 
922 end subroutine make_efg_el

ABINIT/make_efg_ion [ Functions ]

[ Top ] [ Functions ]

NAME

 make_efg_ion

FUNCTION

 compute the electric field gradient due to ionic cores

COPYRIGHT

 Copyright (C) 2005-2018 ABINIT group (JWZ)
 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

 natom, number of atoms in the unit cell
 nsym=number of symmetries in space group
 ntypat, the number of types of atoms in the unit cell
 rprimd(3,3), the matrix giving the transformation from crystal to cartesian coordinates
 symrel(3,3,nsym)=symmetry operators in terms of action on primitive translations
 tnons(3,nsym) = nonsymmorphic translations
 typat(natom), the type of each atom in the unit cell
 ucvol, the volume of the unit cell in atomic units
 xred(3,natom) the location of each atom in the cell in crystallographic coordinates
 zion(ntypat) the net charge on each type of atom

OUTPUT

 efg(3,3,natom), the 3x3 efg tensors at each atomic site

SIDE EFFECTS

NOTES

 This routine computes the electric field gradient, specifically the components
 $\partial^2 V/\partial x_\alpha \partial x_\beta$ of the potential generated by the ionic cores,
 at each atomic site in the unit cell.
 Key references:
 Profeta, Mauri, and Pickard, ``Accurate first principles prediction of $^{17}$O NMR parameters in
 SiO$_2$: Assignment of the zeolite ferrierite spectrum'', J. Am. Chem. Soc. 125, 541--548 (2003) [[cite:Profeta2003]];
 A. Honma, ``Dipolar lattice-sums with applications to the exciton bands of anthracene crystal and
 the crystal field due to point charges'', J. Phys. Soc. Jpn. 42, 1129--1135 (1977) [[cite:Honma1977]];
 and Kresse and Joubert, ``From ultrasoft pseudopotentials to the projector augmented wave method'',
 Phys. Rev. B. 59, 1758--1775 (1999) [[cite:Kresse1999]]. In Kresse and Joubert's notation, the ionic cores are $n_{Zc}$;
 these charges are given by the net core charges on the pseudoatoms. Due to otherwise slow convergence,
 the sum over atoms is carried out by an Ewald method as detailed in the Honma reference, specifically
 his Eq. 4.8.

PARENTS

CHILDREN

SOURCE

539 subroutine make_efg_ion(efg,natom,nsym,ntypat,rprimd,symrel,tnons,typat,ucvol,xred,zion)
540 
541 
542 !This section has been created automatically by the script Abilint (TD).
543 !Do not modify the following lines by hand.
544 #undef ABI_FUNC
545 #define ABI_FUNC 'make_efg_ion'
546 !End of the abilint section
547 
548   implicit none
549 
550   !Arguments ------------------------------------
551   !scalars
552   integer,intent(in) :: natom,nsym,ntypat
553   real(dp) :: ucvol
554   !arrays
555   integer,intent(in) :: symrel(3,3,nsym),typat(natom)
556   real(dp),intent(in) :: rprimd(3,3),tnons(3,nsym)
557   real(dp),intent(in) :: zion(ntypat)
558   real(dp),intent(inout) :: xred(3,natom)
559   real(dp),intent(out) :: efg(3,3,natom)
560   !Local variables-------------------------------
561   !scalars
562   integer :: iatom,ishell,ii,jatom,jj,nshell,sx,sy,sz
563   real(dp) :: cph,dampfac,derfc_karg,derivs,gsq,karg
564   real(dp) :: lenrho,phase,qk,rlkcut,trace,xi0
565   real(dp) :: glkcut
566   !arrays
567   real(dp) :: cvec(3),gvec(3),gpl(3),gprimd(3,3)
568   real(dp) :: rhok(3),rhored(3),rpl(3)
569   real(dp),allocatable :: efg_g(:,:,:),efg_r(:,:,:)
570   real(dp),allocatable :: xcart(:,:)
571 
572   ! ************************************************************************
573 
574   !DEBUG
575   !write(std_out,*)' make_efg_ion : enter'
576   !ENDDEBUG
577 
578   ABI_ALLOCATE(efg_g,(3,3,natom))
579   ABI_ALLOCATE(efg_r,(3,3,natom))
580   ABI_ALLOCATE(xcart,(3,natom))
581   efg(:,:,:) = zero ! final efg tensor
582   efg_g(:,:,:) = zero ! part of tensor accumulated in G space
583   efg_r(:,:,:) = zero ! part of tensor accumulated in R space
584 
585   call xred2xcart(natom,rprimd,xcart,xred) ! get atomic locations in cartesian coords
586 
587   do ii = 1, 3 ! generate the lengths of the unit cell edges in atomic units
588      rpl(ii) = sqrt(rprimd(1,ii)**2+rprimd(2,ii)**2+rprimd(3,ii)**2)
589   end do
590   xi0 = sqrt(pi/(maxval(rpl)*minval(rpl))) ! this estimate for xi0 is from Honma's paper
591 
592   call matr3inv(rprimd,gprimd) ! gprimd holds the inverse transpose of rprimd
593   !remember ordering: rprimd( (x_comp,y_comp,z_comp), (edge 1, edge 2, edge 3) )
594   !while gprimd( (edge 1, edge 2, edge 3),(x_comp, y_comp, z_comp) )
595   do ii = 1, 3 ! generate the lengths of the reciprocal cell edges
596      gpl(ii) = sqrt(gprimd(ii,1)**2+gprimd(ii,2)**2+gprimd(ii,3)**2)
597   end do
598 
599   !go out enough shells such that g**2/4*xi0**2 is of order 30
600   nshell = int(anint(sqrt(30.0)*xi0/(pi*minval(gpl))))
601   glkcut = (0.95*nshell*two*pi*minval(gpl))**2
602 
603   do ishell = 0, nshell ! loop over shells
604      do sx = -ishell, ishell
605         do sy = -ishell, ishell
606            do sz = -ishell, ishell
607               if ( .not. (sx==0 .and. sy==0 .and. sz==0) ) then ! avoid origin
608                  !          constrain to be on shell surface, not interior
609                  if ( abs(sx)==ishell .or. abs(sy)==ishell .or. abs(sz)==ishell ) then
610                     cvec(1)=sx;cvec(2)=sy;cvec(3)=sz
611                     !            make the g vector in cartesian coords
612                     gvec(:) = zero
613                     do ii = 1, 3
614                        do jj = 1, 3
615                           gvec(ii) = gvec(ii) + gprimd(ii,jj)*cvec(jj)*two*pi
616                        end do
617                     end do
618                     gsq = dot_product(gvec,gvec)
619                     if(gsq < glkcut) then
620                        dampfac = exp(-gsq/(4.0*xi0*xi0)) ! see Honma eq. 4.8
621                        do iatom = 1, natom
622                           do jatom = 1, natom
623                              qk = zion(typat(jatom)) ! charge on neighbor atom
624                              rhok = xcart(:,jatom)-xcart(:,iatom)
625                              phase = dot_product(gvec,rhok)
626                              cph = cos(phase)
627                              do ii = 1, 3
628                                 do jj = 1, 3
629                                    derivs = -3.0*gvec(ii)*gvec(jj)/gsq
630                                    if (ii == jj) derivs = 1.0 + derivs
631                                    efg_g(ii,jj,iatom) = efg_g(ii,jj,iatom) + &
632                                         &                       qk*cph*derivs*dampfac
633                                 end do ! end loop over jj
634                              end do ! end loop over ii
635                           end do ! end loop over jatom
636                        end do ! end loop over iatom
637                     end if ! constrain to gsq < glkcut
638                  end if ! end selection on shell edge
639               end if ! end avoidance of origin
640            end do ! end loop over sz
641         end do ! end loop over sy
642      end do ! end loop over sx
643   end do ! end loop over ishell
644 
645   !sum in real space begins here
646 
647   !go out enough shells such that (r*xi0)**2 is of order 30
648   nshell = int(anint(sqrt(30.)/(minval(rpl)*xi0)))
649   rlkcut = nshell*minval(rpl)*0.95
650 !
651   !go out enough shells so that rlkcut is of order 30 bohr
652   !nshell=int(anint(30.0/minval(rpl)))
653   !rlkcut = 0.95*nshell*minval(rpl)
654 
655   do ishell = 0, nshell ! total set of cells to loop over
656      do sx = -ishell, ishell ! loop over all cells in each dimension
657         do sy = -ishell, ishell
658            do sz = -ishell, ishell
659               !        constrain to shell surface, not interior
660               if ( abs(sx)==ishell .or. abs(sy)==ishell .or. abs(sz)==ishell ) then
661                  do jatom = 1, natom ! loop over atoms in shell cell
662                     do iatom = 1, natom ! loop over atoms in central unit cell
663                        if (.NOT. (jatom == iatom .AND. sx == 0 .AND. sy == 0 .AND. sz == 0)) then ! avoid self term
664                           qk = zion(typat(jatom)) ! charge on each neighbor atom
665                           !                ! rhored is the vector in crystal coords from neighbor to target
666                           rhored(1) = xred(1,jatom) + sx - xred(1,iatom)
667                           rhored(2) = xred(2,jatom) + sy - xred(2,iatom)
668                           rhored(3) = xred(3,jatom) + sz - xred(3,iatom)
669                           !                !  rhok is rhored in cartesian coords
670                           rhok(1) = rprimd(1,1)*rhored(1)+rprimd(1,2)*rhored(2)+rprimd(1,3)*rhored(3)
671                           rhok(2) = rprimd(2,1)*rhored(1)+rprimd(2,2)*rhored(2)+rprimd(2,3)*rhored(3)
672                           rhok(3) = rprimd(3,1)*rhored(1)+rprimd(3,2)*rhored(2)+rprimd(3,3)*rhored(3)
673                           trace = dot_product(rhok,rhok)
674                           lenrho = sqrt(trace)
675                           if (lenrho < rlkcut) then ! this restriction is critical as it ensures
676                              !                  ! that we sum over a sphere of atoms in real space
677                              !                  ! no matter what shape the unit cell has
678                              karg = xi0*lenrho
679                              derfc_karg = abi_derfc(karg)
680                              !                  see Honma eq. 2.10 for derivation of the following damping factor
681                              dampfac = (1.0+3.0/(2.0*karg*karg))*exp(-karg*karg)+3.0*sqrt(pi)*derfc_karg/(4.0*karg**3)
682                              do ii = 1, 3 ! loop over tensor elements
683                                 do jj = 1, 3 ! loop over tensor elements
684                                    derivs = -3.0*rhok(ii)*rhok(jj)/trace
685                                    if(ii == jj) derivs = derivs + 1.0 ! see Honma eq 4.8 re: sign
686                                    !                      accumulate real space tensor element,
687                                    !                      weighted by charge of neighbor and Ewald damping factor
688                                    efg_r(ii,jj,iatom) = efg_r(ii,jj,iatom) + qk*derivs*dampfac
689                                 end do ! end loop over jj in efg(ii,jj,iatom)
690                              end do ! end loop over ii in efg(ii,jj,iatom)
691                           end if ! end if statement restricting to a sphere of radius rlkcut
692                        end if ! end if statement avoiding the self atom term
693                     end do ! end loop over i atoms in cell
694                  end do ! end loop over j atoms in cell
695               end if ! end selection on outer shell of cells only
696            end do ! end loop over sz cells
697         end do ! end loop over sy cells
698      end do ! end loop over sx cells
699   end do ! end loop over shells
700 
701   !now combine the g-space and r-space parts, properly weighted (see Honma)
702   do iatom = 1, natom
703      do ii = 1, 3
704         do jj = 1, 3
705            efg(ii,jj,iatom) = four_pi*efg_g(ii,jj,iatom)/(three*ucvol)-&
706                 &       four*xi0**3*efg_r(ii,jj,iatom)/(three*sqrt(pi))
707            !      note extra factor of two: compare Honma eq. 4.6
708         end do
709      end do
710   end do
711 
712   ! symmetrize tensor at each atomic site using point symmetry operations
713   do iatom = 1, natom
714      call matpointsym(iatom,efg(:,:,iatom),natom,nsym,rprimd,symrel,tnons,xred)
715   end do
716 
717   ABI_DEALLOCATE(efg_g)
718   ABI_DEALLOCATE(efg_r)
719   ABI_DEALLOCATE(xcart)
720 
721   !DEBUG
722   !write(std_out,*)' make_efg_ion : exit '
723   !stop
724   !ENDDEBUG
725 
726 end subroutine make_efg_ion