TABLE OF CONTENTS

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

      vtorho

CHILDREN

SOURCE

 42 #if defined HAVE_CONFIG_H
 43 #include "config.h"
 44 #endif
 45 
 46 #include "abi_common.h"
 47 
 48 
 49 subroutine mknucdipmom_k(gmet,kg,kpt,natom,nucdipmom,nucdipmom_k,npw,rprimd,ucvol,xred)
 50 
 51  use defs_basis
 52  use m_profiling_abi
 53 
 54 !This section has been created automatically by the script Abilint (TD).
 55 !Do not modify the following lines by hand.
 56 #undef ABI_FUNC
 57 #define ABI_FUNC 'mknucdipmom_k'
 58 !End of the abilint section
 59 
 60  implicit none
 61 
 62 !Arguments ------------------------------------
 63  !scalars
 64  integer,intent(in) :: natom,npw
 65  real(dp),intent(in) :: ucvol
 66 
 67  !arrays
 68  integer,intent(in) :: kg(3,npw)
 69  real(dp),intent(in) :: gmet(3,3),kpt(3),nucdipmom(3,natom),rprimd(3,3),xred(3,natom)
 70  complex(dpc),intent(out) :: nucdipmom_k(npw*(npw+1)/2)
 71  
 72 !Local variables-------------------------------
 73 !scalars
 74  integer :: atom_nd_tot,col,iatom,ndp_index,row
 75  real(dp) :: crossfac,dg2,permeability,phasefac
 76  complex(dpc) :: cpermfac,cphasefac
 77  !arrays
 78  integer :: atom_nd(natom)
 79  real(dp) :: cprod(3),cprod_cart(3),dgp_red(3), gpk_red(3)
 80 
 81 ! *************************************************************************
 82  !
 83 
 84 ! magnetic permeability mu_0/four_pi in atomic units
 85 ! this constant is also used in m_pawdij.F90/pawdijnd, if you change it here,
 86 ! change it there also for consistency
 87  permeability=5.325135453D-5
 88  ! will need 4*pi*i*(\mu_0/four\pi)
 89  cpermfac = CMPLX(zero,four_pi*permeability)
 90 
 91  ! make list of atoms with non-zero nuclear magnetic dipoles
 92  atom_nd_tot = 0
 93  do iatom = 1, natom
 94    if(any(abs(nucdipmom(:,iatom))>tol8)) then
 95      atom_nd_tot = atom_nd_tot + 1
 96      atom_nd(atom_nd_tot) = iatom
 97    end if
 98  end do
 99 
100  ndp_index = 0
101  do col=1,npw ! enumerate plane waves G
102     ! form k + G at this k point for current plane wave (this is the ket |k+G> )
103     ! in reduced coordinates
104    gpk_red(1)=dble(kg(1,col))+kpt(1)
105    gpk_red(2)=dble(kg(2,col))+kpt(2)
106    gpk_red(3)=dble(kg(3,col))+kpt(3)
107 
108    do row=col,npw ! enumerate lower diagonal from 1 to G
109       ! index of the current matrix element, in lower triangular packed storage
110       ! "packed sequentially, column by column"
111      ndp_index = ndp_index + 1
112      nucdipmom_k(ndp_index) = czero
113      
114       ! form G-G' = \Delta G at this k pt (this is the bra <k+G'| )
115       ! in reduced coordinates
116      dgp_red(1)=dble(kg(1,col)-kg(1,row))
117      dgp_red(2)=dble(kg(2,col)-kg(2,row))
118      dgp_red(3)=dble(kg(3,col)-kg(3,row))
119 
120       ! compute |\Delta G|^2
121       ! must use gmet metric because G's are in reduced coords in reciprocal space
122      dg2 = DOT_PRODUCT(dgp_red,MATMUL(gmet,dgp_red))
123       ! if \Delta G = 0, Hamiltonian term is zero and move on to next one
124      if (abs(dg2)<tol8) then
125        nucdipmom_k(ndp_index)=czero
126        cycle
127      end if
128 
129       ! compute cross product \Delta G \times (k + G)
130       ! notice that \Delta G and (k + G) are in reduced coords in reciprocal space
131      cprod(1) = dgp_red(2)*gpk_red(3) - dgp_red(3)*gpk_red(2)
132      cprod(2) = dgp_red(3)*gpk_red(1) - dgp_red(1)*gpk_red(3)
133      cprod(3) = dgp_red(1)*gpk_red(2) - dgp_red(2)*gpk_red(1)
134 
135       ! proper cross product must account for reduced coords as follows:
136       ! gprimd*dgp \times gprimd*gpk = (det gprimd)*(gprimd^{-1,T})*(dgp \times gpk)
137       ! = rprimd * (dgp \times gpk)/ucvol
138       ! final vector also includes the division by |\Delta G|^2
139      cprod_cart = MATMUL(rprimd,cprod)/(ucvol*dg2)
140 
141       ! loop over the atoms with non-zero nuclear dipoles
142       ! phase factors exp(i*\Delta G*I) where I is ion position,
143       ! might be retrievable from ph1d, need to check
144      do iatom = 1, atom_nd_tot
145        phasefac = two_pi*DOT_PRODUCT(dgp_red,xred(:,atom_nd(iatom)))
146        cphasefac = CMPLX(cos(phasefac),sin(phasefac))
147        crossfac = DOT_PRODUCT(nucdipmom(:,iatom),cprod_cart)
148        nucdipmom_k(ndp_index) = nucdipmom_k(ndp_index) + cpermfac*cphasefac*crossfac
149      end do ! end loop over atoms with nonzero dipoles
150 
151    end do ! end loop over G' = G to npw
152 
153  end do ! end loop over G = 1 to npw
154 
155 end subroutine mknucdipmom_k