(git:1f285aa)
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 ! **************************************************************************************************
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  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 
939 END MODULE cp_ddapc_methods
fac
static GRID_HOST_DEVICE double fac(const int i)
Factorial function, e.g. fac(5) = 5! = 120.
Definition: grid_common.h:48
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
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::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::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_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::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::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::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::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::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
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