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
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)
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 : cell=cell, &
95 : pw_env=pw_env, &
96 : subsys=subsys, &
97 110 : v_hartree_rspace=v_hartree_rspace)
98 :
99 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
100 110 : pw_pools=pw_pools)
101 110 : CALL auxbas_pw_pool%create_pw(wf_r)
102 110 : CALL auxbas_pw_pool%create_pw(vdip_r)
103 :
104 110 : IF (dft_control%qs_control%gapw) THEN
105 0 : IF (dft_control%qs_control%gapw_control%nopaw_as_gpw) THEN
106 0 : CALL pw_axpy(rho_core, rho0_s_gs)
107 0 : CALL pw_transfer(rho0_s_gs, wf_r)
108 0 : CALL pw_axpy(rho_core, rho0_s_gs, -1.0_dp)
109 : ELSE
110 0 : CALL pw_transfer(rho0_s_gs, wf_r)
111 : END IF
112 : ELSE
113 110 : CALL pw_transfer(rho_core, wf_r)
114 : END IF
115 110 : CALL qs_rho_get(rho, rho_r=rho_r)
116 240 : DO ispin = 1, dft_control%nspins
117 240 : CALL pw_axpy(rho_r(ispin), wf_r)
118 : END DO
119 :
120 440 : ngrid(1:3) = wf_r%pw_grid%npts(1:3)
121 110 : idir_surfdip = dft_control%dir_surf_dip
122 :
123 110 : width = 4
124 :
125 440 : DO i = 1, 3
126 440 : IF (i /= idir_surfdip) THEN
127 220 : IF (ABS(wf_r%pw_grid%dh(idir_surfdip, i)) > 1.E-7_dp) THEN
128 : ! stop surface dipole defined only for slab perpendigular to one of the Cartesian axis
129 : CALL cp_abort(__LOCATION__, &
130 0 : "Dipole correction only for surface perpendicular to one Cartesian axis")
131 : ! not properly general, we assume that vectors A, B, and C are along x y and z respectively,
132 : ! in the ortorhombic cell, but in principle it does not need to be this way, importan
133 : ! is that the cell angles are 90 degrees.
134 : END IF
135 : END IF
136 : END DO
137 :
138 110 : ilow = wf_r%pw_grid%bounds(1, idir_surfdip)
139 110 : iup = wf_r%pw_grid%bounds(2, idir_surfdip)
140 :
141 330 : ALLOCATE (rhoavsurf(ilow:iup))
142 27438 : rhoavsurf = 0.0_dp
143 :
144 110 : bo => wf_r%pw_grid%bounds_local
145 1430 : dh = wf_r%pw_grid%dh
146 :
147 110 : CALL pw_scale(wf_r, wf_r%pw_grid%vol)
148 110 : IF (idir_surfdip == 3) THEN
149 56 : isurf = 1
150 56 : jsurf = 2
151 :
152 13784 : DO i = bo(1, 3), bo(2, 3)
153 13784 : rhoavsurf(i) = accurate_sum(wf_r%array(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), i))
154 : END DO
155 :
156 54 : ELSEIF (idir_surfdip == 2) THEN
157 0 : isurf = 3
158 0 : jsurf = 1
159 :
160 0 : DO i = bo(1, 2), bo(2, 2)
161 0 : rhoavsurf(i) = accurate_sum(wf_r%array(bo(1, 1):bo(2, 1), i, bo(1, 3):bo(2, 3)))
162 : END DO
163 : ELSE
164 54 : isurf = 2
165 54 : jsurf = 3
166 :
167 6854 : DO i = bo(1, 1), bo(2, 1)
168 6854 : rhoavsurf(i) = accurate_sum(wf_r%array(i, bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
169 : END DO
170 : END IF
171 110 : CALL pw_scale(wf_r, 1.0_dp/wf_r%pw_grid%vol)
172 27438 : rhoavsurf = rhoavsurf/wf_r%pw_grid%vol
173 :
174 : surfarea = cell%hmat(isurf, isurf)*cell%hmat(jsurf, jsurf) - &
175 110 : cell%hmat(isurf, jsurf)*cell%hmat(jsurf, isurf)
176 110 : dsurf = surfarea/REAL(ngrid(isurf)*ngrid(jsurf), dp)
177 :
178 110 : IF (wf_r%pw_grid%para%mode /= PW_MODE_LOCAL) THEN
179 110 : CALL wf_r%pw_grid%para%group%sum(rhoavsurf)
180 : END IF
181 27438 : rhoavsurf(ilow:iup) = dsurf*rhoavsurf(ilow:iup)
182 :
183 : ! locate where the vacuum is, and set the reference point for the calculation of the dipole
184 27438 : rhoavsurf(ilow:iup) = rhoavsurf(ilow:iup)/surfarea
185 : ! Note: rhosurf has the same dimension as rho
186 110 : IF (dft_control%pos_dir_surf_dip < 0.0_dp) THEN
187 3200 : ilayer_min = ilow - 1 + MINLOC(ABS(rhoavsurf(ilow:iup)), 1)
188 : ELSE
189 78 : pos_surf_dip = dft_control%pos_dir_surf_dip*bohr
190 78 : ilayer_min = ilow - 1 + NINT(pos_surf_dip/dh(idir_surfdip, idir_surfdip)) + 1
191 : END IF
192 110 : rhoav_min = ABS(rhoavsurf(ilayer_min))
193 110 : IF (rhoav_min >= 1.E-5_dp) THEN
194 0 : CPABORT(" Dipole correction needs more vacuum space above the surface ")
195 : END IF
196 :
197 110 : height_min = REAL((ilayer_min - ilow), dp)*dh(idir_surfdip, idir_surfdip)
198 :
199 : ! surface dipole form average rhoavsurf
200 : ! \sum_i NjdjNkdkdi rhoav_i (i-imin)di
201 110 : dip_hh = 0.0_dp
202 110 : dip_fac = wf_r%pw_grid%vol*dh(idir_surfdip, idir_surfdip)/REAL(ngrid(idir_surfdip), dp)
203 :
204 27438 : DO i = ilayer_min + 1, ilayer_min + ngrid(idir_surfdip)
205 27328 : hh = REAL((i - ilayer_min), dp)
206 27328 : IF (i > iup) THEN
207 19702 : irho = i - ngrid(idir_surfdip)
208 : ELSE
209 : irho = i
210 : END IF
211 : ! introduce a cutoff function to smoothen the edges
212 27328 : IF (ABS(irho - ilayer_min) > width) THEN
213 : cutoff = 1.0_dp
214 : ELSE
215 990 : cutoff = ABS(SIN(0.5_dp*pi*REAL(ABS(irho - ilayer_min), dp)/REAL(width, dp)))
216 : END IF
217 27438 : dip_hh = dip_hh + rhoavsurf(irho)*hh*dip_fac*cutoff
218 : END DO
219 :
220 110 : DEALLOCATE (rhoavsurf)
221 : ! for printing purposes [SGh]
222 110 : qs_env%surface_dipole_moment = dip_hh/bohr
223 :
224 : ! Calculation of the dipole potential as a function of the perpendicular coordinate
225 110 : CALL pw_zero(vdip_r)
226 110 : vdip_fac = dip_hh*4.0_dp*pi
227 :
228 27438 : DO i = ilayer_min + 1, ilayer_min + ngrid(idir_surfdip)
229 27328 : hh = REAL((i - ilayer_min), dp)*dh(idir_surfdip, idir_surfdip)
230 : vdip = vdip_fac*(-0.5_dp + (hh/cell%hmat(idir_surfdip, idir_surfdip)))* &
231 27328 : v_hartree_rspace%pw_grid%dvol/surfarea
232 27328 : IF (i > iup) THEN
233 19702 : irho = i - ngrid(idir_surfdip)
234 : ELSE
235 : irho = i
236 : END IF
237 : ! introduce a cutoff function to smoothen the edges
238 27328 : IF (ABS(irho - ilayer_min) > width) THEN
239 : cutoff = 1.0_dp
240 : ELSE
241 990 : cutoff = ABS(SIN(0.5_dp*pi*REAL(ABS(irho - ilayer_min), dp)/REAL(width, dp)))
242 : END IF
243 27328 : vdip = vdip*cutoff
244 :
245 27438 : IF (idir_surfdip == 3) THEN
246 : vdip_r%array(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), irho) = &
247 17343264 : vdip_r%array(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), irho) + vdip
248 13600 : ELSEIF (idir_surfdip == 2) THEN
249 0 : IF (irho >= bo(1, 2) .AND. irho <= bo(2, 2)) THEN
250 : vdip_r%array(bo(1, 1):bo(2, 1), irho, bo(1, 3):bo(2, 3)) = &
251 0 : vdip_r%array(bo(1, 1):bo(2, 1), irho, bo(1, 3):bo(2, 3)) + vdip
252 : END IF
253 : ELSE
254 13600 : IF (irho >= bo(1, 1) .AND. irho <= bo(2, 1)) THEN
255 : vdip_r%array(irho, bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)) = &
256 26378320 : vdip_r%array(irho, bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)) + vdip
257 : END IF
258 : END IF
259 :
260 : END DO
261 :
262 : ! Dipole correction contribution to the energy
263 110 : energy%surf_dipole = 0.5_dp*pw_integral_ab(vdip_r, wf_r, just_sum=.TRUE.)
264 :
265 : ! Add the dipole potential to the hartree potential on the realspace grid
266 110 : CALL pw_axpy(vdip_r, v_hartree_rspace)
267 :
268 110 : CALL auxbas_pw_pool%give_back_pw(wf_r)
269 110 : CALL auxbas_pw_pool%give_back_pw(vdip_r)
270 :
271 110 : CALL timestop(handle)
272 :
273 110 : END SUBROUTINE calc_dipsurf_potential
274 :
275 : END MODULE surface_dipole
|
