TABLE OF CONTENTS
ABINIT/make_efg_onsite [ Functions ]
NAME
make_efg_onsite
FUNCTION
Compute the electric field gradient due to onsite densities
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 my_natom=number of atoms treated by current processor natom=number of atoms in cell. nsym=number of symmetries in space group ntypat=number of atom types 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 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 site due to nhat
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). See in particular
Eq. 11 and 12 of Profeta et al., but note that their sum over occupied states times 2 for occupation number is
replaced in the Kresse and Joubert formulation by the sum over $\rho_{ij}$ occupations for each basis element pair.
PARENTS
calc_efg
CHILDREN
free_my_atmtab,get_my_atmtab,matpointsym,pawdensities,pawrad_deducer0 simp_gen,xmpi_sum
SOURCE
54 #if defined HAVE_CONFIG_H 55 #include "config.h" 56 #endif 57 58 #include "abi_common.h" 59 60 subroutine make_efg_onsite(efg,my_natom,natom,nsym,ntypat,paw_an,pawang,pawrhoij,pawrad,pawtab, & 61 & rprimd,symrel,tnons,xred,& 62 & mpi_atmtab,comm_atom) ! optional arguments (parallelism) 63 64 use defs_basis 65 use m_errors 66 use m_profiling_abi 67 use m_xmpi, only : xmpi_comm_self,xmpi_sum 68 69 use m_pawang, only : pawang_type 70 use m_pawtab, only : pawtab_type 71 use m_pawrad, only : pawrad_type, pawrad_deducer0, simp_gen 72 use m_pawtab, only : pawtab_type 73 use m_paw_an, only : paw_an_type 74 use m_pawrhoij, only : pawrhoij_type 75 use m_paral_atom, only : get_my_atmtab, free_my_atmtab 76 77 !This section has been created automatically by the script Abilint (TD). 78 !Do not modify the following lines by hand. 79 #undef ABI_FUNC 80 #define ABI_FUNC 'make_efg_onsite' 81 use interfaces_45_geomoptim 82 use interfaces_65_paw, except_this_one => make_efg_onsite 83 !End of the abilint section 84 85 implicit none 86 87 !Arguments ------------------------------------ 88 !scalars 89 integer,intent(in) :: my_natom,natom,nsym,ntypat 90 integer,optional,intent(in) :: comm_atom 91 type(pawang_type),intent(in) :: pawang 92 !arrays 93 integer,intent(in) :: symrel(3,3,nsym) 94 integer,optional,target,intent(in) :: mpi_atmtab(:) 95 real(dp),intent(in) :: rprimd(3,3),tnons(3,nsym),xred(3,natom) 96 real(dp),intent(out) :: efg(3,3,natom) 97 type(paw_an_type),intent(in) :: paw_an(my_natom) 98 type(pawrad_type),intent(in) :: pawrad(ntypat) 99 type(pawrhoij_type),intent(in) :: pawrhoij(my_natom) 100 type(pawtab_type),intent(in) :: pawtab(ntypat) 101 102 !Local variables------------------------------- 103 !scalars 104 integer :: cplex,iatom,iatom_tot,ictr,ierr,imesh_size,ispden,itypat 105 integer :: lm,lm_size,lnspden,mesh_size,my_comm_atom,nzlmopt,nspden 106 integer :: opt_compch,opt_dens,opt_l,opt_print 107 logical :: my_atmtab_allocated,paral_atom 108 real(dp) :: c1,c2,c3,compch_sph,intg 109 !arrays 110 integer,pointer :: my_atmtab(:) 111 logical,allocatable :: lmselectin(:),lmselectout(:) 112 real(dp),allocatable :: ff(:),nhat1(:,:,:),rho1(:,:,:),trho1(:,:,:) 113 114 ! ************************************************************************ 115 116 DBG_ENTER("COLL") 117 118 efg(:,:,:) = zero 119 120 !Set up parallelism over atoms 121 paral_atom=(present(comm_atom).and.(my_natom/=natom)) 122 nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab 123 my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom 124 call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,& 125 & my_natom_ref=my_natom) 126 127 !the following factors arise in expanding the angular dependence of the electric field gradient tensor in 128 !terms of real spherical harmonics. The real spherical harmonics are as in the routine initylmr.F90; see 129 !in particular also http://www.unioviedo.es/qcg/art/Theochem419-19-ov-BF97-rotation-matrices.pdf 130 c1 = sqrt(16.0*pi/5.0) 131 c2 = sqrt(4.0*pi/5.0) 132 c3 = sqrt(12.0*pi/5.0) 133 134 !loop over atoms in cell 135 do iatom = 1, my_natom 136 iatom_tot=iatom;if (paral_atom) iatom_tot=my_atmtab(iatom) 137 itypat=pawrhoij(iatom)%itypat 138 139 lm_size = paw_an(iatom)%lm_size 140 if (lm_size < 5) cycle ! if lm_size < 5 then the on-site densities for this atom have no L=2 component 141 ! and therefore nothing to contribute to the on-site electric field gradient 142 143 mesh_size=pawtab(itypat)%mesh_size 144 ABI_ALLOCATE(ff,(mesh_size)) 145 146 cplex = pawrhoij(iatom)%cplex 147 nspden = pawrhoij(iatom)%nspden 148 ABI_ALLOCATE(lmselectin,(lm_size)) 149 ABI_ALLOCATE(lmselectout,(lm_size)) 150 lmselectin = .true. ! compute all moments of densities 151 nzlmopt = -1 152 opt_compch = 0 153 compch_sph = zero 154 opt_dens = 0 ! compute all densities 155 opt_l = -1 ! all moments contribute 156 opt_print = 0 ! do not print out moments 157 158 ABI_ALLOCATE(nhat1,(cplex*mesh_size,lm_size,nspden)) 159 ABI_ALLOCATE(rho1,(cplex*mesh_size,lm_size,nspden)) 160 ABI_ALLOCATE(trho1,(cplex*mesh_size,lm_size,nspden)) 161 162 ! loop over spin components 163 ! nspden = 1: just a single component 164 ! nspden = 2: two components, loop over them 165 ! nspden = 4: total is in component 1, only one of interest 166 if ( nspden == 2 ) then 167 lnspden = 2 168 else 169 lnspden = 1 170 end if 171 do ispden=1,lnspden 172 173 ! construct multipole expansion of on-site charge densities for this atom 174 call pawdensities(compch_sph,cplex,iatom_tot,lmselectin,lmselectout,lm_size,& 175 & nhat1,nspden,nzlmopt,opt_compch,opt_dens,opt_l,opt_print,& 176 & pawang,0,pawrad(itypat),pawrhoij(iatom),pawtab(itypat),& 177 & rho1,trho1) 178 179 do lm = 5, 9 ! loop on L=2 components of multipole expansion 180 181 if(.not. lmselectout(lm)) cycle ! skip moments that contributes zero 182 183 ! the following is r^2*(n1-tn1-nhat)/r^3 for this multipole moment 184 ! use only the real part of the density in case of cplex==2 185 do imesh_size = 2, mesh_size 186 ictr = cplex*(imesh_size - 1) + 1 187 ff(imesh_size)=(rho1(ictr,lm,ispden)-trho1(ictr,lm,ispden)-& 188 & nhat1(ictr,lm,ispden))/& 189 & pawrad(itypat)%rad(imesh_size) 190 end do 191 call pawrad_deducer0(ff,mesh_size,pawrad(itypat)) 192 call simp_gen(intg,ff,pawrad(itypat)) 193 select case (lm) 194 case (5) ! S_{2,-2} 195 efg(1,2,iatom_tot) = efg(1,2,iatom_tot) - c3*intg ! xy case 196 case (6) ! S_{2,-1} 197 efg(2,3,iatom_tot) = efg(2,3,iatom_tot) - c3*intg ! yz case 198 case (7) ! S_{2, 0} 199 efg(1,1,iatom_tot) = efg(1,1,iatom_tot) + c2*intg ! xx case 200 efg(2,2,iatom_tot) = efg(2,2,iatom_tot) + c2*intg ! yy case 201 efg(3,3,iatom_tot) = efg(3,3,iatom_tot) - c1*intg ! zz case 202 case (8) ! S_{2,+1} 203 efg(1,3,iatom_tot) = efg(1,3,iatom_tot) - c3*intg ! xz case 204 case (9) ! S_{2,+2} 205 efg(1,1,iatom_tot) = efg(1,1,iatom_tot) - c3*intg ! xx case 206 efg(2,2,iatom_tot) = efg(2,2,iatom_tot) + c3*intg ! yy case 207 end select 208 209 end do ! end loop over LM components with L=2 210 211 end do ! Loop on spin components 212 213 ! Symmetrization of EFG 214 efg(2,1,iatom_tot) = efg(1,2,iatom_tot) 215 efg(3,1,iatom_tot) = efg(1,3,iatom_tot) 216 efg(3,2,iatom_tot) = efg(2,3,iatom_tot) 217 218 ABI_DEALLOCATE(lmselectin) 219 ABI_DEALLOCATE(lmselectout) 220 ABI_DEALLOCATE(ff) 221 ABI_DEALLOCATE(nhat1) 222 ABI_DEALLOCATE(rho1) 223 ABI_DEALLOCATE(trho1) 224 225 end do ! Loop on atoms 226 227 !Reduction in case of parallelisation over atoms 228 if (paral_atom) then 229 call xmpi_sum(efg,my_comm_atom,ierr) 230 end if 231 232 ! symmetrize tensor at each atomic site using point symmetry operations 233 do iatom = 1, natom 234 call matpointsym(iatom,efg(:,:,iatom),natom,nsym,rprimd,symrel,tnons,xred) 235 end do 236 237 !Destroy atom table used for parallelism 238 call free_my_atmtab(my_atmtab,my_atmtab_allocated) 239 240 DBG_EXIT("COLL") 241 242 end subroutine make_efg_onsite