Fixes to general sorting algorithms for hoomd and lammps #762

CalCraven · 2023-09-11T12:19:28Z

HOOMD functionality right now leads to bugs in writing out the forcefields.

There are two functions that need to be in sync
gmso/external/convert_hoomd.py to_hoomd_snapshot
gmso/external/convert_hoomd.py to_hoomd_forcefield

These are not compatible, as discussed in #758.

Additionally, the sorting algorithms in gmso.utils.sorting.py have some bugs with dihedrals, which need to be flipped if the two central atoms are equivalent, but the first and fourth atoms should be in reversed order.

This PR will restructure some of the sorting algorithms and migrate some of the ones that were used in formats/lammpsdata.py from gmso/core/views.py to gmso/utils/sorting.py

TODO List

codecov · 2023-09-11T12:43:56Z

Codecov Report

Patch coverage: 77.86% and project coverage change: -0.18% ⚠️

Comparison is base (24a6101) 91.89% compared to head (06e9ddd) 91.71%.

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #762      +/-   ##
==========================================
- Coverage   91.89%   91.71%   -0.18%     
==========================================
  Files          66       66              
  Lines        6733     6780      +47     
==========================================
+ Hits         6187     6218      +31     
- Misses        546      562      +16

Files Changed	Coverage Δ
gmso/external/convert_hoomd.py	`91.44% <76.36%> (-0.79%)`	⬇️
gmso/utils/sorting.py	`80.89% <77.77%> (-5.38%)`	⬇️
gmso/core/views.py	`88.46% <100.00%> (-0.79%)`	⬇️
gmso/formats/lammpsdata.py	`95.89% <100.00%> (+<0.01%)`	⬆️

... and 1 file with indirect coverage changes

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

…types

gmso/external/convert_hoomd.py

