Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 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 144 : 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, idim, ip1, ip2, iparticle1, &
91 : iparticle2, k1, k2, k3, n_rep, r1, r2, &
92 : r3
93 : INTEGER, DIMENSION(3) :: gmax, image_cell, rmax, rmin
94 : LOGICAL :: analyt
95 : REAL(KIND=dp) :: alpha, eps, fac, frac_radius, fs, &
96 : galpha, gsq, gsqi, ij_fac, q1t, q2t, &
97 : r, r2tmp, rcut, rcut2, t1, t2, tol, &
98 : tol1
99 : REAL(KIND=dp), DIMENSION(3) :: drvec, g_index, gvec, ra, rvec, svec
100 144 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: d_el, M
101 : TYPE(mp_para_env_type), POINTER :: para_env
102 :
103 144 : NULLIFY (d_el, M, para_env)
104 144 : CALL timeset(routineN, handle)
105 144 : CALL get_qs_env(qs_env, para_env=para_env)
106 144 : CPASSERT(PRESENT(charges))
107 144 : CPASSERT(ASSOCIATED(radii))
108 1440 : rcut = MIN(NORM2(cell%hmat(:, 1)), NORM2(cell%hmat(:, 2)), NORM2(cell%hmat(:, 3)))/2.0_dp
109 144 : CALL section_vals_val_get(multipole_section, "RCUT", n_rep_val=n_rep)
110 144 : IF (n_rep == 1) CALL section_vals_val_get(multipole_section, "RCUT", r_val=rcut)
111 144 : CALL section_vals_val_get(multipole_section, "EWALD_PRECISION", r_val=eps)
112 144 : CALL section_vals_val_get(multipole_section, "ANALYTICAL_GTERM", l_val=analyt)
113 : ! The spline interpolation path is only valid for orthorhombic grids.
114 144 : analyt = analyt .OR. .NOT. cell%orthorhombic .OR. .NOT. ASSOCIATED(coeff)
115 144 : rcut2 = rcut**2
116 : !
117 : ! Setting-up parameters for Ewald summation
118 : !
119 144 : eps = MIN(ABS(eps), 0.5_dp)
120 144 : tol = SQRT(ABS(LOG(eps*rcut)))
121 144 : alpha = SQRT(ABS(LOG(eps*rcut*tol)))/rcut
122 144 : galpha = 1.0_dp/(4.0_dp*alpha*alpha)
123 144 : tol1 = SQRT(-LOG(eps*rcut*(2.0_dp*tol*alpha)**2))
124 144 : IF (cell%orthorhombic) THEN
125 552 : DO idim = 1, 3
126 552 : gmax(idim) = NINT(0.25_dp + cell%hmat(idim, idim)*alpha*tol1/pi)
127 : END DO
128 : ELSE
129 24 : DO idim = 1, 3
130 78 : gmax(idim) = CEILING(alpha*tol1*NORM2(cell%hmat(:, idim))/pi)
131 : END DO
132 : END IF
133 :
134 432 : ALLOCATE (d_el(3, SIZE(particle_set)))
135 1808 : d_el = 0.0_dp
136 144 : fac = 1.e0_dp/cell%deth
137 : !
138 560 : DO iparticle1 = 1, SIZE(particle_set)
139 : !NB parallelization
140 416 : IF (MOD(iparticle1, para_env%num_pe) /= para_env%mepos) CYCLE
141 224 : ip1 = (iparticle1 - 1)*SIZE(radii)
142 896 : q1t = SUM(charges(ip1 + 1:ip1 + SIZE(radii)))
143 898 : DO iparticle2 = 1, iparticle1
144 530 : ij_fac = 1.0_dp
145 530 : IF (iparticle1 == iparticle2) ij_fac = 0.5_dp
146 :
147 530 : ip2 = (iparticle2 - 1)*SIZE(radii)
148 2120 : q2t = SUM(charges(ip2 + 1:ip2 + SIZE(radii)))
149 : !
150 : ! Real-Space Contribution
151 : !
152 2120 : rvec = particle_set(iparticle1)%r - particle_set(iparticle2)%r
153 530 : IF (iparticle1 /= iparticle2) THEN
154 306 : ra = rvec
155 1224 : r2tmp = DOT_PRODUCT(ra, ra)
156 306 : IF (r2tmp <= rcut2) THEN
157 194 : r = SQRT(r2tmp)
158 194 : t1 = erfc(alpha*r)/r
159 776 : drvec = ra/r*q1t*q2t*factor
160 194 : t2 = -2.0_dp*alpha*EXP(-alpha*alpha*r*r)/(r*rootpi) - t1/r
161 776 : d_el(1:3, iparticle1) = d_el(1:3, iparticle1) - t2*drvec
162 776 : d_el(1:3, iparticle2) = d_el(1:3, iparticle2) + t2*drvec
163 : END IF
164 : END IF
165 6890 : svec = MATMUL(cell%h_inv, rvec)
166 2120 : DO idim = 1, 3
167 6360 : frac_radius = rcut*NORM2(cell%h_inv(idim, :))
168 1590 : rmin(idim) = FLOOR(-svec(idim) - frac_radius)
169 2120 : rmax(idim) = CEILING(-svec(idim) + frac_radius)
170 : END DO
171 2200 : DO r1 = rmin(1), rmax(1)
172 7646 : DO r2 = rmin(2), rmax(2)
173 25587 : DO r3 = rmin(3), rmax(3)
174 18471 : IF ((r1 == 0) .AND. (r2 == 0) .AND. (r3 == 0)) CYCLE
175 71764 : image_cell = [r1, r2, r3]
176 340879 : ra = rvec + MATMUL(cell%hmat, REAL(image_cell, KIND=dp))
177 71764 : r2tmp = DOT_PRODUCT(ra, ra)
178 23387 : IF (r2tmp <= rcut2) THEN
179 1020 : r = SQRT(r2tmp)
180 1020 : t1 = erfc(alpha*r)/r
181 4080 : drvec = ra/r*q1t*q2t*factor*ij_fac
182 1020 : t2 = -2.0_dp*alpha*EXP(-alpha*alpha*r*r)/(r*rootpi) - t1/r
183 1020 : d_el(1, iparticle1) = d_el(1, iparticle1) - t2*drvec(1)
184 1020 : d_el(2, iparticle1) = d_el(2, iparticle1) - t2*drvec(2)
185 1020 : d_el(3, iparticle1) = d_el(3, iparticle1) - t2*drvec(3)
186 1020 : d_el(1, iparticle2) = d_el(1, iparticle2) + t2*drvec(1)
187 1020 : d_el(2, iparticle2) = d_el(2, iparticle2) + t2*drvec(2)
188 1020 : d_el(3, iparticle2) = d_el(3, iparticle2) + t2*drvec(3)
189 : END IF
190 : END DO
191 : END DO
192 : END DO
193 : !
194 : ! G-space Contribution
195 : !
196 530 : IF (analyt) THEN
197 2869 : DO k1 = 0, gmax(1)
198 72444 : DO k2 = -gmax(2), gmax(2)
199 2094321 : DO k3 = -gmax(3), gmax(3)
200 2022107 : IF (k1 == 0 .AND. k2 == 0 .AND. k3 == 0) CYCLE
201 2021877 : fs = 2.0_dp; IF (k1 == 0) fs = 1.0_dp
202 8087508 : g_index = [REAL(k1, KIND=dp), REAL(k2, KIND=dp), REAL(k3, KIND=dp)]
203 10109385 : gvec = twopi*MATMUL(TRANSPOSE(cell%h_inv), g_index)
204 8087508 : gsq = DOT_PRODUCT(gvec, gvec)
205 2021877 : gsqi = fs/gsq
206 2021877 : t1 = fac*gsqi*EXP(-galpha*gsq)
207 8087508 : t2 = -SIN(DOT_PRODUCT(gvec, rvec))*t1*q1t*q2t*factor*fourpi
208 8087508 : d_el(1:3, iparticle1) = d_el(1:3, iparticle1) - t2*gvec
209 8157083 : d_el(1:3, iparticle2) = d_el(1:3, iparticle2) + t2*gvec
210 : END DO
211 : END DO
212 : END DO
213 : ELSE
214 1200 : gvec = Eval_d_Interp_Spl3_pbc(rvec, coeff)*q1t*q2t*factor*fourpi
215 1200 : d_el(1:3, iparticle1) = d_el(1:3, iparticle1) - gvec
216 1200 : d_el(1:3, iparticle2) = d_el(1:3, iparticle2) + gvec
217 : END IF
218 946 : IF (iparticle1 /= iparticle2) THEN
219 306 : ra = rvec
220 1224 : r = SQRT(DOT_PRODUCT(ra, ra))
221 306 : t2 = -1.0_dp/(r*r)*factor
222 1224 : drvec = ra/r*q1t*q2t
223 1224 : d_el(1:3, iparticle1) = d_el(1:3, iparticle1) + t2*drvec
224 1224 : d_el(1:3, iparticle2) = d_el(1:3, iparticle2) - t2*drvec
225 : END IF
226 : END DO ! iparticle2
227 : END DO ! iparticle1
228 : !NB parallelization
229 3472 : CALL para_env%sum(d_el)
230 144 : M => qs_env%cp_ddapc_env%Md
231 144 : IF (apply_qmmm_periodic) M => qs_env%cp_ddapc_env%Mr
232 144 : CALL cp_decpl_ddapc_forces(qs_env, M, charges, dq, d_el, particle_set)
233 144 : DEALLOCATE (d_el)
234 144 : CALL timestop(handle)
235 432 : END SUBROUTINE ewald_ddapc_force
236 :
237 : ! **************************************************************************************************
238 : !> \brief Evaluation of the pulay forces due to the fitted charge density
239 : !> \param qs_env ...
240 : !> \param M ...
241 : !> \param charges ...
242 : !> \param dq ...
243 : !> \param d_el ...
244 : !> \param particle_set ...
245 : !> \par History
246 : !> 08.2005 created [tlaino]
247 : !> \author Teodoro Laino
248 : ! **************************************************************************************************
249 158 : SUBROUTINE cp_decpl_ddapc_forces(qs_env, M, charges, dq, d_el, particle_set)
250 : TYPE(qs_environment_type), POINTER :: qs_env
251 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: M
252 : REAL(KIND=dp), DIMENSION(:), POINTER :: charges
253 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: dq
254 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: d_el
255 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
256 :
257 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_decpl_ddapc_forces'
258 :
259 : INTEGER :: handle, i, iatom, ikind, j, k, natom
260 158 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
261 158 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: uv
262 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: chf
263 158 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
264 : TYPE(mp_para_env_type), POINTER :: para_env
265 158 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
266 :
267 158 : NULLIFY (para_env)
268 158 : CALL timeset(routineN, handle)
269 158 : natom = SIZE(particle_set)
270 : CALL get_qs_env(qs_env=qs_env, &
271 : atomic_kind_set=atomic_kind_set, &
272 : para_env=para_env, &
273 158 : force=force)
274 474 : ALLOCATE (chf(3, natom))
275 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
276 : atom_of_kind=atom_of_kind, &
277 158 : kind_of=kind_of)
278 :
279 474 : ALLOCATE (uv(SIZE(M, 1)))
280 47906 : uv(:) = MATMUL(M, charges)
281 612 : DO k = 1, natom
282 1974 : DO j = 1, 3
283 17278 : chf(j, k) = DOT_PRODUCT(uv, dq(:, k, j))
284 : END DO
285 : END DO
286 : !NB now that get_ddapc returns dq that's spread over nodes, must reduce chf here
287 158 : CALL para_env%sum(chf)
288 612 : DO iatom = 1, natom
289 454 : ikind = kind_of(iatom)
290 454 : i = atom_of_kind(iatom)
291 3790 : force(ikind)%ch_pulay(1:3, i) = force(ikind)%ch_pulay(1:3, i) + chf(1:3, iatom) + d_el(1:3, iatom)
292 : END DO
293 158 : DEALLOCATE (chf)
294 158 : DEALLOCATE (uv)
295 158 : CALL timestop(handle)
296 316 : END SUBROUTINE cp_decpl_ddapc_forces
297 :
298 : ! **************************************************************************************************
299 : !> \brief Evaluation of the pulay forces due to the fitted charge density
300 : !> \param qs_env ...
301 : !> \par History
302 : !> 08.2005 created [tlaino]
303 : !> \author Teodoro Laino
304 : ! **************************************************************************************************
305 128 : SUBROUTINE reset_ch_pulay(qs_env)
306 : TYPE(qs_environment_type), POINTER :: qs_env
307 :
308 : CHARACTER(len=*), PARAMETER :: routineN = 'reset_ch_pulay'
309 :
310 : INTEGER :: handle, ind
311 128 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
312 :
313 128 : CALL timeset(routineN, handle)
314 : CALL get_qs_env(qs_env=qs_env, &
315 128 : force=force)
316 348 : DO ind = 1, SIZE(force)
317 1732 : force(ind)%ch_pulay = 0.0_dp
318 : END DO
319 128 : CALL timestop(handle)
320 128 : END SUBROUTINE reset_ch_pulay
321 :
322 : ! **************************************************************************************************
323 : !> \brief computes energy and derivatives given a set of charges
324 : !> \param ddapc_restraint_control ...
325 : !> \param n_gauss ...
326 : !> \param uv derivate of energy wrt the corresponding charge entry
327 : !> \param charges current value of the charges (one number for each gaussian used)
328 : !>
329 : !> \param energy_res energy due to the restraint
330 : !> \par History
331 : !> 02.2006 [Joost VandeVondele]
332 : !> modified [Teo]
333 : !> \note
334 : !> should be easy to adapt for other specialized cases
335 : ! **************************************************************************************************
336 872 : SUBROUTINE evaluate_restraint_functional(ddapc_restraint_control, n_gauss, uv, &
337 : charges, energy_res)
338 : TYPE(ddapc_restraint_type), INTENT(INOUT) :: ddapc_restraint_control
339 : INTEGER, INTENT(in) :: n_gauss
340 : REAL(KIND=dp), DIMENSION(:) :: uv
341 : REAL(KIND=dp), DIMENSION(:), POINTER :: charges
342 : REAL(KIND=dp), INTENT(INOUT) :: energy_res
343 :
344 : INTEGER :: I, ind
345 : REAL(KIND=dp) :: dE, order_p
346 :
347 : ! order parameter (i.e. the sum of the charges of the selected atoms)
348 :
349 872 : order_p = 0.0_dp
350 1882 : DO I = 1, ddapc_restraint_control%natoms
351 1010 : ind = (ddapc_restraint_control%atoms(I) - 1)*n_gauss
352 3532 : order_p = order_p + ddapc_restraint_control%coeff(I)*SUM(charges(ind + 1:ind + n_gauss))
353 : END DO
354 872 : ddapc_restraint_control%ddapc_order_p = order_p
355 :
356 1140 : SELECT CASE (ddapc_restraint_control%functional_form)
357 : CASE (do_ddapc_restraint)
358 : ! the restraint energy
359 268 : energy_res = ddapc_restraint_control%strength*(order_p - ddapc_restraint_control%target)**2.0_dp
360 :
361 : ! derivative of the energy
362 268 : dE = 2.0_dp*ddapc_restraint_control%strength*(order_p - ddapc_restraint_control%target)
363 536 : DO I = 1, ddapc_restraint_control%natoms
364 268 : ind = (ddapc_restraint_control%atoms(I) - 1)*n_gauss
365 1340 : uv(ind + 1:ind + n_gauss) = dE*ddapc_restraint_control%coeff(I)
366 : END DO
367 : CASE (do_ddapc_constraint)
368 604 : energy_res = ddapc_restraint_control%strength*(order_p - ddapc_restraint_control%target)
369 :
370 : ! derivative of the energy
371 1346 : DO I = 1, ddapc_restraint_control%natoms
372 742 : ind = (ddapc_restraint_control%atoms(I) - 1)*n_gauss
373 2192 : uv(ind + 1:ind + n_gauss) = ddapc_restraint_control%strength*ddapc_restraint_control%coeff(I)
374 : END DO
375 :
376 : CASE DEFAULT
377 872 : CPABORT("")
378 : END SELECT
379 :
380 872 : END SUBROUTINE evaluate_restraint_functional
381 :
382 : ! **************************************************************************************************
383 : !> \brief computes derivatives for DDAPC restraint
384 : !> \param qs_env ...
385 : !> \param ddapc_restraint_control ...
386 : !> \param dq ...
387 : !> \param charges ...
388 : !> \param n_gauss ...
389 : !> \param particle_set ...
390 : !> \par History
391 : !> 02.2006 [Joost VandeVondele]
392 : !> modified [Teo]
393 : !> \note
394 : !> should be easy to adapt for other specialized cases
395 : ! **************************************************************************************************
396 36 : SUBROUTINE restraint_functional_force(qs_env, ddapc_restraint_control, dq, charges, &
397 : n_gauss, particle_set)
398 :
399 : TYPE(qs_environment_type), POINTER :: qs_env
400 : TYPE(ddapc_restraint_type), INTENT(INOUT) :: ddapc_restraint_control
401 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: dq
402 : REAL(KIND=dp), DIMENSION(:), POINTER :: charges
403 : INTEGER, INTENT(in) :: n_gauss
404 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
405 :
406 : CHARACTER(len=*), PARAMETER :: routineN = 'restraint_functional_force'
407 :
408 : INTEGER :: handle, i, iatom, ikind, j, k, natom
409 36 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
410 : REAL(kind=dp) :: dum
411 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: uv
412 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: chf
413 36 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
414 : TYPE(mp_para_env_type), POINTER :: para_env
415 36 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
416 :
417 36 : NULLIFY (para_env)
418 36 : CALL timeset(routineN, handle)
419 36 : natom = SIZE(particle_set)
420 : CALL get_qs_env(qs_env=qs_env, &
421 : atomic_kind_set=atomic_kind_set, &
422 : para_env=para_env, &
423 36 : force=force)
424 108 : ALLOCATE (chf(3, natom))
425 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
426 : atom_of_kind=atom_of_kind, &
427 36 : kind_of=kind_of)
428 :
429 108 : ALLOCATE (uv(SIZE(dq, 1)))
430 36 : uv = 0.0_dp
431 : CALL evaluate_restraint_functional(ddapc_restraint_control, n_gauss, uv, &
432 36 : charges, dum)
433 118 : DO k = 1, natom
434 364 : DO j = 1, 3
435 1318 : chf(j, k) = DOT_PRODUCT(uv, dq(:, k, j))
436 : END DO
437 : END DO
438 : !NB now that get_ddapc returns dq that's spread over nodes, must reduce chf here
439 36 : CALL para_env%sum(chf)
440 118 : DO iatom = 1, natom
441 82 : ikind = kind_of(iatom)
442 82 : i = atom_of_kind(iatom)
443 364 : force(ikind)%ch_pulay(1:3, i) = force(ikind)%ch_pulay(1:3, i) + chf(1:3, iatom)
444 : END DO
445 36 : DEALLOCATE (chf)
446 36 : DEALLOCATE (uv)
447 36 : CALL timestop(handle)
448 :
449 72 : END SUBROUTINE restraint_functional_force
450 :
451 : ! **************************************************************************************************
452 : !> \brief Evaluates the electrostatic potential due to a simple solvation model
453 : !> Spherical cavity in a dieletric medium
454 : !> \param qs_env ...
455 : !> \param solvation_section ...
456 : !> \param particle_set ...
457 : !> \param radii ...
458 : !> \param dq ...
459 : !> \param charges ...
460 : !> \par History
461 : !> 08.2006 created [tlaino]
462 : !> \author Teodoro Laino
463 : ! **************************************************************************************************
464 14 : SUBROUTINE solvation_ddapc_force(qs_env, solvation_section, particle_set, &
465 : radii, dq, charges)
466 : TYPE(qs_environment_type), POINTER :: qs_env
467 : TYPE(section_vals_type), POINTER :: solvation_section
468 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
469 : REAL(KIND=dp), DIMENSION(:), POINTER :: radii
470 : REAL(KIND=dp), DIMENSION(:, :, :), OPTIONAL, &
471 : POINTER :: dq
472 : REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: charges
473 :
474 : INTEGER :: i, ip1, ip2, iparticle1, iparticle2, l, &
475 : lmax, n_rep1, n_rep2, weight
476 14 : INTEGER, DIMENSION(:), POINTER :: list
477 : LOGICAL :: fixed_center
478 : REAL(KIND=dp) :: center(3), dcos1(3), dcos2(3), dpos1(3), dpos2(3), eps_in, eps_out, &
479 : factor1(3), factor2(3), lr, mass, mycos, pos1, pos1i, pos2, pos2i, ptcos, q1t, q2t, &
480 : r1(3), r1s, r2(3), r2s, Rs, rvec(3)
481 14 : REAL(KIND=dp), DIMENSION(:), POINTER :: pos, R0
482 14 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: d_el, LocP, M
483 :
484 14 : fixed_center = .FALSE.
485 14 : NULLIFY (d_el, M)
486 14 : eps_in = 1.0_dp
487 14 : CALL section_vals_val_get(solvation_section, "EPS_OUT", r_val=eps_out)
488 14 : CALL section_vals_val_get(solvation_section, "LMAX", i_val=lmax)
489 14 : CALL section_vals_val_get(solvation_section, "SPHERE%RADIUS", r_val=Rs)
490 14 : CALL section_vals_val_get(solvation_section, "SPHERE%CENTER%XYZ", n_rep_val=n_rep1)
491 14 : IF (n_rep1 /= 0) THEN
492 14 : CALL section_vals_val_get(solvation_section, "SPHERE%CENTER%XYZ", r_vals=R0)
493 56 : center = R0
494 : ELSE
495 : CALL section_vals_val_get(solvation_section, "SPHERE%CENTER%ATOM_LIST", &
496 0 : n_rep_val=n_rep2)
497 0 : IF (n_rep2 /= 0) THEN
498 0 : CALL section_vals_val_get(solvation_section, "SPHERE%CENTER%ATOM_LIST", i_vals=list)
499 0 : CALL section_vals_val_get(solvation_section, "SPHERE%CENTER%WEIGHT_TYPE", i_val=weight)
500 0 : ALLOCATE (R0(SIZE(list)))
501 : SELECT CASE (weight)
502 : CASE (weight_type_unit)
503 0 : R0 = 0.0_dp
504 0 : DO i = 1, SIZE(list)
505 0 : R0 = R0 + particle_set(list(i))%r
506 : END DO
507 0 : R0 = R0/REAL(SIZE(list), KIND=dp)
508 : CASE (weight_type_mass)
509 0 : R0 = 0.0_dp
510 0 : mass = 0.0_dp
511 0 : DO i = 1, SIZE(list)
512 0 : R0 = R0 + particle_set(list(i))%r*particle_set(list(i))%atomic_kind%mass
513 0 : mass = mass + particle_set(list(i))%atomic_kind%mass
514 : END DO
515 0 : R0 = R0/mass
516 : END SELECT
517 0 : center = R0
518 0 : DEALLOCATE (R0)
519 : END IF
520 : END IF
521 14 : CPASSERT(n_rep1 /= 0 .OR. n_rep2 /= 0)
522 : ! Potential calculation
523 56 : ALLOCATE (LocP(0:lmax, SIZE(particle_set)))
524 42 : ALLOCATE (pos(SIZE(particle_set)))
525 42 : ALLOCATE (d_el(3, SIZE(particle_set)))
526 166 : d_el = 0.0_dp
527 : ! Determining the single atomic contribution to the dielectric dipole
528 52 : DO i = 1, SIZE(particle_set)
529 152 : rvec = particle_set(i)%r - center
530 152 : r2s = DOT_PRODUCT(rvec, rvec)
531 38 : r1s = SQRT(r2s)
532 190 : LocP(:, i) = 0.0_dp
533 38 : IF (r1s /= 0.0_dp) THEN
534 170 : DO l = 0, lmax
535 : LocP(l, i) = (r1s**l*REAL(l + 1, KIND=dp)*(eps_in - eps_out))/ &
536 170 : (Rs**(2*l + 1)*eps_in*(REAL(l, KIND=dp)*eps_in + REAL(l + 1, KIND=dp)*eps_out))
537 : END DO
538 : END IF
539 52 : pos(i) = r1s
540 : END DO
541 : ! Computes the full derivatives of the interaction energy
542 52 : DO iparticle1 = 1, SIZE(particle_set)
543 38 : ip1 = (iparticle1 - 1)*SIZE(radii)
544 152 : q1t = SUM(charges(ip1 + 1:ip1 + SIZE(radii)))
545 126 : DO iparticle2 = 1, iparticle1
546 74 : ip2 = (iparticle2 - 1)*SIZE(radii)
547 296 : q2t = SUM(charges(ip2 + 1:ip2 + SIZE(radii)))
548 : !
549 296 : r1 = particle_set(iparticle1)%r - center
550 296 : r2 = particle_set(iparticle2)%r - center
551 74 : pos1 = pos(iparticle1)
552 74 : pos2 = pos(iparticle2)
553 74 : factor1 = 0.0_dp
554 74 : factor2 = 0.0_dp
555 74 : IF (pos1*pos2 /= 0.0_dp) THEN
556 62 : pos1i = 1.0_dp/pos1
557 62 : pos2i = 1.0_dp/pos2
558 248 : dpos1 = pos1i*r1
559 248 : dpos2 = pos2i*r2
560 248 : ptcos = DOT_PRODUCT(r1, r2)
561 62 : mycos = ptcos/(pos1*pos2)
562 62 : IF (ABS(mycos) > 1.0_dp) mycos = SIGN(1.0_dp, mycos)
563 248 : dcos1 = (r2*(pos1*pos2) - pos2*dpos1*ptcos)/(pos1*pos2)**2
564 248 : dcos2 = (r1*(pos1*pos2) - pos1*dpos2*ptcos)/(pos1*pos2)**2
565 :
566 248 : DO l = 1, lmax
567 186 : lr = REAL(l, KIND=dp)
568 : factor1 = factor1 + lr*LocP(l, iparticle2)*pos1**(l - 1)*legendre(mycos, l, 0)*dpos1 &
569 744 : + LocP(l, iparticle2)*pos1**l*dlegendre(mycos, l, 0)*dcos1
570 : factor2 = factor2 + lr*LocP(l, iparticle1)*pos2**(l - 1)*legendre(mycos, l, 0)*dpos2 &
571 806 : + LocP(l, iparticle1)*pos2**l*dlegendre(mycos, l, 0)*dcos2
572 : END DO
573 : END IF
574 296 : factor1 = factor1*q1t*q2t
575 296 : factor2 = factor2*q1t*q2t
576 296 : d_el(1:3, iparticle1) = d_el(1:3, iparticle1) + 0.5_dp*factor1
577 334 : d_el(1:3, iparticle2) = d_el(1:3, iparticle2) + 0.5_dp*factor2
578 : END DO
579 : END DO
580 14 : DEALLOCATE (pos)
581 14 : DEALLOCATE (LocP)
582 14 : M => qs_env%cp_ddapc_env%Ms
583 14 : CALL cp_decpl_ddapc_forces(qs_env, M, charges, dq, d_el, particle_set)
584 14 : DEALLOCATE (d_el)
585 14 : END SUBROUTINE solvation_ddapc_force
586 :
587 : END MODULE cp_ddapc_forces
|
