Merge branch 'develop' into opt_topology_serialization

.github/actions/setup-deps/action.yaml

@@ -54,7 +54,7 @@ inputs:
   clustalw:
     default: 'clustalw=2.1'
   h5py:
-    default: 'h5py'
+    default: 'h5py>=2.10'
   joblib:
     default: 'joblib>=0.12'
   netcdf4:

.github/workflows/gh-ci.yaml

@@ -57,7 +57,7 @@ jobs:
               full-deps: false
               install_hole: false
               codecov: false
-              numpy: numpy=1.19.0
+              numpy: numpy=1.20.0
             - name: asv_check
               os: ubuntu-latest
               python-version: 3.8

.travis.yml

@@ -44,7 +44,7 @@ matrix:
       group: edge
       if: type = cron
       before_install:
-        - python -m pip install cython "numpy>=1.19.2" scipy
+        - python -m pip install cython "numpy>=1.20.0" scipy
         - python -m pip install --no-build-isolation hypothesis matplotlib packaging pytest pytest-cov pytest-xdist tqdm threadpoolctl fasteners
       install:
         - cd package

azure-pipelines.yml

@@ -40,25 +40,25 @@ jobs:
           PYTHON_VERSION: '3.9'
           PYTHON_ARCH: 'x64'
           BUILD_TYPE: 'wheel'
-          NUMPY_MIN: '1.19.3'
+          NUMPY_MIN: '1.20.0'
           imageName: 'windows-2019'
         Win-Python38-64bit-full-wheel:
           PYTHON_VERSION: '3.8'
           PYTHON_ARCH: 'x64'
           BUILD_TYPE: 'wheel'
-          NUMPY_MIN: '1.19.0'
+          NUMPY_MIN: '1.20.0'
           imageName: 'windows-2019'
         Linux-Python39-64bit-full-wheel:
           PYTHON_VERSION: '3.9'
           PYTHON_ARCH: 'x64'
           BUILD_TYPE: 'wheel'
-          NUMPY_MIN: '1.19.3'
+          NUMPY_MIN: '1.20.0'
           imageName: 'ubuntu-latest'
         Linux-Python38-64bit-full-wheel:
           PYTHON_VERSION: '3.8'
           PYTHON_ARCH: 'x64'
           BUILD_TYPE: 'wheel'
-          NUMPY_MIN: '1.19.3'
+          NUMPY_MIN: '1.20.0'
           imageName: 'ubuntu-latest'
   pool:
     vmImage: $(imageName)
@@ -92,7 +92,7 @@ jobs:
       pytest-xdist
       scikit-learn
       scipy
-      h5py
+      h5py>=2.10
       tqdm
       threadpoolctl
       fasteners

maintainer/conda/environment.yml

@@ -9,14 +9,14 @@ dependencies:
   - fasteners
   - griddataformats
   - gsd
-  - h5py
+  - h5py>=2.10
   - hypothesis
   - joblib>=0.12
   - matplotlib==3.2.2
   - mmtf-python
   - mock
   - networkx
-  - numpy>=1.19
+  - numpy>=1.20
   - pytest
   - python==3.8
   - scikit-learn

package/AUTHORS

@@ -196,6 +196,7 @@ Chronological list of authors
   - Ricky Sexton
   - Rafael R. Pappalardo
   - Tengyu Xie
+  - Raymond Zhao
 
 
 External code

package/CHANGELOG

@@ -13,15 +13,22 @@ The rules for this file:
   * release numbers follow "Semantic Versioning" http://semver.org
 
 ------------------------------------------------------------------------------
