Added HelloEthane.

examples/HelloArgon.cpp

@@ -19,18 +19,18 @@ void writePdb(const OpenMMContext& context) {
 int main() {
     Platform::loadPluginsFromDirectory(Platform::getDefaultPluginsDirectory());
     System system;
-    NonbondedForce* nonbond = new NonbondedForce();
+    NonbondedForce* nonbond = new NonbondedForce(); 
     system.addForce(nonbond);
 
     // Create atoms
-    int numAtoms = 2;
+    int numAtoms = 3;
     for (int a = 0; a < numAtoms; ++a) {
         system.addParticle(39.95); // mass
         nonbond->addParticle(0.0, 0.3350, 0.001603); // charge, diameter, well depth
     }
 
     // Large step size may be required for stability with small forces
-    VerletIntegrator integrator(0.002); // step size in picoseconds
+    VerletIntegrator integrator(0.020); // step size in picoseconds
     // Let OpenMM Context choose best platform.
     OpenMMContext context(system, integrator);
     printf( "REMARK  Using OpenMM platform %s\n", context.getPlatform().getName().c_str() );

examples/HelloEthane.cpp

+/* -----------------------------------------------------------------------------
+ *                    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) 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;
+
+// These are missing from the current version of OpenMM so we're adding them
+// temporarily here.
+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]);}
+static const double SigmaPerVdwRadius = 2*std::pow(2., -1./6.);
+
+// -----------------------------------------------------------------------------
+//              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* pdbSymbol; Vec3 initPosInAngstroms;
+} atoms[] = {/*0*/C, "C1", Vec3( -.765,   0,   0 ),
+             /*1*/C, "C2", Vec3(  .765,   0,   0 ),
+             /*2*/H, "H1", Vec3(-1.135, 1.03,  0 ), // bonded to C1
+             /*3*/H, "H2", Vec3(-1.135, -.51, .89),
+             /*4*/H, "H3", Vec3(-1.135, -.51,-.89),
+             /*5*/H, "H4", Vec3( 1.135, 1.03,  0 ), // bonded to C2
+             /*6*/H, "H5", Vec3( 1.135, -.51, .89),
+             /*7*/H, "H6", Vec3( 1.135, -.51,-.89),
+             EndOfList};
+
+struct {int type; int a[2];} bonds[]    = {CC,0,1,CH,0,2,CH,0,3,CH,0,4,
+                                                  CH,1,5,CH,1,6,CH,1,7, 
+                                           EndOfList};
+struct {int type; int a[3];} angles[]   = {HCC,2,0,1,HCC,3,0,1,HCC,4,0,1,
+                                           HCC,5,1,0,HCC,6,1,0,HCC,7,1,0,
+                                           HCH,2,0,3,HCH,2,0,4,HCH,3,0,4,
+                                           HCH,5,1,6,HCH,5,1,7,HCH,6,1,7, 
+                                           EndOfList};
+struct {int type; int a[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};
+
+// -----------------------------------------------------------------------------
+//              MODELING AND SIMULATION PARAMETERS
+// -----------------------------------------------------------------------------
+const bool   UseConstraints      = false;   // Should we constrain C-H bonds?
+const double Temperature         = 300;     // bath temperature in Kelvins
+const double FrictionInPs        = 1./91.;  // picoseconds between collisions
+
+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)
+
+// PDB file writer; see below.
+static void writePDB(const OpenMMContext&);
+
+// -----------------------------------------------------------------------------
+//                                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());
+
+    // -------------------------------------------------------------------------
+    // Create a System and Force objects within the System. Retain a reference
+    // to each force object so we can fill in the forces.
+    // -------------------------------------------------------------------------
+    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*      atoms = bonds[i].a;
+        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(atoms[0], atoms[1],
+                            bond.nominalLengthInAngstroms    * NmPerAngstrom,
+                            bond.stiffnessInKcalPerAngstrom2 * 2 * KJPerKcal * AngstromsPerNm * AngstromsPerNm);
+        if (UseConstraints && bond.canConstrain)
+            system.addConstraint(atoms[0], atoms[1],
+                                 bond.nominalLengthInAngstroms * NmPerAngstrom);
+        bondPairs.push_back(std::make_pair(atoms[0], atoms[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*       atoms = angles[i].a;
+        const AngleType& angle = angleType[angles[i].type];
+
+        // See note under bond stretch above regarding the factor of 2 here.
+        bondBend.addAngle(atoms[0],atoms[1],atoms[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*         atoms = torsions[i].a;
+        const TorsionType& torsion = torsionType[torsions[i].type];
+        bondTorsion.addTorsion(atoms[0],atoms[1],atoms[2],atoms[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.
+    // -------------------------------------------------------------------------
+    //LangevinIntegrator integrator(Temperature, FrictionInPs, StepSizeInFs * PsPerFs);
+    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);
+
+    // -------------------------------------------------------------------------
+    // Normal return from main.
+    // -------------------------------------------------------------------------
+    return 0;
+
+  // 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;
+  }
+}
+
+
+// -----------------------------------------------------------------------------
+//                               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 %2s   ETH     1    %8.3f%8.3f%8.3f  1.00  0.00          %2s\n", 
+            i+1, atoms[i].pdbSymbol, pos[0], pos[1], pos[2], atoms[i].pdbSymbol);
+    }
+    printf("ENDMDL\n");
+}
+

examples/HelloSodiumChloride.cpp

@@ -47,7 +47,7 @@ static const double SigmaPerVdwRadius = 2*std::pow(2., -1./6.);
 static void writePDB(const OpenMMContext&);
 
 int main() {
-try {
+  try {
     // Load all available OpenMM plugins from their default location.
     Platform::loadPluginsFromDirectory(Platform::getDefaultPluginsDirectory());
 
@@ -92,12 +92,13 @@ try {
         writePDB(context);
     } while (context.getTime() < SimulationTimePs);
 
-    } catch(const std::exception& e) {
-        std::cout << "EXCEPTION: " << e.what() << std::endl;
-        return 1;
-    }
-
     return 0;
+
+  } catch(const std::exception& e) {
+    std::cout << "EXCEPTION: " << e.what() << std::endl;
+    return 1;
+  }
+
 }
 
 static void