Branch data 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 : : !!***
|
