Skip to content
GitLab
Menu
Projects
Groups
Snippets
Loading...
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Contribute to GitLab
Sign in
Toggle navigation
Menu
Open sidebar
tsoc
openmm
Commits
260b9598
Commit
260b9598
authored
Aug 19, 2013
by
Robert McGibbon
Browse files
2to3 fixes for argon-chemical-potential
parent
443dceb8
Changes
1
Show whitespace changes
Inline
Side-by-side
Showing
1 changed file
with
20 additions
and
19 deletions
+20
-19
examples/argon-chemical-potential.py
examples/argon-chemical-potential.py
+20
-19
No files found.
examples/argon-chemical-potential.py
View file @
260b9598
...
@@ -15,6 +15,7 @@
...
@@ -15,6 +15,7 @@
# [1] Michael R. Shirts and John D. Chodera. Statistically optimal analysis of samples from multiple equilibrium states.
# [1] Michael R. Shirts and John D. Chodera. Statistically optimal analysis of samples from multiple equilibrium states.
# J. Chem. Phys. 129:124105 (2008) http://dx.doi.org/10.1063/1.2978177
# J. Chem. Phys. 129:124105 (2008) http://dx.doi.org/10.1063/1.2978177
from
__future__
import
print_function
from
simtk.openmm
import
*
from
simtk.openmm
import
*
from
simtk.unit
import
*
from
simtk.unit
import
*
import
numpy
import
numpy
...
@@ -51,15 +52,15 @@ nlambda = len(lambda_values)
...
@@ -51,15 +52,15 @@ nlambda = len(lambda_values)
volume
=
nparticles
*
(
sigma
**
3
)
/
reduced_density
volume
=
nparticles
*
(
sigma
**
3
)
/
reduced_density
box_edge
=
volume
**
(
1.0
/
3.0
)
box_edge
=
volume
**
(
1.0
/
3.0
)
cutoff
=
min
(
box_edge
*
0.49
,
2.5
*
sigma
)
# Compute cutoff
cutoff
=
min
(
box_edge
*
0.49
,
2.5
*
sigma
)
# Compute cutoff
print
"sigma = %s"
%
sigma
print
(
"sigma = %s"
%
sigma
)
print
"box_edge = %s"
%
box_edge
print
(
"box_edge = %s"
%
box_edge
)
print
"cutoff = %s"
%
cutoff
print
(
"cutoff = %s"
%
cutoff
)
# =============================================================================
# =============================================================================
# Build systems at each alchemical lambda value.
# Build systems at each alchemical lambda value.
# =============================================================================
# =============================================================================
print
"Building alchemically-modified systems..."
print
(
"Building alchemically-modified systems..."
)
alchemical_systems
=
list
()
# alchemical_systems[i] is the alchemically-modified System object corresponding to lambda_values[i]
alchemical_systems
=
list
()
# alchemical_systems[i] is the alchemically-modified System object corresponding to lambda_values[i]
for
lambda_value
in
lambda_values
:
for
lambda_value
in
lambda_values
:
# Create argon system where first particle is alchemically modified by lambda_value.
# Create argon system where first particle is alchemically modified by lambda_value.
...
@@ -112,17 +113,17 @@ for (lambda_index, lambda_value) in enumerate(lambda_values):
...
@@ -112,17 +113,17 @@ for (lambda_index, lambda_value) in enumerate(lambda_values):
context
.
setPositions
(
positions
)
context
.
setPositions
(
positions
)
# Minimize energy from coordinates.
# Minimize energy from coordinates.
print
"Lambda %d/%d : minimizing..."
%
(
lambda_index
+
1
,
nlambda
)
print
(
"Lambda %d/%d : minimizing..."
%
(
lambda_index
+
1
,
nlambda
)
)
LocalEnergyMinimizer
.
minimize
(
context
)
LocalEnergyMinimizer
.
minimize
(
context
)
# Equilibrate.
# Equilibrate.
print
"Lambda %d/%d : equilibrating..."
%
(
lambda_index
+
1
,
nlambda
)
print
(
"Lambda %d/%d : equilibrating..."
%
(
lambda_index
+
1
,
nlambda
)
)
integrator
.
step
(
nequil_steps
)
integrator
.
step
(
nequil_steps
)
# Sample.
# Sample.
position_history
=
list
()
# position_history[i] is the set of positions after iteration i
position_history
=
list
()
# position_history[i] is the set of positions after iteration i
for
iteration
in
range
(
nprod_iterations
):
for
iteration
in
range
(
nprod_iterations
):
print
"Lambda %d/%d : production iteration %d/%d"
%
(
lambda_index
+
1
,
nlambda
,
iteration
+
1
,
nprod_iterations
)
print
(
"Lambda %d/%d : production iteration %d/%d"
%
(
lambda_index
+
1
,
nlambda
,
iteration
+
1
,
nprod_iterations
)
)
# Run dynamics.
# Run dynamics.
integrator
.
step
(
nprod_steps
)
integrator
.
step
(
nprod_steps
)
...
@@ -138,7 +139,7 @@ for (lambda_index, lambda_value) in enumerate(lambda_values):
...
@@ -138,7 +139,7 @@ for (lambda_index, lambda_value) in enumerate(lambda_values):
del
context
,
integrator
del
context
,
integrator
# Compute reduced potentials of all snapshots at all alchemical states for MBAR.
# Compute reduced potentials of all snapshots at all alchemical states for MBAR.
print
"Lambda %d/%d : computing energies at all states..."
%
(
lambda_index
+
1
,
nlambda
)
print
(
"Lambda %d/%d : computing energies at all states..."
%
(
lambda_index
+
1
,
nlambda
)
)
beta
=
1.0
/
(
BOLTZMANN_CONSTANT_kB
*
AVOGADRO_CONSTANT_NA
*
temperature
)
# inverse temperature
beta
=
1.0
/
(
BOLTZMANN_CONSTANT_kB
*
AVOGADRO_CONSTANT_NA
*
temperature
)
# inverse temperature
for
l
in
range
(
nlambda
):
for
l
in
range
(
nlambda
):
# Set up Context just to evaluate energies.
# Set up Context just to evaluate energies.
...
@@ -159,13 +160,13 @@ try:
...
@@ -159,13 +160,13 @@ try:
from
timeseries
import
subsampleCorrelatedData
from
timeseries
import
subsampleCorrelatedData
from
pymbar
import
MBAR
from
pymbar
import
MBAR
except
:
except
:
raise
"pymbar [https://simtk.org/home/pymbar] must be installed to complete analysis of free energies."
raise
ImportError
(
"pymbar [https://simtk.org/home/pymbar] must be installed to complete analysis of free energies."
)
# =============================================================================
# =============================================================================
# Subsample correlated samples to generate uncorrelated subsample.
# Subsample correlated samples to generate uncorrelated subsample.
# =============================================================================
# =============================================================================
print
"Subsampling data to remove correlation..."
print
(
"Subsampling data to remove correlation..."
)
K
=
nlambda
# number of states
K
=
nlambda
# number of states
N_k
=
nprod_iterations
*
numpy
.
ones
([
K
],
numpy
.
int32
)
# N_k[k] is the number of uncorrelated samples at state k
N_k
=
nprod_iterations
*
numpy
.
ones
([
K
],
numpy
.
int32
)
# N_k[k] is the number of uncorrelated samples at state k
u_kln_subsampled
=
numpy
.
zeros
([
K
,
K
,
nprod_iterations
],
numpy
.
float64
)
# subsampled data
u_kln_subsampled
=
numpy
.
zeros
([
K
,
K
,
nprod_iterations
],
numpy
.
float64
)
# subsampled data
...
@@ -176,24 +177,24 @@ for k in range(K):
...
@@ -176,24 +177,24 @@ for k in range(K):
N_k
[
k
]
=
len
(
indices
)
N_k
[
k
]
=
len
(
indices
)
for
l
in
range
(
K
):
for
l
in
range
(
K
):
u_kln_subsampled
[
k
,
l
,
0
:
len
(
indices
)]
=
u_kln
[
k
,
l
,
indices
]
u_kln_subsampled
[
k
,
l
,
0
:
len
(
indices
)]
=
u_kln
[
k
,
l
,
indices
]
print
"Number of uncorrelated samples per state:"
print
(
"Number of uncorrelated samples per state:"
)
print
N_k
print
(
N_k
)
# =============================================================================
# =============================================================================
# Analyze with MBAR to compute free energy differences and statistical errors.
# Analyze with MBAR to compute free energy differences and statistical errors.
# =============================================================================
# =============================================================================
print
"Analyzing with MBAR..."
print
(
"Analyzing with MBAR..."
)
mbar
=
MBAR
(
u_kln_subsampled
,
N_k
)
mbar
=
MBAR
(
u_kln_subsampled
,
N_k
)
[
Deltaf_ij
,
dDeltaf_ij
]
=
mbar
.
getFreeEnergyDifferences
()
Deltaf_ij
,
dDeltaf_ij
=
mbar
.
getFreeEnergyDifferences
()
print
"Free energy differences (in kT)"
print
(
"Free energy differences (in kT)"
)
print
Deltaf_ij
print
(
Deltaf_ij
)
print
"Statistical errors (in kT)"
print
(
"Statistical errors (in kT)"
)
print
dDeltaf_ij
print
(
dDeltaf_ij
)
# =============================================================================
# =============================================================================
# Report result.
# Report result.
# =============================================================================
# =============================================================================
print
"Free energy of inserting argon particle: %.3f +- %.3f kT"
%
(
Deltaf_ij
[
0
,
K
-
1
],
dDeltaf_ij
[
0
,
K
-
1
])
print
(
"Free energy of inserting argon particle: %.3f +- %.3f kT"
%
(
Deltaf_ij
[
0
,
K
-
1
],
dDeltaf_ij
[
0
,
K
-
1
])
)
Write
Preview
Markdown
is supported
0%
Try again
or
attach a new file
.
Attach a file
Cancel
You are about to add
0
people
to the discussion. Proceed with caution.
Finish editing this message first!
Cancel
Please
register
or
sign in
to comment