parmed/amber/_amberparm.py

@@ -54,6 +54,8 @@ from warnings import warn
 
 # ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
+BE_DIFFICULT = False
+
 class AmberParm(AmberFormat, Structure):
     """
     Amber Topology (parm7 format) class. Gives low, and some high, level access
@@ -310,19 +312,45 @@ class AmberParm(AmberFormat, Structure):
         if struct.unknown_functional:
             raise TypeError('Cannot instantiate an AmberParm from unknown '
                             'functional')
-        if (struct.urey_bradleys or struct.impropers or struct.rb_torsions or
-                struct.cmaps or struct.trigonal_angles or struct.pi_torsions or
-                struct.out_of_plane_bends or struct.stretch_bends or
-                struct.torsion_torsions or struct.multipole_frames):
-            if (struct.rb_torsions or struct.trigonal_angles or
-                    struct.pi_torsions or struct.out_of_plane_bends or
-                    struct.torsion_torsions or struct.multipole_frames or
-                    struct.stretch_bends):
-                raise TypeError('AmberParm does not support all of the '
-                                'parameters defined in the input Structure')
-            # Maybe it just has CHARMM parameters?
-            raise TypeError('AmberParm does not support all of the parameters '
-                            'defined in the input Structure. Try ChamberParm')
+        UNSUPPORTED_ATTRIBUTES = (
+            'urey_bradleys',
+            'impropers',
+            'rb_torsions',
+            'cmaps',
+            'trigonal_angles',
+            'pi_torsions',
+            'out_of_plane_bends',
+            'stretch_bends',
+            'torsion_torsions',
+            'multipole_frames',
+        )
+
+        UNSUPPORTED_AMBER_ATTRIBUTES = (
+            'rb_torsions',
+            'trigonal_angles',
+            'pi_torsions',
+            'out_of_plane_bends',
+            'torsion_torsions',
+            'multipole_frames',
+            'stretch_bends',
+        )
+
+        if BE_DIFFICULT:
+            unsupported_features = [struct.__dict__[attribute] for attribute in UNSUPPORTED_ATTRIBUTES]
+            if any(unsupported_features):
+                unsupported_amber_features = [struct.__dict__[attribute] for attribute in UNSUPPORTED_AMBER_ATTRIBUTES]
+                if any(unsupported_amber_features):
+                    raise TypeError(
+                        'AmberParm does not support all of the parameters defined in the input Structure: {0}'.format(
+                            [attribute for (attribute, value) in zip(UNSUPPORTED_AMBER_ATTRIBUTES, unsupported_amber_features) if value],
+                        ),
+                    )
+                # Maybe it just has CHARMM parameters?
+                raise TypeError(
+                    'AmberParm does not support all of the parameters defined in the input Structure: {0}. Try ChamberParm'.format(
+                        [attribute for (attribute, value) in zip(UNSUPPORTED_ATTRIBUTES, unsupported_features) if value],
+                    ),
+                )
         inst = struct.copy(cls, split_dihedrals=True)
         inst.update_dihedral_exclusions()
         inst._add_missing_13_14()

parmed/amber/_chamberparm.py

@@ -182,7 +182,8 @@ class ChamberParm(AmberParm):
         nbfixes = inst.atoms.assign_nbidx_from_types()
         # Give virtual sites a name that Amber understands
         for atom in inst.atoms:
-            if isinstance(atom, ExtraPoint): atom.type = 'EP'
+            pass
+            #if isinstance(atom, ExtraPoint): atom.type = 'EP'
         # Fill the Lennard-Jones arrays/dicts
         ntyp = 0
         for atom in inst.atoms:

parmed/gromacs/gromacstop.py

@@ -364,16 +364,11 @@ class GromacsTopologyFile(Structure):
                     attype, typ = self._parse_atomtypes(line)
                     params.atom_types[attype] = typ
                 elif current_section == 'nonbond_params':
-                    words = line.split()
-                    a1, a2 = words[:2]
-#                   func = int(words[2]) #... unused
-                    sig, eps = (float(x) for x in words[3:5])
-                    sig *= 10 # Convert to Angstroms
-                    eps *= u.kilojoule.conversion_factor_to(u.kilocalorie)
-                    params.nbfix_types[(a1, a2)] = (eps, sig*2**(1/6))
-                    params.nbfix_types[(a2, a1)] = (eps, sig*2**(1/6))
-                    params.atom_types[a1].add_nbfix(a2, sig*2**(1/6), eps)
-                    params.atom_types[a2].add_nbfix(a1, sig*2**(1/6), eps)
+                    eps, rmin, a1, a2 = eps_rmin_a1_a2_for_line(line)
+                    params.nbfix_types[(a1, a2)] = (eps, rmin)
+                    params.nbfix_types[(a2, a1)] = (eps, rmin)
+                    params.atom_types[a1].add_nbfix(a2, rmin, eps)
+                    params.atom_types[a2].add_nbfix(a1, rmin, eps)
                 elif current_section == 'bondtypes':
                     a, b, t = self._parse_bondtypes(line)
                     params.bond_types[(a, b)] = t
@@ -479,21 +474,29 @@ class GromacsTopologyFile(Structure):
         words = line.split()
         i, j = int(words[0])-1, int(words[1])-1
         funct = int(words[2])
-        if funct != 1:
-            warnings.warn('bond funct != 1; unknown functional',
-                          GromacsWarning)
-            self.unknown_functional = True
         bond = Bond(atoms[i], atoms[j])
         bond.funct = funct
         bond_type = None
-        if len(words) >= 5 and funct == 1:
+        if len(words) >= 5:
             req, k = (float(x) for x in words[3:5])
             if (req, k) in bond_types:
                 bond.type = bond_types[(req, k)]
             else:
+                if funct == 1:
+                    k_harmonic = k / 2
+                elif funct == 2:
+                    # Convert quartic to harmonic
+                    bond.funct = 1
+                    k_harmonic = k * req**2
+                else:
+                    warnings.warn(
+                        'bond funct != 1 or 2; unknown functional',
+                        GromacsWarning,
+                    )
+                    self.unknown_functional = True
                 bond_type = BondType(
-                        k*u.kilojoule_per_mole/u.nanometer**2/2,
-                        req*u.nanometer
+                        k_harmonic * u.kilojoule_per_mole / u.nanometer**2,
+                        req*u.nanometer,
                 )
                 bond_types[(req, k)] = bond.type = bond_type
         return bond, bond_type
@@ -530,22 +533,32 @@ class GromacsTopologyFile(Structure):
         words = line.split()
         i, j, k = [int(w)-1 for w in words[:3]]
         funct = int(words[3])
-        if funct not in (1, 5):
-            warnings.warn('angles funct != 1 or 5; unknown '
-                          'functional', GromacsWarning)
+        if funct not in (1, 2, 5):
+            warnings.warn(
+                'angles funct != 1 or 2 or 5; unknown functional',
+                GromacsWarning,
+            )
             self.unknown_functional = True
         angt = ub = ubt = None
         ang = Angle(atoms[i], atoms[j], atoms[k])
         ang.funct = funct
         if funct == 5:
             ub = UreyBradley(atoms[i], atoms[k])
-        if (funct == 1 and len(words) >= 6) or (funct == 5 and len(words) >= 8):
+        if (funct == 1 and len(words) >= 6) or (funct == 5 and len(words) >= 8) or (funct == 2):
             theteq, k = (float(x) for x in words[4:6])
             if (theteq, k) in angle_types:
                 ang.type = angle_types[(theteq, k)]
             else:
-                angt = AngleType(k*u.kilojoule_per_mole/u.radian**2/2,
-                                 theteq*u.degree)
+                if funct == 1:
+                    k_harmonic = k / 2
+                elif funct == 2:
+                    # Convert quartic to harmonic
+                    ang.funct == 1
+                    k_harmonic = k * math.sin(theteq)**2
+                angt = AngleType(
+                    k_harmonic * u.kilojoule_per_mole/u.radian**2,
+                    theteq*u.degree,
+                )
                 angle_types[(theteq, k)] = ang.type = angt
         if funct == 5 and len(words) >= 8:
             ubreq, ubk = (float(x) for x in words[6:8])
@@ -748,27 +761,35 @@ class GromacsTopologyFile(Structure):
     def _parse_bondtypes(self, line):
         """ Parse bondtypes line. Returns str, str, BondType """
         words = line.split()
-        r = float(words[3]) * u.nanometers
-        k = (float(words[4]) / 2) * (u.kilojoules_per_mole / u.nanometers**2)
-        if words[2] != '1':
-            warnings.warn('bondtypes funct != 1; unknown functional',
+        r = float(words[3])
+        funct = int(words[2])
+        if funct == 1:
+            k_harmonic = (float(words[4]) / 2)
+        elif funct == 2:
+            k_harmonic = float(words[4]) * r**2
+        else:
+            warnings.warn('bondtypes funct != 1 or 2; unknown functional',
                           GromacsWarning)
             self.unknown_functional = True
-        return words[0], words[1], BondType(k, r)
+        return words[0], words[1], BondType(k_harmonic * (u.kilojoules_per_mole / u.nanometers**2), r * u.nanometers)
 
     def _parse_angletypes(self, line):
         """
         Parses angletypes line. Returns str, str, str, AngleType, BondType/None
         """
         words = line.split()
-        theta = float(words[4]) * u.degrees
-        k = (float(words[5]) / 2) * (u.kilojoules_per_mole / u.radians**2)
-        if words[3] != '1' and words[3] != '5':
-            warnings.warn('angletypes funct != 1 or 5; unknown functional',
+        theta = float(words[4])
+        funct = int(words[3])
+        if funct == 1:
+            k_harmonic = (float(words[5]) / 2)
+        elif funct == 2:
+            k_harmonic = float(words[5]) * math.sin(theta)**2
+        if funct not in (1, 2, 5):
+            warnings.warn('angletypes funct != 1 or 2 or 5; unknown functional',
                           GromacsWarning)
             self.unknown_functional = True
         ub = None
-        if words[3] == '5':
+        if funct == '5':
             # Contains the angle with urey-bradley
             ub0 = float(words[6])
             cub = float(words[7]) / 2
@@ -778,7 +799,7 @@ class GromacsTopologyFile(Structure):
                 ub0 *= u.nanometers
                 cub *= u.kilojoules_per_mole / u.nanometers**2
                 ub = BondType(cub, ub0)
-        return words[0], words[1], words[2], AngleType(k, theta), ub
+        return words[0], words[1], words[2], AngleType(k_harmonic * (u.kilojoules_per_mole / u.radians**2), theta * u.degrees), ub
 
     def _parse_dihedraltypes(self, line):
         """ Parse dihedraltypes, returns (str,str,str,str), str, Type, bool """
@@ -842,13 +863,8 @@ class GromacsTopologyFile(Structure):
         return a1, a2, a3, a4, a5, CmapType(res1, grid)
 
     def _parse_pairtypes(self, line):
-        words = line.split()
-        a1, a2 = words[:2]
-#       funct = int(words[2]) # ... unused
-        cs6, cs12 = (float(x) for x in words[3:5])
-        cs6 *= u.nanometers * 2**(1/6)
-        cs12 *= u.kilojoules_per_mole
-        return a1, a2, NonbondedExceptionType(cs6, cs12, self.defaults.fudgeQQ)
+        eps, rmin, a1, a2 = eps_rmin_a1_a2_for_line(line)
+        return a1, a2, NonbondedExceptionType(rmin, eps, self.defaults.fudgeQQ)
 
     #===================================================
 
@@ -2053,3 +2069,23 @@ def _gettype(atom):
     if atom.atom_type not in (None, UnassignedAtomType):
         return atom.atom_type.bond_type
     return atom.type
+
+def eps_rmin_a1_a2_for_line(line):
+    words = line.split()
+    a1, a2 = words[:2]
+    func = int(words[2])
+    if func == 1:
+        C6_gromacs, C12_gromacs = (float(x) for x in words[3:5])
+        C6, C12 = C6_gromacs, C12_gromacs
+        if any(x == 0. for x in (C6, C12)):
+            eps, sig = (0., 0.)
+        else:
+            eps, sig = (C6**2 / (4 * C12), (C12 / C6)**(1./6))
+    elif func in (2, 3):
+        sig, eps = (float(x) for x in words[3:5])
+    else:
+        raise Exception('Invalid nonbond_params func={0}'.format(func))
+    sig *= 10 # Convert to Angstroms
+    eps *= u.kilojoule.conversion_factor_to(u.kilocalorie)
+    rmin = sig*2**(1./6)
+    return (eps, rmin, a1, a2)