Refactor HelloArgon toward better first example program

examples/HelloArgon.cpp

@@ -3,49 +3,79 @@
 
 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");
-}
+// forward declaration of subroutine defined later in this file.
+void writePdb(const OpenMMContext& context);
+
+void simulateArgon() 
+{
+    // Load any shared libraries containing GPU implementations
+    Platform::loadPluginsFromDirectory(
+        Platform::getDefaultPluginsDirectory());
 
-int main() {
-    Platform::loadPluginsFromDirectory(Platform::getDefaultPluginsDirectory());
     System system;
     NonbondedForce* nonbond = new NonbondedForce(); 
     system.addForce(nonbond);
 
-    // Create atoms
-    int numAtoms = 3;
-    for (int a = 0; a < numAtoms; ++a) {
+    // Create three atoms
+    std::vector<Vec3> initialPositions(3);
+    for (int a = 0; a < 3; ++a) 
+    {
         system.addParticle(39.95); // mass
-        nonbond->addParticle(0.0, 0.3350, 0.001603); // charge, diameter, well depth
+        nonbond->addParticle(0.0, 0.3350, 0.001603); // charge, sigma, well depth
+        initialPositions[a] = Vec3(0.5*a,0,0); // location
     }
 
-    // Large step size may be required for stability with small forces
     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() );
+    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);
+    // Set the starting positions of the atoms. Velocities will be zero.
+    context.setPositions(initialPositions);
 
+    // Simulate
     while(context.getTime() < 500.0) { // picoseconds
-        writePdb(context);
-        integrator.step(100); // number of steps between reports
+        writePdb(context); // output coordinates
+        // Run 100 steps at a time, for efficient use of OpenMM
+        integrator.step(100);
+    }
+    writePdb(context); // output final coordinates
+}
+
+int main() 
+{
+    try {
+        simulateArgon();
+        return 0; // success!
     }
-    writePdb(context);
+    // 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; // failure!
+    }
+}
+
+// writePdb() subroutine for quick-and-dirty trajectory output.
+void writePdb(const OpenMMContext& context) 
+{
+    // Request atomic positions from OpenMM
+    const State state = context.getState(State::Positions);
+    const std::vector<Vec3>& pos = state.getPositions();
 
-    return 0;
+    // write out in PDB format
+
+    // Use PDB MODEL cards to number trajectory frames
+    static int modelFrameNumber = 0; 
+    modelFrameNumber++;
+    printf("MODEL     %d\n", modelFrameNumber); // start of frame
+    for (int a = 0; a < context.getSystem().getNumParticles(); ++a)
+    {
+        printf("ATOM  %5d AR    AR     1    ", a+1); // atom number
+        printf("%8.3f%8.3f%8.3f  1.00  0.00          AR\n",
+            // notice "*10" converts nanometers to Angstroms
+            pos[a][0]*10, pos[a][1]*10, pos[a][2]*10);
+    }
+    printf("ENDMDL\n"); // end of frame
 }