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