Skip to content

GitLab

Commit a8783de3 authored by Michael Sherman's avatar Michael Sherman
Browse files

Reworked the Ethane example to use the same structure as the NaCl one, and... 

Reworked the Ethane example to use the same structure as the NaCl one, and made a few minor changes to NaCl too.
parent 138c1c88
Hide whitespace changes
Inline Side-by-side
Showing with 226 additions and 122 deletions
examples/HelloEthane.cpp
View file @ a8783de3
/* -----------------------------------------------------------------------------
* OpenMM(tm) HelloEthane example (May 2009)
* OpenMM(tm) HelloEthane example in C++ (June 2009)
* -----------------------------------------------------------------------------
* This is a complete, self-contained "hello world" example demonstrating
* GPU-accelerated simulation of a system with both bonded and nonbonded forces,
......@@ -13,59 +13,33 @@
* 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 <string>
#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 bool UseConstraints = false; // Should we constrain C-H bonds?
const double StepSizeInFs = 2; // integration step size (fs)
const double ReportIntervalInFs = 10; // how often to generate PDB frame (fs)
const double SimulationTimeInPs = 100; // total simulation time (ps)
static void simulateEthane();
static void writePDB(const OpenMM::OpenMMContext&); // PDB file writer; see below.
// -----------------------------------------------------------------------------
// MAIN PROGRAM
// -----------------------------------------------------------------------------
int main() {
// ALWAYS enclose all OpenMM calls with a try/catch block to make sure that
// usage and runtime errors are caught and reported.
try {
// Load all available OpenMM plugins from their default location.
OpenMM::Platform::loadPluginsFromDirectory
(OpenMM::Platform::getDefaultPluginsDirectory());
simulateEthane();
return 0; // Normal return from main.
}
// Currently energy calculation is not available in the GPU kernels so asking
// for it requires slow Reference Platform computation at reporting intervals.
static const bool WantEnergy = true;
// Catch and report usage and runtime errors detected by OpenMM and fail.
catch(const std::exception& e) {
printf("EXCEPTION: %s\n", e.what());
return 1;
}
}
// -----------------------------------------------------------------------------
// 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.
......@@ -101,9 +75,7 @@ struct TorsionType {
} torsionType[] = {/*0 HCCH*/3, 0., 0.150};
const int HCCH = 0;
// -----------------------------------------------------------------------------
// MOLECULE DATA
// -----------------------------------------------------------------------------
// MOLECULE DATA
// Now describe an ethane molecule by listing its atoms, bonds, angles, and
// torsions. We'll provide a default configuration which centers the molecule
// at (0,0,0) with the C-C bond along the x axis.
......@@ -112,16 +84,16 @@ const int HCCH = 0;
// computer do that!
const int EndOfList=-1;
struct AtomInfo
{ int type; const char* pdb; Vec3 initPosInAngstroms;} atoms[] =
{/*0*/C, " C1 ", Vec3( -.7605, 0, 0 ),
/*1*/C, " C2 ", Vec3( .7605, 0, 0 ),
/*2*/H, "1H1 ", Vec3(-1.135, 1.03, 0 ), // bonded to C1
/*3*/H, "2H1 ", Vec3(-1.135, -.51, .89),
/*4*/H, "3H1 ", Vec3(-1.135, -.51,-.89),
/*5*/H, "1H2 ", Vec3( 1.135, 1.03, 0 ), // bonded to C2
/*6*/H, "2H2 ", Vec3( 1.135, -.51, .89),
/*7*/H, "3H2 ", Vec3( 1.135, -.51,-.89),
struct MyAtomInfo
{ int type; const char* pdb; double initPosInAng[3]; double posInAng[3];} atoms[] =
{/*0*/C, " C1 ", { -.7605, 0, 0 }, {0,0,0},
/*1*/C, " C2 ", { .7605, 0, 0 }, {0,0,0},
/*2*/H, "1H1 ", {-1.135, 1.03, 0 }, {0,0,0}, // bonded to C1
/*3*/H, "2H1 ", {-1.135, -.51, .89}, {0,0,0},
/*4*/H, "3H1 ", {-1.135, -.51,-.89}, {0,0,0},
/*5*/H, "1H2 ", { 1.135, 1.03, 0 }, {0,0,0}, // bonded to C2
/*6*/H, "2H2 ", { 1.135, -.51, .89}, {0,0,0},
/*7*/H, "3H2 ", { 1.135, -.51,-.89}, {0,0,0},
EndOfList};
static struct {int type; int atoms[2];} bonds[] =
......@@ -142,25 +114,153 @@ static struct {int type; int atoms[4];} torsions[] =
EndOfList};
// PDB FILE WRITER
// Given state data, output a single frame (pdb "model") of the trajectory.
static void
myWritePDBFrame(int frameNum, double timeInPs, double energyInKcal,
const MyAtomInfo atoms[])
{
// Write out in PDB format -- printf is so much more compact than formatted cout.
printf("MODEL %d\n", frameNum);
printf("REMARK 250 time=%.3f ps; energy=%.3f kcal/mole\n",
timeInPs, energyInKcal);
for (int n=0; atoms[n].type != EndOfList; ++n)
printf("ATOM %5d %4s ETH 1 %8.3f%8.3f%8.3f 1.00 0.00\n",
n+1, atoms[n].pdb,
atoms[n].posInAng[0], atoms[n].posInAng[1], atoms[n].posInAng[2]);
printf("ENDMDL\n");
}
// -----------------------------------------------------------------------------
// 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.
struct MyOpenMMData;
static MyOpenMMData* myInitializeOpenMM(const MyAtomInfo atoms[],
double stepSizeInFs,
std::string& platformName);
static void myStepWithOpenMM(MyOpenMMData*, int numSteps);
static void myGetOpenMMState(MyOpenMMData*, bool wantEnergy,
double& time, double& energy,
MyAtomInfo atoms[]);
static void myTerminateOpenMM(MyOpenMMData*);
// -----------------------------------------------------------------------------
// MAIN PROGRAM
// -----------------------------------------------------------------------------
int main() {
// ALWAYS enclose all OpenMM calls with a try/catch block to make sure that
// usage and runtime errors are caught and reported.
try {
std::string platformName;
// Set up OpenMM data structures; returns OpenMM Platform name.
MyOpenMMData* omm = myInitializeOpenMM(atoms, StepSizeInFs, platformName);
// 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());
const int NumSilentSteps = (int)(ReportIntervalInFs / StepSizeInFs + 0.5);
for (int frame=1; ; ++frame) {
double time, energy;
myGetOpenMMState(omm, WantEnergy, time, energy, atoms);
myWritePDBFrame(frame, time, energy, atoms);
if (time >= SimulationTimeInPs)
break;
myStepWithOpenMM(omm, NumSilentSteps);
}
// Clean up OpenMM data structures.
myTerminateOpenMM(omm);
return 0; // Normal return from main.
}
// Catch and report usage and runtime errors detected by OpenMM and fail.
catch(const std::exception& e) {
printf("EXCEPTION: %s\n", e.what());
return 1;
}
}
// Add missing scalar product operators for Vec3.
Vec3 operator*(const Vec3& v, double r) {return Vec3(v[0]*r, v[1]*r, v[2]*r);}
Vec3 operator*(double r, const Vec3& v) {return Vec3(r*v[0], r*v[1], r*v[2]);}
// 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.);
// -----------------------------------------------------------------------------
// ETHANE SIMULATION
// OpenMM-USING CODE
// -----------------------------------------------------------------------------
static void simulateEthane() {
// -------------------------------------------------------------------------
// 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.
// -----------------------------------------------------------------------------
// 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"
using OpenMM::Vec3; // so we can just say "Vec3" below
// Add definitions 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;
}
struct MyOpenMMData {
MyOpenMMData() : system(0), context(0), integrator(0) {}
~MyOpenMMData() {delete system; delete context; delete integrator;}
OpenMM::System* system;
OpenMM::OpenMMContext* context;
OpenMM::Integrator* integrator;
};
// -----------------------------------------------------------------------------
// INITIALIZE OpenMM DATA STRUCTURES
// -----------------------------------------------------------------------------
// 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( const MyAtomInfo atoms[],
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
// to each force object so we can fill in the forces. Note: OpenMM owns
// the objects and will take care of deleting them; don't do it yourself!
// -------------------------------------------------------------------------
OpenMM::System system;
// to each force object so we can fill in the forces. Note: the System owns
// the force objects and will take care of deleting them; don't do it yourself!
OpenMM::System& system = *(omm->system = new OpenMM::System());
OpenMM::NonbondedForce& nonbond = *new OpenMM::NonbondedForce();
OpenMM::HarmonicBondForce& bondStretch = *new OpenMM::HarmonicBondForce();
OpenMM::HarmonicAngleForce& bondBend = *new OpenMM::HarmonicAngleForce();
......@@ -170,29 +270,29 @@ static void simulateEthane() {
system.addForce(&bondBend);
system.addForce(&bondTorsion);
// -------------------------------------------------------------------------
// 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.
// -------------------------------------------------------------------------
std::vector<Vec3> initialPositions;
std::vector<Vec3> initialPosInNm;
for (int n=0; atoms[n].type != EndOfList; ++n) {
const AtomType& atype = atomType[atoms[n].type];
system.addParticle(atype.mass);
nonbond.addParticle(atype.charge,
atype.vdwRadiusInAngstroms * OpenMM::NmPerAngstrom
* SigmaPerVdwRadius,
* OpenMM::SigmaPerVdwRadius,
atype.vdwEnergyInKcal * OpenMM::KJPerKcal);
initialPositions.push_back(atoms[n].initPosInAngstroms * OpenMM::NmPerAngstrom);
// Convert the initial position to nm and append to the array.
const Vec3 posInNm(atoms[n].initPosInAng[0] * OpenMM::NmPerAngstrom,
atoms[n].initPosInAng[1] * OpenMM::NmPerAngstrom,
atoms[n].initPosInAng[2] * OpenMM::NmPerAngstrom);
initialPosInNm.push_back(posInNm);
}
// -------------------------------------------------------------------------
// Process the bonds:
// (1) HarmonicBondForce needs bond stretch parameters (in MD units!).
// (2) If we're using constraints, tell System about constrainable bonds.
// (3) Create a list of bonds for generating nonbond exclusions.
// -------------------------------------------------------------------------
std::vector< std::pair<int,int> > bondPairs;
for (int i=0; bonds[i].type != EndOfList; ++i) {
const int* atom = bonds[i].atoms;
......@@ -215,9 +315,7 @@ static void simulateEthane() {
// Exclude 1-2, 1-3 bonded atoms from nonbonded forces, and scale down 1-4 bonded atoms.
nonbond.createExceptionsFromBonds(bondPairs, Coulomb14Scale, LennardJones14Scale);
// -------------------------------------------------------------------------
// Create the 1-2-3 bond angle harmonic terms.
// -------------------------------------------------------------------------
for (int i=0; angles[i].type != EndOfList; ++i) {
const int* atom = angles[i].atoms;
const AngleType& angle = angleType[angles[i].type];
......@@ -228,9 +326,7 @@ static void simulateEthane() {
angle.stiffnessInKcalPerRadian2 * 2 * OpenMM::KJPerKcal);
}
// -------------------------------------------------------------------------
// Create the 1-2-3-4 bond torsion (dihedral) terms.
// -------------------------------------------------------------------------
for (int i=0; torsions[i].type != EndOfList; ++i) {
const int* atom = torsions[i].atoms;
const TorsionType& torsion = torsionType[torsions[i].type];
......@@ -240,57 +336,65 @@ static void simulateEthane() {
torsion.amplitudeInKcal * OpenMM::KJPerKcal);
}
// -------------------------------------------------------------------------
// Choose an Integrator for advancing time, and a Context connecting the
// System with the Integrator for simulation. Let the Context choose the
// best available Platform. Initialize the configuration from the default
// positions we collected above. Initial velocities will be zero.
// -------------------------------------------------------------------------
OpenMM::VerletIntegrator integrator(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);
}
omm->integrator = new OpenMM::VerletIntegrator(StepSizeInFs * OpenMM::PsPerFs);
omm->context = new OpenMM::OpenMMContext(*omm->system, *omm->integrator);
omm->context->setPositions(initialPosInNm);
platformName = omm->context->getPlatform().getName();
return omm;
}
// -----------------------------------------------------------------------------
// PDB FILE WRITER
// COPY STATE BACK TO CPU FROM OPENMM
// -----------------------------------------------------------------------------
static void
writePDB(const OpenMM::OpenMMContext& context) {
// Caution: at the moment asking for energy requires use of slow Reference
// platform calculation.
const OpenMM::State state = context.getState( OpenMM::State::Positions
| OpenMM::State::Velocities
| OpenMM::State::Energy);
const double energy = state.getPotentialEnergy() + state.getKineticEnergy();
const std::vector<Vec3>& positions = state.getPositions();
static int modelFrameNumber = 0; // numbering for MODEL records in pdb output
modelFrameNumber++;
printf("MODEL %d\n", modelFrameNumber);
printf("REMARK 250 time=%.3f picoseconds; Energy = %.3f kilojoules/mole\n",
state.getTime(), energy);
for (unsigned i=0; i < positions.size(); ++i) {
const Vec3 pos = positions[i] * OpenMM::AngstromsPerNm;
printf("ATOM %5d %4s ETH 1 %8.3f%8.3f%8.3f 1.00 0.00 \n",
i+1, atoms[i].pdb, pos[0], pos[1], pos[2]);
myGetOpenMMState(MyOpenMMData* omm, bool wantEnergy,
double& timeInPs, double& energyInKcal,
MyAtomInfo atoms[])
{
int infoMask = 0;
infoMask = OpenMM::State::Positions;
if (wantEnergy) {
infoMask += OpenMM::State::Velocities; // for kinetic energy (cheap)
infoMask += OpenMM::State::Energy; // for pot. energy (expensive)
}
printf("ENDMDL\n");
// Forces are also available (and cheap).
const OpenMM::State state = omm->context->getState(infoMask);
timeInPs = state.getTime(); // OpenMM time is in ps already
// Copy OpenMM positions into atoms array and change units from nm to Angstroms.
const std::vector<Vec3>& positionsInNm = state.getPositions();
for (int i=0; i < (int)positionsInNm.size(); ++i)
for (int j=0; j < 3; ++j)
atoms[i].posInAng[j] = positionsInNm[i][j] * OpenMM::AngstromsPerNm;
// If energy has been requested, obtain it and convert from kJ to kcal.
energyInKcal = 0;
if (wantEnergy)
energyInKcal = (state.getPotentialEnergy() + state.getKineticEnergy())
* OpenMM::KcalPerKJ;
}
// -----------------------------------------------------------------------------
// 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;
}
examples/HelloSodiumChloride.cpp
View file @ a8783de3
......@@ -129,9 +129,9 @@ int main() {
// (3) Write a PDB frame when the time comes.
printf("REMARK Using OpenMM platform %s\n", platformName.c_str());
myGetOpenMMState(omm, WantEnergy, time, energy, atoms);
myWritePDBFrame(0, time, energy, atoms);
myWritePDBFrame(1, time, energy, atoms);
for (int frame=1; frame <= NumReports; ++frame) {
for (int frame=2; frame <= NumReports; ++frame) {
myStepWithOpenMM(omm, NumSilentSteps);
myGetOpenMMState(omm, WantEnergy, time, energy, atoms);
myWritePDBFrame(frame, time, energy, atoms);
......
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