Skip to content

GitLab

Commit bb42a5dd authored by Mark Friedrichs's avatar Mark Friedrichs
Browse files

Update for getState()

parent aff7f7bd
Hide whitespace changes
Inline Side-by-side
Showing with 19 additions and 451 deletions
examples/HelloArgonInC.c
View file @ bb42a5dd
/* ----------------------------------------------------------------------------- Vim: Warning: Output is not to a terminal
* OpenMM(tm) HelloArgon example in C (June 2009) [?1049h[?1h=[?12;25h[?12l[?25h[?25l"svn-commit.2.tmp" 5L, 119C 1
* ----------------------------------------------------------------------------- 2 --This line, and those below, will be ignored--
* This program demonstrates a simple molecular simulation using the OpenMM  3
* API for GPU-accelerated molecular dynamics simulation. The primary goal is 4 M examples/HelloArgonInC.c
* to make sure you can compile, link, and run with OpenMM and view the output.  5 M examples/HelloSodiumChlorideInC.c
* The example is available in C++, C, and Fortran 95. ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 1,0-1All[?12l[?25h[?25l2,1 [?12l[?25h[?25l1,0-1[?12l[?25h[?25l-- INSERT --1,1All[?12l[?25h[?25lM2[?12l[?25h[?25l1[?12l[?25h[?25lU2[?12l[?25h[?25lp3[?12l[?25h[?25ld4[?12l[?25h[?25la5[?12l[?25h[?25lt6[?12l[?25h[?25le7[?12l[?25h[?25l8[?12l[?25h[?25lf9[?12l[?25h[?25lo10[?12l[?25h[?25lr1[?12l[?25h[?25l2[?12l[?25h[?25lg3[?12l[?25h[?25le4[?12l[?25h[?25lt5[?12l[?25h[?25lS6[?12l[?25h[?25lt7[?12l[?25h[?25la8[?12l[?25h[?25lt9[?12l[?25h[?25le20[?12l[?25h[?25l(1[?12l[?25h[?25l)0[?12l[?25h[?25l()1,21All[?12l[?25h[?25l:[?12l[?25hwq[?25l"svn-commit.2.tmp" 5L, 140C written
* [?1l>[?12l[?25h[?1049lSending examples/HelloArgonInC.c
* The system modeled here is a small number of argon atoms in a vacuum. Sending examples/HelloSodiumChlorideInC.c
* 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.
* -------------------------------------------------------------------------- */
#include "OpenMMCWrapper.h"
#include <stdio.h>
/* Forward declaration of routine for printing one frame of the
trajectory, defined later in this source file. */
void writePdbFrame(int frameNum, const OpenMM_State*);
void simulateArgon()
{
OpenMM_System* system;
OpenMM_Integrator* integrator;
OpenMM_Context* context;
OpenMM_Platform* platform;
OpenMM_NonbondedForce* nonbond;
OpenMM_Vec3Array* initPosInNm;
OpenMM_StringArray* pluginList;
int a, frameNum;
/* Load any shared libraries containing GPU implementations. */
pluginList = OpenMM_Platform_loadPluginsFromDirectory(
OpenMM_Platform_getDefaultPluginsDirectory());
OpenMM_StringArray_destroy(pluginList);
/* Create a system with nonbonded forces. System takes ownership
of Force; don't destroy it yourself. */
system = OpenMM_System_create();
nonbond = OpenMM_NonbondedForce_create();
OpenMM_System_addForce(system, (OpenMM_Force*)nonbond);
/* Create three atoms. */
initPosInNm = OpenMM_Vec3Array_create(3);
for (a = 0; a < 3; ++a)
{
const OpenMM_Vec3 posNm = {0.5*a, 0, 0}; /*location, nm*/
OpenMM_Vec3Array_set(initPosInNm, a, posNm);
OpenMM_System_addParticle(system, 39.95); /*mass of Ar, grams/mole*/
/* charge, L-J sigma (nm), well depth (kJ) */
OpenMM_NonbondedForce_addParticle(nonbond, 0.0, 0.3350, 0.996); /*vdWRad(Ar)=.188 nm*/
}
/* Create particular integrator, and recast to generic one. */
integrator = (OpenMM_Integrator*)OpenMM_VerletIntegrator_create(0.004); /*step size in ps*/
/* Let OpenMM Context choose best platform. */
context = OpenMM_Context_create(system, integrator);
platform = OpenMM_Context_getPlatform(context);
printf( "REMARK Using OpenMM platform %s\n",
OpenMM_Platform_getName(platform));
/* Set starting positions of the atoms. Leave time and velocity zero. */
OpenMM_Context_setPositions(context, initPosInNm);
/* Simulate. */
for (frameNum=1; ;++frameNum) {
/* Output current state information. */
OpenMM_State* state = OpenMM_Context_getState(context, OpenMM_State_Positions, 0);
const double timeInPs = OpenMM_State_getTime(state);
writePdbFrame(frameNum, state); /*output coordinates*/
OpenMM_State_destroy(state);
if (timeInPs >= 10.)
break;
/* Advance state many steps at a time, for efficient use of OpenMM. */
OpenMM_Integrator_step(integrator, 10); /*(use a lot more than this normally)*/
}
/* Free heap space for all the objects created above. */
OpenMM_Vec3Array_destroy(initPosInNm);
OpenMM_Context_destroy(context);
OpenMM_Integrator_destroy(integrator);
OpenMM_System_destroy(system);
}
int main()
{
simulateArgon();
return 0;
}
/* Handy homebrew PDB writer for quick-and-dirty trajectory output. */
void writePdbFrame(int frameNum, const OpenMM_State* state)
{
int a;
/* Reference atomic positions in the OpenMM State. */
const OpenMM_Vec3Array* posInNm = OpenMM_State_getPositions(state);
/* Use PDB MODEL cards to number trajectory frames. */
printf("MODEL %d\n", frameNum); /*start of frame*/
for (a = 0; a < OpenMM_Vec3Array_getSize(posInNm); ++a)
{
OpenMM_Vec3 posInAng;
/* "10" here converts nanometers to Angstroms */
posInAng = OpenMM_Vec3_scale(*OpenMM_Vec3Array_get(posInNm, a), 10.);
printf("ATOM %5d AR AR 1 ", a+1); /*atom number*/
printf("%8.3f%8.3f%8.3f 1.00 0.00\n", /*coordinates*/
posInAng.x, posInAng.y, posInAng.z);
}
printf("ENDMDL\n"); /*end of frame*/
}
......
examples/HelloSodiumChlorideInC.c
View file @ bb42a5dd
/* ----------------------------------------------------------------------------- Vim: Warning: Output is not to a terminal
* OpenMM(tm) HelloSodiumChloride example in C (June 2009) [?1049h[?1h=[?12;25h[?12l[?25h[?25l"svn-commit.2.tmp" 5L, 119C 1
* ----------------------------------------------------------------------------- 2 --This line, and those below, will be ignored--
* This is a complete, self-contained "hello world" example demonstrating  3
* GPU-accelerated constant temperature simulation of a very simple system with 4 M examples/HelloArgonInC.c
* just nonbonded forces, consisting of several sodium (Na+) and chloride (Cl-)  5 M examples/HelloSodiumChlorideInC.c
* ions in implicit solvent. A multi-frame PDB file is written to stdout which ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 1,0-1All[?12l[?25h[?25l2,1 [?12l[?25h[?25l1,0-1[?12l[?25h[?25l-- INSERT --1,1All[?12l[?25h[?25lM2[?12l[?25h[?25l1[?12l[?25h[?25lU2[?12l[?25h[?25lp3[?12l[?25h[?25ld4[?12l[?25h[?25la5[?12l[?25h[?25lt6[?12l[?25h[?25le7[?12l[?25h[?25l8[?12l[?25h[?25lf9[?12l[?25h[?25lo10[?12l[?25h[?25lr1[?12l[?25h[?25l2[?12l[?25h[?25lg3[?12l[?25h[?25le4[?12l[?25h[?25lt5[?12l[?25h[?25lS6[?12l[?25h[?25lt7[?12l[?25h[?25la8[?12l[?25h[?25lt9[?12l[?25h[?25le20[?12l[?25h[?25l(1[?12l[?25h[?25l)0[?12l[?25h[?25l()1,21All[?12l[?25h[?25l:[?12l[?25hwq[?25l"svn-commit.2.tmp" 5L, 140C written
* can be read by VMD or other visualization tool to produce an animation of the [?1l>[?12l[?25h[?1049lSending examples/HelloArgonInC.c
* resulting trajectory. Sending examples/HelloSodiumChlorideInC.c
* Transmitting file data ..
* Pay particular attention to the handling of units in this example. Incorrect \ No newline at end of file
* 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.
*
* This example is written entirely in ANSI C, using the OpenMM C bindings.
* -------------------------------------------------------------------------- */
#include <stdio.h>
#include <stdlib.h>
/* --------------------------------------------------------------------------
* MODELING AND SIMULATION PARAMETERS
* -------------------------------------------------------------------------- */
static const double Temperature = 300; /*Kelvins */
static const double FrictionInPerPs = 91.; /*collisions per ps*/
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 = 50; /*how often for 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 int WantEnergy = 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.
*/
typedef struct MyAtomInfo_s {
const char* pdb;
double mass, charge, vdwRadiusInAng, vdwEnergyInKcal,
gbsaRadiusInAng, gbsaScaleFactor;
double initPosInAng[3];
double posInAng[3]; /*leave room for runtime state info*/
} MyAtomInfo;
static MyAtomInfo 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.
*/
typedef struct MyOpenMMData_s MyOpenMMData;
static MyOpenMMData* myInitializeOpenMM(const MyAtomInfo atoms[],
double temperature,
double frictionInPs,
double solventDielectric,
double soluteDielectric,
double stepSizeInFs,
const char** platformName);
static void myStepWithOpenMM(MyOpenMMData*, int numSteps);
static void myGetOpenMMState(MyOpenMMData*, int 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[])
{
int n;
/* Write out in PDB format. */
printf("MODEL %d\n", frameNum);
printf("REMARK 250 time=%.3f ps; energy=%.3f kcal/mole\n",
timeInPs, energyInKcal);
for (n=0; *atoms[n].pdb; ++n)
printf("ATOM %5d %4s SLT 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");
}
/* --------------------------------------------------------------------------
* MAIN PROGRAM
* -------------------------------------------------------------------------- */
int main() {
const int NumReports = (int)(SimulationTimeInPs*1000 / ReportIntervalInFs + 0.5);
const int NumSilentSteps = (int)(ReportIntervalInFs / StepSizeInFs + 0.5);
int frame;
/* TODO: what about thrown exceptions? */
double time, energy;
const char* platformName;
/* Set up OpenMM data structures; returns OpenMM Platform name. */
MyOpenMMData* omm = myInitializeOpenMM(atoms, Temperature, FrictionInPerPs,
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);
myGetOpenMMState(omm, WantEnergy, &time, &energy, atoms);
myWritePDBFrame(0, time, energy, atoms);
for (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. */
}
/* --------------------------------------------------------------------------
* OpenMM-USING CODE
* --------------------------------------------------------------------------
* The OpenMM C-wrapped 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. We suggest that you write them in C++ if possible; in fact you
* can use the implementation from the C++ version of this example if you
* want. However, the methods are reimplemented in C below in case you prefer.
*/
#include "OpenMMCWrapper.h"
struct MyOpenMMData_s {
OpenMM_System* system;
OpenMM_Context* 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 an opaque pointer to 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 frictionInPerPs,
double solventDielectric,
double soluteDielectric,
double stepSizeInFs,
const char** platformName)
{
/* Allocate space to hold OpenMM objects while we're using them. */
MyOpenMMData* omm = (MyOpenMMData*)malloc(sizeof(struct MyOpenMMData_s));
/* These are temporary OpenMM objects used and discarded here. */
OpenMM_Vec3Array* initialPosInNm;
OpenMM_StringArray* pluginList;
OpenMM_NonbondedForce* nonbond;
OpenMM_GBSAOBCForce* gbsa;
OpenMM_Platform* platform;
int n;
/* Load all available OpenMM plugins from their default location. */
pluginList = OpenMM_Platform_loadPluginsFromDirectory
(OpenMM_Platform_getDefaultPluginsDirectory());
OpenMM_StringArray_destroy(pluginList);
/* 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 = OpenMM_System_create();
nonbond = OpenMM_NonbondedForce_create();
gbsa = OpenMM_GBSAOBCForce_create();
OpenMM_System_addForce(omm->system, (OpenMM_Force*)nonbond);
OpenMM_System_addForce(omm->system, (OpenMM_Force*)gbsa);
/* Specify dielectrics for GBSA implicit solvation. */
OpenMM_GBSAOBCForce_setSolventDielectric(gbsa, solventDielectric);
OpenMM_GBSAOBCForce_setSoluteDielectric(gbsa, 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.
* (4) Collect default positions for initializing the simulation later. */
initialPosInNm = OpenMM_Vec3Array_create(0);
for (n=0; *atoms[n].pdb; ++n) {
const MyAtomInfo* atom = &atoms[n];
OpenMM_Vec3 posInNm;
OpenMM_System_addParticle(omm->system, atom->mass);
OpenMM_NonbondedForce_addParticle(nonbond,
atom->charge,
atom->vdwRadiusInAng * OpenMM_NmPerAngstrom
* OpenMM_SigmaPerVdwRadius,
atom->vdwEnergyInKcal * OpenMM_KJPerKcal);
OpenMM_GBSAOBCForce_addParticle(gbsa,
atom->charge,
atom->gbsaRadiusInAng * OpenMM_NmPerAngstrom,
atom->gbsaScaleFactor);
/* Convert the initial position to nm and append to the array. */
posInNm = OpenMM_Vec3_scale(*(const OpenMM_Vec3*)atom->initPosInAng, OpenMM_NmPerAngstrom);
OpenMM_Vec3Array_append(initialPosInNm, 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 = (OpenMM_Integrator*)OpenMM_LangevinIntegrator_create(
temperature, frictionInPerPs,
stepSizeInFs * OpenMM_PsPerFs);
omm->context = OpenMM_Context_create(omm->system, omm->integrator);
OpenMM_Context_setPositions(omm->context, initialPosInNm);
platform = OpenMM_Context_getPlatform(omm->context);
*platformName = OpenMM_Platform_getName(platform);
return omm;
}
/* --------------------------------------------------------------------------
* COPY STATE BACK TO CPU FROM OPENMM
* -------------------------------------------------------------------------- */
static void
myGetOpenMMState(MyOpenMMData* omm, int wantEnergy,
double* timeInPs, double* energyInKcal,
MyAtomInfo atoms[])
{
OpenMM_State* state;
const OpenMM_Vec3Array* posArrayInNm;
int infoMask;
int n;
infoMask = OpenMM_State_Positions;
if (wantEnergy) {
infoMask += OpenMM_State_Velocities; /*for kinetic energy (cheap)*/
infoMask += OpenMM_State_Energy; /*for pot. energy (expensive)*/
}
/* Forces are also available (and cheap). */
/* State object is created here and must be explicitly destroyed below. */
state = OpenMM_Context_getState(omm->context, infoMask, 0);
*timeInPs = OpenMM_State_getTime(state); /* OpenMM time is in ps already. */
/* 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. */
posArrayInNm = OpenMM_State_getPositions(state);
for (n=0; *atoms[n].pdb; ++n)
/* Sets atoms[n].pos = posArray[n] * Angstroms/nm. */
*(OpenMM_Vec3*)atoms[n].posInAng =
OpenMM_Vec3_scale(*OpenMM_Vec3Array_get(posArrayInNm, n),
OpenMM_AngstromsPerNm);
/* If energy has been requested, obtain it and convert from kJ to kcal. */
*energyInKcal = 0;
if (wantEnergy)
*energyInKcal = ( OpenMM_State_getPotentialEnergy(state)
+ OpenMM_State_getKineticEnergy(state))
* OpenMM_KcalPerKJ;
OpenMM_State_destroy(state);
}
// -----------------------------------------------------------------------------
// TAKE MULTIPLE STEPS USING OpenMM
// -----------------------------------------------------------------------------
static void
myStepWithOpenMM(MyOpenMMData* omm, int numSteps) {
OpenMM_Integrator_step(omm->integrator, numSteps);
}
// -----------------------------------------------------------------------------
// DEALLOCATE OpenMM OBJECTS
// -----------------------------------------------------------------------------
static void
myTerminateOpenMM(MyOpenMMData* omm) {
/* Clean up top-level heap allocated objects that we're done with now. */
OpenMM_Context_destroy(omm->context);
OpenMM_Integrator_destroy(omm->integrator);
OpenMM_System_destroy(omm->system);
free(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