Add support for more CHARMM force fields (#1157)

* Add features to CHARMM files

* Adds support for LP (lone pairs) in CHARMM PSF files
* Adds a LocalCoordinatesFrame to support the CHARMM LP types
* Starts to add Drude parsing to CHARMM psf parser

* Avoid possibly duplicating constraints

This is especially possible in water where a H-H bond may be added for
constraints, but that may exist alongside a H-O-H angle. If both exist,
the current code will add the same constraint twice which could cause
problems. Now the angle constraint only adds a distance constraint if
those two atoms are not *already* constrained.

This fix is inspired from an OpenMM commit for the CHARMM files.

* Add test files for drude CHARMM systems

* Add a test case for parsing a DRUDE PSF and fix some issues with it

* Add some protection for changeradii against index errors

* Checkpoint commit

* More progress on Drude support

* More progress towards Drude force fields in CHARMM

Properly parses the drude anisotropies with relevant tests.

* Clean up a bunch of stuff and implement CHARMM PSF writing with drude and lone pairs

* Fix bug

.github/workflows/Test.yaml

@@ -28,6 +28,8 @@ jobs:
         python-version: ${{ matrix.python-version }}
         environment-file: devtools/environment-dev.yaml
         activate-environment: parmed-dev
+        channels: conda-forge,bioconda
+        mamba-version: "*"
 
     - name: Environment Information
       shell: bash -l {0}

devtools/ci/install.sh

@@ -11,5 +11,15 @@ python -c "import parmed; print(parmed.__version__)"
 echo "Using ParmEd version `parmed --version`"
 cd test
 ./run_scripts.sh
-py.test -n 4 --cov-branch --cov=parmed --durations-min=5 --disable-warnings --cov-append --cov-report=xml .
+pytest \
+    -n 4 \
+    --cov-branch \
+    --cov=parmed \
+    --cov-append \
+    --cov-report=xml \
+    --durations-min=5 \
+    -W ignore::UserWarning \
+    -W ignore::parmed.exceptions.ParmedWarning \
+    .
+
 echo "Done!"

devtools/environment-dev.yaml

@@ -1,7 +1,7 @@
 name: parmed-dev
 channels:
-  - conda-forge
-  - bioconda
+  - https://conda.anaconda.org/conda-forge
+  - https://conda.anaconda.org/bioconda
 dependencies:
   - pandas
   - nose

parmed/charmm/parameters.py

@@ -415,7 +415,7 @@ def read_parameter_file(self, pfile, comments=None):
                 section = None
                 continue
             if line.upper().startswith('THOLE'):
-                section = None
+                section = 'NBTHOLE'
                 continue
             # It seems like files? sections? can be terminated with 'END'
             if line[:3].upper() == 'END':
@@ -726,8 +726,22 @@ def read_parameter_file(self, pfile, comments=None):
                 except IndexError:
                     raise CharmmError('Could not parse NBFIX terms.')
                 self.nbfix_types[(min(at1, at2), max(at1, at2))] = (emin, rmin)
-        # If we had any CMAP terms, then the last one will not have been added
-        # yet. Add it here
+                continue
+            if section.upper() == 'NBTHOLE':
+                words = line.split()
+                try:
+                    at1 = words[0]
+                    at2 = words[1]
+                    nbt = abs(conv(words[2], float, 'NBTHOLE a'))
+                    try:
+                        self.atom_types_str[at1].add_nbthole(at2, nbt)
+                        self.atom_types_str[at2].add_nbthole(at1, nbt)
+                    except KeyError:
+                        pass
+                except IndexError as err:
+                    raise CharmmError("Could not parse NBTHOLE terms.") from err
+                self.nbthole_types[(at1, at2)] = self.nbthole_types[(at2, at1)] = nbt
+        # If we had any CMAP terms, then the last one will not have been added yet. Add it here
         if current_cmap is not None:
             typ = CmapType(current_cmap_res, current_cmap_data)
             self.cmap_types[current_cmap] = typ
@@ -800,7 +814,7 @@ def read_topology_file(self, tfile):
                 elif line[:4].upper() == 'DEFA':
                     words = line.split()
                     if len(words) < 5:
-                        warnings.warn('DEFA line has %d tokens; expected 5' % len(words))
+                        warnings.warn('DEFA line has %d tokens; expected 5' % len(words), ParameterWarning)
                     else:
                         it = iter(words[1:5])
                         for tok, val in zip(it, it):
@@ -811,7 +825,7 @@ def read_topology_file(self, tfile):
                             elif tok.upper() == 'LAST':
                                 tpatch = val
                             else:
-                                warnings.warn(f'DEFA patch {val} unknown')
+                                warnings.warn(f'DEFA patch {val} unknown', ParameterWarning)
                 elif line[:4].upper() in ('RESI', 'PRES'):
                     restype = line[:4].upper()
                     # Get the residue definition
@@ -825,7 +839,7 @@ def read_topology_file(self, tfile):
                     try:
                         charge = float(words[2])
                     except (IndexError, ValueError):
-                        warnings.warn(f'No charge for {resname}')
+                        warnings.warn(f'No charge for {resname}', ParameterWarning)
                     if restype == 'RESI':
                         res = ResidueTemplate(resname)
                     elif restype == 'PRES':
@@ -879,7 +893,7 @@ def read_topology_file(self, tfile):
                             else:
                                 warnings.warn(
                                     f'WARNING: Ignoring "{line.strip()}" because entity type '
-                                    f'{entity_type} not used.'
+                                    f'{entity_type} not used.', ParameterWarning
                                 )
                         elif line.strip().upper() and line.split()[0].upper() in ('BOND', 'DOUBLE'):
                             it = iter([w.upper() for w in line.split()[1:]])
@@ -919,7 +933,7 @@ def read_topology_file(self, tfile):
                             if not skip_adding_residue and lptype_keyword not in ['BISE', 'RELA']:
                                 warnings.warn(
                                     f'LONEPAIR type {words[1]} not supported; only BISEctor and '
-                                     'RELAtive supported'
+                                     'RELAtive supported', ParameterWarning
                                 )
                                 skip_adding_residue = True
                                 break
@@ -993,14 +1007,14 @@ def read_topology_file(self, tfile):
                 try:
                     res.first_patch = self.patches[patch_name]
                 except KeyError:
-                    warnings.warn(f'Patch {patch_name} not found')
+                    warnings.warn(f'Patch {patch_name} not found', ParameterWarning)
 
             patch_name = tpatches[resname]
             if patch_name is not None:
                 try:
                     res.last_patch = self.patches[patch_name]
                 except KeyError:
-                    warnings.warn(f'Patch {patch_name} not found')
+                    warnings.warn(f'Patch {patch_name} not found', ParameterWarning)
         # Now update the residues and patches with the ones we parsed here
         self.residues.update(residues)
 
@@ -1032,7 +1046,7 @@ def read_stream_file(self, sfile):
                 self.read_parameter_file(section, comments)
             title, section, comments = f.next_section()
 
-    def write(self, top=None, par=None, stream=None, **kwargs):
+    def write(self, top=None, par=None, stream=None):
         """ Write a CHARMM parameter set to a file
 
         Parameters
@@ -1043,7 +1057,7 @@ def write(self, top=None, par=None, stream=None, **kwargs):
         par : str or file-like object, optional
             If provided, the parameters will be written to this file in PAR
             format. Either this or the ``str`` argument *must* be provided
-        str : str or file-like object, optional
+        stream : str or file-like object, optional
             If provided, the atom types and parameters will be written to this
             file as separate RTF and PAR cards that can be read as a CHARMM
             stream file. Either this or the ``par`` argument *must* be provided
@@ -1052,17 +1066,6 @@ def write(self, top=None, par=None, stream=None, **kwargs):
         ------
         ValueError if both par and str are None
         """
-        if stream is None and kwargs.get("str") is not None:
-            warnings.warn(
-                "Use 'stream' instead of 'str' to write CHARMM parameters", DeprecationWarning
-            )
-            stream = kwargs["str"]
-        if stream is not None and kwargs.get("str") != stream:
-            warnings.warn(
-                f"Ignoring 'str' ({kwargs.get('str')}), using 'stream' ({stream}) instead",
-                DeprecationWarning,
-            )
-
         if par is None and stream is None:
             raise ValueError('Must specify either par *or* str')