(git:ccc2433)
cp_ddapc_forces.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 Density Derived atomic point charges from a QM calculation
10 !> (see J. Chem. Phys. Vol. 103 pp. 7422-7428)
11 !> \par History
12 !> 08.2005 created [tlaino]
13 !> \author Teodoro Laino
14 ! **************************************************************************************************
15 MODULE cp_ddapc_forces
16 
17  USE atomic_kind_types, ONLY: atomic_kind_type,&
18  get_atomic_kind_set
19  USE cell_types, ONLY: cell_type
20  USE cp_control_types, ONLY: ddapc_restraint_type
21  USE input_constants, ONLY: do_ddapc_constraint,&
22  do_ddapc_restraint,&
23  weight_type_mass,&
24  weight_type_unit
25  USE input_section_types, ONLY: section_vals_type,&
26  section_vals_val_get
27  USE kinds, ONLY: dp
28  USE mathconstants, ONLY: fourpi,&
29  pi,&
30  rootpi,&
31  twopi
32  USE message_passing, ONLY: mp_para_env_type
33  USE particle_types, ONLY: particle_type
34  USE pw_spline_utils, ONLY: eval_d_interp_spl3_pbc
35  USE pw_types, ONLY: pw_r3d_rs_type
36  USE qs_environment_types, ONLY: get_qs_env,&
37  qs_environment_type
38  USE qs_force_types, ONLY: qs_force_type
39  USE spherical_harmonics, ONLY: dlegendre,&
40  legendre
41 !NB for reducing results of calculations that use dq, which is now spread over nodes
42 #include "./base/base_uses.f90"
43 
44  IMPLICIT NONE
45  PRIVATE
46 
47  LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .true.
48  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cp_ddapc_forces'
49  PUBLIC :: ewald_ddapc_force, &
50  reset_ch_pulay, &
51  evaluate_restraint_functional, &
52  restraint_functional_force, &
53  solvation_ddapc_force
54 
55 CONTAINS
56 
57 ! **************************************************************************************************
58 !> \brief Evaluates the Ewald term E2 and E3 energy term for the decoupling/coupling
59 !> of periodic images
60 !> \param qs_env ...
61 !> \param coeff ...
62 !> \param apply_qmmm_periodic ...
63 !> \param factor ...
64 !> \param multipole_section ...
65 !> \param cell ...
66 !> \param particle_set ...
67 !> \param radii ...
68 !> \param dq ...
69 !> \param charges ...
70 !> \par History
71 !> 08.2005 created [tlaino]
72 !> \author Teodoro Laino
73 ! **************************************************************************************************
74  RECURSIVE SUBROUTINE ewald_ddapc_force(qs_env, coeff, apply_qmmm_periodic, &
75  factor, multipole_section, cell, particle_set, radii, dq, charges)
76  TYPE(qs_environment_type), POINTER :: qs_env
77  TYPE(pw_r3d_rs_type), INTENT(IN), POINTER :: coeff
78  LOGICAL, INTENT(IN) :: apply_qmmm_periodic
79  REAL(kind=dp), INTENT(IN) :: factor
80  TYPE(section_vals_type), POINTER :: multipole_section
81  TYPE(cell_type), POINTER :: cell
82  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
83  REAL(kind=dp), DIMENSION(:), POINTER :: radii
84  REAL(kind=dp), DIMENSION(:, :, :), OPTIONAL, &
85  POINTER :: dq
86  REAL(kind=dp), DIMENSION(:), OPTIONAL, POINTER :: charges
87 
88  CHARACTER(len=*), PARAMETER :: routinen = 'ewald_ddapc_force'
89 
90  INTEGER :: handle, ip1, ip2, iparticle1, &
91  iparticle2, k1, k2, k3, n_rep, nmax1, &
92  nmax2, nmax3, r1, r2, r3, rmax1, &
93  rmax2, rmax3
94  LOGICAL :: analyt
95  REAL(kind=dp) :: alpha, eps, fac, fac3, fs, galpha, gsq, &
96  gsqi, ij_fac, q1t, q2t, r, r2tmp, &
97  rcut, rcut2, t1, t2, tol, tol1
98  REAL(kind=dp), DIMENSION(3) :: drvec, fvec, gvec, ra, rvec
99  REAL(kind=dp), DIMENSION(:, :), POINTER :: d_el, m
100  TYPE(mp_para_env_type), POINTER :: para_env
101 
102  NULLIFY (d_el, m, para_env)
103  CALL timeset(routinen, handle)
104  CALL get_qs_env(qs_env, para_env=para_env)
105  cpassert(PRESENT(charges))
106  cpassert(ASSOCIATED(radii))
107  cpassert(cell%orthorhombic)
108  rcut = min(cell%hmat(1, 1), cell%hmat(2, 2), cell%hmat(3, 3))/2.0_dp
109  CALL section_vals_val_get(multipole_section, "RCUT", n_rep_val=n_rep)
110  IF (n_rep == 1) CALL section_vals_val_get(multipole_section, "RCUT", r_val=rcut)
111  CALL section_vals_val_get(multipole_section, "EWALD_PRECISION", r_val=eps)
112  CALL section_vals_val_get(multipole_section, "ANALYTICAL_GTERM", l_val=analyt)
113  rcut2 = rcut**2
114  !
115  ! Setting-up parameters for Ewald summation
116  !
117  eps = min(abs(eps), 0.5_dp)
118  tol = sqrt(abs(log(eps*rcut)))
119  alpha = sqrt(abs(log(eps*rcut*tol)))/rcut
120  galpha = 1.0_dp/(4.0_dp*alpha*alpha)
121  tol1 = sqrt(-log(eps*rcut*(2.0_dp*tol*alpha)**2))
122  nmax1 = nint(0.25_dp + cell%hmat(1, 1)*alpha*tol1/pi)
123  nmax2 = nint(0.25_dp + cell%hmat(2, 2)*alpha*tol1/pi)
124  nmax3 = nint(0.25_dp + cell%hmat(3, 3)*alpha*tol1/pi)
125 
126  rmax1 = ceiling(rcut/cell%hmat(1, 1))
127  rmax2 = ceiling(rcut/cell%hmat(2, 2))
128  rmax3 = ceiling(rcut/cell%hmat(3, 3))
129 
130  ALLOCATE (d_el(3, SIZE(particle_set)))
131  d_el = 0.0_dp
132  fac = 1.e0_dp/cell%deth
133  fac3 = fac/8.0_dp
134  fvec = twopi/(/cell%hmat(1, 1), cell%hmat(2, 2), cell%hmat(3, 3)/)
135  !
136  DO iparticle1 = 1, SIZE(particle_set)
137  !NB parallelization
138  IF (mod(iparticle1, para_env%num_pe) /= para_env%mepos) cycle
139  ip1 = (iparticle1 - 1)*SIZE(radii)
140  q1t = sum(charges(ip1 + 1:ip1 + SIZE(radii)))
141  DO iparticle2 = 1, iparticle1
142  ij_fac = 1.0_dp
143  IF (iparticle1 == iparticle2) ij_fac = 0.5_dp
144 
145  ip2 = (iparticle2 - 1)*SIZE(radii)
146  q2t = sum(charges(ip2 + 1:ip2 + SIZE(radii)))
147  !
148  ! Real-Space Contribution
149  !
150  rvec = particle_set(iparticle1)%r - particle_set(iparticle2)%r
151  IF (iparticle1 /= iparticle2) THEN
152  ra = rvec
153  r2tmp = dot_product(ra, ra)
154  IF (r2tmp <= rcut2) THEN
155  r = sqrt(r2tmp)
156  t1 = erfc(alpha*r)/r
157  drvec = ra/r*q1t*q2t*factor
158  t2 = -2.0_dp*alpha*exp(-alpha*alpha*r*r)/(r*rootpi) - t1/r
159  d_el(1:3, iparticle1) = d_el(1:3, iparticle1) - t2*drvec
160  d_el(1:3, iparticle2) = d_el(1:3, iparticle2) + t2*drvec
161  END IF
162  END IF
163  DO r1 = -rmax1, rmax1
164  DO r2 = -rmax2, rmax2
165  DO r3 = -rmax3, rmax3
166  IF ((r1 == 0) .AND. (r2 == 0) .AND. (r3 == 0)) cycle
167  ra(1) = rvec(1) + cell%hmat(1, 1)*r1
168  ra(2) = rvec(2) + cell%hmat(2, 2)*r2
169  ra(3) = rvec(3) + cell%hmat(3, 3)*r3
170  r2tmp = dot_product(ra, ra)
171  IF (r2tmp <= rcut2) THEN
172  r = sqrt(r2tmp)
173  t1 = erfc(alpha*r)/r
174  drvec = ra/r*q1t*q2t*factor*ij_fac
175  t2 = -2.0_dp*alpha*exp(-alpha*alpha*r*r)/(r*rootpi) - t1/r
176  d_el(1, iparticle1) = d_el(1, iparticle1) - t2*drvec(1)
177  d_el(2, iparticle1) = d_el(2, iparticle1) - t2*drvec(2)
178  d_el(3, iparticle1) = d_el(3, iparticle1) - t2*drvec(3)
179  d_el(1, iparticle2) = d_el(1, iparticle2) + t2*drvec(1)
180  d_el(2, iparticle2) = d_el(2, iparticle2) + t2*drvec(2)
181  d_el(3, iparticle2) = d_el(3, iparticle2) + t2*drvec(3)
182  END IF
183  END DO
184  END DO
185  END DO
186  !
187  ! G-space Contribution
188  !
189  IF (analyt) THEN
190  DO k1 = 0, nmax1
191  DO k2 = -nmax2, nmax2
192  DO k3 = -nmax3, nmax3
193  IF (k1 == 0 .AND. k2 == 0 .AND. k3 == 0) cycle
194  fs = 2.0_dp; IF (k1 == 0) fs = 1.0_dp
195  gvec = fvec*(/real(k1, kind=dp), real(k2, kind=dp), real(k3, kind=dp)/)
196  gsq = dot_product(gvec, gvec)
197  gsqi = fs/gsq
198  t1 = fac*gsqi*exp(-galpha*gsq)
199  t2 = -sin(dot_product(gvec, rvec))*t1*q1t*q2t*factor*fourpi
200  d_el(1:3, iparticle1) = d_el(1:3, iparticle1) - t2*gvec
201  d_el(1:3, iparticle2) = d_el(1:3, iparticle2) + t2*gvec
202  END DO
203  END DO
204  END DO
205  ELSE
206  gvec = eval_d_interp_spl3_pbc(rvec, coeff)*q1t*q2t*factor*fourpi
207  d_el(1:3, iparticle1) = d_el(1:3, iparticle1) - gvec
208  d_el(1:3, iparticle2) = d_el(1:3, iparticle2) + gvec
209  END IF
210  IF (iparticle1 /= iparticle2) THEN
211  ra = rvec
212  r = sqrt(dot_product(ra, ra))
213  t2 = -1.0_dp/(r*r)*factor
214  drvec = ra/r*q1t*q2t
215  d_el(1:3, iparticle1) = d_el(1:3, iparticle1) + t2*drvec
216  d_el(1:3, iparticle2) = d_el(1:3, iparticle2) - t2*drvec
217  END IF
218  END DO ! iparticle2
219  END DO ! iparticle1
220  !NB parallelization
221  CALL para_env%sum(d_el)
222  m => qs_env%cp_ddapc_env%Md
223  IF (apply_qmmm_periodic) m => qs_env%cp_ddapc_env%Mr
224  CALL cp_decpl_ddapc_forces(qs_env, m, charges, dq, d_el, particle_set)
225  DEALLOCATE (d_el)
226  CALL timestop(handle)
227  END SUBROUTINE ewald_ddapc_force
228 
229 ! **************************************************************************************************
230 !> \brief Evaluation of the pulay forces due to the fitted charge density
231 !> \param qs_env ...
232 !> \param M ...
233 !> \param charges ...
234 !> \param dq ...
235 !> \param d_el ...
236 !> \param particle_set ...
237 !> \par History
238 !> 08.2005 created [tlaino]
239 !> \author Teodoro Laino
240 ! **************************************************************************************************
241  SUBROUTINE cp_decpl_ddapc_forces(qs_env, M, charges, dq, d_el, particle_set)
242  TYPE(qs_environment_type), POINTER :: qs_env
243  REAL(kind=dp), DIMENSION(:, :), POINTER :: m
244  REAL(kind=dp), DIMENSION(:), POINTER :: charges
245  REAL(kind=dp), DIMENSION(:, :, :), POINTER :: dq
246  REAL(kind=dp), DIMENSION(:, :), POINTER :: d_el
247  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
248 
249  CHARACTER(len=*), PARAMETER :: routinen = 'cp_decpl_ddapc_forces'
250 
251  INTEGER :: handle, i, iatom, ikind, j, k, natom
252  INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
253  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: uv
254  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: chf
255  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
256  TYPE(mp_para_env_type), POINTER :: para_env
257  TYPE(qs_force_type), DIMENSION(:), POINTER :: force
258 
259  NULLIFY (para_env)
260  CALL timeset(routinen, handle)
261  natom = SIZE(particle_set)
262  CALL get_qs_env(qs_env=qs_env, &
263  atomic_kind_set=atomic_kind_set, &
264  para_env=para_env, &
265  force=force)
266  ALLOCATE (chf(3, natom))
267  CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
268  atom_of_kind=atom_of_kind, &
269  kind_of=kind_of)
270 
271  ALLOCATE (uv(SIZE(m, 1)))
272  uv(:) = matmul(m, charges)
273  DO k = 1, natom
274  DO j = 1, 3
275  chf(j, k) = dot_product(uv, dq(:, k, j))
276  END DO
277  END DO
278  !NB now that get_ddapc returns dq that's spread over nodes, must reduce chf here
279  CALL para_env%sum(chf)
280  DO iatom = 1, natom
281  ikind = kind_of(iatom)
282  i = atom_of_kind(iatom)
283  force(ikind)%ch_pulay(1:3, i) = force(ikind)%ch_pulay(1:3, i) + chf(1:3, iatom) + d_el(1:3, iatom)
284  END DO
285  DEALLOCATE (chf)
286  DEALLOCATE (uv)
287  CALL timestop(handle)
288  END SUBROUTINE cp_decpl_ddapc_forces
289 
290 ! **************************************************************************************************
291 !> \brief Evaluation of the pulay forces due to the fitted charge density
292 !> \param qs_env ...
293 !> \par History
294 !> 08.2005 created [tlaino]
295 !> \author Teodoro Laino
296 ! **************************************************************************************************
297  SUBROUTINE reset_ch_pulay(qs_env)
298  TYPE(qs_environment_type), POINTER :: qs_env
299 
300  CHARACTER(len=*), PARAMETER :: routinen = 'reset_ch_pulay'
301 
302  INTEGER :: handle, ind
303  TYPE(qs_force_type), DIMENSION(:), POINTER :: force
304 
305  CALL timeset(routinen, handle)
306  CALL get_qs_env(qs_env=qs_env, &
307  force=force)
308  DO ind = 1, SIZE(force)
309  force(ind)%ch_pulay = 0.0_dp
310  END DO
311  CALL timestop(handle)
312  END SUBROUTINE reset_ch_pulay
313 
314 ! **************************************************************************************************
315 !> \brief computes energy and derivatives given a set of charges
316 !> \param ddapc_restraint_control ...
317 !> \param n_gauss ...
318 !> \param uv derivate of energy wrt the corresponding charge entry
319 !> \param charges current value of the charges (one number for each gaussian used)
320 !>
321 !> \param energy_res energy due to the restraint
322 !> \par History
323 !> 02.2006 [Joost VandeVondele]
324 !> modified [Teo]
325 !> \note
326 !> should be easy to adapt for other specialized cases
327 ! **************************************************************************************************
328  SUBROUTINE evaluate_restraint_functional(ddapc_restraint_control, n_gauss, uv, &
329  charges, energy_res)
330  TYPE(ddapc_restraint_type), INTENT(INOUT) :: ddapc_restraint_control
331  INTEGER, INTENT(in) :: n_gauss
332  REAL(kind=dp), DIMENSION(:) :: uv
333  REAL(kind=dp), DIMENSION(:), POINTER :: charges
334  REAL(kind=dp), INTENT(INOUT) :: energy_res
335 
336  INTEGER :: i, ind
337  REAL(kind=dp) :: de, order_p
338 
339 ! order parameter (i.e. the sum of the charges of the selected atoms)
340 
341  order_p = 0.0_dp
342  DO i = 1, ddapc_restraint_control%natoms
343  ind = (ddapc_restraint_control%atoms(i) - 1)*n_gauss
344  order_p = order_p + ddapc_restraint_control%coeff(i)*sum(charges(ind + 1:ind + n_gauss))
345  END DO
346  ddapc_restraint_control%ddapc_order_p = order_p
347 
348  SELECT CASE (ddapc_restraint_control%functional_form)
349  CASE (do_ddapc_restraint)
350  ! the restraint energy
351  energy_res = ddapc_restraint_control%strength*(order_p - ddapc_restraint_control%target)**2.0_dp
352 
353  ! derivative of the energy
354  de = 2.0_dp*ddapc_restraint_control%strength*(order_p - ddapc_restraint_control%target)
355  DO i = 1, ddapc_restraint_control%natoms
356  ind = (ddapc_restraint_control%atoms(i) - 1)*n_gauss
357  uv(ind + 1:ind + n_gauss) = de*ddapc_restraint_control%coeff(i)
358  END DO
359  CASE (do_ddapc_constraint)
360  energy_res = ddapc_restraint_control%strength*(order_p - ddapc_restraint_control%target)
361 
362  ! derivative of the energy
363  DO i = 1, ddapc_restraint_control%natoms
364  ind = (ddapc_restraint_control%atoms(i) - 1)*n_gauss
365  uv(ind + 1:ind + n_gauss) = ddapc_restraint_control%strength*ddapc_restraint_control%coeff(i)
366  END DO
367 
368  CASE DEFAULT
369  cpabort("")
370  END SELECT
371 
372  END SUBROUTINE evaluate_restraint_functional
373 
374 ! **************************************************************************************************
375 !> \brief computes derivatives for DDAPC restraint
376 !> \param qs_env ...
377 !> \param ddapc_restraint_control ...
378 !> \param dq ...
379 !> \param charges ...
380 !> \param n_gauss ...
381 !> \param particle_set ...
382 !> \par History
383 !> 02.2006 [Joost VandeVondele]
384 !> modified [Teo]
385 !> \note
386 !> should be easy to adapt for other specialized cases
387 ! **************************************************************************************************
388  SUBROUTINE restraint_functional_force(qs_env, ddapc_restraint_control, dq, charges, &
389  n_gauss, particle_set)
390 
391  TYPE(qs_environment_type), POINTER :: qs_env
392  TYPE(ddapc_restraint_type), INTENT(INOUT) :: ddapc_restraint_control
393  REAL(kind=dp), DIMENSION(:, :, :), POINTER :: dq
394  REAL(kind=dp), DIMENSION(:), POINTER :: charges
395  INTEGER, INTENT(in) :: n_gauss
396  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
397 
398  CHARACTER(len=*), PARAMETER :: routinen = 'restraint_functional_force'
399 
400  INTEGER :: handle, i, iatom, ikind, j, k, natom
401  INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
402  REAL(kind=dp) :: dum
403  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: uv
404  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: chf
405  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
406  TYPE(mp_para_env_type), POINTER :: para_env
407  TYPE(qs_force_type), DIMENSION(:), POINTER :: force
408 
409  NULLIFY (para_env)
410  CALL timeset(routinen, handle)
411  natom = SIZE(particle_set)
412  CALL get_qs_env(qs_env=qs_env, &
413  atomic_kind_set=atomic_kind_set, &
414  para_env=para_env, &
415  force=force)
416  ALLOCATE (chf(3, natom))
417  CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
418  atom_of_kind=atom_of_kind, &
419  kind_of=kind_of)
420 
421  ALLOCATE (uv(SIZE(dq, 1)))
422  uv = 0.0_dp
423  CALL evaluate_restraint_functional(ddapc_restraint_control, n_gauss, uv, &
424  charges, dum)
425  DO k = 1, natom
426  DO j = 1, 3
427  chf(j, k) = dot_product(uv, dq(:, k, j))
428  END DO
429  END DO
430  !NB now that get_ddapc returns dq that's spread over nodes, must reduce chf here
431  CALL para_env%sum(chf)
432  DO iatom = 1, natom
433  ikind = kind_of(iatom)
434  i = atom_of_kind(iatom)
435  force(ikind)%ch_pulay(1:3, i) = force(ikind)%ch_pulay(1:3, i) + chf(1:3, iatom)
436  END DO
437  DEALLOCATE (chf)
438  DEALLOCATE (uv)
439  CALL timestop(handle)
440 
441  END SUBROUTINE restraint_functional_force
442 
443 ! **************************************************************************************************
444 !> \brief Evaluates the electrostatic potential due to a simple solvation model
445 !> Spherical cavity in a dieletric medium
446 !> \param qs_env ...
447 !> \param solvation_section ...
448 !> \param particle_set ...
449 !> \param radii ...
450 !> \param dq ...
451 !> \param charges ...
452 !> \par History
453 !> 08.2006 created [tlaino]
454 !> \author Teodoro Laino
455 ! **************************************************************************************************
456  SUBROUTINE solvation_ddapc_force(qs_env, solvation_section, particle_set, &
457  radii, dq, charges)
458  TYPE(qs_environment_type), POINTER :: qs_env
459  TYPE(section_vals_type), POINTER :: solvation_section
460  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
461  REAL(kind=dp), DIMENSION(:), POINTER :: radii
462  REAL(kind=dp), DIMENSION(:, :, :), OPTIONAL, &
463  POINTER :: dq
464  REAL(kind=dp), DIMENSION(:), OPTIONAL, POINTER :: charges
465 
466  INTEGER :: i, ip1, ip2, iparticle1, iparticle2, l, &
467  lmax, n_rep1, n_rep2, weight
468  INTEGER, DIMENSION(:), POINTER :: list
469  LOGICAL :: fixed_center
470  REAL(kind=dp) :: center(3), dcos1(3), dcos2(3), dpos1(3), dpos2(3), eps_in, eps_out, &
471  factor1(3), factor2(3), lr, mass, mycos, pos1, pos1i, pos2, pos2i, ptcos, q1t, q2t, &
472  r1(3), r1s, r2(3), r2s, rs, rvec(3)
473  REAL(kind=dp), DIMENSION(:), POINTER :: pos, r0
474  REAL(kind=dp), DIMENSION(:, :), POINTER :: d_el, locp, m
475 
476  fixed_center = .false.
477  NULLIFY (d_el, m)
478  eps_in = 1.0_dp
479  CALL section_vals_val_get(solvation_section, "EPS_OUT", r_val=eps_out)
480  CALL section_vals_val_get(solvation_section, "LMAX", i_val=lmax)
481  CALL section_vals_val_get(solvation_section, "SPHERE%RADIUS", r_val=rs)
482  CALL section_vals_val_get(solvation_section, "SPHERE%CENTER%XYZ", n_rep_val=n_rep1)
483  IF (n_rep1 /= 0) THEN
484  CALL section_vals_val_get(solvation_section, "SPHERE%CENTER%XYZ", r_vals=r0)
485  center = r0
486  ELSE
487  CALL section_vals_val_get(solvation_section, "SPHERE%CENTER%ATOM_LIST", &
488  n_rep_val=n_rep2)
489  IF (n_rep2 /= 0) THEN
490  CALL section_vals_val_get(solvation_section, "SPHERE%CENTER%ATOM_LIST", i_vals=list)
491  CALL section_vals_val_get(solvation_section, "SPHERE%CENTER%WEIGHT_TYPE", i_val=weight)
492  ALLOCATE (r0(SIZE(list)))
493  SELECT CASE (weight)
494  CASE (weight_type_unit)
495  r0 = 0.0_dp
496  DO i = 1, SIZE(list)
497  r0 = r0 + particle_set(list(i))%r
498  END DO
499  r0 = r0/real(SIZE(list), kind=dp)
500  CASE (weight_type_mass)
501  r0 = 0.0_dp
502  mass = 0.0_dp
503  DO i = 1, SIZE(list)
504  r0 = r0 + particle_set(list(i))%r*particle_set(list(i))%atomic_kind%mass
505  mass = mass + particle_set(list(i))%atomic_kind%mass
506  END DO
507  r0 = r0/mass
508  END SELECT
509  center = r0
510  DEALLOCATE (r0)
511  END IF
512  END IF
513  cpassert(n_rep1 /= 0 .OR. n_rep2 /= 0)
514  ! Potential calculation
515  ALLOCATE (locp(0:lmax, SIZE(particle_set)))
516  ALLOCATE (pos(SIZE(particle_set)))
517  ALLOCATE (d_el(3, SIZE(particle_set)))
518  d_el = 0.0_dp
519  ! Determining the single atomic contribution to the dielectric dipole
520  DO i = 1, SIZE(particle_set)
521  rvec = particle_set(i)%r - center
522  r2s = dot_product(rvec, rvec)
523  r1s = sqrt(r2s)
524  locp(:, i) = 0.0_dp
525  IF (r1s /= 0.0_dp) THEN
526  DO l = 0, lmax
527  locp(l, i) = (r1s**l*real(l + 1, kind=dp)*(eps_in - eps_out))/ &
528  (rs**(2*l + 1)*eps_in*(real(l, kind=dp)*eps_in + real(l + 1, kind=dp)*eps_out))
529  END DO
530  END IF
531  pos(i) = r1s
532  END DO
533  ! Computes the full derivatives of the interaction energy
534  DO iparticle1 = 1, SIZE(particle_set)
535  ip1 = (iparticle1 - 1)*SIZE(radii)
536  q1t = sum(charges(ip1 + 1:ip1 + SIZE(radii)))
537  DO iparticle2 = 1, iparticle1
538  ip2 = (iparticle2 - 1)*SIZE(radii)
539  q2t = sum(charges(ip2 + 1:ip2 + SIZE(radii)))
540  !
541  r1 = particle_set(iparticle1)%r - center
542  r2 = particle_set(iparticle2)%r - center
543  pos1 = pos(iparticle1)
544  pos2 = pos(iparticle2)
545  factor1 = 0.0_dp
546  factor2 = 0.0_dp
547  IF (pos1*pos2 /= 0.0_dp) THEN
548  pos1i = 1.0_dp/pos1
549  pos2i = 1.0_dp/pos2
550  dpos1 = pos1i*r1
551  dpos2 = pos2i*r2
552  ptcos = dot_product(r1, r2)
553  mycos = ptcos/(pos1*pos2)
554  IF (abs(mycos) > 1.0_dp) mycos = sign(1.0_dp, mycos)
555  dcos1 = (r2*(pos1*pos2) - pos2*dpos1*ptcos)/(pos1*pos2)**2
556  dcos2 = (r1*(pos1*pos2) - pos1*dpos2*ptcos)/(pos1*pos2)**2
557 
558  DO l = 1, lmax
559  lr = real(l, kind=dp)
560  factor1 = factor1 + lr*locp(l, iparticle2)*pos1**(l - 1)*legendre(mycos, l, 0)*dpos1 &
561  + locp(l, iparticle2)*pos1**l*dlegendre(mycos, l, 0)*dcos1
562  factor2 = factor2 + lr*locp(l, iparticle1)*pos2**(l - 1)*legendre(mycos, l, 0)*dpos2 &
563  + locp(l, iparticle1)*pos2**l*dlegendre(mycos, l, 0)*dcos2
564  END DO
565  END IF
566  factor1 = factor1*q1t*q2t
567  factor2 = factor2*q1t*q2t
568  d_el(1:3, iparticle1) = d_el(1:3, iparticle1) + 0.5_dp*factor1
569  d_el(1:3, iparticle2) = d_el(1:3, iparticle2) + 0.5_dp*factor2
570  END DO
571  END DO
572  DEALLOCATE (pos)
573  DEALLOCATE (locp)
574  m => qs_env%cp_ddapc_env%Ms
575  CALL cp_decpl_ddapc_forces(qs_env, m, charges, dq, d_el, particle_set)
576  DEALLOCATE (d_el)
577  END SUBROUTINE solvation_ddapc_force
578 
579 END MODULE cp_ddapc_forces
fac
static GRID_HOST_DEVICE double fac(const int i)
Factorial function, e.g. fac(5) = 5! = 120.
Definition: grid_common.h:48
atomic_kind_types
Define the atomic kind types and their sub types.
Definition: atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind_set
subroutine, public get_atomic_kind_set(atomic_kind_set, atom_of_kind, kind_of, natom_of_kind, maxatom, natom, nshell, fist_potential_present, shell_present, shell_adiabatic, shell_check_distance, damping_present)
Get attributes of an atomic kind set.
Definition: atomic_kind_types.F:223
cell_types
Handles all functions related to the CELL.
Definition: cell_types.F:15
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition: cp_control_types.F:12
cp_ddapc_forces
Density Derived atomic point charges from a QM calculation (see J. Chem. Phys. Vol....
Definition: cp_ddapc_forces.F:15
cp_ddapc_forces::solvation_ddapc_force
subroutine, public solvation_ddapc_force(qs_env, solvation_section, particle_set, radii, dq, charges)
Evaluates the electrostatic potential due to a simple solvation model Spherical cavity in a dieletric...
Definition: cp_ddapc_forces.F:458
cp_ddapc_forces::ewald_ddapc_force
recursive subroutine, public ewald_ddapc_force(qs_env, coeff, apply_qmmm_periodic, factor, multipole_section, cell, particle_set, radii, dq, charges)
Evaluates the Ewald term E2 and E3 energy term for the decoupling/coupling of periodic images.
Definition: cp_ddapc_forces.F:76
cp_ddapc_forces::restraint_functional_force
subroutine, public restraint_functional_force(qs_env, ddapc_restraint_control, dq, charges, n_gauss, particle_set)
computes derivatives for DDAPC restraint
Definition: cp_ddapc_forces.F:390
cp_ddapc_forces::evaluate_restraint_functional
subroutine, public evaluate_restraint_functional(ddapc_restraint_control, n_gauss, uv, charges, energy_res)
computes energy and derivatives given a set of charges
Definition: cp_ddapc_forces.F:330
cp_ddapc_forces::reset_ch_pulay
subroutine, public reset_ch_pulay(qs_env)
Evaluation of the pulay forces due to the fitted charge density.
Definition: cp_ddapc_forces.F:298
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition: input_constants.F:17
input_constants::do_ddapc_constraint
integer, parameter, public do_ddapc_constraint
Definition: input_constants.F:714
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_constants::do_ddapc_restraint
integer, parameter, public do_ddapc_restraint
Definition: input_constants.F:714
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_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
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::pi
real(kind=dp), parameter, public pi
Definition: mathconstants.F:110
mathconstants::rootpi
real(kind=dp), parameter, public rootpi
Definition: mathconstants.F:114
mathconstants::fourpi
real(kind=dp), parameter, public fourpi
Definition: mathconstants.F:113
mathconstants::twopi
real(kind=dp), parameter, public twopi
Definition: mathconstants.F:112
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_d_interp_spl3_pbc
real(kind=dp) function, dimension(3), public eval_d_interp_spl3_pbc(vec, pw)
Evaluates the derivatives of the PBC interpolated Spline (pw) function on the generic input vector (v...
Definition: pw_spline_utils.F:3300
pw_types
Definition: pw_types.F:24
qs_environment_types
Definition: qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_RI_aux_kp, matrix_s, matrix_s_RI_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, WannierCentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Get the QUICKSTEP environment.
Definition: qs_environment_types.F:502
qs_force_types
Definition: qs_force_types.F:13
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
spherical_harmonics::dlegendre
real(kind=dp) function, public dlegendre(x, l, m)
...
Definition: spherical_harmonics.F:1337