Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2023 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 helium_common, ONLY: helium_eval_chain,&
16 : helium_eval_expansion,&
17 : helium_pbc,&
18 : helium_spline
19 : USE helium_types, ONLY: e_id_interact,&
20 : e_id_kinetic,&
21 : e_id_potential,&
22 : e_id_thermo,&
23 : e_id_total,&
24 : e_id_virial,&
25 : helium_solvent_p_type,&
26 : helium_solvent_type
27 : USE input_constants, ONLY: helium_sampling_worm,&
28 : helium_solute_intpot_mwater,&
29 : helium_solute_intpot_none
30 : USE kinds, ONLY: dp
31 : USE physcon, ONLY: angstrom,&
32 : kelvin
33 : USE pint_types, ONLY: pint_env_type
34 : USE splines_types, ONLY: spline_data_type
35 : #include "../base/base_uses.f90"
36 :
37 : IMPLICIT NONE
38 :
39 : PRIVATE
40 :
41 : LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
42 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'helium_interactions'
43 :
44 : PUBLIC :: helium_calc_energy
45 : PUBLIC :: helium_total_link_action
46 : PUBLIC :: helium_total_pair_action
47 : PUBLIC :: helium_total_inter_action
48 : PUBLIC :: helium_solute_e_f
49 : PUBLIC :: helium_bead_solute_e_f
50 : PUBLIC :: helium_intpot_scan
51 : PUBLIC :: helium_vij
52 :
53 : CONTAINS
54 :
55 : ! ***************************************************************************
56 : !> \brief Calculate the helium energy (including helium-solute interaction)
57 : !> \param helium helium environment
58 : !> \param pint_env path integral environment
59 : !> \par History
60 : !> 2009-06 moved I/O out from here [lwalewski]
61 : !> \author hforbert
62 : ! **************************************************************************************************
63 7042 : SUBROUTINE helium_calc_energy(helium, pint_env)
64 : TYPE(helium_solvent_type), INTENT(INOUT) :: helium
65 : TYPE(pint_env_type), INTENT(IN) :: pint_env
66 :
67 : INTEGER :: b, bead, i, j, n
68 7042 : INTEGER, DIMENSION(:), POINTER :: perm
69 : LOGICAL :: nperiodic
70 : REAL(KIND=dp) :: a, cell_size, en, interac, kin, pot, &
71 : rmax, rmin, vkin
72 7042 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: work2, work3
73 7042 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: work
74 : REAL(KIND=dp), DIMENSION(3) :: r
75 7042 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: pos
76 : TYPE(spline_data_type), POINTER :: e0
77 :
78 7042 : pos => helium%pos
79 7042 : perm => helium%permutation
80 7042 : e0 => helium%e0
81 7042 : cell_size = 0.5_dp*helium%cell_size
82 7042 : nperiodic = .NOT. helium%periodic
83 7042 : n = helium%atoms
84 7042 : b = helium%beads
85 7042 : en = 0.0_dp
86 7042 : pot = 0.0_dp
87 7042 : rmin = 1.0e20_dp
88 7042 : rmax = 0.0_dp
89 : ALLOCATE (work(3, helium%beads + 1), &
90 : work2(helium%beads + 1), &
91 49294 : work3(SIZE(helium%uoffdiag, 1) + 1))
92 214004 : DO i = 1, n - 1
93 3502716 : DO j = i + 1, n
94 71130344 : DO bead = 1, b
95 274655240 : work(:, bead) = pos(:, i, bead) - pos(:, j, bead)
96 : END DO
97 13154848 : work(:, b + 1) = pos(:, perm(i), 1) - pos(:, perm(j), 1)
98 3288712 : en = en + helium_eval_chain(helium, work, b + 1, work2, work3, energy=.TRUE.)
99 71337306 : DO bead = 1, b
100 67841632 : a = work2(bead)
101 67841632 : IF (a < rmin) rmin = a
102 67841632 : IF (a > rmax) rmax = a
103 71130344 : IF ((a < cell_size) .OR. nperiodic) THEN
104 63517637 : pot = pot + helium_spline(helium%vij, a)
105 : END IF
106 : END DO
107 : END DO
108 : END DO
109 7042 : DEALLOCATE (work, work2, work3)
110 7042 : pot = pot/b
111 7042 : en = en/b
112 :
113 : ! helium-solute interaction energy (all beads of all particles)
114 7042 : interac = 0.0_dp
115 7042 : IF (helium%solute_present) THEN
116 3632 : CALL helium_solute_e(pint_env, helium, interac)
117 : END IF
118 7042 : interac = interac/b
119 :
120 : !TODO:
121 7042 : vkin = 0.0_dp
122 : ! vkin = helium_virial_energy(helium)
123 :
124 7042 : kin = 0.0_dp
125 221046 : DO i = 1, n
126 856016 : r(:) = pos(:, i, b) - pos(:, perm(i), 1)
127 214004 : CALL helium_pbc(helium, r)
128 214004 : kin = kin + r(1)*r(1) + r(2)*r(2) + r(3)*r(3)
129 4413186 : DO bead = 2, b
130 16768560 : r(:) = pos(:, i, bead - 1) - pos(:, i, bead)
131 4192140 : CALL helium_pbc(helium, r)
132 4406144 : kin = kin + r(1)*r(1) + r(2)*r(2) + r(3)*r(3)
133 : END DO
134 : END DO
135 7042 : kin = 1.5_dp*n/helium%tau - 0.5*kin/(b*helium%tau**2*helium%hb2m)
136 :
137 : ! TODO: move printing somewhere else ?
138 : ! print *,"POT = ",(pot/n+helium%e_corr)*kelvin,"K"
139 : ! print *,"INTERAC = ",interac*kelvin,"K"
140 : ! print *,"RMIN= ",rmin*angstrom,"A"
141 : ! print *,"RMAX= ",rmax*angstrom,"A"
142 : ! print *,"EVIRIAL not valid!"
143 : ! print *,"ETHERMO= ",((en+kin)/n+helium%e_corr)*kelvin,"K"
144 : ! print *,"ECORR= ",helium%e_corr*kelvin,"K"
145 : !! kin = helium_total_action(helium)
146 : !! print *,"ACTION= ",kin
147 : ! print *,"WINDING#= ",helium_calc_winding(helium)
148 :
149 7042 : helium%energy_inst(e_id_potential) = pot/n + helium%e_corr
150 7042 : helium%energy_inst(e_id_kinetic) = (en - pot + kin)/n
151 7042 : helium%energy_inst(e_id_interact) = interac
152 7042 : helium%energy_inst(e_id_thermo) = (en + kin)/n + helium%e_corr
153 7042 : helium%energy_inst(e_id_virial) = vkin ! 0.0_dp at the moment
154 7042 : helium%energy_inst(e_id_total) = helium%energy_inst(e_id_thermo)
155 : ! Once vkin is properly implemented, switch to:
156 : ! helium%energy_inst(e_id_total) = (en+vkin)/n+helium%e_corr
157 :
158 14084 : END SUBROUTINE helium_calc_energy
159 :
160 : ! ***************************************************************************
161 : !> \brief Computes the total harmonic link action of the helium
162 : !> \param helium ...
163 : !> \return ...
164 : !> \date 2016-05-03
165 : !> \author Felix Uhl
166 : ! **************************************************************************************************
167 120 : REAL(KIND=dp) FUNCTION helium_total_link_action(helium) RESULT(linkaction)
168 :
169 : TYPE(helium_solvent_type), INTENT(IN) :: helium
170 :
171 : INTEGER :: iatom, ibead
172 120 : INTEGER, DIMENSION(:), POINTER :: perm
173 : REAL(KIND=dp), DIMENSION(3) :: r
174 :
175 120 : perm => helium%permutation
176 120 : linkaction = 0.0_dp
177 :
178 : ! Harmonic Link action
179 : ! (r(m-1) - r(m))**2/(4*lambda*tau)
180 2172 : DO ibead = 1, helium%beads - 1
181 47586 : DO iatom = 1, helium%atoms
182 181656 : r(:) = helium%pos(:, iatom, ibead) - helium%pos(:, iatom, ibead + 1)
183 45414 : CALL helium_pbc(helium, r)
184 47466 : linkaction = linkaction + (r(1)*r(1) + r(2)*r(2) + r(3)*r(3))
185 : END DO
186 : END DO
187 2610 : DO iatom = 1, helium%atoms
188 : ! choose last bead connection according to permutation table
189 9960 : r(:) = helium%pos(:, iatom, helium%beads) - helium%pos(:, perm(iatom), 1)
190 2490 : CALL helium_pbc(helium, r)
191 2610 : linkaction = linkaction + (r(1)*r(1) + r(2)*r(2) + r(3)*r(3))
192 : END DO
193 120 : linkaction = linkaction/(2.0_dp*helium%tau*helium%hb2m)
194 :
195 120 : END FUNCTION helium_total_link_action
196 :
197 : ! ***************************************************************************
198 : !> \brief Computes the total pair action of the helium
199 : !> \param helium ...
200 : !> \return ...
201 : !> \date 2016-05-03
202 : !> \author Felix Uhl
203 : ! **************************************************************************************************
204 120 : REAL(KIND=dp) FUNCTION helium_total_pair_action(helium) RESULT(pairaction)
205 :
206 : TYPE(helium_solvent_type), INTENT(INOUT) :: helium
207 :
208 : INTEGER :: iatom, ibead, jatom, opatom, patom
209 120 : INTEGER, DIMENSION(:), POINTER :: perm
210 120 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: work3
211 : REAL(KIND=dp), DIMENSION(3) :: r, rp
212 :
213 360 : ALLOCATE (work3(SIZE(helium%uoffdiag, 1) + 1))
214 120 : perm => helium%permutation
215 120 : pairaction = 0.0_dp
216 :
217 : ! He-He pair action
218 2172 : DO ibead = 1, helium%beads - 1
219 45534 : DO iatom = 1, helium%atoms - 1
220 698706 : DO jatom = iatom + 1, helium%atoms
221 2613168 : r(:) = helium%pos(:, iatom, ibead) - helium%pos(:, jatom, ibead)
222 2613168 : rp(:) = helium%pos(:, iatom, ibead + 1) - helium%pos(:, jatom, ibead + 1)
223 696654 : pairaction = pairaction + helium_eval_expansion(helium, r, rp, work3)
224 : END DO
225 : END DO
226 : END DO
227 : !Ensure right permutation for pair action of last and first beads.
228 2490 : DO iatom = 1, helium%atoms - 1
229 37710 : DO jatom = iatom + 1, helium%atoms
230 140880 : r(:) = helium%pos(:, iatom, helium%beads) - helium%pos(:, jatom, helium%beads)
231 140880 : rp(:) = helium%pos(:, perm(iatom), 1) - helium%pos(:, perm(jatom), 1)
232 37590 : pairaction = pairaction + helium_eval_expansion(helium, r, rp, work3)
233 : END DO
234 : END DO
235 :
236 : ! correct for open worm configurations
237 120 : IF (.NOT. helium%worm_is_closed) THEN
238 : ! special treatment if double bead is first bead
239 0 : iatom = helium%worm_atom_idx
240 0 : IF (helium%worm_bead_idx == 1) THEN
241 : ! patom is the atom in front of the lone head bead
242 0 : patom = helium%iperm(iatom)
243 : ! go through all atoms
244 0 : DO jatom = 1, helium%atoms
245 0 : IF (jatom == helium%worm_atom_idx) CYCLE
246 0 : opatom = helium%iperm(jatom)
247 : ! subtract pair action for closed link
248 0 : r(:) = helium%pos(:, iatom, 1) - helium%pos(:, jatom, 1)
249 0 : rp(:) = helium%pos(:, patom, helium%beads) - helium%pos(:, opatom, helium%beads)
250 0 : pairaction = pairaction - helium_eval_expansion(helium, r, rp, work3)
251 : ! and add corrected extra link
252 : ! rp stays the same
253 0 : r(:) = helium%worm_xtra_bead(:) - helium%pos(:, jatom, 1)
254 0 : pairaction = pairaction + helium_eval_expansion(helium, r, rp, work3)
255 : END DO
256 : ELSE
257 : ! bead stays constant
258 0 : ibead = helium%worm_bead_idx
259 : ! go through all atoms
260 0 : DO jatom = 1, helium%atoms
261 0 : IF (jatom == helium%worm_atom_idx) CYCLE
262 : ! subtract pair action for closed link
263 0 : r(:) = helium%pos(:, iatom, ibead) - helium%pos(:, jatom, ibead)
264 0 : rp(:) = helium%pos(:, iatom, ibead - 1) - helium%pos(:, jatom, ibead - 1)
265 0 : pairaction = pairaction - helium_eval_expansion(helium, r, rp, work3)
266 : ! and add corrected extra link
267 : ! rp stays the same
268 0 : r(:) = helium%worm_xtra_bead(:) - helium%pos(:, jatom, ibead)
269 0 : pairaction = pairaction + helium_eval_expansion(helium, r, rp, work3)
270 : END DO
271 : END IF
272 : END IF
273 120 : DEALLOCATE (work3)
274 :
275 120 : END FUNCTION helium_total_pair_action
276 :
277 : ! ***************************************************************************
278 : !> \brief Computes the total interaction of the helium with the solute
279 : !> \param pint_env ...
280 : !> \param helium ...
281 : !> \return ...
282 : !> \date 2016-05-03
283 : !> \author Felix Uhl
284 : ! **************************************************************************************************
285 120 : REAL(KIND=dp) FUNCTION helium_total_inter_action(pint_env, helium) RESULT(interaction)
286 :
287 : TYPE(pint_env_type), INTENT(IN) :: pint_env
288 : TYPE(helium_solvent_type), INTENT(IN) :: helium
289 :
290 : INTEGER :: iatom, ibead
291 : REAL(KIND=dp) :: e
292 :
293 120 : interaction = 0.0_dp
294 :
295 : ! InterAction with solute
296 120 : IF (helium%solute_present) THEN
297 1564 : DO ibead = 1, helium%beads
298 27068 : DO iatom = 1, helium%atoms
299 :
300 : CALL helium_bead_solute_e_f(pint_env, helium, &
301 25504 : iatom, ibead, helium%pos(:, iatom, ibead), e)
302 26976 : interaction = interaction + e
303 : END DO
304 : END DO
305 92 : IF (helium%sampling_method == helium_sampling_worm) THEN
306 12 : IF (.NOT. helium%worm_is_closed) THEN
307 : ! subtract half of tail bead interaction again
308 : CALL helium_bead_solute_e_f(pint_env, helium, &
309 : helium%worm_atom_idx, helium%worm_bead_idx, &
310 0 : helium%pos(:, helium%worm_atom_idx, helium%worm_bead_idx), e)
311 0 : interaction = interaction - 0.5_dp*e
312 : ! add half of head bead interaction
313 : CALL helium_bead_solute_e_f(pint_env, helium, &
314 : helium%worm_atom_idx, helium%worm_bead_idx, &
315 0 : helium%worm_xtra_bead, e)
316 0 : interaction = interaction + 0.5_dp*e
317 : END IF
318 : END IF
319 : END IF
320 :
321 120 : interaction = interaction*helium%tau
322 :
323 120 : END FUNCTION helium_total_inter_action
324 :
325 : ! ***************************************************************************
326 : !> \brief Calculate general helium-solute interaction energy (and forces)
327 : !> between one helium bead and the corresponding solute time slice.
328 : !> \param pint_env path integral environment
329 : !> \param helium ...
330 : !> \param helium_part_index helium particle index
331 : !> \param helium_slice_index helium time slice index
332 : !> \param helium_r_opt explicit helium bead coordinates (optional)
333 : !> \param energy calculated energy
334 : !> \param force calculated force (if requested)
335 : !> \par History
336 : !> 2019-09 Added multiple-time striding in imag. time [cschran]
337 : !> \author Lukasz Walewski
338 : ! **************************************************************************************************
339 5509296 : SUBROUTINE helium_bead_solute_e_f(pint_env, helium, helium_part_index, &
340 : helium_slice_index, helium_r_opt, energy, force)
341 :
342 : TYPE(pint_env_type), INTENT(IN) :: pint_env
343 : TYPE(helium_solvent_type), INTENT(IN) :: helium
344 : INTEGER, INTENT(IN) :: helium_part_index, helium_slice_index
345 : REAL(KIND=dp), DIMENSION(3), INTENT(IN), OPTIONAL :: helium_r_opt
346 : REAL(KIND=dp), INTENT(OUT) :: energy
347 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT), &
348 : OPTIONAL, POINTER :: force
349 :
350 : INTEGER :: hbeads, hi, qi, stride
351 : REAL(KIND=dp), DIMENSION(3) :: helium_r
352 5509296 : REAL(KIND=dp), DIMENSION(:), POINTER :: my_force
353 :
354 5509296 : hbeads = helium%beads
355 : ! helium bead index that is invariant wrt the rotations
356 5509296 : hi = MOD(helium_slice_index - 1 + hbeads + helium%relrot, hbeads) + 1
357 : ! solute bead index that belongs to hi helium index
358 5509296 : qi = ((hi - 1)*pint_env%p)/hbeads + 1
359 :
360 : ! coordinates of the helium bead
361 5509296 : IF (PRESENT(helium_r_opt)) THEN
362 1429598 : helium_r(:) = helium_r_opt(:)
363 : ELSE
364 16318792 : helium_r(:) = helium%pos(:, helium_part_index, helium_slice_index)
365 : END IF
366 :
367 11018592 : SELECT CASE (helium%solute_interaction)
368 :
369 : CASE (helium_solute_intpot_mwater)
370 5509296 : IF (PRESENT(force)) THEN
371 54652288 : force(:, :) = 0.0_dp
372 1171264 : my_force => force(qi, :)
373 : CALL helium_intpot_model_water( &
374 : pint_env%x(qi, :), &
375 : helium, &
376 : helium_r, &
377 : energy, &
378 : my_force &
379 1171264 : )
380 : ELSE
381 : CALL helium_intpot_model_water( &
382 : pint_env%x(qi, :), &
383 : helium, &
384 : helium_r, &
385 : energy &
386 4338032 : )
387 : END IF
388 :
389 : CASE (helium_solute_intpot_none)
390 0 : energy = 0.0_dp
391 5509296 : IF (PRESENT(force)) THEN
392 0 : force(:, :) = 0.0_dp
393 : END IF
394 :
395 : CASE DEFAULT
396 :
397 : END SELECT
398 :
399 : ! Account for Imaginary time striding in forces:
400 5509296 : IF (PRESENT(force)) THEN
401 1171264 : IF (hbeads < pint_env%p) THEN
402 3072 : stride = pint_env%p/hbeads
403 915456 : force = force*REAL(stride, dp)
404 : END IF
405 : END IF
406 :
407 5509296 : END SUBROUTINE helium_bead_solute_e_f
408 :
409 : ! ***************************************************************************
410 : !> \brief Calculate total helium-solute interaction energy and forces.
411 : !> \param pint_env path integral environment
412 : !> \param helium ...
413 : !> \param energy calculated interaction energy
414 : !> \author Lukasz Walewski
415 : ! **************************************************************************************************
416 2642 : SUBROUTINE helium_solute_e_f(pint_env, helium, energy)
417 :
418 : TYPE(pint_env_type), INTENT(IN) :: pint_env
419 : TYPE(helium_solvent_type), INTENT(INOUT) :: helium
420 : REAL(KIND=dp), INTENT(OUT) :: energy
421 :
422 : INTEGER :: ia, ib, jb, jc
423 : REAL(KIND=dp) :: my_energy
424 2642 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: force
425 :
426 2642 : NULLIFY (force)
427 2642 : force => helium%force_inst
428 :
429 2642 : energy = 0.0_dp
430 123044 : force(:, :) = 0.0_dp
431 :
432 : ! calculate the total interaction energy and gradients between the
433 : ! solute and the helium, sum over all beads of all He particles
434 75846 : DO ia = 1, helium%atoms
435 1247110 : DO ib = 1, helium%beads
436 : CALL helium_bead_solute_e_f(pint_env, helium, ia, ib, &
437 1171264 : energy=my_energy, force=helium%rtmp_p_ndim_2d)
438 1171264 : energy = energy + my_energy
439 6015540 : DO jb = 1, pint_env%p
440 48881984 : DO jc = 1, pint_env%ndim
441 47710720 : force(jb, jc) = force(jb, jc) + helium%rtmp_p_ndim_2d(jb, jc)
442 : END DO
443 : END DO
444 : END DO
445 : END DO
446 :
447 2642 : END SUBROUTINE helium_solute_e_f
448 :
449 : ! ***************************************************************************
450 : !> \brief Calculate total helium-solute interaction energy.
451 : !> \param pint_env path integral environment
452 : !> \param helium ...
453 : !> \param energy calculated interaction energy
454 : !> \author Lukasz Walewski
455 : ! **************************************************************************************************
456 3632 : SUBROUTINE helium_solute_e(pint_env, helium, energy)
457 :
458 : TYPE(pint_env_type), INTENT(IN) :: pint_env
459 : TYPE(helium_solvent_type), INTENT(IN) :: helium
460 : REAL(KIND=dp), INTENT(OUT) :: energy
461 :
462 : INTEGER :: ia, ib
463 : REAL(KIND=dp) :: my_energy
464 :
465 3632 : energy = 0.0_dp
466 :
467 108516 : DO ia = 1, helium%atoms
468 1786660 : DO ib = 1, helium%beads
469 1678144 : CALL helium_bead_solute_e_f(pint_env, helium, ia, ib, energy=my_energy)
470 1783028 : energy = energy + my_energy
471 : END DO
472 : END DO
473 :
474 3632 : END SUBROUTINE helium_solute_e
475 :
476 : ! ***************************************************************************
477 : !> \brief Scan the helium-solute interaction energy within the periodic cell
478 : !> \param pint_env ...
479 : !> \param helium_env ...
480 : !> \date 2014-01-22
481 : !> \par History
482 : !> 2016-07-14 Modified to work with independent helium_env [cschran]
483 : !> \author Lukasz Walewski
484 : ! **************************************************************************************************
485 0 : SUBROUTINE helium_intpot_scan(pint_env, helium_env)
486 :
487 : TYPE(pint_env_type), INTENT(IN) :: pint_env
488 : TYPE(helium_solvent_p_type), DIMENSION(:), POINTER :: helium_env
489 :
490 : CHARACTER(len=*), PARAMETER :: routineN = 'helium_intpot_scan'
491 :
492 : INTEGER :: handle, ic, ix, iy, iz, k, nbin
493 : LOGICAL :: wrapped
494 : REAL(KIND=dp) :: delr, my_en, ox, oy, oz
495 : REAL(kind=dp), DIMENSION(3) :: pbc1, pbc2, pos
496 :
497 0 : CALL timeset(routineN, handle)
498 :
499 : ! Perform scan only on ionode, since this is only used to output the intpot
500 0 : IF (pint_env%logger%para_env%is_source()) THEN
501 : ! Assume ionode always to have at least one helium_env
502 0 : k = 1
503 0 : helium_env(k)%helium%rho_inst(1, :, :, :) = 0.0_dp
504 0 : nbin = helium_env(k)%helium%rho_nbin
505 0 : delr = helium_env(k)%helium%rho_delr
506 0 : helium_env(k)%helium%center(:) = 0.0_dp
507 0 : ox = helium_env(k)%helium%center(1) - helium_env(k)%helium%rho_maxr/2.0_dp
508 0 : oy = helium_env(k)%helium%center(2) - helium_env(k)%helium%rho_maxr/2.0_dp
509 0 : oz = helium_env(k)%helium%center(3) - helium_env(k)%helium%rho_maxr/2.0_dp
510 :
511 0 : DO ix = 1, nbin
512 0 : DO iy = 1, nbin
513 0 : DO iz = 1, nbin
514 :
515 : ! put the probe in the center of the current voxel
516 0 : pos(:) = (/ox + (ix - 0.5_dp)*delr, oy + (iy - 0.5_dp)*delr, oz + (iz - 0.5_dp)*delr/)
517 :
518 : ! calc interaction energy for the current probe position
519 0 : helium_env(k)%helium%pos(:, 1, 1) = pos(:)
520 0 : CALL helium_bead_solute_e_f(pint_env, helium_env(k)%helium, 1, 1, energy=my_en)
521 :
522 : ! check if the probe fits within the unit cell
523 0 : pbc1(:) = pos(:) - helium_env(k)%helium%center
524 0 : pbc2(:) = pbc1(:)
525 0 : CALL helium_pbc(helium_env(k)%helium, pbc2)
526 0 : wrapped = .FALSE.
527 0 : DO ic = 1, 3
528 0 : IF (ABS(pbc1(ic) - pbc2(ic)) .GT. 10.0_dp*EPSILON(0.0_dp)) THEN
529 0 : wrapped = .TRUE.
530 : END IF
531 : END DO
532 :
533 : ! set the interaction energy value
534 0 : IF (wrapped) THEN
535 0 : helium_env(k)%helium%rho_inst(1, ix, iy, iz) = 0.0_dp
536 : ELSE
537 0 : helium_env(k)%helium%rho_inst(1, ix, iy, iz) = my_en
538 : END IF
539 :
540 : END DO
541 : END DO
542 : END DO
543 : END IF
544 :
545 0 : CALL timestop(handle)
546 0 : END SUBROUTINE helium_intpot_scan
547 :
548 : ! ***************************************************************************
549 : !> \brief Calculate model helium-solute interaction energy and forces
550 : !> between one helium bead and the corresponding solute time
551 : !> slice asuming water solute.
552 : !> \param solute_x solute positions ARR(3*NATOMS)
553 : !> to global atom indices
554 : !> \param helium only needed for helium_pbc call at the moment
555 : !> \param helium_x helium bead position ARR(3)
556 : !> \param energy calculated interaction energy
557 : !> \param force ...
558 : !> \author Felix Uhl
559 : ! **************************************************************************************************
560 5509296 : SUBROUTINE helium_intpot_model_water(solute_x, helium, helium_x, energy, force)
561 :
562 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: solute_x
563 : TYPE(helium_solvent_type), INTENT(IN) :: helium
564 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: helium_x
565 : REAL(KIND=dp), INTENT(OUT) :: energy
566 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT), &
567 : OPTIONAL, POINTER :: force
568 :
569 : INTEGER :: i, ig
570 : REAL(KIND=dp) :: d, d2, dd, ep, eps, s1, s2, sig
571 : REAL(KIND=dp), DIMENSION(3) :: dr, solute_r
572 :
573 5509296 : energy = 0.0_dp
574 5509296 : IF (PRESENT(force)) THEN
575 11712640 : force(:) = 0.0_dp
576 : END IF
577 :
578 5509296 : sig = 2.69_dp ! 1.4 Angstrom
579 5509296 : eps = 60.61e-6_dp ! 19 K
580 5509296 : s1 = 0.0_dp
581 22037184 : DO i = 1, SIZE(helium%solute_element)
582 22037184 : IF (helium%solute_element(i) == "H ") THEN
583 11018592 : ig = i - 1
584 11018592 : solute_r(1) = solute_x(3*ig + 1)
585 11018592 : solute_r(2) = solute_x(3*ig + 2)
586 11018592 : solute_r(3) = solute_x(3*ig + 3)
587 44074368 : dr(:) = solute_r(:) - helium_x(:)
588 11018592 : CALL helium_pbc(helium, dr)
589 11018592 : d2 = dr(1)*dr(1) + dr(2)*dr(2) + dr(3)*dr(3)
590 11018592 : d = SQRT(d2)
591 11018592 : dd = (sig/d)**6
592 11018592 : ep = 4.0_dp*eps*dd*(dd - 1.0_dp)
593 11018592 : s1 = s1 + ep
594 11018592 : s2 = 24.0_dp*eps*dd*(2.0_dp*dd - 1.0_dp)/d2
595 11018592 : IF (PRESENT(force)) THEN
596 2342528 : force(3*ig + 1) = force(3*ig + 1) + s2*dr(1)
597 2342528 : force(3*ig + 2) = force(3*ig + 2) + s2*dr(2)
598 2342528 : force(3*ig + 3) = force(3*ig + 3) + s2*dr(3)
599 : END IF
600 : END IF
601 : END DO ! i = 1, num_hydrogen
602 5509296 : energy = energy + s1
603 :
604 5509296 : sig = 5.01_dp ! 2.6 Angstrom
605 5509296 : eps = 104.5e-6_dp ! 33 K
606 5509296 : s1 = 0.0_dp
607 22037184 : DO i = 1, SIZE(helium%solute_element)
608 22037184 : IF (helium%solute_element(i) == "O ") THEN
609 5509296 : ig = i - 1
610 5509296 : solute_r(1) = solute_x(3*ig + 1)
611 5509296 : solute_r(2) = solute_x(3*ig + 2)
612 5509296 : solute_r(3) = solute_x(3*ig + 3)
613 22037184 : dr(:) = solute_r(:) - helium_x(:)
614 5509296 : CALL helium_pbc(helium, dr)
615 5509296 : d2 = dr(1)*dr(1) + dr(2)*dr(2) + dr(3)*dr(3)
616 5509296 : d = SQRT(d2)
617 5509296 : dd = (sig/d)**6
618 5509296 : ep = 4.0_dp*eps*dd*(dd - 1.0_dp)
619 5509296 : s1 = s1 + ep
620 5509296 : s2 = 24.0_dp*eps*dd*(2.0_dp*dd - 1.0_dp)/d2
621 5509296 : IF (PRESENT(force)) THEN
622 1171264 : force(3*ig + 1) = force(3*ig + 1) + s2*dr(1)
623 1171264 : force(3*ig + 2) = force(3*ig + 2) + s2*dr(2)
624 1171264 : force(3*ig + 3) = force(3*ig + 3) + s2*dr(3)
625 : END IF
626 : END IF
627 : END DO ! i = 1, num_chlorine
628 5509296 : energy = energy + s1
629 :
630 5509296 : END SUBROUTINE helium_intpot_model_water
631 :
632 : ! ***************************************************************************
633 : !> \brief Helium-helium pair interaction potential.
634 : !> \param r ...
635 : !> \return ...
636 : ! **************************************************************************************************
637 1200 : ELEMENTAL FUNCTION helium_vij(r) RESULT(vij)
638 :
639 : REAL(kind=dp), INTENT(IN) :: r
640 : REAL(kind=dp) :: vij
641 :
642 : REAL(kind=dp) :: f, x, x2
643 :
644 1200 : x = angstrom*r/2.9673_dp
645 1200 : IF (x < 1.241314_dp) THEN
646 324 : x2 = 1.241314_dp/x - 1.0_dp
647 324 : f = EXP(-x2*x2)
648 : ELSE
649 : f = 1.0_dp
650 : END IF
651 1200 : x2 = 1.0_dp/(x*x)
652 : vij = 10.8_dp/kelvin*(544850.4_dp*EXP(-13.353384_dp*x) - f* &
653 1200 : ((0.1781_dp*x2 + 0.4253785_dp)*x2 + 1.3732412_dp)*x2*x2*x2)
654 1200 : END FUNCTION helium_vij
655 :
656 : #if 0
657 :
658 : ! this block is currently turned off
659 :
660 : ! ***************************************************************************
661 : !> \brief Helium-helium pair interaction potential's derivative.
662 : !> \param r ...
663 : !> \return ...
664 : ! **************************************************************************************************
665 : ELEMENTAL FUNCTION helium_d_vij(r) RESULT(dvij)
666 :
667 : REAL(kind=dp), INTENT(IN) :: r
668 : REAL(kind=dp) :: dvij
669 :
670 : REAL(kind=dp) :: f, fp, x, x2, y
671 :
672 : x = angstrom*r/2.9673_dp
673 : x = r/2.9673_dp
674 : x2 = 1.0_dp/(x*x)
675 : IF (x < 1.241314_dp) THEN
676 : y = 1.241314_dp/x - 1.0_dp
677 : f = EXP(-y*y)
678 : fp = 2.0_dp*1.241314_dp*f*y* &
679 : ((0.1781_dp*x2 + 0.4253785_dp)*x2 + 1.3732412_dp)*x2*x2*x2*x2
680 : ELSE
681 : f = 1.0_dp
682 : fp = 0.0_dp
683 : END IF
684 :
685 : dvij = angstrom*(10.8_dp/2.9673_dp)*( &
686 : (-13.353384_dp*544850.4_dp)*EXP(-13.353384_dp*x) - fp + &
687 : f*(((10.0_dp*0.1781_dp)*x2 + (8.0_dp*0.4253785_dp))*x2 + (6.0_dp*1.3732412_dp))* &
688 : x2*x2*x2/x)/(r*kelvin)
689 : END FUNCTION helium_d_vij
690 :
691 : ! **************************************************************************************************
692 : !> \brief ...
693 : !> \param helium ...
694 : !> \param n ...
695 : !> \param i ...
696 : !> \return ...
697 : ! **************************************************************************************************
698 : FUNCTION helium_atom_action(helium, n, i) RESULT(res)
699 :
700 : TYPE(helium_solvent_type), INTENT(INOUT) :: helium
701 : INTEGER, INTENT(IN) :: n, i
702 : REAL(KIND=dp) :: res
703 :
704 : INTEGER :: c, j
705 : REAL(KIND=dp) :: r(3), rp(3), s, t
706 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: work3
707 :
708 : ALLOCATE (work3(SIZE(helium%uoffdiag, 1) + 1))
709 : s = 0.0_dp
710 : t = 0.0_dp
711 : IF (n < helium%beads) THEN
712 : DO c = 1, 3
713 : r(c) = helium%pos(c, i, n) - helium%pos(c, i, n + 1)
714 : END DO
715 : CALL helium_pbc(helium, r)
716 : t = r(1)*r(1) + r(2)*r(2) + r(3)*r(3)
717 : DO j = 1, i - 1
718 : DO c = 1, 3
719 : r(c) = helium%pos(c, i, n) - helium%pos(c, j, n)
720 : rp(c) = helium%pos(c, i, n + 1) - helium%pos(c, j, n + 1)
721 : END DO
722 : s = s + helium_eval_expansion(helium, r, rp, work3)
723 : END DO
724 : DO j = i + 1, helium%atoms
725 : DO c = 1, 3
726 : r(c) = helium%pos(c, i, n) - helium%pos(c, j, n)
727 : rp(c) = helium%pos(c, i, n + 1) - helium%pos(c, j, n + 1)
728 : END DO
729 : s = s + helium_eval_expansion(helium, r, rp, work3)
730 : END DO
731 : ELSE
732 : DO c = 1, 3
733 : r(c) = helium%pos(c, i, n) - helium%pos(c, helium%permutation(i), 1)
734 : END DO
735 : CALL helium_pbc(helium, r)
736 : t = r(1)*r(1) + r(2)*r(2) + r(3)*r(3)
737 : DO j = 1, i - 1
738 : DO c = 1, 3
739 : r(c) = helium%pos(c, i, n) - helium%pos(c, j, n)
740 : rp(c) = helium%pos(c, helium%permutation(i), 1) - helium%pos(c, helium%permutation(j), 1)
741 : END DO
742 : s = s + helium_eval_expansion(helium, r, rp, work3)
743 : END DO
744 : DO j = i + 1, helium%atoms
745 : DO c = 1, 3
746 : r(c) = helium%pos(c, i, n) - helium%pos(c, j, n)
747 : rp(c) = helium%pos(c, helium%permutation(i), 1) - helium%pos(c, helium%permutation(j), 1)
748 : END DO
749 : s = s + helium_eval_expansion(helium, r, rp, work3)
750 : END DO
751 : END IF
752 : t = t/(2.0_dp*helium%tau*helium%hb2m)
753 : s = s*0.5_dp
754 : res = s + t
755 : DEALLOCATE (work3)
756 :
757 : END FUNCTION helium_atom_action
758 :
759 : ! **************************************************************************************************
760 : !> \brief ...
761 : !> \param helium ...
762 : !> \param n ...
763 : !> \return ...
764 : ! **************************************************************************************************
765 : FUNCTION helium_link_action(helium, n) RESULT(res)
766 :
767 : TYPE(helium_solvent_type), INTENT(INOUT) :: helium
768 : INTEGER, INTENT(IN) :: n
769 : REAL(KIND=dp) :: res
770 :
771 : INTEGER :: c, i, j
772 : REAL(KIND=dp) :: r(3), rp(3), s, t
773 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: work3
774 :
775 : ALLOCATE (work3(SIZE(helium%uoffdiag, 1) + 1))
776 : s = 0.0_dp
777 : t = 0.0_dp
778 : IF (n < helium%beads) THEN
779 : DO i = 1, helium%atoms
780 : DO c = 1, 3
781 : r(c) = helium%pos(c, i, n) - helium%pos(c, i, n + 1)
782 : END DO
783 : CALL helium_pbc(helium, r)
784 : t = t + r(1)*r(1) + r(2)*r(2) + r(3)*r(3)
785 : DO j = 1, i - 1
786 : DO c = 1, 3
787 : r(c) = helium%pos(c, i, n) - helium%pos(c, j, n)
788 : rp(c) = helium%pos(c, i, n + 1) - helium%pos(c, j, n + 1)
789 : END DO
790 : s = s + helium_eval_expansion(helium, r, rp, work3)
791 : END DO
792 : END DO
793 : ELSE
794 : DO i = 1, helium%atoms
795 : DO c = 1, 3
796 : r(c) = helium%pos(c, i, n) - helium%pos(c, helium%permutation(i), 1)
797 : END DO
798 : CALL helium_pbc(helium, r)
799 : t = t + r(1)*r(1) + r(2)*r(2) + r(3)*r(3)
800 : DO j = 1, i - 1
801 : DO c = 1, 3
802 : r(c) = helium%pos(c, i, n) - helium%pos(c, j, n)
803 : rp(c) = helium%pos(c, helium%permutation(i), 1) - helium%pos(c, helium%permutation(j), 1)
804 : END DO
805 : s = s + helium_eval_expansion(helium, r, rp, work3)
806 : END DO
807 : END DO
808 : END IF
809 : t = t/(2.0_dp*helium%tau*helium%hb2m)
810 : res = s + t
811 : DEALLOCATE (work3)
812 :
813 : END FUNCTION helium_link_action
814 :
815 : ! **************************************************************************************************
816 : !> \brief ...
817 : !> \param helium ...
818 : !> \return ...
819 : ! **************************************************************************************************
820 : FUNCTION helium_total_action(helium) RESULT(res)
821 :
822 : TYPE(helium_solvent_type), INTENT(INOUT) :: helium
823 : REAL(KIND=dp) :: res
824 :
825 : INTEGER :: i
826 : REAL(KIND=dp) :: s
827 :
828 : s = 0.0_dp
829 : DO i = 1, helium%beads
830 : s = s + helium_link_action(helium, i)
831 : END DO
832 : res = s
833 :
834 : END FUNCTION helium_total_action
835 :
836 : ! **************************************************************************************************
837 : !> \brief ...
838 : !> \param helium ...
839 : !> \param part ...
840 : !> \param ref_bead ...
841 : !> \param delta_bead ...
842 : !> \param d ...
843 : ! **************************************************************************************************
844 : SUBROUTINE helium_delta_pos(helium, part, ref_bead, delta_bead, d)
845 :
846 : TYPE(helium_solvent_type), INTENT(INOUT) :: helium
847 : INTEGER, INTENT(IN) :: part, ref_bead, delta_bead
848 : REAL(KIND=dp), DIMENSION(3), INTENT(OUT) :: d
849 :
850 : INTEGER :: b, bead, db, nbead, np, p
851 : REAL(KIND=dp), DIMENSION(3) :: r
852 :
853 : b = helium%beads
854 :
855 : d(:) = 0.0_dp
856 : IF (delta_bead > 0) THEN
857 : bead = ref_bead
858 : p = part
859 : db = delta_bead
860 : DO
861 : IF (db < 1) EXIT
862 : nbead = bead + 1
863 : np = p
864 : IF (nbead > b) THEN
865 : nbead = nbead - b
866 : np = helium%permutation(np)
867 : END IF
868 : r(:) = helium%pos(:, p, bead) - helium%pos(:, np, nbead)
869 : CALL helium_pbc(helium, r)
870 : d(:) = d(:) + r(:)
871 : bead = nbead
872 : p = np
873 : db = db - 1
874 : END DO
875 : ELSEIF (delta_bead < 0) THEN
876 : bead = ref_bead
877 : p = part
878 : db = delta_bead
879 : DO
880 : IF (db >= 0) EXIT
881 : nbead = bead - 1
882 : np = p
883 : IF (nbead < 1) THEN
884 : nbead = nbead + b
885 : np = helium%iperm(np)
886 : END IF
887 : r(:) = helium%pos(:, p, bead) - helium%pos(:, np, nbead)
888 : CALL helium_pbc(helium, r)
889 : d(:) = d(:) + r(:)
890 : bead = nbead
891 : p = np
892 : db = db + 1
893 : END DO
894 : END IF
895 : END SUBROUTINE helium_delta_pos
896 :
897 : #endif
898 :
899 : END MODULE helium_interactions
|
