Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : MODULE surface_dipole
10 :
11 : USE cell_types, ONLY: cell_type
12 : USE cp_control_types, ONLY: dft_control_type
13 : USE kahan_sum, ONLY: accurate_sum
14 : USE kinds, ONLY: dp
15 : USE mathconstants, ONLY: pi
16 : USE physcon, ONLY: bohr
17 : USE pw_env_types, ONLY: pw_env_get,&
18 : pw_env_type
19 : USE pw_grid_types, ONLY: PW_MODE_LOCAL
20 : USE pw_methods, ONLY: pw_axpy,&
21 : pw_integral_ab,&
22 : pw_scale,&
23 : pw_transfer,&
24 : pw_zero
25 : USE pw_pool_types, ONLY: pw_pool_p_type,&
26 : pw_pool_type
27 : USE pw_types, ONLY: pw_c1d_gs_type,&
28 : pw_r3d_rs_type
29 : USE qs_energy_types, ONLY: qs_energy_type
30 : USE qs_environment_types, ONLY: get_qs_env,&
31 : qs_environment_type
32 : USE qs_rho_types, ONLY: qs_rho_get,&
33 : qs_rho_type
34 : USE qs_subsys_types, ONLY: qs_subsys_type
35 : #include "./base/base_uses.f90"
36 :
37 : IMPLICIT NONE
38 :
39 : PRIVATE
40 :
41 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'surface_dipole'
42 :
43 : PUBLIC :: calc_dipsurf_potential
44 :
45 : CONTAINS
46 :
47 : ! **************************************************************************************************
48 : !> \brief compute the surface dipole and the correction to the hartree potential
49 : !> \param qs_env the qs environment
50 : !> \param energy ...
51 : !> \par History
52 : !> 01.2014 created [MI]
53 : !> \author MI
54 : !> \author Soumya Ghosh added SURF_DIP_POS 19.11.2018
55 : ! **************************************************************************************************
56 :
57 110 : SUBROUTINE calc_dipsurf_potential(qs_env, energy)
58 :
59 : TYPE(qs_environment_type), POINTER :: qs_env
60 : TYPE(qs_energy_type), POINTER :: energy
61 :
62 : CHARACTER(len=*), PARAMETER :: routineN = 'calc_dipsurf_potential'
63 :
64 : INTEGER :: handle, i, idir_surfdip, ilayer_min, &
65 : ilow, irho, ispin, isurf, iup, jsurf, &
66 : width
67 : INTEGER, DIMENSION(3) :: ngrid
68 110 : INTEGER, DIMENSION(:, :), POINTER :: bo
69 : REAL(dp) :: cutoff, dh(3, 3), dip_fac, dip_hh, &
70 : dsurf, height_min, hh, pos_surf_dip, &
71 : rhoav_min, surfarea, vdip, vdip_fac
72 110 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: rhoavsurf
73 : TYPE(cell_type), POINTER :: cell
74 : TYPE(dft_control_type), POINTER :: dft_control
75 : TYPE(pw_c1d_gs_type), POINTER :: rho0_s_gs, rho_core, rhoz_cneo_s_gs
76 : TYPE(pw_env_type), POINTER :: pw_env
77 110 : TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools
78 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
79 : TYPE(pw_r3d_rs_type) :: vdip_r, wf_r
80 110 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r
81 : TYPE(pw_r3d_rs_type), POINTER :: v_hartree_rspace
82 : TYPE(qs_rho_type), POINTER :: rho
83 : TYPE(qs_subsys_type), POINTER :: subsys
84 :
85 110 : CALL timeset(routineN, handle)
86 110 : NULLIFY (cell, dft_control, rho, pw_env, auxbas_pw_pool, &
87 110 : pw_pools, subsys, v_hartree_rspace, rho_r, rhoz_cneo_s_gs)
88 :
89 : CALL get_qs_env(qs_env, &
90 : dft_control=dft_control, &
91 : rho=rho, &
92 : rho_core=rho_core, &
93 : rho0_s_gs=rho0_s_gs, &
94 : rhoz_cneo_s_gs=rhoz_cneo_s_gs, &
95 : cell=cell, &
96 : pw_env=pw_env, &
97 : subsys=subsys, &
98 110 : v_hartree_rspace=v_hartree_rspace)
99 :
100 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
101 110 : pw_pools=pw_pools)
102 110 : CALL auxbas_pw_pool%create_pw(wf_r)
103 110 : CALL auxbas_pw_pool%create_pw(vdip_r)
104 :
105 110 : IF (dft_control%qs_control%gapw) THEN
106 0 : IF (dft_control%qs_control%gapw_control%nopaw_as_gpw) THEN
107 0 : CALL pw_axpy(rho_core, rho0_s_gs)
108 0 : IF (ASSOCIATED(rhoz_cneo_s_gs)) THEN
109 0 : CALL pw_axpy(rhoz_cneo_s_gs, rho0_s_gs)
110 : END IF
111 0 : CALL pw_transfer(rho0_s_gs, wf_r)
112 0 : CALL pw_axpy(rho_core, rho0_s_gs, -1.0_dp)
113 0 : IF (ASSOCIATED(rhoz_cneo_s_gs)) THEN
114 0 : CALL pw_axpy(rhoz_cneo_s_gs, rho0_s_gs, -1.0_dp)
115 : END IF
116 : ELSE
117 0 : IF (ASSOCIATED(rhoz_cneo_s_gs)) THEN
118 0 : CALL pw_axpy(rhoz_cneo_s_gs, rho0_s_gs)
119 : END IF
120 0 : CALL pw_transfer(rho0_s_gs, wf_r)
121 0 : IF (ASSOCIATED(rhoz_cneo_s_gs)) THEN
122 0 : CALL pw_axpy(rhoz_cneo_s_gs, rho0_s_gs, -1.0_dp)
123 : END IF
124 : END IF
125 : ELSE
126 110 : CALL pw_transfer(rho_core, wf_r)
127 : END IF
128 110 : CALL qs_rho_get(rho, rho_r=rho_r)
129 240 : DO ispin = 1, dft_control%nspins
130 240 : CALL pw_axpy(rho_r(ispin), wf_r)
131 : END DO
132 :
133 440 : ngrid(1:3) = wf_r%pw_grid%npts(1:3)
134 110 : idir_surfdip = dft_control%dir_surf_dip
135 :
136 110 : width = 4
137 :
138 440 : DO i = 1, 3
139 440 : IF (i /= idir_surfdip) THEN
140 220 : IF (ABS(wf_r%pw_grid%dh(idir_surfdip, i)) > 1.E-7_dp) THEN
141 : ! stop surface dipole defined only for slab perpendigular to one of the Cartesian axis
142 : CALL cp_abort(__LOCATION__, &
143 0 : "Dipole correction only for surface perpendicular to one Cartesian axis")
144 : ! not properly general, we assume that vectors A, B, and C are along x y and z respectively,
145 : ! in the ortorhombic cell, but in principle it does not need to be this way, importan
146 : ! is that the cell angles are 90 degrees.
147 : END IF
148 : END IF
149 : END DO
150 :
151 110 : ilow = wf_r%pw_grid%bounds(1, idir_surfdip)
152 110 : iup = wf_r%pw_grid%bounds(2, idir_surfdip)
153 :
154 330 : ALLOCATE (rhoavsurf(ilow:iup))
155 27438 : rhoavsurf = 0.0_dp
156 :
157 110 : bo => wf_r%pw_grid%bounds_local
158 1430 : dh = wf_r%pw_grid%dh
159 :
160 110 : CALL pw_scale(wf_r, wf_r%pw_grid%vol)
161 110 : IF (idir_surfdip == 3) THEN
162 56 : isurf = 1
163 56 : jsurf = 2
164 :
165 13784 : DO i = bo(1, 3), bo(2, 3)
166 13784 : rhoavsurf(i) = accurate_sum(wf_r%array(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), i))
167 : END DO
168 :
169 54 : ELSEIF (idir_surfdip == 2) THEN
170 0 : isurf = 3
171 0 : jsurf = 1
172 :
173 0 : DO i = bo(1, 2), bo(2, 2)
174 0 : rhoavsurf(i) = accurate_sum(wf_r%array(bo(1, 1):bo(2, 1), i, bo(1, 3):bo(2, 3)))
175 : END DO
176 : ELSE
177 54 : isurf = 2
178 54 : jsurf = 3
179 :
180 6854 : DO i = bo(1, 1), bo(2, 1)
181 6854 : rhoavsurf(i) = accurate_sum(wf_r%array(i, bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
182 : END DO
183 : END IF
184 110 : CALL pw_scale(wf_r, 1.0_dp/wf_r%pw_grid%vol)
185 27438 : rhoavsurf = rhoavsurf/wf_r%pw_grid%vol
186 :
187 : surfarea = cell%hmat(isurf, isurf)*cell%hmat(jsurf, jsurf) - &
188 110 : cell%hmat(isurf, jsurf)*cell%hmat(jsurf, isurf)
189 110 : dsurf = surfarea/REAL(ngrid(isurf)*ngrid(jsurf), dp)
190 :
191 110 : IF (wf_r%pw_grid%para%mode /= PW_MODE_LOCAL) THEN
192 110 : CALL wf_r%pw_grid%para%group%sum(rhoavsurf)
193 : END IF
194 27438 : rhoavsurf(ilow:iup) = dsurf*rhoavsurf(ilow:iup)
195 :
196 : ! locate where the vacuum is, and set the reference point for the calculation of the dipole
197 27438 : rhoavsurf(ilow:iup) = rhoavsurf(ilow:iup)/surfarea
198 : ! Note: rhosurf has the same dimension as rho
199 110 : IF (dft_control%pos_dir_surf_dip < 0.0_dp) THEN
200 3200 : ilayer_min = ilow - 1 + MINLOC(ABS(rhoavsurf(ilow:iup)), 1)
201 : ELSE
202 78 : pos_surf_dip = dft_control%pos_dir_surf_dip*bohr
203 78 : ilayer_min = ilow - 1 + NINT(pos_surf_dip/dh(idir_surfdip, idir_surfdip)) + 1
204 : END IF
205 110 : rhoav_min = ABS(rhoavsurf(ilayer_min))
206 110 : IF (rhoav_min >= 1.E-5_dp) THEN
207 0 : CPABORT(" Dipole correction needs more vacuum space above the surface ")
208 : END IF
209 :
210 110 : height_min = REAL((ilayer_min - ilow), dp)*dh(idir_surfdip, idir_surfdip)
211 :
212 : ! surface dipole form average rhoavsurf
213 : ! \sum_i NjdjNkdkdi rhoav_i (i-imin)di
214 110 : dip_hh = 0.0_dp
215 110 : dip_fac = wf_r%pw_grid%vol*dh(idir_surfdip, idir_surfdip)/REAL(ngrid(idir_surfdip), dp)
216 :
217 27438 : DO i = ilayer_min + 1, ilayer_min + ngrid(idir_surfdip)
218 27328 : hh = REAL((i - ilayer_min), dp)
219 27328 : IF (i > iup) THEN
220 19702 : irho = i - ngrid(idir_surfdip)
221 : ELSE
222 : irho = i
223 : END IF
224 : ! introduce a cutoff function to smoothen the edges
225 27328 : IF (ABS(irho - ilayer_min) > width) THEN
226 : cutoff = 1.0_dp
227 : ELSE
228 990 : cutoff = ABS(SIN(0.5_dp*pi*REAL(ABS(irho - ilayer_min), dp)/REAL(width, dp)))
229 : END IF
230 27438 : dip_hh = dip_hh + rhoavsurf(irho)*hh*dip_fac*cutoff
231 : END DO
232 :
233 110 : DEALLOCATE (rhoavsurf)
234 : ! for printing purposes [SGh]
235 110 : qs_env%surface_dipole_moment = dip_hh/bohr
236 :
237 : ! Calculation of the dipole potential as a function of the perpendicular coordinate
238 110 : CALL pw_zero(vdip_r)
239 110 : vdip_fac = dip_hh*4.0_dp*pi
240 :
241 27438 : DO i = ilayer_min + 1, ilayer_min + ngrid(idir_surfdip)
242 27328 : hh = REAL((i - ilayer_min), dp)*dh(idir_surfdip, idir_surfdip)
243 : vdip = vdip_fac*(-0.5_dp + (hh/cell%hmat(idir_surfdip, idir_surfdip)))* &
244 27328 : v_hartree_rspace%pw_grid%dvol/surfarea
245 27328 : IF (i > iup) THEN
246 19702 : irho = i - ngrid(idir_surfdip)
247 : ELSE
248 : irho = i
249 : END IF
250 : ! introduce a cutoff function to smoothen the edges
251 27328 : IF (ABS(irho - ilayer_min) > width) THEN
252 : cutoff = 1.0_dp
253 : ELSE
254 990 : cutoff = ABS(SIN(0.5_dp*pi*REAL(ABS(irho - ilayer_min), dp)/REAL(width, dp)))
255 : END IF
256 27328 : vdip = vdip*cutoff
257 :
258 27438 : IF (idir_surfdip == 3) THEN
259 : vdip_r%array(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), irho) = &
260 17343264 : vdip_r%array(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), irho) + vdip
261 13600 : ELSEIF (idir_surfdip == 2) THEN
262 0 : IF (irho >= bo(1, 2) .AND. irho <= bo(2, 2)) THEN
263 : vdip_r%array(bo(1, 1):bo(2, 1), irho, bo(1, 3):bo(2, 3)) = &
264 0 : vdip_r%array(bo(1, 1):bo(2, 1), irho, bo(1, 3):bo(2, 3)) + vdip
265 : END IF
266 : ELSE
267 13600 : IF (irho >= bo(1, 1) .AND. irho <= bo(2, 1)) THEN
268 : vdip_r%array(irho, bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)) = &
269 26378320 : vdip_r%array(irho, bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)) + vdip
270 : END IF
271 : END IF
272 :
273 : END DO
274 :
275 : ! Dipole correction contribution to the energy
276 110 : energy%surf_dipole = 0.5_dp*pw_integral_ab(vdip_r, wf_r, just_sum=.TRUE.)
277 :
278 : ! Add the dipole potential to the hartree potential on the realspace grid
279 110 : CALL pw_axpy(vdip_r, v_hartree_rspace)
280 :
281 110 : CALL auxbas_pw_pool%give_back_pw(wf_r)
282 110 : CALL auxbas_pw_pool%give_back_pw(vdip_r)
283 :
284 110 : CALL timestop(handle)
285 :
286 110 : END SUBROUTINE calc_dipsurf_potential
287 :
288 : END MODULE surface_dipole
|
