-
Notifications
You must be signed in to change notification settings - Fork 103
/
diffusion_test.py
executable file
·130 lines (111 loc) · 6.63 KB
/
diffusion_test.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
#!/usr/bin/env python3
# diffusion_test.py
#------------------------------------------------------------------------------------------------#
# This software was written in 2016/17 #
# by Michael P. Allen <[email protected]>/<[email protected]> #
# and Dominic J. Tildesley <[email protected]> ("the authors"), #
# to accompany the book "Computer Simulation of Liquids", second edition, 2017 ("the text"), #
# published by Oxford University Press ("the publishers"). #
# #
# LICENCE #
# Creative Commons CC0 Public Domain Dedication. #
# To the extent possible under law, the authors have dedicated all copyright and related #
# and neighboring rights to this software to the PUBLIC domain worldwide. #
# This software is distributed without any warranty. #
# You should have received a copy of the CC0 Public Domain Dedication along with this software. #
# If not, see <http://creativecommons.org/publicdomain/zero/1.0/>. #
# #
# DISCLAIMER #
# The authors and publishers make no warranties about the software, and disclaim liability #
# for all uses of the software, to the fullest extent permitted by applicable law. #
# The authors and publishers do not recommend use of this software for any purpose. #
# It is made freely available, solely to clarify points made in the text. When using or citing #
# the software, you should not imply endorsement by the authors or publishers. #
#------------------------------------------------------------------------------------------------#
"""Generate test data for diffusion.py."""
def a_propagator ( t, r, v ):
"""A propagator (drift).
t is the time over which to propagate (typically dt/2)
r and v are the current positions and velocities
The function returns the new positions.
"""
return r + v * t
def o_propagator ( t, v ):
"""O propagator (friction and random contributions).
t is the time over which to propagate (typically dt)
v is the current velocity
gamma and temperature are accessed from the calling program
The function returns the new velocities.
"""
import numpy as np
x = gamma*t
c = -np.expm1(-2*x) # 1-exp(-2*x), preserving accuracy for small x
c = np.sqrt(c)
return np.exp(-x) * v + c * np.sqrt(temperature) * np.random.randn(n,3)
# diffusion_test program
import json
import sys
import numpy as np
import math
from platform import python_version
from config_io_module import write_cnf_atoms
print('diffusion_test')
print('Python: '+python_version())
print('NumPy: '+np.__version__)
print()
print('Brownian dynamics without interactions, constant-NVT ensemble')
print('Particle mass m=1 throughout')
# Read parameters in JSON format
try:
nml = json.load(sys.stdin)
except json.JSONDecodeError:
print('Exiting on Invalid JSON format')
sys.exit()
# Set default values, check keys and typecheck values
defaults = {"n":250, "nblock":999, "nstep":25, "dt":0.002, "gamma":1.0, "temperature":1.0, "box":1.0}
for key, val in nml.items():
if key in defaults:
assert type(val) == type(defaults[key]), key+" has the wrong type"
else:
print('Warning', key, 'not in ',list(defaults.keys()))
# Set parameters to input values or defaults
n = nml["n"] if "n" in nml else defaults["n"]
nblock = nml["nblock"] if "nblock" in nml else defaults["nblock"]
nstep = nml["nstep"] if "nstep" in nml else defaults["nstep"]
dt = nml["dt"] if "dt" in nml else defaults["dt"]
gamma = nml["gamma"] if "gamma" in nml else defaults["gamma"]
temperature = nml["temperature"] if "temperature" in nml else defaults["temperature"]
box = nml["box"] if "box" in nml else defaults["box"]
# Write out parameters
print( "{:40}{:15d} ".format('Number of atoms', n) )
print( "{:40}{:15d} ".format('Number of blocks', nblock) )
print( "{:40}{:15d} ".format('Number of steps per block', nstep) )
print( "{:40}{:15.6f}".format('Time step', dt) )
print( "{:40}{:15.6f}".format('Friction coefficient', gamma) )
print( "{:40}{:15.6f}".format('Temperature', temperature) )
print( "{:40}{:15.6f}".format('Ideal diffusion coefft', temperature/gamma) )
np.random.seed()
r = np.random.rand(n,3).astype(np.float64) # Random positions
r = r - 0.5 # Now in range (-1/2,1/2)
r = r * box # Now in range (-box/2,box/2)
v = np.random.randn(n,3).astype(np.float64) # Random velocities
v = v * np.sqrt(temperature) # At desired temperature
cnf_prefix = 'cnf.'
sav_tag = '000'
write_cnf_atoms ( cnf_prefix+sav_tag, n, box, r, v )
for blk in range(nblock): # Loop over blocks
for stp in range(nstep): # Loop over steps
r = a_propagator(dt/2.0,r,v) # A drift half-step
v = o_propagator(dt,v) # O random velocities and friction step
r = a_propagator(dt/2.0,r,v) # A drift half-step
r = r - np.rint(r/box)*box # Periodic boundaries
sav_tag = str(blk+1).zfill(3) if blk<999 else 'sav' # Number configuration by block
write_cnf_atoms ( cnf_prefix+sav_tag, n, box, r, v ) # Save configuration
print('Exact results output to diffusion_exact.out')
with open("diffusion_exact.out","w") as f:
for blk in range(nblock//2): # Loop up to half the run length
t = (blk*nstep)*dt # Time advances block by block
vacf = 3.0*temperature * math.exp(-gamma*t) # Velocity autocorrelation function
rvcf = 3.0*temperature * ( 1.0 - math.exp(-gamma*t) ) / gamma # Velocity-displacement correlation
msd = 6.0*temperature * ( t - ( 1.0 - math.exp(-gamma*t) ) / gamma ) / gamma # Mean-square displacement
print("{:15.6f}{:15.8f}{:15.8f}{:15.8f}".format(t, vacf, rvcf, msd), file=f)