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 : !> \brief contains information regarding the decoupling/recoupling method of Bloechl
10 : !> \author Teodoro Laino
11 : ! **************************************************************************************************
12 : MODULE cp_ddapc_methods
13 : USE cell_types, ONLY: cell_type
14 : USE cp_log_handling, ONLY: cp_logger_get_default_io_unit
15 : USE input_constants, ONLY: weight_type_mass,&
16 : weight_type_unit
17 : USE input_section_types, ONLY: section_vals_type,&
18 : section_vals_val_get,&
19 : section_vals_val_set
20 : USE kahan_sum, ONLY: accurate_sum
21 : USE kinds, ONLY: dp
22 : USE mathconstants, ONLY: fourpi,&
23 : oorootpi,&
24 : pi,&
25 : twopi
26 : USE mathlib, ONLY: diamat_all,&
27 : invert_matrix
28 : USE message_passing, ONLY: mp_para_env_type
29 : USE particle_types, ONLY: particle_type
30 : USE pw_spline_utils, ONLY: Eval_Interp_Spl3_pbc
31 : USE pw_types, ONLY: pw_c1d_gs_type,&
32 : pw_r3d_rs_type
33 : USE spherical_harmonics, ONLY: legendre
34 : #include "./base/base_uses.f90"
35 :
36 : IMPLICIT NONE
37 : PRIVATE
38 : LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .FALSE.
39 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cp_ddapc_methods'
40 : PUBLIC :: ddapc_eval_gfunc, &
41 : build_b_vector, &
42 : build_der_b_vector, &
43 : build_A_matrix, &
44 : build_der_A_matrix_rows, &
45 : prep_g_dot_rvec_sin_cos, &
46 : cleanup_g_dot_rvec_sin_cos, &
47 : ddapc_eval_AmI, &
48 : ewald_ddapc_pot, &
49 : solvation_ddapc_pot
50 :
51 : CONTAINS
52 :
53 : ! **************************************************************************************************
54 : !> \brief ...
55 : !> \param gfunc ...
56 : !> \param w ...
57 : !> \param gcut ...
58 : !> \param rho_tot_g ...
59 : !> \param radii ...
60 : ! **************************************************************************************************
61 296 : SUBROUTINE ddapc_eval_gfunc(gfunc, w, gcut, rho_tot_g, radii)
62 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: gfunc
63 : REAL(kind=dp), DIMENSION(:), POINTER :: w
64 : REAL(KIND=dp), INTENT(IN) :: gcut
65 : TYPE(pw_c1d_gs_type), INTENT(IN) :: rho_tot_g
66 : REAL(kind=dp), DIMENSION(:), POINTER :: radii
67 :
68 : CHARACTER(len=*), PARAMETER :: routineN = 'ddapc_eval_gfunc'
69 :
70 : INTEGER :: e_dim, handle, ig, igauss, s_dim
71 : REAL(KIND=dp) :: g2, gcut2, rc, rc2
72 :
73 296 : CALL timeset(routineN, handle)
74 296 : gcut2 = gcut*gcut
75 : !
76 296 : s_dim = rho_tot_g%pw_grid%first_gne0
77 296 : e_dim = rho_tot_g%pw_grid%ngpts_cut_local
78 1184 : ALLOCATE (gfunc(s_dim:e_dim, SIZE(radii)))
79 888 : ALLOCATE (w(s_dim:e_dim))
80 73665770 : gfunc = 0.0_dp
81 24772666 : w = 0.0_dp
82 1136 : DO igauss = 1, SIZE(radii)
83 840 : rc = radii(igauss)
84 840 : rc2 = rc*rc
85 556592 : DO ig = s_dim, e_dim
86 556296 : g2 = rho_tot_g%pw_grid%gsq(ig)
87 556296 : IF (g2 > gcut2) EXIT
88 556296 : gfunc(ig, igauss) = EXP(-g2*rc2/4.0_dp)
89 : END DO
90 : END DO
91 186764 : DO ig = s_dim, e_dim
92 186764 : g2 = rho_tot_g%pw_grid%gsq(ig)
93 186764 : IF (g2 > gcut2) EXIT
94 186764 : w(ig) = fourpi*(g2 - gcut2)**2/(g2*gcut2)
95 : END DO
96 296 : CALL timestop(handle)
97 296 : END SUBROUTINE ddapc_eval_gfunc
98 :
99 : ! **************************************************************************************************
100 : !> \brief Computes the B vector for the solution of the linear system
101 : !> \param bv ...
102 : !> \param gfunc ...
103 : !> \param w ...
104 : !> \param particle_set ...
105 : !> \param radii ...
106 : !> \param rho_tot_g ...
107 : !> \param gcut ...
108 : !> \par History
109 : !> 08.2005 created [tlaino]
110 : !> \author Teodoro Laino
111 : ! **************************************************************************************************
112 2588 : SUBROUTINE build_b_vector(bv, gfunc, w, particle_set, radii, rho_tot_g, gcut)
113 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: bv
114 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: gfunc
115 : REAL(KIND=dp), DIMENSION(:), POINTER :: w
116 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
117 : REAL(KIND=dp), DIMENSION(:), POINTER :: radii
118 : TYPE(pw_c1d_gs_type), INTENT(IN) :: rho_tot_g
119 : REAL(KIND=dp), INTENT(IN) :: gcut
120 :
121 : CHARACTER(len=*), PARAMETER :: routineN = 'build_b_vector'
122 :
123 : COMPLEX(KIND=dp) :: phase
124 : INTEGER :: e_dim, handle, idim, ig, igauss, igmax, &
125 : iparticle, s_dim
126 : REAL(KIND=dp) :: arg, g2, gcut2, gvec(3), rvec(3)
127 2588 : REAL(KIND=dp), DIMENSION(:), POINTER :: my_bv, my_bvw
128 :
129 2588 : CALL timeset(routineN, handle)
130 2588 : NULLIFY (my_bv, my_bvw)
131 2588 : gcut2 = gcut*gcut
132 2588 : s_dim = rho_tot_g%pw_grid%first_gne0
133 2588 : e_dim = rho_tot_g%pw_grid%ngpts_cut_local
134 2588 : igmax = 0
135 1037996 : DO ig = s_dim, e_dim
136 1037996 : g2 = rho_tot_g%pw_grid%gsq(ig)
137 1037996 : IF (g2 > gcut2) EXIT
138 1037996 : igmax = ig
139 : END DO
140 2588 : IF (igmax >= s_dim) THEN
141 7764 : ALLOCATE (my_bv(s_dim:igmax))
142 5176 : ALLOCATE (my_bvw(s_dim:igmax))
143 : !
144 9884 : DO iparticle = 1, SIZE(particle_set)
145 29184 : rvec = particle_set(iparticle)%r
146 3512084 : my_bv = 0.0_dp
147 3512084 : DO ig = s_dim, igmax
148 14019152 : gvec = rho_tot_g%pw_grid%g(:, ig)
149 14019152 : arg = DOT_PRODUCT(gvec, rvec)
150 3504788 : phase = CMPLX(COS(arg), -SIN(arg), KIND=dp)
151 3512084 : my_bv(ig) = w(ig)*REAL(CONJG(rho_tot_g%array(ig))*phase, KIND=dp)
152 : END DO
153 29396 : DO igauss = 1, SIZE(radii)
154 19512 : idim = (iparticle - 1)*SIZE(radii) + igauss
155 10333500 : DO ig = s_dim, igmax
156 10333500 : my_bvw(ig) = my_bv(ig)*gfunc(ig, igauss)
157 : END DO
158 26808 : bv(idim) = accurate_sum(my_bvw)
159 : END DO
160 : END DO
161 2588 : DEALLOCATE (my_bvw)
162 2588 : DEALLOCATE (my_bv)
163 : ELSE
164 0 : DO iparticle = 1, SIZE(particle_set)
165 0 : DO igauss = 1, SIZE(radii)
166 0 : idim = (iparticle - 1)*SIZE(radii) + igauss
167 0 : bv(idim) = 0.0_dp
168 : END DO
169 : END DO
170 : END IF
171 2588 : CALL timestop(handle)
172 2588 : END SUBROUTINE build_b_vector
173 :
174 : ! **************************************************************************************************
175 : !> \brief Computes the A matrix for the solution of the linear system
176 : !> \param Am ...
177 : !> \param gfunc ...
178 : !> \param w ...
179 : !> \param particle_set ...
180 : !> \param radii ...
181 : !> \param rho_tot_g ...
182 : !> \param gcut ...
183 : !> \param g_dot_rvec_sin ...
184 : !> \param g_dot_rvec_cos ...
185 : !> \par History
186 : !> 08.2005 created [tlaino]
187 : !> \author Teodoro Laino
188 : !> \note NB accept g_dot_rvec_* arrays
189 : ! **************************************************************************************************
190 296 : SUBROUTINE build_A_matrix(Am, gfunc, w, particle_set, radii, rho_tot_g, gcut, g_dot_rvec_sin, g_dot_rvec_cos)
191 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: Am
192 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: gfunc
193 : REAL(KIND=dp), DIMENSION(:), POINTER :: w
194 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
195 : REAL(KIND=dp), DIMENSION(:), POINTER :: radii
196 : TYPE(pw_c1d_gs_type), INTENT(IN) :: rho_tot_g
197 : REAL(KIND=dp), INTENT(IN) :: gcut
198 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: g_dot_rvec_sin, g_dot_rvec_cos
199 :
200 : CHARACTER(len=*), PARAMETER :: routineN = 'build_A_matrix'
201 :
202 : INTEGER :: e_dim, handle, idim1, idim2, ig, &
203 : igauss1, igauss2, igmax, iparticle1, &
204 : iparticle2, istart_g, s_dim
205 : REAL(KIND=dp) :: g2, gcut2, tmp
206 296 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: my_Am, my_Amw
207 296 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: gfunc_sq
208 :
209 : !NB precalculate as many things outside of the innermost loop as possible, in particular w(ig)*gfunc(ig,igauus1)*gfunc(ig,igauss2)
210 :
211 296 : CALL timeset(routineN, handle)
212 296 : gcut2 = gcut*gcut
213 296 : s_dim = rho_tot_g%pw_grid%first_gne0
214 296 : e_dim = rho_tot_g%pw_grid%ngpts_cut_local
215 296 : igmax = 0
216 186764 : DO ig = s_dim, e_dim
217 186764 : g2 = rho_tot_g%pw_grid%gsq(ig)
218 186764 : IF (g2 > gcut2) EXIT
219 186764 : igmax = ig
220 : END DO
221 296 : IF (igmax >= s_dim) THEN
222 888 : ALLOCATE (my_Am(s_dim:igmax))
223 592 : ALLOCATE (my_Amw(s_dim:igmax))
224 1480 : ALLOCATE (gfunc_sq(s_dim:igmax, SIZE(radii), SIZE(radii)))
225 :
226 1136 : DO igauss1 = 1, SIZE(radii)
227 3608 : DO igauss2 = 1, SIZE(radii)
228 1665732 : gfunc_sq(s_dim:igmax, igauss1, igauss2) = w(s_dim:igmax)*gfunc(s_dim:igmax, igauss1)*gfunc(s_dim:igmax, igauss2)
229 : END DO
230 : END DO
231 :
232 1178 : DO iparticle1 = 1, SIZE(particle_set)
233 4912 : DO iparticle2 = iparticle1, SIZE(particle_set)
234 6964954 : DO ig = s_dim, igmax
235 : !NB replace explicit dot product and cosine with cos(A+B) formula - much faster
236 : my_Am(ig) = (g_dot_rvec_cos(ig - s_dim + 1, iparticle1)*g_dot_rvec_cos(ig - s_dim + 1, iparticle2) + &
237 6964954 : g_dot_rvec_sin(ig - s_dim + 1, iparticle1)*g_dot_rvec_sin(ig - s_dim + 1, iparticle2))
238 : END DO
239 15638 : DO igauss1 = 1, SIZE(radii)
240 11022 : idim1 = (iparticle1 - 1)*SIZE(radii) + igauss1
241 11022 : istart_g = 1
242 11022 : IF (iparticle2 == iparticle1) istart_g = igauss1
243 45158 : DO igauss2 = istart_g, SIZE(radii)
244 30402 : idim2 = (iparticle2 - 1)*SIZE(radii) + igauss2
245 60402252 : my_Amw(s_dim:igmax) = my_Am(s_dim:igmax)*gfunc_sq(s_dim:igmax, igauss1, igauss2)
246 : !NB no loss of accuracy in my test cases
247 : !tmp = accurate_sum(my_Amw)
248 60402252 : tmp = SUM(my_Amw)
249 30402 : Am(idim2, idim1) = tmp
250 41424 : Am(idim1, idim2) = tmp
251 : END DO
252 : END DO
253 : END DO
254 : END DO
255 296 : DEALLOCATE (gfunc_sq)
256 296 : DEALLOCATE (my_Amw)
257 296 : DEALLOCATE (my_Am)
258 : END IF
259 296 : CALL timestop(handle)
260 296 : END SUBROUTINE build_A_matrix
261 :
262 : ! **************************************************************************************************
263 : !> \brief Computes the derivative of B vector for the evaluation of the Pulay forces
264 : !> \param dbv ...
265 : !> \param gfunc ...
266 : !> \param w ...
267 : !> \param particle_set ...
268 : !> \param radii ...
269 : !> \param rho_tot_g ...
270 : !> \param gcut ...
271 : !> \param iparticle0 ...
272 : !> \par History
273 : !> 08.2005 created [tlaino]
274 : !> \author Teodoro Laino
275 : ! **************************************************************************************************
276 420 : SUBROUTINE build_der_b_vector(dbv, gfunc, w, particle_set, radii, rho_tot_g, gcut, iparticle0)
277 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: dbv
278 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: gfunc
279 : REAL(KIND=dp), DIMENSION(:), POINTER :: w
280 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
281 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: radii
282 : TYPE(pw_c1d_gs_type), INTENT(IN) :: rho_tot_g
283 : REAL(KIND=dp), INTENT(IN) :: gcut
284 : INTEGER, INTENT(IN) :: iparticle0
285 :
286 : CHARACTER(len=*), PARAMETER :: routineN = 'build_der_b_vector'
287 :
288 : COMPLEX(KIND=dp) :: dphase
289 : INTEGER :: e_dim, handle, idim, ig, igauss, igmax, &
290 : iparticle, s_dim
291 : REAL(KIND=dp) :: arg, g2, gcut2
292 420 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: my_dbvw
293 420 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: my_dbv
294 : REAL(KIND=dp), DIMENSION(3) :: gvec, rvec
295 :
296 420 : CALL timeset(routineN, handle)
297 420 : gcut2 = gcut*gcut
298 420 : s_dim = rho_tot_g%pw_grid%first_gne0
299 420 : e_dim = rho_tot_g%pw_grid%ngpts_cut_local
300 420 : igmax = 0
301 360282 : DO ig = s_dim, e_dim
302 360282 : g2 = rho_tot_g%pw_grid%gsq(ig)
303 360282 : IF (g2 > gcut2) EXIT
304 360282 : igmax = ig
305 : END DO
306 420 : IF (igmax >= s_dim) THEN
307 1260 : ALLOCATE (my_dbv(3, s_dim:igmax))
308 1260 : ALLOCATE (my_dbvw(s_dim:igmax))
309 1908 : DO iparticle = 1, SIZE(particle_set)
310 1488 : IF (iparticle /= iparticle0) CYCLE
311 1680 : rvec = particle_set(iparticle)%r
312 360282 : DO ig = s_dim, igmax
313 1439448 : gvec = rho_tot_g%pw_grid%g(:, ig)
314 1439448 : arg = DOT_PRODUCT(gvec, rvec)
315 359862 : dphase = -CMPLX(SIN(arg), COS(arg), KIND=dp)
316 1439868 : my_dbv(:, ig) = w(ig)*REAL(CONJG(rho_tot_g%array(ig))*dphase, KIND=dp)*gvec(:)
317 : END DO
318 3060 : DO igauss = 1, SIZE(radii)
319 1152 : idim = (iparticle - 1)*SIZE(radii) + igauss
320 1071630 : DO ig = s_dim, igmax
321 1071630 : my_dbvw(ig) = my_dbv(1, ig)*gfunc(ig, igauss)
322 : END DO
323 1152 : dbv(idim, 1) = accurate_sum(my_dbvw)
324 1071630 : DO ig = s_dim, igmax
325 1071630 : my_dbvw(ig) = my_dbv(2, ig)*gfunc(ig, igauss)
326 : END DO
327 1152 : dbv(idim, 2) = accurate_sum(my_dbvw)
328 1071630 : DO ig = s_dim, igmax
329 1071630 : my_dbvw(ig) = my_dbv(3, ig)*gfunc(ig, igauss)
330 : END DO
331 1572 : dbv(idim, 3) = accurate_sum(my_dbvw)
332 : END DO
333 : END DO
334 420 : DEALLOCATE (my_dbvw)
335 420 : DEALLOCATE (my_dbv)
336 : ELSE
337 0 : DO iparticle = 1, SIZE(particle_set)
338 0 : IF (iparticle /= iparticle0) CYCLE
339 0 : DO igauss = 1, SIZE(radii)
340 0 : idim = (iparticle - 1)*SIZE(radii) + igauss
341 0 : dbv(idim, 1:3) = 0.0_dp
342 : END DO
343 : END DO
344 : END IF
345 420 : CALL timestop(handle)
346 420 : END SUBROUTINE build_der_b_vector
347 :
348 : ! **************************************************************************************************
349 : !> \brief Computes the derivative of the A matrix for the evaluation of the
350 : !> Pulay forces
351 : !> \param dAm ...
352 : !> \param gfunc ...
353 : !> \param w ...
354 : !> \param particle_set ...
355 : !> \param radii ...
356 : !> \param rho_tot_g ...
357 : !> \param gcut ...
358 : !> \param iparticle0 ...
359 : !> \param nparticles ...
360 : !> \param g_dot_rvec_sin ...
361 : !> \param g_dot_rvec_cos ...
362 : !> \par History
363 : !> 08.2005 created [tlaino]
364 : !> \author Teodoro Laino
365 : !> \note NB accept g_dot_rvec_* arrays
366 : ! **************************************************************************************************
367 444 : SUBROUTINE build_der_A_matrix_rows(dAm, gfunc, w, particle_set, radii, &
368 148 : rho_tot_g, gcut, iparticle0, nparticles, g_dot_rvec_sin, g_dot_rvec_cos)
369 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: dAm
370 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: gfunc
371 : REAL(KIND=dp), DIMENSION(:), POINTER :: w
372 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
373 : REAL(KIND=dp), DIMENSION(:), POINTER :: radii
374 : TYPE(pw_c1d_gs_type), INTENT(IN) :: rho_tot_g
375 : REAL(KIND=dp), INTENT(IN) :: gcut
376 : INTEGER, INTENT(IN) :: iparticle0, nparticles
377 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: g_dot_rvec_sin, g_dot_rvec_cos
378 :
379 : CHARACTER(len=*), PARAMETER :: routineN = 'build_der_A_matrix_rows'
380 :
381 : INTEGER :: e_dim, handle, ig, igauss2, igmax, &
382 : iparticle1, iparticle2, s_dim
383 : REAL(KIND=dp) :: g2, gcut2
384 :
385 : !NB calculate derivatives for a block of particles, just the row parts (since derivative matrix is symmetric)
386 : !NB Use DGEMM to speed up calculation, can't do accurate_sum() anymore because dgemm does the sum over g
387 :
388 : EXTERNAL DGEMM
389 148 : REAL(KIND=dp), ALLOCATABLE :: lhs(:, :), rhs(:, :)
390 : INTEGER :: Nr, Np, Ng, icomp, ipp
391 :
392 148 : CALL timeset(routineN, handle)
393 148 : gcut2 = gcut*gcut
394 148 : s_dim = rho_tot_g%pw_grid%first_gne0
395 148 : e_dim = rho_tot_g%pw_grid%ngpts_cut_local
396 148 : igmax = 0
397 152182 : DO ig = s_dim, e_dim
398 152182 : g2 = rho_tot_g%pw_grid%gsq(ig)
399 152182 : IF (g2 > gcut2) EXIT
400 152182 : igmax = ig
401 : END DO
402 :
403 148 : Nr = SIZE(radii)
404 148 : Np = SIZE(particle_set)
405 148 : Ng = igmax - s_dim + 1
406 148 : IF (igmax >= s_dim) THEN
407 592 : ALLOCATE (lhs(nparticles*Nr, Ng))
408 592 : ALLOCATE (rhs(Ng, Np*Nr))
409 :
410 : ! rhs with first term of sin(g.(rvec1-rvec2))
411 : ! rhs has all parts that depend on iparticle2
412 568 : DO iparticle2 = 1, Np
413 1720 : DO igauss2 = 1, Nr
414 1072050 : rhs(1:Ng, (iparticle2 - 1)*Nr + igauss2) = g_dot_rvec_sin(1:Ng, iparticle2)*gfunc(s_dim:igmax, igauss2)
415 : END DO
416 : END DO
417 592 : DO icomp = 1, 3
418 : ! create lhs, which has all parts that depend on iparticle1
419 1704 : DO ipp = 1, nparticles
420 1260 : iparticle1 = iparticle0 + ipp - 1
421 1081290 : DO ig = s_dim, igmax
422 : lhs((ipp - 1)*Nr + 1:(ipp - 1)*Nr + Nr, ig - s_dim + 1) = w(ig)*rho_tot_g%pw_grid%g(icomp, ig)* &
423 4292280 : gfunc(ig, 1:Nr)*g_dot_rvec_cos(ig - s_dim + 1, iparticle1)
424 : END DO
425 : END DO ! ipp
426 : ! do main multiply
427 : CALL DGEMM('N', 'N', nparticles*Nr, Np*Nr, Ng, 1.0D0, lhs(1, 1), nparticles*Nr, rhs(1, 1), &
428 444 : Ng, 0.0D0, dAm((iparticle0 - 1)*Nr + 1, 1, icomp), Np*Nr)
429 : ! do extra multiplies to compensate for missing factor of 2
430 1852 : DO ipp = 1, nparticles
431 1260 : iparticle1 = iparticle0 + ipp - 1
432 : CALL DGEMM('N', 'N', Nr, Nr, Ng, 1.0D0, lhs((ipp - 1)*Nr + 1, 1), nparticles*Nr, rhs(1, (iparticle1 - 1)*Nr + 1), &
433 1704 : Ng, 1.0D0, dAm((iparticle1 - 1)*Nr + 1, (iparticle1 - 1)*Nr + 1, icomp), Np*Nr)
434 : END DO
435 : ! now extra columns to account for factor of 2 in some rhs columns
436 : END DO ! icomp
437 :
438 : ! rhs with second term of sin(g.(rvec1-rvec2))
439 : ! rhs has all parts that depend on iparticle2
440 568 : DO iparticle2 = 1, Np
441 1720 : DO igauss2 = 1, Nr
442 1072050 : rhs(1:Ng, (iparticle2 - 1)*Nr + igauss2) = -g_dot_rvec_cos(1:Ng, iparticle2)*gfunc(s_dim:igmax, igauss2)
443 : END DO
444 : END DO
445 592 : DO icomp = 1, 3
446 : ! create lhs, which has all parts that depend on iparticle1
447 1704 : DO ipp = 1, nparticles
448 1260 : iparticle1 = iparticle0 + ipp - 1
449 1081290 : DO ig = s_dim, igmax
450 : lhs((ipp - 1)*Nr + 1:(ipp - 1)*Nr + Nr, ig - s_dim + 1) = w(ig)*rho_tot_g%pw_grid%g(icomp, ig)*gfunc(ig, 1:Nr)* &
451 4292280 : g_dot_rvec_sin(ig - s_dim + 1, iparticle1)
452 : END DO
453 : END DO
454 : ! do main multiply
455 : CALL DGEMM('N', 'N', nparticles*Nr, Np*Nr, Ng, 1.0D0, lhs(1, 1), nparticles*Nr, rhs(1, 1), &
456 444 : Ng, 1.0D0, dAm((iparticle0 - 1)*Nr + 1, 1, icomp), Np*Nr)
457 : ! do extra multiples to compensate for missing factor of 2
458 1852 : DO ipp = 1, nparticles
459 1260 : iparticle1 = iparticle0 + ipp - 1
460 : CALL DGEMM('N', 'N', Nr, Nr, Ng, 1.0D0, lhs((ipp - 1)*Nr + 1, 1), nparticles*Nr, rhs(1, (iparticle1 - 1)*Nr + 1), &
461 1704 : Ng, 1.0D0, dAm((iparticle1 - 1)*Nr + 1, (iparticle1 - 1)*Nr + 1, icomp), Np*Nr)
462 : END DO
463 : END DO
464 :
465 148 : DEALLOCATE (rhs)
466 148 : DEALLOCATE (lhs)
467 : ELSE
468 : ! Some ranks may not own G-vectors below the cutoff.
469 0 : dAm((iparticle0 - 1)*Nr + 1:(iparticle0 + nparticles - 1)*Nr, :, :) = 0.0_dp
470 : END IF
471 148 : CALL timestop(handle)
472 148 : END SUBROUTINE build_der_A_matrix_rows
473 :
474 : ! **************************************************************************************************
475 : !> \brief deallocate g_dot_rvec_* arrays
476 : !> \param g_dot_rvec_sin ...
477 : !> \param g_dot_rvec_cos ...
478 : ! **************************************************************************************************
479 444 : SUBROUTINE cleanup_g_dot_rvec_sin_cos(g_dot_rvec_sin, g_dot_rvec_cos)
480 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: g_dot_rvec_sin, g_dot_rvec_cos
481 :
482 444 : IF (ALLOCATED(g_dot_rvec_sin)) DEALLOCATE (g_dot_rvec_sin)
483 444 : IF (ALLOCATED(g_dot_rvec_cos)) DEALLOCATE (g_dot_rvec_cos)
484 444 : END SUBROUTINE cleanup_g_dot_rvec_sin_cos
485 :
486 : ! **************************************************************************************************
487 : !> \brief precompute sin(g.r) and cos(g.r) for quicker evaluations of sin(g.(r1-r2)) and cos(g.(r1-r2))
488 : !> \param rho_tot_g ...
489 : !> \param particle_set ...
490 : !> \param gcut ...
491 : !> \param g_dot_rvec_sin ...
492 : !> \param g_dot_rvec_cos ...
493 : ! **************************************************************************************************
494 444 : SUBROUTINE prep_g_dot_rvec_sin_cos(rho_tot_g, particle_set, gcut, g_dot_rvec_sin, g_dot_rvec_cos)
495 : TYPE(pw_c1d_gs_type), INTENT(IN) :: rho_tot_g
496 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
497 : REAL(KIND=dp), INTENT(IN) :: gcut
498 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: g_dot_rvec_sin, g_dot_rvec_cos
499 :
500 : INTEGER :: e_dim, ig, igmax, iparticle, s_dim
501 : REAL(KIND=dp) :: g2, g_dot_rvec, gcut2, rvec(3)
502 :
503 444 : gcut2 = gcut*gcut
504 444 : s_dim = rho_tot_g%pw_grid%first_gne0
505 444 : e_dim = rho_tot_g%pw_grid%ngpts_cut_local
506 444 : igmax = 0
507 338946 : DO ig = s_dim, e_dim
508 338946 : g2 = rho_tot_g%pw_grid%gsq(ig)
509 338946 : IF (g2 > gcut2) EXIT
510 338946 : igmax = ig
511 : END DO
512 :
513 444 : IF (igmax >= s_dim) THEN
514 1776 : ALLOCATE (g_dot_rvec_sin(1:igmax - s_dim + 1, SIZE(particle_set)))
515 1332 : ALLOCATE (g_dot_rvec_cos(1:igmax - s_dim + 1, SIZE(particle_set)))
516 :
517 1746 : DO iparticle = 1, SIZE(particle_set)
518 5208 : rvec = particle_set(iparticle)%r
519 1105232 : DO ig = s_dim, igmax
520 4413944 : g_dot_rvec = DOT_PRODUCT(rho_tot_g%pw_grid%g(:, ig), rvec)
521 1103486 : g_dot_rvec_sin(ig - s_dim + 1, iparticle) = SIN(g_dot_rvec)
522 1104788 : g_dot_rvec_cos(ig - s_dim + 1, iparticle) = COS(g_dot_rvec)
523 : END DO
524 : END DO
525 : ELSE
526 : ! Pass valid zero-sized arrays to the downstream routines.
527 0 : ALLOCATE (g_dot_rvec_sin(0, SIZE(particle_set)))
528 0 : ALLOCATE (g_dot_rvec_cos(0, SIZE(particle_set)))
529 : END IF
530 :
531 444 : END SUBROUTINE prep_g_dot_rvec_sin_cos
532 :
533 : ! **************************************************************************************************
534 : !> \brief Computes the inverse AmI of the Am matrix
535 : !> \param GAmI ...
536 : !> \param c0 ...
537 : !> \param gfunc ...
538 : !> \param w ...
539 : !> \param particle_set ...
540 : !> \param gcut ...
541 : !> \param rho_tot_g ...
542 : !> \param radii ...
543 : !> \param iw ...
544 : !> \param Vol ...
545 : !> \par History
546 : !> 12.2005 created [tlaino]
547 : !> \author Teodoro Laino
548 : ! **************************************************************************************************
549 296 : SUBROUTINE ddapc_eval_AmI(GAmI, c0, gfunc, w, particle_set, gcut, &
550 : rho_tot_g, radii, iw, Vol)
551 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: GAmI
552 : REAL(KIND=dp), INTENT(OUT) :: c0
553 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: gfunc
554 : REAL(KIND=dp), DIMENSION(:), POINTER :: w
555 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
556 : REAL(KIND=dp), INTENT(IN) :: gcut
557 : TYPE(pw_c1d_gs_type), INTENT(IN) :: rho_tot_g
558 : REAL(KIND=dp), DIMENSION(:), POINTER :: radii
559 : INTEGER, INTENT(IN) :: iw
560 : REAL(KIND=dp), INTENT(IN) :: Vol
561 :
562 : CHARACTER(len=*), PARAMETER :: routineN = 'ddapc_eval_AmI'
563 :
564 : INTEGER :: handle, ndim
565 : REAL(KIND=dp) :: condition_number, inv_error
566 296 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: AmE, cv
567 296 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: Am, AmI, Amw, g_dot_rvec_cos, &
568 296 : g_dot_rvec_sin
569 :
570 : !NB for precomputation of sin(g.r) and cos(g.r)
571 :
572 296 : CALL timeset(routineN, handle)
573 296 : ndim = SIZE(particle_set)*SIZE(radii)
574 1184 : ALLOCATE (Am(ndim, ndim))
575 888 : ALLOCATE (AmI(ndim, ndim))
576 888 : ALLOCATE (GAmI(ndim, ndim))
577 888 : ALLOCATE (cv(ndim))
578 296 : Am = 0.0_dp
579 296 : AmI = 0.0_dp
580 2834 : cv = 1.0_dp/Vol
581 : !NB precompute sin(g.r) and cos(g.r) for faster evaluation of cos(g.(r1-r2)) in build_A_matrix()
582 296 : CALL prep_g_dot_rvec_sin_cos(rho_tot_g, particle_set, gcut, g_dot_rvec_sin, g_dot_rvec_cos)
583 296 : CALL build_A_matrix(Am, gfunc, w, particle_set, radii, rho_tot_g, gcut, g_dot_rvec_sin, g_dot_rvec_cos)
584 296 : CALL cleanup_g_dot_rvec_sin_cos(g_dot_rvec_sin, g_dot_rvec_cos)
585 61100 : Am(:, :) = Am(:, :)/(Vol*Vol)
586 296 : CALL rho_tot_g%pw_grid%para%group%sum(Am)
587 296 : IF (iw > 0) THEN
588 : ! Checking conditions numbers and eigenvalues
589 0 : ALLOCATE (Amw(ndim, ndim))
590 0 : ALLOCATE (AmE(ndim))
591 0 : Amw(:, :) = Am
592 0 : CALL diamat_all(Amw, AmE)
593 0 : condition_number = MAXVAL(ABS(AmE))/MINVAL(ABS(AmE))
594 0 : WRITE (iw, '(T3,A)') " Eigenvalues of Matrix A:"
595 0 : WRITE (iw, '(T3,4E15.8)') AmE
596 0 : WRITE (iw, '(T3,A,1E15.9)') " Condition number:", condition_number
597 0 : IF (condition_number > 1.0E12_dp) THEN
598 : WRITE (iw, FMT="(/,T2,A)") &
599 0 : "WARNING: high condition number => possibly ill-conditioned matrix"
600 : END IF
601 0 : DEALLOCATE (Amw)
602 0 : DEALLOCATE (AmE)
603 : END IF
604 296 : CALL invert_matrix(Am, AmI, inv_error, "N", improve=.FALSE.)
605 296 : IF (iw > 0) THEN
606 0 : WRITE (iw, '(T3,A,F15.9)') " Error inverting the A matrix: ", inv_error
607 : END IF
608 122496 : c0 = DOT_PRODUCT(cv, MATMUL(AmI, cv))
609 296 : DEALLOCATE (Am)
610 296 : DEALLOCATE (cv)
611 61100 : GAmI = AmI
612 296 : DEALLOCATE (AmI)
613 296 : CALL timestop(handle)
614 592 : END SUBROUTINE ddapc_eval_AmI
615 :
616 : ! **************************************************************************************************
617 : !> \brief Evaluates the Ewald term E2 and E3 energy term for the decoupling/coupling
618 : !> of periodic images
619 : !> \param cp_para_env ...
620 : !> \param coeff ...
621 : !> \param factor ...
622 : !> \param cell ...
623 : !> \param multipole_section ...
624 : !> \param particle_set ...
625 : !> \param M ...
626 : !> \param radii ...
627 : !> \par History
628 : !> 08.2005 created [tlaino]
629 : !> \author Teodoro Laino
630 : !> \note NB receive cp_para_env for parallelization
631 : ! **************************************************************************************************
632 224 : RECURSIVE SUBROUTINE ewald_ddapc_pot(cp_para_env, coeff, factor, cell, multipole_section, &
633 : particle_set, M, radii)
634 : TYPE(mp_para_env_type), INTENT(IN) :: cp_para_env
635 : TYPE(pw_r3d_rs_type), INTENT(IN), POINTER :: coeff
636 : REAL(KIND=dp), INTENT(IN) :: factor
637 : TYPE(cell_type), POINTER :: cell
638 : TYPE(section_vals_type), POINTER :: multipole_section
639 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
640 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: M
641 : REAL(KIND=dp), DIMENSION(:), POINTER :: radii
642 :
643 : CHARACTER(len=*), PARAMETER :: routineN = 'ewald_ddapc_pot'
644 :
645 : INTEGER :: ewmdim, handle, iaxis, idim, idim1, idim2, idimo, igauss1, igauss2, ip1, ip2, &
646 : iparticle1, iparticle2, istart_g, k1, k2, k3, n_rep, ndim, r1, r2, r3
647 : INTEGER, DIMENSION(3) :: gmax, image_cell, rmax, rmin
648 : LOGICAL :: analyt
649 : REAL(KIND=dp) :: alpha, eps, ew_neut, fac, fac3, frac_radius, fs, g_ewald, galpha, gsq, &
650 : gsqi, ij_fac, my_val, r, r2tmp, r_ewald, rc1, rc12, rc2, rc22, rcut, rcut2, t1, tol, tol1
651 : REAL(KIND=dp), DIMENSION(3) :: g_index, gvec, ra, rvec, svec
652 224 : REAL(KIND=dp), DIMENSION(:), POINTER :: EwM
653 :
654 224 : NULLIFY (EwM)
655 224 : CALL timeset(routineN, handle)
656 224 : CPASSERT(.NOT. ASSOCIATED(M))
657 224 : CPASSERT(ASSOCIATED(radii))
658 2240 : rcut = MIN(NORM2(cell%hmat(:, 1)), NORM2(cell%hmat(:, 2)), NORM2(cell%hmat(:, 3)))/2.0_dp
659 224 : CALL section_vals_val_get(multipole_section, "RCUT", n_rep_val=n_rep)
660 224 : IF (n_rep == 1) CALL section_vals_val_get(multipole_section, "RCUT", r_val=rcut)
661 224 : CALL section_vals_val_get(multipole_section, "EWALD_PRECISION", r_val=eps)
662 224 : CALL section_vals_val_get(multipole_section, "ANALYTICAL_GTERM", l_val=analyt)
663 : ! The spline interpolation path is only valid for orthorhombic grids.
664 224 : analyt = analyt .OR. .NOT. cell%orthorhombic .OR. .NOT. ASSOCIATED(coeff)
665 224 : rcut2 = rcut**2
666 : !
667 : ! Setting-up parameters for Ewald summation
668 : !
669 224 : eps = MIN(ABS(eps), 0.5_dp)
670 224 : tol = SQRT(ABS(LOG(eps*rcut)))
671 224 : alpha = SQRT(ABS(LOG(eps*rcut*tol)))/rcut
672 224 : galpha = 1.0_dp/(4.0_dp*alpha*alpha)
673 224 : tol1 = SQRT(-LOG(eps*rcut*(2.0_dp*tol*alpha)**2))
674 224 : IF (cell%orthorhombic) THEN
675 856 : DO iaxis = 1, 3
676 856 : gmax(iaxis) = NINT(0.25_dp + cell%hmat(iaxis, iaxis)*alpha*tol1/pi)
677 : END DO
678 : ELSE
679 40 : DO iaxis = 1, 3
680 130 : gmax(iaxis) = CEILING(alpha*tol1*NORM2(cell%hmat(:, iaxis))/pi)
681 : END DO
682 : END IF
683 224 : fac = 1.e0_dp/cell%deth
684 224 : fac3 = fac*pi
685 224 : ew_neut = -fac*pi/alpha**2
686 : !
687 224 : ewmdim = SIZE(particle_set)*(SIZE(particle_set) + 1)/2
688 224 : ndim = SIZE(particle_set)*SIZE(radii)
689 672 : ALLOCATE (EwM(ewmdim))
690 896 : ALLOCATE (M(ndim, ndim))
691 92216 : M = 0.0_dp
692 : !
693 5584 : idim = 0
694 5584 : EwM = 0.0_dp
695 972 : DO iparticle1 = 1, SIZE(particle_set)
696 6108 : ip1 = (iparticle1 - 1)*SIZE(radii)
697 6332 : DO iparticle2 = 1, iparticle1
698 5360 : ij_fac = 1.0_dp
699 5360 : IF (iparticle1 == iparticle2) ij_fac = 0.5_dp
700 :
701 5360 : ip2 = (iparticle2 - 1)*SIZE(radii)
702 5360 : idim = idim + 1
703 : !NB parallelization, done here so indexing is right
704 5360 : IF (MOD(iparticle1, cp_para_env%num_pe) /= cp_para_env%mepos) CYCLE
705 : !
706 : ! Real-Space Contribution
707 : !
708 2704 : my_val = 0.0_dp
709 10816 : rvec = particle_set(iparticle1)%r - particle_set(iparticle2)%r
710 2704 : r_ewald = 0.0_dp
711 2704 : IF (iparticle1 /= iparticle2) THEN
712 2314 : ra = rvec
713 9256 : r2tmp = DOT_PRODUCT(ra, ra)
714 2314 : IF (r2tmp <= rcut2) THEN
715 2226 : r = SQRT(r2tmp)
716 2226 : t1 = erfc(alpha*r)/r
717 2226 : r_ewald = t1
718 : END IF
719 : END IF
720 35152 : svec = MATMUL(cell%h_inv, rvec)
721 10816 : DO iaxis = 1, 3
722 32448 : frac_radius = rcut*NORM2(cell%h_inv(iaxis, :))
723 8112 : rmin(iaxis) = FLOOR(-svec(iaxis) - frac_radius)
724 10816 : rmax(iaxis) = CEILING(-svec(iaxis) + frac_radius)
725 : END DO
726 16189 : DO r1 = rmin(1), rmax(1)
727 86728 : DO r2 = rmin(2), rmax(2)
728 467857 : DO r3 = rmin(3), rmax(3)
729 383833 : IF ((r1 == 0) .AND. (r2 == 0) .AND. (r3 == 0)) CYCLE
730 1524516 : image_cell = [r1, r2, r3]
731 7241451 : ra = rvec + MATMUL(cell%hmat, REAL(image_cell, KIND=dp))
732 1524516 : r2tmp = DOT_PRODUCT(ra, ra)
733 451668 : IF (r2tmp <= rcut2) THEN
734 47717 : r = SQRT(r2tmp)
735 47717 : t1 = erfc(alpha*r)/r
736 47717 : r_ewald = r_ewald + t1*ij_fac
737 : END IF
738 : END DO
739 : END DO
740 : END DO
741 : !
742 : ! G-space Contribution
743 : !
744 2704 : IF (analyt) THEN
745 278 : g_ewald = 0.0_dp
746 3535 : DO k1 = 0, gmax(1)
747 88788 : DO k2 = -gmax(2), gmax(2)
748 2516283 : DO k3 = -gmax(3), gmax(3)
749 2427773 : IF (k1 == 0 .AND. k2 == 0 .AND. k3 == 0) CYCLE
750 2427495 : fs = 2.0_dp; IF (k1 == 0) fs = 1.0_dp
751 9709980 : g_index = [REAL(k1, KIND=dp), REAL(k2, KIND=dp), REAL(k3, KIND=dp)]
752 12137475 : gvec = twopi*MATMUL(TRANSPOSE(cell%h_inv), g_index)
753 9709980 : gsq = DOT_PRODUCT(gvec, gvec)
754 2427495 : gsqi = fs/gsq
755 2427495 : t1 = fac*gsqi*EXP(-galpha*gsq)
756 9795511 : g_ewald = g_ewald + t1*COS(DOT_PRODUCT(gvec, rvec))
757 : END DO
758 : END DO
759 : END DO
760 : ELSE
761 2426 : g_ewald = Eval_Interp_Spl3_pbc(rvec, coeff)
762 : END IF
763 : !
764 : ! G-EWALD, R-EWALD
765 : !
766 2704 : g_ewald = r_ewald + fourpi*g_ewald
767 : !
768 : ! Self Contribution
769 : !
770 2704 : IF (iparticle1 == iparticle2) THEN
771 390 : g_ewald = g_ewald - 2.0_dp*alpha*oorootpi
772 : END IF
773 : !
774 2704 : IF (iparticle1 /= iparticle2) THEN
775 2314 : ra = rvec
776 9256 : r = NORM2(ra)
777 2314 : my_val = factor/r
778 : END IF
779 6108 : EwM(idim) = my_val - factor*g_ewald
780 : END DO ! iparticle2
781 : END DO ! iparticle1
782 : !NB sum over parallelized contributions of different nodes
783 10944 : CALL cp_para_env%sum(EwM)
784 224 : idim = 0
785 972 : DO iparticle2 = 1, SIZE(particle_set)
786 748 : ip2 = (iparticle2 - 1)*SIZE(radii)
787 748 : idimo = (iparticle2 - 1)
788 748 : idimo = idimo*(idimo + 1)/2
789 3216 : DO igauss2 = 1, SIZE(radii)
790 2244 : idim2 = ip2 + igauss2
791 2244 : rc2 = radii(igauss2)
792 2244 : rc22 = rc2*rc2
793 19072 : DO iparticle1 = 1, iparticle2
794 16080 : ip1 = (iparticle1 - 1)*SIZE(radii)
795 16080 : idim = idimo + iparticle1
796 16080 : istart_g = 1
797 16080 : IF (iparticle1 == iparticle2) istart_g = igauss2
798 64320 : DO igauss1 = istart_g, SIZE(radii)
799 45996 : idim1 = ip1 + igauss1
800 45996 : rc1 = radii(igauss1)
801 45996 : rc12 = rc1*rc1
802 45996 : M(idim1, idim2) = EwM(idim) - factor*ew_neut - factor*fac3*(rc12 + rc22)
803 62076 : M(idim2, idim1) = M(idim1, idim2)
804 : END DO
805 : END DO
806 : END DO ! iparticle2
807 : END DO ! iparticle1
808 224 : DEALLOCATE (EwM)
809 224 : CALL timestop(handle)
810 672 : END SUBROUTINE ewald_ddapc_pot
811 :
812 : ! **************************************************************************************************
813 : !> \brief Evaluates the electrostatic potential due to a simple solvation model
814 : !> Spherical cavity in a dieletric medium
815 : !> \param solvation_section ...
816 : !> \param particle_set ...
817 : !> \param M ...
818 : !> \param radii ...
819 : !> \par History
820 : !> 08.2006 created [tlaino]
821 : !> \author Teodoro Laino
822 : ! **************************************************************************************************
823 26 : SUBROUTINE solvation_ddapc_pot(solvation_section, particle_set, M, radii)
824 : TYPE(section_vals_type), POINTER :: solvation_section
825 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
826 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: M
827 : REAL(KIND=dp), DIMENSION(:), POINTER :: radii
828 :
829 : INTEGER :: i, idim, idim1, idim2, igauss1, igauss2, ip1, ip2, iparticle1, iparticle2, &
830 : istart_g, j, l, lmax, n_rep1, n_rep2, ndim, output_unit, weight
831 26 : INTEGER, DIMENSION(:), POINTER :: list
832 : LOGICAL :: fixed_center
833 : REAL(KIND=dp) :: center(3), eps_in, eps_out, factor, &
834 : mass, mycos, r1, r2, Rs, rvec(3)
835 26 : REAL(KIND=dp), DIMENSION(:), POINTER :: pos, R0
836 26 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: cost, LocP
837 :
838 26 : fixed_center = .FALSE.
839 52 : output_unit = cp_logger_get_default_io_unit()
840 26 : ndim = SIZE(particle_set)*SIZE(radii)
841 104 : ALLOCATE (M(ndim, ndim))
842 1274 : M = 0.0_dp
843 26 : eps_in = 1.0_dp
844 26 : CALL section_vals_val_get(solvation_section, "EPS_OUT", r_val=eps_out)
845 26 : CALL section_vals_val_get(solvation_section, "LMAX", i_val=lmax)
846 26 : CALL section_vals_val_get(solvation_section, "SPHERE%RADIUS", r_val=Rs)
847 26 : CALL section_vals_val_get(solvation_section, "SPHERE%CENTER%XYZ", n_rep_val=n_rep1)
848 26 : IF (n_rep1 /= 0) THEN
849 24 : CALL section_vals_val_get(solvation_section, "SPHERE%CENTER%XYZ", r_vals=R0)
850 96 : center = R0
851 : ELSE
852 : CALL section_vals_val_get(solvation_section, "SPHERE%CENTER%ATOM_LIST", &
853 2 : n_rep_val=n_rep2)
854 2 : IF (n_rep2 /= 0) THEN
855 2 : CALL section_vals_val_get(solvation_section, "SPHERE%CENTER%ATOM_LIST", i_vals=list)
856 2 : CALL section_vals_val_get(solvation_section, "SPHERE%CENTER%WEIGHT_TYPE", i_val=weight)
857 2 : ALLOCATE (R0(3))
858 : SELECT CASE (weight)
859 : CASE (weight_type_unit)
860 8 : R0 = 0.0_dp
861 4 : DO i = 1, SIZE(list)
862 10 : R0 = R0 + particle_set(list(i))%r
863 : END DO
864 8 : R0 = R0/REAL(SIZE(list), KIND=dp)
865 : CASE (weight_type_mass)
866 0 : R0 = 0.0_dp
867 0 : mass = 0.0_dp
868 0 : DO i = 1, SIZE(list)
869 0 : R0 = R0 + particle_set(list(i))%r*particle_set(list(i))%atomic_kind%mass
870 0 : mass = mass + particle_set(list(i))%atomic_kind%mass
871 : END DO
872 2 : R0 = R0/mass
873 : END SELECT
874 8 : center = R0
875 2 : CALL section_vals_val_get(solvation_section, "SPHERE%CENTER%FIXED", l_val=fixed_center)
876 4 : IF (fixed_center) THEN
877 : CALL section_vals_val_set(solvation_section, "SPHERE%CENTER%XYZ", &
878 2 : r_vals_ptr=R0)
879 : ELSE
880 0 : DEALLOCATE (R0)
881 : END IF
882 : END IF
883 : END IF
884 26 : CPASSERT(n_rep1 /= 0 .OR. n_rep2 /= 0)
885 : ! Potential calculation
886 104 : ALLOCATE (LocP(0:lmax, SIZE(particle_set)))
887 78 : ALLOCATE (pos(SIZE(particle_set)))
888 104 : ALLOCATE (cost(SIZE(particle_set), SIZE(particle_set)))
889 : ! Determining the single atomic contribution to the dielectric dipole
890 76 : DO i = 1, SIZE(particle_set)
891 200 : rvec = particle_set(i)%r - center
892 200 : r2 = DOT_PRODUCT(rvec, rvec)
893 50 : r1 = SQRT(r2)
894 50 : IF (r1 >= Rs) THEN
895 0 : IF (output_unit > 0) THEN
896 0 : WRITE (output_unit, '(A,I6,A)') "Atom number :: ", i, " is out of the solvation sphere"
897 0 : WRITE (output_unit, '(2(A,F12.6))') "Distance from the center::", r1, " Radius of the sphere::", rs
898 : END IF
899 0 : CPABORT("Unable to evaluate electrostatic potential in solution")
900 : END IF
901 250 : LocP(:, i) = 0.0_dp
902 50 : IF (r1 /= 0.0_dp) THEN
903 230 : DO l = 0, lmax
904 : LocP(l, i) = (r1**l*REAL(l + 1, KIND=dp)*(eps_in - eps_out))/ &
905 230 : (Rs**(2*l + 1)*eps_in*(REAL(l, KIND=dp)*eps_in + REAL(l + 1, KIND=dp)*eps_out))
906 : END DO
907 : ELSE
908 : ! limit for r->0
909 4 : LocP(0, i) = (eps_in - eps_out)/(Rs*eps_in*eps_out)
910 : END IF
911 76 : pos(i) = r1
912 : END DO
913 : ! Particle-Particle potential energy matrix
914 198 : cost = 0.0_dp
915 76 : DO i = 1, SIZE(particle_set)
916 162 : DO j = 1, i
917 86 : factor = 0.0_dp
918 86 : IF (pos(i)*pos(j) /= 0.0_dp) THEN
919 296 : mycos = DOT_PRODUCT(particle_set(i)%r - center, particle_set(j)%r - center)/(pos(i)*pos(j))
920 74 : IF (ABS(mycos) > 1.0_dp) mycos = SIGN(1.0_dp, mycos)
921 370 : DO l = 0, lmax
922 370 : factor = factor + LocP(l, i)*pos(j)**l*legendre(mycos, l, 0)
923 : END DO
924 : ELSE
925 12 : factor = LocP(0, i)
926 : END IF
927 86 : cost(i, j) = factor
928 136 : cost(j, i) = factor
929 : END DO
930 : END DO
931 : ! Computes the full potential energy matrix
932 26 : idim = 0
933 76 : DO iparticle2 = 1, SIZE(particle_set)
934 50 : ip2 = (iparticle2 - 1)*SIZE(radii)
935 226 : DO igauss2 = 1, SIZE(radii)
936 150 : idim2 = ip2 + igauss2
937 458 : DO iparticle1 = 1, iparticle2
938 258 : ip1 = (iparticle1 - 1)*SIZE(radii)
939 258 : istart_g = 1
940 258 : IF (iparticle1 == iparticle2) istart_g = igauss2
941 1032 : DO igauss1 = istart_g, SIZE(radii)
942 624 : idim1 = ip1 + igauss1
943 624 : M(idim1, idim2) = cost(iparticle1, iparticle2)
944 882 : M(idim2, idim1) = M(idim1, idim2)
945 : END DO
946 : END DO
947 : END DO
948 : END DO
949 26 : DEALLOCATE (cost)
950 26 : DEALLOCATE (pos)
951 26 : DEALLOCATE (LocP)
952 52 : END SUBROUTINE solvation_ddapc_pot
953 :
954 296 : END MODULE cp_ddapc_methods
|
