(git:ed6f26b)
Loading...
Searching...
No Matches
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-2025 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! **************************************************************************************************
15MODULE 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
55CONTAINS
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
579END MODULE cp_ddapc_forces
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:720
input_constants::weight_type_unit
integer, parameter, public weight_type_unit
Definition input_constants.F:70
input_constants::weight_type_mass
integer, parameter, public weight_type_mass
Definition input_constants.F:70
input_constants::do_ddapc_restraint
integer, parameter, public do_ddapc_restraint
Definition input_constants.F:720
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:1047
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
mathconstants::fac
real(kind=dp), dimension(0:maxfac), parameter, public fac
Definition mathconstants.F:37
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_pp, 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, harris_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, eeq, rhs)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:514
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
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
cp_control_types::ddapc_restraint_type
Definition cp_control_types.F:291
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
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_r3d_rs_type
Definition pw_types.F:62
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:217
qs_force_types::qs_force_type
Definition qs_force_types.F:28