(git:374b731)
Loading...
Searching...
No Matches
cp_ddapc_methods.F
Go to the documentation of this file.
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! **************************************************************************************************
12MODULE 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
51CONTAINS
52
53! **************************************************************************************************
54!> \brief ...
55!> \param gfunc ...
56!> \param w ...
57!> \param gcut ...
58!> \param rho_tot_g ...
59!> \param radii ...
60! **************************************************************************************************
61 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 CALL timeset(routinen, handle)
74 gcut2 = gcut*gcut
75 !
76 s_dim = rho_tot_g%pw_grid%first_gne0
77 e_dim = rho_tot_g%pw_grid%ngpts_cut_local
78 ALLOCATE (gfunc(s_dim:e_dim, SIZE(radii)))
79 ALLOCATE (w(s_dim:e_dim))
80 gfunc = 0.0_dp
81 w = 0.0_dp
82 DO igauss = 1, SIZE(radii)
83 rc = radii(igauss)
84 rc2 = rc*rc
85 DO ig = s_dim, e_dim
86 g2 = rho_tot_g%pw_grid%gsq(ig)
87 IF (g2 > gcut2) EXIT
88 gfunc(ig, igauss) = exp(-g2*rc2/4.0_dp)
89 END DO
90 END DO
91 DO ig = s_dim, e_dim
92 g2 = rho_tot_g%pw_grid%gsq(ig)
93 IF (g2 > gcut2) EXIT
94 w(ig) = fourpi*(g2 - gcut2)**2/(g2*gcut2)
95 END DO
96 CALL timestop(handle)
97 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 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 REAL(kind=dp), DIMENSION(:), POINTER :: my_bv, my_bvw
128
129 CALL timeset(routinen, handle)
130 NULLIFY (my_bv, my_bvw)
131 gcut2 = gcut*gcut
132 s_dim = rho_tot_g%pw_grid%first_gne0
133 e_dim = rho_tot_g%pw_grid%ngpts_cut_local
134 igmax = 0
135 DO ig = s_dim, e_dim
136 g2 = rho_tot_g%pw_grid%gsq(ig)
137 IF (g2 > gcut2) EXIT
138 igmax = ig
139 END DO
140 IF (igmax .GE. s_dim) THEN
141 ALLOCATE (my_bv(s_dim:igmax))
142 ALLOCATE (my_bvw(s_dim:igmax))
143 !
144 DO iparticle = 1, SIZE(particle_set)
145 rvec = particle_set(iparticle)%r
146 my_bv = 0.0_dp
147 DO ig = s_dim, igmax
148 gvec = rho_tot_g%pw_grid%g(:, ig)
149 arg = dot_product(gvec, rvec)
150 phase = cmplx(cos(arg), -sin(arg), kind=dp)
151 my_bv(ig) = w(ig)*real(conjg(rho_tot_g%array(ig))*phase, kind=dp)
152 END DO
153 DO igauss = 1, SIZE(radii)
154 idim = (iparticle - 1)*SIZE(radii) + igauss
155 DO ig = s_dim, igmax
156 my_bvw(ig) = my_bv(ig)*gfunc(ig, igauss)
157 END DO
158 bv(idim) = accurate_sum(my_bvw)
159 END DO
160 END DO
161 DEALLOCATE (my_bvw)
162 DEALLOCATE (my_bv)
163 ELSE
164 DO iparticle = 1, SIZE(particle_set)
165 DO igauss = 1, SIZE(radii)
166 idim = (iparticle - 1)*SIZE(radii) + igauss
167 bv(idim) = 0.0_dp
168 END DO
169 END DO
170 END IF
171 CALL timestop(handle)
172 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 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: my_am, my_amw
207 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 CALL timeset(routinen, handle)
212 gcut2 = gcut*gcut
213 s_dim = rho_tot_g%pw_grid%first_gne0
214 e_dim = rho_tot_g%pw_grid%ngpts_cut_local
215 igmax = 0
216 DO ig = s_dim, e_dim
217 g2 = rho_tot_g%pw_grid%gsq(ig)
218 IF (g2 > gcut2) EXIT
219 igmax = ig
220 END DO
221 IF (igmax .GE. s_dim) THEN
222 ALLOCATE (my_am(s_dim:igmax))
223 ALLOCATE (my_amw(s_dim:igmax))
224 ALLOCATE (gfunc_sq(s_dim:igmax, SIZE(radii), SIZE(radii)))
225
226 DO igauss1 = 1, SIZE(radii)
227 DO igauss2 = 1, SIZE(radii)
228 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 DO iparticle1 = 1, SIZE(particle_set)
233 DO iparticle2 = iparticle1, SIZE(particle_set)
234 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 g_dot_rvec_sin(ig - s_dim + 1, iparticle1)*g_dot_rvec_sin(ig - s_dim + 1, iparticle2))
238 END DO
239 DO igauss1 = 1, SIZE(radii)
240 idim1 = (iparticle1 - 1)*SIZE(radii) + igauss1
241 istart_g = 1
242 IF (iparticle2 == iparticle1) istart_g = igauss1
243 DO igauss2 = istart_g, SIZE(radii)
244 idim2 = (iparticle2 - 1)*SIZE(radii) + igauss2
245 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 tmp = sum(my_amw)
249 am(idim2, idim1) = tmp
250 am(idim1, idim2) = tmp
251 END DO
252 END DO
253 END DO
254 END DO
255 DEALLOCATE (gfunc_sq)
256 DEALLOCATE (my_amw)
257 DEALLOCATE (my_am)
258 END IF
259 CALL timestop(handle)
260 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 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: my_dbvw
293 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: my_dbv
294 REAL(kind=dp), DIMENSION(3) :: gvec, rvec
295
296 CALL timeset(routinen, handle)
297 gcut2 = gcut*gcut
298 s_dim = rho_tot_g%pw_grid%first_gne0
299 e_dim = rho_tot_g%pw_grid%ngpts_cut_local
300 igmax = 0
301 DO ig = s_dim, e_dim
302 g2 = rho_tot_g%pw_grid%gsq(ig)
303 IF (g2 > gcut2) EXIT
304 igmax = ig
305 END DO
306 IF (igmax .GE. s_dim) THEN
307 ALLOCATE (my_dbv(3, s_dim:igmax))
308 ALLOCATE (my_dbvw(s_dim:igmax))
309 DO iparticle = 1, SIZE(particle_set)
310 IF (iparticle /= iparticle0) cycle
311 rvec = particle_set(iparticle)%r
312 DO ig = s_dim, igmax
313 gvec = rho_tot_g%pw_grid%g(:, ig)
314 arg = dot_product(gvec, rvec)
315 dphase = -cmplx(sin(arg), cos(arg), kind=dp)
316 my_dbv(:, ig) = w(ig)*real(conjg(rho_tot_g%array(ig))*dphase, kind=dp)*gvec(:)
317 END DO
318 DO igauss = 1, SIZE(radii)
319 idim = (iparticle - 1)*SIZE(radii) + igauss
320 DO ig = s_dim, igmax
321 my_dbvw(ig) = my_dbv(1, ig)*gfunc(ig, igauss)
322 END DO
323 dbv(idim, 1) = accurate_sum(my_dbvw)
324 DO ig = s_dim, igmax
325 my_dbvw(ig) = my_dbv(2, ig)*gfunc(ig, igauss)
326 END DO
327 dbv(idim, 2) = accurate_sum(my_dbvw)
328 DO ig = s_dim, igmax
329 my_dbvw(ig) = my_dbv(3, ig)*gfunc(ig, igauss)
330 END DO
331 dbv(idim, 3) = accurate_sum(my_dbvw)
332 END DO
333 END DO
334 DEALLOCATE (my_dbvw)
335 DEALLOCATE (my_dbv)
336 ELSE
337 DO iparticle = 1, SIZE(particle_set)
338 IF (iparticle /= iparticle0) cycle
339 DO igauss = 1, SIZE(radii)
340 idim = (iparticle - 1)*SIZE(radii) + igauss
341 dbv(idim, 1:3) = 0.0_dp
342 END DO
343 END DO
344 END IF
345 CALL timestop(handle)
346 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 SUBROUTINE build_der_a_matrix_rows(dAm, gfunc, w, particle_set, radii, &
368 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 REAL(kind=dp), ALLOCATABLE :: lhs(:, :), rhs(:, :)
390 INTEGER :: nr, np, ng, icomp, ipp
391
392 CALL timeset(routinen, handle)
393 gcut2 = gcut*gcut
394 s_dim = rho_tot_g%pw_grid%first_gne0
395 e_dim = rho_tot_g%pw_grid%ngpts_cut_local
396 igmax = 0
397 DO ig = s_dim, e_dim
398 g2 = rho_tot_g%pw_grid%gsq(ig)
399 IF (g2 > gcut2) EXIT
400 igmax = ig
401 END DO
402
403 nr = SIZE(radii)
404 np = SIZE(particle_set)
405 ng = igmax - s_dim + 1
406 IF (igmax .GE. s_dim) THEN
407 ALLOCATE (lhs(nparticles*nr, ng))
408 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 DO iparticle2 = 1, np
413 DO igauss2 = 1, nr
414 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 DO icomp = 1, 3
418 ! create lhs, which has all parts that depend on iparticle1
419 DO ipp = 1, nparticles
420 iparticle1 = iparticle0 + ipp - 1
421 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 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 ng, 0.0d0, dam((iparticle0 - 1)*nr + 1, 1, icomp), np*nr)
429 ! do extra multiplies to compensate for missing factor of 2
430 DO ipp = 1, nparticles
431 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 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 DO iparticle2 = 1, np
441 DO igauss2 = 1, nr
442 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 DO icomp = 1, 3
446 ! create lhs, which has all parts that depend on iparticle1
447 DO ipp = 1, nparticles
448 iparticle1 = iparticle0 + ipp - 1
449 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 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 ng, 1.0d0, dam((iparticle0 - 1)*nr + 1, 1, icomp), np*nr)
457 ! do extra multiples to compensate for missing factor of 2
458 DO ipp = 1, nparticles
459 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 ng, 1.0d0, dam((iparticle1 - 1)*nr + 1, (iparticle1 - 1)*nr + 1, icomp), np*nr)
462 END DO
463 END DO
464
465 DEALLOCATE (rhs)
466 DEALLOCATE (lhs)
467 END IF
468 CALL timestop(handle)
469 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 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 IF (ALLOCATED(g_dot_rvec_sin)) DEALLOCATE (g_dot_rvec_sin)
480 IF (ALLOCATED(g_dot_rvec_cos)) DEALLOCATE (g_dot_rvec_cos)
481 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 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 gcut2 = gcut*gcut
501 s_dim = rho_tot_g%pw_grid%first_gne0
502 e_dim = rho_tot_g%pw_grid%ngpts_cut_local
503 igmax = 0
504 DO ig = s_dim, e_dim
505 g2 = rho_tot_g%pw_grid%gsq(ig)
506 IF (g2 > gcut2) EXIT
507 igmax = ig
508 END DO
509
510 IF (igmax .GE. s_dim) THEN
511 ALLOCATE (g_dot_rvec_sin(1:igmax - s_dim + 1, SIZE(particle_set)))
512 ALLOCATE (g_dot_rvec_cos(1:igmax - s_dim + 1, SIZE(particle_set)))
513
514 DO iparticle = 1, SIZE(particle_set)
515 rvec = particle_set(iparticle)%r
516 DO ig = s_dim, igmax
517 g_dot_rvec = dot_product(rho_tot_g%pw_grid%g(:, ig), rvec)
518 g_dot_rvec_sin(ig - s_dim + 1, iparticle) = sin(g_dot_rvec)
519 g_dot_rvec_cos(ig - s_dim + 1, iparticle) = cos(g_dot_rvec)
520 END DO
521 END DO
522 END IF
523
524 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 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: ame, cv
560 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: am, ami, amw, g_dot_rvec_cos, &
561 g_dot_rvec_sin
562
563!NB for precomputation of sin(g.r) and cos(g.r)
564
565 CALL timeset(routinen, handle)
566 ndim = SIZE(particle_set)*SIZE(radii)
567 ALLOCATE (am(ndim, ndim))
568 ALLOCATE (ami(ndim, ndim))
569 ALLOCATE (gami(ndim, ndim))
570 ALLOCATE (cv(ndim))
571 am = 0.0_dp
572 ami = 0.0_dp
573 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 CALL prep_g_dot_rvec_sin_cos(rho_tot_g, particle_set, gcut, g_dot_rvec_sin, g_dot_rvec_cos)
576 CALL build_a_matrix(am, gfunc, w, particle_set, radii, rho_tot_g, gcut, g_dot_rvec_sin, g_dot_rvec_cos)
577 CALL cleanup_g_dot_rvec_sin_cos(g_dot_rvec_sin, g_dot_rvec_cos)
578 am(:, :) = am(:, :)/(vol*vol)
579 CALL rho_tot_g%pw_grid%para%group%sum(am)
580 IF (iw > 0) THEN
581 ! Checking conditions numbers and eigenvalues
582 ALLOCATE (amw(ndim, ndim))
583 ALLOCATE (ame(ndim))
584 amw(:, :) = am
585 CALL diamat_all(amw, ame)
586 condition_number = maxval(abs(ame))/minval(abs(ame))
587 WRITE (iw, '(T3,A)') " Eigenvalues of Matrix A:"
588 WRITE (iw, '(T3,4E15.8)') ame
589 WRITE (iw, '(T3,A,1E15.9)') " Condition number:", condition_number
590 IF (condition_number > 1.0e12_dp) THEN
591 WRITE (iw, fmt="(/,T2,A)") &
592 "WARNING: high condition number => possibly ill-conditioned matrix"
593 END IF
594 DEALLOCATE (amw)
595 DEALLOCATE (ame)
596 END IF
597 CALL invert_matrix(am, ami, inv_error, "N", improve=.false.)
598 IF (iw > 0) THEN
599 WRITE (iw, '(T3,A,F15.9)') " Error inverting the A matrix: ", inv_error
600 END IF
601 c0 = dot_product(cv, matmul(ami, cv))
602 DEALLOCATE (am)
603 DEALLOCATE (cv)
604 gami = ami
605 DEALLOCATE (ami)
606 CALL timestop(handle)
607 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 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 REAL(kind=dp), DIMENSION(:), POINTER :: ewm
646
647 NULLIFY (ewm)
648 CALL timeset(routinen, handle)
649 cpassert(.NOT. ASSOCIATED(m))
650 cpassert(ASSOCIATED(radii))
651 cpassert(cell%orthorhombic)
652 rcut = min(cell%hmat(1, 1), cell%hmat(2, 2), cell%hmat(3, 3))/2.0_dp
653 CALL section_vals_val_get(multipole_section, "RCUT", n_rep_val=n_rep)
654 IF (n_rep == 1) CALL section_vals_val_get(multipole_section, "RCUT", r_val=rcut)
655 CALL section_vals_val_get(multipole_section, "EWALD_PRECISION", r_val=eps)
656 CALL section_vals_val_get(multipole_section, "ANALYTICAL_GTERM", l_val=analyt)
657 rcut2 = rcut**2
658 !
659 ! Setting-up parameters for Ewald summation
660 !
661 eps = min(abs(eps), 0.5_dp)
662 tol = sqrt(abs(log(eps*rcut)))
663 alpha = sqrt(abs(log(eps*rcut*tol)))/rcut
664 galpha = 1.0_dp/(4.0_dp*alpha*alpha)
665 tol1 = sqrt(-log(eps*rcut*(2.0_dp*tol*alpha)**2))
666 nmax1 = nint(0.25_dp + cell%hmat(1, 1)*alpha*tol1/pi)
667 nmax2 = nint(0.25_dp + cell%hmat(2, 2)*alpha*tol1/pi)
668 nmax3 = nint(0.25_dp + cell%hmat(3, 3)*alpha*tol1/pi)
669
670 rmax1 = ceiling(rcut/cell%hmat(1, 1))
671 rmax2 = ceiling(rcut/cell%hmat(2, 2))
672 rmax3 = ceiling(rcut/cell%hmat(3, 3))
673 fac = 1.e0_dp/cell%deth
674 fac3 = fac*pi
675 fvec = twopi/(/cell%hmat(1, 1), cell%hmat(2, 2), cell%hmat(3, 3)/)
676 ew_neut = -fac*pi/alpha**2
677 !
678 ewmdim = SIZE(particle_set)*(SIZE(particle_set) + 1)/2
679 ndim = SIZE(particle_set)*SIZE(radii)
680 ALLOCATE (ewm(ewmdim))
681 ALLOCATE (m(ndim, ndim))
682 m = 0.0_dp
683 !
684 idim = 0
685 ewm = 0.0_dp
686 DO iparticle1 = 1, SIZE(particle_set)
687 ip1 = (iparticle1 - 1)*SIZE(radii)
688 DO iparticle2 = 1, iparticle1
689 ij_fac = 1.0_dp
690 IF (iparticle1 == iparticle2) ij_fac = 0.5_dp
691
692 ip2 = (iparticle2 - 1)*SIZE(radii)
693 idim = idim + 1
694 !NB parallelization, done here so indexing is right
695 IF (mod(iparticle1, cp_para_env%num_pe) /= cp_para_env%mepos) cycle
696 !
697 ! Real-Space Contribution
698 !
699 my_val = 0.0_dp
700 rvec = particle_set(iparticle1)%r - particle_set(iparticle2)%r
701 r_ewald = 0.0_dp
702 IF (iparticle1 /= iparticle2) THEN
703 ra = rvec
704 r2tmp = dot_product(ra, ra)
705 IF (r2tmp <= rcut2) THEN
706 r = sqrt(r2tmp)
707 t1 = erfc(alpha*r)/r
708 r_ewald = t1
709 END IF
710 END IF
711 DO r1 = -rmax1, rmax1
712 DO r2 = -rmax2, rmax2
713 DO r3 = -rmax3, rmax3
714 IF ((r1 == 0) .AND. (r2 == 0) .AND. (r3 == 0)) cycle
715 ra(1) = rvec(1) + cell%hmat(1, 1)*r1
716 ra(2) = rvec(2) + cell%hmat(2, 2)*r2
717 ra(3) = rvec(3) + cell%hmat(3, 3)*r3
718 r2tmp = dot_product(ra, ra)
719 IF (r2tmp <= rcut2) THEN
720 r = sqrt(r2tmp)
721 t1 = erfc(alpha*r)/r
722 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 IF (analyt) THEN
731 g_ewald = 0.0_dp
732 DO k1 = 0, nmax1
733 DO k2 = -nmax2, nmax2
734 DO k3 = -nmax3, nmax3
735 IF (k1 == 0 .AND. k2 == 0 .AND. k3 == 0) cycle
736 fs = 2.0_dp; IF (k1 == 0) fs = 1.0_dp
737 gvec = fvec*(/real(k1, kind=dp), real(k2, kind=dp), real(k3, kind=dp)/)
738 gsq = dot_product(gvec, gvec)
739 gsqi = fs/gsq
740 t1 = fac*gsqi*exp(-galpha*gsq)
741 g_ewald = g_ewald + t1*cos(dot_product(gvec, rvec))
742 END DO
743 END DO
744 END DO
745 ELSE
746 g_ewald = eval_interp_spl3_pbc(rvec, coeff)
747 END IF
748 !
749 ! G-EWALD, R-EWALD
750 !
751 g_ewald = r_ewald + fourpi*g_ewald
752 !
753 ! Self Contribution
754 !
755 IF (iparticle1 == iparticle2) THEN
756 g_ewald = g_ewald - 2.0_dp*alpha*oorootpi
757 END IF
758 !
759 IF (iparticle1 /= iparticle2) THEN
760 ra = rvec
761 r = sqrt(dot_product(ra, ra))
762 my_val = factor/r
763 END IF
764 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 CALL cp_para_env%sum(ewm)
769 idim = 0
770 DO iparticle2 = 1, SIZE(particle_set)
771 ip2 = (iparticle2 - 1)*SIZE(radii)
772 idimo = (iparticle2 - 1)
773 idimo = idimo*(idimo + 1)/2
774 DO igauss2 = 1, SIZE(radii)
775 idim2 = ip2 + igauss2
776 rc2 = radii(igauss2)
777 rc22 = rc2*rc2
778 DO iparticle1 = 1, iparticle2
779 ip1 = (iparticle1 - 1)*SIZE(radii)
780 idim = idimo + iparticle1
781 istart_g = 1
782 IF (iparticle1 == iparticle2) istart_g = igauss2
783 DO igauss1 = istart_g, SIZE(radii)
784 idim1 = ip1 + igauss1
785 rc1 = radii(igauss1)
786 rc12 = rc1*rc1
787 m(idim1, idim2) = ewm(idim) - factor*ew_neut - factor*fac3*(rc12 + rc22)
788 m(idim2, idim1) = m(idim1, idim2)
789 END DO
790 END DO
791 END DO ! iparticle2
792 END DO ! iparticle1
793 DEALLOCATE (ewm)
794 CALL timestop(handle)
795 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 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 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 REAL(kind=dp), DIMENSION(:), POINTER :: pos, r0
821 REAL(kind=dp), DIMENSION(:, :), POINTER :: cost, locp
822
823 fixed_center = .false.
824 output_unit = cp_logger_get_default_io_unit()
825 ndim = SIZE(particle_set)*SIZE(radii)
826 ALLOCATE (m(ndim, ndim))
827 m = 0.0_dp
828 eps_in = 1.0_dp
829 CALL section_vals_val_get(solvation_section, "EPS_OUT", r_val=eps_out)
830 CALL section_vals_val_get(solvation_section, "LMAX", i_val=lmax)
831 CALL section_vals_val_get(solvation_section, "SPHERE%RADIUS", r_val=rs)
832 CALL section_vals_val_get(solvation_section, "SPHERE%CENTER%XYZ", n_rep_val=n_rep1)
833 IF (n_rep1 /= 0) THEN
834 CALL section_vals_val_get(solvation_section, "SPHERE%CENTER%XYZ", r_vals=r0)
835 center = r0
836 ELSE
837 CALL section_vals_val_get(solvation_section, "SPHERE%CENTER%ATOM_LIST", &
838 n_rep_val=n_rep2)
839 IF (n_rep2 /= 0) THEN
840 CALL section_vals_val_get(solvation_section, "SPHERE%CENTER%ATOM_LIST", i_vals=list)
841 CALL section_vals_val_get(solvation_section, "SPHERE%CENTER%WEIGHT_TYPE", i_val=weight)
842 ALLOCATE (r0(3))
843 SELECT CASE (weight)
844 CASE (weight_type_unit)
845 r0 = 0.0_dp
846 DO i = 1, SIZE(list)
847 r0 = r0 + particle_set(list(i))%r
848 END DO
849 r0 = r0/real(SIZE(list), kind=dp)
850 CASE (weight_type_mass)
851 r0 = 0.0_dp
852 mass = 0.0_dp
853 DO i = 1, SIZE(list)
854 r0 = r0 + particle_set(list(i))%r*particle_set(list(i))%atomic_kind%mass
855 mass = mass + particle_set(list(i))%atomic_kind%mass
856 END DO
857 r0 = r0/mass
858 END SELECT
859 center = r0
860 CALL section_vals_val_get(solvation_section, "SPHERE%CENTER%FIXED", l_val=fixed_center)
861 IF (fixed_center) THEN
862 CALL section_vals_val_set(solvation_section, "SPHERE%CENTER%XYZ", &
863 r_vals_ptr=r0)
864 ELSE
865 DEALLOCATE (r0)
866 END IF
867 END IF
868 END IF
869 cpassert(n_rep1 /= 0 .OR. n_rep2 /= 0)
870 ! Potential calculation
871 ALLOCATE (locp(0:lmax, SIZE(particle_set)))
872 ALLOCATE (pos(SIZE(particle_set)))
873 ALLOCATE (cost(SIZE(particle_set), SIZE(particle_set)))
874 ! Determining the single atomic contribution to the dielectric dipole
875 DO i = 1, SIZE(particle_set)
876 rvec = particle_set(i)%r - center
877 r2 = dot_product(rvec, rvec)
878 r1 = sqrt(r2)
879 IF (r1 >= rs) THEN
880 IF (output_unit > 0) THEN
881 WRITE (output_unit, '(A,I6,A)') "Atom number :: ", i, " is out of the solvation sphere"
882 WRITE (output_unit, '(2(A,F12.6))') "Distance from the center::", r1, " Radius of the sphere::", rs
883 END IF
884 cpabort("")
885 END IF
886 locp(:, i) = 0.0_dp
887 IF (r1 /= 0.0_dp) THEN
888 DO l = 0, lmax
889 locp(l, i) = (r1**l*real(l + 1, kind=dp)*(eps_in - eps_out))/ &
890 (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 locp(0, i) = (eps_in - eps_out)/(rs*eps_in*eps_out)
895 END IF
896 pos(i) = r1
897 END DO
898 ! Particle-Particle potential energy matrix
899 cost = 0.0_dp
900 DO i = 1, SIZE(particle_set)
901 DO j = 1, i
902 factor = 0.0_dp
903 IF (pos(i)*pos(j) /= 0.0_dp) THEN
904 mycos = dot_product(particle_set(i)%r - center, particle_set(j)%r - center)/(pos(i)*pos(j))
905 IF (abs(mycos) > 1.0_dp) mycos = sign(1.0_dp, mycos)
906 DO l = 0, lmax
907 factor = factor + locp(l, i)*pos(j)**l*legendre(mycos, l, 0)
908 END DO
909 ELSE
910 factor = locp(0, i)
911 END IF
912 cost(i, j) = factor
913 cost(j, i) = factor
914 END DO
915 END DO
916 ! Computes the full potential energy matrix
917 idim = 0
918 DO iparticle2 = 1, SIZE(particle_set)
919 ip2 = (iparticle2 - 1)*SIZE(radii)
920 DO igauss2 = 1, SIZE(radii)
921 idim2 = ip2 + igauss2
922 DO iparticle1 = 1, iparticle2
923 ip1 = (iparticle1 - 1)*SIZE(radii)
924 istart_g = 1
925 IF (iparticle1 == iparticle2) istart_g = igauss2
926 DO igauss1 = istart_g, SIZE(radii)
927 idim1 = ip1 + igauss1
928 m(idim1, idim2) = cost(iparticle1, iparticle2)
929 m(idim2, idim1) = m(idim1, idim2)
930 END DO
931 END DO
932 END DO
933 END DO
934 DEALLOCATE (cost)
935 DEALLOCATE (pos)
936 DEALLOCATE (locp)
937 END SUBROUTINE solvation_ddapc_pot
938
939END MODULE cp_ddapc_methods
dgemm
static void dgemm(const char transa, const char transb, const int m, const int n, const int k, const double alpha, const double *a, const int lda, const double *b, const int ldb, const double beta, double *c, const int ldc)
Convenient wrapper to hide Fortran nature of dgemm_, swapping a and b.
Definition grid_cpu_task_list.c:195
kahan_sum::accurate_sum
Definition kahan_sum.F:41
mathlib::invert_matrix
Definition mathlib.F:69
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
cp_ddapc_methods
contains information regarding the decoupling/recoupling method of Bloechl
Definition cp_ddapc_methods.F:12
cp_ddapc_methods::cleanup_g_dot_rvec_sin_cos
subroutine, public cleanup_g_dot_rvec_sin_cos(g_dot_rvec_sin, g_dot_rvec_cos)
deallocate g_dot_rvec_* arrays
Definition cp_ddapc_methods.F:477
cp_ddapc_methods::build_der_b_vector
subroutine, public build_der_b_vector(dbv, gfunc, w, particle_set, radii, rho_tot_g, gcut, iparticle0)
Computes the derivative of B vector for the evaluation of the Pulay forces.
Definition cp_ddapc_methods.F:277
cp_ddapc_methods::build_b_vector
subroutine, public build_b_vector(bv, gfunc, w, particle_set, radii, rho_tot_g, gcut)
Computes the B vector for the solution of the linear system.
Definition cp_ddapc_methods.F:113
cp_ddapc_methods::build_a_matrix
subroutine, public build_a_matrix(am, gfunc, w, particle_set, radii, rho_tot_g, gcut, g_dot_rvec_sin, g_dot_rvec_cos)
Computes the A matrix for the solution of the linear system.
Definition cp_ddapc_methods.F:191
cp_ddapc_methods::solvation_ddapc_pot
subroutine, public solvation_ddapc_pot(solvation_section, particle_set, m, radii)
Evaluates the electrostatic potential due to a simple solvation model Spherical cavity in a dieletric...
Definition cp_ddapc_methods.F:809
cp_ddapc_methods::prep_g_dot_rvec_sin_cos
subroutine, public prep_g_dot_rvec_sin_cos(rho_tot_g, particle_set, gcut, g_dot_rvec_sin, g_dot_rvec_cos)
precompute sin(g.r) and cos(g.r) for quicker evaluations of sin(g.(r1-r2)) and cos(g....
Definition cp_ddapc_methods.F:492
cp_ddapc_methods::ewald_ddapc_pot
recursive subroutine, public ewald_ddapc_pot(cp_para_env, coeff, factor, cell, multipole_section, particle_set, m, radii)
Evaluates the Ewald term E2 and E3 energy term for the decoupling/coupling of periodic images.
Definition cp_ddapc_methods.F:627
cp_ddapc_methods::ddapc_eval_ami
subroutine, public ddapc_eval_ami(gami, c0, gfunc, w, particle_set, gcut, rho_tot_g, radii, iw, vol)
Computes the inverse AmI of the Am matrix.
Definition cp_ddapc_methods.F:544
cp_ddapc_methods::build_der_a_matrix_rows
subroutine, public build_der_a_matrix_rows(dam, gfunc, w, particle_set, radii, rho_tot_g, gcut, iparticle0, nparticles, g_dot_rvec_sin, g_dot_rvec_cos)
Computes the derivative of the A matrix for the evaluation of the Pulay forces.
Definition cp_ddapc_methods.F:369
cp_ddapc_methods::ddapc_eval_gfunc
subroutine, public ddapc_eval_gfunc(gfunc, w, gcut, rho_tot_g, radii)
...
Definition cp_ddapc_methods.F:62
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_logger_get_default_io_unit
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
Definition cp_log_handling.F:515
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:17
input_constants::weight_type_unit
integer, parameter, public weight_type_unit
Definition input_constants.F:69
input_constants::weight_type_mass
integer, parameter, public weight_type_mass
Definition input_constants.F:69
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_val_set
subroutine, public section_vals_val_set(section_vals, keyword_name, i_rep_section, i_rep_val, val, l_val, i_val, r_val, c_val, l_vals_ptr, i_vals_ptr, r_vals_ptr, c_vals_ptr)
sets the requested value
Definition input_section_types.F:1216
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition input_section_types.F:1040
kahan_sum
sums arrays of real/complex numbers with much reduced round-off as compared to a naive implementation...
Definition kahan_sum.F:29
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
list
An array-based list which grows on demand. When the internal array is full, a new array of twice the ...
Definition list.F:24
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::oorootpi
real(kind=dp), parameter, public oorootpi
Definition mathconstants.F:115
mathconstants::pi
real(kind=dp), parameter, public pi
Definition mathconstants.F:110
mathconstants::fourpi
real(kind=dp), parameter, public fourpi
Definition mathconstants.F:113
mathconstants::twopi
real(kind=dp), parameter, public twopi
Definition mathconstants.F:112
mathconstants::fac
real(kind=dp), dimension(0:maxfac), parameter, public fac
Definition mathconstants.F:37
mathlib
Collection of simple mathematical functions and subroutines.
Definition mathlib.F:15
mathlib::diamat_all
subroutine, public diamat_all(a, eigval, dac)
Diagonalize the symmetric n by n matrix a using the LAPACK library. Only the upper triangle of matrix...
Definition mathlib.F:372
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
pw_spline_utils
different utils that are useful to manipulate splines on the regular grid of a pw
Definition pw_spline_utils.F:18
pw_spline_utils::eval_interp_spl3_pbc
real(kind=dp) function, public eval_interp_spl3_pbc(vec, pw)
Evaluates the PBC interpolated Spline (pw) function on the generic input vector (vec)
Definition pw_spline_utils.F:3151
pw_types
Definition pw_types.F:24
spherical_harmonics
Calculate spherical harmonics.
Definition spherical_harmonics.F:27
spherical_harmonics::legendre
real(kind=dp) function, public legendre(x, l, m)
...
Definition spherical_harmonics.F:1193
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:126
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:763
particle_types::particle_type
Definition particle_types.F:35
pw_types::pw_c1d_gs_type
Definition pw_types.F:127
pw_types::pw_r3d_rs_type
Definition pw_types.F:62