Skip to content

GitLab

Commit 2c11d35b authored by Michael Sherman's avatar Michael Sherman
Browse files

This is the "final" version of the HelloSodiumChloride example in C++. It has... 

This is the "final" version of the HelloSodiumChloride example in C++. It has now been reorganized to put all the OpenMM calls in a few high-level functions whose interfaces do not use OpenMM. Also, the periodic box is gone with GBSA implicit solvation in place instead.
parent 5b8cade7
Hide whitespace changes
Inline Side-by-side
Showing with 318 additions and 180 deletions
examples/HelloSodiumChloride.cpp
View file @ 2c11d35b
/* -----------------------------------------------------------------------------
* OpenMM(tm) HelloSodiumChloride example (May 2009)
* -----------------------------------------------------------------------------
* This is a complete, self-contained "hello world" example demonstrating
* GPU-accelerated constant energy simulation of a very simple system with just
* nonbonded forces, consisting of several sodium (Na+) and chloride (Cl-) ions.
* 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.
*
* Pay particular attention to the handling of units in this example. Incorrect
* handling of units is a very common error; this example shows how you can
* continue to work with Amber-style units of Angstroms and kCals while correctly
* 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>
using OpenMM::Vec3; // so we can just say "Vec3" below
// -----------------------------------------------------------------------------
// MODELING AND SIMULATION PARAMETERS
// -----------------------------------------------------------------------------
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 simulateNaCl();
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());
simulateNaCl();
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;
}
}
// -----------------------------------------------------------------------------
// ATOM AND FORCE FIELD DATA
// -----------------------------------------------------------------------------
// This is not part of OpenMM; just a struct we can use to collect
// atom parameters for this example. Normally atom parameters would
// come from the force field's parameterization file.
// We're going to use data in Angstrom and Kilocalorie units and
// show how to safely convert to OpenMM's internal unit system
// which uses nanometers and kilojoules.
struct AtomInfo {
const char* pdb;
double mass, charge, vdwRadiusInAng, vdwEnergyInKcal;
Vec3 initPosInAngstroms;
} atoms[] = {
// pdb mass charge vdwRadius vdwEnergy initPos
{" NA ", 22.99, 1, 1.8680, 0.00277, Vec3(8,0,0)},
{" CL ", 35.45, -1, 2.4700, 0.1000, Vec3(-8,0,0)},
{" NA ", 22.99, 1, 1.8680, 0.00277, Vec3(0,9,0)},
{" CL ", 35.45, -1, 2.4700, 0.1000, Vec3(0,-9,0)},
{" NA ", 22.99, 1, 1.8680, 0.00277, Vec3(0,0,-10)},
{" CL ", 35.45, -1, 2.4700, 0.1000, Vec3( 0,0,10)},
{""} // end of list
};
// Add missing scalar product operators for OpenMM::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.);
// -----------------------------------------------------------------------------
// NaCl SIMULATION
// -----------------------------------------------------------------------------
static void simulateNaCl() {
// -------------------------------------------------------------------------
// 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;
OpenMM::NonbondedForce* nonbond = new OpenMM::NonbondedForce();
system.addForce(nonbond);
nonbond->setNonbondedMethod(OpenMM::NonbondedForce::CutoffPeriodic);
nonbond->setCutoffDistance(2);
nonbond->setPeriodicBoxVectors(Vec3(5,0,0), Vec3(0,5,0), Vec3(0,0,5));
// -------------------------------------------------------------------------
// 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;
for (int n=0; *atoms[n].pdb; ++n) {
const AtomInfo& atom = atoms[n];
system.addParticle(atom.mass);
nonbond->addParticle(atom.charge,
atom.vdwRadiusInAng * OpenMM::NmPerAngstrom
* SigmaPerVdwRadius,
atom.vdwEnergyInKcal * OpenMM::KJPerKcal);
initialPositions.push_back(atoms[n].initPosInAngstroms
* OpenMM::NmPerAngstrom);
}
// -------------------------------------------------------------------------
// 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);
}
// -----------------------------------------------------------------------------
// PDB FILE WRITER
// -----------------------------------------------------------------------------
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
// write out in PDB format -- printf is so much more compact than formatted cout
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 SLT 1 %8.3f%8.3f%8.3f 1.00 0.00 \n",
i+1, atoms[i].pdb, pos[0], pos[1], pos[2]);
}
printf("ENDMDL\n");
}
// -----------------------------------------------------------------------------
// OpenMM(tm) HelloSodiumChloride example in C++ (June 2009)
// -----------------------------------------------------------------------------
// This is a complete, self-contained "hello world" example demonstrating
// GPU-accelerated constant temperature simulation of a very simple system with
// just nonbonded forces, consisting of several sodium (Na+) and chloride (Cl-)
// ions in implicit solvent. 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.
//
// Pay particular attention to the handling of units in this example. Incorrect
// handling of units is a very common error; this example shows how you can
// continue to work with Amber-style units like Angstroms, kCals, and van der
// Waals radii while correctly communicating with OpenMM in nm, kJ, and sigma.
// -----------------------------------------------------------------------------
#include <exception>
#include <string>
#include <cstdio>
using std::printf;
// -----------------------------------------------------------------------------
// MODELING AND SIMULATION PARAMETERS
// -----------------------------------------------------------------------------
static const double Temperature = 300; // Kelvins
static const double FrictionInPs = 1./91.; // picoseconds between collisions
static const double SolventDielectric = 80.; // typical for water
static const double SoluteDielectric = 2.; // typical for protein
static const double StepSizeInFs = 2; // integration step size (fs)
static const double ReportIntervalInFs = 10; // how often to issue PDB frame (fs)
static const double SimulationTimeInPs = 100; // total simulation time (ps)
// 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;
// -----------------------------------------------------------------------------
// ATOM AND FORCE FIELD DATA
// -----------------------------------------------------------------------------
// This is not part of OpenMM; just a struct we can use to collect atom
// parameters for this example. Normally atom parameters would come from the
// force field's parameterization file. We're going to use data in Angstrom and
// Kilocalorie units and show how to safely convert to OpenMM's internal unit
// system which uses nanometers and kilojoules.
static struct MyAtomInfo {
const char* pdb;
double mass, charge, vdwRadiusInAng, vdwEnergyInKcal,
gbsaRadiusInAng, gbsaScaleFactor;
double initPosInAng[3];
double posInAng[3]; // leave room for runtime state info
} atoms[] = {
// pdb mass charge vdwRad vdwEnergy gbsaRad gbsaScale initPos
{" NA ", 22.99, 1, 1.8680, 0.00277, 1.992, 0.8, 8, 0, 0},
{" CL ", 35.45, -1, 2.4700, 0.1000, 1.735, 0.8, -8, 0, 0},
{" NA ", 22.99, 1, 1.8680, 0.00277, 1.992, 0.8, 0, 9, 0},
{" CL ", 35.45, -1, 2.4700, 0.1000, 1.735, 0.8, 0,-9, 0},
{" NA ", 22.99, 1, 1.8680, 0.00277, 1.992, 0.8, 0, 0,-10},
{" CL ", 35.45, -1, 2.4700, 0.1000, 1.735, 0.8, 0, 0, 10},
{""} // end of list
};
// -----------------------------------------------------------------------------
// 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.
struct MyOpenMMData;
static MyOpenMMData* myInitializeOpenMM(const MyAtomInfo atoms[],
double temperature,
double frictionInPs,
double solventDielectric,
double soluteDielectric,
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*);
// -----------------------------------------------------------------------------
// 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 i=0; *atoms[i].pdb; ++i)
printf("ATOM %5d %4s SLT 1 %8.3f%8.3f%8.3f 1.00 0.00\n",
i+1, atoms[i].pdb,
atoms[i].posInAng[0], atoms[i].posInAng[1], atoms[i].posInAng[2]);
printf("ENDMDL\n");
}
// -----------------------------------------------------------------------------
// MAIN PROGRAM
// -----------------------------------------------------------------------------
int main() {
const int NumReports = (int)(SimulationTimeInPs*1000 / ReportIntervalInFs + 0.5);
const int NumSilentSteps = (int)(ReportIntervalInFs / StepSizeInFs + 0.5);
// ALWAYS enclose all OpenMM calls with a try/catch block to make sure that
// usage and runtime errors are caught and reported.
try {
double time, energy;
std::string platformName;
// Set up OpenMM data structures; returns OpenMM Platform name.
MyOpenMMData* omm = myInitializeOpenMM(atoms, Temperature, FrictionInPs,
SolventDielectric, SoluteDielectric,
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());
myGetOpenMMState(omm, WantEnergy, time, energy, atoms);
myWritePDBFrame(0, time, energy, atoms);
for (int frame=1; frame <= NumReports; ++frame) {
myStepWithOpenMM(omm, NumSilentSteps);
myGetOpenMMState(omm, WantEnergy, time, energy, atoms);
myWritePDBFrame(frame, time, energy, atoms);
}
// 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;
}
}
// -----------------------------------------------------------------------------
// 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.
// 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 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;
}
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) Allocate a MyOpenMMData structure to hang on to OpenMM data structures
// in a manner which is opaque to the caller.
// (2) Allocate the OpenMM objects which persist from call to call.
// (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 temperature,
double frictionInPs,
double solventDielectric,
double soluteDielectric,
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: the OpenMM
// System takes ownership of the force objects; don't delete them yourself.
omm->system = new OpenMM::System();
OpenMM::NonbondedForce* nonbond = new OpenMM::NonbondedForce();
OpenMM::GBSAOBCForce* gbsa = new OpenMM::GBSAOBCForce();
omm->system->addForce(nonbond);
omm->system->addForce(gbsa);
// Specify dielectrics for GBSA implicit solvation.
gbsa->setSolventDielectric(solventDielectric);
gbsa->setSoluteDielectric(soluteDielectric);
// 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) GBSA needs charge, radius, and scale factor.
// (3) Collect default positions for initializing the simulation later.
std::vector<Vec3> initialPosInNm;
for (int n=0; *atoms[n].pdb; ++n) {
const MyAtomInfo& atom = atoms[n];
omm->system->addParticle(atom.mass);
nonbond->addParticle(atom.charge,
atom.vdwRadiusInAng * OpenMM::NmPerAngstrom
* OpenMM::SigmaPerVdwRadius,
atom.vdwEnergyInKcal * OpenMM::KJPerKcal);
gbsa->addParticle(atom.charge,
atom.gbsaRadiusInAng * OpenMM::NmPerAngstrom,
atom.gbsaScaleFactor);
const Vec3 posInNm(atom.initPosInAng[0] * OpenMM::NmPerAngstrom,
atom.initPosInAng[1] * OpenMM::NmPerAngstrom,
atom.initPosInAng[2] * OpenMM::NmPerAngstrom);
initialPosInNm.push_back(posInNm);
}
// 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 but could
// have been set here.
omm->integrator = new OpenMM::LangevinIntegrator(temperature, frictionInPs,
stepSizeInFs * OpenMM::PsPerFs);
omm->context = new OpenMM::OpenMMContext(*omm->system, *omm->integrator);
omm->context->setPositions(initialPosInNm);
platformName = omm->context->getPlatform().getName();
return omm;
}
// -----------------------------------------------------------------------------
// COPY STATE BACK TO CPU FROM OPENMM
// -----------------------------------------------------------------------------
static void
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
infoMask += OpenMM::State::Energy; // for potential energy (expensive)
}
// Forces are also available (and NOT expensive).
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;
}
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