Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2023 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Methods 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 22 : 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 22 : CALL helium_update_coord_system(helium, pint_env)
62 :
63 22 : ncentratt = 0
64 22 : ncentracc = 0
65 22 : nstagatt = 0
66 22 : nstagacc = 0
67 22 : nopenatt = 0
68 22 : nopenacc = 0
69 22 : ncloseatt = 0
70 22 : ncloseacc = 0
71 22 : nmoveheadatt = 0
72 22 : nmoveheadacc = 0
73 22 : nmovetailatt = 0
74 22 : nmovetailacc = 0
75 22 : ncrawlfwdatt = 0
76 22 : ncrawlfwdacc = 0
77 22 : ncrawlbwdatt = 0
78 22 : ncrawlbwdacc = 0
79 22 : nswapatt = 0
80 22 : npswapacc = 0
81 22 : nswapacc = 0
82 22 : nstat = 0
83 22 : ntot = 0
84 88 : helium%proarea%inst(:) = 0.d0
85 88 : helium%prarea2%inst(:) = 0.d0
86 88 : helium%wnumber%inst(:) = 0.d0
87 88 : helium%wnmber2%inst(:) = 0.d0
88 88 : helium%mominer%inst(:) = 0.d0
89 2074 : IF (helium%solute_present) helium%force_avrg(:, :) = 0.0d0
90 242 : helium%energy_avrg(:) = 0.0d0
91 726 : helium%plength_avrg(:) = 0.0d0
92 22 : IF (helium%worm_max_open_cycles > 0) THEN
93 0 : helium%savepos(:, :, :) = helium%pos(:, :, :)
94 0 : helium%saveiperm(:) = helium%iperm(:)
95 0 : helium%savepermutation(:) = helium%permutation(:)
96 : END IF
97 :
98 22 : nMC = helium%iter_rot
99 22 : ncycle = 0
100 22 : IF (helium%worm_allow_open) THEN
101 : DO ! Exit criterion at the end of the loop
102 19011 : DO iMC = 1, nMC
103 18910 : imove = helium%rng_stream_uniform%next(1, helium%worm_all_limit)
104 18910 : IF (helium%worm_is_closed) THEN
105 4448 : IF ((imove >= helium%worm_centroid_min) .AND. (imove <= helium%worm_centroid_max)) THEN
106 : ! centroid move
107 137 : iatom = helium%rng_stream_uniform%next(1, helium%atoms)
108 : CALL worm_centroid_move(pint_env, helium, &
109 137 : iatom, helium%worm_centroid_drmax, ac)
110 137 : ncentratt = ncentratt + 1
111 137 : 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 4311 : ELSE IF ((imove >= helium%worm_centroid_max + 1) .AND. (imove <= helium%worm_open_close_min - 1)) THEN
115 : ! staging move
116 3766 : iatom = helium%rng_stream_uniform%next(1, helium%atoms)
117 3766 : ibead = helium%rng_stream_uniform%next(1, helium%beads)
118 3766 : staging_l = helium%rng_stream_uniform%next(2, helium%worm_staging_l)
119 : CALL worm_staging_move(pint_env, helium, &
120 3766 : iatom, ibead, staging_l, ac)
121 3766 : nstagatt = nstagatt + 1
122 3766 : nstagacc = nstagacc + ac
123 545 : ELSE IF ((imove >= helium%worm_open_close_min) .AND. (imove <= helium%worm_open_close_max)) THEN
124 : ! attempt opening of worm
125 545 : iatom = helium%rng_stream_uniform%next(1, helium%atoms)
126 545 : ibead = helium%rng_stream_uniform%next(1, helium%beads)
127 545 : staging_l = helium%rng_stream_uniform%next(2, helium%worm_staging_l)
128 : CALL worm_open_move(pint_env, helium, &
129 545 : iatom, ibead, staging_l, ac)
130 545 : nopenatt = nopenatt + 1
131 545 : nopenacc = nopenacc + ac
132 : ELSE
133 : ! this must not occur
134 0 : CPABORT("Undefined move selected in helium worm sampling!")
135 : END IF
136 : ELSE ! worm is open
137 14462 : IF ((imove >= helium%worm_centroid_min) .AND. (imove <= helium%worm_centroid_max)) THEN
138 : ! centroid move
139 366 : iatom = helium%rng_stream_uniform%next(1, helium%atoms)
140 : CALL worm_centroid_move(pint_env, helium, &
141 366 : iatom, helium%worm_centroid_drmax, ac)
142 366 : ncentratt = ncentratt + 1
143 366 : ncentracc = ncentracc + ac
144 14096 : ELSE IF ((imove >= helium%worm_staging_min) .AND. (imove <= helium%worm_staging_max)) THEN
145 : ! staging move
146 871 : iatom = helium%rng_stream_uniform%next(1, helium%atoms)
147 871 : ibead = helium%rng_stream_uniform%next(1, helium%beads)
148 871 : staging_l = helium%rng_stream_uniform%next(2, helium%worm_staging_l)
149 : CALL worm_staging_move(pint_env, helium, &
150 871 : iatom, ibead, staging_l, ac)
151 871 : nstagatt = nstagatt + 1
152 871 : nstagacc = nstagacc + ac
153 13225 : ELSE IF ((imove >= helium%worm_fcrawl_min) .AND. (imove <= helium%worm_fcrawl_max)) THEN
154 : ! crawl forward
155 2190 : DO icrawl = 1, helium%worm_repeat_crawl
156 1460 : staging_l = helium%rng_stream_uniform%next(2, helium%worm_staging_l)
157 : CALL worm_crawl_move_forward(pint_env, helium, &
158 1460 : staging_l, ac)
159 1460 : ncrawlfwdatt = ncrawlfwdatt + 1
160 2190 : ncrawlfwdacc = ncrawlfwdacc + ac
161 : END DO
162 12495 : ELSE IF ((imove >= helium%worm_bcrawl_min) .AND. (imove <= helium%worm_bcrawl_max)) THEN
163 : ! crawl backward
164 2391 : DO icrawl = 1, helium%worm_repeat_crawl
165 1594 : staging_l = helium%rng_stream_uniform%next(2, helium%worm_staging_l)
166 : CALL worm_crawl_move_backward(pint_env, helium, &
167 1594 : staging_l, ac)
168 1594 : ncrawlbwdatt = ncrawlbwdatt + 1
169 2391 : ncrawlbwdacc = ncrawlbwdacc + ac
170 : END DO
171 11698 : ELSE IF ((imove >= helium%worm_head_min) .AND. (imove <= helium%worm_head_max)) THEN
172 : ! move head
173 844 : staging_l = helium%rng_stream_uniform%next(2, helium%worm_staging_l)
174 : CALL worm_head_move(pint_env, helium, &
175 844 : staging_l, ac)
176 844 : nmoveheadatt = nmoveheadatt + 1
177 844 : nmoveheadacc = nmoveheadacc + ac
178 10854 : ELSE IF ((imove >= helium%worm_tail_min) .AND. (imove <= helium%worm_tail_max)) THEN
179 : ! move tail
180 864 : staging_l = helium%rng_stream_uniform%next(2, helium%worm_staging_l)
181 : CALL worm_tail_move(pint_env, helium, &
182 864 : staging_l, ac)
183 864 : nmovetailatt = nmovetailatt + 1
184 864 : nmovetailacc = nmovetailacc + ac
185 9990 : ELSE IF ((imove >= helium%worm_swap_min) .AND. (imove <= helium%worm_swap_max)) THEN
186 8332 : staging_l = helium%rng_stream_uniform%next(2, helium%worm_staging_l)
187 : CALL worm_swap_move(pint_env, helium, &
188 8332 : helium%atoms, staging_l, ac)
189 8332 : npswapacc = npswapacc + ac
190 8332 : nswapacc = nswapacc + ac
191 8332 : nswapatt = nswapatt + 1
192 1658 : ELSE IF ((imove >= helium%worm_open_close_min) .AND. (imove <= helium%worm_open_close_max)) THEN
193 : ! attempt closing of worm
194 1658 : staging_l = helium%rng_stream_uniform%next(2, helium%worm_staging_l)
195 : CALL worm_close_move(pint_env, helium, &
196 1658 : staging_l, ac)
197 1658 : ncloseatt = ncloseatt + 1
198 1658 : ncloseacc = ncloseacc + ac
199 : ELSE
200 : ! this must not occur
201 0 : 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 18910 : IF (helium%worm_is_closed) THEN
207 4448 : nstat = nstat + 1
208 4448 : IF (helium%solute_present) THEN
209 4224 : IF (helium%get_helium_forces == helium_forces_average) THEN
210 : !TODO needs proper averaging!
211 0 : CALL helium_solute_e_f(pint_env, helium, rtmp)
212 0 : helium%force_avrg(:, :) = helium%force_avrg(:, :) + helium%force_inst(:, :)
213 : END IF
214 : END IF
215 : END IF
216 19011 : ntot = ntot + 1
217 : END DO ! MC loop
218 :
219 101 : IF (helium%worm_is_closed) THEN
220 : EXIT
221 : ELSE
222 79 : ncycle = ncycle + 1
223 : ! reset position and permutation and start MC cycle again
224 : ! if stuck in open position for too long
225 79 : IF (helium%worm_max_open_cycles > 0 .AND. ncycle > helium%worm_max_open_cycles) THEN
226 0 : nMC = helium%iter_rot
227 0 : ncycle = 0
228 0 : helium%pos = helium%savepos
229 0 : helium%work = helium%pos
230 0 : helium%permutation = helium%savepermutation
231 0 : helium%iperm = helium%saveiperm
232 0 : CPWARN("Trapped in open state, reset helium to closed state from last MD step.")
233 : ELSE
234 79 : 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 0 : DO iMC = 1, nMC
240 0 : imove = helium%rng_stream_uniform%next(1, helium%worm_all_limit)
241 :
242 0 : IF ((imove >= helium%worm_centroid_min) .AND. (imove <= helium%worm_centroid_max)) THEN
243 : ! centroid move
244 0 : iatom = helium%rng_stream_uniform%next(1, helium%atoms)
245 : CALL worm_centroid_move(pint_env, helium, &
246 0 : iatom, helium%worm_centroid_drmax, ac)
247 0 : ncentratt = ncentratt + 1
248 0 : ncentracc = ncentracc + ac
249 0 : ELSE IF ((imove >= helium%worm_staging_min) .AND. (imove <= helium%worm_staging_max)) THEN
250 : ! staging move
251 0 : iatom = helium%rng_stream_uniform%next(1, helium%atoms)
252 0 : ibead = helium%rng_stream_uniform%next(1, helium%beads)
253 : CALL worm_staging_move(pint_env, helium, &
254 0 : iatom, ibead, helium%worm_staging_l, ac)
255 0 : nstagatt = nstagatt + 1
256 0 : nstagacc = nstagacc + ac
257 : ELSE
258 : ! this must not occur
259 0 : 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 0 : nstat = nstat + 1
264 0 : ntot = ntot + 1
265 0 : IF (helium%solute_present) THEN
266 0 : IF (helium%get_helium_forces == helium_forces_average) THEN
267 : ! TODO: needs proper averaging
268 0 : CALL helium_solute_e_f(pint_env, helium, rtmp)
269 0 : 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 22 : nmoveheadacc + nmovetailacc + ncrawlfwdacc + ncrawlbwdacc
278 : helium%num_accepted(2, 1) = ncentratt + nstagatt + nopenatt + ncloseatt + nswapatt + &
279 22 : nmoveheadatt + nmovetailatt + ncrawlfwdatt + ncrawlbwdatt
280 22 : 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 22 : CALL helium_calc_energy(helium, pint_env)
285 242 : helium%energy_avrg(:) = helium%energy_inst(:)*REAL(helium%worm_nstat, dp)
286 :
287 22 : CALL helium_calc_plength(helium)
288 726 : helium%plength_avrg(:) = helium%plength_inst(:)*REAL(helium%worm_nstat, dp)
289 :
290 : ! Calculate last_force
291 22 : IF (helium%solute_present) THEN
292 12 : IF (helium%get_helium_forces == helium_forces_last) THEN
293 12 : CALL helium_solute_e_f(pint_env, helium, rtmp)
294 2064 : helium%force_avrg(:, :) = helium%force_avrg(:, :) + helium%force_inst(:, :)
295 : END IF
296 : END IF
297 :
298 22 : IF (helium%worm_show_statistics) THEN
299 22 : WRITE (stmp, *) '--------------------------------------------------'
300 22 : CALL helium_write_line(stmp)
301 22 : WRITE (stmp, '(A10,F15.5,2I10)') 'Opening: ', &
302 22 : REAL(nopenacc, dp)/REAL(MAX(1, nopenatt), dp), &
303 44 : nopenacc, nopenatt
304 22 : CALL helium_write_line(stmp)
305 22 : WRITE (stmp, '(A10,F15.5,2I10)') 'Closing: ', &
306 22 : REAL(ncloseacc, dp)/REAL(MAX(1, ncloseatt), dp), &
307 44 : ncloseacc, ncloseatt
308 22 : CALL helium_write_line(stmp)
309 22 : WRITE (stmp, '(A10,F15.5,2I10)') 'Move Head: ', &
310 22 : REAL(nmoveheadacc, dp)/REAL(MAX(1, nmoveheadatt), dp), &
311 44 : nmoveheadacc, nmoveheadatt
312 22 : CALL helium_write_line(stmp)
313 22 : WRITE (stmp, '(A10,F15.5,2I10)') 'Move Tail: ', &
314 22 : REAL(nmovetailacc, dp)/REAL(MAX(1, nmovetailatt), dp), &
315 44 : nmovetailacc, nmovetailatt
316 22 : CALL helium_write_line(stmp)
317 22 : WRITE (stmp, '(A10,F15.5,2I10)') 'Crawl FWD: ', &
318 22 : REAL(ncrawlfwdacc, dp)/REAL(MAX(1, ncrawlfwdatt), dp), &
319 44 : ncrawlfwdacc, ncrawlfwdatt
320 22 : CALL helium_write_line(stmp)
321 22 : WRITE (stmp, '(A10,F15.5,2I10)') 'Crawl BWD: ', &
322 22 : REAL(ncrawlbwdacc, dp)/REAL(MAX(1, ncrawlbwdatt), dp), &
323 44 : ncrawlbwdacc, ncrawlbwdatt
324 22 : CALL helium_write_line(stmp)
325 22 : WRITE (stmp, '(A10,F15.5,2I10)') 'Staging: ', &
326 22 : REAL(nstagacc, dp)/REAL(MAX(1, nstagatt), dp), &
327 44 : nstagacc, nstagatt
328 22 : CALL helium_write_line(stmp)
329 22 : WRITE (stmp, '(A10,F15.5,2I10)') 'Centroid: ', &
330 22 : REAL(ncentracc, dp)/REAL(MAX(1, ncentratt), dp), &
331 44 : ncentracc, ncentratt
332 22 : CALL helium_write_line(stmp)
333 22 : WRITE (stmp, '(A10,F15.5,2I10)') 'Select: ', &
334 22 : REAL(npswapacc, dp)/REAL(MAX(1, nswapatt), dp), &
335 44 : npswapacc, nswapatt
336 22 : CALL helium_write_line(stmp)
337 22 : WRITE (stmp, '(A10,F15.5,2I10)') 'Swapping: ', &
338 22 : REAL(nswapacc, dp)/REAL(MAX(1, nswapatt), dp), &
339 44 : nswapacc, nswapatt
340 22 : CALL helium_write_line(stmp)
341 22 : WRITE (stmp, *) "Open State Probability: ", REAL(ntot - nstat, dp)/REAL(MAX(1, ntot), dp), ntot - nstat, ntot
342 22 : CALL helium_write_line(stmp)
343 22 : WRITE (stmp, *) "Closed State Probability: ", REAL(nstat, dp)/REAL(MAX(1, ntot), dp), nstat, ntot
344 22 : CALL helium_write_line(stmp)
345 : END IF
346 :
347 : !CALL center_pos(helium)
348 :
349 22 : 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 503 : 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 2012 : DO ic = 1, 3
374 1509 : rtmp = helium%rng_stream_uniform%next()
375 2012 : dr(ic) = (2.0_dp*rtmp - 1.0_dp)*drmax
376 : END DO
377 :
378 503 : IF (helium%worm_is_closed) THEN
379 137 : worm_in_moved_cycle = .FALSE.
380 : ! Perform move for first atom
381 2428 : DO ib = 1, helium%beads
382 9301 : helium%work(:, iatom, ib) = helium%work(:, iatom, ib) + dr(:)
383 : END DO
384 : ! move along permutation cycle
385 137 : jatom = helium%permutation(iatom)
386 137 : DO WHILE (jatom /= iatom)
387 0 : DO ib = 1, helium%beads
388 0 : helium%work(:, jatom, ib) = helium%work(:, jatom, ib) + dr(:)
389 : END DO
390 : ! next atom in chain
391 0 : jatom = helium%permutation(jatom)
392 : END DO
393 : ELSE ! worm is open
394 366 : 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 6258 : DO ib = 1, helium%beads
398 23934 : helium%work(:, iatom, ib) = helium%work(:, iatom, ib) + dr(:)
399 : END DO
400 : ! move along permutation cycle
401 366 : jatom = helium%permutation(iatom)
402 366 : DO WHILE (jatom /= iatom)
403 0 : DO ib = 1, helium%beads
404 0 : helium%work(:, jatom, ib) = helium%work(:, jatom, ib) + dr(:)
405 : END DO
406 0 : worm_in_moved_cycle = worm_in_moved_cycle .OR. (helium%worm_atom_idx == jatom)
407 : ! next atom in chain
408 0 : jatom = helium%permutation(jatom)
409 : END DO
410 : ! if atom contains had bead move that as well
411 366 : IF (worm_in_moved_cycle) THEN
412 56 : 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 503 : helium%worm_xtra_bead, worm_in_moved_cycle)
418 :
419 : snew = worm_centroid_move_action(helium, helium%work, iatom, &
420 503 : helium%worm_xtra_bead_work, worm_in_moved_cycle)
421 :
422 503 : IF (helium%solute_present) THEN
423 : sold = sold + worm_centroid_move_inter_action(pint_env, helium, helium%pos, iatom, &
424 488 : helium%worm_xtra_bead, worm_in_moved_cycle)
425 : snew = snew + worm_centroid_move_inter_action(pint_env, helium, helium%work, iatom, &
426 488 : helium%worm_xtra_bead_work, worm_in_moved_cycle)
427 : END IF
428 :
429 : ! Metropolis:
430 503 : sdiff = sold - snew
431 503 : IF (sdiff < 0) THEN
432 246 : should_reject = .FALSE.
433 246 : IF (sdiff < -100.0_dp) THEN ! To protect from exponential underflow
434 : should_reject = .TRUE.
435 : ELSE
436 246 : rtmp = helium%rng_stream_uniform%next()
437 246 : 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 172 : DO ib = 1, helium%beads
445 664 : helium%work(:, iatom, ib) = helium%pos(:, iatom, ib)
446 : END DO
447 : ! move along permutation cycle
448 8 : jatom = helium%permutation(iatom)
449 8 : DO WHILE (jatom /= iatom)
450 0 : DO ib = 1, helium%beads
451 0 : helium%work(:, jatom, ib) = helium%pos(:, jatom, ib)
452 : END DO
453 : ! next atom in chain
454 0 : jatom = helium%permutation(jatom)
455 : END DO
456 32 : helium%worm_xtra_bead_work(:) = helium%worm_xtra_bead(:)
457 8 : ac = 0
458 10 : 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 495 : IF (.NOT. helium%periodic) THEN
467 484 : IF (helium%solute_present) THEN
468 1936 : new_com(:) = helium%center(:)
469 484 : old_com(:) = new_com(:)
470 : ELSE
471 0 : new_com(:) = 0.0_dp
472 0 : DO ia = 1, helium%atoms
473 0 : DO ib = 1, helium%beads
474 0 : new_com(:) = new_com(:) + helium%work(:, ia, ib)
475 : END DO
476 : END DO
477 0 : new_com(:) = new_com(:)/(REAL(helium%atoms*helium%beads, dp))
478 : ! also compute the old center of mass (optimization potential)
479 0 : old_com(:) = 0.0_dp
480 0 : DO ia = 1, helium%atoms
481 0 : DO ib = 1, helium%beads
482 0 : old_com(:) = old_com(:) + helium%pos(:, ia, ib)
483 : END DO
484 : END DO
485 0 : old_com(:) = old_com(:)/(REAL(helium%atoms*helium%beads, dp))
486 : END IF
487 484 : should_reject = .FALSE.
488 15944 : atom: DO ia = 1, helium%atoms
489 15462 : dr(:) = 0.0_dp
490 262854 : DO ib = 1, helium%beads
491 1005030 : dr(:) = dr(:) + helium%work(:, ia, ib) - new_com(:)
492 : END DO
493 61848 : dr(:) = dr(:)/REAL(helium%beads, dp)
494 61848 : rtmp = DOT_PRODUCT(dr, dr)
495 15944 : IF (rtmp >= helium%droplet_radius**2) THEN
496 2 : dro(:) = 0.0_dp
497 34 : DO ib = 1, helium%beads
498 130 : dro(:) = dro(:) + helium%pos(:, ia, ib) - old_com(:)
499 : END DO
500 8 : dro(:) = dro(:)/REAL(helium%beads, dp)
501 8 : rtmpo = DOT_PRODUCT(dro, dro)
502 : ! only reject if it moves away from COM outside the droplet radius
503 2 : 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 484 : IF (should_reject) THEN
511 : ! restore original coordinates
512 : ! write back only changed atoms
513 34 : DO ib = 1, helium%beads
514 130 : helium%work(:, iatom, ib) = helium%pos(:, iatom, ib)
515 : END DO
516 : ! move along permutation cycle
517 2 : jatom = helium%permutation(iatom)
518 2 : DO WHILE (jatom /= iatom)
519 0 : DO ib = 1, helium%beads
520 0 : helium%work(:, jatom, ib) = helium%pos(:, jatom, ib)
521 : END DO
522 : ! next atom in chain
523 0 : jatom = helium%permutation(jatom)
524 : END DO
525 8 : helium%worm_xtra_bead_work(:) = helium%worm_xtra_bead(:)
526 2 : ac = 0
527 2 : RETURN
528 : END IF
529 : END IF
530 :
531 : ! finally accepted
532 493 : ac = 1
533 : ! write changed coordinates to position array
534 8480 : DO ib = 1, helium%beads
535 32441 : helium%pos(:, iatom, ib) = helium%work(:, iatom, ib)
536 : END DO
537 : ! move along permutation cycle
538 493 : jatom = helium%permutation(iatom)
539 493 : DO WHILE (jatom /= iatom)
540 0 : DO ib = 1, helium%beads
541 0 : helium%pos(:, jatom, ib) = helium%work(:, jatom, ib)
542 : END DO
543 : ! next atom in chain
544 0 : jatom = helium%permutation(jatom)
545 : END DO
546 1972 : 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 1006 : 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 1006 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: work2, work3
574 1006 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: work
575 : REAL(KIND=dp), DIMENSION(3) :: r, rp
576 :
577 7042 : 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 1006 : partaction = 0.0_dp
582 :
583 : ! pair of first atom with non-in-cycle atoms
584 1006 : jatom = iatom
585 33198 : DO ia = 1, helium%atoms
586 32192 : IF (ia == jatom) CYCLE
587 31186 : incycle = .FALSE.
588 31186 : katom = helium%permutation(jatom)
589 31186 : DO WHILE (katom /= jatom)
590 0 : IF (katom == ia) THEN
591 : incycle = .TRUE.
592 : EXIT
593 : END IF
594 0 : katom = helium%permutation(katom)
595 : END DO
596 31186 : IF (incycle) CYCLE
597 : ! if not in cycle, compute pair action
598 538532 : DO ib = 1, helium%beads
599 2060570 : work(:, ib) = pos(:, ia, ib) - pos(:, jatom, ib)
600 : END DO
601 : work(:, helium%beads + 1) = pos(:, helium%permutation(ia), 1) - &
602 124744 : pos(:, helium%permutation(jatom), 1)
603 33198 : 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 1006 : jatom = helium%permutation(iatom)
607 1006 : DO WHILE (jatom /= iatom)
608 0 : DO ia = 1, helium%atoms
609 0 : IF (ia == jatom) CYCLE
610 0 : incycle = .FALSE.
611 0 : katom = helium%permutation(jatom)
612 0 : DO WHILE (katom /= jatom)
613 0 : IF (katom == ia) THEN
614 : incycle = .TRUE.
615 : EXIT
616 : END IF
617 0 : katom = helium%permutation(katom)
618 : END DO
619 0 : IF (incycle) CYCLE
620 : ! if not in cycle, compute pair action
621 0 : DO ib = 1, helium%beads
622 0 : work(:, ib) = pos(:, ia, ib) - pos(:, jatom, ib)
623 : END DO
624 : work(:, helium%beads + 1) = pos(:, helium%permutation(ia), 1) - &
625 0 : pos(:, helium%permutation(jatom), 1)
626 0 : partaction = partaction + helium_eval_chain(helium, work, helium%beads + 1, work2, work3)
627 : END DO
628 0 : jatom = helium%permutation(jatom)
629 : END DO
630 : ! correct pair action for open worm configurations
631 1006 : IF (.NOT. helium%worm_is_closed) THEN
632 732 : wbead = helium%worm_bead_idx
633 732 : IF (winc) THEN
634 28 : IF (helium%worm_bead_idx == 1) THEN
635 : ! patom is the atom in front of the lone head bead
636 8 : patom = helium%iperm(helium%worm_atom_idx)
637 : ! go through all atoms, not in the cycle, and correct pair action
638 264 : DO ia = 1, helium%atoms
639 256 : IF (ia == helium%worm_atom_idx) CYCLE
640 248 : incycle = .FALSE.
641 248 : katom = helium%permutation(helium%worm_atom_idx)
642 248 : DO WHILE (katom /= helium%worm_atom_idx)
643 0 : IF (katom == ia) THEN
644 : incycle = .TRUE.
645 : EXIT
646 : END IF
647 0 : katom = helium%permutation(katom)
648 : END DO
649 248 : IF (incycle) CYCLE
650 : ! find other previous atom
651 : ! opatom is the atom in front of the first bead of the second atom
652 248 : opatom = helium%iperm(ia)
653 : ! if not in cycle, compute pair action
654 : ! subtract pair action for closed link
655 992 : r(:) = pos(:, helium%worm_atom_idx, 1) - pos(:, ia, 1)
656 992 : rp(:) = pos(:, patom, helium%beads) - pos(:, opatom, helium%beads)
657 248 : partaction = partaction - helium_eval_expansion(helium, r, rp, work3)
658 : ! and add corrected extra link
659 : ! rp stays the same
660 992 : r(:) = xtrapos(:) - pos(:, ia, 1)
661 264 : 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 660 : DO ia = 1, helium%atoms
667 640 : IF (ia == helium%worm_atom_idx) CYCLE
668 620 : incycle = .FALSE.
669 620 : katom = helium%permutation(helium%worm_atom_idx)
670 620 : DO WHILE (katom /= helium%worm_atom_idx)
671 0 : IF (katom == ia) THEN
672 : incycle = .TRUE.
673 : EXIT
674 : END IF
675 0 : katom = helium%permutation(katom)
676 : END DO
677 620 : IF (incycle) CYCLE
678 : ! if not in cycle, compute pair action
679 : ! subtract pair action for closed link
680 2480 : r(:) = pos(:, helium%worm_atom_idx, wbead) - pos(:, ia, wbead)
681 2480 : rp(:) = pos(:, helium%worm_atom_idx, wbead - 1) - pos(:, ia, wbead - 1)
682 620 : partaction = partaction - helium_eval_expansion(helium, r, rp, work3)
683 : ! and add corrected extra link
684 : ! rp stays the same
685 2480 : r(:) = xtrapos(:) - pos(:, ia, wbead)
686 660 : partaction = partaction + helium_eval_expansion(helium, r, rp, work3)
687 : END DO
688 : END IF
689 : ELSE ! worm is not the moved cycle
690 704 : IF (helium%worm_bead_idx == 1) THEN
691 : ! patom is the atom in front of the lone head bead
692 48 : patom = helium%iperm(helium%worm_atom_idx)
693 48 : opatom = helium%iperm(iatom)
694 : !correct action contribution for first atom in moved cycle
695 192 : r(:) = pos(:, helium%worm_atom_idx, 1) - pos(:, iatom, 1)
696 192 : rp(:) = pos(:, patom, helium%beads) - pos(:, opatom, helium%beads)
697 48 : partaction = partaction - helium_eval_expansion(helium, r, rp, work3)
698 : ! and add corrected extra link
699 : ! rp stays the same
700 192 : r(:) = xtrapos(:) - pos(:, iatom, 1)
701 48 : 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 48 : ia = helium%permutation(iatom)
704 48 : DO WHILE (ia /= iatom)
705 0 : opatom = helium%iperm(ia)
706 : ! subtract pair action for closed link
707 0 : r(:) = pos(:, helium%worm_atom_idx, 1) - pos(:, ia, 1)
708 0 : rp(:) = pos(:, patom, helium%beads) - pos(:, opatom, helium%beads)
709 0 : partaction = partaction - helium_eval_expansion(helium, r, rp, work3)
710 : ! and add corrected extra link
711 : ! rp stays the same
712 0 : r(:) = xtrapos(:) - pos(:, ia, 1)
713 0 : partaction = partaction + helium_eval_expansion(helium, r, rp, work3)
714 48 : 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 2624 : r(:) = pos(:, helium%worm_atom_idx, wbead) - pos(:, iatom, wbead)
720 2624 : rp(:) = pos(:, helium%worm_atom_idx, wbead - 1) - pos(:, iatom, wbead - 1)
721 656 : partaction = partaction - helium_eval_expansion(helium, r, rp, work3)
722 : ! and add corrected extra link
723 : ! rp stays the same
724 2624 : r(:) = xtrapos(:) - pos(:, iatom, wbead)
725 656 : 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 656 : ia = helium%permutation(iatom)
728 656 : DO WHILE (ia /= iatom)
729 : ! subtract pair action for closed link
730 0 : r(:) = pos(:, helium%worm_atom_idx, wbead) - pos(:, ia, wbead)
731 0 : rp(:) = pos(:, helium%worm_atom_idx, wbead - 1) - pos(:, ia, wbead - 1)
732 0 : partaction = partaction - helium_eval_expansion(helium, r, rp, work3)
733 : ! and add corrected extra link
734 : ! rp stays the same
735 0 : r(:) = xtrapos(:) - pos(:, ia, wbead)
736 0 : partaction = partaction + helium_eval_expansion(helium, r, rp, work3)
737 656 : ia = helium%permutation(ia)
738 : END DO
739 : END IF
740 : END IF
741 : END IF
742 1006 : DEALLOCATE (work, work2, work3)
743 :
744 1006 : 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 976 : 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 976 : partaction = 0.0_dp
772 : ! spcial treatment if worm is in moved cycle
773 976 : IF (winc) THEN
774 : ! first atom by hand
775 28 : jatom = iatom
776 : ! if it is worm atom it gets special treatment
777 28 : IF (jatom == helium%worm_atom_idx) THEN
778 : ! up to worm intersection
779 148 : DO jbead = 1, helium%worm_bead_idx - 1
780 : CALL helium_bead_solute_e_f(pint_env, helium, &
781 120 : jatom, jbead, pos(:, jatom, jbead), energy=energy)
782 148 : partaction = partaction + energy
783 : END DO
784 : ! head and tail each with 1/2 weight
785 28 : jbead = helium%worm_bead_idx
786 : ! tail
787 : CALL helium_bead_solute_e_f(pint_env, helium, &
788 28 : jatom, jbead, pos(:, jatom, jbead), energy=energy)
789 28 : partaction = partaction + 0.5_dp*energy
790 : ! head
791 : CALL helium_bead_solute_e_f(pint_env, helium, &
792 28 : jatom, jbead, xtrapos, energy=energy)
793 28 : partaction = partaction + 0.5_dp*energy
794 : ! rest of ring polymer
795 328 : DO jbead = helium%worm_bead_idx + 1, helium%beads
796 : CALL helium_bead_solute_e_f(pint_env, helium, &
797 300 : jatom, jbead, pos(:, jatom, jbead), energy=energy)
798 328 : partaction = partaction + energy
799 : END DO
800 : ELSE
801 0 : DO jbead = 1, helium%beads
802 : CALL helium_bead_solute_e_f(pint_env, helium, &
803 0 : jatom, jbead, pos(:, jatom, jbead), energy=energy)
804 0 : partaction = partaction + energy
805 : END DO
806 : END IF
807 28 : jatom = helium%permutation(iatom)
808 28 : DO WHILE (jatom /= iatom)
809 0 : IF (jatom == helium%worm_atom_idx) THEN
810 : ! up to worm intersection
811 0 : DO jbead = 1, helium%worm_bead_idx - 1
812 : CALL helium_bead_solute_e_f(pint_env, helium, &
813 0 : jatom, jbead, pos(:, jatom, jbead), energy=energy)
814 0 : partaction = partaction + energy
815 : END DO
816 : ! head and tail each with 1/2 weight
817 0 : jbead = helium%worm_bead_idx
818 : ! tail
819 : CALL helium_bead_solute_e_f(pint_env, helium, &
820 0 : jatom, jbead, pos(:, jatom, jbead), energy=energy)
821 0 : partaction = partaction + 0.5_dp*energy
822 : ! head
823 : CALL helium_bead_solute_e_f(pint_env, helium, &
824 0 : jatom, jbead, xtrapos, energy=energy)
825 0 : partaction = partaction + 0.5_dp*energy
826 : ! rest of ring polymer
827 0 : DO jbead = helium%worm_bead_idx + 1, helium%beads
828 : CALL helium_bead_solute_e_f(pint_env, helium, &
829 0 : jatom, jbead, pos(:, jatom, jbead), energy=energy)
830 0 : partaction = partaction + energy
831 : END DO
832 : ELSE
833 0 : DO jbead = 1, helium%beads
834 : CALL helium_bead_solute_e_f(pint_env, helium, &
835 0 : jatom, jbead, pos(:, jatom, jbead), energy=energy)
836 0 : partaction = partaction + energy
837 : END DO
838 : END IF
839 0 : jatom = helium%permutation(jatom)
840 : END DO
841 : ELSE ! worm not in moved cycle
842 : ! first atom by hand
843 948 : jatom = iatom
844 : ! if it is worm atom it gets special treatment
845 16116 : DO jbead = 1, helium%beads
846 : CALL helium_bead_solute_e_f(pint_env, helium, &
847 15168 : jatom, jbead, pos(:, jatom, jbead), energy=energy)
848 16116 : partaction = partaction + energy
849 : END DO
850 948 : jatom = helium%permutation(iatom)
851 948 : DO WHILE (jatom /= iatom)
852 0 : DO jbead = 1, helium%beads
853 : CALL helium_bead_solute_e_f(pint_env, helium, &
854 0 : jatom, jbead, pos(:, jatom, jbead), energy=energy)
855 0 : partaction = partaction + energy
856 : END DO
857 948 : jatom = helium%permutation(jatom)
858 : END DO
859 : END IF
860 :
861 976 : partaction = partaction*helium%tau
862 :
863 976 : 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 14615 : 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 14615 : imass = 1.0_dp/helium%he_mass_au
887 : ! dealing with periodicity
888 14615 : rs(:) = ri(:)
889 58460 : re(:) = rj(:) - rs(:)
890 14615 : CALL helium_pbc(helium, re)
891 58460 : re(:) = re(:) + rs(:)
892 :
893 : ! first construction by hand
894 : ! reusable weight factor 1/(l+1)
895 14615 : rk = REAL(l, dp)
896 14615 : weight = 1.0_dp/(rk + 1.0_dp)
897 : ! staging mass needed for modified variance
898 14615 : invstagemass = rk*weight*imass
899 : ! proposing new positions
900 58460 : DO idim = 1, 3
901 58460 : new_path(idim, 1) = helium%rng_stream_gaussian%next(variance=helium%tau*invstagemass)
902 : END DO
903 58460 : new_path(:, 1) = new_path(:, 1) + weight*(re(:) + rk*rs(:))
904 :
905 49559 : DO istage = 2, l
906 : ! reusable weight factor 1/(k+1)
907 34944 : rk = REAL(l - istage + 1, dp)
908 34944 : weight = 1.0_dp/(rk + 1.0_dp)
909 : ! staging mass needed for modified variance
910 34944 : invstagemass = rk*weight*imass
911 : ! proposing new positions
912 139776 : DO idim = 1, 3
913 139776 : new_path(idim, istage) = helium%rng_stream_gaussian%next(variance=helium%tau*invstagemass)
914 : END DO
915 154391 : new_path(:, istage) = new_path(:, istage) + weight*(rk*new_path(:, istage - 1) + re(:))
916 : END DO
917 :
918 14615 : 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 4637 : 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 4637 : REAL(KIND=dp), DIMENSION(3, l) :: newsection
942 :
943 4637 : ac = 0
944 4637 : endbead = startbead + l + 1
945 : ! Check if the imaginary time section belongs to two atoms
946 4637 : IF (endbead > helium%beads) THEN
947 1277 : endatom = helium%permutation(startatom)
948 1277 : endbead = endbead - helium%beads
949 : ELSE
950 3360 : endatom = startatom
951 : END IF
952 :
953 : ! check if the imaginary time section overlaps with the worm opening
954 4637 : 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 871 : (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 268 : 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 859 : IF (worm_overlap) THEN
976 : sold = worm_path_action_worm_corrected(helium, helium%pos, &
977 : startatom, startbead, endatom, endbead, &
978 256 : helium%worm_xtra_bead, helium%worm_atom_idx, helium%worm_bead_idx)
979 : ELSE
980 : sold = worm_path_action(helium, helium%pos, &
981 4369 : startatom, startbead, endatom, endbead)
982 : END IF
983 :
984 4625 : 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 4422 : 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 4625 : 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 4625 : jbead = 1
1000 18401 : DO ibead = startbead + 1, MIN(helium%beads, startbead + l)
1001 55104 : helium%work(:, startatom, ibead) = newsection(:, jbead)
1002 18401 : jbead = jbead + 1
1003 : END DO
1004 : ! transfer the rest of the beads to coordinates of endatom if necessary
1005 4625 : IF (helium%beads < startbead + l) THEN
1006 3441 : DO ibead = 1, endbead - 1
1007 9780 : helium%work(:, endatom, ibead) = newsection(:, jbead)
1008 3441 : jbead = jbead + 1
1009 : END DO
1010 : END IF
1011 :
1012 : ! action after the move
1013 4625 : IF (worm_overlap) THEN
1014 : snew = worm_path_action_worm_corrected(helium, helium%work, &
1015 : startatom, startbead, endatom, endbead, &
1016 256 : 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 4369 : startatom, startbead, endatom, endbead)
1020 : END IF
1021 :
1022 4625 : 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 4422 : startatom, startbead, endatom, endbead)
1027 : END IF
1028 :
1029 : ! Metropolis:
1030 4625 : sdiff = sold - snew
1031 4625 : IF (sdiff < 0) THEN
1032 2370 : should_reject = .FALSE.
1033 2370 : IF (sdiff < -100.0_dp) THEN ! To protect from exponential underflow
1034 : should_reject = .TRUE.
1035 : ELSE
1036 2369 : rtmp = helium%rng_stream_uniform%next()
1037 2369 : 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 80 : jbead = 1
1045 374 : DO ibead = startbead + 1, MIN(helium%beads, startbead + l)
1046 1176 : helium%work(:, startatom, ibead) = helium%pos(:, startatom, ibead)
1047 374 : jbead = jbead + 1
1048 : END DO
1049 : ! transfer the rest of the beads to coordinates of endatom if necessary
1050 80 : IF (helium%beads < startbead + l) THEN
1051 27 : DO ibead = 1, endbead - 1
1052 72 : helium%work(:, endatom, ibead) = helium%pos(:, endatom, ibead)
1053 27 : jbead = jbead + 1
1054 : END DO
1055 : END IF
1056 80 : ac = 0
1057 80 : 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 4545 : IF (.NOT. helium%periodic) THEN
1066 4380 : IF (helium%solute_present) THEN
1067 17520 : new_com(:) = helium%center(:)
1068 4380 : old_com(:) = new_com(:)
1069 : ELSE
1070 0 : new_com(:) = 0.0_dp
1071 0 : DO ia = 1, helium%atoms
1072 0 : DO ib = 1, helium%beads
1073 0 : new_com(:) = new_com(:) + helium%work(:, ia, ib)
1074 : END DO
1075 : END DO
1076 0 : new_com(:) = new_com(:)/(REAL(helium%atoms*helium%beads, dp))
1077 : ! also compute the old center of mass (optimization potential)
1078 0 : old_com(:) = 0.0_dp
1079 0 : DO ia = 1, helium%atoms
1080 0 : DO ib = 1, helium%beads
1081 0 : old_com(:) = old_com(:) + helium%pos(:, ia, ib)
1082 : END DO
1083 : END DO
1084 0 : old_com(:) = old_com(:)/(REAL(helium%atoms*helium%beads, dp))
1085 : END IF
1086 4380 : should_reject = .FALSE.
1087 144064 : atom: DO ia = 1, helium%atoms
1088 139706 : dr(:) = 0.0_dp
1089 2375002 : DO ib = 1, helium%beads
1090 9080890 : dr(:) = dr(:) + helium%work(:, ia, ib) - new_com(:)
1091 : END DO
1092 558824 : dr(:) = dr(:)/REAL(helium%beads, dp)
1093 558824 : rtmp = DOT_PRODUCT(dr, dr)
1094 144064 : IF (rtmp >= helium%droplet_radius**2) THEN
1095 22 : dro(:) = 0.0_dp
1096 374 : DO ib = 1, helium%beads
1097 1430 : dro(:) = dro(:) + helium%pos(:, ia, ib) - old_com(:)
1098 : END DO
1099 88 : dro(:) = dro(:)/REAL(helium%beads, dp)
1100 88 : rtmpo = DOT_PRODUCT(dro, dro)
1101 : ! only reject if it moves away from COM outside the droplet radius
1102 22 : 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 4380 : IF (should_reject) THEN
1110 : ! restore original coordinates
1111 : ! write back only changed atoms
1112 22 : jbead = 1
1113 80 : DO ibead = startbead + 1, MIN(helium%beads, startbead + l)
1114 232 : helium%work(:, startatom, ibead) = helium%pos(:, startatom, ibead)
1115 80 : jbead = jbead + 1
1116 : END DO
1117 : ! transfer the rest of the beads to coordinates of endatom if necessary
1118 22 : IF (helium%beads < startbead + l) THEN
1119 28 : DO ibead = 1, endbead - 1
1120 80 : helium%work(:, endatom, ibead) = helium%pos(:, endatom, ibead)
1121 28 : jbead = jbead + 1
1122 : END DO
1123 : END IF
1124 22 : ac = 0
1125 22 : RETURN
1126 : END IF
1127 : END IF
1128 :
1129 : ! finally accepted
1130 4523 : ac = 1
1131 : ! write changed coordinates to position array
1132 4523 : jbead = 1
1133 17947 : DO ibead = startbead + 1, MIN(helium%beads, startbead + l)
1134 53696 : helium%pos(:, startatom, ibead) = helium%work(:, startatom, ibead)
1135 17947 : jbead = jbead + 1
1136 : END DO
1137 : ! transfer the rest of the beads to coordinates of endatom if necessary
1138 4523 : IF (helium%beads < startbead + l) THEN
1139 3386 : DO ibead = 1, endbead - 1
1140 9628 : helium%pos(:, endatom, ibead) = helium%work(:, endatom, ibead)
1141 3386 : 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 545 : 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 545 : mass = helium%he_mass_au
1171 :
1172 : ! get index of the atom and bead, where the resampling of the head begins
1173 545 : IF (l < endbead) THEN
1174 : ! startbead belongs to the same atom
1175 438 : startatom = endatom
1176 438 : startbead = endbead - l
1177 : ELSE
1178 : ! startbead belongs to a different atom
1179 : ! find previous atom (assuming l < nbeads)
1180 107 : startatom = helium%iperm(endatom)
1181 107 : startbead = endbead + helium%beads - l
1182 : END IF
1183 : sold = worm_path_action(helium, helium%pos, &
1184 545 : startatom, startbead, endatom, endbead)
1185 :
1186 545 : 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 520 : helium%pos(:, endatom, endbead), endatom, endbead)
1193 : END IF
1194 :
1195 545 : helium%worm_is_closed = .FALSE.
1196 545 : helium%worm_atom_idx_work = endatom
1197 545 : helium%worm_bead_idx_work = endbead
1198 :
1199 : ! alternative grow with consecutive gaussians
1200 545 : IF (startbead < endbead) THEN
1201 : ! everything belongs to the same atom
1202 : ! gro head from startbead
1203 1507 : DO kbead = startbead + 1, endbead - 1
1204 4714 : DO idim = 1, 3
1205 3207 : xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
1206 4276 : helium%work(idim, startatom, kbead) = helium%work(idim, startatom, kbead - 1) + xr
1207 : END DO
1208 : END DO
1209 : ! last grow head bead
1210 1752 : DO idim = 1, 3
1211 1314 : xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
1212 1752 : helium%worm_xtra_bead_work(idim) = helium%work(idim, startatom, endbead - 1) + xr
1213 : END DO
1214 107 : ELSE IF (endbead /= 1) THEN
1215 : ! is distributed among two atoms
1216 : ! grow from startbead
1217 168 : DO kbead = startbead + 1, helium%beads
1218 441 : DO idim = 1, 3
1219 273 : xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
1220 364 : 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 308 : DO idim = 1, 3
1225 231 : xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
1226 308 : helium%work(idim, endatom, 1) = helium%work(idim, startatom, helium%beads) + xr
1227 : END DO
1228 : ! everything on endatom
1229 140 : DO kbead = 2, endbead - 1
1230 329 : DO idim = 1, 3
1231 189 : xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
1232 252 : helium%work(idim, endatom, kbead) = helium%work(idim, endatom, kbead - 1) + xr
1233 : END DO
1234 : END DO
1235 308 : DO idim = 1, 3
1236 231 : xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
1237 308 : 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 110 : DO kbead = startbead + 1, helium%beads
1243 350 : DO idim = 1, 3
1244 240 : xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
1245 320 : 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 120 : DO idim = 1, 3
1250 90 : xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
1251 120 : 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 545 : helium%worm_xtra_bead_work, helium%worm_atom_idx_work, helium%worm_bead_idx_work)
1258 :
1259 545 : IF (helium%solute_present) THEN
1260 : snew = snew + worm_path_inter_action_head(pint_env, helium, helium%work, &
1261 : startatom, startbead, &
1262 520 : 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 2180 : distvec(:) = helium%pos(:, startatom, startbead) - helium%pos(:, endatom, endbead)
1268 545 : CALL helium_pbc(helium, distvec)
1269 2180 : distsq = DOT_PRODUCT(distvec, distvec)
1270 : ! action difference
1271 545 : sdiff = sold - snew
1272 : ! modify action difference due to extra bead
1273 545 : sdiff = sdiff + distsq/(2.0_dp*helium%hb2m*REAL(l, dp)*helium%tau)
1274 545 : sdiff = sdiff + 1.5_dp*LOG(REAL(l, dp)*helium%tau)
1275 545 : sdiff = sdiff + helium%worm_ln_openclose_scale
1276 545 : IF (sdiff < 0) THEN
1277 363 : should_reject = .FALSE.
1278 363 : IF (sdiff < -100.0_dp) THEN ! To protect from exponential underflow
1279 : should_reject = .TRUE.
1280 : ELSE
1281 363 : rtmp = helium%rng_stream_uniform%next()
1282 363 : 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 169 : IF (startbead < endbead) THEN
1291 : ! everything belongs to the same atom
1292 385 : DO kbead = startbead + 1, endbead - 1
1293 1141 : 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 87 : DO kbead = startbead + 1, helium%beads
1299 240 : helium%work(:, startatom, kbead) = helium%pos(:, startatom, kbead)
1300 : END DO
1301 : ! transfer to atom containing the head bead
1302 74 : DO kbead = 1, endbead - 1
1303 188 : helium%work(:, endatom, kbead) = helium%pos(:, endatom, kbead)
1304 : END DO
1305 : END IF
1306 169 : helium%worm_is_closed = .TRUE.
1307 169 : ac = 0
1308 169 : 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 376 : IF (.NOT. helium%periodic) THEN
1317 362 : IF (helium%solute_present) THEN
1318 1448 : new_com(:) = helium%center(:)
1319 362 : old_com(:) = new_com(:)
1320 : ELSE
1321 0 : new_com(:) = 0.0_dp
1322 0 : DO ia = 1, helium%atoms
1323 0 : DO ib = 1, helium%beads
1324 0 : new_com(:) = new_com(:) + helium%work(:, ia, ib)
1325 : END DO
1326 : END DO
1327 0 : new_com(:) = new_com(:)/(REAL(helium%atoms*helium%beads, dp))
1328 : ! also compute the old center of mass (optimization potential)
1329 0 : old_com(:) = 0.0_dp
1330 0 : DO ia = 1, helium%atoms
1331 0 : DO ib = 1, helium%beads
1332 0 : old_com(:) = old_com(:) + helium%pos(:, ia, ib)
1333 : END DO
1334 : END DO
1335 0 : old_com(:) = old_com(:)/(REAL(helium%atoms*helium%beads, dp))
1336 : END IF
1337 362 : should_reject = .FALSE.
1338 11946 : atom: DO ia = 1, helium%atoms
1339 11584 : dr(:) = 0.0_dp
1340 196928 : DO ib = 1, helium%beads
1341 752960 : dr(:) = dr(:) + helium%work(:, ia, ib) - new_com(:)
1342 : END DO
1343 46336 : dr(:) = dr(:)/REAL(helium%beads, dp)
1344 46336 : rtmp = DOT_PRODUCT(dr, dr)
1345 11946 : IF (rtmp >= helium%droplet_radius**2) THEN
1346 0 : dro(:) = 0.0_dp
1347 0 : DO ib = 1, helium%beads
1348 0 : dro(:) = dro(:) + helium%pos(:, ia, ib) - old_com(:)
1349 : END DO
1350 0 : dro(:) = dro(:)/REAL(helium%beads, dp)
1351 0 : rtmpo = DOT_PRODUCT(dro, dro)
1352 : ! only reject if it moves away from COM outside the droplet radius
1353 0 : 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 362 : IF (should_reject) THEN
1361 : ! restore original coordinates
1362 : ! write back only changed atoms
1363 0 : IF (startbead < endbead) THEN
1364 : ! everything belongs to the same atom
1365 0 : DO kbead = startbead + 1, endbead - 1
1366 0 : 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 0 : DO kbead = startbead + 1, helium%beads
1372 0 : helium%work(:, startatom, kbead) = helium%pos(:, startatom, kbead)
1373 : END DO
1374 : ! transfer to atom containing the head bead
1375 0 : DO kbead = 1, endbead - 1
1376 0 : helium%work(:, endatom, kbead) = helium%pos(:, endatom, kbead)
1377 : END DO
1378 : END IF
1379 0 : 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 0 : ac = 0
1383 0 : RETURN
1384 : END IF
1385 : END IF
1386 :
1387 : ! finally accepted
1388 376 : ac = 1
1389 : ! write changed coordinates to position array
1390 376 : IF (startbead < endbead) THEN
1391 : ! everything belongs to the same atom
1392 1122 : DO kbead = startbead + 1, endbead - 1
1393 3573 : 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 191 : DO kbead = startbead + 1, helium%beads
1399 551 : helium%pos(:, startatom, kbead) = helium%work(:, startatom, kbead)
1400 : END DO
1401 : ! transfer to atom containing the head bead
1402 173 : DO kbead = 1, endbead - 1
1403 479 : helium%pos(:, endatom, kbead) = helium%work(:, endatom, kbead)
1404 : END DO
1405 : END IF
1406 1504 : helium%worm_xtra_bead(:) = helium%worm_xtra_bead_work(:)
1407 376 : helium%worm_atom_idx = helium%worm_atom_idx_work
1408 376 : 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 1658 : 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 1658 : REAL(KIND=dp), DIMENSION(3, l-1) :: newsection
1434 :
1435 1658 : mass = helium%he_mass_au
1436 :
1437 1658 : endatom = helium%worm_atom_idx
1438 1658 : endbead = helium%worm_bead_idx
1439 : ! get index of the atom and bead, where the resampling of the head begins
1440 1658 : IF (l < endbead) THEN
1441 : ! startbead belongs to the same atom
1442 1250 : startatom = endatom
1443 1250 : startbead = endbead - l
1444 : ELSE
1445 : ! startbead belongs to a different atom
1446 : ! find previous atom (assuming l < nbeads)
1447 408 : startatom = helium%iperm(endatom)
1448 408 : startbead = endbead + helium%beads - l
1449 : END IF
1450 : sold = worm_path_action_worm_corrected(helium, helium%pos, &
1451 : startatom, startbead, endatom, endbead, &
1452 1658 : helium%worm_xtra_bead, helium%worm_atom_idx, helium%worm_bead_idx)
1453 :
1454 1658 : IF (helium%solute_present) THEN
1455 : sold = sold + worm_path_inter_action_head(pint_env, helium, helium%pos, &
1456 : startatom, startbead, &
1457 1628 : 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 1658 : newsection)
1466 :
1467 : ! transfer the new coordinates to work array
1468 1658 : jbead = 1
1469 1658 : IF (startbead < endbead) THEN
1470 : ! everything belongs to the same atom
1471 4277 : DO kbead = startbead + 1, endbead - 1
1472 12108 : helium%work(:, endatom, kbead) = newsection(:, jbead)
1473 4277 : jbead = jbead + 1
1474 : END DO
1475 : ELSE
1476 : ! is distributed among two atoms
1477 : ! transfer to atom not containing the head bead
1478 954 : DO kbead = startbead + 1, helium%beads
1479 2184 : helium%work(:, startatom, kbead) = newsection(:, jbead)
1480 954 : jbead = jbead + 1
1481 : END DO
1482 : ! transfer to atom containing the head bead
1483 978 : DO kbead = 1, endbead - 1
1484 2280 : helium%work(:, endatom, kbead) = newsection(:, jbead)
1485 978 : jbead = jbead + 1
1486 : END DO
1487 : END IF
1488 :
1489 1658 : helium%worm_is_closed = .TRUE.
1490 :
1491 : snew = worm_path_action(helium, helium%work, &
1492 1658 : startatom, startbead, endatom, endbead)
1493 :
1494 1658 : 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 1628 : helium%work(:, endatom, endbead), endatom, endbead)
1501 : END IF
1502 :
1503 : ! Metropolis:
1504 : ! first compute ln of free density matrix
1505 6632 : distvec(:) = helium%pos(:, startatom, startbead) - helium%pos(:, endatom, endbead)
1506 1658 : CALL helium_pbc(helium, distvec)
1507 6632 : distsq = DOT_PRODUCT(distvec, distvec)
1508 : ! action difference
1509 1658 : sdiff = sold - snew
1510 : ! modify action difference due to extra bead
1511 1658 : sdiff = sdiff - distsq/(2.0_dp*helium%hb2m*REAL(l, dp)*helium%tau)
1512 1658 : sdiff = sdiff - 1.5_dp*LOG(REAL(l, dp)*helium%tau)
1513 1658 : sdiff = sdiff - helium%worm_ln_openclose_scale
1514 1658 : IF (sdiff < 0) THEN
1515 1501 : should_reject = .FALSE.
1516 1501 : IF (sdiff < -100.0_dp) THEN ! To protect from exponential underflow
1517 : should_reject = .TRUE.
1518 : ELSE
1519 1465 : rtmp = helium%rng_stream_uniform%next()
1520 1465 : 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 1282 : IF (startbead < endbead) THEN
1529 : ! everything belongs to the same atom
1530 3205 : DO kbead = startbead + 1, endbead - 1
1531 9952 : 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 746 : DO kbead = startbead + 1, helium%beads
1537 2006 : helium%work(:, startatom, kbead) = helium%pos(:, startatom, kbead)
1538 : END DO
1539 : ! transfer to atom containing the head bead
1540 789 : DO kbead = 1, endbead - 1
1541 2178 : helium%work(:, endatom, kbead) = helium%pos(:, endatom, kbead)
1542 : END DO
1543 : END IF
1544 1282 : helium%worm_is_closed = .FALSE.
1545 1282 : ac = 0
1546 1282 : 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 376 : IF (.NOT. helium%periodic) THEN
1555 362 : IF (helium%solute_present) THEN
1556 1448 : new_com(:) = helium%center(:)
1557 362 : old_com(:) = new_com(:)
1558 : ELSE
1559 0 : new_com(:) = 0.0_dp
1560 0 : DO ia = 1, helium%atoms
1561 0 : DO ib = 1, helium%beads
1562 0 : new_com(:) = new_com(:) + helium%work(:, ia, ib)
1563 : END DO
1564 : END DO
1565 0 : new_com(:) = new_com(:)/(REAL(helium%atoms*helium%beads, dp))
1566 : ! also compute the old center of mass (optimization potential)
1567 0 : old_com(:) = 0.0_dp
1568 0 : DO ia = 1, helium%atoms
1569 0 : DO ib = 1, helium%beads
1570 0 : old_com(:) = old_com(:) + helium%pos(:, ia, ib)
1571 : END DO
1572 : END DO
1573 0 : old_com(:) = old_com(:)/(REAL(helium%atoms*helium%beads, dp))
1574 : END IF
1575 362 : should_reject = .FALSE.
1576 11946 : atom: DO ia = 1, helium%atoms
1577 11584 : dr(:) = 0.0_dp
1578 196928 : DO ib = 1, helium%beads
1579 752960 : dr(:) = dr(:) + helium%work(:, ia, ib) - new_com(:)
1580 : END DO
1581 46336 : dr(:) = dr(:)/REAL(helium%beads, dp)
1582 46336 : rtmp = DOT_PRODUCT(dr, dr)
1583 11946 : IF (rtmp >= helium%droplet_radius**2) THEN
1584 0 : dro(:) = 0.0_dp
1585 0 : DO ib = 1, helium%beads
1586 0 : dro(:) = dro(:) + helium%pos(:, ia, ib) - old_com(:)
1587 : END DO
1588 0 : dro(:) = dro(:)/REAL(helium%beads, dp)
1589 0 : rtmpo = DOT_PRODUCT(dro, dro)
1590 : ! only reject if it moves away from COM outside the droplet radius
1591 0 : 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 362 : IF (should_reject) THEN
1599 : ! restore original coordinates
1600 : ! write back only changed atoms
1601 0 : IF (startbead < endbead) THEN
1602 : ! everything belongs to the same atom
1603 0 : DO kbead = startbead + 1, endbead - 1
1604 0 : 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 0 : DO kbead = startbead + 1, helium%beads
1610 0 : helium%work(:, startatom, kbead) = helium%pos(:, startatom, kbead)
1611 : END DO
1612 : ! transfer to atom containing the head bead
1613 0 : DO kbead = 1, endbead - 1
1614 0 : helium%work(:, endatom, kbead) = helium%pos(:, endatom, kbead)
1615 : END DO
1616 : END IF
1617 0 : helium%worm_is_closed = .FALSE.
1618 0 : ac = 0
1619 0 : RETURN
1620 : END IF
1621 : END IF
1622 :
1623 : ! finally accepted
1624 376 : ac = 1
1625 : ! write changed coordinates to position array
1626 376 : IF (startbead < endbead) THEN
1627 : ! everything belongs to the same atom
1628 1072 : DO kbead = startbead + 1, endbead - 1
1629 3406 : 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 208 : DO kbead = startbead + 1, helium%beads
1635 586 : helium%pos(:, startatom, kbead) = helium%work(:, startatom, kbead)
1636 : END DO
1637 : ! transfer to atom containing the head bead
1638 189 : DO kbead = 1, endbead - 1
1639 510 : 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 844 : 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 844 : mass = helium%he_mass_au
1667 :
1668 : ! get index of the atom and bead, where the resampling of the head begins
1669 844 : endatom = helium%worm_atom_idx
1670 844 : endbead = helium%worm_bead_idx
1671 844 : IF (l < endbead) THEN
1672 : ! startbead belongs to the same atom
1673 663 : startatom = endatom
1674 663 : startbead = endbead - l
1675 : ELSE
1676 : ! startbead belongs to a different atom
1677 : ! find previous atom (assuming l < nbeads)
1678 181 : startatom = helium%iperm(endatom)
1679 181 : 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 844 : helium%worm_xtra_bead, helium%worm_atom_idx, helium%worm_bead_idx)
1685 :
1686 844 : IF (helium%solute_present) THEN
1687 : sold = sold + worm_path_inter_action_head(pint_env, helium, helium%pos, &
1688 : startatom, startbead, &
1689 832 : helium%worm_xtra_bead, helium%worm_atom_idx, helium%worm_bead_idx)
1690 : END IF
1691 :
1692 : ! alternative grow with consecutive gaussians
1693 844 : IF (startbead < endbead) THEN
1694 : ! everything belongs to the same atom
1695 : ! gro head from startbead
1696 2284 : DO kbead = startbead + 1, endbead - 1
1697 7147 : DO idim = 1, 3
1698 4863 : xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
1699 6484 : helium%work(idim, startatom, kbead) = helium%work(idim, startatom, kbead - 1) + xr
1700 : END DO
1701 : END DO
1702 : ! last grow head bead
1703 2652 : DO idim = 1, 3
1704 1989 : xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
1705 2652 : helium%worm_xtra_bead_work(idim) = helium%work(idim, startatom, endbead - 1) + xr
1706 : END DO
1707 181 : ELSE IF (endbead /= 1) THEN
1708 : ! is distributed among two atoms
1709 : ! grow from startbead
1710 242 : DO kbead = startbead + 1, helium%beads
1711 623 : DO idim = 1, 3
1712 381 : xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
1713 508 : 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 460 : DO idim = 1, 3
1718 345 : xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
1719 460 : helium%work(idim, endatom, 1) = helium%work(idim, startatom, helium%beads) + xr
1720 : END DO
1721 : ! everything on endatom
1722 209 : DO kbead = 2, endbead - 1
1723 491 : DO idim = 1, 3
1724 282 : xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
1725 376 : helium%work(idim, endatom, kbead) = helium%work(idim, endatom, kbead - 1) + xr
1726 : END DO
1727 : END DO
1728 460 : DO idim = 1, 3
1729 345 : xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
1730 460 : 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 220 : DO kbead = startbead + 1, helium%beads
1736 682 : DO idim = 1, 3
1737 462 : xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
1738 616 : 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 264 : DO idim = 1, 3
1743 198 : xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
1744 264 : 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 844 : helium%worm_xtra_bead_work, helium%worm_atom_idx_work, helium%worm_bead_idx_work)
1751 :
1752 844 : IF (helium%solute_present) THEN
1753 : snew = snew + worm_path_inter_action_head(pint_env, helium, helium%work, &
1754 : startatom, startbead, &
1755 832 : 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 844 : sdiff = sold - snew
1761 : ! modify action difference due to extra bead
1762 844 : IF (sdiff < 0) THEN
1763 466 : should_reject = .FALSE.
1764 466 : IF (sdiff < -100.0_dp) THEN ! To protect from exponential underflow
1765 : should_reject = .TRUE.
1766 : ELSE
1767 459 : rtmp = helium%rng_stream_uniform%next()
1768 459 : 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 29 : IF (startbead < endbead) THEN
1777 : ! everything belongs to the same atom
1778 107 : DO kbead = startbead + 1, endbead - 1
1779 344 : 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 4 : DO kbead = startbead + 1, helium%beads
1785 13 : helium%work(:, startatom, kbead) = helium%pos(:, startatom, kbead)
1786 : END DO
1787 : ! transfer to atom containing the head bead
1788 2 : DO kbead = 1, endbead - 1
1789 5 : helium%work(:, endatom, kbead) = helium%pos(:, endatom, kbead)
1790 : END DO
1791 : END IF
1792 116 : helium%worm_xtra_bead_work(:) = helium%worm_xtra_bead(:)
1793 29 : ac = 0
1794 29 : 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 815 : IF (.NOT. helium%periodic) THEN
1803 808 : IF (helium%solute_present) THEN
1804 3232 : new_com(:) = helium%center(:)
1805 808 : old_com(:) = new_com(:)
1806 : ELSE
1807 0 : new_com(:) = 0.0_dp
1808 0 : DO ia = 1, helium%atoms
1809 0 : DO ib = 1, helium%beads
1810 0 : new_com(:) = new_com(:) + helium%work(:, ia, ib)
1811 : END DO
1812 : END DO
1813 0 : new_com(:) = new_com(:)/(REAL(helium%atoms*helium%beads, dp))
1814 : ! also compute the old center of mass (optimization potential)
1815 0 : old_com(:) = 0.0_dp
1816 0 : DO ia = 1, helium%atoms
1817 0 : DO ib = 1, helium%beads
1818 0 : old_com(:) = old_com(:) + helium%pos(:, ia, ib)
1819 : END DO
1820 : END DO
1821 0 : old_com(:) = old_com(:)/(REAL(helium%atoms*helium%beads, dp))
1822 : END IF
1823 808 : should_reject = .FALSE.
1824 26664 : atom: DO ia = 1, helium%atoms
1825 25856 : dr(:) = 0.0_dp
1826 439552 : DO ib = 1, helium%beads
1827 1680640 : dr(:) = dr(:) + helium%work(:, ia, ib) - new_com(:)
1828 : END DO
1829 103424 : dr(:) = dr(:)/REAL(helium%beads, dp)
1830 103424 : rtmp = DOT_PRODUCT(dr, dr)
1831 26664 : IF (rtmp >= helium%droplet_radius**2) THEN
1832 0 : dro(:) = 0.0_dp
1833 0 : DO ib = 1, helium%beads
1834 0 : dro(:) = dro(:) + helium%pos(:, ia, ib) - old_com(:)
1835 : END DO
1836 0 : dro(:) = dro(:)/REAL(helium%beads, dp)
1837 0 : rtmpo = DOT_PRODUCT(dro, dro)
1838 : ! only reject if it moves away from COM outside the droplet radius
1839 0 : 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 808 : IF (should_reject) THEN
1847 : ! restore original coordinates
1848 : ! write back only changed atoms
1849 0 : IF (startbead < endbead) THEN
1850 : ! everything belongs to the same atom
1851 0 : DO kbead = startbead + 1, endbead - 1
1852 0 : 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 0 : DO kbead = startbead + 1, helium%beads
1858 0 : helium%work(:, startatom, kbead) = helium%pos(:, startatom, kbead)
1859 : END DO
1860 : ! transfer to atom containing the head bead
1861 0 : DO kbead = 1, endbead - 1
1862 0 : helium%work(:, endatom, kbead) = helium%pos(:, endatom, kbead)
1863 : END DO
1864 : END IF
1865 0 : helium%worm_xtra_bead_work(:) = helium%worm_xtra_bead(:)
1866 0 : ac = 0
1867 0 : RETURN
1868 : END IF
1869 : END IF
1870 :
1871 : ! finally accepted
1872 815 : ac = 1
1873 : ! write changed coordinates to position array
1874 815 : IF (startbead < endbead) THEN
1875 : ! everything belongs to the same atom
1876 2177 : DO kbead = startbead + 1, endbead - 1
1877 6803 : 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 458 : DO kbead = startbead + 1, helium%beads
1883 1292 : helium%pos(:, startatom, kbead) = helium%work(:, startatom, kbead)
1884 : END DO
1885 : ! transfer to atom containing the head bead
1886 388 : DO kbead = 1, endbead - 1
1887 1012 : helium%pos(:, endatom, kbead) = helium%work(:, endatom, kbead)
1888 : END DO
1889 : END IF
1890 3260 : 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 864 : 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 864 : mass = helium%he_mass_au
1916 :
1917 : ! get index of the atom and bead, where the resampling of the tail ends
1918 864 : startatom = helium%worm_atom_idx
1919 864 : startbead = helium%worm_bead_idx
1920 864 : endbead = startbead + l
1921 :
1922 864 : IF (endbead <= helium%beads) THEN
1923 : ! endbead belongs to the same atom
1924 717 : endatom = startatom
1925 : ELSE
1926 : ! endbead belongs to a different atom
1927 : ! find next atom (assuming l < nbeads)
1928 147 : endatom = helium%permutation(startatom)
1929 147 : 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 864 : startatom, startbead, endatom, endbead)
1935 :
1936 864 : IF (helium%solute_present) THEN
1937 : sold = sold + worm_path_inter_action_tail(pint_env, helium, helium%pos, &
1938 : endatom, endbead, &
1939 852 : helium%worm_atom_idx, helium%worm_bead_idx)
1940 : END IF
1941 :
1942 : ! alternative grow with consecutive gaussians
1943 864 : IF (startbead < endbead) THEN
1944 : ! everything belongs to the same atom
1945 : ! gro tail from endbead to startbead (confusing eh?)
1946 3142 : DO kbead = endbead - 1, startbead, -1
1947 10417 : DO idim = 1, 3
1948 7275 : xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
1949 9700 : 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 304 : DO kbead = endbead - 1, 1, -1
1956 775 : DO idim = 1, 3
1957 471 : xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
1958 628 : 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 588 : DO idim = 1, 3
1964 441 : xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
1965 588 : helium%work(idim, startatom, helium%beads) = helium%work(idim, endatom, 1) + xr
1966 : END DO
1967 :
1968 : ! rest on startatom
1969 407 : DO kbead = helium%beads - 1, startbead, -1
1970 1187 : DO idim = 1, 3
1971 780 : xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
1972 1040 : 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 864 : startatom, startbead, endatom, endbead)
1980 :
1981 864 : IF (helium%solute_present) THEN
1982 : snew = snew + worm_path_inter_action_tail(pint_env, helium, helium%work, &
1983 : endatom, endbead, &
1984 852 : helium%worm_atom_idx_work, helium%worm_bead_idx_work)
1985 : END IF
1986 :
1987 : ! Metropolis:
1988 : ! action difference
1989 864 : sdiff = sold - snew
1990 : ! modify action difference due to extra bead
1991 864 : IF (sdiff < 0) THEN
1992 437 : should_reject = .FALSE.
1993 437 : IF (sdiff < -100.0_dp) THEN ! To protect from exponential underflow
1994 : should_reject = .TRUE.
1995 : ELSE
1996 437 : rtmp = helium%rng_stream_uniform%next()
1997 437 : 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 20 : IF (startbead < endbead) THEN
2006 : ! everything belongs to the same atom
2007 84 : DO kbead = startbead, endbead - 1
2008 288 : 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 22 : DO kbead = startbead, helium%beads
2014 76 : helium%work(:, startatom, kbead) = helium%pos(:, startatom, kbead)
2015 : END DO
2016 : ! transfer to atom containing the tail bead
2017 4 : DO kbead = 1, endbead - 1
2018 4 : helium%work(:, endatom, kbead) = helium%pos(:, endatom, kbead)
2019 : END DO
2020 : END IF
2021 80 : helium%worm_xtra_bead_work(:) = helium%worm_xtra_bead(:)
2022 20 : ac = 0
2023 20 : 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 844 : IF (.NOT. helium%periodic) THEN
2032 838 : IF (helium%solute_present) THEN
2033 3352 : new_com(:) = helium%center(:)
2034 838 : old_com(:) = new_com(:)
2035 : ELSE
2036 0 : new_com(:) = 0.0_dp
2037 0 : DO ia = 1, helium%atoms
2038 0 : DO ib = 1, helium%beads
2039 0 : new_com(:) = new_com(:) + helium%work(:, ia, ib)
2040 : END DO
2041 : END DO
2042 0 : new_com(:) = new_com(:)/(REAL(helium%atoms*helium%beads, dp))
2043 : ! also compute the old center of mass (optimization potential)
2044 0 : old_com(:) = 0.0_dp
2045 0 : DO ia = 1, helium%atoms
2046 0 : DO ib = 1, helium%beads
2047 0 : old_com(:) = old_com(:) + helium%pos(:, ia, ib)
2048 : END DO
2049 : END DO
2050 0 : old_com(:) = old_com(:)/(REAL(helium%atoms*helium%beads, dp))
2051 : END IF
2052 838 : should_reject = .FALSE.
2053 27654 : atom: DO ia = 1, helium%atoms
2054 26816 : dr(:) = 0.0_dp
2055 455872 : DO ib = 1, helium%beads
2056 1743040 : dr(:) = dr(:) + helium%work(:, ia, ib) - new_com(:)
2057 : END DO
2058 107264 : dr(:) = dr(:)/REAL(helium%beads, dp)
2059 107264 : rtmp = DOT_PRODUCT(dr, dr)
2060 27654 : IF (rtmp >= helium%droplet_radius**2) THEN
2061 0 : dro(:) = 0.0_dp
2062 0 : DO ib = 1, helium%beads
2063 0 : dro(:) = dro(:) + helium%pos(:, ia, ib) - old_com(:)
2064 : END DO
2065 0 : dro(:) = dro(:)/REAL(helium%beads, dp)
2066 0 : rtmpo = DOT_PRODUCT(dro, dro)
2067 : ! only reject if it moves away from COM outside the droplet radius
2068 0 : 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 838 : IF (should_reject) THEN
2076 : ! restore original coordinates
2077 : ! write back only changed atoms
2078 0 : IF (startbead < endbead) THEN
2079 : ! everything belongs to the same atom
2080 0 : DO kbead = startbead, endbead - 1
2081 0 : 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 0 : DO kbead = startbead, helium%beads
2087 0 : helium%work(:, startatom, kbead) = helium%pos(:, startatom, kbead)
2088 : END DO
2089 : ! transfer to atom containing the tail bead
2090 0 : DO kbead = 1, endbead - 1
2091 0 : helium%work(:, endatom, kbead) = helium%pos(:, endatom, kbead)
2092 : END DO
2093 : END IF
2094 0 : helium%worm_xtra_bead_work(:) = helium%worm_xtra_bead(:)
2095 0 : ac = 0
2096 0 : RETURN
2097 : END IF
2098 : END IF
2099 :
2100 : ! finally accepted
2101 844 : ac = 1
2102 : ! write changed coordinates to position array
2103 844 : IF (startbead < endbead) THEN
2104 : ! everything belongs to the same atom
2105 3058 : DO kbead = startbead, endbead - 1
2106 10129 : 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 532 : DO kbead = startbead, helium%beads
2112 1699 : helium%pos(:, startatom, kbead) = helium%work(:, startatom, kbead)
2113 : END DO
2114 : ! transfer to atom containing the tail bead
2115 300 : DO kbead = 1, endbead - 1
2116 771 : 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 1460 : 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 1460 : mass = helium%he_mass_au
2145 :
2146 : ! determine position of new head in imaginary time
2147 1460 : helium%worm_bead_idx_work = helium%worm_bead_idx + l
2148 1460 : IF (helium%worm_bead_idx_work > helium%beads) THEN
2149 322 : helium%worm_bead_idx_work = helium%worm_bead_idx_work - helium%beads
2150 322 : helium%worm_atom_idx_work = helium%permutation(helium%worm_atom_idx)
2151 : ELSE
2152 1138 : 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 1460 : helium%worm_atom_idx_work, helium%worm_bead_idx_work)
2160 1460 : 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 1428 : 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 1428 : energy=oldtailpot)
2173 1428 : 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 1428 : energy=newtailpot)
2180 1428 : 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 1428 : energy=oldheadpot)
2187 1428 : oldheadpot = oldheadpot*helium%tau
2188 :
2189 : ! new head is not known yet
2190 :
2191 1428 : 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 5840 : helium%work(:, helium%worm_atom_idx, helium%worm_bead_idx) = helium%worm_xtra_bead
2196 :
2197 : ! grow head part with consecutive gaussians
2198 1460 : IF (helium%worm_bead_idx < helium%worm_bead_idx_work) THEN
2199 : ! everything belongs to the same atom
2200 : ! gro head from startbead
2201 3901 : DO kbead = helium%worm_bead_idx + 1, helium%worm_bead_idx_work - 1
2202 12190 : DO idim = 1, 3
2203 8289 : xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
2204 11052 : 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 4552 : DO idim = 1, 3
2209 3414 : xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
2210 4552 : helium%worm_xtra_bead_work(idim) = helium%work(idim, helium%worm_atom_idx, helium%worm_bead_idx_work - 1) + xr
2211 : END DO
2212 322 : ELSE IF (helium%worm_bead_idx_work /= 1) THEN
2213 : ! is distributed among two atoms
2214 : ! grow from startbead
2215 445 : DO kbead = helium%worm_bead_idx + 1, helium%beads
2216 1099 : DO idim = 1, 3
2217 654 : xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
2218 872 : 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 908 : DO idim = 1, 3
2223 681 : xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
2224 908 : 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 435 : DO kbead = 2, helium%worm_bead_idx_work - 1
2228 1059 : DO idim = 1, 3
2229 624 : xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
2230 832 : 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 908 : DO idim = 1, 3
2234 681 : xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
2235 908 : 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 318 : DO kbead = helium%worm_bead_idx + 1, helium%beads
2241 987 : DO idim = 1, 3
2242 669 : xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
2243 892 : 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 380 : DO idim = 1, 3
2248 285 : xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
2249 380 : 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 1460 : helium%worm_xtra_bead_work, helium%worm_atom_idx_work, helium%worm_bead_idx_work)
2257 :
2258 1460 : 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 1428 : helium%worm_xtra_bead_work, helium%worm_atom_idx_work, helium%worm_bead_idx_work)
2262 1428 : snew = snew + 0.5_dp*newtailpot + oldheadpot
2263 : END IF
2264 :
2265 : ! Metropolis:
2266 : ! action difference
2267 1460 : sdiff = sold - snew
2268 : ! modify action difference due to extra bead
2269 1460 : IF (sdiff < 0) THEN
2270 743 : should_reject = .FALSE.
2271 743 : IF (sdiff < -100.0_dp) THEN ! To protect from exponential underflow
2272 : should_reject = .TRUE.
2273 : ELSE
2274 733 : rtmp = helium%rng_stream_uniform%next()
2275 733 : 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 84 : IF (helium%worm_bead_idx < helium%worm_bead_idx_work) THEN
2284 : ! everything belongs to the same atom
2285 287 : DO kbead = helium%worm_bead_idx, helium%worm_bead_idx_work - 1
2286 962 : 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 59 : DO kbead = helium%worm_bead_idx, helium%beads
2292 170 : 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 67 : DO kbead = 1, helium%worm_bead_idx_work - 1
2296 202 : helium%work(:, helium%worm_atom_idx_work, kbead) = helium%pos(:, helium%worm_atom_idx_work, kbead)
2297 : END DO
2298 : END IF
2299 336 : helium%worm_xtra_bead_work(:) = helium%worm_xtra_bead(:)
2300 84 : helium%worm_bead_idx_work = helium%worm_bead_idx
2301 84 : helium%worm_atom_idx_work = helium%worm_atom_idx
2302 84 : ac = 0
2303 90 : 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 1376 : IF (.NOT. helium%periodic) THEN
2312 1364 : IF (helium%solute_present) THEN
2313 5456 : new_com(:) = helium%center(:)
2314 1364 : old_com(:) = new_com(:)
2315 : ELSE
2316 0 : new_com(:) = 0.0_dp
2317 0 : DO ia = 1, helium%atoms
2318 0 : DO ib = 1, helium%beads
2319 0 : new_com(:) = new_com(:) + helium%work(:, ia, ib)
2320 : END DO
2321 : END DO
2322 0 : new_com(:) = new_com(:)/(REAL(helium%atoms*helium%beads, dp))
2323 : ! also compute the old center of mass (optimization potential)
2324 0 : old_com(:) = 0.0_dp
2325 0 : DO ia = 1, helium%atoms
2326 0 : DO ib = 1, helium%beads
2327 0 : old_com(:) = old_com(:) + helium%pos(:, ia, ib)
2328 : END DO
2329 : END DO
2330 0 : old_com(:) = old_com(:)/(REAL(helium%atoms*helium%beads, dp))
2331 : END IF
2332 1364 : should_reject = .FALSE.
2333 44904 : atom: DO ia = 1, helium%atoms
2334 43546 : dr(:) = 0.0_dp
2335 740282 : DO ib = 1, helium%beads
2336 2830490 : dr(:) = dr(:) + helium%work(:, ia, ib) - new_com(:)
2337 : END DO
2338 174184 : dr(:) = dr(:)/REAL(helium%beads, dp)
2339 174184 : rtmp = DOT_PRODUCT(dr, dr)
2340 44904 : IF (rtmp >= helium%droplet_radius**2) THEN
2341 6 : dro(:) = 0.0_dp
2342 102 : DO ib = 1, helium%beads
2343 390 : dro(:) = dro(:) + helium%pos(:, ia, ib) - old_com(:)
2344 : END DO
2345 24 : dro(:) = dro(:)/REAL(helium%beads, dp)
2346 24 : rtmpo = DOT_PRODUCT(dro, dro)
2347 : ! only reject if it moves away from COM outside the droplet radius
2348 6 : 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 1364 : IF (should_reject) THEN
2356 : ! restore original coordinates
2357 : ! write back only changed atoms
2358 6 : IF (helium%worm_bead_idx < helium%worm_bead_idx_work) THEN
2359 : ! everything belongs to the same atom
2360 26 : DO kbead = helium%worm_bead_idx, helium%worm_bead_idx_work - 1
2361 86 : 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 0 : DO kbead = helium%worm_bead_idx, helium%beads
2367 0 : 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 0 : DO kbead = 1, helium%worm_bead_idx_work - 1
2371 0 : helium%work(:, helium%worm_atom_idx_work, kbead) = helium%pos(:, helium%worm_atom_idx_work, kbead)
2372 : END DO
2373 : END IF
2374 24 : helium%worm_xtra_bead_work(:) = helium%worm_xtra_bead(:)
2375 6 : helium%worm_bead_idx_work = helium%worm_bead_idx
2376 6 : helium%worm_atom_idx_work = helium%worm_atom_idx
2377 6 : ac = 0
2378 6 : RETURN
2379 : END IF
2380 : END IF
2381 :
2382 : ! finally accepted
2383 1370 : ac = 1
2384 : ! write changed coordinates to position array
2385 1370 : IF (helium%worm_bead_idx < helium%worm_bead_idx_work) THEN
2386 : ! everything belongs to the same atom
2387 4726 : DO kbead = helium%worm_bead_idx, helium%worm_bead_idx_work - 1
2388 15694 : 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 1026 : DO kbead = helium%worm_bead_idx, helium%beads
2394 3204 : 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 690 : DO kbead = 1, helium%worm_bead_idx_work - 1
2398 1860 : helium%pos(:, helium%worm_atom_idx_work, kbead) = helium%work(:, helium%worm_atom_idx_work, kbead)
2399 : END DO
2400 : END IF
2401 5480 : helium%worm_xtra_bead(:) = helium%worm_xtra_bead_work(:)
2402 1370 : helium%worm_bead_idx = helium%worm_bead_idx_work
2403 1370 : 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 1594 : 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 1594 : mass = helium%he_mass_au
2430 :
2431 : ! determine position of new head in imaginary time
2432 1594 : helium%worm_bead_idx_work = helium%worm_bead_idx - l
2433 1594 : IF (helium%worm_bead_idx_work < 1) THEN
2434 342 : helium%worm_bead_idx_work = helium%worm_bead_idx_work + helium%beads
2435 342 : helium%worm_atom_idx_work = helium%iperm(helium%worm_atom_idx)
2436 : ELSE
2437 1252 : 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 1594 : helium%worm_xtra_bead, helium%worm_atom_idx, helium%worm_bead_idx)
2446 :
2447 1594 : 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 1560 : 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 1560 : energy=oldtailpot)
2460 1560 : 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 1560 : energy=oldheadpot)
2467 1560 : 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 1560 : energy=newheadpot)
2474 1560 : newheadpot = newheadpot*helium%tau
2475 :
2476 : ! new tail not known yet
2477 :
2478 1560 : sold = sold + 0.5_dp*(oldtailpot + oldheadpot) + newheadpot
2479 : END IF
2480 :
2481 : ! copy position to the head bead
2482 6376 : 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 1594 : 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 5568 : DO kbead = helium%worm_bead_idx - 1, helium%worm_bead_idx_work, -1
2489 18516 : DO idim = 1, 3
2490 12948 : xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
2491 17264 : 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 842 : DO kbead = helium%worm_bead_idx - 1, 1, -1
2499 2342 : DO idim = 1, 3
2500 1500 : xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
2501 2000 : 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 1368 : DO idim = 1, 3
2507 1026 : xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
2508 1368 : 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 822 : DO kbead = helium%beads - 1, helium%worm_bead_idx_work, -1
2513 2262 : DO idim = 1, 3
2514 1440 : xr = helium%rng_stream_gaussian%next(variance=helium%hb2m*helium%tau)
2515 1920 : 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 1594 : helium%worm_atom_idx, helium%worm_bead_idx)
2523 :
2524 1594 : 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 1560 : helium%worm_atom_idx_work, helium%worm_bead_idx_work)
2528 1560 : snew = snew + 0.5_dp*newheadpot + oldtailpot
2529 : END IF
2530 :
2531 : ! Metropolis:
2532 : ! action difference
2533 1594 : sdiff = sold - snew
2534 : ! modify action difference due to extra bead
2535 1594 : IF (sdiff < 0) THEN
2536 869 : should_reject = .FALSE.
2537 869 : IF (sdiff < -100.0_dp) THEN ! To protect from exponential underflow
2538 : should_reject = .TRUE.
2539 : ELSE
2540 868 : rtmp = helium%rng_stream_uniform%next()
2541 868 : 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 74 : IF (helium%worm_bead_idx_work < helium%worm_bead_idx) THEN
2550 : ! everything belongs to the same atom
2551 274 : DO kbead = helium%worm_bead_idx_work, helium%worm_bead_idx - 1
2552 922 : 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 44 : DO kbead = helium%worm_bead_idx_work, helium%beads
2558 128 : 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 54 : DO kbead = 1, helium%worm_bead_idx - 1
2562 168 : helium%work(:, helium%worm_atom_idx, kbead) = helium%pos(:, helium%worm_atom_idx, kbead)
2563 : END DO
2564 : END IF
2565 296 : helium%worm_xtra_bead_work(:) = helium%worm_xtra_bead(:)
2566 74 : helium%worm_bead_idx_work = helium%worm_bead_idx
2567 74 : helium%worm_atom_idx_work = helium%worm_atom_idx
2568 74 : ac = 0
2569 112 : 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 1520 : IF (.NOT. helium%periodic) THEN
2578 1506 : IF (helium%solute_present) THEN
2579 6024 : new_com(:) = helium%center(:)
2580 1506 : old_com(:) = new_com(:)
2581 : ELSE
2582 0 : new_com(:) = 0.0_dp
2583 0 : DO ia = 1, helium%atoms
2584 0 : DO ib = 1, helium%beads
2585 0 : new_com(:) = new_com(:) + helium%work(:, ia, ib)
2586 : END DO
2587 : END DO
2588 0 : new_com(:) = new_com(:)/(REAL(helium%atoms*helium%beads, dp))
2589 : ! also compute the old center of mass (optimization potential)
2590 0 : old_com(:) = 0.0_dp
2591 0 : DO ia = 1, helium%atoms
2592 0 : DO ib = 1, helium%beads
2593 0 : old_com(:) = old_com(:) + helium%pos(:, ia, ib)
2594 : END DO
2595 : END DO
2596 0 : old_com(:) = old_com(:)/(REAL(helium%atoms*helium%beads, dp))
2597 : END IF
2598 1506 : should_reject = .FALSE.
2599 49322 : atom: DO ia = 1, helium%atoms
2600 47854 : dr(:) = 0.0_dp
2601 813518 : DO ib = 1, helium%beads
2602 3110510 : dr(:) = dr(:) + helium%work(:, ia, ib) - new_com(:)
2603 : END DO
2604 191416 : dr(:) = dr(:)/REAL(helium%beads, dp)
2605 191416 : rtmp = DOT_PRODUCT(dr, dr)
2606 49322 : IF (rtmp >= helium%droplet_radius**2) THEN
2607 38 : dro(:) = 0.0_dp
2608 646 : DO ib = 1, helium%beads
2609 2470 : dro(:) = dro(:) + helium%pos(:, ia, ib) - old_com(:)
2610 : END DO
2611 152 : dro(:) = dro(:)/REAL(helium%beads, dp)
2612 152 : rtmpo = DOT_PRODUCT(dro, dro)
2613 : ! only reject if it moves away from COM outside the droplet radius
2614 38 : 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 1506 : 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 38 : IF (helium%worm_bead_idx_work < helium%worm_bead_idx) THEN
2627 : ! everything belongs to the same atom
2628 180 : DO kbead = helium%worm_bead_idx_work, helium%worm_bead_idx - 1
2629 606 : 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 0 : DO kbead = helium%worm_bead_idx_work, helium%beads
2635 0 : 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 0 : DO kbead = 1, helium%worm_bead_idx - 1
2639 0 : helium%work(:, helium%worm_atom_idx, kbead) = helium%pos(:, helium%worm_atom_idx, kbead)
2640 : END DO
2641 : END IF
2642 152 : helium%worm_xtra_bead_work(:) = helium%worm_xtra_bead(:)
2643 38 : helium%worm_bead_idx_work = helium%worm_bead_idx
2644 38 : helium%worm_atom_idx_work = helium%worm_atom_idx
2645 38 : ac = 0
2646 38 : RETURN
2647 : END IF
2648 : END IF
2649 :
2650 : ! finally accepted
2651 1482 : ac = 1
2652 : ! write changed coordinates to position array
2653 : ! write back only changed atoms
2654 1482 : IF (helium%worm_bead_idx_work < helium%worm_bead_idx) THEN
2655 : ! everything belongs to the same atom
2656 5114 : DO kbead = helium%worm_bead_idx_work, helium%worm_bead_idx - 1
2657 16988 : 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 1120 : DO kbead = helium%worm_bead_idx_work, helium%beads
2663 3502 : 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 788 : DO kbead = 1, helium%worm_bead_idx - 1
2667 2174 : helium%pos(:, helium%worm_atom_idx, kbead) = helium%work(:, helium%worm_atom_idx, kbead)
2668 : END DO
2669 : END IF
2670 5928 : helium%worm_xtra_bead(:) = helium%worm_xtra_bead_work(:)
2671 1482 : helium%worm_bead_idx = helium%worm_bead_idx_work
2672 1482 : 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 516584 : 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 2066336 : dvec(:) = startpos(:) - endpos(:)
2695 516584 : CALL helium_pbc(helium, dvec)
2696 2066336 : distsq = DOT_PRODUCT(dvec, dvec)
2697 :
2698 516584 : prefac = 0.5_dp/(helium%hb2m*REAL(l, dp)*helium%tau)
2699 516584 : rho = -prefac*distsq
2700 516584 : rho = EXP(rho)
2701 516584 : rho = rho*(SQRT(prefac/pi))**3
2702 :
2703 516584 : 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 8332 : 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 16664 : REAL(KIND=dp), DIMENSION(3, l) :: newsection
2732 8332 : REAL(KIND=dp), DIMENSION(natoms) :: backwarddensmat, forwarddensmat
2733 :
2734 : ! first the endbead of the reconstruction is needed
2735 8332 : startbead = helium%worm_bead_idx
2736 8332 : endbead = helium%worm_bead_idx + l + 1
2737 8332 : fstartatom = helium%worm_atom_idx
2738 8332 : 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 8332 : IF (endbead > helium%beads) THEN
2742 2178 : endbead = endbead - helium%beads
2743 : ! exclude atom is the one not to connect to because it will result in an unnatural state
2744 2178 : excludeatom = helium%permutation(excludeatom)
2745 : END IF
2746 135787 : DO iatom = 1, excludeatom - 1
2747 : forwarddensmat(iatom) = free_density_matrix(helium, helium%worm_xtra_bead(:), &
2748 135787 : helium%pos(:, iatom, endbead), l + 1)
2749 : END DO
2750 8332 : forwarddensmat(excludeatom) = 0.0_dp
2751 139169 : DO iatom = excludeatom + 1, natoms
2752 : forwarddensmat(iatom) = free_density_matrix(helium, helium%worm_xtra_bead(:), &
2753 139169 : helium%pos(:, iatom, endbead), l + 1)
2754 : END DO
2755 274956 : forwarddensmatsum = SUM(forwarddensmat)
2756 8332 : IF (forwarddensmatsum == 0.0_dp) THEN
2757 0 : ac = 0
2758 0 : RETURN
2759 : END IF
2760 :
2761 : ! Select an atom with its corresponding probability
2762 8332 : rtmp = helium%rng_stream_uniform%next()*forwarddensmatsum
2763 8332 : fendatom = 1
2764 148184 : DO WHILE (rtmp >= forwarddensmat(fendatom))
2765 139852 : rtmp = rtmp - forwarddensmat(fendatom)
2766 139852 : fendatom = fendatom + 1
2767 : END DO
2768 : ! just for numerical safety
2769 8332 : fendatom = MIN(fendatom, natoms)
2770 8332 : IF (fendatom == excludeatom) THEN
2771 0 : ac = 0
2772 0 : RETURN
2773 : END IF
2774 :
2775 : ! ensure detailed balance
2776 8332 : IF (endbead > startbead) THEN
2777 6154 : bstartatom = fendatom
2778 : ELSE
2779 2178 : bstartatom = helium%iperm(fendatom)
2780 : END IF
2781 : bendatom = fendatom
2782 :
2783 135787 : DO iatom = 1, excludeatom - 1
2784 : backwarddensmat(iatom) = free_density_matrix(helium, &
2785 : helium%pos(:, bstartatom, startbead), &
2786 135787 : helium%pos(:, iatom, endbead), l + 1)
2787 : END DO
2788 8332 : backwarddensmat(excludeatom) = 0.0_dp
2789 139169 : DO iatom = excludeatom + 1, natoms
2790 : backwarddensmat(iatom) = free_density_matrix(helium, &
2791 : helium%pos(:, bstartatom, startbead), &
2792 139169 : helium%pos(:, iatom, endbead), l + 1)
2793 : END DO
2794 :
2795 274956 : backwarddensmatsum = SUM(backwarddensmat)
2796 :
2797 8332 : mass = helium%he_mass_au
2798 :
2799 : !compute action before the move
2800 : sold = worm_path_action(helium, helium%pos, &
2801 8332 : bstartatom, startbead, fendatom, endbead)
2802 :
2803 8332 : 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 8152 : 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 8152 : energy=oldheadpotential)
2812 8152 : oldheadpotential = oldheadpotential*helium%tau
2813 :
2814 : CALL helium_bead_solute_e_f(pint_env, helium, &
2815 : bstartatom, startbead, helium%pos(:, bstartatom, startbead), &
2816 8152 : energy=newheadpotential)
2817 8152 : newheadpotential = newheadpotential*helium%tau
2818 :
2819 8152 : 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 8332 : 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 8332 : jbead = 1
2832 8332 : IF (startbead < endbead) THEN
2833 : ! everything belongs to the same atom
2834 27123 : DO kbead = startbead + 1, endbead - 1
2835 83876 : helium%work(:, fstartatom, kbead) = newsection(:, jbead)
2836 27123 : jbead = jbead + 1
2837 : END DO
2838 : ELSE
2839 : ! is distributed among two atoms
2840 : ! transfer to atom not containing the head bead
2841 6372 : DO kbead = startbead + 1, helium%beads
2842 16776 : helium%work(:, fstartatom, kbead) = newsection(:, jbead)
2843 6372 : jbead = jbead + 1
2844 : END DO
2845 : ! rest to the second atom
2846 6210 : DO ibead = 1, endbead - 1
2847 16128 : helium%work(:, fendatom, ibead) = newsection(:, jbead)
2848 6210 : jbead = jbead + 1
2849 : END DO
2850 : END IF
2851 :
2852 : !exchange coordinates head coordinate first
2853 33328 : helium%work(:, helium%worm_atom_idx, helium%worm_bead_idx) = helium%worm_xtra_bead(:)
2854 33328 : helium%worm_xtra_bead_work(:) = helium%pos(:, bstartatom, startbead)
2855 :
2856 : ! move coordinates from former worm atom tail to new worm atom tail
2857 81797 : DO ib = startbead, helium%beads
2858 302192 : helium%work(:, bstartatom, ib) = helium%pos(:, fstartatom, ib)
2859 : END DO
2860 : ! need to copy the rest of pselatom to former worm atom
2861 8332 : IF (endbead > startbead) THEN
2862 46124 : DO ib = endbead, helium%beads
2863 166034 : helium%work(:, fstartatom, ib) = helium%pos(:, bstartatom, ib)
2864 : END DO
2865 : END IF
2866 :
2867 : !update permtable
2868 8332 : tmpi = helium%permutation(bstartatom)
2869 8332 : helium%permutation(bstartatom) = helium%permutation(fstartatom)
2870 8332 : helium%permutation(fstartatom) = tmpi
2871 : !update inverse permtable
2872 274956 : DO ib = 1, SIZE(helium%permutation)
2873 274956 : helium%iperm(helium%permutation(ib)) = ib
2874 : END DO
2875 8332 : helium%worm_atom_idx_work = bstartatom
2876 : ! action after the move
2877 :
2878 8332 : IF (endbead > startbead) THEN
2879 : snew = worm_path_action(helium, helium%work, &
2880 6154 : fstartatom, startbead, fstartatom, endbead)
2881 6154 : 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 6060 : fstartatom, startbead, fstartatom, endbead)
2886 :
2887 : ! add the previously computed old and new head actions
2888 6060 : snew = snew + oldheadpotential + 0.5_dp*newheadpotential
2889 : END IF
2890 : ELSE
2891 : snew = worm_path_action(helium, helium%work, &
2892 2178 : fstartatom, startbead, helium%permutation(fstartatom), endbead)
2893 2178 : 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 2092 : fstartatom, startbead, helium%permutation(fstartatom), endbead)
2898 :
2899 : ! add the previously computed old and new head actions
2900 2092 : snew = snew + oldheadpotential + 0.5_dp*newheadpotential
2901 : END IF
2902 : END IF
2903 :
2904 : ! Metropolis:
2905 8332 : sdiff = sold - snew
2906 8332 : sdiff = sdiff + LOG(forwarddensmatsum/backwarddensmatsum)
2907 8332 : IF (sdiff < 0) THEN
2908 8308 : should_reject = .FALSE.
2909 8308 : IF (sdiff < -100.0_dp) THEN ! To protect from exponential underflow
2910 : should_reject = .TRUE.
2911 : ELSE
2912 4041 : rtmp = helium%rng_stream_uniform%next()
2913 4041 : 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 81315 : DO kbead = startbead, helium%beads
2921 300378 : helium%work(:, bstartatom, kbead) = helium%pos(:, bstartatom, kbead)
2922 : END DO
2923 81315 : DO kbead = startbead, helium%beads
2924 300378 : helium%work(:, fstartatom, kbead) = helium%pos(:, fstartatom, kbead)
2925 : END DO
2926 33176 : helium%worm_xtra_bead_work(:) = helium%worm_xtra_bead(:)
2927 8294 : helium%worm_atom_idx_work = helium%worm_atom_idx
2928 8294 : IF (startbead > endbead) THEN
2929 8388 : DO kbead = 1, endbead
2930 27018 : helium%work(:, fendatom, kbead) = helium%pos(:, fendatom, kbead)
2931 : END DO
2932 : END IF
2933 : ! reset permtable
2934 8294 : tmpi = helium%permutation(bstartatom)
2935 8294 : helium%permutation(bstartatom) = helium%permutation(fstartatom)
2936 8294 : helium%permutation(fstartatom) = tmpi
2937 : !update inverse permtable
2938 273702 : DO ib = 1, SIZE(helium%permutation)
2939 273702 : helium%iperm(helium%permutation(ib)) = ib
2940 : END DO
2941 8294 : ac = 0
2942 8294 : 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 38 : IF (.NOT. helium%periodic) THEN
2951 38 : IF (helium%solute_present) THEN
2952 152 : new_com(:) = helium%center(:)
2953 38 : old_com(:) = new_com(:)
2954 : ELSE
2955 0 : new_com(:) = 0.0_dp
2956 0 : DO ia = 1, helium%atoms
2957 0 : DO ib = 1, helium%beads
2958 0 : new_com(:) = new_com(:) + helium%work(:, ia, ib)
2959 : END DO
2960 : END DO
2961 0 : new_com(:) = new_com(:)/(REAL(helium%atoms*helium%beads, dp))
2962 : ! also compute the old center of mass (optimization potential)
2963 0 : old_com(:) = 0.0_dp
2964 0 : DO ia = 1, helium%atoms
2965 0 : DO ib = 1, helium%beads
2966 0 : old_com(:) = old_com(:) + helium%pos(:, ia, ib)
2967 : END DO
2968 : END DO
2969 0 : old_com(:) = old_com(:)/(REAL(helium%atoms*helium%beads, dp))
2970 : END IF
2971 38 : should_reject = .FALSE.
2972 1074 : atom: DO ia = 1, helium%atoms
2973 1066 : dr(:) = 0.0_dp
2974 18122 : DO ib = 1, helium%beads
2975 69290 : dr(:) = dr(:) + helium%work(:, ia, ib) - new_com(:)
2976 : END DO
2977 4264 : dr(:) = dr(:)/REAL(helium%beads, dp)
2978 4264 : rtmp = DOT_PRODUCT(dr, dr)
2979 1074 : IF (rtmp >= helium%droplet_radius**2) THEN
2980 30 : dro(:) = 0.0_dp
2981 510 : DO ib = 1, helium%beads
2982 1950 : dro(:) = dro(:) + helium%pos(:, ia, ib) - old_com(:)
2983 : END DO
2984 120 : dro(:) = dro(:)/REAL(helium%beads, dp)
2985 120 : rtmpo = DOT_PRODUCT(dro, dro)
2986 : ! only reject if it moves away from COM outside the droplet radius
2987 30 : 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 38 : IF (should_reject) THEN
2995 : ! write back only changed atoms
2996 380 : DO kbead = startbead, helium%beads
2997 1430 : helium%work(:, bstartatom, kbead) = helium%pos(:, bstartatom, kbead)
2998 : END DO
2999 380 : DO kbead = startbead, helium%beads
3000 1430 : helium%work(:, fstartatom, kbead) = helium%pos(:, fstartatom, kbead)
3001 : END DO
3002 120 : helium%worm_xtra_bead_work(:) = helium%worm_xtra_bead(:)
3003 30 : helium%worm_atom_idx_work = helium%worm_atom_idx
3004 30 : IF (startbead > endbead) THEN
3005 0 : DO kbead = 1, endbead
3006 0 : helium%work(:, fendatom, kbead) = helium%pos(:, fendatom, kbead)
3007 : END DO
3008 : END IF
3009 :
3010 : ! reset permtable
3011 30 : tmpi = helium%permutation(bstartatom)
3012 30 : helium%permutation(bstartatom) = helium%permutation(fstartatom)
3013 30 : helium%permutation(fstartatom) = tmpi
3014 : !update inverse permtable
3015 990 : DO ib = 1, SIZE(helium%permutation)
3016 990 : helium%iperm(helium%permutation(ib)) = ib
3017 : END DO
3018 30 : ac = 0
3019 30 : RETURN
3020 : END IF
3021 : END IF
3022 :
3023 : ! finally accepted
3024 8 : ac = 1
3025 102 : DO kbead = startbead, helium%beads
3026 384 : helium%pos(:, bstartatom, kbead) = helium%work(:, bstartatom, kbead)
3027 : END DO
3028 102 : DO kbead = startbead, helium%beads
3029 384 : helium%pos(:, fstartatom, kbead) = helium%work(:, fstartatom, kbead)
3030 : END DO
3031 32 : helium%worm_xtra_bead(:) = helium%worm_xtra_bead_work(:)
3032 8 : helium%worm_atom_idx = helium%worm_atom_idx_work
3033 8 : IF (startbead > endbead) THEN
3034 0 : DO kbead = 1, endbead
3035 0 : 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 32387 : 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 32387 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: work2, work3
3062 32387 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: work
3063 :
3064 226709 : ALLOCATE (work(3, helium%beads + 1), work2(helium%beads + 1), work3(SIZE(helium%uoffdiag, 1) + 1))
3065 32387 : partaction = 0.0_dp
3066 : ! do action in two ways, depending
3067 : ! if coordinates are not wrapping
3068 32387 : IF (startbead < endbead) THEN
3069 : ! helium pair action
3070 : ! every atom, with the one (or two) who got a resampling
3071 797016 : DO iatom = 1, helium%atoms
3072 : ! avoid self interaction
3073 772864 : 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 4623619 : DO ibead = startbead, endbead
3077 16248340 : work(:, ibead + 1 - startbead) = pos(:, iatom, ibead) - pos(:, startatom, ibead)
3078 : END DO
3079 797016 : 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 271755 : DO iatom = 1, helium%atoms
3085 : ! avoid self interaction
3086 263520 : 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 978360 : DO ibead = startbead, helium%beads
3090 3147585 : work(:, ibead + 1 - startbead) = pos(:, iatom, ibead) - pos(:, startatom, ibead)
3091 : END DO
3092 : ! wrapping bond
3093 1021140 : work(:, helium%beads + 2 - startbead) = pos(:, helium%permutation(iatom), 1) - pos(:, endatom, 1)
3094 271755 : partaction = partaction + helium_eval_chain(helium, work, helium%beads - startbead + 2, work2, work3)
3095 : END DO
3096 :
3097 : ! ending atom
3098 271755 : DO iatom = 1, helium%atoms
3099 : ! avoid self interaction
3100 263520 : IF (iatom == endatom) CYCLE
3101 : !from first to endbead
3102 263520 : IF (endbead > 1) THEN
3103 843045 : DO ibead = 1, endbead
3104 2785629 : work(:, ibead) = pos(:, iatom, ibead) - pos(:, endatom, ibead)
3105 : END DO
3106 195517 : partaction = partaction + helium_eval_chain(helium, work, endbead, work2, work3)
3107 : END IF
3108 : END DO
3109 :
3110 : END IF
3111 32387 : DEALLOCATE (work, work2, work3)
3112 :
3113 32387 : 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 7457 : 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 7457 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: work2, work3
3141 7457 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: work
3142 : REAL(KIND=dp), DIMENSION(3) :: r, rp
3143 :
3144 52199 : ALLOCATE (work(3, helium%beads + 1), work2(helium%beads + 1), work3(SIZE(helium%uoffdiag, 1) + 1))
3145 7457 : partaction = 0.0_dp
3146 : ! do action in two ways, depending
3147 : ! if coordinates are not wrapping
3148 7457 : IF (startbead < endbead) THEN
3149 : ! helium pair action
3150 : ! every atom, with the one (or two) who got a resampling
3151 190476 : DO iatom = 1, helium%atoms
3152 : ! avoid self interaction
3153 184704 : 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 178932 : IF (worm_bead_idx - startbead > 1) THEN
3157 781355 : DO ibead = startbead, worm_bead_idx - 1
3158 2597180 : work(:, ibead + 1 - startbead) = pos(:, iatom, ibead) - pos(:, startatom, ibead)
3159 : END DO
3160 176080 : partaction = partaction + helium_eval_chain(helium, work, worm_bead_idx - startbead, work2, work3)
3161 : END IF
3162 :
3163 184704 : IF (endbead > worm_bead_idx) THEN
3164 39184 : DO ibead = worm_bead_idx, endbead
3165 129580 : work(:, ibead + 1 - worm_bead_idx) = pos(:, iatom, ibead) - pos(:, startatom, ibead)
3166 : END DO
3167 9052 : partaction = partaction + helium_eval_chain(helium, work, endbead - worm_bead_idx + 1, work2, work3)
3168 : END IF
3169 : END DO
3170 :
3171 5772 : IF (worm_atom_idx /= startatom) THEN
3172 12144 : DO iatom = 1, helium%atoms
3173 11776 : IF (iatom == startatom) CYCLE
3174 11408 : IF (iatom == worm_atom_idx) CYCLE
3175 44160 : r(:) = pos(:, iatom, worm_bead_idx - 1) - pos(:, startatom, worm_bead_idx - 1)
3176 44160 : rp(:) = pos(:, iatom, worm_bead_idx) - pos(:, startatom, worm_bead_idx)
3177 12144 : partaction = partaction + helium_eval_expansion(helium, r, rp, work3)
3178 : END DO
3179 : ! add pair action with worm
3180 1472 : r(:) = pos(:, startatom, worm_bead_idx - 1) - pos(:, worm_atom_idx, worm_bead_idx - 1)
3181 1472 : rp(:) = pos(:, startatom, worm_bead_idx) - xtrapos(:)
3182 368 : partaction = partaction + helium_eval_expansion(helium, r, rp, work3)
3183 : ELSE
3184 : ! worm intersection
3185 178332 : DO iatom = 1, helium%atoms
3186 172928 : IF (iatom == startatom) CYCLE
3187 670096 : r(:) = pos(:, iatom, worm_bead_idx - 1) - pos(:, startatom, worm_bead_idx - 1)
3188 670096 : rp(:) = pos(:, iatom, worm_bead_idx) - xtrapos(:)
3189 178332 : 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 1685 : IF (worm_bead_idx > startbead) THEN
3196 : ! action terms up to worm gap on starting atom
3197 1320 : DO iatom = 1, helium%atoms
3198 : ! avoid self interaction
3199 1280 : 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 1240 : IF (worm_bead_idx - startbead > 1) THEN
3203 2108 : DO ibead = startbead, worm_bead_idx - 1
3204 6572 : work(:, ibead + 1 - startbead) = pos(:, iatom, ibead) - pos(:, startatom, ibead)
3205 : END DO
3206 620 : partaction = partaction + helium_eval_chain(helium, work, worm_bead_idx - startbead, work2, work3)
3207 : END IF
3208 :
3209 : ! up to the wrapping border
3210 3348 : DO ibead = worm_bead_idx, helium%beads
3211 9672 : 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 4960 : pos(:, helium%permutation(startatom), 1)
3216 1320 : partaction = partaction + helium_eval_chain(helium, work, helium%beads - worm_bead_idx + 2, work2, work3)
3217 :
3218 : END DO
3219 :
3220 : ! ending atom
3221 1320 : DO iatom = 1, helium%atoms
3222 : ! avoid self interaction
3223 1280 : IF (iatom == endatom) CYCLE
3224 : !from first to endbead
3225 1280 : IF (endbead > 1) THEN
3226 2480 : DO ibead = 1, endbead
3227 7688 : work(:, ibead) = pos(:, iatom, ibead) - pos(:, endatom, ibead)
3228 : END DO
3229 744 : partaction = partaction + helium_eval_chain(helium, work, endbead, work2, work3)
3230 : END IF
3231 : END DO
3232 :
3233 40 : IF (worm_atom_idx /= startatom) THEN
3234 1320 : DO iatom = 1, helium%atoms
3235 1280 : IF (iatom == startatom) CYCLE
3236 1240 : IF (iatom == worm_atom_idx) CYCLE
3237 4800 : r(:) = pos(:, iatom, worm_bead_idx - 1) - pos(:, startatom, worm_bead_idx - 1)
3238 4800 : rp(:) = pos(:, iatom, worm_bead_idx) - pos(:, startatom, worm_bead_idx)
3239 1320 : partaction = partaction + helium_eval_expansion(helium, r, rp, work3)
3240 : END DO
3241 : ! add pair action with worm
3242 160 : r(:) = pos(:, startatom, worm_bead_idx - 1) - pos(:, worm_atom_idx, worm_bead_idx - 1)
3243 160 : rp(:) = pos(:, startatom, worm_bead_idx) - xtrapos(:)
3244 40 : partaction = partaction + helium_eval_expansion(helium, r, rp, work3)
3245 : ELSE
3246 : ! worm intersection
3247 0 : DO iatom = 1, helium%atoms
3248 0 : IF (iatom == startatom) CYCLE
3249 0 : r(:) = pos(:, iatom, worm_bead_idx - 1) - pos(:, startatom, worm_bead_idx - 1)
3250 0 : rp(:) = pos(:, iatom, worm_bead_idx) - xtrapos(:)
3251 0 : partaction = partaction + helium_eval_expansion(helium, r, rp, work3)
3252 : END DO
3253 : END IF
3254 : ELSE !(worm_bead_idx < endbead)
3255 1645 : IF (worm_bead_idx /= 1) THEN
3256 36993 : DO iatom = 1, helium%atoms
3257 : ! avoid self interaction
3258 35872 : 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 104780 : DO ibead = startbead, helium%beads
3262 314867 : work(:, ibead + 1 - startbead) = pos(:, iatom, ibead) - pos(:, startatom, ibead)
3263 : END DO
3264 : ! wrapping bond
3265 139004 : work(:, helium%beads - startbead + 2) = pos(:, helium%permutation(iatom), 1) - pos(:, helium%permutation(startatom), 1)
3266 36993 : partaction = partaction + helium_eval_chain(helium, work, helium%beads - startbead + 2, work2, work3)
3267 : END DO
3268 :
3269 : ! ending atom
3270 36993 : DO iatom = 1, helium%atoms
3271 : ! avoid self interaction
3272 35872 : IF (iatom == endatom) CYCLE
3273 : !from first to two before the worm gap
3274 34751 : IF (worm_bead_idx > 2) THEN
3275 73284 : DO ibead = 1, worm_bead_idx - 1
3276 232221 : work(:, ibead) = pos(:, iatom, ibead) - pos(:, endatom, ibead)
3277 : END DO
3278 20305 : partaction = partaction + helium_eval_chain(helium, work, worm_bead_idx - 1, work2, work3)
3279 : END IF
3280 :
3281 35872 : IF (endbead > worm_bead_idx) THEN
3282 5332 : DO ibead = worm_bead_idx, endbead
3283 16864 : work(:, ibead - worm_bead_idx + 1) = pos(:, iatom, ibead) - pos(:, endatom, ibead)
3284 : END DO
3285 1488 : partaction = partaction + helium_eval_chain(helium, work, endbead - worm_bead_idx + 1, work2, work3)
3286 : END IF
3287 : END DO
3288 :
3289 1121 : IF (worm_atom_idx /= endatom) THEN
3290 1848 : DO iatom = 1, helium%atoms
3291 1792 : IF (iatom == endatom) CYCLE
3292 1736 : IF (iatom == worm_atom_idx) CYCLE
3293 6720 : r(:) = pos(:, iatom, worm_bead_idx - 1) - pos(:, endatom, worm_bead_idx - 1)
3294 6720 : rp(:) = pos(:, iatom, worm_bead_idx) - pos(:, endatom, worm_bead_idx)
3295 1848 : partaction = partaction + helium_eval_expansion(helium, r, rp, work3)
3296 : END DO
3297 : ! add pair action with worm
3298 224 : r(:) = pos(:, endatom, worm_bead_idx - 1) - pos(:, worm_atom_idx, worm_bead_idx - 1)
3299 224 : rp(:) = pos(:, endatom, worm_bead_idx) - xtrapos(:)
3300 56 : partaction = partaction + helium_eval_expansion(helium, r, rp, work3)
3301 : ELSE
3302 : ! worm intersection
3303 35145 : DO iatom = 1, helium%atoms
3304 34080 : IF (iatom == endatom) CYCLE
3305 132060 : r(:) = pos(:, iatom, worm_bead_idx - 1) - pos(:, endatom, worm_bead_idx - 1)
3306 132060 : rp(:) = pos(:, iatom, worm_bead_idx) - xtrapos(:)
3307 35145 : 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 17292 : DO iatom = 1, helium%atoms
3313 : ! avoid self interaction
3314 16768 : 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 16768 : IF (helium%beads > startbead) THEN
3318 68634 : DO ibead = startbead, helium%beads
3319 227292 : work(:, ibead - startbead + 1) = pos(:, iatom, ibead) - pos(:, startatom, ibead)
3320 : END DO
3321 15748 : partaction = partaction + helium_eval_chain(helium, work, helium%beads - startbead + 1, work2, work3)
3322 : END IF
3323 : END DO
3324 :
3325 : ! ending atom
3326 17292 : DO iatom = 1, helium%atoms
3327 17292 : IF (endbead > 1) THEN
3328 6144 : DO ibead = 1, endbead
3329 20736 : work(:, ibead) = pos(:, iatom, ibead) - pos(:, endatom, ibead)
3330 : END DO
3331 1280 : partaction = partaction + helium_eval_chain(helium, work, endbead, work2, work3)
3332 : END IF
3333 : END DO
3334 :
3335 524 : IF (worm_atom_idx /= endatom) THEN
3336 1584 : DO iatom = 1, helium%atoms
3337 1536 : IF (iatom == endatom) CYCLE
3338 1488 : IF (iatom == worm_atom_idx) CYCLE
3339 5760 : r(:) = pos(:, helium%iperm(iatom), helium%beads) - pos(:, startatom, helium%beads)
3340 5760 : rp(:) = pos(:, iatom, 1) - pos(:, endatom, 1)
3341 1584 : partaction = partaction + helium_eval_expansion(helium, r, rp, work3)
3342 : END DO
3343 : ! add pair action with worm
3344 192 : r(:) = pos(:, startatom, helium%beads) - pos(:, helium%iperm(worm_atom_idx), helium%beads)
3345 192 : rp(:) = pos(:, endatom, 1) - xtrapos(:)
3346 48 : partaction = partaction + helium_eval_expansion(helium, r, rp, work3)
3347 : ELSE
3348 : ! worm intersection
3349 15708 : DO iatom = 1, helium%atoms
3350 15232 : IF (iatom == endatom) CYCLE
3351 59024 : r(:) = pos(:, helium%iperm(iatom), helium%beads) - pos(:, startatom, helium%beads)
3352 59024 : rp(:) = pos(:, iatom, 1) - xtrapos(:)
3353 15708 : 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 7457 : DEALLOCATE (work, work2, work3)
3360 :
3361 7457 : 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 28136 : 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 28136 : 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 28136 : IF (startbead < endbead) THEN
3393 :
3394 : ! interaction is only beadwise due to primitive coupling
3395 : ! startatom and endatom are the same
3396 89486 : DO ibead = startbead + 1, endbead - 1
3397 : CALL helium_bead_solute_e_f(pint_env, helium, &
3398 68646 : startatom, ibead, pos(:, startatom, ibead), energy=energy)
3399 89486 : 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 20942 : DO ibead = startbead + 1, helium%beads
3407 : CALL helium_bead_solute_e_f(pint_env, helium, &
3408 13646 : startatom, ibead, pos(:, startatom, ibead), energy=energy)
3409 20942 : partaction = partaction + energy
3410 : END DO
3411 : ! second atom after imaginary time wrap
3412 20608 : DO ibead = 1, endbead - 1
3413 : CALL helium_bead_solute_e_f(pint_env, helium, &
3414 13312 : endatom, ibead, pos(:, endatom, ibead), energy=energy)
3415 20608 : partaction = partaction + energy
3416 : END DO
3417 : END IF
3418 :
3419 28136 : partaction = partaction*helium%tau
3420 :
3421 28136 : 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 7388 : 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 7388 : partaction = 0.0_dp
3451 :
3452 : ! helium interaction with the solute
3453 : ! if coordinates are not wrapping
3454 7388 : IF (startbead < worm_bead_idx) THEN
3455 19548 : DO ibead = startbead + 1, worm_bead_idx - 1
3456 : CALL helium_bead_solute_e_f(pint_env, helium, &
3457 13846 : startatom, ibead, pos(:, startatom, ibead), energy=energy)
3458 19548 : partaction = partaction + energy
3459 : END DO
3460 : CALL helium_bead_solute_e_f(pint_env, helium, &
3461 5702 : startatom, ibead, xtrapos, energy=energy)
3462 5702 : partaction = partaction + 0.5_dp*energy
3463 : ELSE !(startbead > worm_bead_idx)
3464 4092 : DO ibead = startbead + 1, helium%beads
3465 : CALL helium_bead_solute_e_f(pint_env, helium, &
3466 2406 : startatom, ibead, pos(:, startatom, ibead), energy=energy)
3467 4092 : partaction = partaction + energy
3468 : END DO
3469 3918 : DO ibead = 1, worm_bead_idx - 1
3470 : CALL helium_bead_solute_e_f(pint_env, helium, &
3471 2232 : worm_atom_idx, ibead, pos(:, worm_atom_idx, ibead), energy=energy)
3472 3918 : partaction = partaction + energy
3473 : END DO
3474 : CALL helium_bead_solute_e_f(pint_env, helium, &
3475 1686 : worm_atom_idx, ibead, xtrapos, energy=energy)
3476 1686 : partaction = partaction + 0.5_dp*energy
3477 : END IF
3478 :
3479 7388 : partaction = partaction*helium%tau
3480 :
3481 7388 : 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 3264 : 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 3264 : partaction = 0.0_dp
3509 :
3510 : ! helium interaction with the solute
3511 : ! if coordinates are not wrapping
3512 3264 : IF (worm_bead_idx < endbead) THEN
3513 : CALL helium_bead_solute_e_f(pint_env, helium, &
3514 2642 : worm_atom_idx, worm_bead_idx, pos(:, worm_atom_idx, worm_bead_idx), energy=energy)
3515 2642 : partaction = partaction + 0.5_dp*energy
3516 9018 : DO ibead = worm_bead_idx + 1, endbead - 1
3517 : CALL helium_bead_solute_e_f(pint_env, helium, &
3518 6376 : endatom, ibead, pos(:, endatom, ibead), energy=energy)
3519 9018 : partaction = partaction + energy
3520 : END DO
3521 : ELSE !(worm_bead_idx > endbead)
3522 : CALL helium_bead_solute_e_f(pint_env, helium, &
3523 622 : worm_atom_idx, worm_bead_idx, pos(:, worm_atom_idx, worm_bead_idx), energy=energy)
3524 622 : partaction = partaction + 0.5_dp*energy
3525 1602 : DO ibead = worm_bead_idx + 1, helium%beads
3526 : CALL helium_bead_solute_e_f(pint_env, helium, &
3527 980 : worm_atom_idx, ibead, pos(:, worm_atom_idx, ibead), energy=energy)
3528 1602 : partaction = partaction + energy
3529 : END DO
3530 1418 : DO ibead = 1, endbead - 1
3531 : CALL helium_bead_solute_e_f(pint_env, helium, &
3532 796 : endatom, ibead, pos(:, endatom, ibead), energy=energy)
3533 1418 : partaction = partaction + energy
3534 : END DO
3535 : END IF
3536 :
3537 3264 : partaction = partaction*helium%tau
3538 :
3539 3264 : END FUNCTION worm_path_inter_action_tail
3540 :
3541 : END MODULE helium_worm
|
