HelloSodiumChlorideInC.c

/* -----------------------------------------------------------------------------
 *            OpenMM(tm) HelloSodiumChloride example in C (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.
 *
 * This example is written entirely in ANSI C, using a set of wrappers which
 * are NOT an official part of OpenMM.
 * -------------------------------------------------------------------------- */
#include "wrappers/OpenMM_CWrapper.h"

#include <stdio.h>

/* -----------------------------------------------------------------------------
 *                   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_Context*); // PDB file writer; see below.

/* -----------------------------------------------------------------------------
 *                                MAIN PROGRAM
 * ----------------------------------------------------------------------------- */
int main() {
    /* TODO: what about thrown exceptions? */
    OpenMM_Platform_loadPluginsFromDirectory
       (OpenMM_Platform_getDefaultPluginsDirectory());

    simulateNaCl();

    return 0;
}

/* -----------------------------------------------------------------------------
 *                          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;
    OpenMM_Vec3 initPosInAngstroms;
} atoms[] = {
    /* pdb   mass   charge   vdwRadius  vdwEnergy   initPos */
    {" NA ", 22.99,   1,     1.8680,    0.00277,    { 8, 0,  0}},
    {" CL ", 35.45,  -1,     2.4700,    0.1000,     {-8, 0,  0}},
    {" NA ", 22.99,   1,     1.8680,    0.00277,    { 0, 9,  0}},
    {" CL ", 35.45,  -1,     2.4700,    0.1000,     { 0,-9,  0}},
    {" NA ", 22.99,   1,     1.8680,    0.00277,    { 0, 0,-10}},
    {" CL ", 35.45,  -1,     2.4700,    0.1000,     { 0, 0, 10}},
    {""} /* end of list */
};

/* -----------------------------------------------------------------------------
 *                              NaCl SIMULATION
 * ----------------------------------------------------------------------------- */
void simulateNaCl() {
    OpenMM_Integrator*  integrator;
    OpenMM_Context*     context;
    OpenMM_Vec3Array*   initialPositionsInNm;
    const OpenMM_Vec3   a = {5,0,0}, b = {0,5,0}, c = {0,0,5};
    const int           NumSilentSteps = (int)(ReportIntervalInFs / StepSizeInFs + 0.5);
    int n;
    /* -------------------------------------------------------------------------
     * Create System and Force object. Add the Force to the System. The System
     * takes over ownership, so you should not destroy the Force object yourself.
     * ------------------------------------------------------------------------- */
    OpenMM_System*           system  = OpenMM_System_create();
    OpenMM_NonbondedForce*   nonbond = OpenMM_NonbondedForce_create();
    OpenMM_System_addForce(system, nonbond);

    /* -------------------------------------------------------------------------
     * Specify a periodic box and cutoff distance.
     * ------------------------------------------------------------------------- */
    OpenMM_NonbondedForce_setNonbondedMethod(nonbond, OpenMM_NonbondedForce_CutoffPeriodic);
    OpenMM_NonbondedForce_setCutoffDistance(nonbond, 2);
    OpenMM_NonbondedForce_setPeriodicBoxVectors(nonbond, a, b, c);

    /* -------------------------------------------------------------------------
     * 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.
     * ------------------------------------------------------------------------- */
    initialPositionsInNm = OpenMM_Vec3Array_create(0);
    for (n=0; *atoms[n].pdb; ++n) {
        const struct AtomInfo* atom = &atoms[n];
        OpenMM_Vec3 posInNm;

        OpenMM_System_addParticle(system, atom->mass);
        OpenMM_NonbondedForce_addParticle
           (nonbond, atom->charge,
                     atom->vdwRadiusInAng  * OpenMM_NmPerAngstrom 
                                           * OpenMM_SigmaPerVdwRadius,
                     atom->vdwEnergyInKcal * OpenMM_KJPerKcal);

        /* Convert the initial position to nm and append to the array. */
        OpenMM_Vec3_scale(atom->initPosInAngstroms, OpenMM_NmPerAngstrom, posInNm);
        OpenMM_Vec3Array_append(initialPositionsInNm, 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.
     * ------------------------------------------------------------------------- */
    integrator = (OpenMM_Integrator*)OpenMM_VerletIntegrator_create(StepSizeInFs * OpenMM_PsPerFs);
    context    = OpenMM_Context_create(system, integrator);
    OpenMM_Context_setPositions(context, initialPositionsInNm);

    /* -------------------------------------------------------------------------
     * 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", OpenMM_Context_getPlatform(context));
    writePDB(context);

    do {
        OpenMM_Integrator_step(integrator, NumSilentSteps);
        writePDB(context);
    } while (OpenMM_Context_getTime(context) < SimulationTimeInPs);

    /* Clean up top-level heap allocated objects that we're done with now. */
    OpenMM_Vec3Array_destroy(initialPositionsInNm);
    OpenMM_Context_destroy(context);
    OpenMM_Integrator_destroy(integrator);
}

/* -----------------------------------------------------------------------------
 *                               PDB FILE WRITER
 * ----------------------------------------------------------------------------- */
static void
writePDB(const OpenMM_Context* context) {
    static int modelFrameNumber = 0; /*numbering for MODEL records in pdb output*/

    /* Caution: at the moment asking for energy requires use of slow Reference 
       platform calculation. */

    /* Don't forget to destroy this State when you're done with it. */
    OpenMM_State* state = OpenMM_Context_createState(context,
        OpenMM_State_Positions | OpenMM_State_Velocities | OpenMM_State_Energy);

    const double energy =   OpenMM_State_getPotentialEnergy(state) 
                          + OpenMM_State_getKineticEnergy(state);

    /* 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.
     */
    const OpenMM_Vec3Array* posArray  = OpenMM_State_getPositions(state);
    const OpenMM_Vec3*      positions = OpenMM_Vec3Array_getAsVec3Ptr(posArray);
    const int               npos      = OpenMM_Vec3Array_size(posArray); 
    int i;

    /* write out in PDB format */
    modelFrameNumber++;
    printf("MODEL     %d\n", modelFrameNumber);
    printf("REMARK 250 time=%.3f picoseconds; Energy = %.3f kilojoules/mole\n", 
           OpenMM_State_getTime(state), energy);

    for (i=0; i < npos; ++i) {
        OpenMM_Vec3 pos;
        OpenMM_Vec3_scale(positions[i], OpenMM_AngstromsPerNm, pos);
        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_State_destroy(state);
}