Commit b18c82c8 authored by peastman's avatar peastman
Browse files

Merge pull request #115 from rmcgibbo/py3

Small py3k fixes
parents 7dd882eb 6e89e067
......@@ -15,6 +15,7 @@
# [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
from __future__ import print_function
from simtk.openmm import *
from simtk.unit import *
import numpy
......@@ -51,15 +52,15 @@ nlambda = len(lambda_values)
volume = nparticles*(sigma**3)/reduced_density
box_edge = volume**(1.0/3.0)
cutoff = min(box_edge*0.49, 2.5*sigma) # Compute cutoff
print "sigma = %s" % sigma
print "box_edge = %s" % box_edge
print "cutoff = %s" % cutoff
print("sigma = %s" % sigma)
print("box_edge = %s" % box_edge)
print("cutoff = %s" % cutoff)
# =============================================================================
# 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]
for lambda_value in lambda_values:
# 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):
context.setPositions(positions)
# 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)
# Equilibrate.
print "Lambda %d/%d : equilibrating..." % (lambda_index+1, nlambda)
print("Lambda %d/%d : equilibrating..." % (lambda_index+1, nlambda))
integrator.step(nequil_steps)
# Sample.
position_history = list() # position_history[i] is the set of positions after iteration i
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.
integrator.step(nprod_steps)
......@@ -138,7 +139,7 @@ for (lambda_index, lambda_value) in enumerate(lambda_values):
del context, integrator
# 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
for l in range(nlambda):
# Set up Context just to evaluate energies.
......@@ -159,13 +160,13 @@ try:
from timeseries import subsampleCorrelatedData
from pymbar import MBAR
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.
# =============================================================================
print "Subsampling data to remove correlation..."
print("Subsampling data to remove correlation...")
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
u_kln_subsampled = numpy.zeros([K,K,nprod_iterations], numpy.float64) # subsampled data
......@@ -176,24 +177,24 @@ for k in range(K):
N_k[k] = len(indices)
for l in range(K):
u_kln_subsampled[k,l,0:len(indices)] = u_kln[k,l,indices]
print "Number of uncorrelated samples per state:"
print N_k
print("Number of uncorrelated samples per state:")
print(N_k)
# =============================================================================
# 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)
[Deltaf_ij, dDeltaf_ij] = mbar.getFreeEnergyDifferences()
print "Free energy differences (in kT)"
print Deltaf_ij
print "Statistical errors (in kT)"
print dDeltaf_ij
Deltaf_ij, dDeltaf_ij = mbar.getFreeEnergyDifferences()
print("Free energy differences (in kT)")
print(Deltaf_ij)
print("Statistical errors (in kT)")
print(dDeltaf_ij)
# =============================================================================
# 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]))
from __future__ import print_function
# First make sure OpenMM is installed.
import sys
try:
from simtk.openmm.app import *
from simtk.openmm import *
from simtk.unit import *
except ImportError as err:
print "Failed to import OpenMM packages:", err.message
print "Make sure OpenMM is installed and the library path is set correctly."
print("Failed to import OpenMM packages:", err.message)
print("Make sure OpenMM is installed and the library path is set correctly.")
sys.exit()
# Create a System for the tests.
......@@ -19,28 +21,28 @@ system = forcefield.createSystem(pdb.topology, nonbondedMethod=PME, nonbondedCut
# List all installed platforms and compute forces with each one.
numPlatforms = Platform.getNumPlatforms()
print "There are", numPlatforms, "Platforms available:"
print
print("There are", numPlatforms, "Platforms available:")
print()
forces = [None]*numPlatforms
for i in range(numPlatforms):
platform = Platform.getPlatform(i)
print i+1, platform.getName(),
print(i+1, platform.getName(), end=" ")
integrator = LangevinIntegrator(300*kelvin, 1/picosecond, 0.002*picoseconds)
try:
simulation = Simulation(pdb.topology, system, integrator, platform)
simulation.context.setPositions(pdb.positions)
forces[i] = simulation.context.getState(getForces=True).getForces()
del simulation
print "- Successfully computed forces"
print("- Successfully computed forces")
except:
print "- Error computing forces with", platform.getName(), "platform"
print("- Error computing forces with", platform.getName(), "platform")
# See how well the platforms agree.
if numPlatforms > 1:
print
print "Median difference in forces between platforms:"
print
print()
print("Median difference in forces between platforms:")
print()
for i in range(numPlatforms):
for j in range(i):
if forces[i] is not None and forces[j] is not None:
......@@ -49,4 +51,6 @@ if numPlatforms > 1:
d = f1-f2
error = sqrt((d[0]*d[0]+d[1]*d[1]+d[2]*d[2])/(f1[0]*f1[0]+f1[1]*f1[1]+f1[2]*f1[2]))
errors.append(error)
print "%s vs. %s: %g" % (Platform.getPlatform(j).getName(), Platform.getPlatform(i).getName(), sorted(errors)[len(errors)/2])
print("{} vs. {}: {:g}".format(Platform.getPlatform(j).getName(),
Platform.getPlatform(i).getName(),
sorted(errors)[len(errors)//2]))
......@@ -335,8 +335,8 @@ class PrmtopLoader(object):
% ((bondPointers[ii],
bondPointers[ii+1]),))
iType=int(bondPointers[ii+2])-1
returnList.append((int(bondPointers[ii])/3,
int(bondPointers[ii+1])/3,
returnList.append((int(bondPointers[ii])//3,
int(bondPointers[ii+1])//3,
float(forceConstant[iType])*forceConstConversionFactor,
float(bondEquil[iType])*lengthConversionFactor))
return returnList
......@@ -383,9 +383,9 @@ class PrmtopLoader(object):
anglePointers[ii+1],
anglePointers[ii+2]),))
iType=int(anglePointers[ii+3])-1
self._angleList.append((int(anglePointers[ii])/3,
int(anglePointers[ii+1])/3,
int(anglePointers[ii+2])/3,
self._angleList.append((int(anglePointers[ii])//3,
int(anglePointers[ii+1])//3,
int(anglePointers[ii+2])//3,
float(forceConstant[iType])*forceConstConversionFactor,
float(angleEquil[iType])))
return self._angleList
......@@ -411,10 +411,10 @@ class PrmtopLoader(object):
dihedralPointers[ii+2],
dihedralPointers[ii+3]),))
iType=int(dihedralPointers[ii+4])-1
self._dihedralList.append((int(dihedralPointers[ii])/3,
int(dihedralPointers[ii+1])/3,
abs(int(dihedralPointers[ii+2]))/3,
abs(int(dihedralPointers[ii+3]))/3,
self._dihedralList.append((int(dihedralPointers[ii])//3,
int(dihedralPointers[ii+1])//3,
abs(int(dihedralPointers[ii+2]))//3,
abs(int(dihedralPointers[ii+3]))//3,
float(forceConstant[iType])*forceConstConversionFactor,
float(phase[iType]),
int(0.5+float(periodicity[iType]))))
......@@ -429,8 +429,8 @@ class PrmtopLoader(object):
nonbondTerms = self.getNonbondTerms()
for ii in range(0,len(dihedralPointers),5):
if int(dihedralPointers[ii+2])>0 and int(dihedralPointers[ii+3])>0:
iAtom = int(dihedralPointers[ii])/3
lAtom = int(dihedralPointers[ii+3])/3
iAtom = int(dihedralPointers[ii])//3
lAtom = int(dihedralPointers[ii+3])//3
chargeProd = charges[iAtom]*charges[lAtom]
(rVdwI, epsilonI) = nonbondTerms[iAtom]
(rVdwL, epsilonL) = nonbondTerms[lAtom]
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment