Skip to content

GitLab

Commit 39d85db0 authored by Michael Sherman's avatar Michael Sherman
Browse files

Reworked the NaCl example per dryrun feedback and fixed some warnings in ethane example.

parent 7fccb08e
Hide whitespace changes
Inline Side-by-side
Showing with 422 additions and 366 deletions
examples/HelloEthane.cpp
View file @ 39d85db0
/* -----------------------------------------------------------------------------
* OpenMM(tm) HelloEthane example (May 2009)
* -----------------------------------------------------------------------------
* This is a complete, self-contained "hello world" example demonstrating
* GPU-accelerated simulation of a system with both bonded and nonbonded forces,
* using ethane (H3-C-C-H3) in a vacuum as an example. 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.
/* -----------------------------------------------------------------------------
* OpenMM(tm) HelloEthane example (May 2009)
* -----------------------------------------------------------------------------
* This is a complete, self-contained "hello world" example demonstrating
* GPU-accelerated simulation of a system with both bonded and nonbonded forces,
* using ethane (H3-C-C-H3) in a vacuum as an example. 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 <iostream>
#include <string>
#include <vector>
#include <utility>
using namespace OpenMM;
// -----------------------------------------------------------------------------
// 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 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.
Platform::loadPluginsFromDirectory(Platform::getDefaultPluginsDirectory());
simulateEthane();
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;
}
}
// -----------------------------------------------------------------------------
// 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.
// 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.
const double Coulomb14Scale = 0.5;
const double LennardJones14Scale = 0.5;
struct AtomType {
double mass, charge, vdwRadiusInAngstroms, vdwEnergyInKcal;
} atomType[] = {/*0 H*/ 1.008, 0.0605, 1.4870, 0.0157,
/*1 C*/12.011, -.1815, 1.9080, 0.1094};
const int H = 0, C = 1;
struct BondType {
double nominalLengthInAngstroms, stiffnessInKcalPerAngstrom2;
bool canConstrain;
} bondType[] = {/*0 CC*/1.526, 310., false,
/*1 CH*/1.09 , 340., true};
const int CC = 0, CH = 1;
struct AngleType {
double nominalAngleInDegrees, stiffnessInKcalPerRadian2;
} angleType[] = {/*0 HCC*/109.5, 50.,
/*1 HCH*/109.5, 35.};
const int HCC = 0, HCH = 1;
struct TorsionType {
int periodicity;
double phaseInDegrees, amplitudeInKcal;
} torsionType[] = {/*0 HCCH*/3, 0., 0.150};
const int HCCH = 0;
// -----------------------------------------------------------------------------
// 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.
// Use this as an "end of list" marker so that we do not have to count; let the
// computer do that!
const int EndOfList=-1;
struct AtomInfo
{ int type; char pdb[5]; 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),
EndOfList};
struct {int type; int atoms[2];} bonds[] =
{CC,0,1,
CH,0,2,CH,0,3,CH,0,4, // C1 methyl
CH,1,5,CH,1,6,CH,1,7, // C2 methyl
EndOfList};
struct {int type; int atoms[3];} angles[] =
{HCC,2,0,1,HCC,3,0,1,HCC,4,0,1, // C1 methyl
HCH,2,0,3,HCH,2,0,4,HCH,3,0,4,
HCC,5,1,0,HCC,6,1,0,HCC,7,1,0, // C2 methyl
HCH,5,1,6,HCH,5,1,7,HCH,6,1,7,
EndOfList};
struct {int type; int atoms[4];} torsions[] =
{HCCH,2,0,1,5,HCCH,2,0,1,6,HCCH,2,0,1,7,
HCCH,3,0,1,5,HCCH,3,0,1,6,HCCH,3,0,1,7,
HCCH,4,0,1,5,HCCH,4,0,1,6,HCCH,4,0,1,7,
EndOfList};
// 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
// -----------------------------------------------------------------------------
static void simulateEthane() {
// -------------------------------------------------------------------------
// 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!
// -------------------------------------------------------------------------
System system;
NonbondedForce& nonbond = *new NonbondedForce();
HarmonicBondForce& bondStretch = *new HarmonicBondForce();
HarmonicAngleForce& bondBend = *new HarmonicAngleForce();
PeriodicTorsionForce& bondTorsion = *new PeriodicTorsionForce();
system.addForce(&nonbond);
system.addForce(&bondStretch);
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;
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 * NmPerAngstrom * SigmaPerVdwRadius,
atype.vdwEnergyInKcal * KJPerKcal);
initialPositions.push_back(atoms[n].initPosInAngstroms * NmPerAngstrom);
}
// -------------------------------------------------------------------------
// 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;
const BondType& bond = bondType[bonds[i].type];
// 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
// it as used in the force term kx, with energy kx^2/2.
bondStretch.addBond(atom[0], atom[1],
bond.nominalLengthInAngstroms * NmPerAngstrom,
bond.stiffnessInKcalPerAngstrom2 * 2 * KJPerKcal * AngstromsPerNm * AngstromsPerNm);
if (UseConstraints && bond.canConstrain)
system.addConstraint(atom[0], atom[1],
bond.nominalLengthInAngstroms * NmPerAngstrom);
bondPairs.push_back(std::make_pair(atom[0], atom[1]));
}
// 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];
// See note under bond stretch above regarding the factor of 2 here.
bondBend.addAngle(atom[0],atom[1],atom[2],
angle.nominalAngleInDegrees * RadiansPerDegree,
angle.stiffnessInKcalPerRadian2 * 2 * 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];
bondTorsion.addTorsion(atom[0],atom[1],atom[2],atom[3],
torsion.periodicity,
torsion.phaseInDegrees * RadiansPerDegree,
torsion.amplitudeInKcal * 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.
// -------------------------------------------------------------------------
VerletIntegrator integrator(StepSizeInFs * PsPerFs);
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 OpenMMContext& context) {
// Caution: at the moment asking for energy requires use of slow Reference
// platform calculation.
const State state = context.getState(State::Positions | State::Velocities | 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] * 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]);
}
printf("ENDMDL\n");
}
#endif
#include "OpenMM.h"
#include <iostream>
#include <string>
#include <vector>
#include <utility>
using namespace OpenMM;
// -----------------------------------------------------------------------------
// 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 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.
Platform::loadPluginsFromDirectory(Platform::getDefaultPluginsDirectory());
simulateEthane();
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;
}
}
// -----------------------------------------------------------------------------
// 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.
// 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.
const double Coulomb14Scale = 0.5;
const double LennardJones14Scale = 0.5;
struct AtomType {
double mass, charge, vdwRadiusInAngstroms, vdwEnergyInKcal;
} atomType[] = {/*0 H*/ 1.008, 0.0605, 1.4870, 0.0157,
/*1 C*/12.011, -.1815, 1.9080, 0.1094};
const int H = 0, C = 1;
struct BondType {
double nominalLengthInAngstroms, stiffnessInKcalPerAngstrom2;
bool canConstrain;
} bondType[] = {/*0 CC*/1.526, 310., false,
/*1 CH*/1.09 , 340., true};
const int CC = 0, CH = 1;
struct AngleType {
double nominalAngleInDegrees, stiffnessInKcalPerRadian2;
} angleType[] = {/*0 HCC*/109.5, 50.,
/*1 HCH*/109.5, 35.};
const int HCC = 0, HCH = 1;
struct TorsionType {
int periodicity;
double phaseInDegrees, amplitudeInKcal;
} torsionType[] = {/*0 HCCH*/3, 0., 0.150};
const int HCCH = 0;
// -----------------------------------------------------------------------------
// 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.
// Use this as an "end of list" marker so that we do not have to count; let the
// 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),
EndOfList};
static struct {int type; int atoms[2];} bonds[] =
{CC,0,1,
CH,0,2,CH,0,3,CH,0,4, // C1 methyl
CH,1,5,CH,1,6,CH,1,7, // C2 methyl
EndOfList};
static struct {int type; int atoms[3];} angles[] =
{HCC,2,0,1,HCC,3,0,1,HCC,4,0,1, // C1 methyl
HCH,2,0,3,HCH,2,0,4,HCH,3,0,4,
HCC,5,1,0,HCC,6,1,0,HCC,7,1,0, // C2 methyl
HCH,5,1,6,HCH,5,1,7,HCH,6,1,7,
EndOfList};
static struct {int type; int atoms[4];} torsions[] =
{HCCH,2,0,1,5,HCCH,2,0,1,6,HCCH,2,0,1,7,
HCCH,3,0,1,5,HCCH,3,0,1,6,HCCH,3,0,1,7,
HCCH,4,0,1,5,HCCH,4,0,1,6,HCCH,4,0,1,7,
EndOfList};
// 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
// -----------------------------------------------------------------------------
static void simulateEthane() {
// -------------------------------------------------------------------------
// 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!
// -------------------------------------------------------------------------
System system;
NonbondedForce& nonbond = *new NonbondedForce();
HarmonicBondForce& bondStretch = *new HarmonicBondForce();
HarmonicAngleForce& bondBend = *new HarmonicAngleForce();
PeriodicTorsionForce& bondTorsion = *new PeriodicTorsionForce();
system.addForce(&nonbond);
system.addForce(&bondStretch);
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;
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 * NmPerAngstrom * SigmaPerVdwRadius,
atype.vdwEnergyInKcal * KJPerKcal);
initialPositions.push_back(atoms[n].initPosInAngstroms * NmPerAngstrom);
}
// -------------------------------------------------------------------------
// 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;
const BondType& bond = bondType[bonds[i].type];
// 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
// it as used in the force term kx, with energy kx^2/2.
bondStretch.addBond(atom[0], atom[1],
bond.nominalLengthInAngstroms * NmPerAngstrom,
bond.stiffnessInKcalPerAngstrom2 * 2 * KJPerKcal * AngstromsPerNm * AngstromsPerNm);
if (UseConstraints && bond.canConstrain)
system.addConstraint(atom[0], atom[1],
bond.nominalLengthInAngstroms * NmPerAngstrom);
bondPairs.push_back(std::make_pair(atom[0], atom[1]));
}
// 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];
// See note under bond stretch above regarding the factor of 2 here.
bondBend.addAngle(atom[0],atom[1],atom[2],
angle.nominalAngleInDegrees * RadiansPerDegree,
angle.stiffnessInKcalPerRadian2 * 2 * 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];
bondTorsion.addTorsion(atom[0],atom[1],atom[2],atom[3],
torsion.periodicity,
torsion.phaseInDegrees * RadiansPerDegree,
torsion.amplitudeInKcal * 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.
// -------------------------------------------------------------------------
VerletIntegrator integrator(StepSizeInFs * PsPerFs);
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 OpenMMContext& context) {
// Caution: at the moment asking for energy requires use of slow Reference
// platform calculation.
const State state = context.getState(State::Positions | State::Velocities | 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] * 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]);
}
printf("ENDMDL\n");
}
examples/HelloSodiumChloride.cpp
View file @ 39d85db0
/* -----------------------------------------------------------------------------
* 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 <iostream>
#include <iomanip>
#include <string>
using namespace OpenMM;
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]);
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 {
char* symbol;
double mass, charge, vdwRadiusAng, vdwEnergyKcal;
Vec3 startPosAng;
};
static AtomInfo atoms[] = {
{"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)},
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
};
static const double Temperature = 300; // Kelvins
static const double Friction = 1./91.; // picoseconds between collisions
static const double StepSizeFs = 2; // femtoseconds
static const double ReportIntervalFs = 10;
static const double SimulationTimePs = 100; // total simulation time (ps)
static const double SigmaPerVdwRadius = 2*std::pow(2., -1./6.);
static void writePDB(const OpenMMContext&);
int main() {
try {
// Load all available OpenMM plugins from their default location.
Platform::loadPluginsFromDirectory(Platform::getDefaultPluginsDirectory());
// Create a System and a NonbondedForce object within the System.
System system;
NonbondedForce* nonbond = new NonbondedForce();
// 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);
int numAtoms = 0;
for (; *atoms[numAtoms].symbol; ++numAtoms) {
const AtomInfo& atom = atoms[numAtoms];
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.vdwRadiusAng * NmPerAngstrom * SigmaPerVdwRadius,
atom.vdwEnergyKcal * KJPerKcal);
atom.vdwRadiusInAng * OpenMM::NmPerAngstrom * SigmaPerVdwRadius,
atom.vdwEnergyInKcal * OpenMM::KJPerKcal);
initialPositions.push_back(atoms[n].initPosInAngstroms * OpenMM::NmPerAngstrom);
}
// Create an integrator object for advancing time.
LangevinIntegrator integrator(Temperature, Friction, StepSizeFs * PsPerFs);
//VerletIntegrator integrator(StepSizeFs * PsPerFs);
// Create an OpenMM Context for execution; let it choose best platform.
OpenMMContext context(system, integrator);
const std::string platformName = context.getPlatform().getName();
// -------------------------------------------------------------------------
// 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() );
// Fill in a vector of starting positions, one per atom.
std::vector<Vec3> positions(numAtoms);
for (int i=0; i < numAtoms; ++i)
positions[i] = atoms[i].startPosAng * NmPerAngstrom;
// Set the starting positions in the OpenMM context. Velocities will be zero.
context.setPositions(positions);
// Output the initial state.
writePDB(context);
const int NumSilentSteps = (int)(ReportIntervalFs / StepSizeFs + 0.5);
const int NumSilentSteps = (int)(ReportIntervalInFs / StepSizeInFs + 0.5);
do {
integrator.step(NumSilentSteps);
writePDB(context);
} while (context.getTime() < SimulationTimePs);
return 0;
} catch(const std::exception& e) {
std::cout << "EXCEPTION: " << e.what() << std::endl;
return 1;
}
} while (context.getTime() < SimulationTimeInPs);
}
// -----------------------------------------------------------------------------
// PDB FILE WRITER
// -----------------------------------------------------------------------------
static void
writePDB(const OpenMMContext& context) {
writePDB(const OpenMM::OpenMMContext& context) {
// Caution: at the moment asking for energy requires use of slow reference calculation.
const State state = context.getState(State::Positions | State::Velocities | State::Energy);
const double energy = state.getPotentialEnergy() + state.getKineticEnergy();
const std::vector<Vec3>& positions = state.getPositions();
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);
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] * AngstromsPerNm;
printf("ATOM %3d %2s SLT 1 %8.3f%8.3f%8.3f 1.00 0.00 %2s\n",
i+1, atoms[i].symbol, pos[0], pos[1], pos[2], atoms[i].symbol);
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");
}
......
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