(git:4dc5213)
Loading...
Searching...
No Matches
helium_interactions.F
Go to the documentation of this file.
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! **************************************************************************************************
13MODULE 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
65CONTAINS
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 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 INTEGER, DIMENSION(:), POINTER :: perm
81 LOGICAL :: nperiodic
82 REAL(kind=dp) :: a, cell_size, en, interac, kin, pot, &
83 rmax, rmin, vkin
84 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: work2, work3
85 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: work
86 REAL(kind=dp), DIMENSION(3) :: r
87 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: pos
88 TYPE(spline_data_type), POINTER :: e0
89
90 pos => helium%pos
91 perm => helium%permutation
92 e0 => helium%e0
93 cell_size = 0.5_dp*helium%cell_size
94 nperiodic = .NOT. helium%periodic
95 n = helium%atoms
96 b = helium%beads
97 en = 0.0_dp
98 pot = 0.0_dp
99 rmin = 1.0e20_dp
100 rmax = 0.0_dp
101 ALLOCATE (work(3, helium%beads + 1), &
102 work2(helium%beads + 1), &
103 work3(SIZE(helium%uoffdiag, 1) + 1))
104 DO i = 1, n - 1
105 DO j = i + 1, n
106 DO bead = 1, b
107 work(:, bead) = pos(:, i, bead) - pos(:, j, bead)
108 END DO
109 work(:, b + 1) = pos(:, perm(i), 1) - pos(:, perm(j), 1)
110 en = en + helium_eval_chain(helium, work, b + 1, work2, work3, energy=.true.)
111 DO bead = 1, b
112 a = work2(bead)
113 IF (a < rmin) rmin = a
114 IF (a > rmax) rmax = a
115 IF ((a < cell_size) .OR. nperiodic) THEN
116 pot = pot + helium_spline(helium%vij, a)
117 END IF
118 END DO
119 END DO
120 END DO
121 DEALLOCATE (work, work2, work3)
122 pot = pot/b
123 en = en/b
124
125 ! helium-solute interaction energy (all beads of all particles)
126 interac = 0.0_dp
127 IF (helium%solute_present) THEN
128 CALL helium_solute_e(pint_env, helium, interac)
129 END IF
130 interac = interac/b
131
132!TODO:
133 vkin = 0.0_dp
134! vkin = helium_virial_energy(helium)
135
136 kin = 0.0_dp
137 DO i = 1, n
138 r(:) = pos(:, i, b) - pos(:, perm(i), 1)
139 CALL helium_pbc(helium, r)
140 kin = kin + r(1)*r(1) + r(2)*r(2) + r(3)*r(3)
141 DO bead = 2, b
142 r(:) = pos(:, i, bead - 1) - pos(:, i, bead)
143 CALL helium_pbc(helium, r)
144 kin = kin + r(1)*r(1) + r(2)*r(2) + r(3)*r(3)
145 END DO
146 END DO
147 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 helium%energy_inst(e_id_potential) = pot/n + helium%e_corr
162 helium%energy_inst(e_id_kinetic) = (en - pot + kin)/n
163 helium%energy_inst(e_id_interact) = interac
164 helium%energy_inst(e_id_thermo) = (en + kin)/n + helium%e_corr
165 helium%energy_inst(e_id_virial) = vkin ! 0.0_dp at the moment
166 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 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 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 INTEGER, DIMENSION(:), POINTER :: perm
185 REAL(kind=dp), DIMENSION(3) :: r
186
187 perm => helium%permutation
188 linkaction = 0.0_dp
189
190 ! Harmonic Link action
191 ! (r(m-1) - r(m))**2/(4*lambda*tau)
192 DO ibead = 1, helium%beads - 1
193 DO iatom = 1, helium%atoms
194 r(:) = helium%pos(:, iatom, ibead) - helium%pos(:, iatom, ibead + 1)
195 CALL helium_pbc(helium, r)
196 linkaction = linkaction + (r(1)*r(1) + r(2)*r(2) + r(3)*r(3))
197 END DO
198 END DO
199 DO iatom = 1, helium%atoms
200 ! choose last bead connection according to permutation table
201 r(:) = helium%pos(:, iatom, helium%beads) - helium%pos(:, perm(iatom), 1)
202 CALL helium_pbc(helium, r)
203 linkaction = linkaction + (r(1)*r(1) + r(2)*r(2) + r(3)*r(3))
204 END DO
205 linkaction = linkaction/(2.0_dp*helium%tau*helium%hb2m)
206
207 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 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 INTEGER, DIMENSION(:), POINTER :: perm
222 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: work3
223 REAL(kind=dp), DIMENSION(3) :: r, rp
224
225 ALLOCATE (work3(SIZE(helium%uoffdiag, 1) + 1))
226 perm => helium%permutation
227 pairaction = 0.0_dp
228
229 ! He-He pair action
230 DO ibead = 1, helium%beads - 1
231 DO iatom = 1, helium%atoms - 1
232 DO jatom = iatom + 1, helium%atoms
233 r(:) = helium%pos(:, iatom, ibead) - helium%pos(:, jatom, ibead)
234 rp(:) = helium%pos(:, iatom, ibead + 1) - helium%pos(:, jatom, ibead + 1)
235 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 DO iatom = 1, helium%atoms - 1
241 DO jatom = iatom + 1, helium%atoms
242 r(:) = helium%pos(:, iatom, helium%beads) - helium%pos(:, jatom, helium%beads)
243 rp(:) = helium%pos(:, perm(iatom), 1) - helium%pos(:, perm(jatom), 1)
244 pairaction = pairaction + helium_eval_expansion(helium, r, rp, work3)
245 END DO
246 END DO
247
248 ! correct for open worm configurations
249 IF (.NOT. helium%worm_is_closed) THEN
250 ! special treatment if double bead is first bead
251 iatom = helium%worm_atom_idx
252 IF (helium%worm_bead_idx == 1) THEN
253 ! patom is the atom in front of the lone head bead
254 patom = helium%iperm(iatom)
255 ! go through all atoms
256 DO jatom = 1, helium%atoms
257 IF (jatom == helium%worm_atom_idx) cycle
258 opatom = helium%iperm(jatom)
259 ! subtract pair action for closed link
260 r(:) = helium%pos(:, iatom, 1) - helium%pos(:, jatom, 1)
261 rp(:) = helium%pos(:, patom, helium%beads) - helium%pos(:, opatom, helium%beads)
262 pairaction = pairaction - helium_eval_expansion(helium, r, rp, work3)
263 ! and add corrected extra link
264 ! rp stays the same
265 r(:) = helium%worm_xtra_bead(:) - helium%pos(:, jatom, 1)
266 pairaction = pairaction + helium_eval_expansion(helium, r, rp, work3)
267 END DO
268 ELSE
269 ! bead stays constant
270 ibead = helium%worm_bead_idx
271 ! go through all atoms
272 DO jatom = 1, helium%atoms
273 IF (jatom == helium%worm_atom_idx) cycle
274 ! subtract pair action for closed link
275 r(:) = helium%pos(:, iatom, ibead) - helium%pos(:, jatom, ibead)
276 rp(:) = helium%pos(:, iatom, ibead - 1) - helium%pos(:, jatom, ibead - 1)
277 pairaction = pairaction - helium_eval_expansion(helium, r, rp, work3)
278 ! and add corrected extra link
279 ! rp stays the same
280 r(:) = helium%worm_xtra_bead(:) - helium%pos(:, jatom, ibead)
281 pairaction = pairaction + helium_eval_expansion(helium, r, rp, work3)
282 END DO
283 END IF
284 END IF
285 DEALLOCATE (work3)
286
287 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 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 interaction = 0.0_dp
306
307 ! InterAction with solute
308 IF (helium%solute_present) THEN
309 DO ibead = 1, helium%beads
310 DO iatom = 1, helium%atoms
311
312 CALL helium_bead_solute_e_f(pint_env, helium, &
313 iatom, ibead, helium%pos(:, iatom, ibead), e)
314 interaction = interaction + e
315 END DO
316 END DO
317 IF (helium%sampling_method == helium_sampling_worm) THEN
318 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 helium%pos(:, helium%worm_atom_idx, helium%worm_bead_idx), e)
323 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 helium%worm_xtra_bead, e)
328 interaction = interaction + 0.5_dp*e
329 END IF
330 END IF
331 END IF
332
333 interaction = interaction*helium%tau
334
335 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 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 REAL(kind=dp), DIMENSION(:), POINTER :: my_force
366
367 hbeads = helium%beads
368 ! helium bead index that is invariant wrt the rotations
369 hi = mod(helium_slice_index - 1 + hbeads + helium%relrot, hbeads) + 1
370 ! solute bead index that belongs to hi helium index
371 qi = ((hi - 1)*pint_env%p)/hbeads + 1
372
373 ! coordinates of the helium bead
374 IF (PRESENT(helium_r_opt)) THEN
375 helium_r(:) = helium_r_opt(:)
376 ELSE
377 helium_r(:) = helium%pos(:, helium_part_index, helium_slice_index)
378 END IF
379
380 SELECT CASE (helium%solute_interaction)
381
382 CASE (helium_solute_intpot_mwater)
383 IF (PRESENT(force)) THEN
384 force(:, :) = 0.0_dp
385 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 )
393 ELSE
394 CALL helium_intpot_model_water( &
395 pint_env%x(qi, :), &
396 helium, &
397 helium_r, &
398 energy &
399 )
400 END IF
401
402 CASE (helium_solute_intpot_nnp)
403 IF (PRESENT(force)) THEN
404 force(:, :) = 0.0_dp
405 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 )
413 ELSE
414 CALL helium_intpot_nnp( &
415 pint_env%x(qi, :), &
416 helium, &
417 helium_r, &
418 energy &
419 )
420 END IF
421
422 CASE (helium_solute_intpot_none)
423 energy = 0.0_dp
424 IF (PRESENT(force)) THEN
425 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 IF (PRESENT(force)) THEN
434 IF (hbeads < pint_env%p) THEN
435 stride = pint_env%p/hbeads
436 force = force*real(stride, dp)
437 END IF
438 END IF
439
440 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 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 REAL(kind=dp), DIMENSION(:, :), POINTER :: force
458
459 NULLIFY (force)
460 force => helium%force_inst
461
462 energy = 0.0_dp
463 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 DO ia = 1, helium%atoms
468 DO ib = 1, helium%beads
469 CALL helium_bead_solute_e_f(pint_env, helium, ia, ib, &
470 energy=my_energy, force=helium%rtmp_p_ndim_2d)
471 energy = energy + my_energy
472 DO jb = 1, pint_env%p
473 DO jc = 1, pint_env%ndim
474 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 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 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 energy = 0.0_dp
499
500 DO ia = 1, helium%atoms
501 DO ib = 1, helium%beads
502 CALL helium_bead_solute_e_f(pint_env, helium, ia, ib, energy=my_energy)
503 energy = energy + my_energy
504 END DO
505 END DO
506
507 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 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 CALL timeset(routinen, handle)
531
532 ! Perform scan only on ionode, since this is only used to output the intpot
533 IF (pint_env%logger%para_env%is_source()) THEN
534 ! Assume ionode always to have at least one helium_env
535 k = 1
536 helium_env(k)%helium%rho_inst(1, :, :, :) = 0.0_dp
537 nbin = helium_env(k)%helium%rho_nbin
538 delr = helium_env(k)%helium%rho_delr
539 helium_env(k)%helium%center(:) = 0.0_dp
540 ox = helium_env(k)%helium%center(1) - helium_env(k)%helium%rho_maxr/2.0_dp
541 oy = helium_env(k)%helium%center(2) - helium_env(k)%helium%rho_maxr/2.0_dp
542 oz = helium_env(k)%helium%center(3) - helium_env(k)%helium%rho_maxr/2.0_dp
543
544 DO ix = 1, nbin
545 DO iy = 1, nbin
546 DO iz = 1, nbin
547
548 ! put the probe in the center of the current voxel
549 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 helium_env(k)%helium%pos(:, 1, 1) = pos(:)
553 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 pbc1(:) = pos(:) - helium_env(k)%helium%center
557 pbc2(:) = pbc1(:)
558 CALL helium_pbc(helium_env(k)%helium, pbc2)
559 wrapped = .false.
560 DO ic = 1, 3
561 IF (abs(pbc1(ic) - pbc2(ic)) > 10.0_dp*epsilon(0.0_dp)) THEN
562 wrapped = .true.
563 END IF
564 END DO
565
566 ! set the interaction energy value
567 IF (wrapped) THEN
568 helium_env(k)%helium%rho_inst(1, ix, iy, iz) = 0.0_dp
569 ELSE
570 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 CALL timestop(handle)
579 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 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 energy = 0.0_dp
607 IF (PRESENT(force)) THEN
608 force(:) = 0.0_dp
609 END IF
610
611 sig = 2.69_dp ! 1.4 Angstrom
612 eps = 60.61e-6_dp ! 19 K
613 s1 = 0.0_dp
614 DO i = 1, SIZE(helium%solute_element)
615 IF (helium%solute_element(i) == "H ") THEN
616 ig = i - 1
617 solute_r(1) = solute_x(3*ig + 1)
618 solute_r(2) = solute_x(3*ig + 2)
619 solute_r(3) = solute_x(3*ig + 3)
620 dr(:) = solute_r(:) - helium_x(:)
621 CALL helium_pbc(helium, dr)
622 d2 = dr(1)*dr(1) + dr(2)*dr(2) + dr(3)*dr(3)
623 d = sqrt(d2)
624 dd = (sig/d)**6
625 ep = 4.0_dp*eps*dd*(dd - 1.0_dp)
626 s1 = s1 + ep
627 s2 = 24.0_dp*eps*dd*(2.0_dp*dd - 1.0_dp)/d2
628 IF (PRESENT(force)) THEN
629 force(3*ig + 1) = force(3*ig + 1) + s2*dr(1)
630 force(3*ig + 2) = force(3*ig + 2) + s2*dr(2)
631 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 energy = energy + s1
636
637 sig = 5.01_dp ! 2.6 Angstrom
638 eps = 104.5e-6_dp ! 33 K
639 s1 = 0.0_dp
640 DO i = 1, SIZE(helium%solute_element)
641 IF (helium%solute_element(i) == "O ") THEN
642 ig = i - 1
643 solute_r(1) = solute_x(3*ig + 1)
644 solute_r(2) = solute_x(3*ig + 2)
645 solute_r(3) = solute_x(3*ig + 3)
646 dr(:) = solute_r(:) - helium_x(:)
647 CALL helium_pbc(helium, dr)
648 d2 = dr(1)*dr(1) + dr(2)*dr(2) + dr(3)*dr(3)
649 d = sqrt(d2)
650 dd = (sig/d)**6
651 ep = 4.0_dp*eps*dd*(dd - 1.0_dp)
652 s1 = s1 + ep
653 s2 = 24.0_dp*eps*dd*(2.0_dp*dd - 1.0_dp)/d2
654 IF (PRESENT(force)) THEN
655 force(3*ig + 1) = force(3*ig + 1) + s2*dr(1)
656 force(3*ig + 2) = force(3*ig + 2) + s2*dr(2)
657 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 energy = energy + s1
662
663 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 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 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 NULLIFY (logger)
694 logger => cp_get_default_logger()
695
696 IF (PRESENT(force)) THEN
697 helium%nnp%myforce(:, :, :) = 0.0_dp
698 END IF
699
700 !fill coord array
701 ig = 1
702 DO i = 1, helium%nnp%n_ele
703 IF (helium%nnp%ele(i) == 'He') THEN
704 ind_he = ig
705 DO m = 1, 3
706 helium%nnp%coord(m, ig) = helium_x(m)
707 END DO
708 ig = ig + 1
709 END IF
710 DO j = 1, helium%solute_atoms
711 IF (helium%nnp%ele(i) == helium%solute_element(j)) THEN
712 DO m = 1, 3
713 helium%nnp%coord(m, ig) = solute_x(3*(j - 1) + m)
714 END DO
715 ig = ig + 1
716 END IF
717 END DO
718 END DO
719
720 ! check for hard core condition
721 IF (ASSOCIATED(helium%nnp_sr_cut)) THEN
722 DO i = 1, helium%nnp%num_atoms
723 IF (i == ind_he) cycle
724 rvect(:) = helium%nnp%coord(:, i) - helium%nnp%coord(:, ind_he)
725 CALL helium_pbc(helium, rvect)
726 rsqr = rvect(1)*rvect(1) + rvect(2)*rvect(2) + rvect(3)*rvect(3)
727 IF (rsqr < helium%nnp_sr_cut(helium%nnp%ele_ind(i))) THEN
728 energy = 0.3_dp + 1.0_dp/rsqr
729 IF (PRESENT(force)) THEN
730 force = 0.0_dp
731 END IF
732 RETURN
733 END IF
734 END DO
735 END IF
736
737 ! reset flag if there's an extrapolation to report:
738 helium%nnp%output_expol = .false.
739 nnp => helium%nnp
740 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 DO i = 1, helium%nnp%num_atoms
746
747 !determine index of atom type
748 ind = helium%nnp%ele_ind(i)
749
750 !reset input nodes and grads of ele(ind):
751 helium%nnp%arc(ind)%layer(1)%node(:) = 0.0_dp
752 IF (PRESENT(force)) THEN
753 helium%nnp%arc(ind)%layer(1)%node_grad(:) = 0.0_dp
754 ALLOCATE (denergydsym(helium%nnp%arc(ind)%n_nodes(1)))
755 ! this path is serial and fills the shared nnp%arc(ind) directly
756 CALL nnp_calc_acsf(nnp, i, .true., nnp%arc(ind))
757 ELSE
758 CALL nnp_calc_acsf(nnp, i, .false., nnp%arc(ind))
759 END IF
760
761 ! input nodes filled, perform prediction:
762 DO i_com = 1, helium%nnp%n_committee !loop over committee members
763 ! Predict energy
764 CALL nnp_predict(helium%nnp%arc(ind), helium%nnp, i_com)
765 helium%nnp%atomic_energy(i, i_com) = helium%nnp%arc(ind)%layer(helium%nnp%n_layer)%node(1)
766
767 !Gradients
768 IF (PRESENT(force)) THEN
769
770 denergydsym(:) = 0.0_dp
771 CALL nnp_gradients(helium%nnp%arc(ind), helium%nnp, i_com, denergydsym)
772
773 ! Per-element workspace dGdr scatter; shared with the main NNP
774 ! force path (see nnp_scatter_dgdr_to_forces in nnp_force.F).
775 CALL nnp_scatter_dgdr_to_forces(nnp, ind, i, denergydsym, &
776 helium%nnp%myforce(:, :, i_com))
777
778 END IF
779 END DO ! end loop over committee members
780
781 !deallocate memory
782 IF (PRESENT(force)) THEN
783 DEALLOCATE (denergydsym)
784 END IF
785
786 END DO ! end loop over num_atoms
787
788 ! calculate energy:
789 helium%nnp%committee_energy(:) = sum(helium%nnp%atomic_energy, 1)
790 energy = sum(helium%nnp%committee_energy)/real(helium%nnp%n_committee, dp)
791 helium%nnp%nnp_potential_energy = energy
792
793 IF (PRESENT(force)) THEN
794 ! bring myforce to force array
795 DO j = 1, helium%nnp%num_atoms
796 DO k = 1, 3
797 helium%nnp%committee_forces(k, j, :) = helium%nnp%myforce(k, j, :)
798 END DO
799 END DO
800 helium%nnp%nnp_forces(:, :) = sum(helium%nnp%committee_forces, dim=3)/real(helium%nnp%n_committee, dp)
801 ! project out helium force entry
802 ig = 1
803 DO j = 1, helium%nnp%num_atoms
804 IF (j == ind_he) cycle
805 DO k = 1, 3
806 force(3*(helium%nnp%sort(ig) - 1) + k) = helium%nnp%nnp_forces(k, j)
807 END DO
808 ig = ig + 1
809 END DO
810 END IF
811
812 ! print properties if requested
813 print_section => section_vals_get_subs_vals(helium%nnp%nnp_input, "PRINT")
814 CALL helium_nnp_print(helium%nnp, print_section, ind_he)
815
816 RETURN
817
818 END SUBROUTINE helium_intpot_nnp
819
820! ***************************************************************************
821!> \brief Helium-helium pair interaction potential.
822!> \param r ...
823!> \return ...
824! **************************************************************************************************
825 ELEMENTAL FUNCTION helium_vij(r) RESULT(vij)
826
827 REAL(kind=dp), INTENT(IN) :: r
828 REAL(kind=dp) :: vij
829
830 REAL(kind=dp) :: f, x, x2
831
832 x = angstrom*r/2.9673_dp
833 IF (x < 1.241314_dp) THEN
834 x2 = 1.241314_dp/x - 1.0_dp
835 f = exp(-x2*x2)
836 ELSE
837 f = 1.0_dp
838 END IF
839 x2 = 1.0_dp/(x*x)
840 vij = 10.8_dp/kelvin*(544850.4_dp*exp(-13.353384_dp*x) - f* &
841 ((0.1781_dp*x2 + 0.4253785_dp)*x2 + 1.3732412_dp)*x2*x2*x2)
842 END FUNCTION helium_vij
843
844#if 0
845
846 ! this block is currently turned off
847
848! ***************************************************************************
849!> \brief Helium-helium pair interaction potential's derivative.
850!> \param r ...
851!> \return ...
852! **************************************************************************************************
853 ELEMENTAL FUNCTION helium_d_vij(r) RESULT(dvij)
854
855 REAL(kind=dp), INTENT(IN) :: r
856 REAL(kind=dp) :: dvij
857
858 REAL(kind=dp) :: f, fp, x, x2, y
859
860 x = angstrom*r/2.9673_dp
861 x = r/2.9673_dp
862 x2 = 1.0_dp/(x*x)
863 IF (x < 1.241314_dp) THEN
864 y = 1.241314_dp/x - 1.0_dp
865 f = exp(-y*y)
866 fp = 2.0_dp*1.241314_dp*f*y* &
867 ((0.1781_dp*x2 + 0.4253785_dp)*x2 + 1.3732412_dp)*x2*x2*x2*x2
868 ELSE
869 f = 1.0_dp
870 fp = 0.0_dp
871 END IF
872
873 dvij = angstrom*(10.8_dp/2.9673_dp)*( &
874 (-13.353384_dp*544850.4_dp)*exp(-13.353384_dp*x) - fp + &
875 f*(((10.0_dp*0.1781_dp)*x2 + (8.0_dp*0.4253785_dp))*x2 + (6.0_dp*1.3732412_dp))* &
876 x2*x2*x2/x)/(r*kelvin)
877 END FUNCTION helium_d_vij
878
879! **************************************************************************************************
880!> \brief ...
881!> \param helium ...
882!> \param n ...
883!> \param i ...
884!> \return ...
885! **************************************************************************************************
886 FUNCTION helium_atom_action(helium, n, i) RESULT(res)
887
888 TYPE(helium_solvent_type), INTENT(INOUT) :: helium
889 INTEGER, INTENT(IN) :: n, i
890 REAL(kind=dp) :: res
891
892 INTEGER :: c, j
893 REAL(kind=dp) :: r(3), rp(3), s, t
894 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: work3
895
896 ALLOCATE (work3(SIZE(helium%uoffdiag, 1) + 1))
897 s = 0.0_dp
898 t = 0.0_dp
899 IF (n < helium%beads) THEN
900 DO c = 1, 3
901 r(c) = helium%pos(c, i, n) - helium%pos(c, i, n + 1)
902 END DO
903 CALL helium_pbc(helium, r)
904 t = r(1)*r(1) + r(2)*r(2) + r(3)*r(3)
905 DO j = 1, i - 1
906 DO c = 1, 3
907 r(c) = helium%pos(c, i, n) - helium%pos(c, j, n)
908 rp(c) = helium%pos(c, i, n + 1) - helium%pos(c, j, n + 1)
909 END DO
910 s = s + helium_eval_expansion(helium, r, rp, work3)
911 END DO
912 DO j = i + 1, helium%atoms
913 DO c = 1, 3
914 r(c) = helium%pos(c, i, n) - helium%pos(c, j, n)
915 rp(c) = helium%pos(c, i, n + 1) - helium%pos(c, j, n + 1)
916 END DO
917 s = s + helium_eval_expansion(helium, r, rp, work3)
918 END DO
919 ELSE
920 DO c = 1, 3
921 r(c) = helium%pos(c, i, n) - helium%pos(c, helium%permutation(i), 1)
922 END DO
923 CALL helium_pbc(helium, r)
924 t = r(1)*r(1) + r(2)*r(2) + r(3)*r(3)
925 DO j = 1, i - 1
926 DO c = 1, 3
927 r(c) = helium%pos(c, i, n) - helium%pos(c, j, n)
928 rp(c) = helium%pos(c, helium%permutation(i), 1) - helium%pos(c, helium%permutation(j), 1)
929 END DO
930 s = s + helium_eval_expansion(helium, r, rp, work3)
931 END DO
932 DO j = i + 1, helium%atoms
933 DO c = 1, 3
934 r(c) = helium%pos(c, i, n) - helium%pos(c, j, n)
935 rp(c) = helium%pos(c, helium%permutation(i), 1) - helium%pos(c, helium%permutation(j), 1)
936 END DO
937 s = s + helium_eval_expansion(helium, r, rp, work3)
938 END DO
939 END IF
940 t = t/(2.0_dp*helium%tau*helium%hb2m)
941 s = s*0.5_dp
942 res = s + t
943 DEALLOCATE (work3)
944
945 END FUNCTION helium_atom_action
946
947! **************************************************************************************************
948!> \brief ...
949!> \param helium ...
950!> \param n ...
951!> \return ...
952! **************************************************************************************************
953 FUNCTION helium_link_action(helium, n) RESULT(res)
954
955 TYPE(helium_solvent_type), INTENT(INOUT) :: helium
956 INTEGER, INTENT(IN) :: n
957 REAL(kind=dp) :: res
958
959 INTEGER :: c, i, j
960 REAL(kind=dp) :: r(3), rp(3), s, t
961 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: work3
962
963 ALLOCATE (work3(SIZE(helium%uoffdiag, 1) + 1))
964 s = 0.0_dp
965 t = 0.0_dp
966 IF (n < helium%beads) THEN
967 DO i = 1, helium%atoms
968 DO c = 1, 3
969 r(c) = helium%pos(c, i, n) - helium%pos(c, i, n + 1)
970 END DO
971 CALL helium_pbc(helium, r)
972 t = t + r(1)*r(1) + r(2)*r(2) + r(3)*r(3)
973 DO j = 1, i - 1
974 DO c = 1, 3
975 r(c) = helium%pos(c, i, n) - helium%pos(c, j, n)
976 rp(c) = helium%pos(c, i, n + 1) - helium%pos(c, j, n + 1)
977 END DO
978 s = s + helium_eval_expansion(helium, r, rp, work3)
979 END DO
980 END DO
981 ELSE
982 DO i = 1, helium%atoms
983 DO c = 1, 3
984 r(c) = helium%pos(c, i, n) - helium%pos(c, helium%permutation(i), 1)
985 END DO
986 CALL helium_pbc(helium, r)
987 t = t + r(1)*r(1) + r(2)*r(2) + r(3)*r(3)
988 DO j = 1, i - 1
989 DO c = 1, 3
990 r(c) = helium%pos(c, i, n) - helium%pos(c, j, n)
991 rp(c) = helium%pos(c, helium%permutation(i), 1) - helium%pos(c, helium%permutation(j), 1)
992 END DO
993 s = s + helium_eval_expansion(helium, r, rp, work3)
994 END DO
995 END DO
996 END IF
997 t = t/(2.0_dp*helium%tau*helium%hb2m)
998 res = s + t
999 DEALLOCATE (work3)
1000
1001 END FUNCTION helium_link_action
1002
1003! **************************************************************************************************
1004!> \brief ...
1005!> \param helium ...
1006!> \return ...
1007! **************************************************************************************************
1008 FUNCTION helium_total_action(helium) RESULT(res)
1009
1010 TYPE(helium_solvent_type), INTENT(INOUT) :: helium
1011 REAL(kind=dp) :: res
1012
1013 INTEGER :: i
1014 REAL(kind=dp) :: s
1015
1016 s = 0.0_dp
1017 DO i = 1, helium%beads
1018 s = s + helium_link_action(helium, i)
1019 END DO
1020 res = s
1021
1022 END FUNCTION helium_total_action
1023
1024! **************************************************************************************************
1025!> \brief ...
1026!> \param helium ...
1027!> \param part ...
1028!> \param ref_bead ...
1029!> \param delta_bead ...
1030!> \param d ...
1031! **************************************************************************************************
1032 SUBROUTINE helium_delta_pos(helium, part, ref_bead, delta_bead, d)
1033
1034 TYPE(helium_solvent_type), INTENT(INOUT) :: helium
1035 INTEGER, INTENT(IN) :: part, ref_bead, delta_bead
1036 REAL(kind=dp), DIMENSION(3), INTENT(OUT) :: d
1037
1038 INTEGER :: b, bead, db, nbead, np, p
1039 REAL(kind=dp), DIMENSION(3) :: r
1040
1041 b = helium%beads
1042
1043 d(:) = 0.0_dp
1044 IF (delta_bead > 0) THEN
1045 bead = ref_bead
1046 p = part
1047 db = delta_bead
1048 DO
1049 IF (db < 1) EXIT
1050 nbead = bead + 1
1051 np = p
1052 IF (nbead > b) THEN
1053 nbead = nbead - b
1054 np = helium%permutation(np)
1055 END IF
1056 r(:) = helium%pos(:, p, bead) - helium%pos(:, np, nbead)
1057 CALL helium_pbc(helium, r)
1058 d(:) = d(:) + r(:)
1059 bead = nbead
1060 p = np
1061 db = db - 1
1062 END DO
1063 ELSE IF (delta_bead < 0) THEN
1064 bead = ref_bead
1065 p = part
1066 db = delta_bead
1067 DO
1068 IF (db >= 0) EXIT
1069 nbead = bead - 1
1070 np = p
1071 IF (nbead < 1) THEN
1072 nbead = nbead + b
1073 np = helium%iperm(np)
1074 END IF
1075 r(:) = helium%pos(:, p, bead) - helium%pos(:, np, nbead)
1076 CALL helium_pbc(helium, r)
1077 d(:) = d(:) + r(:)
1078 bead = nbead
1079 p = np
1080 db = db + 1
1081 END DO
1082 END IF
1083 END SUBROUTINE helium_delta_pos
1084
1085#endif
1086
1087END MODULE helium_interactions
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
helium_common
Independent helium subroutines shared with other modules.
Definition helium_common.F:14
helium_common::helium_pbc
subroutine, public helium_pbc(helium, r, enforce)
General PBC routine for helium.
Definition helium_common.F:81
helium_common::helium_eval_chain
real(kind=dp) function, public helium_eval_chain(helium, rchain, nchain, aij, vcoef, energy)
Calculate the pair-product action or energy by evaluating the power series expansion according to Eq....
Definition helium_common.F:936
helium_common::helium_eval_expansion
real(kind=dp) function, public helium_eval_expansion(helium, r, rp, work, action)
Calculate the pair-product action or energy by evaluating the power series expansion according to Eq....
Definition helium_common.F:832
helium_common::helium_spline
real(kind=dp) function, public helium_spline(spl, xx)
...
Definition helium_common.F:1423
helium_interactions
Methods that handle helium-solvent and helium-helium interactions.
Definition helium_interactions.F:13
helium_interactions::helium_solute_e_f
subroutine, public helium_solute_e_f(pint_env, helium, energy)
Calculate total helium-solute interaction energy and forces.
Definition helium_interactions.F:450
helium_interactions::helium_total_pair_action
real(kind=dp) function, public helium_total_pair_action(helium)
Computes the total pair action of the helium.
Definition helium_interactions.F:217
helium_interactions::helium_total_inter_action
real(kind=dp) function, public helium_total_inter_action(pint_env, helium)
Computes the total interaction of the helium with the solute.
Definition helium_interactions.F:298
helium_interactions::helium_intpot_scan
subroutine, public helium_intpot_scan(pint_env, helium_env)
Scan the helium-solute interaction energy within the periodic cell.
Definition helium_interactions.F:519
helium_interactions::helium_bead_solute_e_f
subroutine, public helium_bead_solute_e_f(pint_env, helium, helium_part_index, helium_slice_index, helium_r_opt, energy, force)
Calculate general helium-solute interaction energy (and forces) between one helium bead and the corre...
Definition helium_interactions.F:354
helium_interactions::helium_vij
elemental real(kind=dp) function, public helium_vij(r)
Helium-helium pair interaction potential.
Definition helium_interactions.F:826
helium_interactions::helium_calc_energy
subroutine, public helium_calc_energy(helium, pint_env)
Calculate the helium energy (including helium-solute interaction)
Definition helium_interactions.F:76
helium_interactions::helium_total_link_action
real(kind=dp) function, public helium_total_link_action(helium)
Computes the total harmonic link action of the helium.
Definition helium_interactions.F:180
helium_nnp
Methods dealing with Neural Network interaction potential.
Definition helium_nnp.F:13
helium_nnp::helium_nnp_print
subroutine, public helium_nnp_print(nnp, print_section, ind_he)
Print properties according to the requests in input file.
Definition helium_nnp.F:168
helium_types
Data types representing superfluid helium.
Definition helium_types.F:15
helium_types::e_id_potential
integer, parameter, public e_id_potential
Definition helium_types.F:38
helium_types::e_id_thermo
integer, parameter, public e_id_thermo
Definition helium_types.F:38
helium_types::e_id_virial
integer, parameter, public e_id_virial
Definition helium_types.F:38
helium_types::e_id_interact
integer, parameter, public e_id_interact
Definition helium_types.F:38
helium_types::e_id_kinetic
integer, parameter, public e_id_kinetic
Definition helium_types.F:38
helium_types::e_id_total
integer, parameter, public e_id_total
Energy contributions - symbolic names for indexing energy arrays.
Definition helium_types.F:38
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:18
input_constants::helium_solute_intpot_mwater
integer, parameter, public helium_solute_intpot_mwater
Definition input_constants.F:337
input_constants::helium_solute_intpot_none
integer, parameter, public helium_solute_intpot_none
Definition input_constants.F:337
input_constants::helium_sampling_worm
integer, parameter, public helium_sampling_worm
Definition input_constants.F:362
input_constants::helium_solute_intpot_nnp
integer, parameter, public helium_solute_intpot_nnp
Definition input_constants.F:337
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition input_section_types.F:735
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
nnp_acsf
Functionality for atom centered symmetry functions for neural network potentials.
Definition nnp_acsf.F:16
nnp_acsf::nnp_calc_acsf
subroutine, public nnp_calc_acsf(nnp, i, calc_forces, arc, stress)
Calculate atom centered symmetry functions for given atom i.
Definition nnp_acsf.F:88
nnp_acsf::nnp_prepare_neighbor_cache
subroutine, public nnp_prepare_neighbor_cache(nnp)
Prepare or update the linked-cell / Verlet cache for the current geometry. Lazily allocates nnpcell_l...
Definition nnp_acsf.F:689
nnp_environment_types
Data types for neural network potentials.
Definition nnp_environment_types.F:14
nnp_force
Methods dealing with Neural Network potentials.
Definition nnp_force.F:15
nnp_force::nnp_scatter_dgdr_to_forces
subroutine, public nnp_scatter_dgdr_to_forces(nnp, ind, i, denergydsym, force_xyz)
Scatter the per-neighbour dG/dr arrays held in the nnp_neighbor_workspace into a per-atom Cartesian f...
Definition nnp_force.F:816
nnp_model
Methods dealing with core routines for artificial neural networks.
Definition nnp_model.F:13
nnp_model::nnp_predict
subroutine, public nnp_predict(arc, nnp, i_com)
Predict energy by evaluating neural network.
Definition nnp_model.F:82
nnp_model::nnp_gradients
subroutine, public nnp_gradients(arc, nnp, i_com, denergydsym)
Calculate gradients of neural network.
Definition nnp_model.F:169
physcon
Definition of physical constants:
Definition physcon.F:68
physcon::kelvin
real(kind=dp), parameter, public kelvin
Definition physcon.F:165
physcon::angstrom
real(kind=dp), parameter, public angstrom
Definition physcon.F:144
pint_types
Definition pint_types.F:8
splines_types
routines for handling splines_types
Definition splines_types.F:14
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
helium_types::helium_solvent_p_type
data structure for array of solvent helium environments
Definition helium_types.F:306
helium_types::helium_solvent_type
data structure for solvent helium
Definition helium_types.F:93
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
nnp_environment_types::nnp_type
Main data type collecting all relevant data for neural network potentials.
Definition nnp_environment_types.F:93
pint_types::pint_env_type
environment for a path integral run
Definition pint_types.F:112
splines_types::spline_data_type
Data-structure that holds all needed information about a specific spline interpolation.
Definition splines_types.F:40