mirror of
https://github.com/cp2k/cp2k.git
synced 2026-07-28 14:15:19 -04:00
Useful initialization for helium droplet simulations, as well as refinement on metropolis acceptance for out of droplet steps [Felix Uhl]
svn-origin-rev: 18410
This commit is contained in:
Harald Forbert
parent
137198a7c9
commit
0ceed53270
4 changed files with 80 additions and 29 deletions
|
|
@ -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
|
||||
|
||||
|
|
|
|||
|
|
@ -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.
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue