Lammps charmm dihedrals (#787)

* Adjust molecule handling to be consistently set as 0 index if a default value is provided

* Add kelvin units and charmm dihedrals to lammps

* Fix setting scaling factor with latest numpy 1.25

* Add workaround for reordering multiple charmm params

* fix charge being grabbed from site.atomtype

* [pre-commit.ci] auto fixes from pre-commit.com hooks
for more information, see https://pre-commit.ci

* Add parsing for cvff impropers from PeriodicTorsion

* [pre-commit.ci] auto fixes from pre-commit.com hooks
for more information, see https://pre-commit.ci

* Add back in charmm harmonic improper support

* fix failing tests in handling charges from atom_type and not site

* fix bug in shifting units but not normalizing by a unit array

* Updates to parmed handling

* Revert "fix bug in shifting units but not normalizing by a unit array"

This reverts commit 33a3957.

* remove unwanted imports

* Minor changes to decrease validation time in top compatibility during checking

* Change lammps formatting for handling layered dihedrals so copies of each line are made with multiple dihedral types

* Reorder lammps layered dihedrals, and accurately count duplicated dihedrals in the header n_dihedrals for dihedrals that are counted twice due to multiple layers of the dihedral

* add optional arguments to check_compatibility to provide quicker way to iterate through topology attributes with different potential filters

* Modify readability of some lammps writer functions

* Move lammpswriter functions outside of main writer, and move molecule reindex to sorting.py

* Update gmso/formats/lammpsdata.py

* fix tests with deprecated function

---------

Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Co-authored-by: Co Quach <[email protected]>
Co-authored-by: Co Quach <[email protected]>

gmso/core/topology.py

@@ -815,7 +815,7 @@ def _set_scaling_factor(self, value, molecule_id, interaction, name, index):
         else:
             if interaction not in scaling_interaction_idxes:
                 raise GMSOError(f"Unknown `{name}` interaction `{interaction}`")
-            all_scales[index][scaling_interaction_idxes[interaction]] = value
+            all_scales[index][scaling_interaction_idxes[interaction]] = value[0]
 
     def add_site(self, site, update_types=False):
         """Add a site to the topology.

gmso/external/convert_mbuild.py

@@ -315,7 +315,7 @@ def _parse_molecule_residue(site_map, compound):
         if molecule_tag.name in molecule_tracker:
             molecule_tracker[molecule_tag.name] += 1
         else:
-            molecule_tracker[molecule_tag.name] = 1
+            molecule_tracker[molecule_tag.name] = 0
         molecule_number = molecule_tracker[molecule_tag.name]
         """End of molecule parsing"""
 
@@ -330,7 +330,7 @@ def _parse_molecule_residue(site_map, compound):
                         residue_tracker[residue_tag.name]
                     )
             else:
-                residue_tracker[residue_tag.name] = {residue_tag: 1}
+                residue_tracker[residue_tag.name] = {residue_tag: 0}
 
             residue_number = residue_tracker[residue_tag.name][residue_tag]
             site_map[particle]["residue"] = (residue_tag.name, residue_number)

gmso/external/convert_parmed.py

@@ -67,20 +67,26 @@ def from_parmed(structure, refer_type=True):
             element = (
                 element_by_atomic_number(atom.element) if atom.element else None
             )
+            if residue.number == -1:  # use default value of 0 in GMSO
+                residue_number = 0
+            else:
+                residue_number = residue.number - 1
             site = gmso.Atom(
                 name=atom.name,
                 charge=atom.charge * u.elementary_charge,
                 position=[atom.xx, atom.xy, atom.xz] * u.angstrom,
                 atom_type=None,
-                residue=(residue.name, residue.number),
+                residue=(residue.name, residue_number),
                 element=element,
             )
-            site.molecule = (residue.name, residue.number) if ind_res else None
+            site.molecule = (residue.name, residue_number) if ind_res else None
             site.atom_type = (
                 copy.deepcopy(pmd_top_atomtypes[atom.atom_type])
                 if refer_type and isinstance(atom.atom_type, pmd.AtomType)
                 else None
             )
+            if site.atom_type:
+                site.atom_type.charge = atom.charge * u.elementary_charge
             site_map[atom] = site
             top.add_site(site)
 
@@ -480,7 +486,7 @@ def to_parmed(top, refer_type=True):
             structure.add_atom(
                 pmd_atom,
                 resname=site.residue.name,
-                resnum=site.residue.number,
+                resnum=site.residue.number + 1,
             )
         else:
             structure.add_atom(pmd_atom, resname="RES", resnum=-1)

gmso/formats/gro.py

@@ -64,7 +64,9 @@ def read_gro(filename):
                     "atoms were expected, but at least one fewer was found."
                 )
                 raise ValueError(msg.format(n_atoms))
-            res_id = int(line[:5].strip())
+            res_id = (
+                int(line[:5].strip()) - 1
+            )  # reformat from 1 to 0 index in gmso
             res_name = line[5:10].strip()
             atom_name = line[10:15].strip()
             atom_id = line[15:20].strip()