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 Methods that handle helium-solvent and helium-helium interactions
10 : !> \author Lukasz Walewski
11 : !> \date 2009-06-10
12 : ! **************************************************************************************************
13 : MODULE helium_interactions
14 :
15 : USE cp_log_handling, ONLY: cp_get_default_logger,&
16 : cp_logger_type
17 : USE helium_common, ONLY: helium_eval_chain,&
18 : helium_eval_expansion,&
19 : helium_pbc,&
20 : helium_spline
21 : USE helium_nnp, ONLY: helium_nnp_print
22 : USE helium_types, ONLY: e_id_interact,&
23 : e_id_kinetic,&
24 : e_id_potential,&
25 : e_id_thermo,&
26 : e_id_total,&
27 : e_id_virial,&
28 : helium_solvent_p_type,&
29 : helium_solvent_type
30 : USE input_constants, ONLY: helium_sampling_worm,&
31 : helium_solute_intpot_mwater,&
32 : helium_solute_intpot_nnp,&
33 : helium_solute_intpot_none
34 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
35 : section_vals_type
36 : USE kinds, ONLY: dp
37 : USE nnp_acsf, ONLY: nnp_calc_acsf
38 : USE nnp_environment_types, ONLY: nnp_type
39 : USE nnp_model, ONLY: nnp_gradients,&
40 : nnp_predict
41 : USE physcon, ONLY: angstrom,&
42 : kelvin
43 : USE pint_types, ONLY: pint_env_type
44 : USE splines_types, ONLY: spline_data_type
45 : #include "../base/base_uses.f90"
46 :
47 : IMPLICIT NONE
48 :
49 : PRIVATE
50 :
51 : LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
52 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'helium_interactions'
53 :
54 : PUBLIC :: helium_calc_energy
55 : PUBLIC :: helium_total_link_action
56 : PUBLIC :: helium_total_pair_action
57 : PUBLIC :: helium_total_inter_action
58 : PUBLIC :: helium_solute_e_f
59 : PUBLIC :: helium_bead_solute_e_f
60 : PUBLIC :: helium_intpot_scan
61 : PUBLIC :: helium_vij
62 :
63 : CONTAINS
64 :
65 : ! ***************************************************************************
66 : !> \brief Calculate the helium energy (including helium-solute interaction)
67 : !> \param helium helium environment
68 : !> \param pint_env path integral environment
69 : !> \par History
70 : !> 2009-06 moved I/O out from here [lwalewski]
71 : !> \author hforbert
72 : ! **************************************************************************************************
73 7047 : SUBROUTINE helium_calc_energy(helium, pint_env)
74 : TYPE(helium_solvent_type), INTENT(INOUT) :: helium
75 : TYPE(pint_env_type), INTENT(IN) :: pint_env
76 :
77 : INTEGER :: b, bead, i, j, n
78 7047 : INTEGER, DIMENSION(:), POINTER :: perm
79 : LOGICAL :: nperiodic
80 : REAL(KIND=dp) :: a, cell_size, en, interac, kin, pot, &
81 : rmax, rmin, vkin
82 7047 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: work2, work3
83 7047 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: work
84 : REAL(KIND=dp), DIMENSION(3) :: r
85 7047 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: pos
86 : TYPE(spline_data_type), POINTER :: e0
87 :
88 7047 : pos => helium%pos
89 7047 : perm => helium%permutation
90 7047 : e0 => helium%e0
91 7047 : cell_size = 0.5_dp*helium%cell_size
92 7047 : nperiodic = .NOT. helium%periodic
93 7047 : n = helium%atoms
94 7047 : b = helium%beads
95 7047 : en = 0.0_dp
96 7047 : pot = 0.0_dp
97 7047 : rmin = 1.0e20_dp
98 7047 : rmax = 0.0_dp
99 : ALLOCATE (work(3, helium%beads + 1), &
100 : work2(helium%beads + 1), &
101 49329 : work3(SIZE(helium%uoffdiag, 1) + 1))
102 214104 : DO i = 1, n - 1
103 3503766 : DO j = i + 1, n
104 71146494 : DO bead = 1, b
105 274716990 : work(:, bead) = pos(:, i, bead) - pos(:, j, bead)
106 : END DO
107 13158648 : work(:, b + 1) = pos(:, perm(i), 1) - pos(:, perm(j), 1)
108 3289662 : en = en + helium_eval_chain(helium, work, b + 1, work2, work3, energy=.TRUE.)
109 71353551 : DO bead = 1, b
110 67856832 : a = work2(bead)
111 67856832 : IF (a < rmin) rmin = a
112 67856832 : IF (a > rmax) rmax = a
113 71146494 : IF ((a < cell_size) .OR. nperiodic) THEN
114 63532837 : pot = pot + helium_spline(helium%vij, a)
115 : END IF
116 : END DO
117 : END DO
118 : END DO
119 7047 : DEALLOCATE (work, work2, work3)
120 7047 : pot = pot/b
121 7047 : en = en/b
122 :
123 : ! helium-solute interaction energy (all beads of all particles)
124 7047 : interac = 0.0_dp
125 7047 : IF (helium%solute_present) THEN
126 3637 : CALL helium_solute_e(pint_env, helium, interac)
127 : END IF
128 7047 : interac = interac/b
129 :
130 : !TODO:
131 7047 : vkin = 0.0_dp
132 : ! vkin = helium_virial_energy(helium)
133 :
134 7047 : kin = 0.0_dp
135 221151 : DO i = 1, n
136 856416 : r(:) = pos(:, i, b) - pos(:, perm(i), 1)
137 214104 : CALL helium_pbc(helium, r)
138 214104 : kin = kin + r(1)*r(1) + r(2)*r(2) + r(3)*r(3)
139 4414791 : DO bead = 2, b
140 16774560 : r(:) = pos(:, i, bead - 1) - pos(:, i, bead)
141 4193640 : CALL helium_pbc(helium, r)
142 4407744 : kin = kin + r(1)*r(1) + r(2)*r(2) + r(3)*r(3)
143 : END DO
144 : END DO
145 7047 : kin = 1.5_dp*n/helium%tau - 0.5*kin/(b*helium%tau**2*helium%hb2m)
146 :
147 : ! TODO: move printing somewhere else ?
148 : ! print *,"POT = ",(pot/n+helium%e_corr)*kelvin,"K"
149 : ! print *,"INTERAC = ",interac*kelvin,"K"
150 : ! print *,"RMIN= ",rmin*angstrom,"A"
151 : ! print *,"RMAX= ",rmax*angstrom,"A"
152 : ! print *,"EVIRIAL not valid!"
153 : ! print *,"ETHERMO= ",((en+kin)/n+helium%e_corr)*kelvin,"K"
154 : ! print *,"ECORR= ",helium%e_corr*kelvin,"K"
155 : !! kin = helium_total_action(helium)
156 : !! print *,"ACTION= ",kin
157 : ! print *,"WINDING#= ",helium_calc_winding(helium)
158 :
159 7047 : helium%energy_inst(e_id_potential) = pot/n + helium%e_corr
160 7047 : helium%energy_inst(e_id_kinetic) = (en - pot + kin)/n
161 7047 : helium%energy_inst(e_id_interact) = interac
162 7047 : helium%energy_inst(e_id_thermo) = (en + kin)/n + helium%e_corr
163 7047 : helium%energy_inst(e_id_virial) = vkin ! 0.0_dp at the moment
164 7047 : helium%energy_inst(e_id_total) = helium%energy_inst(e_id_thermo)
165 : ! Once vkin is properly implemented, switch to:
166 : ! helium%energy_inst(e_id_total) = (en+vkin)/n+helium%e_corr
167 :
168 14094 : END SUBROUTINE helium_calc_energy
169 :
170 : ! ***************************************************************************
171 : !> \brief Computes the total harmonic link action of the helium
172 : !> \param helium ...
173 : !> \return ...
174 : !> \date 2016-05-03
175 : !> \author Felix Uhl
176 : ! **************************************************************************************************
177 50 : REAL(KIND=dp) FUNCTION helium_total_link_action(helium) RESULT(linkaction)
178 :
179 : TYPE(helium_solvent_type), INTENT(IN) :: helium
180 :
181 : INTEGER :: iatom, ibead
182 50 : INTEGER, DIMENSION(:), POINTER :: perm
183 : REAL(KIND=dp), DIMENSION(3) :: r
184 :
185 50 : perm => helium%permutation
186 50 : linkaction = 0.0_dp
187 :
188 : ! Harmonic Link action
189 : ! (r(m-1) - r(m))**2/(4*lambda*tau)
190 800 : DO ibead = 1, helium%beads - 1
191 4550 : DO iatom = 1, helium%atoms
192 15000 : r(:) = helium%pos(:, iatom, ibead) - helium%pos(:, iatom, ibead + 1)
193 3750 : CALL helium_pbc(helium, r)
194 4500 : linkaction = linkaction + (r(1)*r(1) + r(2)*r(2) + r(3)*r(3))
195 : END DO
196 : END DO
197 300 : DO iatom = 1, helium%atoms
198 : ! choose last bead connection according to permutation table
199 1000 : r(:) = helium%pos(:, iatom, helium%beads) - helium%pos(:, perm(iatom), 1)
200 250 : CALL helium_pbc(helium, r)
201 300 : linkaction = linkaction + (r(1)*r(1) + r(2)*r(2) + r(3)*r(3))
202 : END DO
203 50 : linkaction = linkaction/(2.0_dp*helium%tau*helium%hb2m)
204 :
205 50 : END FUNCTION helium_total_link_action
206 :
207 : ! ***************************************************************************
208 : !> \brief Computes the total pair action of the helium
209 : !> \param helium ...
210 : !> \return ...
211 : !> \date 2016-05-03
212 : !> \author Felix Uhl
213 : ! **************************************************************************************************
214 50 : REAL(KIND=dp) FUNCTION helium_total_pair_action(helium) RESULT(pairaction)
215 :
216 : TYPE(helium_solvent_type), INTENT(INOUT) :: helium
217 :
218 : INTEGER :: iatom, ibead, jatom, opatom, patom
219 50 : INTEGER, DIMENSION(:), POINTER :: perm
220 50 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: work3
221 : REAL(KIND=dp), DIMENSION(3) :: r, rp
222 :
223 150 : ALLOCATE (work3(SIZE(helium%uoffdiag, 1) + 1))
224 50 : perm => helium%permutation
225 50 : pairaction = 0.0_dp
226 :
227 : ! He-He pair action
228 800 : DO ibead = 1, helium%beads - 1
229 3800 : DO iatom = 1, helium%atoms - 1
230 11250 : DO jatom = iatom + 1, helium%atoms
231 30000 : r(:) = helium%pos(:, iatom, ibead) - helium%pos(:, jatom, ibead)
232 30000 : rp(:) = helium%pos(:, iatom, ibead + 1) - helium%pos(:, jatom, ibead + 1)
233 10500 : pairaction = pairaction + helium_eval_expansion(helium, r, rp, work3)
234 : END DO
235 : END DO
236 : END DO
237 : !Ensure right permutation for pair action of last and first beads.
238 250 : DO iatom = 1, helium%atoms - 1
239 750 : DO jatom = iatom + 1, helium%atoms
240 2000 : r(:) = helium%pos(:, iatom, helium%beads) - helium%pos(:, jatom, helium%beads)
241 2000 : rp(:) = helium%pos(:, perm(iatom), 1) - helium%pos(:, perm(jatom), 1)
242 700 : pairaction = pairaction + helium_eval_expansion(helium, r, rp, work3)
243 : END DO
244 : END DO
245 :
246 : ! correct for open worm configurations
247 50 : IF (.NOT. helium%worm_is_closed) THEN
248 : ! special treatment if double bead is first bead
249 0 : iatom = helium%worm_atom_idx
250 0 : IF (helium%worm_bead_idx == 1) THEN
251 : ! patom is the atom in front of the lone head bead
252 0 : patom = helium%iperm(iatom)
253 : ! go through all atoms
254 0 : DO jatom = 1, helium%atoms
255 0 : IF (jatom == helium%worm_atom_idx) CYCLE
256 0 : opatom = helium%iperm(jatom)
257 : ! subtract pair action for closed link
258 0 : r(:) = helium%pos(:, iatom, 1) - helium%pos(:, jatom, 1)
259 0 : rp(:) = helium%pos(:, patom, helium%beads) - helium%pos(:, opatom, helium%beads)
260 0 : pairaction = pairaction - helium_eval_expansion(helium, r, rp, work3)
261 : ! and add corrected extra link
262 : ! rp stays the same
263 0 : r(:) = helium%worm_xtra_bead(:) - helium%pos(:, jatom, 1)
264 0 : pairaction = pairaction + helium_eval_expansion(helium, r, rp, work3)
265 : END DO
266 : ELSE
267 : ! bead stays constant
268 0 : ibead = helium%worm_bead_idx
269 : ! go through all atoms
270 0 : DO jatom = 1, helium%atoms
271 0 : IF (jatom == helium%worm_atom_idx) CYCLE
272 : ! subtract pair action for closed link
273 0 : r(:) = helium%pos(:, iatom, ibead) - helium%pos(:, jatom, ibead)
274 0 : rp(:) = helium%pos(:, iatom, ibead - 1) - helium%pos(:, jatom, ibead - 1)
275 0 : pairaction = pairaction - helium_eval_expansion(helium, r, rp, work3)
276 : ! and add corrected extra link
277 : ! rp stays the same
278 0 : r(:) = helium%worm_xtra_bead(:) - helium%pos(:, jatom, ibead)
279 0 : pairaction = pairaction + helium_eval_expansion(helium, r, rp, work3)
280 : END DO
281 : END IF
282 : END IF
283 50 : DEALLOCATE (work3)
284 :
285 50 : END FUNCTION helium_total_pair_action
286 :
287 : ! ***************************************************************************
288 : !> \brief Computes the total interaction of the helium with the solute
289 : !> \param pint_env ...
290 : !> \param helium ...
291 : !> \return ...
292 : !> \date 2016-05-03
293 : !> \author Felix Uhl
294 : ! **************************************************************************************************
295 50 : REAL(KIND=dp) FUNCTION helium_total_inter_action(pint_env, helium) RESULT(interaction)
296 :
297 : TYPE(pint_env_type), INTENT(IN) :: pint_env
298 : TYPE(helium_solvent_type), INTENT(IN) :: helium
299 :
300 : INTEGER :: iatom, ibead
301 : REAL(KIND=dp) :: e
302 :
303 50 : interaction = 0.0_dp
304 :
305 : ! InterAction with solute
306 50 : IF (helium%solute_present) THEN
307 850 : DO ibead = 1, helium%beads
308 4850 : DO iatom = 1, helium%atoms
309 :
310 : CALL helium_bead_solute_e_f(pint_env, helium, &
311 4000 : iatom, ibead, helium%pos(:, iatom, ibead), e)
312 4800 : interaction = interaction + e
313 : END DO
314 : END DO
315 50 : IF (helium%sampling_method == helium_sampling_worm) THEN
316 0 : IF (.NOT. helium%worm_is_closed) THEN
317 : ! subtract half of tail bead interaction again
318 : CALL helium_bead_solute_e_f(pint_env, helium, &
319 : helium%worm_atom_idx, helium%worm_bead_idx, &
320 0 : helium%pos(:, helium%worm_atom_idx, helium%worm_bead_idx), e)
321 0 : interaction = interaction - 0.5_dp*e
322 : ! add half of head bead interaction
323 : CALL helium_bead_solute_e_f(pint_env, helium, &
324 : helium%worm_atom_idx, helium%worm_bead_idx, &
325 0 : helium%worm_xtra_bead, e)
326 0 : interaction = interaction + 0.5_dp*e
327 : END IF
328 : END IF
329 : END IF
330 :
331 50 : interaction = interaction*helium%tau
332 :
333 50 : END FUNCTION helium_total_inter_action
334 :
335 : ! ***************************************************************************
336 : !> \brief Calculate general helium-solute interaction energy (and forces)
337 : !> between one helium bead and the corresponding solute time slice.
338 : !> \param pint_env path integral environment
339 : !> \param helium ...
340 : !> \param helium_part_index helium particle index
341 : !> \param helium_slice_index helium time slice index
342 : !> \param helium_r_opt explicit helium bead coordinates (optional)
343 : !> \param energy calculated energy
344 : !> \param force calculated force (if requested)
345 : !> \par History
346 : !> 2019-09 Added multiple-time striding in imag. time [cschran]
347 : !> 2023-07-23 Modified to work with NNP solute-solvent interactions [lduran]
348 : !> \author Lukasz Walewski
349 : ! **************************************************************************************************
350 5542964 : SUBROUTINE helium_bead_solute_e_f(pint_env, helium, helium_part_index, &
351 : helium_slice_index, helium_r_opt, energy, force)
352 :
353 : TYPE(pint_env_type), INTENT(IN) :: pint_env
354 : TYPE(helium_solvent_type), INTENT(IN) :: helium
355 : INTEGER, INTENT(IN) :: helium_part_index, helium_slice_index
356 : REAL(KIND=dp), DIMENSION(3), INTENT(IN), OPTIONAL :: helium_r_opt
357 : REAL(KIND=dp), INTENT(OUT) :: energy
358 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT), &
359 : OPTIONAL, POINTER :: force
360 :
361 : INTEGER :: hbeads, hi, qi, stride
362 : REAL(KIND=dp), DIMENSION(3) :: helium_r
363 5542964 : REAL(KIND=dp), DIMENSION(:), POINTER :: my_force
364 :
365 5542964 : hbeads = helium%beads
366 : ! helium bead index that is invariant wrt the rotations
367 5542964 : hi = MOD(helium_slice_index - 1 + hbeads + helium%relrot, hbeads) + 1
368 : ! solute bead index that belongs to hi helium index
369 5542964 : qi = ((hi - 1)*pint_env%p)/hbeads + 1
370 :
371 : ! coordinates of the helium bead
372 5542964 : IF (PRESENT(helium_r_opt)) THEN
373 1460066 : helium_r(:) = helium_r_opt(:)
374 : ELSE
375 16331592 : helium_r(:) = helium%pos(:, helium_part_index, helium_slice_index)
376 : END IF
377 :
378 11030756 : SELECT CASE (helium%solute_interaction)
379 :
380 : CASE (helium_solute_intpot_mwater)
381 5487792 : IF (PRESENT(force)) THEN
382 54652288 : force(:, :) = 0.0_dp
383 1171264 : my_force => force(qi, :)
384 : CALL helium_intpot_model_water( &
385 : pint_env%x(qi, :), &
386 : helium, &
387 : helium_r, &
388 : energy, &
389 : my_force &
390 1171264 : )
391 : ELSE
392 : CALL helium_intpot_model_water( &
393 : pint_env%x(qi, :), &
394 : helium, &
395 : helium_r, &
396 : energy &
397 4316528 : )
398 : END IF
399 :
400 : CASE (helium_solute_intpot_nnp)
401 55172 : IF (PRESENT(force)) THEN
402 246400 : force(:, :) = 0.0_dp
403 1600 : my_force => force(qi, :)
404 : CALL helium_intpot_nnp( &
405 : pint_env%x(qi, :), &
406 : helium, &
407 : helium_r, &
408 : energy, &
409 : my_force &
410 1600 : )
411 : ELSE
412 : CALL helium_intpot_nnp( &
413 : pint_env%x(qi, :), &
414 : helium, &
415 : helium_r, &
416 : energy &
417 53572 : )
418 : END IF
419 :
420 : CASE (helium_solute_intpot_none)
421 0 : energy = 0.0_dp
422 5542964 : IF (PRESENT(force)) THEN
423 0 : force(:, :) = 0.0_dp
424 : END IF
425 :
426 : CASE DEFAULT
427 :
428 : END SELECT
429 :
430 : ! Account for Imaginary time striding in forces:
431 5542964 : IF (PRESENT(force)) THEN
432 1172864 : IF (hbeads < pint_env%p) THEN
433 3072 : stride = pint_env%p/hbeads
434 915456 : force = force*REAL(stride, dp)
435 : END IF
436 : END IF
437 :
438 5542964 : END SUBROUTINE helium_bead_solute_e_f
439 :
440 : ! ***************************************************************************
441 : !> \brief Calculate total helium-solute interaction energy and forces.
442 : !> \param pint_env path integral environment
443 : !> \param helium ...
444 : !> \param energy calculated interaction energy
445 : !> \author Lukasz Walewski
446 : ! **************************************************************************************************
447 2647 : SUBROUTINE helium_solute_e_f(pint_env, helium, energy)
448 :
449 : TYPE(pint_env_type), INTENT(IN) :: pint_env
450 : TYPE(helium_solvent_type), INTENT(INOUT) :: helium
451 : REAL(KIND=dp), INTENT(OUT) :: energy
452 :
453 : INTEGER :: ia, ib, jb, jc
454 : REAL(KIND=dp) :: my_energy
455 2647 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: force
456 :
457 2647 : NULLIFY (force)
458 2647 : force => helium%force_inst
459 :
460 2647 : energy = 0.0_dp
461 123814 : force(:, :) = 0.0_dp
462 :
463 : ! calculate the total interaction energy and gradients between the
464 : ! solute and the helium, sum over all beads of all He particles
465 75951 : DO ia = 1, helium%atoms
466 1248815 : DO ib = 1, helium%beads
467 : CALL helium_bead_solute_e_f(pint_env, helium, ia, ib, &
468 1172864 : energy=my_energy, force=helium%rtmp_p_ndim_2d)
469 1172864 : energy = energy + my_energy
470 6042840 : DO jb = 1, pint_env%p
471 49139584 : DO jc = 1, pint_env%ndim
472 47966720 : force(jb, jc) = force(jb, jc) + helium%rtmp_p_ndim_2d(jb, jc)
473 : END DO
474 : END DO
475 : END DO
476 : END DO
477 :
478 2647 : END SUBROUTINE helium_solute_e_f
479 :
480 : ! ***************************************************************************
481 : !> \brief Calculate total helium-solute interaction energy.
482 : !> \param pint_env path integral environment
483 : !> \param helium ...
484 : !> \param energy calculated interaction energy
485 : !> \author Lukasz Walewski
486 : ! **************************************************************************************************
487 3637 : SUBROUTINE helium_solute_e(pint_env, helium, energy)
488 :
489 : TYPE(pint_env_type), INTENT(IN) :: pint_env
490 : TYPE(helium_solvent_type), INTENT(IN) :: helium
491 : REAL(KIND=dp), INTENT(OUT) :: energy
492 :
493 : INTEGER :: ia, ib
494 : REAL(KIND=dp) :: my_energy
495 :
496 3637 : energy = 0.0_dp
497 :
498 108621 : DO ia = 1, helium%atoms
499 1788365 : DO ib = 1, helium%beads
500 1679744 : CALL helium_bead_solute_e_f(pint_env, helium, ia, ib, energy=my_energy)
501 1784728 : energy = energy + my_energy
502 : END DO
503 : END DO
504 :
505 3637 : END SUBROUTINE helium_solute_e
506 :
507 : ! ***************************************************************************
508 : !> \brief Scan the helium-solute interaction energy within the periodic cell
509 : !> \param pint_env ...
510 : !> \param helium_env ...
511 : !> \date 2014-01-22
512 : !> \par History
513 : !> 2016-07-14 Modified to work with independent helium_env [cschran]
514 : !> \author Lukasz Walewski
515 : ! **************************************************************************************************
516 0 : SUBROUTINE helium_intpot_scan(pint_env, helium_env)
517 :
518 : TYPE(pint_env_type), INTENT(IN) :: pint_env
519 : TYPE(helium_solvent_p_type), DIMENSION(:), POINTER :: helium_env
520 :
521 : CHARACTER(len=*), PARAMETER :: routineN = 'helium_intpot_scan'
522 :
523 : INTEGER :: handle, ic, ix, iy, iz, k, nbin
524 : LOGICAL :: wrapped
525 : REAL(KIND=dp) :: delr, my_en, ox, oy, oz
526 : REAL(kind=dp), DIMENSION(3) :: pbc1, pbc2, pos
527 :
528 0 : CALL timeset(routineN, handle)
529 :
530 : ! Perform scan only on ionode, since this is only used to output the intpot
531 0 : IF (pint_env%logger%para_env%is_source()) THEN
532 : ! Assume ionode always to have at least one helium_env
533 0 : k = 1
534 0 : helium_env(k)%helium%rho_inst(1, :, :, :) = 0.0_dp
535 0 : nbin = helium_env(k)%helium%rho_nbin
536 0 : delr = helium_env(k)%helium%rho_delr
537 0 : helium_env(k)%helium%center(:) = 0.0_dp
538 0 : ox = helium_env(k)%helium%center(1) - helium_env(k)%helium%rho_maxr/2.0_dp
539 0 : oy = helium_env(k)%helium%center(2) - helium_env(k)%helium%rho_maxr/2.0_dp
540 0 : oz = helium_env(k)%helium%center(3) - helium_env(k)%helium%rho_maxr/2.0_dp
541 :
542 0 : DO ix = 1, nbin
543 0 : DO iy = 1, nbin
544 0 : DO iz = 1, nbin
545 :
546 : ! put the probe in the center of the current voxel
547 0 : pos(:) = [ox + (ix - 0.5_dp)*delr, oy + (iy - 0.5_dp)*delr, oz + (iz - 0.5_dp)*delr]
548 :
549 : ! calc interaction energy for the current probe position
550 0 : helium_env(k)%helium%pos(:, 1, 1) = pos(:)
551 0 : CALL helium_bead_solute_e_f(pint_env, helium_env(k)%helium, 1, 1, energy=my_en)
552 :
553 : ! check if the probe fits within the unit cell
554 0 : pbc1(:) = pos(:) - helium_env(k)%helium%center
555 0 : pbc2(:) = pbc1(:)
556 0 : CALL helium_pbc(helium_env(k)%helium, pbc2)
557 0 : wrapped = .FALSE.
558 0 : DO ic = 1, 3
559 0 : IF (ABS(pbc1(ic) - pbc2(ic)) > 10.0_dp*EPSILON(0.0_dp)) THEN
560 0 : wrapped = .TRUE.
561 : END IF
562 : END DO
563 :
564 : ! set the interaction energy value
565 0 : IF (wrapped) THEN
566 0 : helium_env(k)%helium%rho_inst(1, ix, iy, iz) = 0.0_dp
567 : ELSE
568 0 : helium_env(k)%helium%rho_inst(1, ix, iy, iz) = my_en
569 : END IF
570 :
571 : END DO
572 : END DO
573 : END DO
574 : END IF
575 :
576 0 : CALL timestop(handle)
577 0 : END SUBROUTINE helium_intpot_scan
578 :
579 : ! ***************************************************************************
580 : !> \brief Calculate model helium-solute interaction energy and forces
581 : !> between one helium bead and the corresponding solute time
582 : !> slice asuming water solute.
583 : !> \param solute_x solute positions ARR(3*NATOMS)
584 : !> to global atom indices
585 : !> \param helium only needed for helium_pbc call at the moment
586 : !> \param helium_x helium bead position ARR(3)
587 : !> \param energy calculated interaction energy
588 : !> \param force ...
589 : !> \author Felix Uhl
590 : ! **************************************************************************************************
591 5487792 : SUBROUTINE helium_intpot_model_water(solute_x, helium, helium_x, energy, force)
592 :
593 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: solute_x
594 : TYPE(helium_solvent_type), INTENT(IN) :: helium
595 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: helium_x
596 : REAL(KIND=dp), INTENT(OUT) :: energy
597 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT), &
598 : OPTIONAL, POINTER :: force
599 :
600 : INTEGER :: i, ig
601 : REAL(KIND=dp) :: d, d2, dd, ep, eps, s1, s2, sig
602 : REAL(KIND=dp), DIMENSION(3) :: dr, solute_r
603 :
604 5487792 : energy = 0.0_dp
605 5487792 : IF (PRESENT(force)) THEN
606 11712640 : force(:) = 0.0_dp
607 : END IF
608 :
609 5487792 : sig = 2.69_dp ! 1.4 Angstrom
610 5487792 : eps = 60.61e-6_dp ! 19 K
611 5487792 : s1 = 0.0_dp
612 21951168 : DO i = 1, SIZE(helium%solute_element)
613 21951168 : IF (helium%solute_element(i) == "H ") THEN
614 10975584 : ig = i - 1
615 10975584 : solute_r(1) = solute_x(3*ig + 1)
616 10975584 : solute_r(2) = solute_x(3*ig + 2)
617 10975584 : solute_r(3) = solute_x(3*ig + 3)
618 43902336 : dr(:) = solute_r(:) - helium_x(:)
619 10975584 : CALL helium_pbc(helium, dr)
620 10975584 : d2 = dr(1)*dr(1) + dr(2)*dr(2) + dr(3)*dr(3)
621 10975584 : d = SQRT(d2)
622 10975584 : dd = (sig/d)**6
623 10975584 : ep = 4.0_dp*eps*dd*(dd - 1.0_dp)
624 10975584 : s1 = s1 + ep
625 10975584 : s2 = 24.0_dp*eps*dd*(2.0_dp*dd - 1.0_dp)/d2
626 10975584 : IF (PRESENT(force)) THEN
627 2342528 : force(3*ig + 1) = force(3*ig + 1) + s2*dr(1)
628 2342528 : force(3*ig + 2) = force(3*ig + 2) + s2*dr(2)
629 2342528 : force(3*ig + 3) = force(3*ig + 3) + s2*dr(3)
630 : END IF
631 : END IF
632 : END DO ! i = 1, num_hydrogen
633 5487792 : energy = energy + s1
634 :
635 5487792 : sig = 5.01_dp ! 2.6 Angstrom
636 5487792 : eps = 104.5e-6_dp ! 33 K
637 5487792 : s1 = 0.0_dp
638 21951168 : DO i = 1, SIZE(helium%solute_element)
639 21951168 : IF (helium%solute_element(i) == "O ") THEN
640 5487792 : ig = i - 1
641 5487792 : solute_r(1) = solute_x(3*ig + 1)
642 5487792 : solute_r(2) = solute_x(3*ig + 2)
643 5487792 : solute_r(3) = solute_x(3*ig + 3)
644 21951168 : dr(:) = solute_r(:) - helium_x(:)
645 5487792 : CALL helium_pbc(helium, dr)
646 5487792 : d2 = dr(1)*dr(1) + dr(2)*dr(2) + dr(3)*dr(3)
647 5487792 : d = SQRT(d2)
648 5487792 : dd = (sig/d)**6
649 5487792 : ep = 4.0_dp*eps*dd*(dd - 1.0_dp)
650 5487792 : s1 = s1 + ep
651 5487792 : s2 = 24.0_dp*eps*dd*(2.0_dp*dd - 1.0_dp)/d2
652 5487792 : IF (PRESENT(force)) THEN
653 1171264 : force(3*ig + 1) = force(3*ig + 1) + s2*dr(1)
654 1171264 : force(3*ig + 2) = force(3*ig + 2) + s2*dr(2)
655 1171264 : force(3*ig + 3) = force(3*ig + 3) + s2*dr(3)
656 : END IF
657 : END IF
658 : END DO ! i = 1, num_chlorine
659 5487792 : energy = energy + s1
660 :
661 5487792 : END SUBROUTINE helium_intpot_model_water
662 :
663 : ! ***************************************************************************
664 : !> \brief Calculate helium-solute interaction energy and forces between one
665 : !> helium bead and the corresponding solute time slice using NNP.
666 : !> \param solute_x solute positions ARR(3*NATOMS)
667 : !> to global atom indices
668 : !> \param helium only needed for helium_pbc call at the moment
669 : !> \param helium_x helium bead position ARR(3)
670 : !> \param energy calculated interaction energy
671 : !> \param force (optional) calculated force
672 : !> \date 2023-02-22
673 : !> \author Laura Duran
674 : ! **************************************************************************************************
675 55172 : SUBROUTINE helium_intpot_nnp(solute_x, helium, helium_x, energy, force)
676 :
677 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: solute_x
678 : TYPE(helium_solvent_type), INTENT(IN) :: helium
679 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: helium_x
680 : REAL(KIND=dp), INTENT(OUT) :: energy
681 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT), &
682 : OPTIONAL, POINTER :: force
683 :
684 : INTEGER :: i, i_com, ig, ind, ind_he, j, k, m
685 : LOGICAL :: extrapolate
686 : REAL(KIND=dp) :: rsqr, rvect(3), threshold
687 55172 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: denergydsym
688 55172 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: dsymdxyz
689 : TYPE(cp_logger_type), POINTER :: logger
690 : TYPE(nnp_type), POINTER :: nnp
691 : TYPE(section_vals_type), POINTER :: print_section
692 :
693 55172 : NULLIFY (logger)
694 55172 : logger => cp_get_default_logger()
695 :
696 55172 : IF (PRESENT(force)) THEN
697 28800 : helium%nnp%myforce(:, :, :) = 0.0_dp
698 : END IF
699 :
700 55172 : extrapolate = .FALSE.
701 55172 : threshold = 0.0001d0
702 :
703 : !fill coord array
704 55172 : ig = 1
705 220688 : DO i = 1, helium%nnp%n_ele
706 165516 : IF (helium%nnp%ele(i) == 'He') THEN
707 55172 : ind_he = ig
708 220688 : DO m = 1, 3
709 220688 : helium%nnp%coord(m, ig) = helium_x(m)
710 : END DO
711 55172 : ig = ig + 1
712 : END IF
713 717236 : DO j = 1, helium%solute_atoms
714 662064 : IF (helium%nnp%ele(i) == helium%solute_element(j)) THEN
715 662064 : DO m = 1, 3
716 662064 : helium%nnp%coord(m, ig) = solute_x(3*(j - 1) + m)
717 : END DO
718 165516 : ig = ig + 1
719 : END IF
720 : END DO
721 : END DO
722 :
723 : ! check for hard core condition
724 55172 : IF (ASSOCIATED(helium%nnp_sr_cut)) THEN
725 275860 : DO i = 1, helium%nnp%num_atoms
726 220688 : IF (i == ind_he) CYCLE
727 662064 : rvect(:) = helium%nnp%coord(:, i) - helium%nnp%coord(:, ind_he)
728 165516 : CALL helium_pbc(helium, rvect)
729 165516 : rsqr = rvect(1)*rvect(1) + rvect(2)*rvect(2) + rvect(3)*rvect(3)
730 220688 : IF (rsqr < helium%nnp_sr_cut(helium%nnp%ele_ind(i))) THEN
731 0 : energy = 0.3_dp + 1.0_dp/rsqr
732 0 : IF (PRESENT(force)) THEN
733 0 : force = 0.0_dp
734 : END IF
735 0 : RETURN
736 : END IF
737 : END DO
738 : END IF
739 :
740 : ! reset flag if there's an extrapolation to report:
741 55172 : helium%nnp%output_expol = .FALSE.
742 :
743 : ! calc atomic contribution to energy and force
744 : !NOTE corresponds to nnp_force line with parallelization:
745 : !DO i = istart, istart + mecalc - 1
746 275860 : DO i = 1, helium%nnp%num_atoms
747 :
748 : !determine index of atom type
749 220688 : ind = helium%nnp%ele_ind(i)
750 :
751 : !reset input nodes and grads of ele(ind):
752 4744792 : helium%nnp%arc(ind)%layer(1)%node(:) = 0.0_dp
753 220688 : nnp => helium%nnp ! work around wrong INTENT of nnp_calc_acsf
754 220688 : IF (PRESENT(force)) THEN
755 137600 : helium%nnp%arc(ind)%layer(1)%node_grad(:) = 0.0_dp
756 25600 : ALLOCATE (dsymdxyz(3, helium%nnp%arc(ind)%n_nodes(1), helium%nnp%num_atoms))
757 19200 : ALLOCATE (denergydsym(helium%nnp%arc(ind)%n_nodes(1)))
758 2131200 : dsymdxyz(:, :, :) = 0.0_dp
759 6400 : CALL nnp_calc_acsf(nnp, i, dsymdxyz)
760 : ELSE
761 214288 : CALL nnp_calc_acsf(nnp, i)
762 : END IF
763 :
764 : ! input nodes filled, perform prediction:
765 441376 : DO i_com = 1, helium%nnp%n_committee !loop over committee members
766 : ! Predict energy
767 220688 : CALL nnp_predict(helium%nnp%arc(ind), helium%nnp, i_com)
768 220688 : helium%nnp%atomic_energy(i, i_com) = helium%nnp%arc(ind)%layer(helium%nnp%n_layer)%node(1) ! + helium%nnp%atom_energies(ind)
769 :
770 : !Gradients
771 441376 : IF (PRESENT(force)) THEN
772 :
773 137600 : denergydsym(:) = 0.0_dp
774 :
775 6400 : CALL nnp_gradients(helium%nnp%arc(ind), helium%nnp, i_com, denergydsym)
776 137600 : DO j = 1, helium%nnp%arc(ind)%n_nodes(1)
777 662400 : DO k = 1, helium%nnp%num_atoms
778 2230400 : DO m = 1, 3
779 : helium%nnp%myforce(m, k, i_com) = helium%nnp%myforce(m, k, i_com) &
780 2099200 : - denergydsym(j)*dsymdxyz(m, j, k)
781 : END DO
782 : END DO
783 : END DO
784 :
785 : END IF
786 : END DO ! end loop over committee members
787 :
788 : !deallocate memory
789 275860 : IF (PRESENT(force)) THEN
790 6400 : DEALLOCATE (denergydsym)
791 6400 : DEALLOCATE (dsymdxyz)
792 : END IF
793 :
794 : END DO ! end loop over num_atoms
795 :
796 : ! calculate energy:
797 331032 : helium%nnp%committee_energy(:) = SUM(helium%nnp%atomic_energy, 1)
798 110344 : energy = SUM(helium%nnp%committee_energy)/REAL(helium%nnp%n_committee, dp)
799 55172 : helium%nnp%nnp_potential_energy = energy
800 :
801 55172 : IF (PRESENT(force)) THEN
802 : ! bring myforce to force array
803 8000 : DO j = 1, helium%nnp%num_atoms
804 27200 : DO k = 1, 3
805 44800 : helium%nnp%committee_forces(k, j, :) = helium%nnp%myforce(k, j, :)
806 : END DO
807 : END DO
808 46400 : helium%nnp%nnp_forces(:, :) = SUM(helium%nnp%committee_forces, DIM=3)/REAL(helium%nnp%n_committee, dp)
809 : ! project out helium force entry
810 1600 : ig = 1
811 8000 : DO j = 1, helium%nnp%num_atoms
812 6400 : IF (j == ind_he) CYCLE
813 19200 : DO k = 1, 3
814 19200 : force(3*(helium%nnp%sort(ig) - 1) + k) = helium%nnp%nnp_forces(k, j)
815 : END DO
816 8000 : ig = ig + 1
817 : END DO
818 : END IF
819 :
820 : ! print properties if requested
821 55172 : print_section => section_vals_get_subs_vals(helium%nnp%nnp_input, "PRINT")
822 55172 : CALL helium_nnp_print(helium%nnp, print_section, ind_he)
823 :
824 55172 : RETURN
825 :
826 55172 : END SUBROUTINE helium_intpot_nnp
827 :
828 : ! ***************************************************************************
829 : !> \brief Helium-helium pair interaction potential.
830 : !> \param r ...
831 : !> \return ...
832 : ! **************************************************************************************************
833 1300 : ELEMENTAL FUNCTION helium_vij(r) RESULT(vij)
834 :
835 : REAL(kind=dp), INTENT(IN) :: r
836 : REAL(kind=dp) :: vij
837 :
838 : REAL(kind=dp) :: f, x, x2
839 :
840 1300 : x = angstrom*r/2.9673_dp
841 1300 : IF (x < 1.241314_dp) THEN
842 351 : x2 = 1.241314_dp/x - 1.0_dp
843 351 : f = EXP(-x2*x2)
844 : ELSE
845 : f = 1.0_dp
846 : END IF
847 1300 : x2 = 1.0_dp/(x*x)
848 : vij = 10.8_dp/kelvin*(544850.4_dp*EXP(-13.353384_dp*x) - f* &
849 1300 : ((0.1781_dp*x2 + 0.4253785_dp)*x2 + 1.3732412_dp)*x2*x2*x2)
850 1300 : END FUNCTION helium_vij
851 :
852 : #if 0
853 :
854 : ! this block is currently turned off
855 :
856 : ! ***************************************************************************
857 : !> \brief Helium-helium pair interaction potential's derivative.
858 : !> \param r ...
859 : !> \return ...
860 : ! **************************************************************************************************
861 : ELEMENTAL FUNCTION helium_d_vij(r) RESULT(dvij)
862 :
863 : REAL(kind=dp), INTENT(IN) :: r
864 : REAL(kind=dp) :: dvij
865 :
866 : REAL(kind=dp) :: f, fp, x, x2, y
867 :
868 : x = angstrom*r/2.9673_dp
869 : x = r/2.9673_dp
870 : x2 = 1.0_dp/(x*x)
871 : IF (x < 1.241314_dp) THEN
872 : y = 1.241314_dp/x - 1.0_dp
873 : f = EXP(-y*y)
874 : fp = 2.0_dp*1.241314_dp*f*y* &
875 : ((0.1781_dp*x2 + 0.4253785_dp)*x2 + 1.3732412_dp)*x2*x2*x2*x2
876 : ELSE
877 : f = 1.0_dp
878 : fp = 0.0_dp
879 : END IF
880 :
881 : dvij = angstrom*(10.8_dp/2.9673_dp)*( &
882 : (-13.353384_dp*544850.4_dp)*EXP(-13.353384_dp*x) - fp + &
883 : f*(((10.0_dp*0.1781_dp)*x2 + (8.0_dp*0.4253785_dp))*x2 + (6.0_dp*1.3732412_dp))* &
884 : x2*x2*x2/x)/(r*kelvin)
885 : END FUNCTION helium_d_vij
886 :
887 : ! **************************************************************************************************
888 : !> \brief ...
889 : !> \param helium ...
890 : !> \param n ...
891 : !> \param i ...
892 : !> \return ...
893 : ! **************************************************************************************************
894 : FUNCTION helium_atom_action(helium, n, i) RESULT(res)
895 :
896 : TYPE(helium_solvent_type), INTENT(INOUT) :: helium
897 : INTEGER, INTENT(IN) :: n, i
898 : REAL(KIND=dp) :: res
899 :
900 : INTEGER :: c, j
901 : REAL(KIND=dp) :: r(3), rp(3), s, t
902 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: work3
903 :
904 : ALLOCATE (work3(SIZE(helium%uoffdiag, 1) + 1))
905 : s = 0.0_dp
906 : t = 0.0_dp
907 : IF (n < helium%beads) THEN
908 : DO c = 1, 3
909 : r(c) = helium%pos(c, i, n) - helium%pos(c, i, n + 1)
910 : END DO
911 : CALL helium_pbc(helium, r)
912 : t = r(1)*r(1) + r(2)*r(2) + r(3)*r(3)
913 : DO j = 1, i - 1
914 : DO c = 1, 3
915 : r(c) = helium%pos(c, i, n) - helium%pos(c, j, n)
916 : rp(c) = helium%pos(c, i, n + 1) - helium%pos(c, j, n + 1)
917 : END DO
918 : s = s + helium_eval_expansion(helium, r, rp, work3)
919 : END DO
920 : DO j = i + 1, helium%atoms
921 : DO c = 1, 3
922 : r(c) = helium%pos(c, i, n) - helium%pos(c, j, n)
923 : rp(c) = helium%pos(c, i, n + 1) - helium%pos(c, j, n + 1)
924 : END DO
925 : s = s + helium_eval_expansion(helium, r, rp, work3)
926 : END DO
927 : ELSE
928 : DO c = 1, 3
929 : r(c) = helium%pos(c, i, n) - helium%pos(c, helium%permutation(i), 1)
930 : END DO
931 : CALL helium_pbc(helium, r)
932 : t = r(1)*r(1) + r(2)*r(2) + r(3)*r(3)
933 : DO j = 1, i - 1
934 : DO c = 1, 3
935 : r(c) = helium%pos(c, i, n) - helium%pos(c, j, n)
936 : rp(c) = helium%pos(c, helium%permutation(i), 1) - helium%pos(c, helium%permutation(j), 1)
937 : END DO
938 : s = s + helium_eval_expansion(helium, r, rp, work3)
939 : END DO
940 : DO j = i + 1, helium%atoms
941 : DO c = 1, 3
942 : r(c) = helium%pos(c, i, n) - helium%pos(c, j, n)
943 : rp(c) = helium%pos(c, helium%permutation(i), 1) - helium%pos(c, helium%permutation(j), 1)
944 : END DO
945 : s = s + helium_eval_expansion(helium, r, rp, work3)
946 : END DO
947 : END IF
948 : t = t/(2.0_dp*helium%tau*helium%hb2m)
949 : s = s*0.5_dp
950 : res = s + t
951 : DEALLOCATE (work3)
952 :
953 : END FUNCTION helium_atom_action
954 :
955 : ! **************************************************************************************************
956 : !> \brief ...
957 : !> \param helium ...
958 : !> \param n ...
959 : !> \return ...
960 : ! **************************************************************************************************
961 : FUNCTION helium_link_action(helium, n) RESULT(res)
962 :
963 : TYPE(helium_solvent_type), INTENT(INOUT) :: helium
964 : INTEGER, INTENT(IN) :: n
965 : REAL(KIND=dp) :: res
966 :
967 : INTEGER :: c, i, j
968 : REAL(KIND=dp) :: r(3), rp(3), s, t
969 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: work3
970 :
971 : ALLOCATE (work3(SIZE(helium%uoffdiag, 1) + 1))
972 : s = 0.0_dp
973 : t = 0.0_dp
974 : IF (n < helium%beads) THEN
975 : DO i = 1, helium%atoms
976 : DO c = 1, 3
977 : r(c) = helium%pos(c, i, n) - helium%pos(c, i, n + 1)
978 : END DO
979 : CALL helium_pbc(helium, r)
980 : t = t + r(1)*r(1) + r(2)*r(2) + r(3)*r(3)
981 : DO j = 1, i - 1
982 : DO c = 1, 3
983 : r(c) = helium%pos(c, i, n) - helium%pos(c, j, n)
984 : rp(c) = helium%pos(c, i, n + 1) - helium%pos(c, j, n + 1)
985 : END DO
986 : s = s + helium_eval_expansion(helium, r, rp, work3)
987 : END DO
988 : END DO
989 : ELSE
990 : DO i = 1, helium%atoms
991 : DO c = 1, 3
992 : r(c) = helium%pos(c, i, n) - helium%pos(c, helium%permutation(i), 1)
993 : END DO
994 : CALL helium_pbc(helium, r)
995 : t = t + r(1)*r(1) + r(2)*r(2) + r(3)*r(3)
996 : DO j = 1, i - 1
997 : DO c = 1, 3
998 : r(c) = helium%pos(c, i, n) - helium%pos(c, j, n)
999 : rp(c) = helium%pos(c, helium%permutation(i), 1) - helium%pos(c, helium%permutation(j), 1)
1000 : END DO
1001 : s = s + helium_eval_expansion(helium, r, rp, work3)
1002 : END DO
1003 : END DO
1004 : END IF
1005 : t = t/(2.0_dp*helium%tau*helium%hb2m)
1006 : res = s + t
1007 : DEALLOCATE (work3)
1008 :
1009 : END FUNCTION helium_link_action
1010 :
1011 : ! **************************************************************************************************
1012 : !> \brief ...
1013 : !> \param helium ...
1014 : !> \return ...
1015 : ! **************************************************************************************************
1016 : FUNCTION helium_total_action(helium) RESULT(res)
1017 :
1018 : TYPE(helium_solvent_type), INTENT(INOUT) :: helium
1019 : REAL(KIND=dp) :: res
1020 :
1021 : INTEGER :: i
1022 : REAL(KIND=dp) :: s
1023 :
1024 : s = 0.0_dp
1025 : DO i = 1, helium%beads
1026 : s = s + helium_link_action(helium, i)
1027 : END DO
1028 : res = s
1029 :
1030 : END FUNCTION helium_total_action
1031 :
1032 : ! **************************************************************************************************
1033 : !> \brief ...
1034 : !> \param helium ...
1035 : !> \param part ...
1036 : !> \param ref_bead ...
1037 : !> \param delta_bead ...
1038 : !> \param d ...
1039 : ! **************************************************************************************************
1040 : SUBROUTINE helium_delta_pos(helium, part, ref_bead, delta_bead, d)
1041 :
1042 : TYPE(helium_solvent_type), INTENT(INOUT) :: helium
1043 : INTEGER, INTENT(IN) :: part, ref_bead, delta_bead
1044 : REAL(KIND=dp), DIMENSION(3), INTENT(OUT) :: d
1045 :
1046 : INTEGER :: b, bead, db, nbead, np, p
1047 : REAL(KIND=dp), DIMENSION(3) :: r
1048 :
1049 : b = helium%beads
1050 :
1051 : d(:) = 0.0_dp
1052 : IF (delta_bead > 0) THEN
1053 : bead = ref_bead
1054 : p = part
1055 : db = delta_bead
1056 : DO
1057 : IF (db < 1) EXIT
1058 : nbead = bead + 1
1059 : np = p
1060 : IF (nbead > b) THEN
1061 : nbead = nbead - b
1062 : np = helium%permutation(np)
1063 : END IF
1064 : r(:) = helium%pos(:, p, bead) - helium%pos(:, np, nbead)
1065 : CALL helium_pbc(helium, r)
1066 : d(:) = d(:) + r(:)
1067 : bead = nbead
1068 : p = np
1069 : db = db - 1
1070 : END DO
1071 : ELSEIF (delta_bead < 0) THEN
1072 : bead = ref_bead
1073 : p = part
1074 : db = delta_bead
1075 : DO
1076 : IF (db >= 0) EXIT
1077 : nbead = bead - 1
1078 : np = p
1079 : IF (nbead < 1) THEN
1080 : nbead = nbead + b
1081 : np = helium%iperm(np)
1082 : END IF
1083 : r(:) = helium%pos(:, p, bead) - helium%pos(:, np, nbead)
1084 : CALL helium_pbc(helium, r)
1085 : d(:) = d(:) + r(:)
1086 : bead = nbead
1087 : p = np
1088 : db = db + 1
1089 : END DO
1090 : END IF
1091 : END SUBROUTINE helium_delta_pos
1092 :
1093 : #endif
1094 :
1095 : END MODULE helium_interactions
|