+            [site.atom_type.atomclass for site in dihedral.connection_members]
+        )
+        unique_dihedrals[unique_members] = dihedral
+    unique_dtypes = [


gmso/external/convert_hoomd.py

        # TODO: The range of ks is mismatched (GMSO go from k0 to k5)
        # May need to do a check that k0 == k5 == 0 or raise a warning
-        container.params["-".join(dtype.member_types)] = {
+        member_classes = sort_by_classes(dtype)


daico007

Left some comments about the code. One high-level comment I have (probably need a broader discussion later too), so it seems right now, we are going with atom class for everything (storing bond and bond params in the snapshot and forcefield), which would work for most of our current use case, i.e., when connection types are defined between atom classes. This, however, may throw an issue if there is a connection that is defined between atom types. The latter case is rarer but, IMO, it might happen and should be supported. In other word, I think we should go with atom types for everything instead of atomclass. This will make thing a bit bulkier, but it will allow us to cover those edge case.

daico007 · 2023-09-18T15:34:51Z

gmso/utils/sorting.py

-    else:
-        raise TypeError("Provided connection_type not supported.")
-
-
 def sort_connection_members(connection, sort_by="name"):
    """Sort connection_members of connection."""


I think we need to add a blurp here about what options sort_by can take

daico007 · 2023-09-18T15:44:09Z

gmso/utils/sorting.py

+    elif isinstance(potential, ImproperType):
+        return (
+            potential.member_classes[0],
+            *potential.member_classes[1:],


Suggested change

*potential.member_classes[1:],

*[sorted(potential.member_classes[1:]),

daico007 · 2023-09-18T15:44:37Z

gmso/utils/sorting.py

+    elif isinstance(potential, ImproperType):
+        return (
+            potential.member_types[0],
+            *potential.member_types[1:],


Suggested change

*potential.member_types[1:],

*sorted(potential.member_types[1:]),

daico007 · 2023-09-18T15:47:13Z

gmso/utils/sorting.py

+    elif len(namesList) == 4 and improperBool:
+        return tuple(
+            namesList[0],
+            sorted(*namesList[1:]),


Suggested change

sorted(*namesList[1:]),

*sorted(namesList[1:]),

daico007 · 2023-09-19T12:17:27Z

gmso/utils/sorting.py

    elif isinstance(connection, gmso.Improper):
        site1, site2, site3, site4 = connection.connection_members
        site2, site3, site4 = sorted([site2, site3, site4], key=sorting_key)
        return [site1, site2, site3, site4]
    else:
        raise TypeError("Provided connection not supported.")
+
+
+def sort_by_classes(potential):


Might consider using a more descriptive name for this method, since sort_by_classes is kinda vague (regarding the desired input). Maybe something like sorted_potential_member_classes (sounds a bit clumsy but I will keep brainstorming some other option).

I guess this method is similar for sort_by_types, so probably we need something that would distinct these two (sort_by_types and sort_by_classes) with the sort_connection_strings.

other names I can think sort_potential_member_types and sort_potential_member_classes

Fixed by adding more detailed doc strings. I don't think every function needs to be that detailed. potential is the only argument, so we don't need it to also be part of the function name. I added details about what potential is and how it's sorted.

daico007 · 2023-09-19T12:36:07Z

gmso/tests/test_hoomd.py

+            if (  # TODO: why are these skipped?
+                isinstance(mb_force, hoomd.md.long_range.pppm.Coulomb)
+                or isinstance(mb_force, hoomd.md.pair.pair.LJ)
+                or isinstance(mb_force, hoomd.md.special_pair.LJ)
+                or isinstance(mb_force, hoomd.md.pair.pair.Ewald)
+                or isinstance(mb_force, hoomd.md.special_pair.Coulomb)
+            ):


Suggested change

if ( # TODO: why are these skipped?

isinstance(mb_force, hoomd.md.long_range.pppm.Coulomb)

or isinstance(mb_force, hoomd.md.pair.pair.LJ)

or isinstance(mb_force, hoomd.md.special_pair.LJ)

or isinstance(mb_force, hoomd.md.pair.pair.Ewald)

or isinstance(mb_force, hoomd.md.special_pair.Coulomb)

):

if ( # TODO: why are these skipped?

isinstance(mb_force, [

hoomd.md.long_range.pppm.Coulomb,

hoomd.md.pair.pair.LJ,

hoomd.md.special_pair.LJ,

hoomd.md.pair.pair.Ewald,

hoomd.md.special_pair.Coulomb,

]

):

daico007 · 2023-09-19T12:37:41Z

gmso/tests/test_hoomd.py

+                continue
+            keys = mb_force.params.param_dict.keys()
+            gmso_keys = gmso_force.params.param_dict.keys()
+            print("\n\n", keys, gmso_keys, gmso_force, "\n\n")


Suggested change

print("\n\n", keys, gmso_keys, gmso_force, "\n\n")

daico007 · 2023-09-19T12:59:35Z

gmso/external/convert_hoomd.py

@@ -1413,3 +1430,26 @@ def _convert_params_units(
        potential.parameters = converted_params
        converted_potentials.append(potential)
    return converted_potentials
+
+
+def _convert_single_param_units(


should this belong to utils/units.py?

I agree. It should get moved into using the new UnitsClasses, but maybe we can do that in a later PR.

CalCraven · 2023-09-20T10:33:21Z

One high-level comment I have (probably need a broader discussion later too), so it seems right now, we are going with atom class for everything (storing bond and bond params in the snapshot and forcefield), which would work for most of our current use case, i.e., when connection types are defined between atom classes. This, however, may throw an issue if there is a connection that is defined between atom types. The latter case is rarer but, IMO, it might happen and should be supported. In other word, I think we should go with atom types for everything instead of atomclass.

So I think it's a necessity to use atomclass instead of atomtype. This is related to #765, which is to try to get hoomd to support GAFF. In fact, that PR is updated with the changes to this PR, because both are necessary for GAFF, so in the end I'll probably address these comments/fixes directly in there and close this PR. But I'm getting away from the point. Here's the situation that arises.

A topology has 5 atomtypes, each with unique names, bonded in a line to form 2 unique dihedrals. When writing out the snapshot and the forcefield, the important attribute is the string that is used as the key for these lookup tables. If we use atomtypes, then it'll look like:

Molecule representation -> atom1 - atom2 - atom3 - atom4 - atom5
In snapshots
Dihedral1 = "atype1-atype2-atype3-atype4"
Dihedral2 = "atype2-atype3-atype4-atype5"

in forces
Dihedral_type1 = "atype1-atype2-atype3-atype4"
Dihedral_type2 = "atype2-atype3-atype4-atype5"

But if the classes of dihedrals are all the same, then we're duplicating dihedrals.
Dihedral_type1 = Dihedral_type2. If we used classes, we would reduce the number of dihedrals written. And 90% or more of our forcefields are defined by atomclasses, and with big forcefields we're at risk of writing huge numbers of dihedrals out with all of the duplication.

Now the situation you're describing, where the dihedrals are defined by their atomtypes, not atomclasses. In this scenario, the atomtypes = atomclasses.
atomtype1 = "OPLS_135" atomclass1 = "OPLS_135". In that case, we should still be okay to use the classes as they are equivalent in using their strings.

The strange situation would be:
atomtype1 = "OPLS_135" atomclass1 = "CT"
atomtype2="OPLS_136 atomclass1="CT"
but two different dihedrals are defined as
"OPLS_135-OPLS_136-X-X"
"CT-CT-X-X"
This would write these two dihedrals the same way. But I don't thing writing the atomtypes would be a useful solution here as well. Because then what are the atomtypes for the "CT-CT-X-X" dihedral? In theory it could be any combination of atomtypes that correlate to that class. So I'm not sure we could handle that situation either way, unless we know the original formulation of the dihedral.

for dihedral in top.dihedrals:
        unique_members = tuple(
            [site.atom_type.atomclass for site in dihedral.connection_members]
        )
        unique_dihedrals[unique_members] = dihedral

daico007 · 2023-09-20T14:43:22Z

I am not too worry about the duplication issue (since that will be more or less constrained by how many unique molecules in a systems and how many unique dihedrals in those molecules). There will be duplications, but I don't think that will be to that much of a degree. Consider that against the trade-off of supporting the 10% non-common case, that might be worth it.

What I mean about dihedral define between atomtypes, I mean it's look like:
type1="aa" type2="bb" type3="cc" type4="dd". Making atomclass to take the value of atomtype is just a work-around.

This topic has been brought up a while ago, but in the strange case that you brought up, the dihedral define by type will take precedent to that defined by atomclass.

CalCraven · 2023-09-20T15:54:19Z

I think the duplication issue will get out of hand. For instance, I have simulations with up to 20 different atomtypes in them. But maybe only 10 or so atomclasses. I have the feeling that will result in a large number of dihedrals between the two cases. There's a reason something large like the OPLS forcefield defines the dihedrals by classes, not by types. I will put together an example of this and we can count the number of dihedrals, but run time will be affected by the number of dihedrals in hoomd no doubt about it. I've never actually seen a forcefield with hybrid type and class definitions. It seems a bit silly: if you have a unique site that needs it's own unique dihedral, then give it an atomclass that's the same as its atomtype. We certainly do not have any tested examples of such a forcefield if it is something we would want to support anyways.

CalCraven · 2023-09-20T15:55:08Z

I also am totally spitballing with that 10% uncommon case. I've never actually seen that done. So who knows exactly what that number is.

CalCraven · 2023-09-20T15:56:18Z

if the dihedral is type1="aa" and is given no class1, then it's class should also be "aa".

daico007 · 2023-09-20T16:11:50Z

It seems a bit silly: if you have a unique site that needs it's own unique dihedral, then give it an atomclass that's the same as its atomtype.

I am thinking more about the case where the bond/angle is defined on atomclass and the dihedral on atomtype (for some reason, like you try to develop a forcefield I try to have it very specific to the atom type level). In that case, the atom type still need its original atomclass info.

CalCraven · 2023-09-20T16:12:26Z

Having just checked this, it resulted in 23 unique dihedrals by types and 14 by classes.

CalCraven · 2023-09-20T16:13:57Z

test_dihedrals.txt

CalCraven · 2023-09-20T16:22:57Z

I am thinking more about the case where the bond/angle is defined on atomclass and the dihedral on atomtype (for some reason, like you try to develop a forcefield I try to have it very specific to the atom type level). In that case, the atom type still need its original atomclass info.

I can see this. We need tests for it though, which we don't currently have. I think we could have a check at the beginning of the writer that checks to see if all the classes capture all of the unique dihedral expressions. If not, then we populate the snapshot and forces by types, but that shouldn't be the default case.

I think this PR should get merged tomorrow for CECAM, and I don't have the time to add that feature in. So potentially we can just raise this as a bug in the issues, and address it when we have the time/ a reasonable test case.

daico007 · 2023-09-21T19:18:11Z

included in #765

Fixes to general sorting algorithms for hoomd and lammps

37fdf57

Generate dihedral_types in to_hoomd_forces by dihedral, not dihedral_…

adebaaa

…types

github-advanced-security bot found potential problems Sep 11, 2023

View reviewed changes

Merge branch 'main' into connection-members-classes-only

6fa1518

daico007 reviewed Sep 19, 2023

View reviewed changes

Merge branch 'main' into connection-members-classes-only

06e9ddd

daico007 closed this Sep 21, 2023

daico007 mentioned this pull request Sep 21, 2023

Sort indices in Frame.bond.group #769

Merged

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Fixes to general sorting algorithms for hoomd and lammps #762

Fixes to general sorting algorithms for hoomd and lammps #762

CalCraven commented Sep 11, 2023 •

edited

Loading

codecov bot commented Sep 11, 2023 •

edited

Loading

daico007 left a comment

daico007 Sep 18, 2023

CalCraven Sep 20, 2023

daico007 Sep 18, 2023

daico007 Sep 18, 2023

daico007 Sep 18, 2023

daico007 Sep 19, 2023

daico007 Sep 19, 2023

CalCraven Sep 20, 2023

daico007 Sep 19, 2023

daico007 Sep 19, 2023

daico007 Sep 19, 2023

CalCraven Sep 20, 2023

CalCraven commented Sep 20, 2023

daico007 commented Sep 20, 2023

CalCraven commented Sep 20, 2023

CalCraven commented Sep 20, 2023

CalCraven commented Sep 20, 2023

daico007 commented Sep 20, 2023

CalCraven commented Sep 20, 2023

CalCraven commented Sep 20, 2023

CalCraven commented Sep 20, 2023

daico007 commented Sep 21, 2023

	*potential.member_classes[1:],
	*[sorted(potential.member_classes[1:]),

	*potential.member_types[1:],
	*sorted(potential.member_types[1:]),

Fixes to general sorting algorithms for hoomd and lammps #762

Fixes to general sorting algorithms for hoomd and lammps #762

Conversation

CalCraven commented Sep 11, 2023 • edited Loading

codecov bot commented Sep 11, 2023 • edited Loading

Codecov Report

daico007 left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

CalCraven commented Sep 20, 2023

daico007 commented Sep 20, 2023

CalCraven commented Sep 20, 2023

CalCraven commented Sep 20, 2023

CalCraven commented Sep 20, 2023

daico007 commented Sep 20, 2023

CalCraven commented Sep 20, 2023

CalCraven commented Sep 20, 2023

CalCraven commented Sep 20, 2023

daico007 commented Sep 21, 2023

CalCraven commented Sep 11, 2023 •

edited

Loading

codecov bot commented Sep 11, 2023 •

edited

Loading