Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 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 : evolt
18 : USE pw_env_types, ONLY: pw_env_get,&
19 : pw_env_type
20 : USE pw_grid_types, ONLY: PW_MODE_LOCAL
21 : USE pw_methods, ONLY: pw_axpy,&
22 : pw_integral_ab,&
23 : pw_scale,&
24 : pw_transfer,&
25 : pw_zero
26 : USE pw_pool_types, ONLY: pw_pool_p_type,&
27 : pw_pool_type
28 : USE pw_types, ONLY: pw_c1d_gs_type,&
29 : pw_r3d_rs_type
30 : USE qs_energy_types, ONLY: qs_energy_type
31 : USE qs_environment_types, ONLY: get_qs_env,&
32 : qs_environment_type
33 : USE qs_rho_types, ONLY: qs_rho_get,&
34 : qs_rho_type
35 : USE qs_subsys_types, ONLY: qs_subsys_type
36 : #include "./base/base_uses.f90"
37 :
38 : IMPLICIT NONE
39 :
40 : PRIVATE
41 :
42 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'surface_dipole'
43 :
44 : PUBLIC :: calc_dipsurf_potential
45 :
46 : CONTAINS
47 :
48 : ! **************************************************************************************************
49 : !> \brief compute the surface dipole and the correction to the hartree potential
50 : !> \param qs_env the qs environment
51 : !> \param energy ...
52 : !> \par History
53 : !> 01.2014 created [MI]
54 : !> \author MI
55 : !> \author Soumya Ghosh added SURF_DIP_POS 19.11.2018
56 : ! **************************************************************************************************
57 :
58 110 : SUBROUTINE calc_dipsurf_potential(qs_env, energy)
59 :
60 : TYPE(qs_environment_type), POINTER :: qs_env
61 : TYPE(qs_energy_type), POINTER :: energy
62 :
63 : CHARACTER(len=*), PARAMETER :: routineN = 'calc_dipsurf_potential'
64 :
65 : INTEGER :: handle, i, i_above, i_below, &
66 : idir_surfdip, ilayer_min, ilow, irho, &
67 : ispin, isurf, iup, jsurf, width
68 : INTEGER, DIMENSION(3) :: ngrid
69 110 : INTEGER, DIMENSION(:, :), POINTER :: bo
70 : REAL(dp) :: cutoff, dh(3, 3), dip_fac, dip_hh, dsurf, height_min, hh, pos_surf_dip, &
71 : rhoav_min, surfarea, vac_above, vac_below, vdip, vdip_fac
72 : 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 110 : 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 : ELSE IF (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 110 : qs_env%surface_dipole_ref_pos = height_min/bohr
237 :
238 : ! Calculation of the dipole potential as a function of the perpendicular coordinate
239 110 : CALL pw_zero(vdip_r)
240 110 : vdip_fac = dip_hh*4.0_dp*pi
241 :
242 27438 : DO i = ilayer_min + 1, ilayer_min + ngrid(idir_surfdip)
243 27328 : hh = REAL((i - ilayer_min), dp)*dh(idir_surfdip, idir_surfdip)
244 : vdip = vdip_fac*(-0.5_dp + (hh/cell%hmat(idir_surfdip, idir_surfdip)))* &
245 27328 : v_hartree_rspace%pw_grid%dvol/surfarea
246 27328 : IF (i > iup) THEN
247 19702 : irho = i - ngrid(idir_surfdip)
248 : ELSE
249 : irho = i
250 : END IF
251 : ! introduce a cutoff function to smoothen the edges
252 27328 : IF (ABS(irho - ilayer_min) > width) THEN
253 : cutoff = 1.0_dp
254 : ELSE
255 990 : cutoff = ABS(SIN(0.5_dp*pi*REAL(ABS(irho - ilayer_min), dp)/REAL(width, dp)))
256 : END IF
257 27328 : vdip = vdip*cutoff
258 :
259 27438 : IF (idir_surfdip == 3) THEN
260 : vdip_r%array(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), irho) = &
261 17343264 : vdip_r%array(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), irho) + vdip
262 13600 : ELSE IF (idir_surfdip == 2) THEN
263 0 : IF (irho >= bo(1, 2) .AND. irho <= bo(2, 2)) THEN
264 : vdip_r%array(bo(1, 1):bo(2, 1), irho, bo(1, 3):bo(2, 3)) = &
265 0 : vdip_r%array(bo(1, 1):bo(2, 1), irho, bo(1, 3):bo(2, 3)) + vdip
266 : END IF
267 : ELSE
268 13600 : IF (irho >= bo(1, 1) .AND. irho <= bo(2, 1)) THEN
269 : vdip_r%array(irho, bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)) = &
270 26378320 : vdip_r%array(irho, bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)) + vdip
271 : END IF
272 : END IF
273 :
274 : END DO
275 :
276 : ! Dipole correction contribution to the energy
277 110 : energy%surf_dipole = 0.5_dp*pw_integral_ab(vdip_r, wf_r, just_sum=.TRUE.)
278 :
279 : ! Add the dipole potential to the hartree potential on the realspace grid
280 110 : CALL pw_axpy(vdip_r, v_hartree_rspace)
281 :
282 : ! Vacuum level (plane-averaged, corrected Hartree potential) immediately below and above the
283 : ! dipole correction reference plane. Note: v_hartree_rspace does not carry the GAPW one-center
284 : ! (hard/soft) corrections, but those are localized on the atoms and vanish at these sampling
285 : ! points, which by construction lie in the vacuum.
286 110 : i_below = ilayer_min - 1
287 110 : IF (i_below < ilow) i_below = iup
288 110 : i_above = ilayer_min + 1
289 110 : IF (i_above > iup) i_above = ilow
290 :
291 110 : vac_below = 0.0_dp
292 110 : vac_above = 0.0_dp
293 110 : IF (idir_surfdip == 3) THEN
294 56 : IF (i_below >= bo(1, 3) .AND. i_below <= bo(2, 3)) THEN
295 56 : vac_below = accurate_sum(v_hartree_rspace%array(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), i_below))
296 : END IF
297 56 : IF (i_above >= bo(1, 3) .AND. i_above <= bo(2, 3)) THEN
298 56 : vac_above = accurate_sum(v_hartree_rspace%array(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), i_above))
299 : END IF
300 54 : ELSE IF (idir_surfdip == 2) THEN
301 0 : IF (i_below >= bo(1, 2) .AND. i_below <= bo(2, 2)) THEN
302 0 : vac_below = accurate_sum(v_hartree_rspace%array(bo(1, 1):bo(2, 1), i_below, bo(1, 3):bo(2, 3)))
303 : END IF
304 0 : IF (i_above >= bo(1, 2) .AND. i_above <= bo(2, 2)) THEN
305 0 : vac_above = accurate_sum(v_hartree_rspace%array(bo(1, 1):bo(2, 1), i_above, bo(1, 3):bo(2, 3)))
306 : END IF
307 : ELSE
308 54 : IF (i_below >= bo(1, 1) .AND. i_below <= bo(2, 1)) THEN
309 27 : vac_below = accurate_sum(v_hartree_rspace%array(i_below, bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
310 : END IF
311 54 : IF (i_above >= bo(1, 1) .AND. i_above <= bo(2, 1)) THEN
312 27 : vac_above = accurate_sum(v_hartree_rspace%array(i_above, bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
313 : END IF
314 : END IF
315 :
316 110 : IF (wf_r%pw_grid%para%mode /= PW_MODE_LOCAL) THEN
317 110 : CALL wf_r%pw_grid%para%group%sum(vac_below)
318 110 : CALL wf_r%pw_grid%para%group%sum(vac_above)
319 : END IF
320 :
321 110 : qs_env%vacuum_level_below = vac_below/REAL(ngrid(isurf)*ngrid(jsurf), dp)*evolt
322 110 : qs_env%vacuum_level_above = vac_above/REAL(ngrid(isurf)*ngrid(jsurf), dp)*evolt
323 :
324 110 : CALL auxbas_pw_pool%give_back_pw(wf_r)
325 110 : CALL auxbas_pw_pool%give_back_pw(vdip_r)
326 :
327 110 : CALL timestop(handle)
328 :
329 110 : END SUBROUTINE calc_dipsurf_potential
330 :
331 : END MODULE surface_dipole
|
