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 274 : 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 274 : CALL timeset(routineN, handle)
74 274 : gcut2 = gcut*gcut
75 : !
76 274 : s_dim = rho_tot_g%pw_grid%first_gne0
77 274 : e_dim = rho_tot_g%pw_grid%ngpts_cut_local
78 1096 : ALLOCATE (gfunc(s_dim:e_dim, SIZE(radii)))
79 822 : ALLOCATE (w(s_dim:e_dim))
80 66546619 : gfunc = 0.0_dp
81 22399623 : w = 0.0_dp
82 1048 : DO igauss = 1, SIZE(radii)
83 774 : rc = radii(igauss)
84 774 : rc2 = rc*rc
85 545002 : DO ig = s_dim, e_dim
86 544728 : g2 = rho_tot_g%pw_grid%gsq(ig)
87 544728 : IF (g2 > gcut2) EXIT
88 544728 : gfunc(ig, igauss) = EXP(-g2*rc2/4.0_dp)
89 : END DO
90 : END DO
91 182908 : DO ig = s_dim, e_dim
92 182908 : g2 = rho_tot_g%pw_grid%gsq(ig)
93 182908 : IF (g2 > gcut2) EXIT
94 182908 : w(ig) = fourpi*(g2 - gcut2)**2/(g2*gcut2)
95 : END DO
96 274 : CALL timestop(handle)
97 274 : 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 2128 : 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 2128 : REAL(KIND=dp), DIMENSION(:), POINTER :: my_bv, my_bvw
128 :
129 2128 : CALL timeset(routineN, handle)
130 2128 : NULLIFY (my_bv, my_bvw)
131 2128 : gcut2 = gcut*gcut
132 2128 : s_dim = rho_tot_g%pw_grid%first_gne0
133 2128 : e_dim = rho_tot_g%pw_grid%ngpts_cut_local
134 2128 : igmax = 0
135 965836 : DO ig = s_dim, e_dim
136 965836 : g2 = rho_tot_g%pw_grid%gsq(ig)
137 965836 : IF (g2 > gcut2) EXIT
138 965836 : igmax = ig
139 : END DO
140 2128 : IF (igmax >= s_dim) THEN
141 6384 : ALLOCATE (my_bv(s_dim:igmax))
142 4256 : ALLOCATE (my_bvw(s_dim:igmax))
143 : !
144 8212 : DO iparticle = 1, SIZE(particle_set)
145 24336 : rvec = particle_set(iparticle)%r
146 3337980 : my_bv = 0.0_dp
147 3337980 : DO ig = s_dim, igmax
148 13327584 : gvec = rho_tot_g%pw_grid%g(:, ig)
149 13327584 : arg = DOT_PRODUCT(gvec, rvec)
150 3331896 : phase = CMPLX(COS(arg), -SIN(arg), KIND=dp)
151 3337980 : my_bv(ig) = w(ig)*REAL(CONJG(rho_tot_g%array(ig))*phase, KIND=dp)
152 : END DO
153 24088 : DO igauss = 1, SIZE(radii)
154 15876 : idim = (iparticle - 1)*SIZE(radii) + igauss
155 9811188 : DO ig = s_dim, igmax
156 9811188 : my_bvw(ig) = my_bv(ig)*gfunc(ig, igauss)
157 : END DO
158 21960 : bv(idim) = accurate_sum(my_bvw)
159 : END DO
160 : END DO
161 2128 : DEALLOCATE (my_bvw)
162 2128 : 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 2128 : CALL timestop(handle)
172 2128 : 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 274 : 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 274 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: my_Am, my_Amw
207 274 : 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 274 : CALL timeset(routineN, handle)
212 274 : gcut2 = gcut*gcut
213 274 : s_dim = rho_tot_g%pw_grid%first_gne0
214 274 : e_dim = rho_tot_g%pw_grid%ngpts_cut_local
215 274 : igmax = 0
216 182908 : DO ig = s_dim, e_dim
217 182908 : g2 = rho_tot_g%pw_grid%gsq(ig)
218 182908 : IF (g2 > gcut2) EXIT
219 182908 : igmax = ig
220 : END DO
221 274 : IF (igmax >= s_dim) THEN
222 822 : ALLOCATE (my_Am(s_dim:igmax))
223 548 : ALLOCATE (my_Amw(s_dim:igmax))
224 1370 : ALLOCATE (gfunc_sq(s_dim:igmax, SIZE(radii), SIZE(radii)))
225 :
226 1048 : DO igauss1 = 1, SIZE(radii)
227 3322 : DO igauss2 = 1, SIZE(radii)
228 1630962 : 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 1100 : DO iparticle1 = 1, SIZE(particle_set)
233 4732 : DO iparticle2 = iparticle1, SIZE(particle_set)
234 6949714 : 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 6949714 : g_dot_rvec_sin(ig - s_dim + 1, iparticle1)*g_dot_rvec_sin(ig - s_dim + 1, iparticle2))
238 : END DO
239 15174 : DO igauss1 = 1, SIZE(radii)
240 10716 : idim1 = (iparticle1 - 1)*SIZE(radii) + igauss1
241 10716 : istart_g = 1
242 10716 : IF (iparticle2 == iparticle1) istart_g = igauss1
243 44000 : DO igauss2 = istart_g, SIZE(radii)
244 29652 : idim2 = (iparticle2 - 1)*SIZE(radii) + igauss2
245 60291900 : 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 60291900 : tmp = SUM(my_Amw)
249 29652 : Am(idim2, idim1) = tmp
250 40368 : Am(idim1, idim2) = tmp
251 : END DO
252 : END DO
253 : END DO
254 : END DO
255 274 : DEALLOCATE (gfunc_sq)
256 274 : DEALLOCATE (my_Amw)
257 274 : DEALLOCATE (my_Am)
258 : END IF
259 274 : CALL timestop(handle)
260 274 : 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 : END IF
468 148 : CALL timestop(handle)
469 148 : END SUBROUTINE build_der_A_matrix_rows
470 :
471 : ! **************************************************************************************************
472 : !> \brief deallocate g_dot_rvec_* arrays
473 : !> \param g_dot_rvec_sin ...
474 : !> \param g_dot_rvec_cos ...
475 : ! **************************************************************************************************
476 422 : SUBROUTINE cleanup_g_dot_rvec_sin_cos(g_dot_rvec_sin, g_dot_rvec_cos)
477 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: g_dot_rvec_sin, g_dot_rvec_cos
478 :
479 422 : IF (ALLOCATED(g_dot_rvec_sin)) DEALLOCATE (g_dot_rvec_sin)
480 422 : IF (ALLOCATED(g_dot_rvec_cos)) DEALLOCATE (g_dot_rvec_cos)
481 422 : END SUBROUTINE cleanup_g_dot_rvec_sin_cos
482 :
483 : ! **************************************************************************************************
484 : !> \brief precompute sin(g.r) and cos(g.r) for quicker evaluations of sin(g.(r1-r2)) and cos(g.(r1-r2))
485 : !> \param rho_tot_g ...
486 : !> \param particle_set ...
487 : !> \param gcut ...
488 : !> \param g_dot_rvec_sin ...
489 : !> \param g_dot_rvec_cos ...
490 : ! **************************************************************************************************
491 422 : SUBROUTINE prep_g_dot_rvec_sin_cos(rho_tot_g, particle_set, gcut, g_dot_rvec_sin, g_dot_rvec_cos)
492 : TYPE(pw_c1d_gs_type), INTENT(IN) :: rho_tot_g
493 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
494 : REAL(KIND=dp), INTENT(IN) :: gcut
495 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: g_dot_rvec_sin, g_dot_rvec_cos
496 :
497 : INTEGER :: e_dim, ig, igmax, iparticle, s_dim
498 : REAL(KIND=dp) :: g2, g_dot_rvec, gcut2, rvec(3)
499 :
500 422 : gcut2 = gcut*gcut
501 422 : s_dim = rho_tot_g%pw_grid%first_gne0
502 422 : e_dim = rho_tot_g%pw_grid%ngpts_cut_local
503 422 : igmax = 0
504 335090 : DO ig = s_dim, e_dim
505 335090 : g2 = rho_tot_g%pw_grid%gsq(ig)
506 335090 : IF (g2 > gcut2) EXIT
507 335090 : igmax = ig
508 : END DO
509 :
510 422 : IF (igmax >= s_dim) THEN
511 1688 : ALLOCATE (g_dot_rvec_sin(1:igmax - s_dim + 1, SIZE(particle_set)))
512 1266 : ALLOCATE (g_dot_rvec_cos(1:igmax - s_dim + 1, SIZE(particle_set)))
513 :
514 1668 : DO iparticle = 1, SIZE(particle_set)
515 4984 : rvec = particle_set(iparticle)%r
516 1096274 : DO ig = s_dim, igmax
517 4378424 : g_dot_rvec = DOT_PRODUCT(rho_tot_g%pw_grid%g(:, ig), rvec)
518 1094606 : g_dot_rvec_sin(ig - s_dim + 1, iparticle) = SIN(g_dot_rvec)
519 1095852 : g_dot_rvec_cos(ig - s_dim + 1, iparticle) = COS(g_dot_rvec)
520 : END DO
521 : END DO
522 : END IF
523 :
524 422 : END SUBROUTINE prep_g_dot_rvec_sin_cos
525 :
526 : ! **************************************************************************************************
527 : !> \brief Computes the inverse AmI of the Am matrix
528 : !> \param GAmI ...
529 : !> \param c0 ...
530 : !> \param gfunc ...
531 : !> \param w ...
532 : !> \param particle_set ...
533 : !> \param gcut ...
534 : !> \param rho_tot_g ...
535 : !> \param radii ...
536 : !> \param iw ...
537 : !> \param Vol ...
538 : !> \par History
539 : !> 12.2005 created [tlaino]
540 : !> \author Teodoro Laino
541 : ! **************************************************************************************************
542 274 : SUBROUTINE ddapc_eval_AmI(GAmI, c0, gfunc, w, particle_set, gcut, &
543 : rho_tot_g, radii, iw, Vol)
544 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: GAmI
545 : REAL(KIND=dp), INTENT(OUT) :: c0
546 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: gfunc
547 : REAL(KIND=dp), DIMENSION(:), POINTER :: w
548 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
549 : REAL(KIND=dp), INTENT(IN) :: gcut
550 : TYPE(pw_c1d_gs_type), INTENT(IN) :: rho_tot_g
551 : REAL(KIND=dp), DIMENSION(:), POINTER :: radii
552 : INTEGER, INTENT(IN) :: iw
553 : REAL(KIND=dp), INTENT(IN) :: Vol
554 :
555 : CHARACTER(len=*), PARAMETER :: routineN = 'ddapc_eval_AmI'
556 :
557 : INTEGER :: handle, ndim
558 : REAL(KIND=dp) :: condition_number, inv_error
559 274 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: AmE, cv
560 274 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: Am, AmI, Amw, g_dot_rvec_cos, &
561 274 : g_dot_rvec_sin
562 :
563 : !NB for precomputation of sin(g.r) and cos(g.r)
564 :
565 274 : CALL timeset(routineN, handle)
566 274 : ndim = SIZE(particle_set)*SIZE(radii)
567 1096 : ALLOCATE (Am(ndim, ndim))
568 822 : ALLOCATE (AmI(ndim, ndim))
569 822 : ALLOCATE (GAmI(ndim, ndim))
570 822 : ALLOCATE (cv(ndim))
571 274 : Am = 0.0_dp
572 274 : AmI = 0.0_dp
573 2644 : cv = 1.0_dp/Vol
574 : !NB precompute sin(g.r) and cos(g.r) for faster evaluation of cos(g.(r1-r2)) in build_A_matrix()
575 274 : CALL prep_g_dot_rvec_sin_cos(rho_tot_g, particle_set, gcut, g_dot_rvec_sin, g_dot_rvec_cos)
576 274 : CALL build_A_matrix(Am, gfunc, w, particle_set, radii, rho_tot_g, gcut, g_dot_rvec_sin, g_dot_rvec_cos)
577 274 : CALL cleanup_g_dot_rvec_sin_cos(g_dot_rvec_sin, g_dot_rvec_cos)
578 59578 : Am(:, :) = Am(:, :)/(Vol*Vol)
579 274 : CALL rho_tot_g%pw_grid%para%group%sum(Am)
580 274 : IF (iw > 0) THEN
581 : ! Checking conditions numbers and eigenvalues
582 0 : ALLOCATE (Amw(ndim, ndim))
583 0 : ALLOCATE (AmE(ndim))
584 0 : Amw(:, :) = Am
585 0 : CALL diamat_all(Amw, AmE)
586 0 : condition_number = MAXVAL(ABS(AmE))/MINVAL(ABS(AmE))
587 0 : WRITE (iw, '(T3,A)') " Eigenvalues of Matrix A:"
588 0 : WRITE (iw, '(T3,4E15.8)') AmE
589 0 : WRITE (iw, '(T3,A,1E15.9)') " Condition number:", condition_number
590 0 : IF (condition_number > 1.0E12_dp) THEN
591 : WRITE (iw, FMT="(/,T2,A)") &
592 0 : "WARNING: high condition number => possibly ill-conditioned matrix"
593 : END IF
594 0 : DEALLOCATE (Amw)
595 0 : DEALLOCATE (AmE)
596 : END IF
597 274 : CALL invert_matrix(Am, AmI, inv_error, "N", improve=.FALSE.)
598 274 : IF (iw > 0) THEN
599 0 : WRITE (iw, '(T3,A,F15.9)') " Error inverting the A matrix: ", inv_error
600 : END IF
601 119430 : c0 = DOT_PRODUCT(cv, MATMUL(AmI, cv))
602 274 : DEALLOCATE (Am)
603 274 : DEALLOCATE (cv)
604 59578 : GAmI = AmI
605 274 : DEALLOCATE (AmI)
606 274 : CALL timestop(handle)
607 548 : END SUBROUTINE ddapc_eval_AmI
608 :
609 : ! **************************************************************************************************
610 : !> \brief Evaluates the Ewald term E2 and E3 energy term for the decoupling/coupling
611 : !> of periodic images
612 : !> \param cp_para_env ...
613 : !> \param coeff ...
614 : !> \param factor ...
615 : !> \param cell ...
616 : !> \param multipole_section ...
617 : !> \param particle_set ...
618 : !> \param M ...
619 : !> \param radii ...
620 : !> \par History
621 : !> 08.2005 created [tlaino]
622 : !> \author Teodoro Laino
623 : !> \note NB receive cp_para_env for parallelization
624 : ! **************************************************************************************************
625 202 : RECURSIVE SUBROUTINE ewald_ddapc_pot(cp_para_env, coeff, factor, cell, multipole_section, &
626 : particle_set, M, radii)
627 : TYPE(mp_para_env_type), INTENT(IN) :: cp_para_env
628 : TYPE(pw_r3d_rs_type), INTENT(IN), POINTER :: coeff
629 : REAL(KIND=dp), INTENT(IN) :: factor
630 : TYPE(cell_type), POINTER :: cell
631 : TYPE(section_vals_type), POINTER :: multipole_section
632 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
633 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: M
634 : REAL(KIND=dp), DIMENSION(:), POINTER :: radii
635 :
636 : CHARACTER(len=*), PARAMETER :: routineN = 'ewald_ddapc_pot'
637 :
638 : INTEGER :: ewmdim, handle, iaxis, idim, idim1, idim2, idimo, igauss1, igauss2, ip1, ip2, &
639 : iparticle1, iparticle2, istart_g, k1, k2, k3, n_rep, ndim, r1, r2, r3
640 : INTEGER, DIMENSION(3) :: gmax, image_cell, rmax, rmin
641 : LOGICAL :: analyt
642 : REAL(KIND=dp) :: alpha, eps, ew_neut, fac, fac3, frac_radius, fs, g_ewald, galpha, gsq, &
643 : gsqi, ij_fac, my_val, r, r2tmp, r_ewald, rc1, rc12, rc2, rc22, rcut, rcut2, t1, tol, tol1
644 : REAL(KIND=dp), DIMENSION(3) :: g_index, gvec, ra, rvec, svec
645 202 : REAL(KIND=dp), DIMENSION(:), POINTER :: EwM
646 :
647 202 : NULLIFY (EwM)
648 202 : CALL timeset(routineN, handle)
649 202 : CPASSERT(.NOT. ASSOCIATED(M))
650 202 : CPASSERT(ASSOCIATED(radii))
651 2020 : rcut = MIN(NORM2(cell%hmat(:, 1)), NORM2(cell%hmat(:, 2)), NORM2(cell%hmat(:, 3)))/2.0_dp
652 202 : CALL section_vals_val_get(multipole_section, "RCUT", n_rep_val=n_rep)
653 202 : IF (n_rep == 1) CALL section_vals_val_get(multipole_section, "RCUT", r_val=rcut)
654 202 : CALL section_vals_val_get(multipole_section, "EWALD_PRECISION", r_val=eps)
655 202 : CALL section_vals_val_get(multipole_section, "ANALYTICAL_GTERM", l_val=analyt)
656 : ! The spline interpolation path is only valid for orthorhombic grids.
657 202 : analyt = analyt .OR. .NOT. cell%orthorhombic .OR. .NOT. ASSOCIATED(coeff)
658 202 : rcut2 = rcut**2
659 : !
660 : ! Setting-up parameters for Ewald summation
661 : !
662 202 : eps = MIN(ABS(eps), 0.5_dp)
663 202 : tol = SQRT(ABS(LOG(eps*rcut)))
664 202 : alpha = SQRT(ABS(LOG(eps*rcut*tol)))/rcut
665 202 : galpha = 1.0_dp/(4.0_dp*alpha*alpha)
666 202 : tol1 = SQRT(-LOG(eps*rcut*(2.0_dp*tol*alpha)**2))
667 202 : IF (cell%orthorhombic) THEN
668 768 : DO iaxis = 1, 3
669 768 : gmax(iaxis) = NINT(0.25_dp + cell%hmat(iaxis, iaxis)*alpha*tol1/pi)
670 : END DO
671 : ELSE
672 40 : DO iaxis = 1, 3
673 130 : gmax(iaxis) = CEILING(alpha*tol1*NORM2(cell%hmat(:, iaxis))/pi)
674 : END DO
675 : END IF
676 202 : fac = 1.e0_dp/cell%deth
677 202 : fac3 = fac*pi
678 202 : ew_neut = -fac*pi/alpha**2
679 : !
680 202 : ewmdim = SIZE(particle_set)*(SIZE(particle_set) + 1)/2
681 202 : ndim = SIZE(particle_set)*SIZE(radii)
682 606 : ALLOCATE (EwM(ewmdim))
683 808 : ALLOCATE (M(ndim, ndim))
684 90694 : M = 0.0_dp
685 : !
686 5460 : idim = 0
687 5460 : EwM = 0.0_dp
688 894 : DO iparticle1 = 1, SIZE(particle_set)
689 5950 : ip1 = (iparticle1 - 1)*SIZE(radii)
690 6152 : DO iparticle2 = 1, iparticle1
691 5258 : ij_fac = 1.0_dp
692 5258 : IF (iparticle1 == iparticle2) ij_fac = 0.5_dp
693 :
694 5258 : ip2 = (iparticle2 - 1)*SIZE(radii)
695 5258 : idim = idim + 1
696 : !NB parallelization, done here so indexing is right
697 5258 : IF (MOD(iparticle1, cp_para_env%num_pe) /= cp_para_env%mepos) CYCLE
698 : !
699 : ! Real-Space Contribution
700 : !
701 2653 : my_val = 0.0_dp
702 10612 : rvec = particle_set(iparticle1)%r - particle_set(iparticle2)%r
703 2653 : r_ewald = 0.0_dp
704 2653 : IF (iparticle1 /= iparticle2) THEN
705 2291 : ra = rvec
706 9164 : r2tmp = DOT_PRODUCT(ra, ra)
707 2291 : IF (r2tmp <= rcut2) THEN
708 2203 : r = SQRT(r2tmp)
709 2203 : t1 = erfc(alpha*r)/r
710 2203 : r_ewald = t1
711 : END IF
712 : END IF
713 34489 : svec = MATMUL(cell%h_inv, rvec)
714 10612 : DO iaxis = 1, 3
715 31836 : frac_radius = rcut*NORM2(cell%h_inv(iaxis, :))
716 7959 : rmin(iaxis) = FLOOR(-svec(iaxis) - frac_radius)
717 10612 : rmax(iaxis) = CEILING(-svec(iaxis) + frac_radius)
718 : END DO
719 15985 : DO r1 = rmin(1), rmax(1)
720 86065 : DO r2 = rmin(2), rmax(2)
721 465868 : DO r3 = rmin(3), rmax(3)
722 382456 : IF ((r1 == 0) .AND. (r2 == 0) .AND. (r3 == 0)) CYCLE
723 1519212 : image_cell = [r1, r2, r3]
724 7216257 : ra = rvec + MATMUL(cell%hmat, REAL(image_cell, KIND=dp))
725 1519212 : r2tmp = DOT_PRODUCT(ra, ra)
726 449883 : IF (r2tmp <= rcut2) THEN
727 47717 : r = SQRT(r2tmp)
728 47717 : t1 = erfc(alpha*r)/r
729 47717 : r_ewald = r_ewald + t1*ij_fac
730 : END IF
731 : END DO
732 : END DO
733 : END DO
734 : !
735 : ! G-space Contribution
736 : !
737 2653 : IF (analyt) THEN
738 278 : g_ewald = 0.0_dp
739 3535 : DO k1 = 0, gmax(1)
740 88788 : DO k2 = -gmax(2), gmax(2)
741 2516283 : DO k3 = -gmax(3), gmax(3)
742 2427773 : IF (k1 == 0 .AND. k2 == 0 .AND. k3 == 0) CYCLE
743 2427495 : fs = 2.0_dp; IF (k1 == 0) fs = 1.0_dp
744 9709980 : g_index = [REAL(k1, KIND=dp), REAL(k2, KIND=dp), REAL(k3, KIND=dp)]
745 12137475 : gvec = twopi*MATMUL(TRANSPOSE(cell%h_inv), g_index)
746 9709980 : gsq = DOT_PRODUCT(gvec, gvec)
747 2427495 : gsqi = fs/gsq
748 2427495 : t1 = fac*gsqi*EXP(-galpha*gsq)
749 9795511 : g_ewald = g_ewald + t1*COS(DOT_PRODUCT(gvec, rvec))
750 : END DO
751 : END DO
752 : END DO
753 : ELSE
754 2375 : g_ewald = Eval_Interp_Spl3_pbc(rvec, coeff)
755 : END IF
756 : !
757 : ! G-EWALD, R-EWALD
758 : !
759 2653 : g_ewald = r_ewald + fourpi*g_ewald
760 : !
761 : ! Self Contribution
762 : !
763 2653 : IF (iparticle1 == iparticle2) THEN
764 362 : g_ewald = g_ewald - 2.0_dp*alpha*oorootpi
765 : END IF
766 : !
767 2653 : IF (iparticle1 /= iparticle2) THEN
768 2291 : ra = rvec
769 9164 : r = SQRT(DOT_PRODUCT(ra, ra))
770 2291 : my_val = factor/r
771 : END IF
772 5950 : EwM(idim) = my_val - factor*g_ewald
773 : END DO ! iparticle2
774 : END DO ! iparticle1
775 : !NB sum over parallelized contributions of different nodes
776 10718 : CALL cp_para_env%sum(EwM)
777 202 : idim = 0
778 894 : DO iparticle2 = 1, SIZE(particle_set)
779 692 : ip2 = (iparticle2 - 1)*SIZE(radii)
780 692 : idimo = (iparticle2 - 1)
781 692 : idimo = idimo*(idimo + 1)/2
782 2970 : DO igauss2 = 1, SIZE(radii)
783 2076 : idim2 = ip2 + igauss2
784 2076 : rc2 = radii(igauss2)
785 2076 : rc22 = rc2*rc2
786 18542 : DO iparticle1 = 1, iparticle2
787 15774 : ip1 = (iparticle1 - 1)*SIZE(radii)
788 15774 : idim = idimo + iparticle1
789 15774 : istart_g = 1
790 15774 : IF (iparticle1 == iparticle2) istart_g = igauss2
791 63096 : DO igauss1 = istart_g, SIZE(radii)
792 45246 : idim1 = ip1 + igauss1
793 45246 : rc1 = radii(igauss1)
794 45246 : rc12 = rc1*rc1
795 45246 : M(idim1, idim2) = EwM(idim) - factor*ew_neut - factor*fac3*(rc12 + rc22)
796 61020 : M(idim2, idim1) = M(idim1, idim2)
797 : END DO
798 : END DO
799 : END DO ! iparticle2
800 : END DO ! iparticle1
801 202 : DEALLOCATE (EwM)
802 202 : CALL timestop(handle)
803 606 : END SUBROUTINE ewald_ddapc_pot
804 :
805 : ! **************************************************************************************************
806 : !> \brief Evaluates the electrostatic potential due to a simple solvation model
807 : !> Spherical cavity in a dieletric medium
808 : !> \param solvation_section ...
809 : !> \param particle_set ...
810 : !> \param M ...
811 : !> \param radii ...
812 : !> \par History
813 : !> 08.2006 created [tlaino]
814 : !> \author Teodoro Laino
815 : ! **************************************************************************************************
816 26 : SUBROUTINE solvation_ddapc_pot(solvation_section, particle_set, M, radii)
817 : TYPE(section_vals_type), POINTER :: solvation_section
818 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
819 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: M
820 : REAL(KIND=dp), DIMENSION(:), POINTER :: radii
821 :
822 : INTEGER :: i, idim, idim1, idim2, igauss1, igauss2, ip1, ip2, iparticle1, iparticle2, &
823 : istart_g, j, l, lmax, n_rep1, n_rep2, ndim, output_unit, weight
824 26 : INTEGER, DIMENSION(:), POINTER :: list
825 : LOGICAL :: fixed_center
826 : REAL(KIND=dp) :: center(3), eps_in, eps_out, factor, &
827 : mass, mycos, r1, r2, Rs, rvec(3)
828 26 : REAL(KIND=dp), DIMENSION(:), POINTER :: pos, R0
829 26 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: cost, LocP
830 :
831 26 : fixed_center = .FALSE.
832 52 : output_unit = cp_logger_get_default_io_unit()
833 26 : ndim = SIZE(particle_set)*SIZE(radii)
834 104 : ALLOCATE (M(ndim, ndim))
835 1274 : M = 0.0_dp
836 26 : eps_in = 1.0_dp
837 26 : CALL section_vals_val_get(solvation_section, "EPS_OUT", r_val=eps_out)
838 26 : CALL section_vals_val_get(solvation_section, "LMAX", i_val=lmax)
839 26 : CALL section_vals_val_get(solvation_section, "SPHERE%RADIUS", r_val=Rs)
840 26 : CALL section_vals_val_get(solvation_section, "SPHERE%CENTER%XYZ", n_rep_val=n_rep1)
841 26 : IF (n_rep1 /= 0) THEN
842 24 : CALL section_vals_val_get(solvation_section, "SPHERE%CENTER%XYZ", r_vals=R0)
843 96 : center = R0
844 : ELSE
845 : CALL section_vals_val_get(solvation_section, "SPHERE%CENTER%ATOM_LIST", &
846 2 : n_rep_val=n_rep2)
847 2 : IF (n_rep2 /= 0) THEN
848 2 : CALL section_vals_val_get(solvation_section, "SPHERE%CENTER%ATOM_LIST", i_vals=list)
849 2 : CALL section_vals_val_get(solvation_section, "SPHERE%CENTER%WEIGHT_TYPE", i_val=weight)
850 2 : ALLOCATE (R0(3))
851 : SELECT CASE (weight)
852 : CASE (weight_type_unit)
853 8 : R0 = 0.0_dp
854 4 : DO i = 1, SIZE(list)
855 10 : R0 = R0 + particle_set(list(i))%r
856 : END DO
857 8 : R0 = R0/REAL(SIZE(list), KIND=dp)
858 : CASE (weight_type_mass)
859 0 : R0 = 0.0_dp
860 0 : mass = 0.0_dp
861 0 : DO i = 1, SIZE(list)
862 0 : R0 = R0 + particle_set(list(i))%r*particle_set(list(i))%atomic_kind%mass
863 0 : mass = mass + particle_set(list(i))%atomic_kind%mass
864 : END DO
865 2 : R0 = R0/mass
866 : END SELECT
867 8 : center = R0
868 2 : CALL section_vals_val_get(solvation_section, "SPHERE%CENTER%FIXED", l_val=fixed_center)
869 4 : IF (fixed_center) THEN
870 : CALL section_vals_val_set(solvation_section, "SPHERE%CENTER%XYZ", &
871 2 : r_vals_ptr=R0)
872 : ELSE
873 0 : DEALLOCATE (R0)
874 : END IF
875 : END IF
876 : END IF
877 26 : CPASSERT(n_rep1 /= 0 .OR. n_rep2 /= 0)
878 : ! Potential calculation
879 104 : ALLOCATE (LocP(0:lmax, SIZE(particle_set)))
880 78 : ALLOCATE (pos(SIZE(particle_set)))
881 104 : ALLOCATE (cost(SIZE(particle_set), SIZE(particle_set)))
882 : ! Determining the single atomic contribution to the dielectric dipole
883 76 : DO i = 1, SIZE(particle_set)
884 200 : rvec = particle_set(i)%r - center
885 200 : r2 = DOT_PRODUCT(rvec, rvec)
886 50 : r1 = SQRT(r2)
887 50 : IF (r1 >= Rs) THEN
888 0 : IF (output_unit > 0) THEN
889 0 : WRITE (output_unit, '(A,I6,A)') "Atom number :: ", i, " is out of the solvation sphere"
890 0 : WRITE (output_unit, '(2(A,F12.6))') "Distance from the center::", r1, " Radius of the sphere::", rs
891 : END IF
892 0 : CPABORT("")
893 : END IF
894 250 : LocP(:, i) = 0.0_dp
895 50 : IF (r1 /= 0.0_dp) THEN
896 230 : DO l = 0, lmax
897 : LocP(l, i) = (r1**l*REAL(l + 1, KIND=dp)*(eps_in - eps_out))/ &
898 230 : (Rs**(2*l + 1)*eps_in*(REAL(l, KIND=dp)*eps_in + REAL(l + 1, KIND=dp)*eps_out))
899 : END DO
900 : ELSE
901 : ! limit for r->0
902 4 : LocP(0, i) = (eps_in - eps_out)/(Rs*eps_in*eps_out)
903 : END IF
904 76 : pos(i) = r1
905 : END DO
906 : ! Particle-Particle potential energy matrix
907 198 : cost = 0.0_dp
908 76 : DO i = 1, SIZE(particle_set)
909 162 : DO j = 1, i
910 86 : factor = 0.0_dp
911 86 : IF (pos(i)*pos(j) /= 0.0_dp) THEN
912 296 : mycos = DOT_PRODUCT(particle_set(i)%r - center, particle_set(j)%r - center)/(pos(i)*pos(j))
913 74 : IF (ABS(mycos) > 1.0_dp) mycos = SIGN(1.0_dp, mycos)
914 370 : DO l = 0, lmax
915 370 : factor = factor + LocP(l, i)*pos(j)**l*legendre(mycos, l, 0)
916 : END DO
917 : ELSE
918 12 : factor = LocP(0, i)
919 : END IF
920 86 : cost(i, j) = factor
921 136 : cost(j, i) = factor
922 : END DO
923 : END DO
924 : ! Computes the full potential energy matrix
925 26 : idim = 0
926 76 : DO iparticle2 = 1, SIZE(particle_set)
927 50 : ip2 = (iparticle2 - 1)*SIZE(radii)
928 226 : DO igauss2 = 1, SIZE(radii)
929 150 : idim2 = ip2 + igauss2
930 458 : DO iparticle1 = 1, iparticle2
931 258 : ip1 = (iparticle1 - 1)*SIZE(radii)
932 258 : istart_g = 1
933 258 : IF (iparticle1 == iparticle2) istart_g = igauss2
934 1032 : DO igauss1 = istart_g, SIZE(radii)
935 624 : idim1 = ip1 + igauss1
936 624 : M(idim1, idim2) = cost(iparticle1, iparticle2)
937 882 : M(idim2, idim1) = M(idim1, idim2)
938 : END DO
939 : END DO
940 : END DO
941 : END DO
942 26 : DEALLOCATE (cost)
943 26 : DEALLOCATE (pos)
944 26 : DEALLOCATE (LocP)
945 52 : END SUBROUTINE solvation_ddapc_pot
946 :
947 274 : END MODULE cp_ddapc_methods
|
