Modified the OpenMM wrapper interface to support the strange CUDA platform implementa...

The CUDA platform in OpenMM does not allow the user to generate more
than one instance of a force kernel object. That requires to pack
force field terms that have the same form into one instance,
e.g. Urey-Bradley and Bond/Angle, and the harmonic dihedrals.

diff --git a/src/kernel/openmm_wrapper.cpp b/src/kernel/openmm_wrapper.cpp
index d043f0756f702b46b4ed653203c54cca55d98bd1..3ed2c01360ecc164ebb0031f149b43a3b92269fd 100644 (file)
--- a/src/kernel/openmm_wrapper.cpp
+++ b/src/kernel/openmm_wrapper.cpp
@@ -845,7 +845,7 @@ void* openmm_init(FILE *fplog, const char *platformOptStr,
         /* check wheter Gromacs options compatibility with OpenMM */
         checkGmxOptions(fplog, opt, ir, top, fr, state);
 
-        // Create the system.
+        /* Create the system. */
         const t_idef& idef = top->idef;
         const int numAtoms = top_global->natoms;
         const int numConstraints = idef.il[F_CONSTR].nr/3;
@@ -863,6 +863,13 @@ void* openmm_init(FILE *fplog, const char *platformOptStr,
             sys->addForce(new CMMotionRemover(ir->nstcomm));
 
         /* Set bonded force field terms. */
+
+               /* 
+                * CUDA platform currently doesn't support more than one
+                * instance of a force object, so we pack all forces that
+                * use the same form into one.
+               */
+
         const int* bondAtoms = (int*) idef.il[F_BONDS].iatoms;
         HarmonicBondForce* bondForce = new HarmonicBondForce();
         sys->addForce(bondForce);
@@ -876,25 +883,6 @@ void* openmm_init(FILE *fplog, const char *platformOptStr,
                                idef.iparams[type].harmonic.rA, idef.iparams[type].harmonic.krA);
         }
 
-        /* Urey-Bradley includes both the angle and bond potential for 1-3 interactions */
-        const int* ubAtoms = (int*) idef.il[F_UREY_BRADLEY].iatoms;
-        HarmonicBondForce* ubBondForce = new HarmonicBondForce();
-        HarmonicAngleForce* ubAngleForce = new HarmonicAngleForce();
-        sys->addForce(ubBondForce);
-        sys->addForce(ubAngleForce);
-        offset = 0;
-        for (int i = 0; i < numUB; ++i)
-        {
-            int type = ubAtoms[offset++];
-            int atom1 = ubAtoms[offset++];
-            int atom2 = ubAtoms[offset++];
-            int atom3 = ubAtoms[offset++];
-            ubBondForce->addBond(atom1, atom3,
-                               idef.iparams[type].u_b.r13, idef.iparams[type].u_b.kUB);
-            ubAngleForce->addAngle(atom1, atom2, atom3, 
-                    idef.iparams[type].u_b.theta*M_PI/180.0, idef.iparams[type].u_b.ktheta);
-        }
-
                /* Set the angle force field terms */
         const int* angleAtoms = (int*) idef.il[F_ANGLES].iatoms;
         HarmonicAngleForce* angleForce = new HarmonicAngleForce();
@@ -910,6 +898,27 @@ void* openmm_init(FILE *fplog, const char *platformOptStr,
                     idef.iparams[type].harmonic.rA*M_PI/180.0, idef.iparams[type].harmonic.krA);
         }
 
+        /* Urey-Bradley includes both the angle and bond potential for 1-3 interactions */
+        const int* ubAtoms = (int*) idef.il[F_UREY_BRADLEY].iatoms;
+               /* HarmonicBondForce* ubBondForce = new HarmonicBondForce(); */
+               /*  HarmonicAngleForce* ubAngleForce = new HarmonicAngleForce(); */
+        /* sys->addForce(ubBondForce); */
+        /* sys->addForce(ubAngleForce); */
+        offset = 0;
+        for (int i = 0; i < numUB; ++i)
+        {
+            int type = ubAtoms[offset++];
+            int atom1 = ubAtoms[offset++];
+            int atom2 = ubAtoms[offset++];
+            int atom3 = ubAtoms[offset++];
+            /* ubBondForce->addBond(atom1, atom3, */
+            bondForce->addBond(atom1, atom3,
+                               idef.iparams[type].u_b.r13, idef.iparams[type].u_b.kUB);
+            /* ubAngleForce->addAngle(atom1, atom2, atom3, */ 
+            angleForce->addAngle(atom1, atom2, atom3, 
+                    idef.iparams[type].u_b.theta*M_PI/180.0, idef.iparams[type].u_b.ktheta);
+        }
+
                /* Set proper dihedral terms */
         const int* periodicAtoms = (int*) idef.il[F_PDIHS].iatoms;
         PeriodicTorsionForce* periodicForce = new PeriodicTorsionForce();
@@ -930,8 +939,8 @@ void* openmm_init(FILE *fplog, const char *platformOptStr,
 
                /* Set improper dihedral terms that are represented by a periodic function (as in AMBER FF) */
         const int* periodicImproperAtoms = (int*) idef.il[F_PIDIHS].iatoms;
-        PeriodicTorsionForce* periodicImproperForce = new PeriodicTorsionForce();
-        sys->addForce(periodicImproperForce);
+        /* PeriodicTorsionForce* periodicImproperForce = new PeriodicTorsionForce(); */
+        /* sys->addForce(periodicImproperForce); */
         offset = 0;
         for (int i = 0; i < numPeriodicImproper; ++i)
         {
@@ -940,7 +949,8 @@ void* openmm_init(FILE *fplog, const char *platformOptStr,
             int atom2 = periodicImproperAtoms[offset++];
             int atom3 = periodicImproperAtoms[offset++];
             int atom4 = periodicImproperAtoms[offset++];
-            periodicImproperForce->addTorsion(atom1, atom2, atom3, atom4,
+            /* periodicImproperForce->addTorsion(atom1, atom2, atom3, atom4, */
+            periodicForce->addTorsion(atom1, atom2, atom3, atom4,
                                       idef.iparams[type].pdihs.mult,
                                       idef.iparams[type].pdihs.phiA*M_PI/180.0,
                                       idef.iparams[type].pdihs.cpA);