(git:15c1bfc)
Loading...
Searching...
No Matches
surface_dipole.F
Go to the documentation of this file.
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! **************************************************************************************************
10
11 USE cell_types, ONLY: cell_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,&
20 USE pw_methods, ONLY: pw_axpy,&
22 pw_scale,&
25 USE pw_pool_types, ONLY: pw_pool_p_type,&
27 USE pw_types, ONLY: pw_c1d_gs_type,&
32 USE qs_rho_types, ONLY: qs_rho_get,&
35#include "./base/base_uses.f90"
36
37 IMPLICIT NONE
38
39 PRIVATE
40
41 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'surface_dipole'
42
44
45CONTAINS
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 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 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 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 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 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 CALL timeset(routinen, handle)
86 NULLIFY (cell, dft_control, rho, pw_env, auxbas_pw_pool, &
87 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 v_hartree_rspace=v_hartree_rspace)
99
100 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
101 pw_pools=pw_pools)
102 CALL auxbas_pw_pool%create_pw(wf_r)
103 CALL auxbas_pw_pool%create_pw(vdip_r)
104
105 IF (dft_control%qs_control%gapw) THEN
106 IF (dft_control%qs_control%gapw_control%nopaw_as_gpw) THEN
107 CALL pw_axpy(rho_core, rho0_s_gs)
108 IF (ASSOCIATED(rhoz_cneo_s_gs)) THEN
109 CALL pw_axpy(rhoz_cneo_s_gs, rho0_s_gs)
110 END IF
111 CALL pw_transfer(rho0_s_gs, wf_r)
112 CALL pw_axpy(rho_core, rho0_s_gs, -1.0_dp)
113 IF (ASSOCIATED(rhoz_cneo_s_gs)) THEN
114 CALL pw_axpy(rhoz_cneo_s_gs, rho0_s_gs, -1.0_dp)
115 END IF
116 ELSE
117 IF (ASSOCIATED(rhoz_cneo_s_gs)) THEN
118 CALL pw_axpy(rhoz_cneo_s_gs, rho0_s_gs)
119 END IF
120 CALL pw_transfer(rho0_s_gs, wf_r)
121 IF (ASSOCIATED(rhoz_cneo_s_gs)) THEN
122 CALL pw_axpy(rhoz_cneo_s_gs, rho0_s_gs, -1.0_dp)
123 END IF
124 END IF
125 ELSE
126 CALL pw_transfer(rho_core, wf_r)
127 END IF
128 CALL qs_rho_get(rho, rho_r=rho_r)
129 DO ispin = 1, dft_control%nspins
130 CALL pw_axpy(rho_r(ispin), wf_r)
131 END DO
132
133 ngrid(1:3) = wf_r%pw_grid%npts(1:3)
134 idir_surfdip = dft_control%dir_surf_dip
135
136 width = 4
137
138 DO i = 1, 3
139 IF (i /= idir_surfdip) THEN
140 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 "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 ilow = wf_r%pw_grid%bounds(1, idir_surfdip)
152 iup = wf_r%pw_grid%bounds(2, idir_surfdip)
153
154 ALLOCATE (rhoavsurf(ilow:iup))
155 rhoavsurf = 0.0_dp
156
157 bo => wf_r%pw_grid%bounds_local
158 dh = wf_r%pw_grid%dh
159
160 CALL pw_scale(wf_r, wf_r%pw_grid%vol)
161 IF (idir_surfdip == 3) THEN
162 isurf = 1
163 jsurf = 2
164
165 DO i = bo(1, 3), bo(2, 3)
166 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 ELSEIF (idir_surfdip == 2) THEN
170 isurf = 3
171 jsurf = 1
172
173 DO i = bo(1, 2), bo(2, 2)
174 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 isurf = 2
178 jsurf = 3
179
180 DO i = bo(1, 1), bo(2, 1)
181 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 CALL pw_scale(wf_r, 1.0_dp/wf_r%pw_grid%vol)
185 rhoavsurf = rhoavsurf/wf_r%pw_grid%vol
186
187 surfarea = cell%hmat(isurf, isurf)*cell%hmat(jsurf, jsurf) - &
188 cell%hmat(isurf, jsurf)*cell%hmat(jsurf, isurf)
189 dsurf = surfarea/real(ngrid(isurf)*ngrid(jsurf), dp)
190
191 IF (wf_r%pw_grid%para%mode /= pw_mode_local) THEN
192 CALL wf_r%pw_grid%para%group%sum(rhoavsurf)
193 END IF
194 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 rhoavsurf(ilow:iup) = rhoavsurf(ilow:iup)/surfarea
198 ! Note: rhosurf has the same dimension as rho
199 IF (dft_control%pos_dir_surf_dip < 0.0_dp) THEN
200 ilayer_min = ilow - 1 + minloc(abs(rhoavsurf(ilow:iup)), 1)
201 ELSE
202 pos_surf_dip = dft_control%pos_dir_surf_dip*bohr
203 ilayer_min = ilow - 1 + nint(pos_surf_dip/dh(idir_surfdip, idir_surfdip)) + 1
204 END IF
205 rhoav_min = abs(rhoavsurf(ilayer_min))
206 IF (rhoav_min >= 1.e-5_dp) THEN
207 cpabort(" Dipole correction needs more vacuum space above the surface ")
208 END IF
209
210 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 dip_hh = 0.0_dp
215 dip_fac = wf_r%pw_grid%vol*dh(idir_surfdip, idir_surfdip)/real(ngrid(idir_surfdip), dp)
216
217 DO i = ilayer_min + 1, ilayer_min + ngrid(idir_surfdip)
218 hh = real((i - ilayer_min), dp)
219 IF (i > iup) THEN
220 irho = i - ngrid(idir_surfdip)
221 ELSE
222 irho = i
223 END IF
224! introduce a cutoff function to smoothen the edges
225 IF (abs(irho - ilayer_min) > width) THEN
226 cutoff = 1.0_dp
227 ELSE
228 cutoff = abs(sin(0.5_dp*pi*real(abs(irho - ilayer_min), dp)/real(width, dp)))
229 END IF
230 dip_hh = dip_hh + rhoavsurf(irho)*hh*dip_fac*cutoff
231 END DO
232
233 DEALLOCATE (rhoavsurf)
234! for printing purposes [SGh]
235 qs_env%surface_dipole_moment = dip_hh/bohr
236
237! Calculation of the dipole potential as a function of the perpendicular coordinate
238 CALL pw_zero(vdip_r)
239 vdip_fac = dip_hh*4.0_dp*pi
240
241 DO i = ilayer_min + 1, ilayer_min + ngrid(idir_surfdip)
242 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 v_hartree_rspace%pw_grid%dvol/surfarea
245 IF (i > iup) THEN
246 irho = i - ngrid(idir_surfdip)
247 ELSE
248 irho = i
249 END IF
250! introduce a cutoff function to smoothen the edges
251 IF (abs(irho - ilayer_min) > width) THEN
252 cutoff = 1.0_dp
253 ELSE
254 cutoff = abs(sin(0.5_dp*pi*real(abs(irho - ilayer_min), dp)/real(width, dp)))
255 END IF
256 vdip = vdip*cutoff
257
258 IF (idir_surfdip == 3) THEN
259 vdip_r%array(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), irho) = &
260 vdip_r%array(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), irho) + vdip
261 ELSEIF (idir_surfdip == 2) THEN
262 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 vdip_r%array(bo(1, 1):bo(2, 1), irho, bo(1, 3):bo(2, 3)) + vdip
265 END IF
266 ELSE
267 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 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 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 CALL pw_axpy(vdip_r, v_hartree_rspace)
280
281 CALL auxbas_pw_pool%give_back_pw(wf_r)
282 CALL auxbas_pw_pool%give_back_pw(vdip_r)
283
284 CALL timestop(handle)
285
286 END SUBROUTINE calc_dipsurf_potential
287
288END MODULE surface_dipole
Handles all functions related to the CELL.
Definition cell_types.F:15
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
sums arrays of real/complex numbers with much reduced round-off as compared to a naive implementation...
Definition kahan_sum.F:29
Defines the basic variable types.
Definition kinds.F:23
integer, parameter, public dp
Definition kinds.F:34
Definition of mathematical constants and functions.
real(kind=dp), parameter, public pi
Definition of physical constants:
Definition physcon.F:68
real(kind=dp), parameter, public bohr
Definition physcon.F:147
container for various plainwaves related things
subroutine, public pw_env_get(pw_env, pw_pools, cube_info, gridlevel_info, auxbas_pw_pool, auxbas_grid, auxbas_rs_desc, auxbas_rs_grid, rs_descs, rs_grids, xc_pw_pool, vdw_pw_pool, poisson_env, interp_section)
returns the various attributes of the pw env
integer, parameter, public pw_mode_local
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
superstucture that hold various representations of the density and keeps track of which ones are vali...
subroutine, public qs_rho_get(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
returns info about the density described by this object. If some representation is not available an e...
types that represent a quickstep subsys
subroutine, public calc_dipsurf_potential(qs_env, energy)
compute the surface dipole and the correction to the hartree potential
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
contained for different pw related things
to create arrays of pools
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
keeps the density in various representations, keeping track of which ones are valid.