-??/??/?? IAlibay, PicoCentauri, orbeckst, hmacdope, yuxuanzhuang, rmeli, miss77jun
+??/??/?? IAlibay, PicoCentauri, orbeckst, hmacdope, rmeli, miss77jun, rzhao271,
+         yuxuanzhuang
 
  * 2.3.0
 
 Fixes
   * Fixes awk call in deploy.yaml tests for macos runners (Issue #3693)
   * add a 0.5 for correct midpoints in hole analysis (Issue #3715)
+  * Fix reading error when PDB CONECT records are corrupt. (Issue #988)
+  * Fix writing unusable PDB CONECT records with index>100000. (Issue #988)
 
 Enhancements
+  * Add a new `formalcharge` attribute for storing formal charges (PR #3755)
+  * Add the ability to read & write formal charges from PDB files (PR #3755)
+  * Change functions in `lib.distances` to accept AtomGroups as well as NumPy
+    arrays (CZI Performance track PR #3730)
   * Add `norm` parameter to InterRDF, InterRDF_s to normalize as rdf,
     number density or do not normalize at all. (Issue #3687)
   * Additional logger.info output when per-frame analysis starts (PR #3710)
@@ -31,18 +38,25 @@ Enhancements
   * Optimize topology serialization/construction by lazy-building _RA and _SR.
 
 Changes
+  * To add additional optional packages for different file formats when
+    using `pip` use the new "extra_formats" extras_require target as in
+    `pip install ./package[extra_formats]` (Issue #3701, PR #3711)
+  * Minimum supported NumPy version is now 1.20 (1.21 for macosx-arm64)
+    as per NEP29 (PR #3737)
   * Narrowed variable scope to reduce use of OpenMP `private` clause (PR #3706, PR
     #3728)
 
 Deprecations
+  * The extras_requires target "AMBER" for `pip install ./package[AMBER]`
+    will be removed in 2.4. Use "extra_formats". (Issue #3701, PR #3711)
   * Deprecate `density` parameter in favor of `norm` in InterRDF_s
     (Issue #3687)
 
 
 06/01/22 IAlibay, BFedder, inomag, Agorfa, aya9aladdin, shudipto-amin, cbouy,
          HenokB, umak1106, tamandeeps, Mrqeoqqt, megosato, AnirG, rishu235,
          mtiberti, manishsaini6421, Sukeerti1, robotjellyzone, markvrma, alescoulie,
-         mjtadema, PicoCentauri, Atharva7K, aditi2906, orbeckst,  yuxuanzhuang,
+         mjtadema, PicoCentauri, Atharva7K, aditi2906, orbeckst, yuxuanzhuang,
          rsexton2, rafaelpap, richardjgowers, orioncohen
 
  * 2.2.0

package/MDAnalysis/coordinates/PDB.py

@@ -540,24 +540,27 @@ class PDBWriter(base.WriterBase):
        An indexing issue meant it previously used the first charater (Issue #2224)
 
     .. versionchanged:: 2.0.0
-        Add the `redindex` argument. Setting this keyword to ``True``
-        (the default) preserves the behavior in earlier versions of MDAnalysis.
-        The PDB writer checks for a valid chainID entry instead of using the
-        last character of segid. Should a chainID not be present, or not
-        conform to the PDB standard, the default value of  'X' is used.
-
+       Add the `redindex` argument. Setting this keyword to ``True``
+       (the default) preserves the behavior in earlier versions of MDAnalysis.
+       The PDB writer checks for a valid chainID entry instead of using the
+       last character of segid. Should a chainID not be present, or not
+       conform to the PDB standard, the default value of  'X' is used.
+
+    .. versionchanged:: 2.3.0
+       Do not write unusable conect records when ag index
+       is larger than 100000.
     """
     fmt = {
         'ATOM': (
             "ATOM  {serial:5d} {name:<4s}{altLoc:<1s}{resName:<4s}"
             "{chainID:1s}{resSeq:4d}{iCode:1s}"
             "   {pos[0]:8.3f}{pos[1]:8.3f}{pos[2]:8.3f}{occupancy:6.2f}"
-            "{tempFactor:6.2f}      {segID:<4s}{element:>2s}\n"),
+            "{tempFactor:6.2f}      {segID:<4s}{element:>2s}{charge:2s}\n"),
         'HETATM': (
             "HETATM{serial:5d} {name:<4s}{altLoc:<1s}{resName:<4s}"
             "{chainID:1s}{resSeq:4d}{iCode:1s}"
             "   {pos[0]:8.3f}{pos[1]:8.3f}{pos[2]:8.3f}{occupancy:6.2f}"
-            "{tempFactor:6.2f}      {segID:<4s}{element:>2s}\n"),
+            "{tempFactor:6.2f}      {segID:<4s}{element:>2s}{charge:2s}\n"),
         'REMARK': "REMARK     {0}\n",
         'COMPND': "COMPND    {0}\n",
         'HEADER': "HEADER    {0}\n",
@@ -851,6 +854,10 @@ def _write_pdb_bonds(self):
         for a1, a2 in bonds:
             if not (a1 in mapping and a2 in mapping):
                 continue
+            if mapping[a1] >= 100000 or mapping[a2] >= 100000:
+                warnings.warn("Atom with index >=100000 cannot write "
+                              "bonds to PDB CONECT records.")
+                return
             con[a2].append(a1)
             con[a1].append(a2)
 
@@ -1023,6 +1030,47 @@ def _deduce_PDB_atom_name(self, atomname, resname):
             return '{:<4}'.format(atomname)
         return ' {:<3}'.format(atomname)
 
+    @staticmethod
+    def _format_PDB_charges(charges: np.ndarray) -> np.ndarray:
+        """Format formal charges to match PDB requirements.
+
+        Formal charge entry is set to empty if charge is 0, otherwise the
+        charge is set to a two character ```<charge value><charge sign>``
+        entry, e.g. ``1+`` or ``2-``.
+
+        This method also verifies that formal charges can adhere to the PDB
+        format (i.e. charge cannot be > 9 or < -9).
+
+        Parameters
+        ----------
+        charges: np.ndarray
+            NumPy array of integers representing the formal charges of
+            the atoms being written.
+
+        Returns
+        -------
+        np.ndarray
+            NumPy array of dtype object with strings representing the
+            formal charges of the atoms being written.
+        """
+        if not np.issubdtype(charges.dtype, np.integer):
+            raise ValueError("formal charges array should be of `int` type")
+
+        outcharges = charges.astype(object)
+        outcharges[outcharges == 0] = ''  # empty strings for no charge case
+        # using np.where is more efficient than looping in sparse cases
+        for i in np.where(charges < 0)[0]:
+            if charges[i] < -9:
+                errmsg = "formal charge < -9 is not supported by PDB standard"
+                raise ValueError(errmsg)
+            outcharges[i] = f"{abs(charges[i])}-"
+        for i in np.where(charges > 0)[0]:
+            if charges[i] > 9:
+                errmsg = "formal charge > 9 is not supported by PDB standard"
+                raise ValueError(errmsg)
+            outcharges[i] = f"{charges[i]}+"
+        return outcharges
+
     def _write_timestep(self, ts, multiframe=False):
         """Write a new timestep *ts* to file
 
@@ -1093,6 +1141,7 @@ def get_attr(attrname, default):
         atomnames = get_attr('names', 'X')
         elements = get_attr('elements', ' ')
         record_types = get_attr('record_types', 'ATOM')
+        formal_charges = self._format_PDB_charges(get_attr('formalcharges', 0))
 
         def validate_chainids(chainids, default):
             """Validate each atom's chainID
@@ -1158,6 +1207,7 @@ def validate_chainids(chainids, default):
             vals['segID'] = segids[i][:4]
             vals['chainID'] = chainids[i]
             vals['element'] = elements[i][:2].upper()
+            vals['charge'] = formal_charges[i]
 
             # record_type attribute, if exists, can be ATOM or HETATM
             try:

package/MDAnalysis/core/groups.py

@@ -3041,6 +3041,11 @@ def select_atoms(self, sel, *othersel, periodic=True, rtol=1e-05,
                 ``S`` will be possible options but other values will not raise
                 an error.
 
+            formalcharge *formal-charge*
+                select atoms based on their formal charge, e.g.
+                ``name O and formalcharge -1`` to select all oxygens with a
+                negative 1 formal charge.
+
         **Boolean**
 
             not

package/MDAnalysis/core/topologyattrs.py

@@ -2033,6 +2033,18 @@ def total_charge(group, compound='group'):
         ('total_charge', total_charge))
 
 
+class FormalCharges(AtomAttr):
+    """Formal charge on each atom"""
+    attrname = 'formalcharges'
+    singular = 'formalcharge'
+    per_object = 'atom'
+    dtype = int
+
+    @staticmethod
+    def _gen_initial_values(na, nr, ns):
+        return np.zeros(na)
+
+
 # TODO: update docs to property doc
 class Occupancies(AtomAttr):
     attrname = 'occupancies'