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 : !> \author Ziwei Chai
9 : !> \date 09.04.2026
10 : !> \version 1.0
11 : ! **************************************************************************************************
12 : MODULE gce_methods
13 :
14 : USE cp_control_types, ONLY: dft_control_type,&
15 : pcc_control_type
16 : USE kinds, ONLY: dp
17 : USE message_passing, ONLY: mp_para_env_type
18 : USE pw_methods, ONLY: pw_axpy,&
19 : pw_integrate_function,&
20 : pw_scale,&
21 : pw_transfer,&
22 : pw_zero
23 : USE pw_pool_types, ONLY: pw_pool_type
24 : USE pw_types, ONLY: pw_c1d_gs_type,&
25 : pw_r3d_rs_type
26 : #include "./base/base_uses.f90"
27 :
28 : IMPLICIT NONE
29 :
30 : PRIVATE
31 :
32 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'gce_methods'
33 :
34 : PUBLIC :: planar_averaged_v_hartree_3d, planar_counter_charge
35 :
36 : CONTAINS
37 :
38 : ! **************************************************************************************************
39 : !> \brief map the surface normal to the two in-plane axes and one normal axis
40 : !> \param surf_normal ...
41 : !> \param a ...
42 : !> \param b ...
43 : !> \param c ...
44 : !> \author Ziwei Chai
45 : ! **************************************************************************************************
46 36 : SUBROUTINE get_planar_axes(surf_normal, a, b, c)
47 :
48 : INTEGER, INTENT(IN) :: surf_normal
49 : INTEGER, INTENT(OUT) :: a, b, c
50 :
51 36 : SELECT CASE (surf_normal)
52 : CASE (3)
53 0 : a = 1
54 0 : b = 2
55 0 : c = 3
56 : CASE (1)
57 18 : a = 2
58 18 : b = 3
59 18 : c = 1
60 : CASE (2)
61 18 : a = 3
62 18 : b = 1
63 18 : c = 2
64 : CASE DEFAULT
65 36 : CPABORT("Invalid planar surface normal.")
66 : END SELECT
67 :
68 36 : END SUBROUTINE get_planar_axes
69 :
70 : ! **************************************************************************************************
71 : !> \brief calculate the planar averaged real space potential (e.g. Hartree potential) along the
72 : !> surface normal of a slab model and the corresponding reference (electrostatic) energy
73 : !> near the simulation cell edge.
74 : !> \param v_rspace ...
75 : !> \param dft_control ...
76 : !> \param do_gce ...
77 : !> \param ref_esp ...
78 : !> \param para_env ...
79 : !> \author Ziwei Chai
80 : ! **************************************************************************************************
81 18 : SUBROUTINE planar_averaged_v_hartree_3d(v_rspace, dft_control, do_gce, ref_esp, para_env)
82 :
83 : TYPE(pw_r3d_rs_type), INTENT(IN) :: v_rspace
84 : TYPE(dft_control_type), POINTER :: dft_control
85 : LOGICAL, INTENT(IN) :: do_gce
86 : REAL(KIND=dp), INTENT(INOUT) :: ref_esp
87 : TYPE(mp_para_env_type), POINTER :: para_env
88 :
89 : CHARACTER(LEN=*), PARAMETER :: routineN = 'planar_averaged_v_hartree_3d'
90 :
91 : INTEGER :: a, b, c, handle, x, y, z
92 : REAL(KIND=dp) :: my_ref_esp, ngrids_in_plane
93 18 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: profile, profile_local
94 :
95 18 : CALL timeset(routineN, handle)
96 :
97 18 : IF (do_gce) THEN
98 :
99 : ! a and b cell vectors are parallel to the surface, c vector is not.
100 0 : CALL get_planar_axes(dft_control%pcc_control%surf_normal, a, b, c)
101 :
102 : my_ref_esp = 0.0_dp
103 0 : z = v_rspace%pw_grid%npts(a)*v_rspace%pw_grid%npts(b)
104 0 : ngrids_in_plane = REAL(z, dp)
105 :
106 : !$OMP PARALLEL DO DEFAULT(NONE) &
107 : !$OMP PRIVATE(x,y,z) &
108 : !$OMP SHARED(v_rspace,A,b,c)&
109 0 : !$OMP REDUCTION(+:my_ref_esp)
110 : DO y = v_rspace%pw_grid%bounds_local(1, b), v_rspace%pw_grid%bounds_local(2, b)
111 : DO x = v_rspace%pw_grid%bounds_local(1, a), v_rspace%pw_grid%bounds_local(2, a)
112 : DO z = v_rspace%pw_grid%bounds_local(1, c), v_rspace%pw_grid%bounds_local(2, c)
113 : IF (z == v_rspace%pw_grid%bounds(2, c)) THEN
114 : IF (c == 3) THEN
115 : my_ref_esp = v_rspace%array(x, y, z) + my_ref_esp
116 : ELSE IF (c == 2) THEN
117 : my_ref_esp = v_rspace%array(y, z, x) + my_ref_esp
118 : ELSE
119 : my_ref_esp = v_rspace%array(z, x, y) + my_ref_esp
120 : END IF
121 : END IF
122 : END DO
123 : END DO
124 : END DO
125 : !$OMP END PARALLEL DO
126 :
127 0 : CALL para_env%sum(my_ref_esp)
128 0 : my_ref_esp = my_ref_esp/ngrids_in_plane
129 0 : ref_esp = my_ref_esp
130 :
131 : END IF
132 :
133 18 : IF (dft_control%do_paep) THEN
134 :
135 18 : CALL get_planar_axes(dft_control%paep_control%surf_normal, a, b, c)
136 :
137 18 : z = v_rspace%pw_grid%npts(a)*v_rspace%pw_grid%npts(b)
138 18 : ngrids_in_plane = REAL(z, dp)
139 :
140 54 : ALLOCATE (profile(v_rspace%pw_grid%bounds(1, c):v_rspace%pw_grid%bounds(2, c)))
141 18 : profile(:) = 0.0_dp
142 :
143 : !$OMP PARALLEL DEFAULT(NONE) &
144 : !$OMP PRIVATE(x,y,z,profile_local) &
145 : !$OMP SHARED(v_rspace,A,b,c) &
146 18 : !$OMP SHARED(profile)
147 :
148 : ALLOCATE (profile_local(v_rspace%pw_grid%bounds(1, c):v_rspace%pw_grid%bounds(2, c)))
149 : profile_local(:) = 0.0_dp
150 :
151 : !$OMP DO
152 : DO z = v_rspace%pw_grid%bounds_local(1, c), v_rspace%pw_grid%bounds_local(2, c)
153 : DO y = v_rspace%pw_grid%bounds_local(1, b), v_rspace%pw_grid%bounds_local(2, b)
154 : DO x = v_rspace%pw_grid%bounds_local(1, a), v_rspace%pw_grid%bounds_local(2, a)
155 : IF (c == 3) profile_local(z) = profile_local(z) + v_rspace%array(x, y, z)
156 : IF (c == 2) profile_local(z) = profile_local(z) + v_rspace%array(y, z, x)
157 : IF (c == 1) profile_local(z) = profile_local(z) + v_rspace%array(z, x, y)
158 : END DO
159 : END DO
160 : END DO
161 : !$OMP END DO
162 :
163 : !$OMP CRITICAL
164 : profile(:) = profile(:) + profile_local(:)
165 : !$OMP END CRITICAL
166 :
167 : DEALLOCATE (profile_local)
168 :
169 : !$OMP END PARALLEL
170 :
171 18 : CALL para_env%sum(profile(:))
172 990 : profile(:) = profile(:)/ngrids_in_plane
173 :
174 36 : DEALLOCATE (profile)
175 :
176 : END IF
177 :
178 18 : CALL timestop(handle)
179 :
180 18 : END SUBROUTINE planar_averaged_v_hartree_3d
181 :
182 : ! **************************************************************************************************
183 : !> \brief add the planar counter charge density to the total charge density
184 : !> \param rho_tot_gspace ...
185 : !> \param pcc_env ...
186 : !> \param auxbas_pw_pool ...
187 : !> \author Ziwei Chai
188 : ! **************************************************************************************************
189 36 : SUBROUTINE planar_counter_charge(rho_tot_gspace, pcc_env, auxbas_pw_pool)
190 :
191 : TYPE(pw_c1d_gs_type), INTENT(INOUT) :: rho_tot_gspace
192 : TYPE(pcc_control_type), POINTER :: pcc_env
193 : TYPE(pw_pool_type), INTENT(IN), POINTER :: auxbas_pw_pool
194 :
195 : CHARACTER(LEN=*), PARAMETER :: routineN = 'planar_counter_charge'
196 :
197 : INTEGER :: a, b, c, handle, lcenter, ucenter, x, y, &
198 : z
199 : REAL(KIND=dp) :: dz, gau_a, tot_charge, val
200 : TYPE(pw_c1d_gs_type) :: pcc_density_g
201 : TYPE(pw_r3d_rs_type) :: pcc_density_r
202 :
203 18 : CALL timeset(routineN, handle)
204 :
205 : ! a and b cell vectors are parallel to the surface, c vector is not.
206 18 : CALL get_planar_axes(pcc_env%surf_normal, a, b, c)
207 :
208 18 : CALL auxbas_pw_pool%create_pw(pcc_density_r)
209 18 : CALL auxbas_pw_pool%create_pw(pcc_density_g)
210 18 : CALL pw_zero(pcc_density_r)
211 18 : CALL pw_zero(pcc_density_g)
212 :
213 18 : tot_charge = pw_integrate_function(rho_tot_gspace, isign=-1)
214 :
215 18 : dz = pcc_density_r%pw_grid%dr(c)
216 :
217 18 : ucenter = pcc_density_r%pw_grid%bounds(2, c) - INT(pcc_env%dist_edge/dz)
218 18 : lcenter = pcc_density_r%pw_grid%bounds(1, c) + INT(pcc_env%dist_edge/dz)
219 18 : gau_a = 1.0_dp
220 : !$OMP PARALLEL DO DEFAULT(NONE) &
221 : !$OMP PRIVATE(x,y,z,val) &
222 18 : !$OMP SHARED(pcc_density_r,gau_a,ucenter,lcenter,A,b,c,dz,pcc_env)
223 : DO z = pcc_density_r%pw_grid%bounds_local(1, c), pcc_density_r%pw_grid%bounds_local(2, c)
224 : val = gau_a*EXP(-(dz*REAL(z - ucenter, dp))**2.0_dp/(2.0_dp*pcc_env%gau_c**2.0_dp))
225 : val = val + gau_a*EXP(-(dz*REAL(z - lcenter, dp))**2.0_dp/(2.0_dp*pcc_env%gau_c**2.0_dp))
226 : IF (val < 1.0E-15_dp) val = 0.0_dp
227 : DO y = pcc_density_r%pw_grid%bounds_local(1, b), pcc_density_r%pw_grid%bounds_local(2, b)
228 : DO x = pcc_density_r%pw_grid%bounds_local(1, a), pcc_density_r%pw_grid%bounds_local(2, a)
229 : IF (c == 3) pcc_density_r%array(x, y, z) = val
230 : IF (c == 2) pcc_density_r%array(y, z, x) = val
231 : IF (c == 1) pcc_density_r%array(z, x, y) = val
232 : END DO
233 : END DO
234 : END DO
235 : !$OMP END PARALLEL DO
236 :
237 18 : tot_charge = tot_charge/pw_integrate_function(pcc_density_r, isign=1)
238 18 : CALL pw_scale(pcc_density_r, tot_charge)
239 18 : CALL pw_transfer(pcc_density_r, pcc_density_g)
240 18 : CALL pw_axpy(pcc_density_g, rho_tot_gspace)
241 18 : CALL pw_transfer(rho_tot_gspace, pcc_density_r)
242 :
243 18 : tot_charge = pw_integrate_function(rho_tot_gspace, isign=-1)
244 :
245 18 : CALL auxbas_pw_pool%give_back_pw(pcc_density_r)
246 18 : CALL auxbas_pw_pool%give_back_pw(pcc_density_g)
247 :
248 18 : CALL timestop(handle)
249 :
250 18 : END SUBROUTINE planar_counter_charge
251 :
252 : END MODULE gce_methods
|
