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