TABLE OF CONTENTS
ABINIT/make_efg_el [ Functions ]
NAME
make_efg_el
FUNCTION
compute the electric field gradient due to electron density
COPYRIGHT
Copyright (C) 2005-2018 ABINIT group (JJ) 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),
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). 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
54 #if defined HAVE_CONFIG_H 55 #include "config.h" 56 #endif 57 #include "abi_common.h" 58 59 subroutine make_efg_el(efg,mpi_enreg,natom,nfft,ngfft,nspden,nsym,paral_kgb,rhor,rprimd,symrel,tnons,xred) 60 61 use defs_basis 62 use defs_abitypes 63 use m_profiling_abi 64 65 use m_geometry, only : xred2xcart 66 use m_mpinfo, only : ptabs_fourdp 67 use m_xmpi, only : xmpi_sum 68 69 !This section has been created automatically by the script Abilint (TD). 70 !Do not modify the following lines by hand. 71 #undef ABI_FUNC 72 #define ABI_FUNC 'make_efg_el' 73 use interfaces_32_util 74 use interfaces_45_geomoptim 75 use interfaces_53_ffts 76 !End of the abilint section 77 78 implicit none 79 80 !Arguments ------------------------------------ 81 !scalars 82 integer,intent(in) :: natom,nfft,nspden,nsym,paral_kgb 83 type(MPI_type),intent(in) :: mpi_enreg 84 !arrays 85 integer,intent(in) :: ngfft(18),symrel(3,3,nsym) 86 real(dp),intent(in) :: rhor(nfft,nspden),rprimd(3,3),tnons(3,nsym),xred(3,natom) 87 real(dp),intent(out) :: efg(3,3,natom) 88 89 !Local variables------------------------------- 90 !scalars 91 integer :: cplex,fftdir,fofg_index,iatom,i1,i2,i2_local,i23,i3,id1,id2,id3 92 integer :: ierr,ig,ig2,ig3,ii,ii1,ing,jj 93 integer :: me_fft,n1,n2,n3,nproc_fft,tim_fourdp 94 real(dp) :: cph,derivs,phase,sph,trace 95 ! type(MPI_type) :: mpi_enreg_seq 96 !arrays 97 integer :: id(3) 98 integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:) 99 integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:) 100 real(dp) :: gprimd(3,3),gred(3),gvec(3),ratom(3) 101 real(dp),allocatable :: fofg(:,:),fofr(:),gq(:,:),xcart(:,:) 102 103 ! ************************************************************************ 104 105 !DEBUG 106 !write(std_out,*)' make_efg_el : enter' 107 !ENDDEBUG 108 109 ABI_ALLOCATE(fofg,(2,nfft)) 110 ABI_ALLOCATE(fofr,(nfft)) 111 ABI_ALLOCATE(xcart,(3,natom)) 112 113 efg(:,:,:) = zero 114 call xred2xcart(natom,rprimd,xcart,xred) ! get atomic locations in cartesian coords 115 call matr3inv(rprimd,gprimd) 116 117 tim_fourdp = 0 ! timing code, not using 118 fftdir = -1 ! FT from R to G 119 cplex = 1 ! fofr is real 120 !here we are only interested in the total charge density including nhat, which is rhor(:,1) 121 !regardless of the value of nspden. This may change in the future depending on 122 !developments with noncollinear magnetization and so forth. Such a change will 123 !require an additional loop over nspden. 124 !Multiply by -1 to convert the electron particle density to the charge density 125 fofr(:) = -rhor(:,1) 126 127 ! Get the distrib associated with this fft_grid See hartre.F90 for another example where 128 ! this is done 129 n1=ngfft(1); n2=ngfft(2); n3=ngfft(3) 130 nproc_fft = mpi_enreg%nproc_fft; me_fft = mpi_enreg%me_fft 131 call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local) 132 133 call fourdp(cplex,fofg,fofr,fftdir,mpi_enreg,nfft,ngfft,paral_kgb,tim_fourdp) ! construct charge density in G space 134 135 ! the following loops over G vectors has been copied from hartre.F90 in order to be compatible with 136 ! possible FFT parallelism 137 138 ! In order to speed the routine, precompute the components of g 139 ! Also check if the booked space was large enough... 140 ABI_ALLOCATE(gq,(3,max(n1,n2,n3))) 141 do ii=1,3 142 id(ii)=ngfft(ii)/2+2 143 do ing=1,ngfft(ii) 144 ig=ing-(ing/id(ii))*ngfft(ii)-1 145 gq(ii,ing)=ig 146 end do 147 end do 148 id1=n1/2+2;id2=n2/2+2;id3=n3/2+2 149 150 ! Triple loop on each dimension 151 do i3=1,n3 152 ig3=i3-(i3/id3)*n3-1 153 gred(3) = gq(3,i3) 154 155 do i2=1,n2 156 ig2=i2-(i2/id2)*n2-1 157 if (fftn2_distrib(i2) == me_fft) then 158 159 gred(2) = gq(2,i2) 160 i2_local = ffti2_local(i2) 161 i23=n1*(i2_local-1 +(n2/nproc_fft)*(i3-1)) 162 ! Do the test that eliminates the Gamma point outside of the inner loop 163 ii1=1 164 if(i23==0 .and. ig2==0 .and. ig3==0) ii1=2 165 166 ! Final inner loop on the first dimension (note the lower limit) 167 do i1=ii1,n1 168 ! gs=gs2+ gq(1,i1)*(gq(1,i1)*gmet(1,1)+gqg2p3) 169 gred(1) = gq(1,i1) 170 gvec(1:3) = MATMUL(gprimd,gred) 171 fofg_index=i1+i23 172 trace = dot_product(gvec,gvec) 173 do ii = 1, 3 ! sum over components of efg tensor 174 do jj = 1, 3 ! sum over components of efg tensor 175 derivs = gvec(ii)*gvec(jj) ! This term is $G_\alpha G_\beta$ 176 if (ii == jj) derivs = derivs - trace/three 177 do iatom = 1, natom ! sum over atoms in unit cell 178 ratom(:) = xcart(:,iatom) ! extract location of atom iatom 179 phase = two_pi*dot_product(gvec,ratom) ! argument of $e^{2\pi i G\cdot R}$ 180 cph = cos(phase) 181 sph = sin(phase) 182 efg(ii,jj,iatom) = efg(ii,jj,iatom) - & 183 & four_pi*derivs*(fofg(1,fofg_index)*cph-fofg(2,fofg_index)*sph)/trace ! real part of efg tensor 184 end do ! end loop over atoms in cell 185 end do ! end loop over jj in V_ij 186 end do ! end loop over ii in V_ij 187 end do ! End loop on i1 188 end if 189 end do ! End loop on i2 190 end do ! End loop on i3 191 192 call xmpi_sum(efg,mpi_enreg%comm_fft,ierr) 193 194 ! symmetrize tensor at each atomic site using point symmetry operations 195 do iatom = 1, natom 196 call matpointsym(iatom,efg(:,:,iatom),natom,nsym,rprimd,symrel,tnons,xred) 197 end do 198 199 ABI_DEALLOCATE(fofg) 200 ABI_DEALLOCATE(fofr) 201 ABI_DEALLOCATE(xcart) 202 ABI_DEALLOCATE(gq) 203 204 !DEBUG 205 !write(std_out,*)' make_efg_el : exit ' 206 !stop 207 !ENDDEBUG 208 209 end subroutine make_efg_el