diff --git a/src/motion/helium_methods.F b/src/motion/helium_methods.F index 70378db9b3..3279033755 100644 --- a/src/motion/helium_methods.F +++ b/src/motion/helium_methods.F @@ -1066,6 +1066,7 @@ CONTAINS !> \date 2009-11-09 !> \par History !> 2016-07-14 Modified to work with independent helium_env [cschran] +!> 2018-04-30 Usefull initialization for droplets [fuhl] !> \author Lukasz Walewski ! ************************************************************************************************** SUBROUTINE helium_coord_init(helium_env) @@ -1076,28 +1077,60 @@ CONTAINS INTEGER :: ia, ib, ic, k REAL(KIND=dp) :: r1, r2 + REAL(KIND=dp), DIMENSION(3) :: rvek DO k = 1, SIZE(helium_env) - DO ia = 1, helium_env(k)%helium%atoms - DO ic = 1, 3 - r1 = next_random_number(helium_env(k)%helium%rng_stream_uniform) - r1 = r1*helium_env(k)%helium%cell_size + IF (helium_env(k)%helium%periodic) THEN + DO ia = 1, helium_env(k)%helium%atoms + DO ic = 1, 3 + r1 = next_random_number(helium_env(k)%helium%rng_stream_uniform) + r1 = r1*helium_env(k)%helium%cell_size + DO ib = 1, helium_env(k)%helium%beads +!TODO use thermal gaussian instead + r2 = next_random_number(helium_env(k)%helium%rng_stream_uniform) + helium_env(k)%helium%pos(ic, ia, ib) = r1+0.1_dp*r2 + END DO + END DO DO ib = 1, helium_env(k)%helium%beads -!TODO use thermal gaussian instead - r2 = next_random_number(helium_env(k)%helium%rng_stream_uniform) - helium_env(k)%helium%pos(ic, ia, ib) = r1+0.1_dp*r2 + CALL helium_pbc(helium_env(k)%helium, helium_env(k)%helium%pos(:, ia, ib)) END DO END DO - DO ib = 1, helium_env(k)%helium%beads - CALL helium_pbc(helium_env(k)%helium, helium_env(k)%helium%pos(:, ia, ib)) + ELSE + DO ia = 1, helium_env(k)%helium%atoms + r2 = 2.0_dp*helium_env(k)%helium%droplet_radius**2 + DO WHILE (r2 > helium_env(k)%helium%droplet_radius**2) + DO ic = 1, 3 + rvek(ic) = next_random_number(helium_env(k)%helium%rng_stream_uniform) + rvek(ic) = 2.0_dp*rvek(ic)-1.0_dp + rvek(ic) = rvek(ic)*helium_env(k)%helium%droplet_radius + END DO + r2 = DOT_PRODUCT(rvek, rvek) + ! make sure the helium does not collide with the solute + IF (helium_env(k)%helium%solute_present) THEN + IF (r2 <= 2.0_dp) r2 = 2.0_dp*helium_env(k)%helium%droplet_radius**2 + END IF + END DO + DO ib = 1, helium_env(k)%helium%beads + DO ic = 1, 3 + r2 = next_random_number(helium_env(k)%helium%rng_stream_uniform) + helium_env(k)%helium%pos(ic, ia, ib) = rvek(ic)+0.1_dp*r2 + END DO + ! Make sure, that nothing lies outside the droplet radius + r1 = DOT_PRODUCT(helium_env(k)%helium%pos(:, ia, ib), & + helium_env(k)%helium%pos(:, ia, ib)) + IF (r1 >= helium_env(k)%helium%droplet_radius**2) THEN + r1 = SQRT(r1) + helium_env(k)%helium%pos(:, ia, ib) = & + helium_env(k)%helium%droplet_radius/r1* & + helium_env(k)%helium%pos(:, ia, ib) + END IF + END DO END DO - END DO - + END IF ! store positions at time slice nbeads+1 (rperez): DO ia = 1, helium_env(k)%helium%atoms helium_env(k)%helium%pos(:, ia, helium_env(k)%helium%beads+1) = helium_env(k)%helium%pos(:, ia, 1) END DO - ! initialize work array (rperez): helium_env(k)%helium%work = helium_env(k)%helium%pos END DO diff --git a/src/motion/helium_sampling.F b/src/motion/helium_sampling.F index c000ef8017..2f57616106 100644 --- a/src/motion/helium_sampling.F +++ b/src/motion/helium_sampling.F @@ -761,12 +761,13 @@ CONTAINS routineP = moduleN//':'//routineN CHARACTER(len=default_string_length) :: err_str, stmp - INTEGER :: c, ifix, j, k, l, level, pk1, pk2, stride + INTEGER :: c, ia, ib, ic, ifix, j, k, l, level, & + pk1, pk2, stride INTEGER, DIMENSION(:), POINTER :: p, perm LOGICAL :: nperiodic, should_reject REAL(KIND=dp) :: cell_size, ds, dtk, e1, e2, pds, & - prev_ds, r, rtmp, sigma, x - REAL(KIND=dp), DIMENSION(3) :: bis, new_com, rm1, rm2, tmp1, tmp2 + prev_ds, r, rtmp, rtmpo, sigma, x + REAL(KIND=dp), DIMENSION(3) :: bis, biso, new_com, rm1, rm2, tmp1, tmp2 REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: pos, work TYPE(spline_data_p_type), DIMENSION(:, :), POINTER :: uij @@ -1158,8 +1159,9 @@ CONTAINS END IF ! accepted. - ! rejection of the whole move if any of the new positions is farther away + ! rejection of the whole move if any centroid is farther away ! from the current center of gravity than HELIUM%DROPLET_RADIUS [lwalewski] + ! AND ist not moving towards the center IF (helium%droplet_radius .LT. HUGE(1.0_dp)) THEN IF (helium%solute_present) THEN new_com(:) = helium%center(:) @@ -1172,18 +1174,33 @@ CONTAINS END DO new_com(:) = new_com(:)/helium%atoms/helium%beads END IF - ! cycle through new positions + ! Cycle through atoms (ignore connectivity) should_reject = .FALSE. - outer: DO c = ifix, helium%beads - DO k = 1, cyclen - bis(:) = work(:, p(k), c)-new_com(:) - rtmp = SQRT(bis(1)*bis(1)+bis(2)*bis(2)+bis(3)*bis(3)) - IF (rtmp .GE. helium%droplet_radius) THEN + outer: DO ia = 1, helium%atoms + bis(:) = 0.0_dp + DO ib = 1, helium%beads + DO ic = 1, 3 + bis(ic) = bis(ic)+work(ic, ia, ib)-new_com(ic) + END DO + END DO + bis(:) = bis(:)/helium%beads + rtmp = DOT_PRODUCT(bis, bis) + IF (rtmp .GE. helium%droplet_radius**2) THEN + biso(:) = 0.0_dp + DO ib = 1, helium%beads + DO ic = 1, 3 + biso(ic) = biso(ic)+pos(ic, ia, ib)-helium%center(ic) + END DO + END DO + biso(:) = biso(:)/helium%beads + rtmpo = DOT_PRODUCT(biso, biso) + ! only reject if it moves away from COM outside the droplet radius + IF (rtmpo < rtmp) THEN ! found - this one does not fit within R from the center should_reject = .TRUE. EXIT outer END IF - END DO + END IF END DO outer IF (should_reject) THEN ! restore work and perm, then return @@ -1200,6 +1217,7 @@ CONTAINS perm(p(1)) = c END IF RETURN + ELSE END IF END IF ! accepted. diff --git a/tests/Pimd/TEST_FILES b/tests/Pimd/TEST_FILES index 9c83b3c477..a76ca07368 100644 --- a/tests/Pimd/TEST_FILES +++ b/tests/Pimd/TEST_FILES @@ -13,7 +13,7 @@ he32_only_restart.inp 40 2e-14 -7.039 water-in-helium.inp 9 7e-14 1.0713913726284391E-003 water-in-helium-restart.inp 9 1e-13 1.0717908679185048E-003 water-in-helium-lastforce.inp 9 7e-14 1.0714865845104339E-003 -water-micro-helium.inp 9 1e-13 1.0713355696395585E-003 +water-micro-helium.inp 9 1e-13 1.0713363424216917E-003 w512_pint_nose.inp 9 8e-14 9.5780711530210407 w512_pint_gle.inp 9 3e-13 -3.5646497002020765 w512_pint_pile.inp 9 3e-12 0.86626194213382846 diff --git a/tests/Pimd/water-micro-helium.inp b/tests/Pimd/water-micro-helium.inp index 1c03fed306..58c45005fb 100644 --- a/tests/Pimd/water-micro-helium.inp +++ b/tests/Pimd/water-micro-helium.inp @@ -21,18 +21,18 @@ PRESAMPLE T SOLUTE_INTERACTION MWATER POTENTIAL_FILE_NAME lj-test-pot.dat - NATOMS 2 + NATOMS 5 NBEADS 16 - NUM_ENV 7 - INOROT 200 - IROT 100 + NUM_ENV 5 + N_INNER 50 + N_OUTER 7 SAMPLING_METHOD CEPERLEY &CEPERLEY BISECTION 8 MAX_PERM_CYCLE 1 &END CEPERLEY PERIODIC F - DROPLET_RADIUS 15 + DROPLET_RADIUS 10 &PRINT &COORDINATES ON &EACH