Fix: LAMMPS to Gromacs Conversion Issues

.gitignore

@@ -37,3 +37,9 @@ Thumbs.db
 
 # python #
 *.pyc
+
+build/
+dist/
+intermol.egg-info/
+
+

README.md

@@ -1,8 +1,11 @@
+
+Hi!
+This is a privately maintained version of Intermol.
+
+
 InterMol: a conversion tool for molecular dynamics simulations
 ==============================================================
 
-[![Linux Build Status](https://travis-ci.org/shirtsgroup/InterMol.svg?branch=develop)](https://travis-ci.org/shirtsgroup/InterMol)
-[![Coverage Status](https://coveralls.io/repos/shirtsgroup/InterMol/badge.svg?branch=develop)](https://coveralls.io/r/shirtsgroup/InterMol)
 
 We are currently in beta testing phase. Desmond<=>Gromacs<=>Lammps conversions are carried out natively in InterMol.  AMBER->X is carried out by converting AMBER to GROMACS, then to other programs using ParmEd. AMBER->CHARMM is carried out by ParmEd directly.
 

__conda_version__.txt

@@ -0,0 +1 @@
+0.1.0

intermol/forces/__init__.py

@@ -61,7 +61,7 @@
 from intermol.forces.trig_dihedral_type import TrigDihedralType, TrigDihedral
 from intermol.forces.restricted_bending_dihedral_type import RestrictedBendingDihedralType, RestrictedBendingDihedral
 from intermol.forces.bending_torsion_dihedral_type import BendingTorsionDihedralType, BendingTorsionDihedral
-
+from intermol.forces.improper_cvff_dihedral_type import ImproperCvffDihedralType, ImproperCvffDihedral
 
 # virtual_sites
 from intermol.forces.two_virtual_type import TwoVirtualType, TwoVirtual
@@ -80,5 +80,5 @@
 from intermol.forces.convert_dihedrals import convert_dihedral_from_proper_to_trig
 from intermol.forces.convert_dihedrals import convert_dihedral_from_fourier_to_trig
 from intermol.forces.convert_dihedrals import convert_dihedral_from_trig_to_fourier
-
+from intermol.forces.convert_dihedrals import convert_dihedral_from_improper_cvff_to_trig
 

intermol/forces/convert_dihedrals.py

@@ -34,6 +34,27 @@ def convert_dihedral_from_proper_to_trig(p):
     fcs[fk] = k
     return fcs
 
+def convert_dihedral_from_improper_cvff_to_trig(p):
+
+    k = p['k']
+    multiplicity = p['multiplicity']
+    sign = p['sign']
+    zu = 0*k.unit
+    fcs = {
+        'phi': 0 * units.degrees,
+        'fc0': k,
+        'fc1': zu,
+        'fc2': zu,
+        'fc3': zu,
+        'fc4': zu,
+        'fc5': zu,
+        'fc6': zu
+        }
+
+    k # which force constant is nonzero because of the multiplicity?
+    fk = "fc%d" % (multiplicity._value)
+    fcs[fk] = sign * k
+    return fcs
 
 def convert_dihedral_from_fourier_to_trig(f):
 
@@ -212,7 +233,10 @@ def convert_dihedral_from_RB_to_trig(c):
     c2 = c['C2']
     c3 = c['C3']
     c4 = c['C4']
-    c5 = c['C5']
+    if 'C5' in c:  # program might not define this one, need to check it exists.
+        c5 = c['C5']
+    else:
+        c5 = 0*c0.unit
     if 'C6' in c:  # program might not define this one, need to check it exists.
         c6 = c['C6']
     else:

intermol/forces/forcedata.py

@@ -362,6 +362,10 @@
 # ========
 # dihedrals
 # ========
+lammps_paramlist['improper_cvff_dihedral'] = ['k', 'sign', 'multiplicity']
+master_paramlist['improper_cvff_dihedral'] = ['k', 'sign', 'multiplicity']
+master_unitlist['improper_cvff_dihedral']  = ['energy', 'units.dimensionless', 'units.dimensionless']
+
 doclist['improper_harmonic_dihedral'] = 'stub documentation\n'
 master_paramlist['improper_harmonic_dihedral'] = ['xi', 'k']
 master_unitlist['improper_harmonic_dihedral'] = ['angleD',

intermol/forces/forcefunctions.py

@@ -136,7 +136,10 @@ def create_kwds_from_entries(unitvars, paramlist, entries, force_type, offset=0)
     u = unitvars[typename]
     params = paramlist[typename]
     for i, p in enumerate(params):
-        kwds[p] = float(entries[offset+i]) * u[i]
+        if (len(entries)<=(offset+i)):
+            kwds[p] = float(0) * u[i]
+        else:
+            kwds[p] = float(entries[offset+i]) * u[i]
     return kwds
 
 

intermol/forces/improper_cvff_dihedral_type.py

@@ -0,0 +1,40 @@
+import parmed.unit as units
+
+from intermol.decorators import accepts_compatible_units
+from intermol.forces.abstract_dihedral_type import AbstractDihedralType
+
+
+class ImproperCvffDihedralType(AbstractDihedralType):
+    __slots__ = ['k', 'sign', 'multiplicity']
+
+    @accepts_compatible_units(None, None, None, None, 
+                              k=units.kilojoules_per_mole,
+                              sign=units.dimensionless,
+                              multiplicity=units.dimensionless,)
+    def __init__(self, bondingtype1, bondingtype2, bondingtype3, bondingtype4, 
+                 k=0.0 * units.kilojoules_per_mole,
+                 sign=1.0 * units.dimensionless,
+                 multiplicity=0.0 * units.dimensionless,
+                 ):
+        AbstractDihedralType.__init__(self, bondingtype1, bondingtype2, bondingtype3, bondingtype4)
+        self.k = k
+        self.sign = sign
+        self.multiplicity = multiplicity
+
+
+class ImproperCvffDihedral(ImproperCvffDihedralType):
+    """
+    stub documentation
+    """
+    def __init__(self, atom1, atom2, atom3, atom4, bondingtype1=None, bondingtype2=None, bondingtype3=None, bondingtype4=None, 
+                 k=0.0 * units.kilojoules_per_mole,
+                 sign=1.0 * units.dimensionless,
+                 multiplicity=0.0 * units.dimensionless):
+        self.atom1 = atom1
+        self.atom2 = atom2
+        self.atom3 = atom3
+        self.atom4 = atom4
+        ImproperCvffDihedralType.__init__(self, bondingtype1, bondingtype2, bondingtype3, bondingtype4, 
+                k=k,
+                sign=sign,
+                multiplicity=multiplicity)

intermol/gromacs/grofile_parser.py

@@ -102,10 +102,10 @@ def write(self, system):
                 gro.write('{0:5d}{1:<5s}{2:5s}{3:5d}'.format(
                         atom.residue_index%100000, atom.residue_name, atom.name, (n + 1)%100000))
                 for pos in atom.position:
-                    gro.write('{0:17.12f}'.format(pos.value_in_unit(nanometers)))
+                    gro.write('{0:8.3f}'.format(pos.value_in_unit(nanometers)))
                 if np.any(atom.velocity):
                     for vel in atom.velocity:
-                        gro.write('{0:17.12f}'.format(vel.value_in_unit(nanometers / picoseconds)))
+                        gro.write('{0:8.4f}'.format(vel.value_in_unit(nanometers / picoseconds)))
                 gro.write('\n')
 
             # Check for rectangular; should be symmetric, so we don't have to

intermol/lammps/lammps_parser.py

@@ -73,26 +73,25 @@ class LammpsParser(object):
     lammps_bond_types = dict(
         (k, eval(v.__name__ + 'Type')) for k, v in lammps_bonds.items())
 
-    def canonical_bond(self, params, bond, direction='into'):
+    def canonical_bond(self, params, bond, direction):
         """Convert to/from the canonical form of this interaction. """
         # TODO: Gromacs says harmonic potential bonds do not have constraints or
         #       exclusions. Check that this logic is supported.
         if direction == 'into':
             canonical_force_scale = self.SCALE_INTO
+            bondtest = bond
         else:
             try:
                 typename = self.lookup_lammps_bonds[bond.__class__]
             except KeyError:
-                if bond.__class__.__name__ in ['FeneBond', 'ConnectionBond']:
-                    raise UnimplementedFunctional(bond, ENGINE)
-                else:
-                    raise UnsupportedFunctional(bond, ENGINE)
+                raise UnsupportedFunctional(bond, ENGINE)
+            bondtest = bond.__class__
             canonical_force_scale = self.SCALE_FROM
 
-        if bond.__class__ in [HarmonicBond, HarmonicPotentialBond]:
+        if bondtest in [HarmonicBond, HarmonicPotentialBond]:
             params['k'] *= canonical_force_scale
 
-        if bond.__class__ == HarmonicPotentialBond:
+        if bondtest == HarmonicPotentialBond:
             typename = 'harmonic'
 
         if direction == 'into':
@@ -153,7 +152,7 @@ def canonical_angle(self, params, angle, direction):
 
     lammps_impropers = {
         'harmonic': ImproperHarmonicDihedral,
-        'cvff': TrigDihedral,
+        'cvff': ImproperCvffDihedral,
         }
     lookup_lammps_impropers = dict((v, k) for k, v in lammps_impropers.items())
     lammps_improper_types = dict(
@@ -169,18 +168,22 @@ def canonical_dihedral(self, params, dihedral, direction='into'):
                 converted_dihedral = TrigDihedral
             elif dihedral == ImproperHarmonicDihedral:
                 convertfunc = convert_nothing
+                converted_dihedral = ImproperHarmonicDihedral
             elif dihedral == RbDihedral:
                 convertfunc = convert_dihedral_from_RB_to_trig
                 converted_dihedral = TrigDihedral
             elif dihedral == FourierDihedral:
                 convertfunc = convert_dihedral_from_fourier_to_trig
                 converted_dihedral = TrigDihedral
+            elif dihedral == ImproperCvffDihedral:
+                convertfunc = convert_dihedral_from_improper_cvff_to_trig
+                converted_dihedral = TrigDihedral
                 # Now actually convert the dihedral.
             params = convertfunc(params)
 
             # Adjust scaling conventions.
             canonical_force_scale = self.SCALE_INTO
-            if converted_dihedral == ImproperHarmonicDihedralType:
+            if converted_dihedral == ImproperHarmonicDihedral:
                 params['k'] *= canonical_force_scale
 
             return converted_dihedral, params
@@ -706,7 +709,6 @@ def parse_force_coeffs(self, data_lines, force_name, force_classes,
             force_classes[int(fields[0])] = [force_class, kwds]
 
     def parse_bond_coeffs(self, data_lines):
-
         self.bond_classes = dict()
         self.parse_force_coeffs(data_lines, "Bond",
                                 self.bond_classes, self.bond_style,