Added Hello Argon, made some mods to the other examples.

examples/HelloArgon.cpp

+#include "OpenMM.h"
+#include <cstdio>
+
+using namespace OpenMM;
+
+void writePdb(const OpenMMContext& context) {
+    const State state = context.getState(State::Positions);
+    const std::vector<Vec3>& pos = state.getPositions();
+    static int modelFrameNumber = 0; // numbering for MODEL records in pdb output
+    // write out in PDB format
+    modelFrameNumber++;
+    printf("MODEL     %d\n", modelFrameNumber);
+    for (int a = 0; a < context.getSystem().getNumParticles(); ++a)
+        printf("ATOM  %5d AR    AR     1    %8.3f%8.3f%8.3f  1.00  0.00          AR\n",
+            a+1, pos[a][0]*10, pos[a][1]*10, pos[a][2]*10);
+    printf("ENDMDL\n");
+}
+
+int main() {
+    Platform::loadPluginsFromDirectory(Platform::getDefaultPluginsDirectory());
+    System system;
+    NonbondedForce* nonbond = new NonbondedForce();
+    system.addForce(nonbond);
+
+    // Create atoms
+    int numAtoms = 2;
+    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
+    // Let OpenMM Context choose best platform.
+    OpenMMContext context(system, integrator);
+    printf( "REMARK  Using OpenMM platform %s\n", context.getPlatform().getName().c_str() );
+
+    // Set the starting positions in the OpenMM context. Velocities will be zero.
+    std::vector<Vec3> positions(numAtoms);
+    for (int a = 0; a < numAtoms; ++a)
+        positions[a] = Vec3(0.5*a,0,0);
+    context.setPositions(positions);
+
+    while(context.getTime() < 500.0) { // picoseconds
+        writePdb(context);
+        integrator.step(100); // number of steps between reports
+    }
+    writePdb(context);
+
+    return 0;
+}

examples/HelloSodiumChloride.cpp

@@ -27,20 +27,24 @@ struct AtomInfo {
 };
 
 static AtomInfo atoms[] = {
-    {"Na+", 22.99,  1, 1.8680, 0.00277, Vec3(-3,0,0)},
-    {"Cl-", 35.45, -1, 2.4700, 0.1000,  Vec3( 3,0,0)},
+    {"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         = 100;     // Kelvins
+static const double Temperature         = 300;    // Kelvins
 static const double Friction            = 1./91.; // picoseconds between collisions
-static const double StepSizeFs          = .1;     // femtoseconds
-static const double ReportIntervalFs    = 1000;
-static const double SimulationTimePs    = 1000;  // total simulation time (ps)
+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 showState(const OpenMMContext&);
+static void writePDB(const OpenMMContext&);
 
 int main() {
 try {
@@ -62,14 +66,14 @@ try {
     }
 
     // Create an integrator object for advancing time.
-    //LangevinIntegrator integrator(Temperature, Friction, StepSizeFs * PsPerFs);
-    VerletIntegrator integrator(StepSizeFs * PsPerFs);
+    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();
-    std::cout << "Will use OpenMM platform " << platformName << std::endl;
+    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);
@@ -80,13 +84,13 @@ try {
     context.setPositions(positions);
 
     // Output the initial state.
-    showState(context);
+    writePDB(context);
 
     const int NumSilentSteps = (int)(ReportIntervalFs / StepSizeFs + 0.5);
     do {
         integrator.step(NumSilentSteps);
-        showState(context);
-    } while (context.getTime() <= SimulationTimePs);
+        writePDB(context);
+    } while (context.getTime() < SimulationTimePs);
 
     } catch(const std::exception& e) {
         std::cout << "EXCEPTION: " << e.what() << std::endl;
@@ -97,7 +101,7 @@ try {
 }
 
 static void
-showState(const OpenMMContext& context) {
+writePDB(const 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();
@@ -108,10 +112,11 @@ showState(const OpenMMContext& context) {
     modelFrameNumber++;
     printf("MODEL     %d\n", modelFrameNumber);
     printf("REMARK 250 time=%.3f picoseconds; Energy = %.3f kilojoules/mole\n", state.getTime(), energy);
-    printf("ATOM      1 NA   SLT     1    %8.3f%8.3f%8.3f  1.00  0.00          NA\n", 
-        positions[0][0]*AngstromsPerNm, positions[0][1]*AngstromsPerNm, positions[0][2]*AngstromsPerNm);
-    printf("ATOM      2 CL   SLT     1    %8.3f%8.3f%8.3f  1.00  0.00          CL\n", 
-        positions[1][0]*AngstromsPerNm, positions[1][1]*AngstromsPerNm, positions[1][2]*AngstromsPerNm);
+    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);
+    }
     printf("ENDMDL\n");
 }
 

examples/HelloWorld.cpp

@@ -2,6 +2,7 @@
 
 #include <iostream>
 #include <string>
+#include <limits>
 
 using namespace OpenMM;
 
@@ -33,7 +34,7 @@ static AtomInfo atoms[] = {
 
 static const double Temperature         = 100;     // Kelvins
 static const double Friction            = 1./91.;  // picoseconds between collisions
-static const double StepSizeFs          = 0.1;     // femtoseconds
+static const double StepSizeFs          = 2;     // femtoseconds
 static const double ReportIntervalFs    = 1000;
 static const double SimulationTimePs    = 1000;  // total simulation time (ps)
 
@@ -68,6 +69,7 @@ try {
     OpenMMContext context(system, integrator);
 
     std::cout << "Will use OpenMM platform " << context.getPlatform().getName() << std::endl;
+    std::cout << "eps(float)=" << std::numeric_limits<float>::epsilon() << std::endl;
 
     // Fill in a vector of starting positions, one per atom.
     std::vector<Vec3> positions(numAtoms);
@@ -84,7 +86,7 @@ try {
     do {
         integrator.step(NumSilentSteps);
         showState(context);
-    } while (context.getTime() <= SimulationTimePs);
+    } while (context.getTime() < SimulationTimePs);
 
     } catch(const std::exception& e) {
         std::cout << "EXCEPTION: " << e.what() << std::endl;
@@ -97,11 +99,17 @@ try {
 static void
 showState(const 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 State                 state       = context.getState(  State::Positions | State::Velocities 
+                                                               | State::Forces    | State::Energy);
     const double                energy      = state.getPotentialEnergy() + state.getKineticEnergy();
     const std::vector<Vec3>&    positions   = state.getPositions();
+    const std::vector<Vec3>&    forces      = state.getForces();
 
     std::cout << "t=" << state.getTime() << " E=" << energy << std::endl;
-    for (unsigned i=0; i < positions.size(); ++i)
-        std::cout << i << " " << positions[i] << std::endl;
+    const double dt = StepSizeFs*PsPerFs;
+    for (unsigned i=0; i < positions.size(); ++i) {
+        const Vec3 dq = dt*dt/(2*atoms[i].mass) * forces[i];
+        std::cout << i << " " << positions[i] * AngstromsPerNm << std::endl;
+        std::cout << "f=" << forces[i] << " dq=" << dq << std::endl;
+    }
 }
\ No newline at end of file