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