TABLE OF CONTENTS
ABINIT/m_vcoul/m_cutoff_slab [ Modules ]
NAME
m_cutoff_slab
FUNCTION
SOURCE
9 #if defined HAVE_CONFIG_H 10 #include "config.h" 11 #endif 12 13 #include "abi_common.h" 14 15 module m_cutoff_slab 16 17 use defs_basis 18 use m_abicore 19 use m_errors 20 21 use m_fstrings, only : sjoin, itoa 22 23 implicit none 24 25 private
m_vcoul/cutoff_slab [ Functions ]
[ Top ] [ m_vcoul ] [ Functions ]
NAME
cutoff_slab
FUNCTION
Calculate the Fourier components of an effective Coulomb interaction within a slab of thickness 2*rcut which is symmetric with respect to the xy plane. In this implementation rcut=L_z/2 where L_z is the periodicity along z
INPUTS
qpt(3)=q-point ng=Number of G vectors. gvec(3,ng)=G vectors in reduced coordinates. gprimd(3,3)=Dimensional primitive translations in reciprocal space ($\textrm{bohr}^{-1}$). gmet(3,3)=Metric in reciprocal space.
OUTPUT
vc_cut(ng)=Fourier components of the effective Coulomb interaction.
NOTES
The Fourier expression for an interaction truncated along the z-direction (i.e non-zero only if |z|<R) is: vc(q.G) = 4pi/|q+G|^2 * [ 1 + e^{-((q+G)_xy)*R} * ( (q_z+G_z)/(q+G)_xy * sin((q_z+G_z)R) - - cos((q_z+G_Z)R)) ] (1) Equation (1) diverges when q_xy+G_xy --> 0 for any non zero q_z+G_z However if we choose R=L/2, where L defines the periodicity along z, and we limit ourselves to consider q-points with q_z==0, then sin((q_z+G_z)R)=sin(G_Z 2pi/L)=0 for every G. Under these assumptions we obtain v(q,G) = 4pi/|q+G|^2 [1-e^{-(q+G)_xy*L/2}\cos((q_z+G_z)R)] which is always finite when G_z /=0 while it diverges as 4piR/(q+G)_xy as (q+G)_xy -->0 but only in the x-y plane.
SOURCE
75 subroutine cutoff_slab(qpt, ng, gvec, gprimd, rcut, boxcenter, pdir, alpha, vc_cut, method) 76 77 !Arguments ------------------------------------ 78 !scalars 79 integer,intent(in) :: method,ng 80 real(dp),intent(in) :: rcut 81 !arrays 82 integer,intent(in) :: gvec(3,ng),pdir(3) 83 real(dp),intent(in) :: alpha(3),boxcenter(3),gprimd(3,3),qpt(3) 84 real(dp),intent(out) :: vc_cut(ng) 85 86 !Local variables------------------------------- 87 !scalars 88 integer :: ig,igs 89 real(dp),parameter :: SMALL=tol4 !@WC: was tol6 90 real(dp) :: qpg2,qpg_para,qpg_perp 91 character(len=500) :: msg 92 !arrays 93 real(dp) :: b1(3),b2(3),b3(3),gcart(3),qc(3),qpg(3) 94 95 ! ************************************************************************* 96 97 ABI_UNUSED(pdir) 98 ABI_UNUSED(boxcenter) 99 100 ! From reduced to cartesian coordinates. 101 b1(:)=two_pi*gprimd(:,1) 102 b2(:)=two_pi*gprimd(:,2) 103 b3(:)=two_pi*gprimd(:,3) 104 105 qc = b1*qpt(1) + b2*qpt(2) + b3*qpt(3) 106 107 ! Different approaches according to method 108 vc_cut = zero 109 110 select case (method) 111 112 case (1) 113 ! Beigi's expression. 114 ! q-points with non-zero component along the z-axis are not allowed if 115 ! the simplified Eq.1 for the Coulomb interaction is used. 116 if (ANY(ABS(qc) > SMALL)) then 117 write(std_out,*)qc 118 write(msg,'(5a)')& 119 'Found q-points with non-zero component along non-periodic direction ',ch10,& 120 'This is not allowed, see Notes in cutoff_slab.F90 ',ch10,& 121 'ACTION: Modify the q-point sampling ' 122 ABI_ERROR(msg) 123 end if 124 125 ! Calculate truncated Coulomb interaction for a infinite surface 126 ! supposing input q-points are different from zero. 127 igs=1; if (SQRT(DOT_PRODUCT(qc,qc))<tol16) igs=2 ! avoid (q=0, G=0) 128 do ig=igs,ng 129 gcart(:) = b1(:)*gvec(1,ig)+b2(:)*gvec(2,ig)+b3(:)*gvec(3,ig) 130 qpg(:) = qc(:) + gcart(:) 131 qpg2 = DOT_PRODUCT(qpg(:),qpg(:)) 132 qpg_para = SQRT(qpg(1)**2+qpg(2)**2) ; qpg_perp=qpg(3) 133 vc_cut(ig) = four_pi/qpg2*(one-EXP(-qpg_para*rcut)*COS(qpg_perp*rcut)) 134 end do 135 136 case (2) 137 ! Rozzi's method 138 ABI_ERROR("Work in progress") 139 ABI_UNUSED(alpha) ! just to keep alpha as an argument 140 !alpha=?? ; ap1sqrt=SQRT(one+alpha**2) 141 do ig=1,ng 142 gcart(:) = b1(:)*gvec(1,ig)+b2(:)*gvec(2,ig)+b3(:)*gvec(3,ig) 143 qpg(:) = qc(:) + gcart(:) 144 qpg2 =DOT_PRODUCT(qpg(:),qpg(:)) 145 qpg_para=SQRT(qpg(1)**2+qpg(2)**2) ; qpg_perp =qpg(3) 146 if (qpg_para>SMALL) then 147 vc_cut(ig)=four_pi/qpg2*(one+EXP(-qpg_para*rcut)*(qpg_perp/qpg_para*SIN(qpg_perp*rcut)-COS(qpg_perp*rcut))) 148 else 149 if (ABS(qpg_perp)>SMALL) then 150 vc_cut(ig)=four_pi/qpg_perp**2*(one-COS(qpg_perp*rcut)-qpg_perp*rcut*SIN(qpg_perp*rcut)) ! & 151 ! contribution due to finite slab 152 ! + 8*rcut*SIN(qpg_perp*rcut)/qpg_perp*LOG((alpha+ap1sqrt)*(one+ap1sqrt)/alpha) 153 else 154 vc_cut(ig)=-two_pi*rcut**2 155 end if 156 end if 157 end do !ig 158 159 case default 160 ABI_BUG(sjoin('Wrong value for method:', itoa(method))) 161 end select 162 163 end subroutine cutoff_slab