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

Minor tweaks and comments.

parent f337005b
...@@ -8,87 +8,87 @@ ...@@ -8,87 +8,87 @@
// other visualization tool to produce an animation of the resulting trajectory. // other visualization tool to produce an animation of the resulting trajectory.
// ----------------------------------------------------------------------------- // -----------------------------------------------------------------------------
#include "OpenMM.h" #include "OpenMM.h"
#include <cstdio> #include <cstdio>
// Forward declaration of writePdb() subroutine for printing atomic // Forward declaration of writePdb() subroutine for printing atomic
// coordinates, defined later in this source file. // coordinates, defined later in this source file.
void writePdb(const OpenMM::OpenMMContext& context); void writePdb(const OpenMM::OpenMMContext& context);
void simulateArgon() void simulateArgon()
{ {
// Load any shared libraries containing GPU implementations // Load any shared libraries containing GPU implementations
OpenMM::Platform::loadPluginsFromDirectory( OpenMM::Platform::loadPluginsFromDirectory(
OpenMM::Platform::getDefaultPluginsDirectory()); OpenMM::Platform::getDefaultPluginsDirectory());
// Create a system with nonbonded forces // Create a system with nonbonded forces
OpenMM::System system; OpenMM::System system;
OpenMM::NonbondedForce* nonbond = new OpenMM::NonbondedForce(); OpenMM::NonbondedForce* nonbond = new OpenMM::NonbondedForce();
system.addForce(nonbond); system.addForce(nonbond);
// Create three atoms // Create three atoms
std::vector<OpenMM::Vec3> initialPositions(3); std::vector<OpenMM::Vec3> initialPositions(3);
for (int a = 0; a < 3; ++a) for (int a = 0; a < 3; ++a)
{ {
system.addParticle(39.95); // mass, grams per mole system.addParticle(39.95); // mass, grams per mole
// charge, sigma, well depth // charge, sigma, well depth
nonbond->addParticle(0.0, 0.3350, 0.996); nonbond->addParticle(0.0, 0.3350, 0.996);
initialPositions[a] = OpenMM::Vec3(0.5*a,0,0); // location, nanometers initialPositions[a] = OpenMM::Vec3(0.5*a,0,0); // location, nanometers
} }
OpenMM::VerletIntegrator integrator(0.004); // step size in picoseconds OpenMM::VerletIntegrator integrator(0.004); // step size in picoseconds
// Let OpenMM Context choose best platform. // Let OpenMM Context choose best platform.
OpenMM::OpenMMContext context(system, integrator); OpenMM::OpenMMContext context(system, integrator);
printf( "REMARK Using OpenMM platform %s\n", printf( "REMARK Using OpenMM platform %s\n",
context.getPlatform().getName().c_str() ); context.getPlatform().getName().c_str() );
// Set the starting positions of the atoms. Velocities will be zero. // Set the starting positions of the atoms. Velocities will be zero.
context.setPositions(initialPositions); context.setPositions(initialPositions);
// Simulate // Simulate
while(context.getTime() < 10.0) { // picoseconds while(context.getTime() < 10.0) { // picoseconds
writePdb(context); // output coordinates writePdb(context); // output coordinates
// Run 100 steps at a time, for efficient use of OpenMM // Run 100 steps at a time, for efficient use of OpenMM
integrator.step(100); integrator.step(100);
} }
writePdb(context); // output final coordinates writePdb(context); // output final coordinates
} }
int main() int main()
{ {
try { try {
simulateArgon(); simulateArgon();
return 0; // success! return 0; // success!
} }
// Catch and report usage and runtime errors detected by OpenMM and fail. // Catch and report usage and runtime errors detected by OpenMM and fail.
catch(const std::exception& e) { catch(const std::exception& e) {
printf("EXCEPTION: %s\n", e.what()); printf("EXCEPTION: %s\n", e.what());
return 1; // failure! return 1; // failure!
} }
} }
// writePdb() subroutine for quick-and-dirty trajectory output. // writePdb() subroutine for quick-and-dirty trajectory output.
void writePdb(const OpenMM::OpenMMContext& context) void writePdb(const OpenMM::OpenMMContext& context)
{ {
// Request atomic positions from OpenMM // Request atomic positions from OpenMM
const OpenMM::State state = context.getState(OpenMM::State::Positions); const OpenMM::State state = context.getState(OpenMM::State::Positions);
const std::vector<OpenMM::Vec3>& pos = state.getPositions(); const std::vector<OpenMM::Vec3>& pos = state.getPositions();
// write out in PDB format // write out in PDB format
// Use PDB MODEL cards to number trajectory frames // Use PDB MODEL cards to number trajectory frames
static int modelFrameNumber = 0; static int modelFrameNumber = 0;
modelFrameNumber++; modelFrameNumber++;
printf("MODEL %d\n", modelFrameNumber); // start of frame printf("MODEL %d\n", modelFrameNumber); // start of frame
for (int a = 0; a < context.getSystem().getNumParticles(); ++a) for (int a = 0; a < context.getSystem().getNumParticles(); ++a)
{ {
printf("ATOM %5d AR AR 1 ", a+1); // atom number printf("ATOM %5d AR AR 1 ", a+1); // atom number
printf("%8.3f%8.3f%8.3f 1.00 0.00 AR\n", // coordinates printf("%8.3f%8.3f%8.3f 1.00 0.00 AR\n", // coordinates
// "*10" converts nanometers to Angstroms // "*10" converts nanometers to Angstroms
pos[a][0]*10, pos[a][1]*10, pos[a][2]*10); pos[a][0]*10, pos[a][1]*10, pos[a][2]*10);
} }
printf("ENDMDL\n"); // end of frame printf("ENDMDL\n"); // end of frame
} }
...@@ -187,9 +187,9 @@ struct MyOpenMMData { ...@@ -187,9 +187,9 @@ struct MyOpenMMData {
// INITIALIZE OpenMM DATA STRUCTURES // INITIALIZE OpenMM DATA STRUCTURES
// ----------------------------------------------------------------------------- // -----------------------------------------------------------------------------
// We take these actions here: // We take these actions here:
// (1) Allocate a MyOpenMMData structure to hang on to OpenMM data structures // (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. // 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 // (3) Fill the OpenMM::System with the force field parameters we want to
// use and the particular set of atoms to be simulated. // use and the particular set of atoms to be simulated.
// (4) Create an Integrator and a Context associating the Integrator with // (4) Create an Integrator and a Context associating the Integrator with
...@@ -232,7 +232,7 @@ myInitializeOpenMM( const MyAtomInfo atoms[], ...@@ -232,7 +232,7 @@ myInitializeOpenMM( const MyAtomInfo atoms[],
// (1) System needs to know the masses. // (1) System needs to know the masses.
// (2) NonbondedForce needs charges,van der Waals properties (in MD units!). // (2) NonbondedForce needs charges,van der Waals properties (in MD units!).
// (3) GBSA needs charge, radius, and scale factor. // (3) GBSA needs charge, radius, and scale factor.
// (3) Collect default positions for initializing the simulation later. // (4) Collect default positions for initializing the simulation later.
std::vector<Vec3> initialPosInNm; std::vector<Vec3> initialPosInNm;
for (int n=0; *atoms[n].pdb; ++n) { for (int n=0; *atoms[n].pdb; ++n) {
const MyAtomInfo& atom = atoms[n]; const MyAtomInfo& atom = atoms[n];
......
...@@ -184,7 +184,7 @@ struct MyOpenMMData_s { ...@@ -184,7 +184,7 @@ struct MyOpenMMData_s {
* *
* Note that this function must understand the calling MD code's molecule and * 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. * force field data structures so will need to be customized for each MD code.
*/ */
static MyOpenMMData* static MyOpenMMData*
myInitializeOpenMM( const MyAtomInfo atoms[], myInitializeOpenMM( const MyAtomInfo atoms[],
double temperature, double temperature,
...@@ -274,7 +274,7 @@ myGetOpenMMState(MyOpenMMData* omm, int wantEnergy, ...@@ -274,7 +274,7 @@ myGetOpenMMState(MyOpenMMData* omm, int wantEnergy,
MyAtomInfo atoms[]) MyAtomInfo atoms[])
{ {
OpenMM_State* state; OpenMM_State* state;
const OpenMM_Vec3Array* posArray; const OpenMM_Vec3Array* posArrayInNm;
int infoMask; int infoMask;
int n; int n;
...@@ -291,12 +291,11 @@ myGetOpenMMState(MyOpenMMData* omm, int wantEnergy, ...@@ -291,12 +291,11 @@ myGetOpenMMState(MyOpenMMData* omm, int wantEnergy,
/* Positions are maintained as a Vec3Array inside the State. This will give /* Positions are maintained as a Vec3Array inside the State. This will give
* us access, but don't destroy it yourself -- it will go away with the State. */ * us access, but don't destroy it yourself -- it will go away with the State. */
posArray = OpenMM_State_getPositions(state); posArrayInNm = OpenMM_State_getPositions(state);
for (n=0; *atoms[n].pdb; ++n) { for (n=0; *atoms[n].pdb; ++n)
double posInNm[3]; /* Sets atoms[n].pos = posArray[n] * Angstroms/nm. */
OpenMM_Vec3Array_get(posArray, n, posInNm); OpenMM_Vec3Array_getScaled(posArrayInNm, n, OpenMM_AngstromsPerNm,
OpenMM_Vec3_scale(posInNm, OpenMM_AngstromsPerNm, atoms[n].posInAng); atoms[n].posInAng);
}
/* If energy has been requested, obtain it and convert from kJ to kcal. */ /* If energy has been requested, obtain it and convert from kJ to kcal. */
*energyInKcal = 0; *energyInKcal = 0;
......
...@@ -299,7 +299,7 @@ void OpenMM_GBSAOBCForce_destroy(OpenMM_GBSAOBCForce* doomed) ...@@ -299,7 +299,7 @@ void OpenMM_GBSAOBCForce_destroy(OpenMM_GBSAOBCForce* doomed)
void openmm_gbsaobcforce_destroy_(OpenMM_GBSAOBCForce*& doomed) void openmm_gbsaobcforce_destroy_(OpenMM_GBSAOBCForce*& doomed)
{ OpenMM_GBSAOBCForce_destroy(doomed); doomed = 0;} { OpenMM_GBSAOBCForce_destroy(doomed); doomed = 0;}
// Fortran only: recast NonbondedForce as a Force. // Fortran only: recast GBSAOBCForce as a Force.
void openmm_gbsaobcforce_asforce_(OpenMM_GBSAOBCForce* const& gbsa, void openmm_gbsaobcforce_asforce_(OpenMM_GBSAOBCForce* const& gbsa,
OpenMM_Force*& force) OpenMM_Force*& force)
{ force = (OpenMM_Force*)gbsa; } { force = (OpenMM_Force*)gbsa; }
......
...@@ -129,7 +129,7 @@ extern void OpenMM_Vec3Array_get(const OpenMM_Vec3Array*, int i, ...@@ -129,7 +129,7 @@ extern void OpenMM_Vec3Array_get(const OpenMM_Vec3Array*, int i,
extern void OpenMM_Vec3Array_getScaled(const OpenMM_Vec3Array*, int i, double s, double[3]); extern void OpenMM_Vec3Array_getScaled(const OpenMM_Vec3Array*, int i, double s, double[3]);
extern void OpenMM_Vec3Array_set(OpenMM_Vec3Array*, int i, const double[3]); extern void OpenMM_Vec3Array_set(OpenMM_Vec3Array*, int i, const double[3]);
extern void OpenMM_Vec3Array_setScaled(OpenMM_Vec3Array*, int i, const double[3], double s); extern void OpenMM_Vec3Array_setScaled(OpenMM_Vec3Array*, int i, const double[3], double s);
extern void OpenMM_Vec3_scale(const double[3], double s, double[3]); extern void OpenMM_Vec3_scale(const double in[3], double s, double out[3]);
/* OpenMM_String */ /* OpenMM_String */
extern OpenMM_String* OpenMM_String_create(const char* init); extern OpenMM_String* OpenMM_String_create(const char* init);
......
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