(git:e7e05ae)
helium_worm.F
Go to the documentation of this file.
1 !--------------------------------------------------------------------------------------------------!
2 ! CP2K: A general program to perform molecular dynamics simulations !
3 ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 ! !
5 ! SPDX-License-Identifier: GPL-2.0-or-later !
6 !--------------------------------------------------------------------------------------------------!
7 
8 ! **************************************************************************************************
9 !> \brief Methods dealing with the canonical worm alogrithm
10 !> \author Felix Uhl [fuhl]
11 !> \date 2018-10-10
12 ! **************************************************************************************************
13 MODULE helium_worm
14 
15  USE helium_common, ONLY: helium_calc_plength,&
16  helium_eval_chain,&
17  helium_eval_expansion,&
18  helium_pbc,&
19  helium_update_coord_system
20  USE helium_interactions, ONLY: helium_bead_solute_e_f,&
21  helium_calc_energy,&
22  helium_solute_e_f
23  USE helium_io, ONLY: helium_write_line
24  USE helium_types, ONLY: helium_solvent_type
25  USE input_constants, ONLY: helium_forces_average,&
26  helium_forces_last
27  USE kinds, ONLY: default_string_length,&
28  dp
29  USE mathconstants, ONLY: pi
30  USE pint_types, ONLY: pint_env_type
31 #include "../base/base_uses.f90"
32 
33  IMPLICIT NONE
34 
35  PRIVATE
36  LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .true.
37  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'helium_worm'
38 
39  PUBLIC :: helium_sample_worm
40 
41 CONTAINS
42 
43 ! **************************************************************************************************
44 !> \brief Main worm sampling routine
45 !> \param helium ...
46 !> \param pint_env ...
47 !> \author fuhl
48 ! **************************************************************************************************
49  SUBROUTINE helium_sample_worm(helium, pint_env)
50 
51  TYPE(helium_solvent_type), INTENT(INOUT) :: helium
52  TYPE(pint_env_type), INTENT(IN) :: pint_env
53 
54  CHARACTER(len=default_string_length) :: stmp
55  INTEGER :: ac, iatom, ibead, icrawl, imc, imove, ncentracc, ncentratt, ncloseacc, ncloseatt, &
56  ncrawlbwdacc, ncrawlbwdatt, ncrawlfwdacc, ncrawlfwdatt, ncycle, nmc, nmoveheadacc, &
57  nmoveheadatt, nmovetailacc, nmovetailatt, nopenacc, nopenatt, npswapacc, nstagacc, &
58  nstagatt, nstat, nswapacc, nswapatt, ntot, staging_l
59  REAL(kind=dp) :: rtmp
60 
61  CALL helium_update_coord_system(helium, pint_env)
62 
63  ncentratt = 0
64  ncentracc = 0
65  nstagatt = 0
66  nstagacc = 0
67  nopenatt = 0
68  nopenacc = 0
69  ncloseatt = 0
70  ncloseacc = 0
71  nmoveheadatt = 0
72  nmoveheadacc = 0
73  nmovetailatt = 0
74  nmovetailacc = 0
75  ncrawlfwdatt = 0
76  ncrawlfwdacc = 0
77  ncrawlbwdatt = 0
78  ncrawlbwdacc = 0
79  nswapatt = 0
80  npswapacc = 0
81  nswapacc = 0
82  nstat = 0
83  ntot = 0
84  helium%proarea%inst(:) = 0.d0
85  helium%prarea2%inst(:) = 0.d0
86  helium%wnumber%inst(:) = 0.d0
87  helium%wnmber2%inst(:) = 0.d0
88  helium%mominer%inst(:) = 0.d0
89  IF (helium%solute_present) helium%force_avrg(:, :) = 0.0d0
90  helium%energy_avrg(:) = 0.0d0
91  helium%plength_avrg(:) = 0.0d0
92  IF (helium%worm_max_open_cycles > 0) THEN
93  helium%savepos(:, :, :) = helium%pos(:, :, :)
94  helium%saveiperm(:) = helium%iperm(:)
95  helium%savepermutation(:) = helium%permutation(:)
96  END IF
97 
98  nmc = helium%iter_rot
99  ncycle = 0
100  IF (helium%worm_allow_open) THEN
101  DO ! Exit criterion at the end of the loop
102  DO imc = 1, nmc
103  imove = helium%rng_stream_uniform%next(1, helium%worm_all_limit)
104  IF (helium%worm_is_closed) THEN
105  IF ((imove >= helium%worm_centroid_min) .AND. (imove <= helium%worm_centroid_max)) THEN
106  ! centroid move
107  iatom = helium%rng_stream_uniform%next(1, helium%atoms)
108  CALL worm_centroid_move(pint_env, helium, &
109  iatom, helium%worm_centroid_drmax, ac)
110  ncentratt = ncentratt + 1
111  ncentracc = ncentracc + ac
112  ! Note: weights for open and centroid move are taken from open sampling
113  ! staging is adjusted to conserve these weights
114  ELSE IF ((imove >= helium%worm_centroid_max + 1) .AND. (imove <= helium%worm_open_close_min - 1)) THEN
115  ! staging move
116  iatom = helium%rng_stream_uniform%next(1, helium%atoms)
117  ibead = helium%rng_stream_uniform%next(1, helium%beads)
118  staging_l = helium%rng_stream_uniform%next(2, helium%worm_staging_l)
119  CALL worm_staging_move(pint_env, helium, &
120  iatom, ibead, staging_l, ac)
121  nstagatt = nstagatt + 1
122  nstagacc = nstagacc + ac
123  ELSE IF ((imove >= helium%worm_open_close_min) .AND. (imove <= helium%worm_open_close_max)) THEN
124  ! attempt opening of worm
125  iatom = helium%rng_stream_uniform%next(1, helium%atoms)
126  ibead = helium%rng_stream_uniform%next(1, helium%beads)
127  staging_l = helium%rng_stream_uniform%next(2, helium%worm_staging_l)
128  CALL worm_open_move(pint_env, helium, &
129  iatom, ibead, staging_l, ac)
130  nopenatt = nopenatt + 1
131  nopenacc = nopenacc + ac
132  ELSE
133  ! this must not occur
134  cpabort("Undefined move selected in helium worm sampling!")
135  END IF
136  ELSE ! worm is open
137  IF ((imove >= helium%worm_centroid_min) .AND. (imove <= helium%worm_centroid_max)) THEN
138  ! centroid move
139  iatom = helium%rng_stream_uniform%next(1, helium%atoms)
140  CALL worm_centroid_move(pint_env, helium, &
141  iatom, helium%worm_centroid_drmax, ac)
142  ncentratt = ncentratt + 1
143  ncentracc = ncentracc + ac
144  ELSE IF ((imove >= helium%worm_staging_min) .AND. (imove <= helium%worm_staging_max)) THEN
145  ! staging move
146  iatom = helium%rng_stream_uniform%next(1, helium%atoms)
147  ibead = helium%rng_stream_uniform%next(1, helium%beads)
148  staging_l = helium%rng_stream_uniform%next(2, helium%worm_staging_l)
149  CALL worm_staging_move(pint_env, helium, &
150  iatom, ibead, staging_l, ac)
151  nstagatt = nstagatt + 1
152  nstagacc = nstagacc + ac
153  ELSE IF ((imove >= helium%worm_fcrawl_min) .AND. (imove <= helium%worm_fcrawl_max)) THEN
154  ! crawl forward
155  DO icrawl = 1, helium%worm_repeat_crawl
156  staging_l = helium%rng_stream_uniform%next(2, helium%worm_staging_l)
157  CALL worm_crawl_move_forward(pint_env, helium, &
158  staging_l, ac)
159  ncrawlfwdatt = ncrawlfwdatt + 1
160  ncrawlfwdacc = ncrawlfwdacc + ac
161  END DO
162  ELSE IF ((imove >= helium%worm_bcrawl_min) .AND. (imove <= helium%worm_bcrawl_max)) THEN
163  ! crawl backward
164  DO icrawl = 1, helium%worm_repeat_crawl
165  staging_l = helium%rng_stream_uniform%next(2, helium%worm_staging_l)
166  CALL worm_crawl_move_backward(pint_env, helium, &
167  staging_l, ac)
168  ncrawlbwdatt = ncrawlbwdatt + 1
169  ncrawlbwdacc = ncrawlbwdacc + ac
170  END DO
171  ELSE IF ((imove >= helium%worm_head_min) .AND. (imove <= helium%worm_head_max)) THEN
172  ! move head
173  staging_l = helium%rng_stream_uniform%next(2, helium%worm_staging_l)
174  CALL worm_head_move(pint_env, helium, &
175  staging_l, ac)
176  nmoveheadatt = nmoveheadatt + 1
177  nmoveheadacc = nmoveheadacc + ac
178  ELSE IF ((imove >= helium%worm_tail_min) .AND. (imove <= helium%worm_tail_max)) THEN
179  ! move tail
180  staging_l = helium%rng_stream_uniform%next(2, helium%worm_staging_l)
181  CALL worm_tail_move(pint_env, helium, &
182  staging_l, ac)
183  nmovetailatt = nmovetailatt + 1
184  nmovetailacc = nmovetailacc + ac
185  ELSE IF ((imove >= helium%worm_swap_min) .AND. (imove <= helium%worm_swap_max)) THEN
186  staging_l = helium%rng_stream_uniform%next(2, helium%worm_staging_l)
187  CALL worm_swap_move(pint_env, helium, &
188  helium%atoms, staging_l, ac)
189  npswapacc = npswapacc + ac
190  nswapacc = nswapacc + ac
191  nswapatt = nswapatt + 1
192  ELSE IF ((imove >= helium%worm_open_close_min) .AND. (imove <= helium%worm_open_close_max)) THEN
193  ! attempt closing of worm
194  staging_l = helium%rng_stream_uniform%next(2, helium%worm_staging_l)
195  CALL worm_close_move(pint_env, helium, &
196  staging_l, ac)
197  ncloseatt = ncloseatt + 1
198  ncloseacc = ncloseacc + ac
199  ELSE
200  ! this must not occur
201  cpabort("Undefined move selected in helium worm sampling!")
202  END IF
203  END IF
204 
205  ! Accumulate statistics if we are in the Z-sector:
206  IF (helium%worm_is_closed) THEN
207  nstat = nstat + 1
208  IF (helium%solute_present) THEN
209  IF (helium%get_helium_forces == helium_forces_average) THEN
210  !TODO needs proper averaging!
211  CALL helium_solute_e_f(pint_env, helium, rtmp)
212  helium%force_avrg(:, :) = helium%force_avrg(:, :) + helium%force_inst(:, :)
213  END IF
214  END IF
215  END IF
216  ntot = ntot + 1
217  END DO ! MC loop
218 
219  IF (helium%worm_is_closed) THEN
220  EXIT
221  ELSE
222  ncycle = ncycle + 1
223  ! reset position and permutation and start MC cycle again
224  ! if stuck in open position for too long
225  IF (helium%worm_max_open_cycles > 0 .AND. ncycle > helium%worm_max_open_cycles) THEN
226  nmc = helium%iter_rot
227  ncycle = 0
228  helium%pos = helium%savepos
229  helium%work = helium%pos
230  helium%permutation = helium%savepermutation
231  helium%iperm = helium%saveiperm
232  cpwarn("Trapped in open state, reset helium to closed state from last MD step.")
233  ELSE
234  nmc = max(helium%iter_rot/10, 10)
235  END IF
236  END IF
237  END DO !attempts loop
238  ELSE ! only closed configurations allowed
239  DO imc = 1, nmc
240  imove = helium%rng_stream_uniform%next(1, helium%worm_all_limit)
241 
242  IF ((imove >= helium%worm_centroid_min) .AND. (imove <= helium%worm_centroid_max)) THEN
243  ! centroid move
244  iatom = helium%rng_stream_uniform%next(1, helium%atoms)
245  CALL worm_centroid_move(pint_env, helium, &
246  iatom, helium%worm_centroid_drmax, ac)
247  ncentratt = ncentratt + 1
248  ncentracc = ncentracc + ac
249  ELSE IF ((imove >= helium%worm_staging_min) .AND. (imove <= helium%worm_staging_max)) THEN
250  ! staging move
251  iatom = helium%rng_stream_uniform%next(1, helium%atoms)
252  ibead = helium%rng_stream_uniform%next(1, helium%beads)
253  CALL worm_staging_move(pint_env, helium, &
254  iatom, ibead, helium%worm_staging_l, ac)
255  nstagatt = nstagatt + 1
256  nstagacc = nstagacc + ac
257  ELSE
258  ! this must not occur
259  cpabort("Undefined move selected in helium worm sampling!")
260  END IF
261 
262  ! Accumulate statistics if we are in closed configurations (which we always are)
263  nstat = nstat + 1
264  ntot = ntot + 1
265  IF (helium%solute_present) THEN
266  IF (helium%get_helium_forces == helium_forces_average) THEN
267  ! TODO: needs proper averaging
268  CALL helium_solute_e_f(pint_env, helium, rtmp)
269  helium%force_avrg(:, :) = helium%force_avrg(:, :) + helium%force_inst(:, :)
270  END IF
271  END IF
272  END DO ! MC loop
273  END IF
274 
275  ! Save naccepted and ntot
276  helium%num_accepted(1, 1) = ncentracc + nstagacc + nopenacc + ncloseacc + nswapacc + &
277  nmoveheadacc + nmovetailacc + ncrawlfwdacc + ncrawlbwdacc
278  helium%num_accepted(2, 1) = ncentratt + nstagatt + nopenatt + ncloseatt + nswapatt + &
279  nmoveheadatt + nmovetailatt + ncrawlfwdatt + ncrawlbwdatt
280  helium%worm_nstat = nstat
281 
282  ! Calculate energy and permutation path length
283  ! Multiply times helium%worm_nstat for consistent averaging in helium_sampling
284  CALL helium_calc_energy(helium, pint_env)
285  helium%energy_avrg(:) = helium%energy_inst(:)*real(helium%worm_nstat, dp)
286 
287  CALL helium_calc_plength(helium)
288  helium%plength_avrg(:) = helium%plength_inst(:)*real(helium%worm_nstat, dp)
289 
290  ! Calculate last_force
291  IF (helium%solute_present) THEN
292  IF (helium%get_helium_forces == helium_forces_last) THEN
293  CALL helium_solute_e_f(pint_env, helium, rtmp)
294  helium%force_avrg(:, :) = helium%force_avrg(:, :) + helium%force_inst(:, :)
295  END IF
296  END IF
297 
298  IF (helium%worm_show_statistics) THEN
299  WRITE (stmp, *) '--------------------------------------------------'
300  CALL helium_write_line(stmp)
301  WRITE (stmp, '(A10,F15.5,2I10)') 'Opening: ', &
302  REAL(nopenacc, dp)/REAL(MAX(1, nopenatt), dp), &
303  nopenacc, nopenatt
304  CALL helium_write_line(stmp)
305  WRITE (stmp, '(A10,F15.5,2I10)') 'Closing: ', &
306  REAL(ncloseacc, dp)/REAL(MAX(1, ncloseatt), dp), &
307  ncloseacc, ncloseatt
308  CALL helium_write_line(stmp)
309  WRITE (stmp, '(A10,F15.5,2I10)') 'Move Head: ', &
310  REAL(nmoveheadacc, dp)/REAL(MAX(1, nmoveheadatt), dp), &
311  nmoveheadacc, nmoveheadatt
312  CALL helium_write_line(stmp)
313  WRITE (stmp, '(A10,F15.5,2I10)') 'Move Tail: ', &
314  REAL(nmovetailacc, dp)/REAL(MAX(1, nmovetailatt), dp), &
315  nmovetailacc, nmovetailatt
316  CALL helium_write_line(stmp)
317  WRITE (stmp, '(A10,F15.5,2I10)') 'Crawl FWD: ', &
318  REAL(ncrawlfwdacc, dp)/REAL(MAX(1, ncrawlfwdatt), dp), &
319  ncrawlfwdacc, ncrawlfwdatt
320  CALL helium_write_line(stmp)
321  WRITE (stmp, '(A10,F15.5,2I10)') 'Crawl BWD: ', &
322  REAL(ncrawlbwdacc, dp)/REAL(MAX(1, ncrawlbwdatt), dp), &
323  ncrawlbwdacc, ncrawlbwdatt
324  CALL helium_write_line(stmp)
325  WRITE (stmp, '(A10,F15.5,2I10)') 'Staging: ', &
326  REAL(nstagacc, dp)/REAL(MAX(1, nstagatt), dp), &
327  nstagacc, nstagatt
328  CALL helium_write_line(stmp)
329  WRITE (stmp, '(A10,F15.5,2I10)') 'Centroid: ', &
330  REAL(ncentracc, dp)/REAL(MAX(1, ncentratt), dp), &
331  ncentracc, ncentratt
332  CALL helium_write_line(stmp)
333  WRITE (stmp, '(A10,F15.5,2I10)') 'Select: ', &
334  REAL(npswapacc, dp)/REAL(MAX(1, nswapatt), dp), &
335  npswapacc, nswapatt
336  CALL helium_write_line(stmp)
337  WRITE (stmp, '(A10,F15.5,2I10)') 'Swapping: ', &
338  REAL(nswapacc, dp)/REAL(MAX(1, nswapatt), dp), &
339  nswapacc, nswapatt
340  CALL helium_write_line(stmp)
341  WRITE (stmp, *) "Open State Probability: ", real(ntot - nstat, dp)/real(max(1, ntot), dp), ntot - nstat, ntot
342  CALL helium_write_line(stmp)
343  WRITE (stmp, *) "Closed State Probability: ", real(nstat, dp)/real(max(1, ntot), dp), nstat, ntot
344  CALL helium_write_line(stmp)
345  END IF
346 
347  !CALL center_pos(helium)
348 
349  END SUBROUTINE helium_sample_worm
350 
351 ! **************************************************************************************************
352 !> \brief Centroid Move (accounting for exchanged particles)
353 !> \param pint_env ...
354 !> \param helium ...
355 !> \param iatom ...
356 !> \param drmax ...
357 !> \param ac ...
358 !> \author fuhl
359 ! **************************************************************************************************
360  SUBROUTINE worm_centroid_move(pint_env, helium, iatom, drmax, ac)
361 
362  TYPE(pint_env_type), INTENT(IN) :: pint_env
363  TYPE(helium_solvent_type), INTENT(INOUT) :: helium
364  INTEGER, INTENT(IN) :: iatom
365  REAL(kind=dp), INTENT(IN) :: drmax
366  INTEGER, INTENT(OUT) :: ac
367 
368  INTEGER :: ia, ib, ic, jatom
369  LOGICAL :: should_reject, worm_in_moved_cycle
370  REAL(kind=dp) :: rtmp, rtmpo, sdiff, snew, sold
371  REAL(kind=dp), DIMENSION(3) :: dr, dro, new_com, old_com
372 
373  DO ic = 1, 3
374  rtmp = helium%rng_stream_uniform%next()
375  dr(ic) = (2.0_dp*rtmp - 1.0_dp)*drmax
376  END DO
377 
378  IF (helium%worm_is_closed) THEN
379  worm_in_moved_cycle = .false.
380  ! Perform move for first atom
381  DO ib = 1, helium%beads
382  helium%work(:, iatom, ib) = helium%work(:, iatom, ib) + dr(:)
383  END DO
384  ! move along permutation cycle
385  jatom = helium%permutation(iatom)
386  DO WHILE (jatom /= iatom)
387  DO ib = 1, helium%beads
388  helium%work(:, jatom, ib) = helium%work(:, jatom, ib) + dr(:)
389  END DO
390  ! next atom in chain
391  jatom = helium%permutation(jatom)
392  END DO
393  ELSE ! worm is open
394  worm_in_moved_cycle = (helium%worm_atom_idx == iatom)
395  ! while moving, check if worm is in moved cycle
396  ! Perform move for first atom
397  DO ib = 1, helium%beads
398  helium%work(:, iatom, ib) = helium%work(:, iatom, ib) + dr(:)
399  END DO
400  ! move along permutation cycle
401  jatom = helium%permutation(iatom)
402  DO WHILE (jatom /= iatom)
403  DO ib = 1, helium%beads
404  helium%work(:, jatom, ib) = helium%work(:, jatom, ib) + dr(:)
405  END DO
406  worm_in_moved_cycle = worm_in_moved_cycle .OR. (helium%worm_atom_idx == jatom)
407  ! next atom in chain
408  jatom = helium%permutation(jatom)
409  END DO
410  ! if atom contains had bead move that as well
411  IF (worm_in_moved_cycle) THEN
412  helium%worm_xtra_bead_work(:) = helium%worm_xtra_bead(:) + dr(:)
413  END IF
414  END IF
415 
416  sold = worm_centroid_move_action(helium, helium%pos, iatom, &
417  helium%worm_xtra_bead, worm_in_moved_cycle)
418 
419  snew = worm_centroid_move_action(helium, helium%work, iatom, &
420  helium%worm_xtra_bead_work, worm_in_moved_cycle)
421 
422  IF (helium%solute_present) THEN
423  sold = sold + worm_centroid_move_inter_action(pint_env, helium, helium%pos, iatom, &
424  helium%worm_xtra_bead, worm_in_moved_cycle)
425  snew = snew + worm_centroid_move_inter_action(pint_env, helium, helium%work, iatom, &
426  helium%worm_xtra_bead_work, worm_in_moved_cycle)
427  END IF
428 
429  ! Metropolis:
430  sdiff = sold - snew
431  IF (sdiff < 0) THEN
432  should_reject = .false.
433  IF (sdiff < -100.0_dp) THEN ! To protect from exponential underflow
434  should_reject = .true.
435  ELSE
436  rtmp = helium%rng_stream_uniform%next()
437  IF (exp(sdiff) < rtmp) THEN
438  should_reject = .true.
439  END IF
440  END IF
441  IF (should_reject) THEN
442  ! rejected !
443  ! write back only changed atoms
444  DO ib = 1, helium%beads
445  helium%work(:, iatom, ib) = helium%pos(:, iatom, ib)
446  END DO
447  ! move along permutation cycle
448  jatom = helium%permutation(iatom)
449  DO WHILE (jatom /= iatom)
450  DO ib = 1, helium%beads
451  helium%work(:, jatom, ib) = helium%pos(:, jatom, ib)
452  END DO
453  ! next atom in chain
454  jatom = helium%permutation(jatom)
455  END DO
456  helium%worm_xtra_bead_work(:) = helium%worm_xtra_bead(:)
457  ac = 0
458  RETURN
459  END IF
460  END IF
461 
462  ! for now accepted
463  ! rejection of the whole move if any centroid is farther away
464  ! from the current center of gravity than HELIUM%DROPLET_RADIUS [lwalewski]
465  ! AND ist not moving towards the center
466  IF (.NOT. helium%periodic) THEN
467  IF (helium%solute_present) THEN
468  new_com(:) = helium%center(:)
469  old_com(:) = new_com(:)
470  ELSE
471  new_com(:) = 0.0_dp
472  DO ia = 1, helium%atoms
473  DO ib = 1, helium%beads
474  new_com(:) = new_com(:) + helium%work(:, ia, ib)
475  END DO
476  END DO
477  new_com(:) = new_com(:)/(real(helium%atoms*helium%beads, dp))
478  ! also compute the old center of mass (optimization potential)
479  old_com(:) = 0.0_dp
480  DO ia = 1, helium%atoms
481  DO ib = 1, helium%beads
482  old_com(:) = old_com(:) + helium%pos(:, ia, ib)
483  END DO
484  END DO
485  old_com(:) = old_com(:)/(real(helium%atoms*helium%beads, dp))
486  END IF
487  should_reject = .false.
488  atom: DO ia = 1, helium%atoms
489  dr(:) = 0.0_dp
490  DO ib = 1, helium%beads
491  dr(:) = dr(:) + helium%work(:, ia, ib) - new_com(:)
492  END DO
493  dr(:) = dr(:)/real(helium%beads, dp)
494  rtmp = dot_product(dr, dr)
495  IF (rtmp >= helium%droplet_radius**2) THEN
496  dro(:) = 0.0_dp
497  DO ib = 1, helium%beads
498  dro(:) = dro(:) + helium%pos(:, ia, ib) - old_com(:)
499  END DO
500  dro(:) = dro(:)/real(helium%beads, dp)
501  rtmpo = dot_product(dro, dro)
502  ! only reject if it moves away from COM outside the droplet radius
503  IF (rtmpo < rtmp) THEN
504  ! found - this one does not fit within R from the center
505  should_reject = .true.
506  EXIT atom
507  END IF
508  END IF
509  END DO atom
510  IF (should_reject) THEN
511  ! restore original coordinates
512  ! write back only changed atoms
513  DO ib = 1, helium%beads
514  helium%work(:, iatom, ib) = helium%pos(:, iatom, ib)
515  END DO
516  ! move along permutation cycle
517  jatom = helium%permutation(iatom)
518  DO WHILE (jatom /= iatom)
519  DO ib = 1, helium%beads
520  helium%work(:, jatom, ib) = helium%pos(:, jatom, ib)
521  END DO
522  ! next atom in chain
523  jatom = helium%permutation(jatom)
524  END DO
525  helium%worm_xtra_bead_work(:) = helium%worm_xtra_bead(:)
526  ac = 0
527  RETURN
528  END IF
529  END IF
530 
531  ! finally accepted
532  ac = 1
533  ! write changed coordinates to position array
534  DO ib = 1, helium%beads
535  helium%pos(:, iatom, ib) = helium%work(:, iatom, ib)
536  END DO
537  ! move along permutation cycle
538  jatom = helium%permutation(iatom)
539  DO WHILE (jatom /= iatom)
540  DO ib = 1, helium%beads
541  helium%pos(:, jatom, ib) = helium%work(:, jatom, ib)
542  END DO
543  ! next atom in chain
544  jatom = helium%permutation(jatom)
545  END DO
546  helium%worm_xtra_bead(:) = helium%worm_xtra_bead_work(:)
547 
548  END SUBROUTINE worm_centroid_move
549 
550 ! **************************************************************************************************
551 !> \brief Action for centroid move
552 !> \param helium ...
553 !> \param pos ...
554 !> \param iatom ...
555 !> \param xtrapos ...
556 !> \param winc ...
557 !> \return ...
558 !> \author fuhl
559 ! **************************************************************************************************
560  REAL(kind=dp) FUNCTION worm_centroid_move_action(helium, pos, iatom, xtrapos, winc) &
561  result(partaction)
562 
563  TYPE(helium_solvent_type), INTENT(INOUT) :: helium
564  REAL(kind=dp), DIMENSION(:, :, :), INTENT(IN), &
565  POINTER :: pos
566  INTEGER, INTENT(IN) :: iatom
567  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: xtrapos
568  LOGICAL, INTENT(IN) :: winc
569 
570  INTEGER :: ia, ib, jatom, katom, opatom, patom, &
571  wbead
572  LOGICAL :: incycle
573  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: work2, work3
574  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: work
575  REAL(kind=dp), DIMENSION(3) :: r, rp
576 
577  ALLOCATE (work(3, helium%beads + 1), work2(helium%beads + 1), work3(SIZE(helium%uoffdiag, 1) + 1))
578 
579  ! compute action difference
580  ! action before the move from pos array
581  partaction = 0.0_dp
582 
583  ! pair of first atom with non-in-cycle atoms
584  jatom = iatom
585  DO ia = 1, helium%atoms
586  IF (ia == jatom) cycle
587  incycle = .false.
588  katom = helium%permutation(jatom)
589  DO WHILE (katom /= jatom)
590  IF (katom == ia) THEN
591  incycle = .true.
592  EXIT
593  END IF
594  katom = helium%permutation(katom)
595  END DO
596  IF (incycle) cycle
597  ! if not in cycle, compute pair action
598  DO ib = 1, helium%beads
599  work(:, ib) = pos(:, ia, ib) - pos(:, jatom, ib)
600  END DO
601  work(:, helium%beads + 1) = pos(:, helium%permutation(ia), 1) - &
602  pos(:, helium%permutation(jatom), 1)
603  partaction = partaction + helium_eval_chain(helium, work, helium%beads + 1, work2, work3)
604  END DO
605  ! all other cycle atoms with non-in-cycle atoms
606  jatom = helium%permutation(iatom)
607  DO WHILE (jatom /= iatom)
608  DO ia = 1, helium%atoms
609  IF (ia == jatom) cycle
610  incycle = .false.
611  katom = helium%permutation(jatom)
612  DO WHILE (katom /= jatom)
613  IF (katom == ia) THEN
614  incycle = .true.
615  EXIT
616  END IF
617  katom = helium%permutation(katom)
618  END DO
619  IF (incycle) cycle
620  ! if not in cycle, compute pair action
621  DO ib = 1, helium%beads
622  work(:, ib) = pos(:, ia, ib) - pos(:, jatom, ib)
623  END DO
624  work(:, helium%beads + 1) = pos(:, helium%permutation(ia), 1) - &
625  pos(:, helium%permutation(jatom), 1)
626  partaction = partaction + helium_eval_chain(helium, work, helium%beads + 1, work2, work3)
627  END DO
628  jatom = helium%permutation(jatom)
629  END DO
630  ! correct pair action for open worm configurations
631  IF (.NOT. helium%worm_is_closed) THEN
632  wbead = helium%worm_bead_idx
633  IF (winc) THEN
634  IF (helium%worm_bead_idx == 1) THEN
635  ! patom is the atom in front of the lone head bead
636  patom = helium%iperm(helium%worm_atom_idx)
637  ! go through all atoms, not in the cycle, and correct pair action
638  DO ia = 1, helium%atoms
639  IF (ia == helium%worm_atom_idx) cycle
640  incycle = .false.
641  katom = helium%permutation(helium%worm_atom_idx)
642  DO WHILE (katom /= helium%worm_atom_idx)
643  IF (katom == ia) THEN
644  incycle = .true.
645  EXIT
646  END IF
647  katom = helium%permutation(katom)
648  END DO
649  IF (incycle) cycle
650  ! find other previous atom
651  ! opatom is the atom in front of the first bead of the second atom
652  opatom = helium%iperm(ia)
653  ! if not in cycle, compute pair action
654  ! subtract pair action for closed link
655  r(:) = pos(:, helium%worm_atom_idx, 1) - pos(:, ia, 1)
656  rp(:) = pos(:, patom, helium%beads) - pos(:, opatom, helium%beads)
657  partaction = partaction - helium_eval_expansion(helium, r, rp, work3)
658  ! and add corrected extra link
659  ! rp stays the same
660  r(:) = xtrapos(:) - pos(:, ia, 1)
661  partaction = partaction + helium_eval_expansion(helium, r, rp, work3)
662  END DO
663  ELSE ! bead index /= 1
664  ! atom index is constant
665  ! go through all atoms, not in the cycle, and correct pair action
666  DO ia = 1, helium%atoms
667  IF (ia == helium%worm_atom_idx) cycle
668  incycle = .false.
669  katom = helium%permutation(helium%worm_atom_idx)
670  DO WHILE (katom /= helium%worm_atom_idx)
671  IF (katom == ia) THEN
672  incycle = .true.
673  EXIT
674  END IF
675  katom = helium%permutation(katom)
676  END DO
677  IF (incycle) cycle
678  ! if not in cycle, compute pair action
679  ! subtract pair action for closed link
680  r(:) = pos(:, helium%worm_atom_idx, wbead) - pos(:, ia, wbead)
681  rp(:) = pos(:, helium%worm_atom_idx, wbead - 1) - pos(:, ia, wbead - 1)
682  partaction = partaction - helium_eval_expansion(helium, r, rp, work3)
683  ! and add corrected extra link
684  ! rp stays the same
685  r(:) = xtrapos(:) - pos(:, ia, wbead)
686  partaction = partaction + helium_eval_expansion(helium, r, rp, work3)
687  END DO
688  END IF
689  ELSE ! worm is not the moved cycle
690  IF (helium%worm_bead_idx == 1) THEN
691  ! patom is the atom in front of the lone head bead
692  patom = helium%iperm(helium%worm_atom_idx)
693  opatom = helium%iperm(iatom)
694  !correct action contribution for first atom in moved cycle
695  r(:) = pos(:, helium%worm_atom_idx, 1) - pos(:, iatom, 1)
696  rp(:) = pos(:, patom, helium%beads) - pos(:, opatom, helium%beads)
697  partaction = partaction - helium_eval_expansion(helium, r, rp, work3)
698  ! and add corrected extra link
699  ! rp stays the same
700  r(:) = xtrapos(:) - pos(:, iatom, 1)
701  partaction = partaction + helium_eval_expansion(helium, r, rp, work3)
702  ! go through all other atoms, not in the exchange cycle, and correct pair action
703  ia = helium%permutation(iatom)
704  DO WHILE (ia /= iatom)
705  opatom = helium%iperm(ia)
706  ! subtract pair action for closed link
707  r(:) = pos(:, helium%worm_atom_idx, 1) - pos(:, ia, 1)
708  rp(:) = pos(:, patom, helium%beads) - pos(:, opatom, helium%beads)
709  partaction = partaction - helium_eval_expansion(helium, r, rp, work3)
710  ! and add corrected extra link
711  ! rp stays the same
712  r(:) = xtrapos(:) - pos(:, ia, 1)
713  partaction = partaction + helium_eval_expansion(helium, r, rp, work3)
714  ia = helium%permutation(ia)
715  END DO
716  ELSE ! bead index /= 1
717  ! patom is the atom in front of the lone head bead
718  !correct action contribution for first atom in moved cycle
719  r(:) = pos(:, helium%worm_atom_idx, wbead) - pos(:, iatom, wbead)
720  rp(:) = pos(:, helium%worm_atom_idx, wbead - 1) - pos(:, iatom, wbead - 1)
721  partaction = partaction - helium_eval_expansion(helium, r, rp, work3)
722  ! and add corrected extra link
723  ! rp stays the same
724  r(:) = xtrapos(:) - pos(:, iatom, wbead)
725  partaction = partaction + helium_eval_expansion(helium, r, rp, work3)
726  ! go through all other atoms, not in the exchange cycle, and correct pair action
727  ia = helium%permutation(iatom)
728  DO WHILE (ia /= iatom)
729  ! subtract pair action for closed link
730  r(:) = pos(:, helium%worm_atom_idx, wbead) - pos(:, ia, wbead)
731  rp(:) = pos(:, helium%worm_atom_idx, wbead - 1) - pos(:, ia, wbead - 1)
732  partaction = partaction - helium_eval_expansion(helium, r, rp, work3)
733  ! and add corrected extra link
734  ! rp stays the same
735  r(:) = xtrapos(:) - pos(:, ia, wbead)
736  partaction = partaction + helium_eval_expansion(helium, r, rp, work3)
737  ia = helium%permutation(ia)
738  END DO
739  END IF
740  END IF
741  END IF
742  DEALLOCATE (work, work2, work3)
743 
744  END FUNCTION worm_centroid_move_action
745 
746 ! **************************************************************************************************
747 !> \brief Centroid interaction
748 !> \param pint_env ...
749 !> \param helium ...
750 !> \param pos ...
751 !> \param iatom ...
752 !> \param xtrapos ...
753 !> \param winc ...
754 !> \return ...
755 !> \author fuhl
756 ! **************************************************************************************************
757  REAL(kind=dp) FUNCTION worm_centroid_move_inter_action(pint_env, helium, pos, iatom, xtrapos, winc) &
758  result(partaction)
759 
760  TYPE(pint_env_type), INTENT(IN) :: pint_env
761  TYPE(helium_solvent_type), INTENT(IN) :: helium
762  REAL(kind=dp), DIMENSION(:, :, :), INTENT(IN), &
763  POINTER :: pos
764  INTEGER, INTENT(IN) :: iatom
765  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: xtrapos
766  LOGICAL, INTENT(IN) :: winc
767 
768  INTEGER :: jatom, jbead
769  REAL(kind=dp) :: energy
770 
771  partaction = 0.0_dp
772  ! spcial treatment if worm is in moved cycle
773  IF (winc) THEN
774  ! first atom by hand
775  jatom = iatom
776  ! if it is worm atom it gets special treatment
777  IF (jatom == helium%worm_atom_idx) THEN
778  ! up to worm intersection
779  DO jbead = 1, helium%worm_bead_idx - 1
780  CALL helium_bead_solute_e_f(pint_env, helium, &
781  jatom, jbead, pos(:, jatom, jbead), energy=energy)
782  partaction = partaction + energy
783  END DO
784  ! head and tail each with 1/2 weight
785  jbead = helium%worm_bead_idx
786  ! tail
787  CALL helium_bead_solute_e_f(pint_env, helium, &
788  jatom, jbead, pos(:, jatom, jbead), energy=energy)
789  partaction = partaction + 0.5_dp*energy
790  ! head
791  CALL helium_bead_solute_e_f(pint_env, helium, &
792  jatom, jbead, xtrapos, energy=energy)
793  partaction = partaction + 0.5_dp*energy
794  ! rest of ring polymer
795  DO jbead = helium%worm_bead_idx + 1, helium%beads
796  CALL helium_bead_solute_e_f(pint_env, helium, &
797  jatom, jbead, pos(:, jatom, jbead), energy=energy)
798  partaction = partaction + energy
799  END DO
800  ELSE
801  DO jbead = 1, helium%beads
802  CALL helium_bead_solute_e_f(pint_env, helium, &
803  jatom, jbead, pos(:, jatom, jbead), energy=energy)
804  partaction = partaction + energy
805  END DO
806  END IF
807  jatom = helium%permutation(iatom)
808  DO WHILE (jatom /= iatom)
809  IF (jatom == helium%worm_atom_idx) THEN
810  ! up to worm intersection
811  DO jbead = 1, helium%worm_bead_idx - 1
812  CALL helium_bead_solute_e_f(pint_env, helium, &
813  jatom, jbead, pos(:, jatom, jbead), energy=energy)
814  partaction = partaction + energy
815  END DO
816  ! head and tail each with 1/2 weight
817  jbead = helium%worm_bead_idx
818  ! tail
819  CALL helium_bead_solute_e_f(pint_env, helium, &
820  jatom, jbead, pos(:, jatom, jbead), energy=energy)
821  partaction = partaction + 0.5_dp*energy
822  ! head
823  CALL helium_bead_solute_e_f(pint_env, helium, &
824  jatom, jbead, xtrapos, energy=energy)
825  partaction = partaction + 0.5_dp*energy
826  ! rest of ring polymer
827  DO jbead = helium%worm_bead_idx + 1, helium%beads
828  CALL helium_bead_solute_e_f(pint_env, helium, &
829  jatom, jbead, pos(:, jatom, jbead), energy=energy)
830  partaction = partaction + energy
831  END DO
832  ELSE
833  DO jbead = 1, helium%beads
834  CALL helium_bead_solute_e_f(pint_env, helium, &
835  jatom, jbead, pos(:, jatom, jbead), energy=energy)
836  partaction = partaction + energy
837  END DO
838  END IF
839  jatom = helium%permutation(jatom)
840  END DO
841  ELSE ! worm not in moved cycle
842  ! first atom by hand
843  jatom = iatom
844  ! if it is worm atom it gets special treatment
845  DO jbead = 1, helium%beads
846  CALL helium_bead_solute_e_f(pint_env, helium, &
847  jatom, jbead, pos(:, jatom, jbead), energy=energy)
848  partaction = partaction + energy
849  END DO
850  jatom = helium%permutation(iatom)
851  DO WHILE (jatom /= iatom)
852  DO jbead = 1, helium%beads
853  CALL helium_bead_solute_e_f(pint_env, helium, &
854  jatom, jbead, pos(:, jatom, jbead), energy=energy)
855  partaction = partaction + energy
856  END DO
857  jatom = helium%permutation(jatom)
858  END DO
859  END IF
860 
861  partaction = partaction*helium%tau
862 
863  END FUNCTION worm_centroid_move_inter_action
864 
865 ! **************************************************************************************************
866 !> \brief General path construct based on staging moves
867 !> \param helium ...
868 !> \param ri ...
869 !> \param rj ...
870 !> \param l ...
871 !> \param new_path ...
872 !> \author fuhl
873 ! **************************************************************************************************
874  SUBROUTINE path_construct(helium, ri, rj, l, new_path)
875 
876  !constructs a path of length l between the positions ri and rj
877  TYPE(helium_solvent_type), INTENT(INOUT) :: helium
878  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: ri, rj
879  INTEGER, INTENT(IN) :: l
880  REAL(kind=dp), DIMENSION(3, l), INTENT(OUT) :: new_path
881 
882  INTEGER :: idim, istage
883  REAL(kind=dp) :: imass, invstagemass, rk, weight
884  REAL(kind=dp), DIMENSION(3) :: re, rs
885 
886  imass = 1.0_dp/helium%he_mass_au
887  ! dealing with periodicity
888  rs(:) = ri(:)
889  re(:) = rj(:) - rs(:)
890  CALL helium_pbc(helium, re)
891  re(:) = re(:) + rs(:)
892 
893  ! first construction by hand
894  ! reusable weight factor 1/(l+1)
895  rk = real(l, dp)
896  weight = 1.0_dp/(rk + 1.0_dp)
897  ! staging mass needed for modified variance
898  invstagemass = rk*weight*imass
899  ! proposing new positions
900  DO idim = 1, 3
901  new_path(idim, 1) = helium%rng_stream_gaussian%next(variance=helium%tau*invstagemass)
902  END DO
903  new_path(:, 1) = new_path(:, 1) + weight*(re(:) + rk*rs(:))
904 
905  DO istage = 2, l
906  ! reusable weight factor 1/(k+1)
907  rk = real(l - istage + 1, dp)
908  weight = 1.0_dp/(rk + 1.0_dp)
909  ! staging mass needed for modified variance
910  invstagemass = rk*weight*imass
911  ! proposing new positions
912  DO idim = 1, 3
913  new_path(idim, istage) = helium%rng_stream_gaussian%next(variance=helium%tau*invstagemass)
914  END DO
915  new_path(:, istage) = new_path(:, istage) + weight*(rk*new_path(:, istage - 1) + re(:))
916  END DO
917 
918  END SUBROUTINE path_construct
919 
920 ! **************************************************************************************************
921 !> \brief Staging move
922 !> \param pint_env ...
923 !> \param helium ...
924 !> \param startatom ...
925 !> \param startbead ...
926 !> \param l ...
927 !> \param ac ...
928 !> \author fuhl
929 ! **************************************************************************************************
930  SUBROUTINE worm_staging_move(pint_env, helium, startatom, startbead, l, ac)
931 
932  TYPE(pint_env_type), INTENT(IN) :: pint_env
933  TYPE(helium_solvent_type), INTENT(INOUT) :: helium
934  INTEGER, INTENT(IN) :: startatom, startbead, l
935  INTEGER, INTENT(OUT) :: ac
936 
937  INTEGER :: endatom, endbead, ia, ib, ibead, jbead
938  LOGICAL :: should_reject, worm_overlap
939  REAL(kind=dp) :: rtmp, rtmpo, sdiff, snew, sold
940  REAL(kind=dp), DIMENSION(3) :: dr, dro, new_com, old_com
941  REAL(kind=dp), DIMENSION(3, l) :: newsection
942 
943  ac = 0
944  endbead = startbead + l + 1
945  ! Check if the imaginary time section belongs to two atoms
946  IF (endbead > helium%beads) THEN
947  endatom = helium%permutation(startatom)
948  endbead = endbead - helium%beads
949  ELSE
950  endatom = startatom
951  END IF
952 
953  ! check if the imaginary time section overlaps with the worm opening
954  IF (helium%worm_is_closed) THEN
955  worm_overlap = .false.
956  ELSE
957  ! if it does give it special treatment during action evaluation
958  worm_overlap = ((startbead < endbead) .AND. &
959  (helium%worm_bead_idx > startbead) .AND. &
960  (helium%worm_bead_idx <= endbead)) &
961  .OR. &
962  ((startbead > endbead) .AND. &
963  ((helium%worm_bead_idx > startbead) .OR. &
964  (helium%worm_bead_idx <= endbead)))
965  IF (worm_overlap) THEN
966  ! if in addition the worm end is IN the reconstructed section reject immediately
967  IF (((helium%worm_atom_idx == startatom) .AND. (helium%worm_bead_idx >= startbead)) .OR. &
968  ((helium%worm_atom_idx == endatom) .AND. (helium%worm_bead_idx <= endbead))) THEN
969  ac = 0
970  RETURN
971  END IF
972  END IF
973  END IF
974  ! action before the move
975  IF (worm_overlap) THEN
976  sold = worm_path_action_worm_corrected(helium, helium%pos, &
977  startatom, startbead, endatom, endbead, &
978  helium%worm_xtra_bead, helium%worm_atom_idx, helium%worm_bead_idx)
979  ELSE
980  sold = worm_path_action(helium, helium%pos, &
981  startatom, startbead, endatom, endbead)
982  END IF
983 
984  IF (helium%solute_present) THEN
985  ! no special head treatment needed, because a swap can't go over
986  ! the worm gap and due to primitive coupling no cross bead terms are considered
987  sold = sold + worm_path_inter_action(pint_env, helium, helium%pos, &
988  startatom, startbead, endatom, endbead)
989  END IF
990 
991  ! construct a new path connecting the start and endbead
992  CALL path_construct(helium, &
993  helium%pos(:, startatom, startbead), &
994  helium%pos(:, endatom, endbead), l, &
995  newsection)
996 
997  ! write new path segment to work array
998  ! first the part that is guaranteed to fit on the coorinates of startatom
999  jbead = 1
1000  DO ibead = startbead + 1, min(helium%beads, startbead + l)
1001  helium%work(:, startatom, ibead) = newsection(:, jbead)
1002  jbead = jbead + 1
1003  END DO
1004  ! transfer the rest of the beads to coordinates of endatom if necessary
1005  IF (helium%beads < startbead + l) THEN
1006  DO ibead = 1, endbead - 1
1007  helium%work(:, endatom, ibead) = newsection(:, jbead)
1008  jbead = jbead + 1
1009  END DO
1010  END IF
1011 
1012  ! action after the move
1013  IF (worm_overlap) THEN
1014  snew = worm_path_action_worm_corrected(helium, helium%work, &
1015  startatom, startbead, endatom, endbead, &
1016  helium%worm_xtra_bead_work, helium%worm_atom_idx_work, helium%worm_bead_idx_work)
1017  ELSE
1018  snew = worm_path_action(helium, helium%work, &
1019  startatom, startbead, endatom, endbead)
1020  END IF
1021 
1022  IF (helium%solute_present) THEN
1023  ! no special head treatment needed, because a swap can't go over
1024  ! the worm gap and due to primitive coupling no cross bead terms are considered
1025  snew = snew + worm_path_inter_action(pint_env, helium, helium%work, &
1026  startatom, startbead, endatom, endbead)
1027  END IF
1028 
1029  ! Metropolis:
1030  sdiff = sold - snew
1031  IF (sdiff < 0) THEN
1032  should_reject = .false.
1033  IF (sdiff < -100.0_dp) THEN ! To protect from exponential underflow
1034  should_reject = .true.
1035  ELSE
1036  rtmp = helium%rng_stream_uniform%next()
1037  IF (exp(sdiff) < rtmp) THEN
1038  should_reject = .true.
1039  END IF
1040  END IF
1041  IF (should_reject) THEN
1042  ! rejected !
1043  ! write back only changed atoms
1044  jbead = 1
1045  DO ibead = startbead + 1, min(helium%beads, startbead + l)
1046  helium%work(:, startatom, ibead) = helium%pos(:, startatom, ibead)
1047  jbead = jbead + 1
1048  END DO
1049  ! transfer the rest of the beads to coordinates of endatom if necessary
1050  IF (helium%beads < startbead + l) THEN
1051  DO ibead = 1, endbead - 1
1052  helium%work(:, endatom, ibead) = helium%pos(:, endatom, ibead)
1053  jbead = jbead + 1
1054  END DO
1055  END IF
1056  ac = 0
1057  RETURN
1058  END IF
1059  END IF
1060 
1061  ! for now accepted
1062  ! rejection of the whole move if any centroid is farther away
1063  ! from the current center of gravity than HELIUM%DROPLET_RADIUS [lwalewski]
1064  ! AND ist not moving towards the center
1065  IF (.NOT. helium%periodic) THEN
1066  IF (helium%solute_present) THEN
1067  new_com(:) = helium%center(:)
1068  old_com(:) = new_com(:)
1069  ELSE
1070  new_com(:) = 0.0_dp
1071  DO ia = 1, helium%atoms
1072  DO ib = 1, helium%beads
1073  new_com(:) = new_com(:) + helium%work(:, ia, ib)
1074  END DO
1075  END DO
1076  new_com(:) = new_com(:)/(real(helium%atoms*helium%beads, dp))
1077  ! also compute the old center of mass (optimization potential)
1078  old_com(:) = 0.0_dp
1079  DO ia = 1, helium%atoms
1080  DO ib = 1, helium%beads
1081  old_com(:) = old_com(:) + helium%pos(:, ia, ib)
1082  END DO
1083  END DO
1084  old_com(:) = old_com(:)/(real(helium%atoms*helium%beads, dp))
1085  END IF
1086  should_reject = .false.
1087  atom: DO ia = 1, helium%atoms
1088  dr(:) = 0.0_dp
1089  DO ib = 1, helium%beads
1090  dr(:) = dr(:) + helium%work(:, ia, ib) - new_com(:)
1091  END DO
1092  dr(:) = dr(:)/real(helium%beads, dp)
1093  rtmp = dot_product(dr, dr)
1094  IF (rtmp >= helium%droplet_radius**2) THEN
1095  dro(:) = 0.0_dp
1096  DO ib = 1, helium%beads
1097  dro(:) = dro(:) + helium%pos(:, ia, ib) - old_com(:)
1098  END DO
1099  dro(:) = dro(:)/real(helium%beads, dp)
1100  rtmpo = dot_product(dro, dro)
1101  ! only reject if it moves away from COM outside the droplet radius
1102  IF (rtmpo < rtmp) THEN
1103  ! found - this one does not fit within R from the center
1104  should_reject = .true.
1105  EXIT atom
1106  END IF
1107  END IF
1108  END DO atom
1109  IF (should_reject) THEN
1110  ! restore original coordinates
1111  ! write back only changed atoms
1112  jbead = 1
1113  DO ibead = startbead + 1, min(helium%beads, startbead + l)
1114  helium%work(:, startatom, ibead) = helium%pos(:, startatom, ibead)
1115  jbead = jbead + 1
1116  END DO
1117  ! transfer the rest of the beads to coordinates of endatom if necessary
1118  IF (helium%beads < startbead + l) THEN
1119  DO ibead = 1, endbead - 1
1120  helium%work(:, endatom, ibead) = helium%pos(:, endatom, ibead)
1121  jbead = jbead + 1
1122  END DO
1123  END IF
1124  ac = 0
1125  RETURN
1126  END IF
1127  END IF
1128 
1129  ! finally accepted
1130  ac = 1
1131  ! write changed coordinates to position array
1132  jbead = 1
1133  DO ibead = startbead + 1, min(helium%beads, startbead + l)
1134  helium%pos(:, startatom, ibead) = helium%work(:, startatom, ibead)
1135  jbead = jbead + 1
1136  END DO
1137  ! transfer the rest of the beads to coordinates of endatom if necessary
1138  IF (helium%beads < startbead + l) THEN
1139  DO ibead = 1, endbead - 1
1140  helium%pos(:, endatom, ibead) = helium%work(:, endatom, ibead)
1141  jbead = jbead + 1
1142  END DO
1143  END IF
1144 
1145  END SUBROUTINE worm_staging_move
1146 
1147 ! **************************************************************************************************
1148 !> \brief Open move to off-diagonal elements of the density matrix, allows to sample permutations
1149 !> \param pint_env ...
1150 !> \param helium ...
1151 !> \param endatom ...
1152 !> \param endbead ...
1153 !> \param l ...
1154 !> \param ac ...
1155 !> \author fuhl
1156 ! **************************************************************************************************
1157  SUBROUTINE worm_open_move(pint_env, helium, endatom, endbead, l, ac)
1158 
1159  TYPE(pint_env_type), INTENT(IN) :: pint_env
1160  TYPE(helium_solvent_type), INTENT(INOUT) :: helium
1161  INTEGER, INTENT(IN) :: endatom, endbead, l
1162  INTEGER, INTENT(OUT) :: ac
1163 
1164  INTEGER :: ia, ib, idim, kbead, startatom, startbead
1165  LOGICAL :: should_reject
1166  REAL(kind=dp) :: distsq, mass, rtmp, rtmpo, sdiff, snew, &
1167  sold, xr
1168  REAL(kind=dp), DIMENSION(3) :: distvec, dr, dro, new_com, old_com
1169 
1170  mass = helium%he_mass_au
1171 
1172  ! get index of the atom and bead, where the resampling of the head begins
1173  IF (l < endbead) THEN
1174  ! startbead belongs to the same atom
1175  startatom = endatom
1176  startbead = endbead - l
1177  ELSE
1178  ! startbead belongs to a different atom
1179  ! find previous atom (assuming l < nbeads)
1180  startatom = helium%iperm(endatom)
1181  startbead = endbead + helium%beads - l
1182  END IF
1183  sold = worm_path_action(helium, helium%pos, &
1184  startatom, startbead, endatom, endbead)
1185 
1186  IF (helium%solute_present) THEN
1187  ! yes this is correct, as the bead, that splits into tail and head only changes half
1188  ! therefore only half of its action needs to be considred
1189  ! and is cheated in here by passing it as head bead
1190  sold = sold + worm_path_inter_action_head(pint_env, helium, helium%pos, &
1191  startatom, startbead, &
1192  helium%pos(:, endatom, endbead), endatom, endbead)
1193  END IF
1194 
1195  helium%worm_is_closed = .false.
1196  helium%worm_atom_idx_work = endatom
1197  helium%worm_bead_idx_work = endbead
1198 
1199  ! alternative grow with consecutive gaussians
1200  IF (startbead < endbead) THEN
1201  ! everything belongs to the same atom
1202  ! gro head from startbead
1203  DO kbead = startbead + 1, endbead - 1
1204  DO idim = 1, 3
1205  xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
1206  helium%work(idim, startatom, kbead) = helium%work(idim, startatom, kbead - 1) + xr
1207  END DO
1208  END DO
1209  ! last grow head bead
1210  DO idim = 1, 3
1211  xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
1212  helium%worm_xtra_bead_work(idim) = helium%work(idim, startatom, endbead - 1) + xr
1213  END DO
1214  ELSE IF (endbead /= 1) THEN
1215  ! is distributed among two atoms
1216  ! grow from startbead
1217  DO kbead = startbead + 1, helium%beads
1218  DO idim = 1, 3
1219  xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
1220  helium%work(idim, startatom, kbead) = helium%work(idim, startatom, kbead - 1) + xr
1221  END DO
1222  END DO
1223  ! bead one of endatom relative to last on startatom
1224  DO idim = 1, 3
1225  xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
1226  helium%work(idim, endatom, 1) = helium%work(idim, startatom, helium%beads) + xr
1227  END DO
1228  ! everything on endatom
1229  DO kbead = 2, endbead - 1
1230  DO idim = 1, 3
1231  xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
1232  helium%work(idim, endatom, kbead) = helium%work(idim, endatom, kbead - 1) + xr
1233  END DO
1234  END DO
1235  DO idim = 1, 3
1236  xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
1237  helium%worm_xtra_bead_work(idim) = helium%work(idim, endatom, endbead - 1) + xr
1238  END DO
1239  ELSE ! imagtimewrap and headbead = 1
1240  ! is distributed among two atoms
1241  ! grow from startbead
1242  DO kbead = startbead + 1, helium%beads
1243  DO idim = 1, 3
1244  xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
1245  helium%work(idim, startatom, kbead) = helium%work(idim, startatom, kbead - 1) + xr
1246  END DO
1247  END DO
1248  ! bead one of endatom relative to last on startatom
1249  DO idim = 1, 3
1250  xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
1251  helium%worm_xtra_bead_work(idim) = helium%work(idim, startatom, helium%beads) + xr
1252  END DO
1253  END IF
1254 
1255  snew = worm_path_action_worm_corrected(helium, helium%work, &
1256  startatom, startbead, endatom, endbead, &
1257  helium%worm_xtra_bead_work, helium%worm_atom_idx_work, helium%worm_bead_idx_work)
1258 
1259  IF (helium%solute_present) THEN
1260  snew = snew + worm_path_inter_action_head(pint_env, helium, helium%work, &
1261  startatom, startbead, &
1262  helium%worm_xtra_bead_work, helium%worm_atom_idx_work, helium%worm_bead_idx_work)
1263  END IF
1264 
1265  ! Metropolis:
1266  ! first compute ln of free density matrix
1267  distvec(:) = helium%pos(:, startatom, startbead) - helium%pos(:, endatom, endbead)
1268  CALL helium_pbc(helium, distvec)
1269  distsq = dot_product(distvec, distvec)
1270  ! action difference
1271  sdiff = sold - snew
1272  ! modify action difference due to extra bead
1273  sdiff = sdiff + distsq/(2.0_dp*helium%hb2m*real(l, dp)*helium%tau)
1274  sdiff = sdiff + 1.5_dp*log(real(l, dp)*helium%tau)
1275  sdiff = sdiff + helium%worm_ln_openclose_scale
1276  IF (sdiff < 0) THEN
1277  should_reject = .false.
1278  IF (sdiff < -100.0_dp) THEN ! To protect from exponential underflow
1279  should_reject = .true.
1280  ELSE
1281  rtmp = helium%rng_stream_uniform%next()
1282  IF (exp(sdiff) < rtmp) THEN
1283  should_reject = .true.
1284  END IF
1285  END IF
1286  IF (should_reject) THEN
1287  ! rejected !
1288  ! write back only changed atoms
1289  ! transfer the new coordinates to work array
1290  IF (startbead < endbead) THEN
1291  ! everything belongs to the same atom
1292  DO kbead = startbead + 1, endbead - 1
1293  helium%work(:, endatom, kbead) = helium%pos(:, endatom, kbead)
1294  END DO
1295  ELSE
1296  ! is distributed among two atoms
1297  ! transfer to atom not containing the head bead
1298  DO kbead = startbead + 1, helium%beads
1299  helium%work(:, startatom, kbead) = helium%pos(:, startatom, kbead)
1300  END DO
1301  ! transfer to atom containing the head bead
1302  DO kbead = 1, endbead - 1
1303  helium%work(:, endatom, kbead) = helium%pos(:, endatom, kbead)
1304  END DO
1305  END IF
1306  helium%worm_is_closed = .true.
1307  ac = 0
1308  RETURN
1309  END IF
1310  END IF
1311 
1312  ! for now accepted
1313  ! rejection of the whole move if any centroid is farther away
1314  ! from the current center of gravity than HELIUM%DROPLET_RADIUS [lwalewski]
1315  ! AND ist not moving towards the center
1316  IF (.NOT. helium%periodic) THEN
1317  IF (helium%solute_present) THEN
1318  new_com(:) = helium%center(:)
1319  old_com(:) = new_com(:)
1320  ELSE
1321  new_com(:) = 0.0_dp
1322  DO ia = 1, helium%atoms
1323  DO ib = 1, helium%beads
1324  new_com(:) = new_com(:) + helium%work(:, ia, ib)
1325  END DO
1326  END DO
1327  new_com(:) = new_com(:)/(real(helium%atoms*helium%beads, dp))
1328  ! also compute the old center of mass (optimization potential)
1329  old_com(:) = 0.0_dp
1330  DO ia = 1, helium%atoms
1331  DO ib = 1, helium%beads
1332  old_com(:) = old_com(:) + helium%pos(:, ia, ib)
1333  END DO
1334  END DO
1335  old_com(:) = old_com(:)/(real(helium%atoms*helium%beads, dp))
1336  END IF
1337  should_reject = .false.
1338  atom: DO ia = 1, helium%atoms
1339  dr(:) = 0.0_dp
1340  DO ib = 1, helium%beads
1341  dr(:) = dr(:) + helium%work(:, ia, ib) - new_com(:)
1342  END DO
1343  dr(:) = dr(:)/real(helium%beads, dp)
1344  rtmp = dot_product(dr, dr)
1345  IF (rtmp >= helium%droplet_radius**2) THEN
1346  dro(:) = 0.0_dp
1347  DO ib = 1, helium%beads
1348  dro(:) = dro(:) + helium%pos(:, ia, ib) - old_com(:)
1349  END DO
1350  dro(:) = dro(:)/real(helium%beads, dp)
1351  rtmpo = dot_product(dro, dro)
1352  ! only reject if it moves away from COM outside the droplet radius
1353  IF (rtmpo < rtmp) THEN
1354  ! found - this one does not fit within R from the center
1355  should_reject = .true.
1356  EXIT atom
1357  END IF
1358  END IF
1359  END DO atom
1360  IF (should_reject) THEN
1361  ! restore original coordinates
1362  ! write back only changed atoms
1363  IF (startbead < endbead) THEN
1364  ! everything belongs to the same atom
1365  DO kbead = startbead + 1, endbead - 1
1366  helium%work(:, endatom, kbead) = helium%pos(:, endatom, kbead)
1367  END DO
1368  ELSE
1369  ! is distributed among two atoms
1370  ! transfer to atom not containing the head bead
1371  DO kbead = startbead + 1, helium%beads
1372  helium%work(:, startatom, kbead) = helium%pos(:, startatom, kbead)
1373  END DO
1374  ! transfer to atom containing the head bead
1375  DO kbead = 1, endbead - 1
1376  helium%work(:, endatom, kbead) = helium%pos(:, endatom, kbead)
1377  END DO
1378  END IF
1379  helium%worm_is_closed = .true.
1380  !helium%worm_atom_idx_work = helium%worm_atom_idx
1381  !helium%worm_bead_idx_work = helium%worm_bead_idx
1382  ac = 0
1383  RETURN
1384  END IF
1385  END IF
1386 
1387  ! finally accepted
1388  ac = 1
1389  ! write changed coordinates to position array
1390  IF (startbead < endbead) THEN
1391  ! everything belongs to the same atom
1392  DO kbead = startbead + 1, endbead - 1
1393  helium%pos(:, endatom, kbead) = helium%work(:, endatom, kbead)
1394  END DO
1395  ELSE
1396  ! is distributed among two atoms
1397  ! transfer to atom not containing the head bead
1398  DO kbead = startbead + 1, helium%beads
1399  helium%pos(:, startatom, kbead) = helium%work(:, startatom, kbead)
1400  END DO
1401  ! transfer to atom containing the head bead
1402  DO kbead = 1, endbead - 1
1403  helium%pos(:, endatom, kbead) = helium%work(:, endatom, kbead)
1404  END DO
1405  END IF
1406  helium%worm_xtra_bead(:) = helium%worm_xtra_bead_work(:)
1407  helium%worm_atom_idx = helium%worm_atom_idx_work
1408  helium%worm_bead_idx = helium%worm_bead_idx_work
1409 
1410  END SUBROUTINE worm_open_move
1411 
1412 ! **************************************************************************************************
1413 !> \brief Close move back to diagonal elements of density matrix, permutation fixed in closed state
1414 !> \param pint_env ...
1415 !> \param helium ...
1416 !> \param l ...
1417 !> \param ac ...
1418 !> \author fuhl
1419 ! **************************************************************************************************
1420  SUBROUTINE worm_close_move(pint_env, helium, l, ac)
1421 
1422  TYPE(pint_env_type), INTENT(IN) :: pint_env
1423  TYPE(helium_solvent_type), INTENT(INOUT) :: helium
1424  INTEGER, INTENT(IN) :: l
1425  INTEGER, INTENT(OUT) :: ac
1426 
1427  INTEGER :: endatom, endbead, ia, ib, jbead, kbead, &
1428  startatom, startbead
1429  LOGICAL :: should_reject
1430  REAL(kind=dp) :: distsq, mass, rtmp, rtmpo, sdiff, snew, &
1431  sold
1432  REAL(kind=dp), DIMENSION(3) :: distvec, dr, dro, new_com, old_com
1433  REAL(kind=dp), DIMENSION(3, l-1) :: newsection
1434 
1435  mass = helium%he_mass_au
1436 
1437  endatom = helium%worm_atom_idx
1438  endbead = helium%worm_bead_idx
1439  ! get index of the atom and bead, where the resampling of the head begins
1440  IF (l < endbead) THEN
1441  ! startbead belongs to the same atom
1442  startatom = endatom
1443  startbead = endbead - l
1444  ELSE
1445  ! startbead belongs to a different atom
1446  ! find previous atom (assuming l < nbeads)
1447  startatom = helium%iperm(endatom)
1448  startbead = endbead + helium%beads - l
1449  END IF
1450  sold = worm_path_action_worm_corrected(helium, helium%pos, &
1451  startatom, startbead, endatom, endbead, &
1452  helium%worm_xtra_bead, helium%worm_atom_idx, helium%worm_bead_idx)
1453 
1454  IF (helium%solute_present) THEN
1455  sold = sold + worm_path_inter_action_head(pint_env, helium, helium%pos, &
1456  startatom, startbead, &
1457  helium%worm_xtra_bead, helium%worm_atom_idx, helium%worm_bead_idx)
1458  END IF
1459 
1460  ! close between head and tail
1461  ! only l-1 beads need to be reconstructed
1462  CALL path_construct(helium, &
1463  helium%pos(:, startatom, startbead), &
1464  helium%pos(:, endatom, endbead), l - 1, &
1465  newsection)
1466 
1467  ! transfer the new coordinates to work array
1468  jbead = 1
1469  IF (startbead < endbead) THEN
1470  ! everything belongs to the same atom
1471  DO kbead = startbead + 1, endbead - 1
1472  helium%work(:, endatom, kbead) = newsection(:, jbead)
1473  jbead = jbead + 1
1474  END DO
1475  ELSE
1476  ! is distributed among two atoms
1477  ! transfer to atom not containing the head bead
1478  DO kbead = startbead + 1, helium%beads
1479  helium%work(:, startatom, kbead) = newsection(:, jbead)
1480  jbead = jbead + 1
1481  END DO
1482  ! transfer to atom containing the head bead
1483  DO kbead = 1, endbead - 1
1484  helium%work(:, endatom, kbead) = newsection(:, jbead)
1485  jbead = jbead + 1
1486  END DO
1487  END IF
1488 
1489  helium%worm_is_closed = .true.
1490 
1491  snew = worm_path_action(helium, helium%work, &
1492  startatom, startbead, endatom, endbead)
1493 
1494  IF (helium%solute_present) THEN
1495  ! yes this is correct, as the bead, that was split into tail and head only changes half
1496  ! therefore only half of its action needs to be considred
1497  ! and is cheated in here by passing it as head bead
1498  snew = snew + worm_path_inter_action_head(pint_env, helium, helium%work, &
1499  startatom, startbead, &
1500  helium%work(:, endatom, endbead), endatom, endbead)
1501  END IF
1502 
1503  ! Metropolis:
1504  ! first compute ln of free density matrix
1505  distvec(:) = helium%pos(:, startatom, startbead) - helium%pos(:, endatom, endbead)
1506  CALL helium_pbc(helium, distvec)
1507  distsq = dot_product(distvec, distvec)
1508  ! action difference
1509  sdiff = sold - snew
1510  ! modify action difference due to extra bead
1511  sdiff = sdiff - distsq/(2.0_dp*helium%hb2m*real(l, dp)*helium%tau)
1512  sdiff = sdiff - 1.5_dp*log(real(l, dp)*helium%tau)
1513  sdiff = sdiff - helium%worm_ln_openclose_scale
1514  IF (sdiff < 0) THEN
1515  should_reject = .false.
1516  IF (sdiff < -100.0_dp) THEN ! To protect from exponential underflow
1517  should_reject = .true.
1518  ELSE
1519  rtmp = helium%rng_stream_uniform%next()
1520  IF (exp(sdiff) < rtmp) THEN
1521  should_reject = .true.
1522  END IF
1523  END IF
1524  IF (should_reject) THEN
1525  ! rejected !
1526  ! write back only changed atoms
1527  ! transfer the new coordinates to work array
1528  IF (startbead < endbead) THEN
1529  ! everything belongs to the same atom
1530  DO kbead = startbead + 1, endbead - 1
1531  helium%work(:, endatom, kbead) = helium%pos(:, endatom, kbead)
1532  END DO
1533  ELSE
1534  ! is distributed among two atoms
1535  ! transfer to atom not containing the head bead
1536  DO kbead = startbead + 1, helium%beads
1537  helium%work(:, startatom, kbead) = helium%pos(:, startatom, kbead)
1538  END DO
1539  ! transfer to atom containing the head bead
1540  DO kbead = 1, endbead - 1
1541  helium%work(:, endatom, kbead) = helium%pos(:, endatom, kbead)
1542  END DO
1543  END IF
1544  helium%worm_is_closed = .false.
1545  ac = 0
1546  RETURN
1547  END IF
1548  END IF
1549 
1550  ! for now accepted
1551  ! rejection of the whole move if any centroid is farther away
1552  ! from the current center of gravity than HELIUM%DROPLET_RADIUS [lwalewski]
1553  ! AND ist not moving towards the center
1554  IF (.NOT. helium%periodic) THEN
1555  IF (helium%solute_present) THEN
1556  new_com(:) = helium%center(:)
1557  old_com(:) = new_com(:)
1558  ELSE
1559  new_com(:) = 0.0_dp
1560  DO ia = 1, helium%atoms
1561  DO ib = 1, helium%beads
1562  new_com(:) = new_com(:) + helium%work(:, ia, ib)
1563  END DO
1564  END DO
1565  new_com(:) = new_com(:)/(real(helium%atoms*helium%beads, dp))
1566  ! also compute the old center of mass (optimization potential)
1567  old_com(:) = 0.0_dp
1568  DO ia = 1, helium%atoms
1569  DO ib = 1, helium%beads
1570  old_com(:) = old_com(:) + helium%pos(:, ia, ib)
1571  END DO
1572  END DO
1573  old_com(:) = old_com(:)/(real(helium%atoms*helium%beads, dp))
1574  END IF
1575  should_reject = .false.
1576  atom: DO ia = 1, helium%atoms
1577  dr(:) = 0.0_dp
1578  DO ib = 1, helium%beads
1579  dr(:) = dr(:) + helium%work(:, ia, ib) - new_com(:)
1580  END DO
1581  dr(:) = dr(:)/real(helium%beads, dp)
1582  rtmp = dot_product(dr, dr)
1583  IF (rtmp >= helium%droplet_radius**2) THEN
1584  dro(:) = 0.0_dp
1585  DO ib = 1, helium%beads
1586  dro(:) = dro(:) + helium%pos(:, ia, ib) - old_com(:)
1587  END DO
1588  dro(:) = dro(:)/real(helium%beads, dp)
1589  rtmpo = dot_product(dro, dro)
1590  ! only reject if it moves away from COM outside the droplet radius
1591  IF (rtmpo < rtmp) THEN
1592  ! found - this one does not fit within R from the center
1593  should_reject = .true.
1594  EXIT atom
1595  END IF
1596  END IF
1597  END DO atom
1598  IF (should_reject) THEN
1599  ! restore original coordinates
1600  ! write back only changed atoms
1601  IF (startbead < endbead) THEN
1602  ! everything belongs to the same atom
1603  DO kbead = startbead + 1, endbead - 1
1604  helium%work(:, endatom, kbead) = helium%pos(:, endatom, kbead)
1605  END DO
1606  ELSE
1607  ! is distributed among two atoms
1608  ! transfer to atom not containing the head bead
1609  DO kbead = startbead + 1, helium%beads
1610  helium%work(:, startatom, kbead) = helium%pos(:, startatom, kbead)
1611  END DO
1612  ! transfer to atom containing the head bead
1613  DO kbead = 1, endbead - 1
1614  helium%work(:, endatom, kbead) = helium%pos(:, endatom, kbead)
1615  END DO
1616  END IF
1617  helium%worm_is_closed = .false.
1618  ac = 0
1619  RETURN
1620  END IF
1621  END IF
1622 
1623  ! finally accepted
1624  ac = 1
1625  ! write changed coordinates to position array
1626  IF (startbead < endbead) THEN
1627  ! everything belongs to the same atom
1628  DO kbead = startbead + 1, endbead - 1
1629  helium%pos(:, endatom, kbead) = helium%work(:, endatom, kbead)
1630  END DO
1631  ELSE
1632  ! is distributed among two atoms
1633  ! transfer to atom not containing the head bead
1634  DO kbead = startbead + 1, helium%beads
1635  helium%pos(:, startatom, kbead) = helium%work(:, startatom, kbead)
1636  END DO
1637  ! transfer to atom containing the head bead
1638  DO kbead = 1, endbead - 1
1639  helium%pos(:, endatom, kbead) = helium%work(:, endatom, kbead)
1640  END DO
1641  END IF
1642 
1643  END SUBROUTINE worm_close_move
1644 
1645 ! **************************************************************************************************
1646 !> \brief Move head in open state
1647 !> \param pint_env ...
1648 !> \param helium ...
1649 !> \param l ...
1650 !> \param ac ...
1651 !> \author fuhl
1652 ! **************************************************************************************************
1653  SUBROUTINE worm_head_move(pint_env, helium, l, ac)
1654 
1655  TYPE(pint_env_type), INTENT(IN) :: pint_env
1656  TYPE(helium_solvent_type), INTENT(INOUT) :: helium
1657  INTEGER, INTENT(IN) :: l
1658  INTEGER, INTENT(OUT) :: ac
1659 
1660  INTEGER :: endatom, endbead, ia, ib, idim, kbead, &
1661  startatom, startbead
1662  LOGICAL :: should_reject
1663  REAL(kind=dp) :: mass, rtmp, rtmpo, sdiff, snew, sold, xr
1664  REAL(kind=dp), DIMENSION(3) :: dr, dro, new_com, old_com
1665 
1666  mass = helium%he_mass_au
1667 
1668  ! get index of the atom and bead, where the resampling of the head begins
1669  endatom = helium%worm_atom_idx
1670  endbead = helium%worm_bead_idx
1671  IF (l < endbead) THEN
1672  ! startbead belongs to the same atom
1673  startatom = endatom
1674  startbead = endbead - l
1675  ELSE
1676  ! startbead belongs to a different atom
1677  ! find previous atom (assuming l < nbeads)
1678  startatom = helium%iperm(endatom)
1679  startbead = endbead + helium%beads - l
1680  END IF
1681 
1682  sold = worm_path_action_worm_corrected(helium, helium%pos, &
1683  startatom, startbead, endatom, endbead, &
1684  helium%worm_xtra_bead, helium%worm_atom_idx, helium%worm_bead_idx)
1685 
1686  IF (helium%solute_present) THEN
1687  sold = sold + worm_path_inter_action_head(pint_env, helium, helium%pos, &
1688  startatom, startbead, &
1689  helium%worm_xtra_bead, helium%worm_atom_idx, helium%worm_bead_idx)
1690  END IF
1691 
1692  ! alternative grow with consecutive gaussians
1693  IF (startbead < endbead) THEN
1694  ! everything belongs to the same atom
1695  ! gro head from startbead
1696  DO kbead = startbead + 1, endbead - 1
1697  DO idim = 1, 3
1698  xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
1699  helium%work(idim, startatom, kbead) = helium%work(idim, startatom, kbead - 1) + xr
1700  END DO
1701  END DO
1702  ! last grow head bead
1703  DO idim = 1, 3
1704  xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
1705  helium%worm_xtra_bead_work(idim) = helium%work(idim, startatom, endbead - 1) + xr
1706  END DO
1707  ELSE IF (endbead /= 1) THEN
1708  ! is distributed among two atoms
1709  ! grow from startbead
1710  DO kbead = startbead + 1, helium%beads
1711  DO idim = 1, 3
1712  xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
1713  helium%work(idim, startatom, kbead) = helium%work(idim, startatom, kbead - 1) + xr
1714  END DO
1715  END DO
1716  ! bead one of endatom relative to last on startatom
1717  DO idim = 1, 3
1718  xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
1719  helium%work(idim, endatom, 1) = helium%work(idim, startatom, helium%beads) + xr
1720  END DO
1721  ! everything on endatom
1722  DO kbead = 2, endbead - 1
1723  DO idim = 1, 3
1724  xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
1725  helium%work(idim, endatom, kbead) = helium%work(idim, endatom, kbead - 1) + xr
1726  END DO
1727  END DO
1728  DO idim = 1, 3
1729  xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
1730  helium%worm_xtra_bead_work(idim) = helium%work(idim, endatom, endbead - 1) + xr
1731  END DO
1732  ELSE ! imagtimewrap and headbead = 1
1733  ! is distributed among two atoms
1734  ! grow from startbead
1735  DO kbead = startbead + 1, helium%beads
1736  DO idim = 1, 3
1737  xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
1738  helium%work(idim, startatom, kbead) = helium%work(idim, startatom, kbead - 1) + xr
1739  END DO
1740  END DO
1741  ! bead one of endatom relative to last on startatom
1742  DO idim = 1, 3
1743  xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
1744  helium%worm_xtra_bead_work(idim) = helium%work(idim, startatom, helium%beads) + xr
1745  END DO
1746  END IF
1747 
1748  snew = worm_path_action_worm_corrected(helium, helium%work, &
1749  startatom, startbead, endatom, endbead, &
1750  helium%worm_xtra_bead_work, helium%worm_atom_idx_work, helium%worm_bead_idx_work)
1751 
1752  IF (helium%solute_present) THEN
1753  snew = snew + worm_path_inter_action_head(pint_env, helium, helium%work, &
1754  startatom, startbead, &
1755  helium%worm_xtra_bead_work, helium%worm_atom_idx_work, helium%worm_bead_idx_work)
1756  END IF
1757 
1758  ! Metropolis:
1759  ! action difference
1760  sdiff = sold - snew
1761  ! modify action difference due to extra bead
1762  IF (sdiff < 0) THEN
1763  should_reject = .false.
1764  IF (sdiff < -100.0_dp) THEN ! To protect from exponential underflow
1765  should_reject = .true.
1766  ELSE
1767  rtmp = helium%rng_stream_uniform%next()
1768  IF (exp(sdiff) < rtmp) THEN
1769  should_reject = .true.
1770  END IF
1771  END IF
1772  IF (should_reject) THEN
1773  ! rejected !
1774  ! write back only changed atoms
1775  ! transfer the new coordinates to work array
1776  IF (startbead < endbead) THEN
1777  ! everything belongs to the same atom
1778  DO kbead = startbead + 1, endbead - 1
1779  helium%work(:, endatom, kbead) = helium%pos(:, endatom, kbead)
1780  END DO
1781  ELSE
1782  ! is distributed among two atoms
1783  ! transfer to atom not containing the head bead
1784  DO kbead = startbead + 1, helium%beads
1785  helium%work(:, startatom, kbead) = helium%pos(:, startatom, kbead)
1786  END DO
1787  ! transfer to atom containing the head bead
1788  DO kbead = 1, endbead - 1
1789  helium%work(:, endatom, kbead) = helium%pos(:, endatom, kbead)
1790  END DO
1791  END IF
1792  helium%worm_xtra_bead_work(:) = helium%worm_xtra_bead(:)
1793  ac = 0
1794  RETURN
1795  END IF
1796  END IF
1797 
1798  ! for now accepted
1799  ! rejection of the whole move if any centroid is farther away
1800  ! from the current center of gravity than HELIUM%DROPLET_RADIUS [lwalewski]
1801  ! AND ist not moving towards the center
1802  IF (.NOT. helium%periodic) THEN
1803  IF (helium%solute_present) THEN
1804  new_com(:) = helium%center(:)
1805  old_com(:) = new_com(:)
1806  ELSE
1807  new_com(:) = 0.0_dp
1808  DO ia = 1, helium%atoms
1809  DO ib = 1, helium%beads
1810  new_com(:) = new_com(:) + helium%work(:, ia, ib)
1811  END DO
1812  END DO
1813  new_com(:) = new_com(:)/(real(helium%atoms*helium%beads, dp))
1814  ! also compute the old center of mass (optimization potential)
1815  old_com(:) = 0.0_dp
1816  DO ia = 1, helium%atoms
1817  DO ib = 1, helium%beads
1818  old_com(:) = old_com(:) + helium%pos(:, ia, ib)
1819  END DO
1820  END DO
1821  old_com(:) = old_com(:)/(real(helium%atoms*helium%beads, dp))
1822  END IF
1823  should_reject = .false.
1824  atom: DO ia = 1, helium%atoms
1825  dr(:) = 0.0_dp
1826  DO ib = 1, helium%beads
1827  dr(:) = dr(:) + helium%work(:, ia, ib) - new_com(:)
1828  END DO
1829  dr(:) = dr(:)/real(helium%beads, dp)
1830  rtmp = dot_product(dr, dr)
1831  IF (rtmp >= helium%droplet_radius**2) THEN
1832  dro(:) = 0.0_dp
1833  DO ib = 1, helium%beads
1834  dro(:) = dro(:) + helium%pos(:, ia, ib) - old_com(:)
1835  END DO
1836  dro(:) = dro(:)/real(helium%beads, dp)
1837  rtmpo = dot_product(dro, dro)
1838  ! only reject if it moves away from COM outside the droplet radius
1839  IF (rtmpo < rtmp) THEN
1840  ! found - this one does not fit within R from the center
1841  should_reject = .true.
1842  EXIT atom
1843  END IF
1844  END IF
1845  END DO atom
1846  IF (should_reject) THEN
1847  ! restore original coordinates
1848  ! write back only changed atoms
1849  IF (startbead < endbead) THEN
1850  ! everything belongs to the same atom
1851  DO kbead = startbead + 1, endbead - 1
1852  helium%work(:, endatom, kbead) = helium%pos(:, endatom, kbead)
1853  END DO
1854  ELSE
1855  ! is distributed among two atoms
1856  ! transfer to atom not containing the head bead
1857  DO kbead = startbead + 1, helium%beads
1858  helium%work(:, startatom, kbead) = helium%pos(:, startatom, kbead)
1859  END DO
1860  ! transfer to atom containing the head bead
1861  DO kbead = 1, endbead - 1
1862  helium%work(:, endatom, kbead) = helium%pos(:, endatom, kbead)
1863  END DO
1864  END IF
1865  helium%worm_xtra_bead_work(:) = helium%worm_xtra_bead(:)
1866  ac = 0
1867  RETURN
1868  END IF
1869  END IF
1870 
1871  ! finally accepted
1872  ac = 1
1873  ! write changed coordinates to position array
1874  IF (startbead < endbead) THEN
1875  ! everything belongs to the same atom
1876  DO kbead = startbead + 1, endbead - 1
1877  helium%pos(:, endatom, kbead) = helium%work(:, endatom, kbead)
1878  END DO
1879  ELSE
1880  ! is distributed among two atoms
1881  ! transfer to atom not containing the head bead
1882  DO kbead = startbead + 1, helium%beads
1883  helium%pos(:, startatom, kbead) = helium%work(:, startatom, kbead)
1884  END DO
1885  ! transfer to atom containing the head bead
1886  DO kbead = 1, endbead - 1
1887  helium%pos(:, endatom, kbead) = helium%work(:, endatom, kbead)
1888  END DO
1889  END IF
1890  helium%worm_xtra_bead(:) = helium%worm_xtra_bead_work(:)
1891 
1892  END SUBROUTINE worm_head_move
1893 
1894 ! **************************************************************************************************
1895 !> \brief Move tail in open state
1896 !> \param pint_env ...
1897 !> \param helium ...
1898 !> \param l ...
1899 !> \param ac ...
1900 !> \author fuhl
1901 ! **************************************************************************************************
1902  SUBROUTINE worm_tail_move(pint_env, helium, l, ac)
1903 
1904  TYPE(pint_env_type), INTENT(IN) :: pint_env
1905  TYPE(helium_solvent_type), INTENT(INOUT) :: helium
1906  INTEGER, INTENT(IN) :: l
1907  INTEGER, INTENT(OUT) :: ac
1908 
1909  INTEGER :: endatom, endbead, ia, ib, idim, kbead, &
1910  startatom, startbead
1911  LOGICAL :: should_reject
1912  REAL(kind=dp) :: mass, rtmp, rtmpo, sdiff, snew, sold, xr
1913  REAL(kind=dp), DIMENSION(3) :: dr, dro, new_com, old_com
1914 
1915  mass = helium%he_mass_au
1916 
1917  ! get index of the atom and bead, where the resampling of the tail ends
1918  startatom = helium%worm_atom_idx
1919  startbead = helium%worm_bead_idx
1920  endbead = startbead + l
1921 
1922  IF (endbead <= helium%beads) THEN
1923  ! endbead belongs to the same atom
1924  endatom = startatom
1925  ELSE
1926  ! endbead belongs to a different atom
1927  ! find next atom (assuming l < nbeads)
1928  endatom = helium%permutation(startatom)
1929  endbead = endbead - helium%beads
1930  END IF
1931 
1932  !yes this is correct, as the head does not play any role here
1933  sold = worm_path_action(helium, helium%pos, &
1934  startatom, startbead, endatom, endbead)
1935 
1936  IF (helium%solute_present) THEN
1937  sold = sold + worm_path_inter_action_tail(pint_env, helium, helium%pos, &
1938  endatom, endbead, &
1939  helium%worm_atom_idx, helium%worm_bead_idx)
1940  END IF
1941 
1942  ! alternative grow with consecutive gaussians
1943  IF (startbead < endbead) THEN
1944  ! everything belongs to the same atom
1945  ! gro tail from endbead to startbead (confusing eh?)
1946  DO kbead = endbead - 1, startbead, -1
1947  DO idim = 1, 3
1948  xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
1949  helium%work(idim, startatom, kbead) = helium%work(idim, startatom, kbead + 1) + xr
1950  END DO
1951  END DO
1952  ELSE
1953  ! is distributed among two atoms
1954  ! grow from endbead
1955  DO kbead = endbead - 1, 1, -1
1956  DO idim = 1, 3
1957  xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
1958  helium%work(idim, endatom, kbead) = helium%work(idim, endatom, kbead + 1) + xr
1959  END DO
1960  END DO
1961 
1962  ! over imaginary time boundary
1963  DO idim = 1, 3
1964  xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
1965  helium%work(idim, startatom, helium%beads) = helium%work(idim, endatom, 1) + xr
1966  END DO
1967 
1968  ! rest on startatom
1969  DO kbead = helium%beads - 1, startbead, -1
1970  DO idim = 1, 3
1971  xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
1972  helium%work(idim, startatom, kbead) = helium%work(idim, startatom, kbead + 1) + xr
1973  END DO
1974  END DO
1975  END IF
1976 
1977  !yes this is correct, as the head does not play any role here
1978  snew = worm_path_action(helium, helium%work, &
1979  startatom, startbead, endatom, endbead)
1980 
1981  IF (helium%solute_present) THEN
1982  snew = snew + worm_path_inter_action_tail(pint_env, helium, helium%work, &
1983  endatom, endbead, &
1984  helium%worm_atom_idx_work, helium%worm_bead_idx_work)
1985  END IF
1986 
1987  ! Metropolis:
1988  ! action difference
1989  sdiff = sold - snew
1990  ! modify action difference due to extra bead
1991  IF (sdiff < 0) THEN
1992  should_reject = .false.
1993  IF (sdiff < -100.0_dp) THEN ! To protect from exponential underflow
1994  should_reject = .true.
1995  ELSE
1996  rtmp = helium%rng_stream_uniform%next()
1997  IF (exp(sdiff) < rtmp) THEN
1998  should_reject = .true.
1999  END IF
2000  END IF
2001  IF (should_reject) THEN
2002  ! rejected !
2003  ! write back only changed atoms
2004  ! transfer the new coordinates to work array
2005  IF (startbead < endbead) THEN
2006  ! everything belongs to the same atom
2007  DO kbead = startbead, endbead - 1
2008  helium%work(:, endatom, kbead) = helium%pos(:, endatom, kbead)
2009  END DO
2010  ELSE
2011  ! is distributed among two atoms
2012  ! transfer to atom not containing the tail bead
2013  DO kbead = startbead, helium%beads
2014  helium%work(:, startatom, kbead) = helium%pos(:, startatom, kbead)
2015  END DO
2016  ! transfer to atom containing the tail bead
2017  DO kbead = 1, endbead - 1
2018  helium%work(:, endatom, kbead) = helium%pos(:, endatom, kbead)
2019  END DO
2020  END IF
2021  helium%worm_xtra_bead_work(:) = helium%worm_xtra_bead(:)
2022  ac = 0
2023  RETURN
2024  END IF
2025  END IF
2026 
2027  ! for now accepted
2028  ! rejection of the whole move if any centroid is farther away
2029  ! from the current center of gravity than HELIUM%DROPLET_RADIUS [lwalewski]
2030  ! AND ist not moving towards the center
2031  IF (.NOT. helium%periodic) THEN
2032  IF (helium%solute_present) THEN
2033  new_com(:) = helium%center(:)
2034  old_com(:) = new_com(:)
2035  ELSE
2036  new_com(:) = 0.0_dp
2037  DO ia = 1, helium%atoms
2038  DO ib = 1, helium%beads
2039  new_com(:) = new_com(:) + helium%work(:, ia, ib)
2040  END DO
2041  END DO
2042  new_com(:) = new_com(:)/(real(helium%atoms*helium%beads, dp))
2043  ! also compute the old center of mass (optimization potential)
2044  old_com(:) = 0.0_dp
2045  DO ia = 1, helium%atoms
2046  DO ib = 1, helium%beads
2047  old_com(:) = old_com(:) + helium%pos(:, ia, ib)
2048  END DO
2049  END DO
2050  old_com(:) = old_com(:)/(real(helium%atoms*helium%beads, dp))
2051  END IF
2052  should_reject = .false.
2053  atom: DO ia = 1, helium%atoms
2054  dr(:) = 0.0_dp
2055  DO ib = 1, helium%beads
2056  dr(:) = dr(:) + helium%work(:, ia, ib) - new_com(:)
2057  END DO
2058  dr(:) = dr(:)/real(helium%beads, dp)
2059  rtmp = dot_product(dr, dr)
2060  IF (rtmp >= helium%droplet_radius**2) THEN
2061  dro(:) = 0.0_dp
2062  DO ib = 1, helium%beads
2063  dro(:) = dro(:) + helium%pos(:, ia, ib) - old_com(:)
2064  END DO
2065  dro(:) = dro(:)/real(helium%beads, dp)
2066  rtmpo = dot_product(dro, dro)
2067  ! only reject if it moves away from COM outside the droplet radius
2068  IF (rtmpo < rtmp) THEN
2069  ! found - this one does not fit within R from the center
2070  should_reject = .true.
2071  EXIT atom
2072  END IF
2073  END IF
2074  END DO atom
2075  IF (should_reject) THEN
2076  ! restore original coordinates
2077  ! write back only changed atoms
2078  IF (startbead < endbead) THEN
2079  ! everything belongs to the same atom
2080  DO kbead = startbead, endbead - 1
2081  helium%work(:, endatom, kbead) = helium%pos(:, endatom, kbead)
2082  END DO
2083  ELSE
2084  ! is distributed among two atoms
2085  ! transfer to atom not containing the tail bead
2086  DO kbead = startbead, helium%beads
2087  helium%work(:, startatom, kbead) = helium%pos(:, startatom, kbead)
2088  END DO
2089  ! transfer to atom containing the tail bead
2090  DO kbead = 1, endbead - 1
2091  helium%work(:, endatom, kbead) = helium%pos(:, endatom, kbead)
2092  END DO
2093  END IF
2094  helium%worm_xtra_bead_work(:) = helium%worm_xtra_bead(:)
2095  ac = 0
2096  RETURN
2097  END IF
2098  END IF
2099 
2100  ! finally accepted
2101  ac = 1
2102  ! write changed coordinates to position array
2103  IF (startbead < endbead) THEN
2104  ! everything belongs to the same atom
2105  DO kbead = startbead, endbead - 1
2106  helium%pos(:, endatom, kbead) = helium%work(:, endatom, kbead)
2107  END DO
2108  ELSE
2109  ! is distributed among two atoms
2110  ! transfer to atom not containing the tail bead
2111  DO kbead = startbead, helium%beads
2112  helium%pos(:, startatom, kbead) = helium%work(:, startatom, kbead)
2113  END DO
2114  ! transfer to atom containing the tail bead
2115  DO kbead = 1, endbead - 1
2116  helium%pos(:, endatom, kbead) = helium%work(:, endatom, kbead)
2117  END DO
2118  END IF
2119 
2120  END SUBROUTINE worm_tail_move
2121 
2122 ! **************************************************************************************************
2123 !> \brief Crawl forward in open state, can avoid being trapped in open state
2124 !> \param pint_env ...
2125 !> \param helium ...
2126 !> \param l ...
2127 !> \param ac ...
2128 !> \author fuhl
2129 ! **************************************************************************************************
2130  SUBROUTINE worm_crawl_move_forward(pint_env, helium, l, ac)
2131 
2132  TYPE(pint_env_type), INTENT(IN) :: pint_env
2133  TYPE(helium_solvent_type), INTENT(INOUT) :: helium
2134  INTEGER, INTENT(IN) :: l
2135  INTEGER, INTENT(OUT) :: ac
2136 
2137  INTEGER :: ia, ib, idim, kbead
2138  LOGICAL :: should_reject
2139  REAL(kind=dp) :: mass, newtailpot, oldheadpot, &
2140  oldtailpot, rtmp, rtmpo, sdiff, snew, &
2141  sold, xr
2142  REAL(kind=dp), DIMENSION(3) :: dr, dro, new_com, old_com
2143 
2144  mass = helium%he_mass_au
2145 
2146  ! determine position of new head in imaginary time
2147  helium%worm_bead_idx_work = helium%worm_bead_idx + l
2148  IF (helium%worm_bead_idx_work > helium%beads) THEN
2149  helium%worm_bead_idx_work = helium%worm_bead_idx_work - helium%beads
2150  helium%worm_atom_idx_work = helium%permutation(helium%worm_atom_idx)
2151  ELSE
2152  helium%worm_atom_idx_work = helium%worm_atom_idx
2153  END IF
2154 
2155  ! compute action before move
2156  ! head is not involved in pair action
2157  sold = worm_path_action(helium, helium%pos, &
2158  helium%worm_atom_idx, helium%worm_bead_idx, &
2159  helium%worm_atom_idx_work, helium%worm_bead_idx_work)
2160  IF (helium%solute_present) THEN
2161  !this will leave out the old and new tail bead
2162  ! due to efficiency reasons they are treated separately
2163  sold = sold + worm_path_inter_action(pint_env, helium, helium%pos, &
2164  helium%worm_atom_idx, helium%worm_bead_idx, &
2165  helium%worm_atom_idx_work, helium%worm_bead_idx_work)
2166 
2167  ! compute old/new head/tail interactions
2168  ! old tail
2169  CALL helium_bead_solute_e_f(pint_env, helium, &
2170  helium%worm_atom_idx, helium%worm_bead_idx, &
2171  helium%pos(:, helium%worm_atom_idx, helium%worm_bead_idx), &
2172  energy=oldtailpot)
2173  oldtailpot = oldtailpot*helium%tau
2174 
2175  ! new tail
2176  CALL helium_bead_solute_e_f(pint_env, helium, &
2177  helium%worm_atom_idx_work, helium%worm_bead_idx_work, &
2178  helium%pos(:, helium%worm_atom_idx_work, helium%worm_bead_idx_work), &
2179  energy=newtailpot)
2180  newtailpot = newtailpot*helium%tau
2181 
2182  ! old head
2183  CALL helium_bead_solute_e_f(pint_env, helium, &
2184  helium%worm_atom_idx, helium%worm_bead_idx, &
2185  helium%worm_xtra_bead, &
2186  energy=oldheadpot)
2187  oldheadpot = oldheadpot*helium%tau
2188 
2189  ! new head is not known yet
2190 
2191  sold = sold + 0.5_dp*(oldtailpot + oldheadpot) + newtailpot
2192  END IF
2193 
2194  ! copy over old head position to working array and grow from there
2195  helium%work(:, helium%worm_atom_idx, helium%worm_bead_idx) = helium%worm_xtra_bead
2196 
2197  ! grow head part with consecutive gaussians
2198  IF (helium%worm_bead_idx < helium%worm_bead_idx_work) THEN
2199  ! everything belongs to the same atom
2200  ! gro head from startbead
2201  DO kbead = helium%worm_bead_idx + 1, helium%worm_bead_idx_work - 1
2202  DO idim = 1, 3
2203  xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
2204  helium%work(idim, helium%worm_atom_idx, kbead) = helium%work(idim, helium%worm_atom_idx, kbead - 1) + xr
2205  END DO
2206  END DO
2207  ! last grow head bead
2208  DO idim = 1, 3
2209  xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
2210  helium%worm_xtra_bead_work(idim) = helium%work(idim, helium%worm_atom_idx, helium%worm_bead_idx_work - 1) + xr
2211  END DO
2212  ELSE IF (helium%worm_bead_idx_work /= 1) THEN
2213  ! is distributed among two atoms
2214  ! grow from startbead
2215  DO kbead = helium%worm_bead_idx + 1, helium%beads
2216  DO idim = 1, 3
2217  xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
2218  helium%work(idim, helium%worm_atom_idx, kbead) = helium%work(idim, helium%worm_atom_idx, kbead - 1) + xr
2219  END DO
2220  END DO
2221  ! bead one of endatom relative to last on helium%worm_atom_idx
2222  DO idim = 1, 3
2223  xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
2224  helium%work(idim, helium%worm_atom_idx_work, 1) = helium%work(idim, helium%worm_atom_idx, helium%beads) + xr
2225  END DO
2226  ! everything on endatom
2227  DO kbead = 2, helium%worm_bead_idx_work - 1
2228  DO idim = 1, 3
2229  xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
2230  helium%work(idim, helium%worm_atom_idx_work, kbead) = helium%work(idim, helium%worm_atom_idx_work, kbead - 1) + xr
2231  END DO
2232  END DO
2233  DO idim = 1, 3
2234  xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
2235  helium%worm_xtra_bead_work(idim) = helium%work(idim, helium%worm_atom_idx_work, helium%worm_bead_idx_work - 1) + xr
2236  END DO
2237  ELSE ! imagtimewrap and headbead = 1
2238  ! is distributed among two atoms
2239  ! grow from startbead
2240  DO kbead = helium%worm_bead_idx + 1, helium%beads
2241  DO idim = 1, 3
2242  xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
2243  helium%work(idim, helium%worm_atom_idx, kbead) = helium%work(idim, helium%worm_atom_idx, kbead - 1) + xr
2244  END DO
2245  END DO
2246  ! bead one of endatom relative to last on helium%worm_atom_idx
2247  DO idim = 1, 3
2248  xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
2249  helium%worm_xtra_bead_work(idim) = helium%work(idim, helium%worm_atom_idx, helium%beads) + xr
2250  END DO
2251  END IF
2252 
2253  snew = worm_path_action_worm_corrected(helium, helium%work, &
2254  helium%worm_atom_idx, helium%worm_bead_idx, &
2255  helium%worm_atom_idx_work, helium%worm_bead_idx_work, &
2256  helium%worm_xtra_bead_work, helium%worm_atom_idx_work, helium%worm_bead_idx_work)
2257 
2258  IF (helium%solute_present) THEN
2259  snew = snew + worm_path_inter_action_head(pint_env, helium, helium%work, &
2260  helium%worm_atom_idx, helium%worm_bead_idx, &
2261  helium%worm_xtra_bead_work, helium%worm_atom_idx_work, helium%worm_bead_idx_work)
2262  snew = snew + 0.5_dp*newtailpot + oldheadpot
2263  END IF
2264 
2265  ! Metropolis:
2266  ! action difference
2267  sdiff = sold - snew
2268  ! modify action difference due to extra bead
2269  IF (sdiff < 0) THEN
2270  should_reject = .false.
2271  IF (sdiff < -100.0_dp) THEN ! To protect from exponential underflow
2272  should_reject = .true.
2273  ELSE
2274  rtmp = helium%rng_stream_uniform%next()
2275  IF (exp(sdiff) < rtmp) THEN
2276  should_reject = .true.
2277  END IF
2278  END IF
2279  IF (should_reject) THEN
2280  ! rejected !
2281  ! write back only changed atoms
2282  ! transfer the new coordinates to work array
2283  IF (helium%worm_bead_idx < helium%worm_bead_idx_work) THEN
2284  ! everything belongs to the same atom
2285  DO kbead = helium%worm_bead_idx, helium%worm_bead_idx_work - 1
2286  helium%work(:, helium%worm_atom_idx_work, kbead) = helium%pos(:, helium%worm_atom_idx_work, kbead)
2287  END DO
2288  ELSE
2289  ! is distributed among two atoms
2290  ! transfer to atom not containing the head bead
2291  DO kbead = helium%worm_bead_idx, helium%beads
2292  helium%work(:, helium%worm_atom_idx, kbead) = helium%pos(:, helium%worm_atom_idx, kbead)
2293  END DO
2294  ! transfer to atom containing the head bead
2295  DO kbead = 1, helium%worm_bead_idx_work - 1
2296  helium%work(:, helium%worm_atom_idx_work, kbead) = helium%pos(:, helium%worm_atom_idx_work, kbead)
2297  END DO
2298  END IF
2299  helium%worm_xtra_bead_work(:) = helium%worm_xtra_bead(:)
2300  helium%worm_bead_idx_work = helium%worm_bead_idx
2301  helium%worm_atom_idx_work = helium%worm_atom_idx
2302  ac = 0
2303  RETURN
2304  END IF
2305  END IF
2306 
2307  ! for now accepted
2308  ! rejection of the whole move if any centroid is farther away
2309  ! from the current center of gravity than HELIUM%DROPLET_RADIUS [lwalewski]
2310  ! AND ist not moving towards the center
2311  IF (.NOT. helium%periodic) THEN
2312  IF (helium%solute_present) THEN
2313  new_com(:) = helium%center(:)
2314  old_com(:) = new_com(:)
2315  ELSE
2316  new_com(:) = 0.0_dp
2317  DO ia = 1, helium%atoms
2318  DO ib = 1, helium%beads
2319  new_com(:) = new_com(:) + helium%work(:, ia, ib)
2320  END DO
2321  END DO
2322  new_com(:) = new_com(:)/(real(helium%atoms*helium%beads, dp))
2323  ! also compute the old center of mass (optimization potential)
2324  old_com(:) = 0.0_dp
2325  DO ia = 1, helium%atoms
2326  DO ib = 1, helium%beads
2327  old_com(:) = old_com(:) + helium%pos(:, ia, ib)
2328  END DO
2329  END DO
2330  old_com(:) = old_com(:)/(real(helium%atoms*helium%beads, dp))
2331  END IF
2332  should_reject = .false.
2333  atom: DO ia = 1, helium%atoms
2334  dr(:) = 0.0_dp
2335  DO ib = 1, helium%beads
2336  dr(:) = dr(:) + helium%work(:, ia, ib) - new_com(:)
2337  END DO
2338  dr(:) = dr(:)/real(helium%beads, dp)
2339  rtmp = dot_product(dr, dr)
2340  IF (rtmp >= helium%droplet_radius**2) THEN
2341  dro(:) = 0.0_dp
2342  DO ib = 1, helium%beads
2343  dro(:) = dro(:) + helium%pos(:, ia, ib) - old_com(:)
2344  END DO
2345  dro(:) = dro(:)/real(helium%beads, dp)
2346  rtmpo = dot_product(dro, dro)
2347  ! only reject if it moves away from COM outside the droplet radius
2348  IF (rtmpo < rtmp) THEN
2349  ! found - this one does not fit within R from the center
2350  should_reject = .true.
2351  EXIT atom
2352  END IF
2353  END IF
2354  END DO atom
2355  IF (should_reject) THEN
2356  ! restore original coordinates
2357  ! write back only changed atoms
2358  IF (helium%worm_bead_idx < helium%worm_bead_idx_work) THEN
2359  ! everything belongs to the same atom
2360  DO kbead = helium%worm_bead_idx, helium%worm_bead_idx_work - 1
2361  helium%work(:, helium%worm_atom_idx_work, kbead) = helium%pos(:, helium%worm_atom_idx_work, kbead)
2362  END DO
2363  ELSE
2364  ! is distributed among two atoms
2365  ! transfer to atom not containing the head bead
2366  DO kbead = helium%worm_bead_idx, helium%beads
2367  helium%work(:, helium%worm_atom_idx, kbead) = helium%pos(:, helium%worm_atom_idx, kbead)
2368  END DO
2369  ! transfer to atom containing the head bead
2370  DO kbead = 1, helium%worm_bead_idx_work - 1
2371  helium%work(:, helium%worm_atom_idx_work, kbead) = helium%pos(:, helium%worm_atom_idx_work, kbead)
2372  END DO
2373  END IF
2374  helium%worm_xtra_bead_work(:) = helium%worm_xtra_bead(:)
2375  helium%worm_bead_idx_work = helium%worm_bead_idx
2376  helium%worm_atom_idx_work = helium%worm_atom_idx
2377  ac = 0
2378  RETURN
2379  END IF
2380  END IF
2381 
2382  ! finally accepted
2383  ac = 1
2384  ! write changed coordinates to position array
2385  IF (helium%worm_bead_idx < helium%worm_bead_idx_work) THEN
2386  ! everything belongs to the same atom
2387  DO kbead = helium%worm_bead_idx, helium%worm_bead_idx_work - 1
2388  helium%pos(:, helium%worm_atom_idx_work, kbead) = helium%work(:, helium%worm_atom_idx_work, kbead)
2389  END DO
2390  ELSE
2391  ! is distributed among two atoms
2392  ! transfer to atom not containing the head bead
2393  DO kbead = helium%worm_bead_idx, helium%beads
2394  helium%pos(:, helium%worm_atom_idx, kbead) = helium%work(:, helium%worm_atom_idx, kbead)
2395  END DO
2396  ! transfer to atom containing the head bead
2397  DO kbead = 1, helium%worm_bead_idx_work - 1
2398  helium%pos(:, helium%worm_atom_idx_work, kbead) = helium%work(:, helium%worm_atom_idx_work, kbead)
2399  END DO
2400  END IF
2401  helium%worm_xtra_bead(:) = helium%worm_xtra_bead_work(:)
2402  helium%worm_bead_idx = helium%worm_bead_idx_work
2403  helium%worm_atom_idx = helium%worm_atom_idx_work
2404 
2405  END SUBROUTINE worm_crawl_move_forward
2406 
2407 ! **************************************************************************************************
2408 !> \brief Crawl backwards in open state, can avoid being trapped in open state
2409 !> \param pint_env ...
2410 !> \param helium ...
2411 !> \param l ...
2412 !> \param ac ...
2413 !> \author fuhl
2414 ! **************************************************************************************************
2415  SUBROUTINE worm_crawl_move_backward(pint_env, helium, l, ac)
2416 
2417  TYPE(pint_env_type), INTENT(IN) :: pint_env
2418  TYPE(helium_solvent_type), INTENT(INOUT) :: helium
2419  INTEGER, INTENT(IN) :: l
2420  INTEGER, INTENT(OUT) :: ac
2421 
2422  INTEGER :: ia, ib, idim, kbead
2423  LOGICAL :: should_reject
2424  REAL(kind=dp) :: mass, newheadpot, oldheadpot, &
2425  oldtailpot, rtmp, rtmpo, sdiff, snew, &
2426  sold, xr
2427  REAL(kind=dp), DIMENSION(3) :: dr, dro, new_com, old_com
2428 
2429  mass = helium%he_mass_au
2430 
2431  ! determine position of new head in imaginary time
2432  helium%worm_bead_idx_work = helium%worm_bead_idx - l
2433  IF (helium%worm_bead_idx_work < 1) THEN
2434  helium%worm_bead_idx_work = helium%worm_bead_idx_work + helium%beads
2435  helium%worm_atom_idx_work = helium%iperm(helium%worm_atom_idx)
2436  ELSE
2437  helium%worm_atom_idx_work = helium%worm_atom_idx
2438  END IF
2439 
2440  ! compute action before move
2441  ! head is not involved in pair action
2442  sold = worm_path_action_worm_corrected(helium, helium%work, &
2443  helium%worm_atom_idx_work, helium%worm_bead_idx_work, &
2444  helium%worm_atom_idx, helium%worm_bead_idx, &
2445  helium%worm_xtra_bead, helium%worm_atom_idx, helium%worm_bead_idx)
2446 
2447  IF (helium%solute_present) THEN
2448  !this will leave out the old and new tail bead
2449  ! due to efficiency reasons they are treated separately
2450  sold = sold + worm_path_inter_action(pint_env, helium, helium%pos, &
2451  helium%worm_atom_idx_work, helium%worm_bead_idx_work, &
2452  helium%worm_atom_idx, helium%worm_bead_idx)
2453 
2454  ! compute old/new head/tail interactions
2455  ! old tail
2456  CALL helium_bead_solute_e_f(pint_env, helium, &
2457  helium%worm_atom_idx, helium%worm_bead_idx, &
2458  helium%pos(:, helium%worm_atom_idx, helium%worm_bead_idx), &
2459  energy=oldtailpot)
2460  oldtailpot = oldtailpot*helium%tau
2461 
2462  ! old head
2463  CALL helium_bead_solute_e_f(pint_env, helium, &
2464  helium%worm_atom_idx, helium%worm_bead_idx, &
2465  helium%worm_xtra_bead, &
2466  energy=oldheadpot)
2467  oldheadpot = oldheadpot*helium%tau
2468 
2469  ! new head
2470  CALL helium_bead_solute_e_f(pint_env, helium, &
2471  helium%worm_atom_idx_work, helium%worm_bead_idx_work, &
2472  helium%pos(:, helium%worm_atom_idx_work, helium%worm_bead_idx_work), &
2473  energy=newheadpot)
2474  newheadpot = newheadpot*helium%tau
2475 
2476  ! new tail not known yet
2477 
2478  sold = sold + 0.5_dp*(oldtailpot + oldheadpot) + newheadpot
2479  END IF
2480 
2481  ! copy position to the head bead
2482  helium%worm_xtra_bead_work = helium%pos(:, helium%worm_atom_idx_work, helium%worm_bead_idx_work)
2483 
2484  ! alternative grow with consecutive gaussians
2485  IF (helium%worm_bead_idx_work < helium%worm_bead_idx) THEN
2486  ! everything belongs to the same atom
2487  ! gro tail from endbead to startbead (confusing eh?)
2488  DO kbead = helium%worm_bead_idx - 1, helium%worm_bead_idx_work, -1
2489  DO idim = 1, 3
2490  xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
2491  helium%work(idim, helium%worm_atom_idx, kbead) = helium%work(idim, helium%worm_atom_idx, kbead + 1) + xr
2492  END DO
2493 
2494  END DO
2495  ELSE
2496  ! is distributed among two atoms
2497  ! grow from endbead
2498  DO kbead = helium%worm_bead_idx - 1, 1, -1
2499  DO idim = 1, 3
2500  xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
2501  helium%work(idim, helium%worm_atom_idx, kbead) = helium%work(idim, helium%worm_atom_idx, kbead + 1) + xr
2502  END DO
2503  END DO
2504 
2505  ! over imaginary time boundary
2506  DO idim = 1, 3
2507  xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
2508  helium%work(idim, helium%worm_atom_idx_work, helium%beads) = helium%work(idim, helium%worm_atom_idx, 1) + xr
2509  END DO
2510 
2511  ! rest on startatom
2512  DO kbead = helium%beads - 1, helium%worm_bead_idx_work, -1
2513  DO idim = 1, 3
2514  xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
2515  helium%work(idim, helium%worm_atom_idx_work, kbead) = helium%work(idim, helium%worm_atom_idx_work, kbead + 1) + xr
2516  END DO
2517  END DO
2518  END IF
2519 
2520  snew = worm_path_action(helium, helium%work, &
2521  helium%worm_atom_idx_work, helium%worm_bead_idx_work, &
2522  helium%worm_atom_idx, helium%worm_bead_idx)
2523 
2524  IF (helium%solute_present) THEN
2525  snew = snew + worm_path_inter_action_tail(pint_env, helium, helium%work, &
2526  helium%worm_atom_idx, helium%worm_bead_idx, &
2527  helium%worm_atom_idx_work, helium%worm_bead_idx_work)
2528  snew = snew + 0.5_dp*newheadpot + oldtailpot
2529  END IF
2530 
2531  ! Metropolis:
2532  ! action difference
2533  sdiff = sold - snew
2534  ! modify action difference due to extra bead
2535  IF (sdiff < 0) THEN
2536  should_reject = .false.
2537  IF (sdiff < -100.0_dp) THEN ! To protect from exponential underflow
2538  should_reject = .true.
2539  ELSE
2540  rtmp = helium%rng_stream_uniform%next()
2541  IF (exp(sdiff) < rtmp) THEN
2542  should_reject = .true.
2543  END IF
2544  END IF
2545  IF (should_reject) THEN
2546  ! rejected !
2547  ! write back only changed atoms
2548  ! transfer the new coordinates to work array
2549  IF (helium%worm_bead_idx_work < helium%worm_bead_idx) THEN
2550  ! everything belongs to the same atom
2551  DO kbead = helium%worm_bead_idx_work, helium%worm_bead_idx - 1
2552  helium%work(:, helium%worm_atom_idx, kbead) = helium%pos(:, helium%worm_atom_idx, kbead)
2553  END DO
2554  ELSE
2555  ! is distributed among two atoms
2556  ! transfer to atom not containing the tail bead
2557  DO kbead = helium%worm_bead_idx_work, helium%beads
2558  helium%work(:, helium%worm_atom_idx_work, kbead) = helium%pos(:, helium%worm_atom_idx_work, kbead)
2559  END DO
2560  ! transfer to atom containing the tail bead
2561  DO kbead = 1, helium%worm_bead_idx - 1
2562  helium%work(:, helium%worm_atom_idx, kbead) = helium%pos(:, helium%worm_atom_idx, kbead)
2563  END DO
2564  END IF
2565  helium%worm_xtra_bead_work(:) = helium%worm_xtra_bead(:)
2566  helium%worm_bead_idx_work = helium%worm_bead_idx
2567  helium%worm_atom_idx_work = helium%worm_atom_idx
2568  ac = 0
2569  RETURN
2570  END IF
2571  END IF
2572 
2573  ! for now accepted
2574  ! rejection of the whole move if any centroid is farther away
2575  ! from the current center of gravity than HELIUM%DROPLET_RADIUS [lwalewski]
2576  ! AND ist not moving towards the center
2577  IF (.NOT. helium%periodic) THEN
2578  IF (helium%solute_present) THEN
2579  new_com(:) = helium%center(:)
2580  old_com(:) = new_com(:)
2581  ELSE
2582  new_com(:) = 0.0_dp
2583  DO ia = 1, helium%atoms
2584  DO ib = 1, helium%beads
2585  new_com(:) = new_com(:) + helium%work(:, ia, ib)
2586  END DO
2587  END DO
2588  new_com(:) = new_com(:)/(real(helium%atoms*helium%beads, dp))
2589  ! also compute the old center of mass (optimization potential)
2590  old_com(:) = 0.0_dp
2591  DO ia = 1, helium%atoms
2592  DO ib = 1, helium%beads
2593  old_com(:) = old_com(:) + helium%pos(:, ia, ib)
2594  END DO
2595  END DO
2596  old_com(:) = old_com(:)/(real(helium%atoms*helium%beads, dp))
2597  END IF
2598  should_reject = .false.
2599  atom: DO ia = 1, helium%atoms
2600  dr(:) = 0.0_dp
2601  DO ib = 1, helium%beads
2602  dr(:) = dr(:) + helium%work(:, ia, ib) - new_com(:)
2603  END DO
2604  dr(:) = dr(:)/real(helium%beads, dp)
2605  rtmp = dot_product(dr, dr)
2606  IF (rtmp >= helium%droplet_radius**2) THEN
2607  dro(:) = 0.0_dp
2608  DO ib = 1, helium%beads
2609  dro(:) = dro(:) + helium%pos(:, ia, ib) - old_com(:)
2610  END DO
2611  dro(:) = dro(:)/real(helium%beads, dp)
2612  rtmpo = dot_product(dro, dro)
2613  ! only reject if it moves away from COM outside the droplet radius
2614  IF (rtmpo < rtmp) THEN
2615  ! found - this one does not fit within R from the center
2616  should_reject = .true.
2617  EXIT atom
2618  END IF
2619  END IF
2620  END DO atom
2621  IF (should_reject) THEN
2622  ! restore original coordinates
2623  ! write back only changed atoms
2624  ! write back only changed atoms
2625  ! transfer the new coordinates to work array
2626  IF (helium%worm_bead_idx_work < helium%worm_bead_idx) THEN
2627  ! everything belongs to the same atom
2628  DO kbead = helium%worm_bead_idx_work, helium%worm_bead_idx - 1
2629  helium%work(:, helium%worm_atom_idx, kbead) = helium%pos(:, helium%worm_atom_idx, kbead)
2630  END DO
2631  ELSE
2632  ! is distributed among two atoms
2633  ! transfer to atom not containing the tail bead
2634  DO kbead = helium%worm_bead_idx_work, helium%beads
2635  helium%work(:, helium%worm_atom_idx_work, kbead) = helium%pos(:, helium%worm_atom_idx_work, kbead)
2636  END DO
2637  ! transfer to atom containing the tail bead
2638  DO kbead = 1, helium%worm_bead_idx - 1
2639  helium%work(:, helium%worm_atom_idx, kbead) = helium%pos(:, helium%worm_atom_idx, kbead)
2640  END DO
2641  END IF
2642  helium%worm_xtra_bead_work(:) = helium%worm_xtra_bead(:)
2643  helium%worm_bead_idx_work = helium%worm_bead_idx
2644  helium%worm_atom_idx_work = helium%worm_atom_idx
2645  ac = 0
2646  RETURN
2647  END IF
2648  END IF
2649 
2650  ! finally accepted
2651  ac = 1
2652  ! write changed coordinates to position array
2653  ! write back only changed atoms
2654  IF (helium%worm_bead_idx_work < helium%worm_bead_idx) THEN
2655  ! everything belongs to the same atom
2656  DO kbead = helium%worm_bead_idx_work, helium%worm_bead_idx - 1
2657  helium%pos(:, helium%worm_atom_idx, kbead) = helium%work(:, helium%worm_atom_idx, kbead)
2658  END DO
2659  ELSE
2660  ! is distributed among two atoms
2661  ! transfer to atom not containing the tail bead
2662  DO kbead = helium%worm_bead_idx_work, helium%beads
2663  helium%pos(:, helium%worm_atom_idx_work, kbead) = helium%work(:, helium%worm_atom_idx_work, kbead)
2664  END DO
2665  ! transfer to atom containing the tail bead
2666  DO kbead = 1, helium%worm_bead_idx - 1
2667  helium%pos(:, helium%worm_atom_idx, kbead) = helium%work(:, helium%worm_atom_idx, kbead)
2668  END DO
2669  END IF
2670  helium%worm_xtra_bead(:) = helium%worm_xtra_bead_work(:)
2671  helium%worm_bead_idx = helium%worm_bead_idx_work
2672  helium%worm_atom_idx = helium%worm_atom_idx_work
2673 
2674  END SUBROUTINE worm_crawl_move_backward
2675 
2676 ! **************************************************************************************************
2677 !> \brief Free particle density matrix
2678 !> \param helium ...
2679 !> \param startpos ...
2680 !> \param endpos ...
2681 !> \param l ...
2682 !> \return ...
2683 !> \author fuhl
2684 ! **************************************************************************************************
2685  REAL(kind=dp) FUNCTION free_density_matrix(helium, startpos, endpos, l) RESULT(rho)
2686 
2687  TYPE(helium_solvent_type), INTENT(IN) :: helium
2688  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: startpos, endpos
2689  INTEGER, INTENT(IN) :: l
2690 
2691  REAL(kind=dp) :: distsq, prefac
2692  REAL(kind=dp), DIMENSION(3) :: dvec
2693 
2694  dvec(:) = startpos(:) - endpos(:)
2695  CALL helium_pbc(helium, dvec)
2696  distsq = dot_product(dvec, dvec)
2697 
2698  prefac = 0.5_dp/(helium%hb2m*real(l, dp)*helium%tau)
2699  rho = -prefac*distsq
2700  rho = exp(rho)
2701  rho = rho*(sqrt(prefac/pi))**3
2702 
2703  END FUNCTION free_density_matrix
2704 
2705 ! **************************************************************************************************
2706 !> \brief Swap move in open state to change permutation state
2707 !> \param pint_env ...
2708 !> \param helium ...
2709 !> \param natoms ...
2710 !> \param l ...
2711 !> \param ac ...
2712 !> \author fuhl
2713 ! **************************************************************************************************
2714  SUBROUTINE worm_swap_move(pint_env, helium, natoms, l, ac)
2715 
2716  TYPE(pint_env_type), INTENT(IN) :: pint_env
2717  TYPE(helium_solvent_type), INTENT(INOUT) :: helium
2718  INTEGER, INTENT(IN) :: natoms, l
2719  INTEGER, INTENT(OUT) :: ac
2720 
2721  INTEGER :: bendatom, bstartatom, endbead, &
2722  excludeatom, fendatom, fstartatom, ia, &
2723  iatom, ib, ibead, jbead, kbead, &
2724  startbead, tmpi
2725  LOGICAL :: should_reject
2726  REAL(kind=dp) :: backwarddensmatsum, forwarddensmatsum, &
2727  mass, newheadpotential, &
2728  oldheadpotential, rtmp, rtmpo, sdiff, &
2729  snew, sold
2730  REAL(kind=dp), DIMENSION(3) :: dr, dro, new_com, old_com
2731  REAL(kind=dp), DIMENSION(3, l) :: newsection
2732  REAL(kind=dp), DIMENSION(natoms) :: backwarddensmat, forwarddensmat
2733 
2734  ! first the endbead of the reconstruction is needed
2735  startbead = helium%worm_bead_idx
2736  endbead = helium%worm_bead_idx + l + 1
2737  fstartatom = helium%worm_atom_idx
2738  excludeatom = fstartatom
2739  ! compute the atomwise probabilities to be the worms swap partner
2740  ! Check if the imaginary time section belongs to two atoms
2741  IF (endbead > helium%beads) THEN
2742  endbead = endbead - helium%beads
2743  ! exclude atom is the one not to connect to because it will result in an unnatural state
2744  excludeatom = helium%permutation(excludeatom)
2745  END IF
2746  DO iatom = 1, excludeatom - 1
2747  forwarddensmat(iatom) = free_density_matrix(helium, helium%worm_xtra_bead(:), &
2748  helium%pos(:, iatom, endbead), l + 1)
2749  END DO
2750  forwarddensmat(excludeatom) = 0.0_dp
2751  DO iatom = excludeatom + 1, natoms
2752  forwarddensmat(iatom) = free_density_matrix(helium, helium%worm_xtra_bead(:), &
2753  helium%pos(:, iatom, endbead), l + 1)
2754  END DO
2755  forwarddensmatsum = sum(forwarddensmat)
2756  IF (forwarddensmatsum == 0.0_dp) THEN
2757  ac = 0
2758  RETURN
2759  END IF
2760 
2761  ! Select an atom with its corresponding probability
2762  rtmp = helium%rng_stream_uniform%next()*forwarddensmatsum
2763  fendatom = 1
2764  DO WHILE (rtmp >= forwarddensmat(fendatom))
2765  rtmp = rtmp - forwarddensmat(fendatom)
2766  fendatom = fendatom + 1
2767  END DO
2768  ! just for numerical safety
2769  fendatom = min(fendatom, natoms)
2770  IF (fendatom == excludeatom) THEN
2771  ac = 0
2772  RETURN
2773  END IF
2774 
2775  ! ensure detailed balance
2776  IF (endbead > startbead) THEN
2777  bstartatom = fendatom
2778  ELSE
2779  bstartatom = helium%iperm(fendatom)
2780  END IF
2781  bendatom = fendatom
2782 
2783  DO iatom = 1, excludeatom - 1
2784  backwarddensmat(iatom) = free_density_matrix(helium, &
2785  helium%pos(:, bstartatom, startbead), &
2786  helium%pos(:, iatom, endbead), l + 1)
2787  END DO
2788  backwarddensmat(excludeatom) = 0.0_dp
2789  DO iatom = excludeatom + 1, natoms
2790  backwarddensmat(iatom) = free_density_matrix(helium, &
2791  helium%pos(:, bstartatom, startbead), &
2792  helium%pos(:, iatom, endbead), l + 1)
2793  END DO
2794 
2795  backwarddensmatsum = sum(backwarddensmat)
2796 
2797  mass = helium%he_mass_au
2798 
2799  !compute action before the move
2800  sold = worm_path_action(helium, helium%pos, &
2801  bstartatom, startbead, fendatom, endbead)
2802 
2803  IF (helium%solute_present) THEN
2804  ! no special head treatment needed, as it will change due to swapping
2805  ! the worm gap and due to primitive coupling no cross bead terms are considered
2806  sold = sold + worm_path_inter_action(pint_env, helium, helium%pos, &
2807  bstartatom, startbead, fendatom, endbead)
2808  ! compute potential of old and new head here (only once, as later only a rescaling is necessary)
2809  CALL helium_bead_solute_e_f(pint_env, helium, &
2810  fstartatom, startbead, helium%worm_xtra_bead, &
2811  energy=oldheadpotential)
2812  oldheadpotential = oldheadpotential*helium%tau
2813 
2814  CALL helium_bead_solute_e_f(pint_env, helium, &
2815  bstartatom, startbead, helium%pos(:, bstartatom, startbead), &
2816  energy=newheadpotential)
2817  newheadpotential = newheadpotential*helium%tau
2818 
2819  sold = sold + 0.5_dp*oldheadpotential + newheadpotential
2820  END IF
2821 
2822  ! construct a new path connecting the start and endbead
2823  ! need to be the old coordinates due to reordering of the work coordinates
2824  CALL path_construct(helium, &
2825  helium%worm_xtra_bead(:), &
2826  helium%pos(:, fendatom, endbead), l, &
2827  newsection)
2828 
2829  !write new path segment to work array
2830  !first the part that is guaranteed to fit on the coorinates of startatom of the
2831  jbead = 1
2832  IF (startbead < endbead) THEN
2833  ! everything belongs to the same atom
2834  DO kbead = startbead + 1, endbead - 1
2835  helium%work(:, fstartatom, kbead) = newsection(:, jbead)
2836  jbead = jbead + 1
2837  END DO
2838  ELSE
2839  ! is distributed among two atoms
2840  ! transfer to atom not containing the head bead
2841  DO kbead = startbead + 1, helium%beads
2842  helium%work(:, fstartatom, kbead) = newsection(:, jbead)
2843  jbead = jbead + 1
2844  END DO
2845  ! rest to the second atom
2846  DO ibead = 1, endbead - 1
2847  helium%work(:, fendatom, ibead) = newsection(:, jbead)
2848  jbead = jbead + 1
2849  END DO
2850  END IF
2851 
2852  !exchange coordinates head coordinate first
2853  helium%work(:, helium%worm_atom_idx, helium%worm_bead_idx) = helium%worm_xtra_bead(:)
2854  helium%worm_xtra_bead_work(:) = helium%pos(:, bstartatom, startbead)
2855 
2856  ! move coordinates from former worm atom tail to new worm atom tail
2857  DO ib = startbead, helium%beads
2858  helium%work(:, bstartatom, ib) = helium%pos(:, fstartatom, ib)
2859  END DO
2860  ! need to copy the rest of pselatom to former worm atom
2861  IF (endbead > startbead) THEN
2862  DO ib = endbead, helium%beads
2863  helium%work(:, fstartatom, ib) = helium%pos(:, bstartatom, ib)
2864  END DO
2865  END IF
2866 
2867  !update permtable
2868  tmpi = helium%permutation(bstartatom)
2869  helium%permutation(bstartatom) = helium%permutation(fstartatom)
2870  helium%permutation(fstartatom) = tmpi
2871  !update inverse permtable
2872  DO ib = 1, SIZE(helium%permutation)
2873  helium%iperm(helium%permutation(ib)) = ib
2874  END DO
2875  helium%worm_atom_idx_work = bstartatom
2876  ! action after the move
2877 
2878  IF (endbead > startbead) THEN
2879  snew = worm_path_action(helium, helium%work, &
2880  fstartatom, startbead, fstartatom, endbead)
2881  IF (helium%solute_present) THEN
2882  ! no special head treatment needed, because a swap can't go over
2883  ! the worm gap and due to primitive coupling no cross bead terms are considered
2884  snew = snew + worm_path_inter_action(pint_env, helium, helium%work, &
2885  fstartatom, startbead, fstartatom, endbead)
2886 
2887  ! add the previously computed old and new head actions
2888  snew = snew + oldheadpotential + 0.5_dp*newheadpotential
2889  END IF
2890  ELSE
2891  snew = worm_path_action(helium, helium%work, &
2892  fstartatom, startbead, helium%permutation(fstartatom), endbead)
2893  IF (helium%solute_present) THEN
2894  ! no special head treatment needed, because a swap can't go over
2895  ! the worm gap and due to primitive coupling no cross bead terms are considered
2896  snew = snew + worm_path_inter_action(pint_env, helium, helium%work, &
2897  fstartatom, startbead, helium%permutation(fstartatom), endbead)
2898 
2899  ! add the previously computed old and new head actions
2900  snew = snew + oldheadpotential + 0.5_dp*newheadpotential
2901  END IF
2902  END IF
2903 
2904  ! Metropolis:
2905  sdiff = sold - snew
2906  sdiff = sdiff + log(forwarddensmatsum/backwarddensmatsum)
2907  IF (sdiff < 0) THEN
2908  should_reject = .false.
2909  IF (sdiff < -100.0_dp) THEN ! To protect from exponential underflow
2910  should_reject = .true.
2911  ELSE
2912  rtmp = helium%rng_stream_uniform%next()
2913  IF (exp(sdiff) < rtmp) THEN
2914  should_reject = .true.
2915  END IF
2916  END IF
2917  IF (should_reject) THEN
2918  ! rejected !
2919  ! write back only changed atoms
2920  DO kbead = startbead, helium%beads
2921  helium%work(:, bstartatom, kbead) = helium%pos(:, bstartatom, kbead)
2922  END DO
2923  DO kbead = startbead, helium%beads
2924  helium%work(:, fstartatom, kbead) = helium%pos(:, fstartatom, kbead)
2925  END DO
2926  helium%worm_xtra_bead_work(:) = helium%worm_xtra_bead(:)
2927  helium%worm_atom_idx_work = helium%worm_atom_idx
2928  IF (startbead > endbead) THEN
2929  DO kbead = 1, endbead
2930  helium%work(:, fendatom, kbead) = helium%pos(:, fendatom, kbead)
2931  END DO
2932  END IF
2933  ! reset permtable
2934  tmpi = helium%permutation(bstartatom)
2935  helium%permutation(bstartatom) = helium%permutation(fstartatom)
2936  helium%permutation(fstartatom) = tmpi
2937  !update inverse permtable
2938  DO ib = 1, SIZE(helium%permutation)
2939  helium%iperm(helium%permutation(ib)) = ib
2940  END DO
2941  ac = 0
2942  RETURN
2943  END IF
2944  END IF
2945 
2946  ! for now accepted
2947  ! rejection of the whole move if any centroid is farther away
2948  ! from the current center of gravity than HELIUM%DROPLET_RADIUS [lwalewski]
2949  ! AND ist not moving towards the center
2950  IF (.NOT. helium%periodic) THEN
2951  IF (helium%solute_present) THEN
2952  new_com(:) = helium%center(:)
2953  old_com(:) = new_com(:)
2954  ELSE
2955  new_com(:) = 0.0_dp
2956  DO ia = 1, helium%atoms
2957  DO ib = 1, helium%beads
2958  new_com(:) = new_com(:) + helium%work(:, ia, ib)
2959  END DO
2960  END DO
2961  new_com(:) = new_com(:)/(real(helium%atoms*helium%beads, dp))
2962  ! also compute the old center of mass (optimization potential)
2963  old_com(:) = 0.0_dp
2964  DO ia = 1, helium%atoms
2965  DO ib = 1, helium%beads
2966  old_com(:) = old_com(:) + helium%pos(:, ia, ib)
2967  END DO
2968  END DO
2969  old_com(:) = old_com(:)/(real(helium%atoms*helium%beads, dp))
2970  END IF
2971  should_reject = .false.
2972  atom: DO ia = 1, helium%atoms
2973  dr(:) = 0.0_dp
2974  DO ib = 1, helium%beads
2975  dr(:) = dr(:) + helium%work(:, ia, ib) - new_com(:)
2976  END DO
2977  dr(:) = dr(:)/real(helium%beads, dp)
2978  rtmp = dot_product(dr, dr)
2979  IF (rtmp >= helium%droplet_radius**2) THEN
2980  dro(:) = 0.0_dp
2981  DO ib = 1, helium%beads
2982  dro(:) = dro(:) + helium%pos(:, ia, ib) - old_com(:)
2983  END DO
2984  dro(:) = dro(:)/real(helium%beads, dp)
2985  rtmpo = dot_product(dro, dro)
2986  ! only reject if it moves away from COM outside the droplet radius
2987  IF (rtmpo < rtmp) THEN
2988  ! found - this one does not fit within R from the center
2989  should_reject = .true.
2990  EXIT atom
2991  END IF
2992  END IF
2993  END DO atom
2994  IF (should_reject) THEN
2995  ! write back only changed atoms
2996  DO kbead = startbead, helium%beads
2997  helium%work(:, bstartatom, kbead) = helium%pos(:, bstartatom, kbead)
2998  END DO
2999  DO kbead = startbead, helium%beads
3000  helium%work(:, fstartatom, kbead) = helium%pos(:, fstartatom, kbead)
3001  END DO
3002  helium%worm_xtra_bead_work(:) = helium%worm_xtra_bead(:)
3003  helium%worm_atom_idx_work = helium%worm_atom_idx
3004  IF (startbead > endbead) THEN
3005  DO kbead = 1, endbead
3006  helium%work(:, fendatom, kbead) = helium%pos(:, fendatom, kbead)
3007  END DO
3008  END IF
3009 
3010  ! reset permtable
3011  tmpi = helium%permutation(bstartatom)
3012  helium%permutation(bstartatom) = helium%permutation(fstartatom)
3013  helium%permutation(fstartatom) = tmpi
3014  !update inverse permtable
3015  DO ib = 1, SIZE(helium%permutation)
3016  helium%iperm(helium%permutation(ib)) = ib
3017  END DO
3018  ac = 0
3019  RETURN
3020  END IF
3021  END IF
3022 
3023  ! finally accepted
3024  ac = 1
3025  DO kbead = startbead, helium%beads
3026  helium%pos(:, bstartatom, kbead) = helium%work(:, bstartatom, kbead)
3027  END DO
3028  DO kbead = startbead, helium%beads
3029  helium%pos(:, fstartatom, kbead) = helium%work(:, fstartatom, kbead)
3030  END DO
3031  helium%worm_xtra_bead(:) = helium%worm_xtra_bead_work(:)
3032  helium%worm_atom_idx = helium%worm_atom_idx_work
3033  IF (startbead > endbead) THEN
3034  DO kbead = 1, endbead
3035  helium%pos(:, fendatom, kbead) = helium%work(:, fendatom, kbead)
3036  END DO
3037  END IF
3038 
3039  END SUBROUTINE worm_swap_move
3040 
3041 ! **************************************************************************************************
3042 !> \brief Action along path
3043 !> \param helium ...
3044 !> \param pos ...
3045 !> \param startatom ...
3046 !> \param startbead ...
3047 !> \param endatom ...
3048 !> \param endbead ...
3049 !> \return ...
3050 !> \author fuhl
3051 ! **************************************************************************************************
3052  REAL(kind=dp) FUNCTION worm_path_action(helium, pos, &
3053  startatom, startbead, endatom, endbead) RESULT(partaction)
3054 
3055  TYPE(helium_solvent_type), INTENT(INOUT) :: helium
3056  REAL(kind=dp), DIMENSION(:, :, :), INTENT(IN), &
3057  POINTER :: pos
3058  INTEGER, INTENT(IN) :: startatom, startbead, endatom, endbead
3059 
3060  INTEGER :: iatom, ibead
3061  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: work2, work3
3062  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: work
3063 
3064  ALLOCATE (work(3, helium%beads + 1), work2(helium%beads + 1), work3(SIZE(helium%uoffdiag, 1) + 1))
3065  partaction = 0.0_dp
3066  ! do action in two ways, depending
3067  ! if coordinates are not wrapping
3068  IF (startbead < endbead) THEN
3069  ! helium pair action
3070  ! every atom, with the one (or two) who got a resampling
3071  DO iatom = 1, helium%atoms
3072  ! avoid self interaction
3073  IF (iatom == startatom) cycle
3074  ! first the section up to the worm gap
3075  ! two less, because we need to work on the worm intersection separately
3076  DO ibead = startbead, endbead
3077  work(:, ibead + 1 - startbead) = pos(:, iatom, ibead) - pos(:, startatom, ibead)
3078  END DO
3079  partaction = partaction + helium_eval_chain(helium, work, endbead - startbead + 1, work2, work3)
3080  END DO
3081  ELSE !(startbead > endbead)
3082  ! helium pair action
3083  ! every atom, with the one (or two) who got a resampling
3084  DO iatom = 1, helium%atoms
3085  ! avoid self interaction
3086  IF (iatom == startatom) cycle
3087  ! first the section up to the worm gap
3088  ! two less, because we need to work on the worm intersection separately
3089  DO ibead = startbead, helium%beads
3090  work(:, ibead + 1 - startbead) = pos(:, iatom, ibead) - pos(:, startatom, ibead)
3091  END DO
3092  ! wrapping bond
3093  work(:, helium%beads + 2 - startbead) = pos(:, helium%permutation(iatom), 1) - pos(:, endatom, 1)
3094  partaction = partaction + helium_eval_chain(helium, work, helium%beads - startbead + 2, work2, work3)
3095  END DO
3096 
3097  ! ending atom
3098  DO iatom = 1, helium%atoms
3099  ! avoid self interaction
3100  IF (iatom == endatom) cycle
3101  !from first to endbead
3102  IF (endbead > 1) THEN
3103  DO ibead = 1, endbead
3104  work(:, ibead) = pos(:, iatom, ibead) - pos(:, endatom, ibead)
3105  END DO
3106  partaction = partaction + helium_eval_chain(helium, work, endbead, work2, work3)
3107  END IF
3108  END DO
3109 
3110  END IF
3111  DEALLOCATE (work, work2, work3)
3112 
3113  END FUNCTION worm_path_action
3114 
3115 ! **************************************************************************************************
3116 !> \brief Corrected path action for worm
3117 !> \param helium ...
3118 !> \param pos ...
3119 !> \param startatom ...
3120 !> \param startbead ...
3121 !> \param endatom ...
3122 !> \param endbead ...
3123 !> \param xtrapos ...
3124 !> \param worm_atom_idx ...
3125 !> \param worm_bead_idx ...
3126 !> \return ...
3127 !> \author fuhl
3128 ! **************************************************************************************************
3129  REAL(kind=dp) FUNCTION worm_path_action_worm_corrected(helium, pos, &
3130  startatom, startbead, endatom, endbead, xtrapos, worm_atom_idx, worm_bead_idx) RESULT(partaction)
3131 
3132  TYPE(helium_solvent_type), INTENT(INOUT) :: helium
3133  REAL(kind=dp), DIMENSION(:, :, :), INTENT(IN), &
3134  POINTER :: pos
3135  INTEGER, INTENT(IN) :: startatom, startbead, endatom, endbead
3136  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: xtrapos
3137  INTEGER, INTENT(IN) :: worm_atom_idx, worm_bead_idx
3138 
3139  INTEGER :: iatom, ibead
3140  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: work2, work3
3141  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: work
3142  REAL(kind=dp), DIMENSION(3) :: r, rp
3143 
3144  ALLOCATE (work(3, helium%beads + 1), work2(helium%beads + 1), work3(SIZE(helium%uoffdiag, 1) + 1))
3145  partaction = 0.0_dp
3146  ! do action in two ways, depending
3147  ! if coordinates are not wrapping
3148  IF (startbead < endbead) THEN
3149  ! helium pair action
3150  ! every atom, with the one (or two) who got a resampling
3151  DO iatom = 1, helium%atoms
3152  ! avoid self interaction
3153  IF (iatom == startatom) cycle
3154  ! first the section up to the worm gap
3155  ! two less, because we need to work on the worm intersection separately
3156  IF (worm_bead_idx - startbead > 1) THEN
3157  DO ibead = startbead, worm_bead_idx - 1
3158  work(:, ibead + 1 - startbead) = pos(:, iatom, ibead) - pos(:, startatom, ibead)
3159  END DO
3160  partaction = partaction + helium_eval_chain(helium, work, worm_bead_idx - startbead, work2, work3)
3161  END IF
3162 
3163  IF (endbead > worm_bead_idx) THEN
3164  DO ibead = worm_bead_idx, endbead
3165  work(:, ibead + 1 - worm_bead_idx) = pos(:, iatom, ibead) - pos(:, startatom, ibead)
3166  END DO
3167  partaction = partaction + helium_eval_chain(helium, work, endbead - worm_bead_idx + 1, work2, work3)
3168  END IF
3169  END DO
3170 
3171  IF (worm_atom_idx /= startatom) THEN
3172  DO iatom = 1, helium%atoms
3173  IF (iatom == startatom) cycle
3174  IF (iatom == worm_atom_idx) cycle
3175  r(:) = pos(:, iatom, worm_bead_idx - 1) - pos(:, startatom, worm_bead_idx - 1)
3176  rp(:) = pos(:, iatom, worm_bead_idx) - pos(:, startatom, worm_bead_idx)
3177  partaction = partaction + helium_eval_expansion(helium, r, rp, work3)
3178  END DO
3179  ! add pair action with worm
3180  r(:) = pos(:, startatom, worm_bead_idx - 1) - pos(:, worm_atom_idx, worm_bead_idx - 1)
3181  rp(:) = pos(:, startatom, worm_bead_idx) - xtrapos(:)
3182  partaction = partaction + helium_eval_expansion(helium, r, rp, work3)
3183  ELSE
3184  ! worm intersection
3185  DO iatom = 1, helium%atoms
3186  IF (iatom == startatom) cycle
3187  r(:) = pos(:, iatom, worm_bead_idx - 1) - pos(:, startatom, worm_bead_idx - 1)
3188  rp(:) = pos(:, iatom, worm_bead_idx) - xtrapos(:)
3189  partaction = partaction + helium_eval_expansion(helium, r, rp, work3)
3190  END DO
3191  END IF
3192  ELSE !(startbead > endbead)
3193  ! section wraps around the atom
3194  ! two cases: worm gap is before or after wrap
3195  IF (worm_bead_idx > startbead) THEN
3196  ! action terms up to worm gap on starting atom
3197  DO iatom = 1, helium%atoms
3198  ! avoid self interaction
3199  IF (iatom == startatom) cycle
3200  ! first the section up to the worm gap
3201  ! two less, because we need to work on the worm intersection separately
3202  IF (worm_bead_idx - startbead > 1) THEN
3203  DO ibead = startbead, worm_bead_idx - 1
3204  work(:, ibead + 1 - startbead) = pos(:, iatom, ibead) - pos(:, startatom, ibead)
3205  END DO
3206  partaction = partaction + helium_eval_chain(helium, work, worm_bead_idx - startbead, work2, work3)
3207  END IF
3208 
3209  ! up to the wrapping border
3210  DO ibead = worm_bead_idx, helium%beads
3211  work(:, ibead + 1 - worm_bead_idx) = pos(:, iatom, ibead) - pos(:, startatom, ibead)
3212  END DO
3213  ! wrapping bond
3214  work(:, helium%beads - worm_bead_idx + 2) = pos(:, helium%permutation(iatom), 1) - &
3215  pos(:, helium%permutation(startatom), 1)
3216  partaction = partaction + helium_eval_chain(helium, work, helium%beads - worm_bead_idx + 2, work2, work3)
3217 
3218  END DO
3219 
3220  ! ending atom
3221  DO iatom = 1, helium%atoms
3222  ! avoid self interaction
3223  IF (iatom == endatom) cycle
3224  !from first to endbead
3225  IF (endbead > 1) THEN
3226  DO ibead = 1, endbead
3227  work(:, ibead) = pos(:, iatom, ibead) - pos(:, endatom, ibead)
3228  END DO
3229  partaction = partaction + helium_eval_chain(helium, work, endbead, work2, work3)
3230  END IF
3231  END DO
3232 
3233  IF (worm_atom_idx /= startatom) THEN
3234  DO iatom = 1, helium%atoms
3235  IF (iatom == startatom) cycle
3236  IF (iatom == worm_atom_idx) cycle
3237  r(:) = pos(:, iatom, worm_bead_idx - 1) - pos(:, startatom, worm_bead_idx - 1)
3238  rp(:) = pos(:, iatom, worm_bead_idx) - pos(:, startatom, worm_bead_idx)
3239  partaction = partaction + helium_eval_expansion(helium, r, rp, work3)
3240  END DO
3241  ! add pair action with worm
3242  r(:) = pos(:, startatom, worm_bead_idx - 1) - pos(:, worm_atom_idx, worm_bead_idx - 1)
3243  rp(:) = pos(:, startatom, worm_bead_idx) - xtrapos(:)
3244  partaction = partaction + helium_eval_expansion(helium, r, rp, work3)
3245  ELSE
3246  ! worm intersection
3247  DO iatom = 1, helium%atoms
3248  IF (iatom == startatom) cycle
3249  r(:) = pos(:, iatom, worm_bead_idx - 1) - pos(:, startatom, worm_bead_idx - 1)
3250  rp(:) = pos(:, iatom, worm_bead_idx) - xtrapos(:)
3251  partaction = partaction + helium_eval_expansion(helium, r, rp, work3)
3252  END DO
3253  END IF
3254  ELSE !(worm_bead_idx < endbead)
3255  IF (worm_bead_idx /= 1) THEN
3256  DO iatom = 1, helium%atoms
3257  ! avoid self interaction
3258  IF (iatom == startatom) cycle
3259  ! first the section up to the end of the atom
3260  ! one less, because we need to work on the wrapping separately
3261  DO ibead = startbead, helium%beads
3262  work(:, ibead + 1 - startbead) = pos(:, iatom, ibead) - pos(:, startatom, ibead)
3263  END DO
3264  ! wrapping bond
3265  work(:, helium%beads - startbead + 2) = pos(:, helium%permutation(iatom), 1) - pos(:, helium%permutation(startatom), 1)
3266  partaction = partaction + helium_eval_chain(helium, work, helium%beads - startbead + 2, work2, work3)
3267  END DO
3268 
3269  ! ending atom
3270  DO iatom = 1, helium%atoms
3271  ! avoid self interaction
3272  IF (iatom == endatom) cycle
3273  !from first to two before the worm gap
3274  IF (worm_bead_idx > 2) THEN
3275  DO ibead = 1, worm_bead_idx - 1
3276  work(:, ibead) = pos(:, iatom, ibead) - pos(:, endatom, ibead)
3277  END DO
3278  partaction = partaction + helium_eval_chain(helium, work, worm_bead_idx - 1, work2, work3)
3279  END IF
3280 
3281  IF (endbead > worm_bead_idx) THEN
3282  DO ibead = worm_bead_idx, endbead
3283  work(:, ibead - worm_bead_idx + 1) = pos(:, iatom, ibead) - pos(:, endatom, ibead)
3284  END DO
3285  partaction = partaction + helium_eval_chain(helium, work, endbead - worm_bead_idx + 1, work2, work3)
3286  END IF
3287  END DO
3288 
3289  IF (worm_atom_idx /= endatom) THEN
3290  DO iatom = 1, helium%atoms
3291  IF (iatom == endatom) cycle
3292  IF (iatom == worm_atom_idx) cycle
3293  r(:) = pos(:, iatom, worm_bead_idx - 1) - pos(:, endatom, worm_bead_idx - 1)
3294  rp(:) = pos(:, iatom, worm_bead_idx) - pos(:, endatom, worm_bead_idx)
3295  partaction = partaction + helium_eval_expansion(helium, r, rp, work3)
3296  END DO
3297  ! add pair action with worm
3298  r(:) = pos(:, endatom, worm_bead_idx - 1) - pos(:, worm_atom_idx, worm_bead_idx - 1)
3299  rp(:) = pos(:, endatom, worm_bead_idx) - xtrapos(:)
3300  partaction = partaction + helium_eval_expansion(helium, r, rp, work3)
3301  ELSE
3302  ! worm intersection
3303  DO iatom = 1, helium%atoms
3304  IF (iatom == endatom) cycle
3305  r(:) = pos(:, iatom, worm_bead_idx - 1) - pos(:, endatom, worm_bead_idx - 1)
3306  rp(:) = pos(:, iatom, worm_bead_idx) - xtrapos(:)
3307  partaction = partaction + helium_eval_expansion(helium, r, rp, work3)
3308  END DO
3309  END IF
3310  ELSE !(worm_bead_idx == 1)
3311  ! special treatment if first bead is worm bead
3312  DO iatom = 1, helium%atoms
3313  ! avoid self interaction
3314  IF (iatom == startatom) cycle
3315  ! first the section up to the end of the atom
3316  ! one less, because we need to work on the wrapping separately
3317  IF (helium%beads > startbead) THEN
3318  DO ibead = startbead, helium%beads
3319  work(:, ibead - startbead + 1) = pos(:, iatom, ibead) - pos(:, startatom, ibead)
3320  END DO
3321  partaction = partaction + helium_eval_chain(helium, work, helium%beads - startbead + 1, work2, work3)
3322  END IF
3323  END DO
3324 
3325  ! ending atom
3326  DO iatom = 1, helium%atoms
3327  IF (endbead > 1) THEN
3328  DO ibead = 1, endbead
3329  work(:, ibead) = pos(:, iatom, ibead) - pos(:, endatom, ibead)
3330  END DO
3331  partaction = partaction + helium_eval_chain(helium, work, endbead, work2, work3)
3332  END IF
3333  END DO
3334 
3335  IF (worm_atom_idx /= endatom) THEN
3336  DO iatom = 1, helium%atoms
3337  IF (iatom == endatom) cycle
3338  IF (iatom == worm_atom_idx) cycle
3339  r(:) = pos(:, helium%iperm(iatom), helium%beads) - pos(:, startatom, helium%beads)
3340  rp(:) = pos(:, iatom, 1) - pos(:, endatom, 1)
3341  partaction = partaction + helium_eval_expansion(helium, r, rp, work3)
3342  END DO
3343  ! add pair action with worm
3344  r(:) = pos(:, startatom, helium%beads) - pos(:, helium%iperm(worm_atom_idx), helium%beads)
3345  rp(:) = pos(:, endatom, 1) - xtrapos(:)
3346  partaction = partaction + helium_eval_expansion(helium, r, rp, work3)
3347  ELSE
3348  ! worm intersection
3349  DO iatom = 1, helium%atoms
3350  IF (iatom == endatom) cycle
3351  r(:) = pos(:, helium%iperm(iatom), helium%beads) - pos(:, startatom, helium%beads)
3352  rp(:) = pos(:, iatom, 1) - xtrapos(:)
3353  partaction = partaction + helium_eval_expansion(helium, r, rp, work3)
3354  END DO
3355  END IF
3356  END IF
3357  END IF
3358  END IF
3359  DEALLOCATE (work, work2, work3)
3360 
3361  END FUNCTION worm_path_action_worm_corrected
3362 
3363 ! **************************************************************************************************
3364 !> \brief Path interaction
3365 !> \param pint_env ...
3366 !> \param helium ...
3367 !> \param pos ...
3368 !> \param startatom ...
3369 !> \param startbead ...
3370 !> \param endatom ...
3371 !> \param endbead ...
3372 !> \return ...
3373 !> \author fuhl
3374 ! **************************************************************************************************
3375  REAL(kind=dp) FUNCTION worm_path_inter_action(pint_env, helium, pos, &
3376  startatom, startbead, endatom, endbead) RESULT(partaction)
3377 
3378  TYPE(pint_env_type), INTENT(IN) :: pint_env
3379  TYPE(helium_solvent_type), INTENT(IN) :: helium
3380  REAL(kind=dp), DIMENSION(:, :, :), INTENT(IN), &
3381  POINTER :: pos
3382  INTEGER, INTENT(IN) :: startatom, startbead, endatom, endbead
3383 
3384  INTEGER :: ibead
3385  REAL(kind=dp) :: energy
3386 
3387  partaction = 0.0_dp
3388 
3389  ! helium interaction with the solute
3390  ! do action in two ways, depending
3391  ! if coordinates are not wrapping
3392  IF (startbead < endbead) THEN
3393 
3394  ! interaction is only beadwise due to primitive coupling
3395  ! startatom and endatom are the same
3396  DO ibead = startbead + 1, endbead - 1
3397  CALL helium_bead_solute_e_f(pint_env, helium, &
3398  startatom, ibead, pos(:, startatom, ibead), energy=energy)
3399  partaction = partaction + energy
3400  END DO
3401 
3402  ELSE !(startbead > endbead)
3403 
3404  ! interaction is only beadwise due to primitive coupling
3405  ! startatom and endatom are different
3406  DO ibead = startbead + 1, helium%beads
3407  CALL helium_bead_solute_e_f(pint_env, helium, &
3408  startatom, ibead, pos(:, startatom, ibead), energy=energy)
3409  partaction = partaction + energy
3410  END DO
3411  ! second atom after imaginary time wrap
3412  DO ibead = 1, endbead - 1
3413  CALL helium_bead_solute_e_f(pint_env, helium, &
3414  endatom, ibead, pos(:, endatom, ibead), energy=energy)
3415  partaction = partaction + energy
3416  END DO
3417  END IF
3418 
3419  partaction = partaction*helium%tau
3420 
3421  END FUNCTION worm_path_inter_action
3422 
3423 ! **************************************************************************************************
3424 !> \brief Path interaction for head
3425 !> \param pint_env ...
3426 !> \param helium ...
3427 !> \param pos ...
3428 !> \param startatom ...
3429 !> \param startbead ...
3430 !> \param xtrapos ...
3431 !> \param worm_atom_idx ...
3432 !> \param worm_bead_idx ...
3433 !> \return ...
3434 !> \author fuhl
3435 ! **************************************************************************************************
3436  REAL(kind=dp) FUNCTION worm_path_inter_action_head(pint_env, helium, pos, &
3437  startatom, startbead, xtrapos, worm_atom_idx, worm_bead_idx) RESULT(partaction)
3438 
3439  TYPE(pint_env_type), INTENT(IN) :: pint_env
3440  TYPE(helium_solvent_type), INTENT(IN) :: helium
3441  REAL(kind=dp), DIMENSION(:, :, :), INTENT(IN), &
3442  POINTER :: pos
3443  INTEGER, INTENT(IN) :: startatom, startbead
3444  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: xtrapos
3445  INTEGER, INTENT(IN) :: worm_atom_idx, worm_bead_idx
3446 
3447  INTEGER :: ibead
3448  REAL(kind=dp) :: energy
3449 
3450  partaction = 0.0_dp
3451 
3452  ! helium interaction with the solute
3453  ! if coordinates are not wrapping
3454  IF (startbead < worm_bead_idx) THEN
3455  DO ibead = startbead + 1, worm_bead_idx - 1
3456  CALL helium_bead_solute_e_f(pint_env, helium, &
3457  startatom, ibead, pos(:, startatom, ibead), energy=energy)
3458  partaction = partaction + energy
3459  END DO
3460  CALL helium_bead_solute_e_f(pint_env, helium, &
3461  startatom, ibead, xtrapos, energy=energy)
3462  partaction = partaction + 0.5_dp*energy
3463  ELSE !(startbead > worm_bead_idx)
3464  DO ibead = startbead + 1, helium%beads
3465  CALL helium_bead_solute_e_f(pint_env, helium, &
3466  startatom, ibead, pos(:, startatom, ibead), energy=energy)
3467  partaction = partaction + energy
3468  END DO
3469  DO ibead = 1, worm_bead_idx - 1
3470  CALL helium_bead_solute_e_f(pint_env, helium, &
3471  worm_atom_idx, ibead, pos(:, worm_atom_idx, ibead), energy=energy)
3472  partaction = partaction + energy
3473  END DO
3474  CALL helium_bead_solute_e_f(pint_env, helium, &
3475  worm_atom_idx, ibead, xtrapos, energy=energy)
3476  partaction = partaction + 0.5_dp*energy
3477  END IF
3478 
3479  partaction = partaction*helium%tau
3480 
3481  END FUNCTION worm_path_inter_action_head
3482 
3483 ! **************************************************************************************************
3484 !> \brief Path interaction for tail
3485 !> \param pint_env ...
3486 !> \param helium ...
3487 !> \param pos ...
3488 !> \param endatom ...
3489 !> \param endbead ...
3490 !> \param worm_atom_idx ...
3491 !> \param worm_bead_idx ...
3492 !> \return ...
3493 !> \author fuhl
3494 ! **************************************************************************************************
3495  REAL(kind=dp) FUNCTION worm_path_inter_action_tail(pint_env, helium, pos, &
3496  endatom, endbead, worm_atom_idx, worm_bead_idx) RESULT(partaction)
3497 
3498  TYPE(pint_env_type), INTENT(IN) :: pint_env
3499  TYPE(helium_solvent_type), INTENT(IN) :: helium
3500  REAL(kind=dp), DIMENSION(:, :, :), INTENT(IN), &
3501  POINTER :: pos
3502  INTEGER, INTENT(IN) :: endatom, endbead, worm_atom_idx, &
3503  worm_bead_idx
3504 
3505  INTEGER :: ibead
3506  REAL(kind=dp) :: energy
3507 
3508  partaction = 0.0_dp
3509 
3510  ! helium interaction with the solute
3511  ! if coordinates are not wrapping
3512  IF (worm_bead_idx < endbead) THEN
3513  CALL helium_bead_solute_e_f(pint_env, helium, &
3514  worm_atom_idx, worm_bead_idx, pos(:, worm_atom_idx, worm_bead_idx), energy=energy)
3515  partaction = partaction + 0.5_dp*energy
3516  DO ibead = worm_bead_idx + 1, endbead - 1
3517  CALL helium_bead_solute_e_f(pint_env, helium, &
3518  endatom, ibead, pos(:, endatom, ibead), energy=energy)
3519  partaction = partaction + energy
3520  END DO
3521  ELSE !(worm_bead_idx > endbead)
3522  CALL helium_bead_solute_e_f(pint_env, helium, &
3523  worm_atom_idx, worm_bead_idx, pos(:, worm_atom_idx, worm_bead_idx), energy=energy)
3524  partaction = partaction + 0.5_dp*energy
3525  DO ibead = worm_bead_idx + 1, helium%beads
3526  CALL helium_bead_solute_e_f(pint_env, helium, &
3527  worm_atom_idx, ibead, pos(:, worm_atom_idx, ibead), energy=energy)
3528  partaction = partaction + energy
3529  END DO
3530  DO ibead = 1, endbead - 1
3531  CALL helium_bead_solute_e_f(pint_env, helium, &
3532  endatom, ibead, pos(:, endatom, ibead), energy=energy)
3533  partaction = partaction + energy
3534  END DO
3535  END IF
3536 
3537  partaction = partaction*helium%tau
3538 
3539  END FUNCTION worm_path_inter_action_tail
3540 
3541 END MODULE helium_worm
atom
Definition: atom.F:9
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_calc_plength
pure subroutine, public helium_calc_plength(helium)
Calculate probability distribution of the permutation lengths.
Definition: helium_common.F:763
helium_common::helium_update_coord_system
pure subroutine, public helium_update_coord_system(helium, pint_env)
Set coordinate system, e.g. for RHO calculations.
Definition: helium_common.F:2335
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:448
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:352
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:74
helium_io
I/O subroutines for helium.
Definition: helium_io.F:13
helium_io::helium_write_line
subroutine, public helium_write_line(line)
Writes out a line of text to the default output unit.
Definition: helium_io.F:447
helium_types
Data types representing superfluid helium.
Definition: helium_types.F:15
helium_worm
Methods dealing with the canonical worm alogrithm.
Definition: helium_worm.F:13
helium_worm::helium_sample_worm
subroutine, public helium_sample_worm(helium, pint_env)
Main worm sampling routine.
Definition: helium_worm.F:50
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition: input_constants.F:17
input_constants::helium_forces_average
integer, parameter, public helium_forces_average
Definition: input_constants.F:330
input_constants::helium_forces_last
integer, parameter, public helium_forces_last
Definition: input_constants.F:330
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition: kinds.F:57
mathconstants
Definition of mathematical constants and functions.
Definition: mathconstants.F:16
mathconstants::pi
real(kind=dp), parameter, public pi
Definition: mathconstants.F:110
pint_types
Definition: pint_types.F:8