Skip to content

GitLab

Commit 35f0c21c authored by Michael Sherman's avatar Michael Sherman
Browse files

Reorganized the HelloWaterBox example to look like a clone of the HelloEthane example.

parent 2763eed1
Hide whitespace changes
Inline Side-by-side
Showing with 310 additions and 189 deletions
examples/HelloEthane.cpp
View file @ 35f0c21c
...@@ -16,7 +16,6 @@ ...@@ -16,7 +16,6 @@
#include <cstdio> #include <cstdio>
#include <string> #include <string>
#include <vector> #include <vector>
#include <utility>
// ----------------------------------------------------------------------------- // -----------------------------------------------------------------------------
// MOCK MD CODE // MOCK MD CODE
...@@ -40,8 +39,6 @@ const double SimulationTimeInPs = 100; // total simulation time (ps) ...@@ -40,8 +39,6 @@ const double SimulationTimeInPs = 100; // total simulation time (ps)
static const bool WantEnergy = true; static const bool WantEnergy = true;
// FORCE FIELD DATA // FORCE FIELD DATA
// These data structures are not part of OpenMM; they are a model of the kinds
// of data structures an MD code uses to hold a set of force field parameters.
// For this example we're using a tiny subset of the Amber99 force field. // For this example we're using a tiny subset of the Amber99 force field.
// We want to keep the data in the original unit system to avoid conversion // We want to keep the data in the original unit system to avoid conversion
// bugs; this requires conversion on the way in and out of OpenMM. // bugs; this requires conversion on the way in and out of OpenMM.
...@@ -220,12 +217,16 @@ namespace OpenMM { ...@@ -220,12 +217,16 @@ namespace OpenMM {
static const double SigmaPerVdwRadius = 1.78179743628068; static const double SigmaPerVdwRadius = 1.78179743628068;
} }
// This is our opaque "handle" class containing all the OpenMM objects that
// must persist from call to call during a simulation. The main program gets
// a pointer to one of these but sees it as essentially a void* since it
// doesn't know the definition of this class.
struct MyOpenMMData { struct MyOpenMMData {
MyOpenMMData() : system(0), context(0), integrator(0) {} MyOpenMMData() : system(0), context(0), integrator(0) {}
~MyOpenMMData() {delete system; delete context; delete integrator;} ~MyOpenMMData() {delete context; delete integrator; delete system;}
OpenMM::System* system; OpenMM::System* system;
OpenMM::OpenMMContext* context;
OpenMM::Integrator* integrator; OpenMM::Integrator* integrator;
OpenMM::OpenMMContext* context;
}; };
......
examples/HelloWaterBox.cpp
View file @ 35f0c21c
/* ----------------------------------------------------------------------------- /* -----------------------------------------------------------------------------
* OpenMM(tm) HelloWaterBox example (May 2009) * OpenMM(tm) HelloWaterBox example in C++ (June 2009)
* ----------------------------------------------------------------------------- * -----------------------------------------------------------------------------
* This is a complete, self-contained "hello world" example demonstrating * This is a complete, self-contained "hello world" example demonstrating
* GPU-accelerated simulation of a system with both bonded and nonbonded forces, * GPU-accelerated simulation of a system with both bonded and nonbonded forces,
* using water (H-O-H) in a periodic box as an example. A multi-frame PDB file is * using water (H-O-H) in a periodic box as an example. This is a constant-
* written to stdout which can be read by VMD or other visualization tool to * temperature simulation using an Andersen thermostat. A multi-frame PDB file
* is written to stdout which can be read by VMD or other visualization tool to
* produce an animation of the resulting trajectory. * produce an animation of the resulting trajectory.
* *
* Pay particular attention to the handling of units in this example. Incorrect * Pay particular attention to the handling of units in this example. Incorrect
...@@ -13,44 +14,139 @@ ...@@ -13,44 +14,139 @@
* communicating with OpenMM in nanometers and kJoules. * communicating with OpenMM in nanometers and kJoules.
* -------------------------------------------------------------------------- */ * -------------------------------------------------------------------------- */
// Suppress irrelevant warnings from Microsoft's compiler.
#ifdef _MSC_VER
#pragma warning(disable:4996) // sprintf is unsafe
#pragma warning(disable:4251) // no dll interface for some classes
#endif
#include "OpenMM.h"
#include <cstdio> #include <cstdio>
#include <cstdlib>
#include <string> #include <string>
#include <vector> #include <vector>
#include <utility>
using OpenMM::Vec3; // so we can just say "Vec3" below
// ----------------------------------------------------------------------------- // -----------------------------------------------------------------------------
// MODELING AND SIMULATION PARAMETERS // MOCK MD CODE
// ----------------------------------------------------------------------------- // -----------------------------------------------------------------------------
// The code starting here and through main() below is meant to represent in
// simplified form some pre-existing molecular dynamics code, which defines its
// own data structures for force fields, the atoms in this simulation, and the
// simulation parameters, and takes care of recording the trajectory. All this
// has nothing to do with OpenMM; the OpenMM-dependent code comes later and is
// clearly marked below.
// -----------------------------------------------------------------------------
// MODELING AND SIMULATION PARAMETERS
const int NumWatersAlongEdge = 10; // Size of box is NxNxN waters.
const double Temperature = 300; // Kelvins
const double FrictionInPerPs = 91.; // collisions per picosecond
const double CutoffDistanceInAng = 10.; // Angstroms
const bool UseConstraints = true; // Should we constrain O-H bonds? const bool UseConstraints = true; // Should we constrain O-H bonds?
const double StepSizeInFs = 2; // integration step size (fs) const double StepSizeInFs = 2; // integration step size (fs)
const double ReportIntervalInFs = 100; // how often to generate PDB frame (fs) const double ReportIntervalInFs = 100; // how often to generate PDB frame (fs)
const double SimulationTimeInPs = 10; // total simulation time (ps) const double SimulationTimeInPs = 10; // total simulation time (ps)
// FORCE FIELD DATA
// For this example we're using a tiny subset of the Amber99 force field.
// We want to keep the data in the original unit system to avoid conversion
// bugs; this requires conversion on the way in and out of OpenMM.
// Amber reduces nonbonded forces between 1-4 bonded atoms. (These won't be
// used in our all-water simulation.)
const double Coulomb14Scale = 0.5;
const double LennardJones14Scale = 0.5;
// We only need force field parameters for water here.
const double O_mass = 15.9994; // Daltons
const double O_charge = -0.834; // e
const double O_vdwRadInAng = 1.7683; // Angstroms
const double O_vdwEnergyInKcal = 0.1520; // kcal per mole
const double H_mass = 1.00794;
const double H_charge = 0.417;
const double H_vdwRadInAng = 0.0001;
const double H_vdwEnergyInKcal = 0.0000;
// Parameters for the O-H bonds.
const double OH_nominalLengthInAng = 0.9572;
const double OH_stiffnessInKcalPerAng2 = 553.0; // that is, e=k(x-x0)^2
// Parameters for the H-O-H angle.
const double HOH_nominalAngleInDeg = 104.52;
const double HOH_stiffnessInKcalPerRad2 = 100.; // that is e=k(a-a0)^2
// PDB FILE WRITER
// This is a PDB writer that only knows how to write out water molecules. It is
// just here for this example and has nothing to do with OpenMM!
static void
myWritePDBFrame(int frameNum, double timeInPs, const std::vector<double>& atomPosInAng)
{
const char* atomNames[] = {" O ", " H1 ", " H2 "}; // cycle through these
printf("MODEL %d\n", frameNum);
printf("REMARK 250 time=%.3f picoseconds\n", timeInPs);
for (int atom=0; atom < (int)atomPosInAng.size()/3; ++atom)
{
printf("HETATM%5d %4s HOH %4d ", // start of pdb HETATM line
atom+1, atomNames[atom%3], 1 + atom/3); // atom number, name, residue #
printf("%8.3f%8.3f%8.3f", // middle of pdb HETATM line
atomPosInAng[3*atom+0], atomPosInAng[3*atom+1], atomPosInAng[3*atom+2]);
printf(" 1.00 0.00 \n"); // end of pdb HETATM line
}
printf("ENDMDL\n"); // end of trajectory frame
}
// -----------------------------------------------------------------------------
// INTERFACE TO OpenMM
// -----------------------------------------------------------------------------
// These four functions and an opaque structure are used to interface our main
// program with OpenMM without the main program having any direct interaction
// with the OpenMM API. This is a clean approach for interfacing with any MD
// code, although the details of the interface routines will differ. This is
// still just "locally written" code and is not required by OpenMM. Normally
// these would be in another compilation unit but they are defined later in
// this file.
struct MyOpenMMData;
static MyOpenMMData* myInitializeOpenMM(int numWatersAlongEdge,
double temperature,
double frictionInPerPs,
double stepSizeInFs,
std::string& platformName);
static void myStepWithOpenMM(MyOpenMMData*, int numSteps);
static void myGetOpenMMState(MyOpenMMData*, double& time,
std::vector<double>& atomPosInAng);
static void myTerminateOpenMM(MyOpenMMData*);
static void simulateWaterBox();
static void writePDB(const OpenMM::OpenMMContext&); // PDB file writer; see below.
// ----------------------------------------------------------------------------- // -----------------------------------------------------------------------------
// MAIN PROGRAM // WATER BOX MAIN PROGRAM
// ----------------------------------------------------------------------------- // -----------------------------------------------------------------------------
int main() { int main() {
// ALWAYS enclose all OpenMM calls with a try/catch block to make sure that // ALWAYS enclose all OpenMM calls with a try/catch block to make sure that
// usage and runtime errors are caught and reported. // usage and runtime errors are caught and reported.
try { try {
// Load all available OpenMM plugins from their default location. std::string platformName;
OpenMM::Platform::loadPluginsFromDirectory
(OpenMM::Platform::getDefaultPluginsDirectory()); // Set up OpenMM data structures; return handle and OpenMM Platform name.
simulateWaterBox(); MyOpenMMData* omm = myInitializeOpenMM(NumWatersAlongEdge, Temperature,
FrictionInPerPs, StepSizeInFs,
platformName); // output
// Run the simulation:
// (1) Write the first line of the PDB file and the initial configuration.
// (2) Run silently entirely within OpenMM between reporting intervals.
// (3) Write a PDB frame when the time comes.
printf("REMARK Using OpenMM platform %s\n", platformName.c_str());
std::vector<double> atomPositionsInAng; // x,y,z,x,y,z, ...
const int NumSilentSteps = (int)(ReportIntervalInFs / StepSizeInFs + 0.5);
for (int frame=1; ; ++frame) {
double time;
myGetOpenMMState(omm, time, atomPositionsInAng);
myWritePDBFrame(frame, time, atomPositionsInAng);
if (time >= SimulationTimeInPs)
break;
myStepWithOpenMM(omm, NumSilentSteps);
}
// Clean up OpenMM data structures.
myTerminateOpenMM(omm);
return 0; // Normal return from main. return 0; // Normal return from main.
} }
...@@ -64,33 +160,76 @@ int main() { ...@@ -64,33 +160,76 @@ int main() {
// ----------------------------------------------------------------------------- // -----------------------------------------------------------------------------
// FORCE FIELD DATA // OpenMM-USING CODE
// -----------------------------------------------------------------------------
// The OpenMM API is visible only at this point and below. Normally this would
// be in a separate compilation module; we're including it here for simplicity.
// ----------------------------------------------------------------------------- // -----------------------------------------------------------------------------
// These data structures are not part of OpenMM; they are a model of the kinds
// of data structures an MD code uses to hold a set of force field parameters.
// For this example we're using a tiny subset of the Amber99 force field.
// We want to keep the data in the original unit system to avoid conversion
// bugs; this requires conversion on the way in and out of OpenMM.
// Amber reduces nonbonded forces between 1-4 bonded atoms. // Suppress irrelevant warnings from Microsoft's compiler.
const double Coulomb14Scale = 0.5; #ifdef _MSC_VER
const double LennardJones14Scale = 0.5; #pragma warning(disable:4996) // sprintf is unsafe
#pragma warning(disable:4251) // no dll interface for some classes
#endif
#include "OpenMM.h"
using OpenMM::Vec3; // so we can just say "Vec3" below
// Add definition missing from the current version of OpenMM.
namespace OpenMM {
// This conversion factor takes you from a van der Waals radius (defined as 1/2
// the minimum energy separation) to the related Lennard Jones "sigma" parameter
// (defined as the zero crossing separation). The value is 2*pow(2, -1/6).
static const double SigmaPerVdwRadius = 1.78179743628068;
}
// This is our opaque "handle" class containing all the OpenMM objects that
// must persist from call to call during a simulation. The main program gets
// a pointer to one of these but sees it as essentially a void* since it
// doesn't know the definition of this class.
struct MyOpenMMData {
MyOpenMMData() : system(0), context(0), integrator(0) {}
~MyOpenMMData() {delete context; delete integrator; delete system;}
OpenMM::System* system;
OpenMM::Integrator* integrator;
OpenMM::OpenMMContext* context;
};
// This is the conversion factor that takes you from a van der Waals radius
// (defined as 1/2 the minimum energy separation) to the related Lennard Jones
// "sigma" parameter (defined as the zero crossing separation).
static const double SigmaPerVdwRadius = 2*std::pow(2., -1./6.);
// ----------------------------------------------------------------------------- // -----------------------------------------------------------------------------
// WATER SIMULATION // INITIALIZE OpenMM DATA STRUCTURES
// ----------------------------------------------------------------------------- // -----------------------------------------------------------------------------
static void simulateWaterBox() { // We take these actions here:
// ------------------------------------------------------------------------- // (1) Load any available OpenMM plugins, e.g. Cuda and Brook.
// (2) Allocate a MyOpenMMData structure to hang on to OpenMM data structures
// in a manner which is opaque to the caller.
// (3) Fill the OpenMM::System with the force field parameters we want to
// use and the particular set of atoms to be simulated.
// (4) Create an Integrator and a Context associating the Integrator with
// the System.
// (5) Select the OpenMM platform to be used.
// (6) Return the MyOpenMMData struct and the name of the Platform in use.
//
// Note that this function must understand the calling MD code's molecule and
// force field data structures so will need to be customized for each MD code.
static MyOpenMMData*
myInitializeOpenMM( int numWatersAlongEdge,
double temperature,
double frictionInPerPs,
double stepSizeInFs,
std::string& platformName)
{
// Load all available OpenMM plugins from their default location.
OpenMM::Platform::loadPluginsFromDirectory
(OpenMM::Platform::getDefaultPluginsDirectory());
// Allocate space to hold OpenMM objects while we're using them.
MyOpenMMData* omm = new MyOpenMMData();
// Create a System and Force objects within the System. Retain a reference // Create a System and Force objects within the System. Retain a reference
// to each force object so we can fill in the forces. Note: OpenMM owns // to each force object so we can fill in the forces. Note: the System owns
// the objects and will take care of deleting them; don't do it yourself! // the force objects and will take care of deleting them; don't do it yourself!
// ------------------------------------------------------------------------- OpenMM::System& system = *(omm->system = new OpenMM::System());
OpenMM::System system;
OpenMM::NonbondedForce& nonbond = *new OpenMM::NonbondedForce(); OpenMM::NonbondedForce& nonbond = *new OpenMM::NonbondedForce();
system.addForce(&nonbond); system.addForce(&nonbond);
...@@ -102,172 +241,153 @@ static void simulateWaterBox() { ...@@ -102,172 +241,153 @@ static void simulateWaterBox() {
system.addForce(&bondBend); system.addForce(&bondBend);
OpenMM::AndersenThermostat& thermostat = *new OpenMM::AndersenThermostat( OpenMM::AndersenThermostat& thermostat = *new OpenMM::AndersenThermostat(
300, // temperature in kelvins temperature, // kelvins
91.0); // collision frequency in 1/picoseconds frictionInPerPs); // collision frequency in 1/picoseconds
system.addForce(&thermostat);
// -------------------------------------------------------------------------
// Specify the atoms and their properties:
// (1) System needs to know the masses.
// (2) NonbondedForce needs charges,van der Waals properties (in MD units!).
// (3) Collect default positions for initializing the simulation later.
// -------------------------------------------------------------------------
// Volume of one water is 30 Angstroms cubed; // Volume of one water is 30 Angstroms cubed;
// Thus length in one dimension is cube-root of 30, // Thus length in one dimension is cube-root of 30,
// or 3.107 Angstroms or 0.3107 nanometers // or 3.107 Angstroms or 0.3107 nanometers
double waterSize = 0.3107; // edge of cube containing one water, in nanometers const double WaterSizeInNm = 0.3107; // edge of cube containing one water, in nanometers
// Place water molecules one at a time in an NxNxN rectilinear grid
// Place 1000 water molecules one at a time in a 10x10x10 rectilinear grid const double boxEdgeLengthInNm = WaterSizeInNm * numWatersAlongEdge;
const int watersPerEdge = 10;
double boxEdgeLength = waterSize * watersPerEdge;
// Create periodic box // Create periodic box
nonbond.setNonbondedMethod(OpenMM::NonbondedForce::CutoffPeriodic); nonbond.setNonbondedMethod(OpenMM::NonbondedForce::CutoffPeriodic);
nonbond.setCutoffDistance(2); nonbond.setCutoffDistance(CutoffDistanceInAng * OpenMM::NmPerAngstrom);
nonbond.setPeriodicBoxVectors(Vec3(boxEdgeLength,0,0), Vec3(0,boxEdgeLength,0), Vec3(0,0,boxEdgeLength)); nonbond.setPeriodicBoxVectors(Vec3(boxEdgeLengthInNm,0,0),
Vec3(0,boxEdgeLengthInNm,0),
Vec3(0,0,boxEdgeLengthInNm));
// Specify the atoms and their properties:
// (1) System needs to know the masses and constraints (if any).
// (2) NonbondedForce needs charges,van der Waals properties (in MD units!).
// (3) Collect starting positions for initializing the simulation later.
// Create data structures to hold lists of initial positions and bonds // Create data structures to hold lists of initial positions and bonds
std::vector<Vec3> initialPositions; std::vector<Vec3> initialPosInNm;
std::vector< std::pair<int,int> > bondPairs; std::vector< std::pair<int,int> > bondPairs;
// Add water molecules one at a time in the 10x10x10 cubic lattice // Add water molecules one at a time in the 10x10x10 cubic lattice
for (int latticeX = 0; latticeX < watersPerEdge; ++latticeX) for (int latticeX = 0; latticeX < numWatersAlongEdge; ++latticeX)
for (int latticeY = 0; latticeY < watersPerEdge; ++latticeY) for (int latticeY = 0; latticeY < numWatersAlongEdge; ++latticeY)
for (int latticeZ = 0; latticeZ < watersPerEdge; ++latticeZ) for (int latticeZ = 0; latticeZ < numWatersAlongEdge; ++latticeZ)
{ {
// Add parameters for one water molecule // Add parameters for one water molecule
// Add atom masses to system // Add atom masses to system
int oIndex = system.addParticle(15.9994); // O int oIndex = system.addParticle(O_mass); // O
int h1Index = system.addParticle(1.00794); // H1 int h1Index = system.addParticle(H_mass); // H1
int h2Index = system.addParticle(1.00794); // H2 int h2Index = system.addParticle(H_mass); // H2
// Add atom charge, sigma, and stiffness to nonbonded force // Add atom charge, sigma, and stiffness to nonbonded force
nonbond.addParticle( // Oxygen nonbond.addParticle( // Oxygen
-0.834, // charge O_charge,
1.7683 * OpenMM::NmPerAngstrom * SigmaPerVdwRadius, O_vdwRadInAng * OpenMM::NmPerAngstrom * OpenMM::SigmaPerVdwRadius,
0.1520 * OpenMM::KJPerKcal); O_vdwEnergyInKcal * OpenMM::KJPerKcal);
nonbond.addParticle( // Hydrogen1 nonbond.addParticle( // Hydrogen1
0.417, // charge H_charge,
0.0001 * OpenMM::NmPerAngstrom * SigmaPerVdwRadius, H_vdwRadInAng * OpenMM::NmPerAngstrom * OpenMM::SigmaPerVdwRadius,
0.0000 * OpenMM::KJPerKcal); H_vdwEnergyInKcal * OpenMM::KJPerKcal);
nonbond.addParticle( // Hydrogen2 nonbond.addParticle( // Hydrogen2
0.417, // charge H_charge,
0.0001 * OpenMM::NmPerAngstrom * SigmaPerVdwRadius, H_vdwRadInAng * OpenMM::NmPerAngstrom * OpenMM::SigmaPerVdwRadius,
0.0000 * OpenMM::KJPerKcal); H_vdwEnergyInKcal * OpenMM::KJPerKcal);
// constrain O-H bond lengths // Constrain O-H bond lengths or use harmonic forces.
if (UseConstraints) { if (UseConstraints) {
system.addConstraint(oIndex, h1Index, system.addConstraint(oIndex, h1Index,
0.9572 * OpenMM::NmPerAngstrom); OH_nominalLengthInAng * OpenMM::NmPerAngstrom);
system.addConstraint(oIndex, h2Index, system.addConstraint(oIndex, h2Index,
0.9572 * OpenMM::NmPerAngstrom); OH_nominalLengthInAng * OpenMM::NmPerAngstrom);
} else { } else {
// Add stretch parameters for two covalent bonds // Add stretch parameters for two covalent bonds
// Note factor of 2 for stiffness below because Amber specifies the constant // Note factor of 2 for stiffness below because Amber specifies the constant
// as it is used in the harmonic energy term kx^2 with force 2kx; OpenMM wants // as it is used in the harmonic energy term kx^2 with force 2kx; OpenMM wants
// it as used in the force term kx, with energy kx^2/2. // it as used in the force term kx, with energy kx^2/2.
bondStretch.addBond(oIndex, h1Index, bondStretch.addBond(oIndex, h1Index,
0.9572 * OpenMM::NmPerAngstrom, OH_nominalLengthInAng * OpenMM::NmPerAngstrom,
553.0 * 2 * OpenMM::KJPerKcal OH_stiffnessInKcalPerAng2 * 2 * OpenMM::KJPerKcal
* OpenMM::AngstromsPerNm * OpenMM::AngstromsPerNm); * OpenMM::AngstromsPerNm * OpenMM::AngstromsPerNm);
bondStretch.addBond(oIndex, h2Index, bondStretch.addBond(oIndex, h2Index,
0.9572 * OpenMM::NmPerAngstrom, OH_nominalLengthInAng * OpenMM::NmPerAngstrom,
553.0 * 2 * OpenMM::KJPerKcal OH_stiffnessInKcalPerAng2 * 2 * OpenMM::KJPerKcal
* OpenMM::AngstromsPerNm * OpenMM::AngstromsPerNm); * OpenMM::AngstromsPerNm * OpenMM::AngstromsPerNm);
} }
// Store bonds for exclusion list // Store bonds for exclusion list
bondPairs.push_back(std::make_pair(oIndex, h1Index)); bondPairs.push_back(std::make_pair(oIndex, h1Index));
bondPairs.push_back(std::make_pair(oIndex, h2Index)); bondPairs.push_back(std::make_pair(oIndex, h2Index));
// Add bond bend parameters for one angle // Add bond bend parameters for one angle.
// See note under bond stretch above regarding the factor of 2 here. // See note under bond stretch above regarding the factor of 2 here.
bondBend.addAngle(h1Index, oIndex, h2Index, bondBend.addAngle(h1Index, oIndex, h2Index,
104.52 * OpenMM::RadiansPerDegree, HOH_nominalAngleInDeg * OpenMM::RadiansPerDegree,
100.00 * 2 * OpenMM::KJPerKcal); HOH_stiffnessInKcalPerRad2 * 2 * OpenMM::KJPerKcal);
// Location of this molecule in the lattice // Location of this molecule in the lattice
Vec3 latticeVec(waterSize * latticeX, Vec3 latticeVec(WaterSizeInNm * latticeX,
waterSize * latticeY, WaterSizeInNm * latticeY,
waterSize * latticeZ); WaterSizeInNm * latticeZ);
// flip half the waters to prevent giant dipole // flip half the waters to prevent giant dipole
int flip = (rand() % 100) > 49 ? 1 : -1; int flip = (rand() % 100) > 49 ? 1 : -1;
// place this water // place this water
initialPositions.push_back(Vec3(0,0,0) + latticeVec); // O initialPosInNm.push_back(Vec3(0,0,0) + latticeVec); // O
initialPositions.push_back(Vec3(0.09572*flip,0,0) + latticeVec); // H1 initialPosInNm.push_back(Vec3(0.09572*flip,0,0) + latticeVec); // H1
initialPositions.push_back(Vec3(-0.02397*flip,0.09267*flip,0) + latticeVec); // H2 initialPosInNm.push_back(Vec3(-0.02397*flip,0.09267*flip,0) + latticeVec); // H2
} }
// Populate nonbonded exclusions // Populate nonbonded exclusions
nonbond.createExceptionsFromBonds(bondPairs, Coulomb14Scale, LennardJones14Scale); nonbond.createExceptionsFromBonds(bondPairs, Coulomb14Scale, LennardJones14Scale);
// -------------------------------------------------------------------------
// Choose an Integrator for advancing time, and a Context connecting the // Choose an Integrator for advancing time, and a Context connecting the
// System with the Integrator for simulation. Let the Context choose the // System with the Integrator for simulation. Let the Context choose the
// best available Platform. Initialize the configuration from the default // best available Platform. Initialize the configuration from the default
// positions we collected above. Initial velocities will be zero. // positions we collected above. Initial velocities will be zero.
// ------------------------------------------------------------------------- omm->integrator = new OpenMM::VerletIntegrator(StepSizeInFs * OpenMM::PsPerFs);
omm->context = new OpenMM::OpenMMContext(*omm->system, *omm->integrator);
OpenMM::VerletIntegrator integrator( omm->context->setPositions(initialPosInNm);
StepSizeInFs * OpenMM::PsPerFs);
OpenMM::OpenMMContext context(system, integrator);
context.setPositions(initialPositions);
// -------------------------------------------------------------------------
// Run the simulation:
// (1) Write the first line of the PDB file and the initial configuration.
// (2) Run silently entirely within OpenMM between reporting intervals.
// (3) Write a PDB frame when the time comes.
// -------------------------------------------------------------------------
printf("REMARK Using OpenMM platform %s\n", context.getPlatform().getName().c_str() );
writePDB(context);
const int NumSilentSteps = (int)(ReportIntervalInFs / StepSizeInFs + 0.5);
do {
integrator.step(NumSilentSteps);
writePDB(context);
} while (context.getTime() < SimulationTimeInPs);
}
platformName = omm->context->getPlatform().getName();
return omm;
}
// ----------------------------------------------------------------------------- // -----------------------------------------------------------------------------
// PDB FILE WRITER // COPY STATE BACK TO CPU FROM OPENMM
// ----------------------------------------------------------------------------- // -----------------------------------------------------------------------------
static void static void
writePDB(const OpenMM::OpenMMContext& context) { myGetOpenMMState(MyOpenMMData* omm, double& timeInPs,
std::vector<double>& atomPositionsInAng)
{
// We don't calculate energy in this example because it would take most of the time // We don't calculate energy in this example because it would take most of the time
const OpenMM::State state = context.getState(OpenMM::State::Positions); const OpenMM::State state = omm->context->getState(OpenMM::State::Positions);
const std::vector<Vec3>& positions = state.getPositions(); timeInPs = state.getTime(); // OpenMM time is in ps already
static int modelFrameNumber = 0; // numbering for MODEL records in pdb output
// Copy OpenMM positions into output array and change units from nm to Angstroms.
modelFrameNumber++; const std::vector<Vec3>& positionsInNm = state.getPositions();
printf("MODEL %d\n", modelFrameNumber); atomPositionsInAng.resize(3*positionsInNm.size());
printf("REMARK 250 time=%.3f picoseconds\n", state.getTime()); for (int i=0; i < (int)positionsInNm.size(); ++i)
int residueNumber = 1; for (int j=0; j < 3; ++j)
char* atomName = " O "; atomPositionsInAng[3*i+j] = positionsInNm[i][j] * OpenMM::AngstromsPerNm;
for (unsigned int atomIndex=0; atomIndex < positions.size(); ++atomIndex) }
{
// cycle through same three atom names
if (0 == atomIndex%3) {
atomName = " O ";
// increment residue number once every three atoms
++residueNumber;
}
else if (1 == atomIndex%3) atomName = " H1 ";
else atomName = " H2 ";
printf("HETATM%5d %4s HOH %4d ", // start of pdb ATOM line
atomIndex, atomName, residueNumber);
printf("%8.3f%8.3f%8.3f", // middle of pdb ATOM line
positions[atomIndex][0] * OpenMM::AngstromsPerNm, // x
positions[atomIndex][1] * OpenMM::AngstromsPerNm, // y
positions[atomIndex][2] * OpenMM::AngstromsPerNm); // z
printf(" 1.00 0.00 \n"); // end of pdb ATOM line
} // -----------------------------------------------------------------------------
printf("ENDMDL\n"); // end of trajectory frame // TAKE MULTIPLE STEPS USING OpenMM
// -----------------------------------------------------------------------------
static void
myStepWithOpenMM(MyOpenMMData* omm, int numSteps) {
omm->integrator->step(numSteps);
} }
// -----------------------------------------------------------------------------
// DEALLOCATE OpenMM OBJECTS
// -----------------------------------------------------------------------------
static void
myTerminateOpenMM(MyOpenMMData* omm) {
delete omm;
}
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