Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 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 246 : 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 246 : CALL timeset(routineN, handle)
74 246 : gcut2 = gcut*gcut
75 : !
76 246 : s_dim = rho_tot_g%pw_grid%first_gne0
77 246 : e_dim = rho_tot_g%pw_grid%ngpts_cut_local
78 984 : ALLOCATE (gfunc(s_dim:e_dim, SIZE(radii)))
79 738 : ALLOCATE (w(s_dim:e_dim))
80 54440205 : gfunc = 0.0_dp
81 18364161 : w = 0.0_dp
82 936 : DO igauss = 1, SIZE(radii)
83 690 : rc = radii(igauss)
84 690 : rc2 = rc*rc
85 521214 : DO ig = s_dim, e_dim
86 520968 : g2 = rho_tot_g%pw_grid%gsq(ig)
87 520968 : IF (g2 > gcut2) EXIT
88 520968 : gfunc(ig, igauss) = EXP(-g2*rc2/4.0_dp)
89 : END DO
90 : END DO
91 174988 : DO ig = s_dim, e_dim
92 174988 : g2 = rho_tot_g%pw_grid%gsq(ig)
93 174988 : IF (g2 > gcut2) EXIT
94 174988 : w(ig) = fourpi*(g2 - gcut2)**2/(g2*gcut2)
95 : END DO
96 246 : CALL timestop(handle)
97 246 : 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 1342 : 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 1342 : REAL(KIND=dp), DIMENSION(:), POINTER :: my_bv, my_bvw
128 :
129 1342 : CALL timeset(routineN, handle)
130 1342 : NULLIFY (my_bv, my_bvw)
131 1342 : gcut2 = gcut*gcut
132 1342 : s_dim = rho_tot_g%pw_grid%first_gne0
133 1342 : e_dim = rho_tot_g%pw_grid%ngpts_cut_local
134 1342 : igmax = 0
135 821400 : DO ig = s_dim, e_dim
136 821400 : g2 = rho_tot_g%pw_grid%gsq(ig)
137 821400 : IF (g2 > gcut2) EXIT
138 821400 : igmax = ig
139 : END DO
140 1342 : IF (igmax .GE. s_dim) THEN
141 4026 : ALLOCATE (my_bv(s_dim:igmax))
142 4026 : ALLOCATE (my_bvw(s_dim:igmax))
143 : !
144 5196 : DO iparticle = 1, SIZE(particle_set)
145 15416 : rvec = particle_set(iparticle)%r
146 2916828 : my_bv = 0.0_dp
147 2916828 : DO ig = s_dim, igmax
148 11651896 : gvec = rho_tot_g%pw_grid%g(:, ig)
149 11651896 : arg = DOT_PRODUCT(gvec, rvec)
150 2912974 : phase = CMPLX(COS(arg), -SIN(arg), KIND=dp)
151 2916828 : my_bv(ig) = w(ig)*REAL(CONJG(rho_tot_g%array(ig))*phase, KIND=dp)
152 : END DO
153 15354 : DO igauss = 1, SIZE(radii)
154 10158 : idim = (iparticle - 1)*SIZE(radii) + igauss
155 8630676 : DO ig = s_dim, igmax
156 8630676 : my_bvw(ig) = my_bv(ig)*gfunc(ig, igauss)
157 : END DO
158 14012 : bv(idim) = accurate_sum(my_bvw)
159 : END DO
160 : END DO
161 1342 : DEALLOCATE (my_bvw)
162 1342 : 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 1342 : CALL timestop(handle)
172 1342 : 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 246 : 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 246 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: my_Am, my_Amw
207 246 : 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 246 : CALL timeset(routineN, handle)
212 246 : gcut2 = gcut*gcut
213 246 : s_dim = rho_tot_g%pw_grid%first_gne0
214 246 : e_dim = rho_tot_g%pw_grid%ngpts_cut_local
215 246 : igmax = 0
216 174988 : DO ig = s_dim, e_dim
217 174988 : g2 = rho_tot_g%pw_grid%gsq(ig)
218 174988 : IF (g2 > gcut2) EXIT
219 174988 : igmax = ig
220 : END DO
221 246 : IF (igmax .GE. s_dim) THEN
222 738 : ALLOCATE (my_Am(s_dim:igmax))
223 738 : ALLOCATE (my_Amw(s_dim:igmax))
224 1230 : ALLOCATE (gfunc_sq(s_dim:igmax, SIZE(radii), SIZE(radii)))
225 :
226 936 : DO igauss1 = 1, SIZE(radii)
227 2958 : DO igauss2 = 1, SIZE(radii)
228 1559598 : 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 992 : DO iparticle1 = 1, SIZE(particle_set)
233 4468 : DO iparticle2 = iparticle1, SIZE(particle_set)
234 6903904 : 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 6903904 : g_dot_rvec_sin(ig - s_dim + 1, iparticle1)*g_dot_rvec_sin(ig - s_dim + 1, iparticle2))
238 : END DO
239 14470 : DO igauss1 = 1, SIZE(radii)
240 10248 : idim1 = (iparticle1 - 1)*SIZE(radii) + igauss1
241 10248 : istart_g = 1
242 10248 : IF (iparticle2 == iparticle1) istart_g = igauss1
243 42212 : DO igauss2 = istart_g, SIZE(radii)
244 28488 : idim2 = (iparticle2 - 1)*SIZE(radii) + igauss2
245 59949180 : 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 59949180 : tmp = SUM(my_Amw)
249 28488 : Am(idim2, idim1) = tmp
250 38736 : Am(idim1, idim2) = tmp
251 : END DO
252 : END DO
253 : END DO
254 : END DO
255 246 : DEALLOCATE (gfunc_sq)
256 246 : DEALLOCATE (my_Amw)
257 246 : DEALLOCATE (my_Am)
258 : END IF
259 246 : CALL timestop(handle)
260 246 : 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 396 : 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 396 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: my_dbvw
293 396 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: my_dbv
294 : REAL(KIND=dp), DIMENSION(3) :: gvec, rvec
295 :
296 396 : CALL timeset(routineN, handle)
297 396 : gcut2 = gcut*gcut
298 396 : s_dim = rho_tot_g%pw_grid%first_gne0
299 396 : e_dim = rho_tot_g%pw_grid%ngpts_cut_local
300 396 : igmax = 0
301 350556 : DO ig = s_dim, e_dim
302 350556 : g2 = rho_tot_g%pw_grid%gsq(ig)
303 350556 : IF (g2 > gcut2) EXIT
304 350556 : igmax = ig
305 : END DO
306 396 : IF (igmax .GE. s_dim) THEN
307 1188 : ALLOCATE (my_dbv(3, s_dim:igmax))
308 1188 : ALLOCATE (my_dbvw(s_dim:igmax))
309 1812 : DO iparticle = 1, SIZE(particle_set)
310 1416 : IF (iparticle /= iparticle0) CYCLE
311 1584 : rvec = particle_set(iparticle)%r
312 350556 : DO ig = s_dim, igmax
313 1400640 : gvec = rho_tot_g%pw_grid%g(:, ig)
314 1400640 : arg = DOT_PRODUCT(gvec, rvec)
315 350160 : dphase = -CMPLX(SIN(arg), COS(arg), KIND=dp)
316 1401036 : my_dbv(:, ig) = w(ig)*REAL(CONJG(rho_tot_g%array(ig))*dphase, KIND=dp)*gvec(:)
317 : END DO
318 2892 : DO igauss = 1, SIZE(radii)
319 1080 : idim = (iparticle - 1)*SIZE(radii) + igauss
320 1042452 : DO ig = s_dim, igmax
321 1042452 : my_dbvw(ig) = my_dbv(1, ig)*gfunc(ig, igauss)
322 : END DO
323 1080 : dbv(idim, 1) = accurate_sum(my_dbvw)
324 1042452 : DO ig = s_dim, igmax
325 1042452 : my_dbvw(ig) = my_dbv(2, ig)*gfunc(ig, igauss)
326 : END DO
327 1080 : dbv(idim, 2) = accurate_sum(my_dbvw)
328 1042452 : DO ig = s_dim, igmax
329 1042452 : my_dbvw(ig) = my_dbv(3, ig)*gfunc(ig, igauss)
330 : END DO
331 1476 : dbv(idim, 3) = accurate_sum(my_dbvw)
332 : END DO
333 : END DO
334 396 : DEALLOCATE (my_dbvw)
335 396 : 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 396 : CALL timestop(handle)
346 396 : 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 420 : SUBROUTINE build_der_A_matrix_rows(dAm, gfunc, w, particle_set, radii, &
368 140 : 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 140 : REAL(KIND=dp), ALLOCATABLE :: lhs(:, :), rhs(:, :)
390 : INTEGER :: Nr, Np, Ng, icomp, ipp
391 :
392 140 : CALL timeset(routineN, handle)
393 140 : gcut2 = gcut*gcut
394 140 : s_dim = rho_tot_g%pw_grid%first_gne0
395 140 : e_dim = rho_tot_g%pw_grid%ngpts_cut_local
396 140 : igmax = 0
397 148940 : DO ig = s_dim, e_dim
398 148940 : g2 = rho_tot_g%pw_grid%gsq(ig)
399 148940 : IF (g2 > gcut2) EXIT
400 148940 : igmax = ig
401 : END DO
402 :
403 140 : Nr = SIZE(radii)
404 140 : Np = SIZE(particle_set)
405 140 : Ng = igmax - s_dim + 1
406 140 : IF (igmax .GE. s_dim) THEN
407 560 : ALLOCATE (lhs(nparticles*Nr, Ng))
408 560 : 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 536 : DO iparticle2 = 1, Np
413 1616 : DO igauss2 = 1, Nr
414 1042848 : 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 560 : DO icomp = 1, 3
418 : ! create lhs, which has all parts that depend on iparticle1
419 1608 : DO ipp = 1, nparticles
420 1188 : iparticle1 = iparticle0 + ipp - 1
421 1052088 : 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 4175784 : 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 420 : Ng, 0.0D0, dAm((iparticle0 - 1)*Nr + 1, 1, icomp), Np*Nr)
429 : ! do extra multiplies to compensate for missing factor of 2
430 1748 : DO ipp = 1, nparticles
431 1188 : 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 1608 : 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 536 : DO iparticle2 = 1, Np
441 1616 : DO igauss2 = 1, Nr
442 1042848 : 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 560 : DO icomp = 1, 3
446 : ! create lhs, which has all parts that depend on iparticle1
447 1608 : DO ipp = 1, nparticles
448 1188 : iparticle1 = iparticle0 + ipp - 1
449 1052088 : 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 4175784 : 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 420 : Ng, 1.0D0, dAm((iparticle0 - 1)*Nr + 1, 1, icomp), Np*Nr)
457 : ! do extra multiples to compensate for missing factor of 2
458 1748 : DO ipp = 1, nparticles
459 1188 : 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 1608 : Ng, 1.0D0, dAm((iparticle1 - 1)*Nr + 1, (iparticle1 - 1)*Nr + 1, icomp), Np*Nr)
462 : END DO
463 : END DO
464 :
465 140 : DEALLOCATE (rhs)
466 140 : DEALLOCATE (lhs)
467 : END IF
468 140 : CALL timestop(handle)
469 140 : 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 386 : 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 386 : IF (ALLOCATED(g_dot_rvec_sin)) DEALLOCATE (g_dot_rvec_sin)
480 386 : IF (ALLOCATED(g_dot_rvec_cos)) DEALLOCATE (g_dot_rvec_cos)
481 386 : 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 386 : 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 386 : gcut2 = gcut*gcut
501 386 : s_dim = rho_tot_g%pw_grid%first_gne0
502 386 : e_dim = rho_tot_g%pw_grid%ngpts_cut_local
503 386 : igmax = 0
504 323928 : DO ig = s_dim, e_dim
505 323928 : g2 = rho_tot_g%pw_grid%gsq(ig)
506 323928 : IF (g2 > gcut2) EXIT
507 323928 : igmax = ig
508 : END DO
509 :
510 386 : IF (igmax .GE. s_dim) THEN
511 1544 : ALLOCATE (g_dot_rvec_sin(1:igmax - s_dim + 1, SIZE(particle_set)))
512 1544 : ALLOCATE (g_dot_rvec_cos(1:igmax - s_dim + 1, SIZE(particle_set)))
513 :
514 1528 : DO iparticle = 1, SIZE(particle_set)
515 4568 : rvec = particle_set(iparticle)%r
516 1063322 : DO ig = s_dim, igmax
517 4247176 : g_dot_rvec = DOT_PRODUCT(rho_tot_g%pw_grid%g(:, ig), rvec)
518 1061794 : g_dot_rvec_sin(ig - s_dim + 1, iparticle) = SIN(g_dot_rvec)
519 1062936 : g_dot_rvec_cos(ig - s_dim + 1, iparticle) = COS(g_dot_rvec)
520 : END DO
521 : END DO
522 : END IF
523 :
524 386 : 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 246 : 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 246 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: AmE, cv
560 246 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: Am, AmI, Amw, g_dot_rvec_cos, &
561 246 : g_dot_rvec_sin
562 :
563 : !NB for precomputation of sin(g.r) and cos(g.r)
564 :
565 246 : CALL timeset(routineN, handle)
566 246 : ndim = SIZE(particle_set)*SIZE(radii)
567 984 : ALLOCATE (Am(ndim, ndim))
568 984 : ALLOCATE (AmI(ndim, ndim))
569 984 : ALLOCATE (GAmI(ndim, ndim))
570 738 : ALLOCATE (cv(ndim))
571 57222 : Am = 0.0_dp
572 57222 : AmI = 0.0_dp
573 2376 : 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 246 : CALL prep_g_dot_rvec_sin_cos(rho_tot_g, particle_set, gcut, g_dot_rvec_sin, g_dot_rvec_cos)
576 246 : CALL build_A_matrix(Am, gfunc, w, particle_set, radii, rho_tot_g, gcut, g_dot_rvec_sin, g_dot_rvec_cos)
577 246 : CALL cleanup_g_dot_rvec_sin_cos(g_dot_rvec_sin, g_dot_rvec_cos)
578 57222 : Am(:, :) = Am(:, :)/(Vol*Vol)
579 246 : CALL rho_tot_g%pw_grid%para%group%sum(Am)
580 246 : 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 246 : CALL invert_matrix(Am, AmI, inv_error, "N", improve=.FALSE.)
598 246 : IF (iw > 0) THEN
599 0 : WRITE (iw, '(T3,A,F15.9)') " Error inverting the A matrix: ", inv_error
600 : END IF
601 117066 : c0 = DOT_PRODUCT(cv, MATMUL(AmI, cv))
602 246 : DEALLOCATE (Am)
603 246 : DEALLOCATE (cv)
604 57222 : GAmI = AmI
605 246 : DEALLOCATE (AmI)
606 246 : CALL timestop(handle)
607 492 : 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 174 : 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, idim, idim1, idim2, idimo, igauss1, igauss2, ip1, ip2, &
639 : iparticle1, iparticle2, istart_g, k1, k2, k3, n_rep, ndim, nmax1, nmax2, nmax3, r1, r2, &
640 : r3, rmax1, rmax2, rmax3
641 : LOGICAL :: analyt
642 : REAL(KIND=dp) :: alpha, eps, ew_neut, fac, fac3, fs, g_ewald, galpha, gsq, gsqi, ij_fac, &
643 : my_val, r, r2tmp, r_ewald, rc1, rc12, rc2, rc22, rcut, rcut2, t1, tol, tol1
644 : REAL(KIND=dp), DIMENSION(3) :: fvec, gvec, ra, rvec
645 174 : REAL(KIND=dp), DIMENSION(:), POINTER :: EwM
646 :
647 174 : NULLIFY (EwM)
648 174 : CALL timeset(routineN, handle)
649 174 : CPASSERT(.NOT. ASSOCIATED(M))
650 174 : CPASSERT(ASSOCIATED(radii))
651 174 : CPASSERT(cell%orthorhombic)
652 174 : rcut = MIN(cell%hmat(1, 1), cell%hmat(2, 2), cell%hmat(3, 3))/2.0_dp
653 174 : CALL section_vals_val_get(multipole_section, "RCUT", n_rep_val=n_rep)
654 174 : IF (n_rep == 1) CALL section_vals_val_get(multipole_section, "RCUT", r_val=rcut)
655 174 : CALL section_vals_val_get(multipole_section, "EWALD_PRECISION", r_val=eps)
656 174 : CALL section_vals_val_get(multipole_section, "ANALYTICAL_GTERM", l_val=analyt)
657 174 : rcut2 = rcut**2
658 : !
659 : ! Setting-up parameters for Ewald summation
660 : !
661 174 : eps = MIN(ABS(eps), 0.5_dp)
662 174 : tol = SQRT(ABS(LOG(eps*rcut)))
663 174 : alpha = SQRT(ABS(LOG(eps*rcut*tol)))/rcut
664 174 : galpha = 1.0_dp/(4.0_dp*alpha*alpha)
665 174 : tol1 = SQRT(-LOG(eps*rcut*(2.0_dp*tol*alpha)**2))
666 174 : nmax1 = NINT(0.25_dp + cell%hmat(1, 1)*alpha*tol1/pi)
667 174 : nmax2 = NINT(0.25_dp + cell%hmat(2, 2)*alpha*tol1/pi)
668 174 : nmax3 = NINT(0.25_dp + cell%hmat(3, 3)*alpha*tol1/pi)
669 :
670 174 : rmax1 = CEILING(rcut/cell%hmat(1, 1))
671 174 : rmax2 = CEILING(rcut/cell%hmat(2, 2))
672 174 : rmax3 = CEILING(rcut/cell%hmat(3, 3))
673 174 : fac = 1.e0_dp/cell%deth
674 174 : fac3 = fac*pi
675 696 : fvec = twopi/(/cell%hmat(1, 1), cell%hmat(2, 2), cell%hmat(3, 3)/)
676 174 : ew_neut = -fac*pi/alpha**2
677 : !
678 174 : ewmdim = SIZE(particle_set)*(SIZE(particle_set) + 1)/2
679 174 : ndim = SIZE(particle_set)*SIZE(radii)
680 522 : ALLOCATE (EwM(ewmdim))
681 696 : ALLOCATE (M(ndim, ndim))
682 88338 : M = 0.0_dp
683 : !
684 5276 : idim = 0
685 5276 : EwM = 0.0_dp
686 786 : DO iparticle1 = 1, SIZE(particle_set)
687 5714 : ip1 = (iparticle1 - 1)*SIZE(radii)
688 5888 : DO iparticle2 = 1, iparticle1
689 5102 : ij_fac = 1.0_dp
690 5102 : IF (iparticle1 == iparticle2) ij_fac = 0.5_dp
691 :
692 5102 : ip2 = (iparticle2 - 1)*SIZE(radii)
693 5102 : idim = idim + 1
694 : !NB parallelization, done here so indexing is right
695 5102 : IF (MOD(iparticle1, cp_para_env%num_pe) /= cp_para_env%mepos) CYCLE
696 : !
697 : ! Real-Space Contribution
698 : !
699 2575 : my_val = 0.0_dp
700 10300 : rvec = particle_set(iparticle1)%r - particle_set(iparticle2)%r
701 2575 : r_ewald = 0.0_dp
702 2575 : IF (iparticle1 /= iparticle2) THEN
703 2253 : ra = rvec
704 9012 : r2tmp = DOT_PRODUCT(ra, ra)
705 2253 : IF (r2tmp <= rcut2) THEN
706 2165 : r = SQRT(r2tmp)
707 2165 : t1 = erfc(alpha*r)/r
708 2165 : r_ewald = t1
709 : END IF
710 : END IF
711 14670 : DO r1 = -rmax1, rmax1
712 73181 : DO r2 = -rmax2, rmax2
713 360381 : DO r3 = -rmax3, rmax3
714 289775 : IF ((r1 == 0) .AND. (r2 == 0) .AND. (r3 == 0)) CYCLE
715 287200 : ra(1) = rvec(1) + cell%hmat(1, 1)*r1
716 287200 : ra(2) = rvec(2) + cell%hmat(2, 2)*r2
717 287200 : ra(3) = rvec(3) + cell%hmat(3, 3)*r3
718 1148800 : r2tmp = DOT_PRODUCT(ra, ra)
719 345711 : IF (r2tmp <= rcut2) THEN
720 47717 : r = SQRT(r2tmp)
721 47717 : t1 = erfc(alpha*r)/r
722 47717 : r_ewald = r_ewald + t1*ij_fac
723 : END IF
724 : END DO
725 : END DO
726 : END DO
727 : !
728 : ! G-space Contribution
729 : !
730 2575 : IF (analyt) THEN
731 : g_ewald = 0.0_dp
732 2917 : DO k1 = 0, nmax1
733 74424 : DO k2 = -nmax2, nmax2
734 2157813 : DO k3 = -nmax3, nmax3
735 2083619 : IF (k1 == 0 .AND. k2 == 0 .AND. k3 == 0) CYCLE
736 2083389 : fs = 2.0_dp; IF (k1 == 0) fs = 1.0_dp
737 8333556 : gvec = fvec*(/REAL(k1, KIND=dp), REAL(k2, KIND=dp), REAL(k3, KIND=dp)/)
738 8333556 : gsq = DOT_PRODUCT(gvec, gvec)
739 2083389 : gsqi = fs/gsq
740 2083389 : t1 = fac*gsqi*EXP(-galpha*gsq)
741 8405293 : g_ewald = g_ewald + t1*COS(DOT_PRODUCT(gvec, rvec))
742 : END DO
743 : END DO
744 : END DO
745 : ELSE
746 2345 : g_ewald = Eval_Interp_Spl3_pbc(rvec, coeff)
747 : END IF
748 : !
749 : ! G-EWALD, R-EWALD
750 : !
751 2575 : g_ewald = r_ewald + fourpi*g_ewald
752 : !
753 : ! Self Contribution
754 : !
755 2575 : IF (iparticle1 == iparticle2) THEN
756 322 : g_ewald = g_ewald - 2.0_dp*alpha*oorootpi
757 : END IF
758 : !
759 2575 : IF (iparticle1 /= iparticle2) THEN
760 2253 : ra = rvec
761 9012 : r = SQRT(DOT_PRODUCT(ra, ra))
762 2253 : my_val = factor/r
763 : END IF
764 5714 : EwM(idim) = my_val - factor*g_ewald
765 : END DO ! iparticle2
766 : END DO ! iparticle1
767 : !NB sum over parallelized contributions of different nodes
768 10378 : CALL cp_para_env%sum(EwM)
769 174 : idim = 0
770 786 : DO iparticle2 = 1, SIZE(particle_set)
771 612 : ip2 = (iparticle2 - 1)*SIZE(radii)
772 612 : idimo = (iparticle2 - 1)
773 612 : idimo = idimo*(idimo + 1)/2
774 2622 : DO igauss2 = 1, SIZE(radii)
775 1836 : idim2 = ip2 + igauss2
776 1836 : rc2 = radii(igauss2)
777 1836 : rc22 = rc2*rc2
778 17754 : DO iparticle1 = 1, iparticle2
779 15306 : ip1 = (iparticle1 - 1)*SIZE(radii)
780 15306 : idim = idimo + iparticle1
781 15306 : istart_g = 1
782 15306 : IF (iparticle1 == iparticle2) istart_g = igauss2
783 61224 : DO igauss1 = istart_g, SIZE(radii)
784 44082 : idim1 = ip1 + igauss1
785 44082 : rc1 = radii(igauss1)
786 44082 : rc12 = rc1*rc1
787 44082 : M(idim1, idim2) = EwM(idim) - factor*ew_neut - factor*fac3*(rc12 + rc22)
788 59388 : M(idim2, idim1) = M(idim1, idim2)
789 : END DO
790 : END DO
791 : END DO ! iparticle2
792 : END DO ! iparticle1
793 174 : DEALLOCATE (EwM)
794 174 : CALL timestop(handle)
795 522 : END SUBROUTINE ewald_ddapc_pot
796 :
797 : ! **************************************************************************************************
798 : !> \brief Evaluates the electrostatic potential due to a simple solvation model
799 : !> Spherical cavity in a dieletric medium
800 : !> \param solvation_section ...
801 : !> \param particle_set ...
802 : !> \param M ...
803 : !> \param radii ...
804 : !> \par History
805 : !> 08.2006 created [tlaino]
806 : !> \author Teodoro Laino
807 : ! **************************************************************************************************
808 26 : SUBROUTINE solvation_ddapc_pot(solvation_section, particle_set, M, radii)
809 : TYPE(section_vals_type), POINTER :: solvation_section
810 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
811 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: M
812 : REAL(KIND=dp), DIMENSION(:), POINTER :: radii
813 :
814 : INTEGER :: i, idim, idim1, idim2, igauss1, igauss2, ip1, ip2, iparticle1, iparticle2, &
815 : istart_g, j, l, lmax, n_rep1, n_rep2, ndim, output_unit, weight
816 26 : INTEGER, DIMENSION(:), POINTER :: list
817 : LOGICAL :: fixed_center
818 : REAL(KIND=dp) :: center(3), eps_in, eps_out, factor, &
819 : mass, mycos, r1, r2, Rs, rvec(3)
820 26 : REAL(KIND=dp), DIMENSION(:), POINTER :: pos, R0
821 26 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: cost, LocP
822 :
823 26 : fixed_center = .FALSE.
824 52 : output_unit = cp_logger_get_default_io_unit()
825 26 : ndim = SIZE(particle_set)*SIZE(radii)
826 104 : ALLOCATE (M(ndim, ndim))
827 1274 : M = 0.0_dp
828 26 : eps_in = 1.0_dp
829 26 : CALL section_vals_val_get(solvation_section, "EPS_OUT", r_val=eps_out)
830 26 : CALL section_vals_val_get(solvation_section, "LMAX", i_val=lmax)
831 26 : CALL section_vals_val_get(solvation_section, "SPHERE%RADIUS", r_val=Rs)
832 26 : CALL section_vals_val_get(solvation_section, "SPHERE%CENTER%XYZ", n_rep_val=n_rep1)
833 26 : IF (n_rep1 /= 0) THEN
834 24 : CALL section_vals_val_get(solvation_section, "SPHERE%CENTER%XYZ", r_vals=R0)
835 96 : center = R0
836 : ELSE
837 : CALL section_vals_val_get(solvation_section, "SPHERE%CENTER%ATOM_LIST", &
838 2 : n_rep_val=n_rep2)
839 2 : IF (n_rep2 /= 0) THEN
840 2 : CALL section_vals_val_get(solvation_section, "SPHERE%CENTER%ATOM_LIST", i_vals=list)
841 2 : CALL section_vals_val_get(solvation_section, "SPHERE%CENTER%WEIGHT_TYPE", i_val=weight)
842 2 : ALLOCATE (R0(3))
843 : SELECT CASE (weight)
844 : CASE (weight_type_unit)
845 8 : R0 = 0.0_dp
846 4 : DO i = 1, SIZE(list)
847 10 : R0 = R0 + particle_set(list(i))%r
848 : END DO
849 8 : R0 = R0/REAL(SIZE(list), KIND=dp)
850 : CASE (weight_type_mass)
851 0 : R0 = 0.0_dp
852 0 : mass = 0.0_dp
853 0 : DO i = 1, SIZE(list)
854 0 : R0 = R0 + particle_set(list(i))%r*particle_set(list(i))%atomic_kind%mass
855 0 : mass = mass + particle_set(list(i))%atomic_kind%mass
856 : END DO
857 2 : R0 = R0/mass
858 : END SELECT
859 8 : center = R0
860 2 : CALL section_vals_val_get(solvation_section, "SPHERE%CENTER%FIXED", l_val=fixed_center)
861 4 : IF (fixed_center) THEN
862 : CALL section_vals_val_set(solvation_section, "SPHERE%CENTER%XYZ", &
863 2 : r_vals_ptr=R0)
864 : ELSE
865 0 : DEALLOCATE (R0)
866 : END IF
867 : END IF
868 : END IF
869 26 : CPASSERT(n_rep1 /= 0 .OR. n_rep2 /= 0)
870 : ! Potential calculation
871 104 : ALLOCATE (LocP(0:lmax, SIZE(particle_set)))
872 78 : ALLOCATE (pos(SIZE(particle_set)))
873 104 : ALLOCATE (cost(SIZE(particle_set), SIZE(particle_set)))
874 : ! Determining the single atomic contribution to the dielectric dipole
875 76 : DO i = 1, SIZE(particle_set)
876 200 : rvec = particle_set(i)%r - center
877 200 : r2 = DOT_PRODUCT(rvec, rvec)
878 50 : r1 = SQRT(r2)
879 50 : IF (r1 >= Rs) THEN
880 0 : IF (output_unit > 0) THEN
881 0 : WRITE (output_unit, '(A,I6,A)') "Atom number :: ", i, " is out of the solvation sphere"
882 0 : WRITE (output_unit, '(2(A,F12.6))') "Distance from the center::", r1, " Radius of the sphere::", rs
883 : END IF
884 0 : CPABORT("")
885 : END IF
886 250 : LocP(:, i) = 0.0_dp
887 50 : IF (r1 /= 0.0_dp) THEN
888 230 : DO l = 0, lmax
889 : LocP(l, i) = (r1**l*REAL(l + 1, KIND=dp)*(eps_in - eps_out))/ &
890 230 : (Rs**(2*l + 1)*eps_in*(REAL(l, KIND=dp)*eps_in + REAL(l + 1, KIND=dp)*eps_out))
891 : END DO
892 : ELSE
893 : ! limit for r->0
894 4 : LocP(0, i) = (eps_in - eps_out)/(Rs*eps_in*eps_out)
895 : END IF
896 76 : pos(i) = r1
897 : END DO
898 : ! Particle-Particle potential energy matrix
899 198 : cost = 0.0_dp
900 76 : DO i = 1, SIZE(particle_set)
901 162 : DO j = 1, i
902 86 : factor = 0.0_dp
903 86 : IF (pos(i)*pos(j) /= 0.0_dp) THEN
904 296 : mycos = DOT_PRODUCT(particle_set(i)%r - center, particle_set(j)%r - center)/(pos(i)*pos(j))
905 74 : IF (ABS(mycos) > 1.0_dp) mycos = SIGN(1.0_dp, mycos)
906 370 : DO l = 0, lmax
907 370 : factor = factor + LocP(l, i)*pos(j)**l*legendre(mycos, l, 0)
908 : END DO
909 : ELSE
910 12 : factor = LocP(0, i)
911 : END IF
912 86 : cost(i, j) = factor
913 136 : cost(j, i) = factor
914 : END DO
915 : END DO
916 : ! Computes the full potential energy matrix
917 26 : idim = 0
918 76 : DO iparticle2 = 1, SIZE(particle_set)
919 50 : ip2 = (iparticle2 - 1)*SIZE(radii)
920 226 : DO igauss2 = 1, SIZE(radii)
921 150 : idim2 = ip2 + igauss2
922 458 : DO iparticle1 = 1, iparticle2
923 258 : ip1 = (iparticle1 - 1)*SIZE(radii)
924 258 : istart_g = 1
925 258 : IF (iparticle1 == iparticle2) istart_g = igauss2
926 1032 : DO igauss1 = istart_g, SIZE(radii)
927 624 : idim1 = ip1 + igauss1
928 624 : M(idim1, idim2) = cost(iparticle1, iparticle2)
929 882 : M(idim2, idim1) = M(idim1, idim2)
930 : END DO
931 : END DO
932 : END DO
933 : END DO
934 26 : DEALLOCATE (cost)
935 26 : DEALLOCATE (pos)
936 26 : DEALLOCATE (LocP)
937 52 : END SUBROUTINE solvation_ddapc_pot
938 :
939 246 : END MODULE cp_ddapc_methods
|
