Line data Source code
1 : !!****m* ABINIT/m_vcoul/m_cutoff_slab
2 : !! NAME
3 : !! m_cutoff_slab
4 : !!
5 : !! FUNCTION
6 : !!
7 : !! SOURCE
8 :
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
26 : !!***
27 :
28 : public :: cutoff_slab
29 :
30 : !integer,public,parameter :: SURFACE_BEIGI = 1
31 : !integer,public,parameter :: SURFACE_ROZZI = 2
32 :
33 : CONTAINS !========================================================================================
34 : !!***
35 :
36 : !----------------------------------------------------------------------
37 :
38 : !!****f* m_vcoul/cutoff_slab
39 : !! NAME
40 : !! cutoff_slab
41 : !!
42 : !! FUNCTION
43 : !! Calculate the Fourier components of an effective Coulomb interaction
44 : !! within a slab of thickness 2*rcut which is symmetric with respect to the xy plane.
45 : !! In this implementation rcut=L_z/2 where L_z is the periodicity along z
46 : !!
47 : !! INPUTS
48 : !! qpt(3)=q-point
49 : !! ng=Number of G vectors.
50 : !! gvec(3,ng)=G vectors in reduced coordinates.
51 : !! gprimd(3,3)=Dimensional primitive translations in reciprocal space ($\textrm{bohr}^{-1}$).
52 : !! gmet(3,3)=Metric in reciprocal space.
53 : !!
54 : !! OUTPUT
55 : !! vc_cut(ng)=Fourier components of the effective Coulomb interaction.
56 : !!
57 : !! NOTES
58 : !! The Fourier expression for an interaction truncated along the z-direction (i.e non-zero only if |z|<R) is:
59 : !!
60 : !! 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) -
61 : !! - cos((q_z+G_Z)R)) ] (1)
62 : !!
63 : !! Equation (1) diverges when q_xy+G_xy --> 0 for any non zero q_z+G_z
64 : !! However if we choose R=L/2, where L defines the periodicity along z,
65 : !! 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.
66 : !! Under these assumptions we obtain
67 : !!
68 : !! v(q,G) = 4pi/|q+G|^2 [1-e^{-(q+G)_xy*L/2}\cos((q_z+G_z)R)]
69 : !!
70 : !! which is always finite when G_z /=0 while it diverges as 4piR/(q+G)_xy as (q+G)_xy -->0
71 : !! but only in the x-y plane.
72 : !!
73 : !! SOURCE
74 :
75 0 : 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 0 : b1(:)=two_pi*gprimd(:,1)
102 0 : b2(:)=two_pi*gprimd(:,2)
103 0 : b3(:)=two_pi*gprimd(:,3)
104 :
105 0 : qc = b1*qpt(1) + b2*qpt(2) + b3*qpt(3)
106 :
107 : ! Different approaches according to method
108 0 : vc_cut = zero
109 :
110 0 : 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 0 : if (ANY(ABS(qc) > SMALL)) then
117 0 : write(std_out,*)qc
118 : write(msg,'(5a)')&
119 0 : 'Found q-points with non-zero component along non-periodic direction ',ch10,&
120 0 : 'This is not allowed, see Notes in cutoff_slab.F90 ',ch10,&
121 0 : 'ACTION: Modify the q-point sampling '
122 0 : 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 0 : igs=1; if (NORM2(qc)<tol16) igs=2 ! avoid (q=0, G=0)
128 0 : do ig=igs,ng
129 0 : gcart(:) = b1(:)*gvec(1,ig)+b2(:)*gvec(2,ig)+b3(:)*gvec(3,ig)
130 0 : qpg(:) = qc(:) + gcart(:)
131 0 : qpg2 = DOT_PRODUCT(qpg(:),qpg(:))
132 0 : qpg_para = SQRT(qpg(1)**2+qpg(2)**2) ; qpg_perp=qpg(3)
133 0 : 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 0 : ABI_ERROR("Work in progress")
139 : ABI_UNUSED(alpha) ! just to keep alpha as an argument
140 : !alpha=?? ; ap1sqrt=SQRT(one+alpha**2)
141 0 : do ig=1,ng
142 0 : gcart(:) = b1(:)*gvec(1,ig)+b2(:)*gvec(2,ig)+b3(:)*gvec(3,ig)
143 0 : qpg(:) = qc(:) + gcart(:)
144 0 : qpg2 =DOT_PRODUCT(qpg(:),qpg(:))
145 0 : qpg_para=SQRT(qpg(1)**2+qpg(2)**2) ; qpg_perp =qpg(3)
146 0 : if (qpg_para>SMALL) then
147 0 : 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 0 : if (ABS(qpg_perp)>SMALL) then
150 0 : 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 0 : vc_cut(ig)=-two_pi*rcut**2
155 : end if
156 : end if
157 : end do !ig
158 :
159 : case default
160 0 : ABI_BUG(sjoin('Wrong value for method:', itoa(method)))
161 : end select
162 :
163 0 : end subroutine cutoff_slab
164 : !!***
165 :
166 : end module m_cutoff_slab
167 : !!***
|
