Line data Source code
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 : ! **************************************************************************************************
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 136 : 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 136 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: d_el, M
100 : TYPE(mp_para_env_type), POINTER :: para_env
101 :
102 136 : NULLIFY (d_el, M, para_env)
103 136 : CALL timeset(routineN, handle)
104 136 : CALL get_qs_env(qs_env, para_env=para_env)
105 136 : CPASSERT(PRESENT(charges))
106 136 : CPASSERT(ASSOCIATED(radii))
107 136 : CPASSERT(cell%orthorhombic)
108 136 : rcut = MIN(cell%hmat(1, 1), cell%hmat(2, 2), cell%hmat(3, 3))/2.0_dp
109 136 : CALL section_vals_val_get(multipole_section, "RCUT", n_rep_val=n_rep)
110 136 : IF (n_rep == 1) CALL section_vals_val_get(multipole_section, "RCUT", r_val=rcut)
111 136 : CALL section_vals_val_get(multipole_section, "EWALD_PRECISION", r_val=eps)
112 136 : CALL section_vals_val_get(multipole_section, "ANALYTICAL_GTERM", l_val=analyt)
113 136 : rcut2 = rcut**2
114 : !
115 : ! Setting-up parameters for Ewald summation
116 : !
117 136 : eps = MIN(ABS(eps), 0.5_dp)
118 136 : tol = SQRT(ABS(LOG(eps*rcut)))
119 136 : alpha = SQRT(ABS(LOG(eps*rcut*tol)))/rcut
120 136 : galpha = 1.0_dp/(4.0_dp*alpha*alpha)
121 136 : tol1 = SQRT(-LOG(eps*rcut*(2.0_dp*tol*alpha)**2))
122 136 : nmax1 = NINT(0.25_dp + cell%hmat(1, 1)*alpha*tol1/pi)
123 136 : nmax2 = NINT(0.25_dp + cell%hmat(2, 2)*alpha*tol1/pi)
124 136 : nmax3 = NINT(0.25_dp + cell%hmat(3, 3)*alpha*tol1/pi)
125 :
126 136 : rmax1 = CEILING(rcut/cell%hmat(1, 1))
127 136 : rmax2 = CEILING(rcut/cell%hmat(2, 2))
128 136 : rmax3 = CEILING(rcut/cell%hmat(3, 3))
129 :
130 408 : ALLOCATE (d_el(3, SIZE(particle_set)))
131 1704 : d_el = 0.0_dp
132 136 : fac = 1.e0_dp/cell%deth
133 136 : fac3 = fac/8.0_dp
134 544 : fvec = twopi/(/cell%hmat(1, 1), cell%hmat(2, 2), cell%hmat(3, 3)/)
135 : !
136 528 : DO iparticle1 = 1, SIZE(particle_set)
137 : !NB parallelization
138 392 : IF (MOD(iparticle1, para_env%num_pe) /= para_env%mepos) CYCLE
139 212 : ip1 = (iparticle1 - 1)*SIZE(radii)
140 848 : q1t = SUM(charges(ip1 + 1:ip1 + SIZE(radii)))
141 854 : DO iparticle2 = 1, iparticle1
142 506 : ij_fac = 1.0_dp
143 506 : IF (iparticle1 == iparticle2) ij_fac = 0.5_dp
144 :
145 506 : ip2 = (iparticle2 - 1)*SIZE(radii)
146 2024 : q2t = SUM(charges(ip2 + 1:ip2 + SIZE(radii)))
147 : !
148 : ! Real-Space Contribution
149 : !
150 2024 : rvec = particle_set(iparticle1)%r - particle_set(iparticle2)%r
151 506 : IF (iparticle1 /= iparticle2) THEN
152 294 : ra = rvec
153 1176 : r2tmp = DOT_PRODUCT(ra, ra)
154 294 : IF (r2tmp <= rcut2) THEN
155 182 : r = SQRT(r2tmp)
156 182 : t1 = erfc(alpha*r)/r
157 728 : drvec = ra/r*q1t*q2t*factor
158 182 : t2 = -2.0_dp*alpha*EXP(-alpha*alpha*r*r)/(r*rootpi) - t1/r
159 728 : d_el(1:3, iparticle1) = d_el(1:3, iparticle1) - t2*drvec
160 728 : d_el(1:3, iparticle2) = d_el(1:3, iparticle2) + t2*drvec
161 : END IF
162 : END IF
163 2096 : DO r1 = -rmax1, rmax1
164 7226 : DO r2 = -rmax2, rmax2
165 23910 : DO r3 = -rmax3, rmax3
166 17190 : IF ((r1 == 0) .AND. (r2 == 0) .AND. (r3 == 0)) CYCLE
167 16684 : ra(1) = rvec(1) + cell%hmat(1, 1)*r1
168 16684 : ra(2) = rvec(2) + cell%hmat(2, 2)*r2
169 16684 : ra(3) = rvec(3) + cell%hmat(3, 3)*r3
170 66736 : r2tmp = DOT_PRODUCT(ra, ra)
171 21814 : IF (r2tmp <= rcut2) THEN
172 1020 : r = SQRT(r2tmp)
173 1020 : t1 = erfc(alpha*r)/r
174 4080 : drvec = ra/r*q1t*q2t*factor*ij_fac
175 1020 : t2 = -2.0_dp*alpha*EXP(-alpha*alpha*r*r)/(r*rootpi) - t1/r
176 1020 : d_el(1, iparticle1) = d_el(1, iparticle1) - t2*drvec(1)
177 1020 : d_el(2, iparticle1) = d_el(2, iparticle1) - t2*drvec(2)
178 1020 : d_el(3, iparticle1) = d_el(3, iparticle1) - t2*drvec(3)
179 1020 : d_el(1, iparticle2) = d_el(1, iparticle2) + t2*drvec(1)
180 1020 : d_el(2, iparticle2) = d_el(2, iparticle2) + t2*drvec(2)
181 1020 : 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 506 : IF (analyt) THEN
190 2599 : DO k1 = 0, nmax1
191 66864 : DO k2 = -nmax2, nmax2
192 1964583 : DO k3 = -nmax3, nmax3
193 1897925 : IF (k1 == 0 .AND. k2 == 0 .AND. k3 == 0) CYCLE
194 1897719 : fs = 2.0_dp; IF (k1 == 0) fs = 1.0_dp
195 7590876 : gvec = fvec*(/REAL(k1, KIND=dp), REAL(k2, KIND=dp), REAL(k3, KIND=dp)/)
196 7590876 : gsq = DOT_PRODUCT(gvec, gvec)
197 1897719 : gsqi = fs/gsq
198 1897719 : t1 = fac*gsqi*EXP(-galpha*gsq)
199 7590876 : t2 = -SIN(DOT_PRODUCT(gvec, rvec))*t1*q1t*q2t*factor*fourpi
200 7590876 : d_el(1:3, iparticle1) = d_el(1:3, iparticle1) - t2*gvec
201 7655141 : d_el(1:3, iparticle2) = d_el(1:3, iparticle2) + t2*gvec
202 : END DO
203 : END DO
204 : END DO
205 : ELSE
206 1200 : gvec = Eval_d_Interp_Spl3_pbc(rvec, coeff)*q1t*q2t*factor*fourpi
207 1200 : d_el(1:3, iparticle1) = d_el(1:3, iparticle1) - gvec
208 1200 : d_el(1:3, iparticle2) = d_el(1:3, iparticle2) + gvec
209 : END IF
210 898 : IF (iparticle1 /= iparticle2) THEN
211 294 : ra = rvec
212 1176 : r = SQRT(DOT_PRODUCT(ra, ra))
213 294 : t2 = -1.0_dp/(r*r)*factor
214 1176 : drvec = ra/r*q1t*q2t
215 1176 : d_el(1:3, iparticle1) = d_el(1:3, iparticle1) + t2*drvec
216 1176 : 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 3272 : CALL para_env%sum(d_el)
222 136 : M => qs_env%cp_ddapc_env%Md
223 136 : IF (apply_qmmm_periodic) M => qs_env%cp_ddapc_env%Mr
224 136 : CALL cp_decpl_ddapc_forces(qs_env, M, charges, dq, d_el, particle_set)
225 136 : DEALLOCATE (d_el)
226 136 : CALL timestop(handle)
227 408 : 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 150 : 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 150 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
253 150 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: uv
254 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: chf
255 150 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
256 : TYPE(mp_para_env_type), POINTER :: para_env
257 150 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
258 :
259 150 : NULLIFY (para_env)
260 150 : CALL timeset(routineN, handle)
261 150 : 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 150 : force=force)
266 450 : ALLOCATE (chf(3, natom))
267 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
268 : atom_of_kind=atom_of_kind, &
269 150 : kind_of=kind_of)
270 :
271 450 : ALLOCATE (uv(SIZE(M, 1)))
272 47172 : uv(:) = MATMUL(M, charges)
273 580 : DO k = 1, natom
274 1870 : DO j = 1, 3
275 16534 : 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 150 : CALL para_env%sum(chf)
280 580 : DO iatom = 1, natom
281 430 : ikind = kind_of(iatom)
282 430 : i = atom_of_kind(iatom)
283 3590 : 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 150 : DEALLOCATE (chf)
286 150 : DEALLOCATE (uv)
287 150 : CALL timestop(handle)
288 300 : 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 120 : 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 120 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
304 :
305 120 : CALL timeset(routineN, handle)
306 : CALL get_qs_env(qs_env=qs_env, &
307 120 : force=force)
308 324 : DO ind = 1, SIZE(force)
309 1612 : force(ind)%ch_pulay = 0.0_dp
310 : END DO
311 120 : CALL timestop(handle)
312 120 : 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 872 : 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 872 : order_p = 0.0_dp
342 1882 : DO I = 1, ddapc_restraint_control%natoms
343 1010 : ind = (ddapc_restraint_control%atoms(I) - 1)*n_gauss
344 3532 : order_p = order_p + ddapc_restraint_control%coeff(I)*SUM(charges(ind + 1:ind + n_gauss))
345 : END DO
346 872 : ddapc_restraint_control%ddapc_order_p = order_p
347 :
348 1140 : SELECT CASE (ddapc_restraint_control%functional_form)
349 : CASE (do_ddapc_restraint)
350 : ! the restraint energy
351 268 : energy_res = ddapc_restraint_control%strength*(order_p - ddapc_restraint_control%target)**2.0_dp
352 :
353 : ! derivative of the energy
354 268 : dE = 2.0_dp*ddapc_restraint_control%strength*(order_p - ddapc_restraint_control%target)
355 536 : DO I = 1, ddapc_restraint_control%natoms
356 268 : ind = (ddapc_restraint_control%atoms(I) - 1)*n_gauss
357 1340 : uv(ind + 1:ind + n_gauss) = dE*ddapc_restraint_control%coeff(I)
358 : END DO
359 : CASE (do_ddapc_constraint)
360 604 : energy_res = ddapc_restraint_control%strength*(order_p - ddapc_restraint_control%target)
361 :
362 : ! derivative of the energy
363 1346 : DO I = 1, ddapc_restraint_control%natoms
364 742 : ind = (ddapc_restraint_control%atoms(I) - 1)*n_gauss
365 2192 : uv(ind + 1:ind + n_gauss) = ddapc_restraint_control%strength*ddapc_restraint_control%coeff(I)
366 : END DO
367 :
368 : CASE DEFAULT
369 872 : CPABORT("")
370 : END SELECT
371 :
372 872 : 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 36 : 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 36 : 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 36 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
406 : TYPE(mp_para_env_type), POINTER :: para_env
407 36 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
408 :
409 36 : NULLIFY (para_env)
410 36 : CALL timeset(routineN, handle)
411 36 : 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 36 : force=force)
416 108 : ALLOCATE (chf(3, natom))
417 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
418 : atom_of_kind=atom_of_kind, &
419 36 : kind_of=kind_of)
420 :
421 108 : ALLOCATE (uv(SIZE(dq, 1)))
422 174 : uv = 0.0_dp
423 : CALL evaluate_restraint_functional(ddapc_restraint_control, n_gauss, uv, &
424 36 : charges, dum)
425 118 : DO k = 1, natom
426 364 : DO j = 1, 3
427 1318 : 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 36 : CALL para_env%sum(chf)
432 118 : DO iatom = 1, natom
433 82 : ikind = kind_of(iatom)
434 82 : i = atom_of_kind(iatom)
435 364 : force(ikind)%ch_pulay(1:3, i) = force(ikind)%ch_pulay(1:3, i) + chf(1:3, iatom)
436 : END DO
437 36 : DEALLOCATE (chf)
438 36 : DEALLOCATE (uv)
439 36 : CALL timestop(handle)
440 :
441 72 : 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 14 : 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 14 : 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 14 : REAL(KIND=dp), DIMENSION(:), POINTER :: pos, R0
474 14 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: d_el, LocP, M
475 :
476 14 : fixed_center = .FALSE.
477 14 : NULLIFY (d_el, M)
478 14 : eps_in = 1.0_dp
479 14 : CALL section_vals_val_get(solvation_section, "EPS_OUT", r_val=eps_out)
480 14 : CALL section_vals_val_get(solvation_section, "LMAX", i_val=lmax)
481 14 : CALL section_vals_val_get(solvation_section, "SPHERE%RADIUS", r_val=Rs)
482 14 : CALL section_vals_val_get(solvation_section, "SPHERE%CENTER%XYZ", n_rep_val=n_rep1)
483 14 : IF (n_rep1 /= 0) THEN
484 14 : CALL section_vals_val_get(solvation_section, "SPHERE%CENTER%XYZ", r_vals=R0)
485 56 : center = R0
486 : ELSE
487 : CALL section_vals_val_get(solvation_section, "SPHERE%CENTER%ATOM_LIST", &
488 0 : n_rep_val=n_rep2)
489 0 : IF (n_rep2 /= 0) THEN
490 0 : CALL section_vals_val_get(solvation_section, "SPHERE%CENTER%ATOM_LIST", i_vals=list)
491 0 : CALL section_vals_val_get(solvation_section, "SPHERE%CENTER%WEIGHT_TYPE", i_val=weight)
492 0 : ALLOCATE (R0(SIZE(list)))
493 : SELECT CASE (weight)
494 : CASE (weight_type_unit)
495 0 : R0 = 0.0_dp
496 0 : DO i = 1, SIZE(list)
497 0 : R0 = R0 + particle_set(list(i))%r
498 : END DO
499 0 : R0 = R0/REAL(SIZE(list), KIND=dp)
500 : CASE (weight_type_mass)
501 0 : R0 = 0.0_dp
502 0 : mass = 0.0_dp
503 0 : DO i = 1, SIZE(list)
504 0 : R0 = R0 + particle_set(list(i))%r*particle_set(list(i))%atomic_kind%mass
505 0 : mass = mass + particle_set(list(i))%atomic_kind%mass
506 : END DO
507 0 : R0 = R0/mass
508 : END SELECT
509 0 : center = R0
510 0 : DEALLOCATE (R0)
511 : END IF
512 : END IF
513 14 : CPASSERT(n_rep1 /= 0 .OR. n_rep2 /= 0)
514 : ! Potential calculation
515 56 : ALLOCATE (LocP(0:lmax, SIZE(particle_set)))
516 42 : ALLOCATE (pos(SIZE(particle_set)))
517 42 : ALLOCATE (d_el(3, SIZE(particle_set)))
518 166 : d_el = 0.0_dp
519 : ! Determining the single atomic contribution to the dielectric dipole
520 52 : DO i = 1, SIZE(particle_set)
521 152 : rvec = particle_set(i)%r - center
522 152 : r2s = DOT_PRODUCT(rvec, rvec)
523 38 : r1s = SQRT(r2s)
524 190 : LocP(:, i) = 0.0_dp
525 38 : IF (r1s /= 0.0_dp) THEN
526 170 : DO l = 0, lmax
527 : LocP(l, i) = (r1s**l*REAL(l + 1, KIND=dp)*(eps_in - eps_out))/ &
528 170 : (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 52 : pos(i) = r1s
532 : END DO
533 : ! Computes the full derivatives of the interaction energy
534 52 : DO iparticle1 = 1, SIZE(particle_set)
535 38 : ip1 = (iparticle1 - 1)*SIZE(radii)
536 152 : q1t = SUM(charges(ip1 + 1:ip1 + SIZE(radii)))
537 126 : DO iparticle2 = 1, iparticle1
538 74 : ip2 = (iparticle2 - 1)*SIZE(radii)
539 296 : q2t = SUM(charges(ip2 + 1:ip2 + SIZE(radii)))
540 : !
541 296 : r1 = particle_set(iparticle1)%r - center
542 296 : r2 = particle_set(iparticle2)%r - center
543 74 : pos1 = pos(iparticle1)
544 74 : pos2 = pos(iparticle2)
545 74 : factor1 = 0.0_dp
546 74 : factor2 = 0.0_dp
547 74 : IF (pos1*pos2 /= 0.0_dp) THEN
548 62 : pos1i = 1.0_dp/pos1
549 62 : pos2i = 1.0_dp/pos2
550 248 : dpos1 = pos1i*r1
551 248 : dpos2 = pos2i*r2
552 248 : ptcos = DOT_PRODUCT(r1, r2)
553 62 : mycos = ptcos/(pos1*pos2)
554 62 : IF (ABS(mycos) > 1.0_dp) mycos = SIGN(1.0_dp, mycos)
555 248 : dcos1 = (r2*(pos1*pos2) - pos2*dpos1*ptcos)/(pos1*pos2)**2
556 248 : dcos2 = (r1*(pos1*pos2) - pos1*dpos2*ptcos)/(pos1*pos2)**2
557 :
558 248 : DO l = 1, lmax
559 186 : lr = REAL(l, KIND=dp)
560 : factor1 = factor1 + lr*LocP(l, iparticle2)*pos1**(l - 1)*legendre(mycos, l, 0)*dpos1 &
561 744 : + 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 806 : + LocP(l, iparticle1)*pos2**l*dlegendre(mycos, l, 0)*dcos2
564 : END DO
565 : END IF
566 296 : factor1 = factor1*q1t*q2t
567 296 : factor2 = factor2*q1t*q2t
568 296 : d_el(1:3, iparticle1) = d_el(1:3, iparticle1) + 0.5_dp*factor1
569 334 : d_el(1:3, iparticle2) = d_el(1:3, iparticle2) + 0.5_dp*factor2
570 : END DO
571 : END DO
572 14 : DEALLOCATE (pos)
573 14 : DEALLOCATE (LocP)
574 14 : M => qs_env%cp_ddapc_env%Ms
575 14 : CALL cp_decpl_ddapc_forces(qs_env, M, charges, dq, d_el, particle_set)
576 14 : DEALLOCATE (d_el)
577 14 : END SUBROUTINE solvation_ddapc_force
578 :
579 : END MODULE cp_ddapc_forces
|
