Single molecule force fields

[rsdefever]

Our group has many custom force fields for individual molecules that have been developed over the years. We would like to use gmso to store and use these force fields in our simulations. This results in a slightly different workflow than the typical workflow, so this thread will be a reference for the features we need in gmso to accomplish our goals. Some of these things may already be possible in gmso, but I am including everything for the sake of brainstorming and discussion.

Workflow:

1. Load one or more force field(s) for a fully typed single molecule(s) from disk. This should include all atoms, atomtypes, partial charges, atomic positions, bonds, angles, dihedrals, impropers, and all the associated parameters and functional forms, for each molecule. I imagine this will be a serialized (on disk) representation of a gmso.Topology for a single molecule.
2. Create a system with some number of each molecule (e.g., 300 BMIM and 300 PF6). This step would most likely be done with the fill_box functionality of mbuild.
3. Copy or replicate the single molecule force field files so that we have a complete description of our system as a gmso.Topology.
4. Write force field input files (LAMMPS, GROMACS, Cassandra) to disk.

Needs:

1. Serialized representation of gmso.Topology.
2. Ability to replicate a gmso.Topology many times.
3. Ability to ensure there are no conflicts. E.g., if molecule_1 has an atom type ca and molecule_2 also has a ca atom type, we need to either:
   • Ensure they have the same parameters, or
   • Ensure that they will be written to disk as different atom types with different parameters if the engine requires it
4. Check that the specification for e.g., mixing rules, 1-4 interactions, etc, are consistent

That is what I have for now. If anyone has additional thoughts on the needs/workflow/etc, please pitch in below.

[umesh-timalsina]

@rsdefever Thanks for the issue. I think this workflow in on itself would be a huge milestone for gmso to meet.

Workflow:

1. Load one or more force field(s) for a single molecule(s) from disk. This should include all atoms, atomtypes, partial charges, atomic positions, bonds, angles, dihedrals, impropers, and all the associated parameters and functional forms, for each molecule. I imagine this will be a serialized (on disk) representation of a gmso.Topology for a single molecule.

What formats are the molecules represented currently? The serialized on disk representation of a Topology depends on #494, which I assume is almost resolved. Also, separating potentials from structures in a Topology can be done which of which a rough outline has been discussed here. Also, if we don't follow the serialization route an alternative would be to write a reader for the Molecules (from disk) to create a typed Topology. Then, we can extend the ForceField class to implement

@classmethod
ForceField.from_topology(top: Topology, *args, **kwargs) -> 'ForceField':

2. Create a system with some number of each molecule (e.g., 300 BMIM and 300 PF6). This step would most likely be done with the fill_box functionality of mbuild.

This would be great, but I am assuming would require an extra step for conversion to mbuild <-> gmso. I am not entirely sure if this is possible with some refactoring of this function to incorporate this function to take in Topology.

3. Copy or replicate the single molecule force field files so that we have a complete description of our system as a gmso.Topology.

This part is not clear to me. Can we discuss this further?

4. Write force field input files (LAMMPS, GROMACS, Cassandra) to disk.

We have some support for these writers pinned in #6

Needs:

1. Serialized representation of gmso.Topology.

Overall, could be achieved and in development.

2. Ability to replicate a gmso.Topology many times.

See, #32 and #463

3. Ability to ensure there are no conflicts. E.g., if molecule_1 has an atom type ca and molecule_2 also has a ca atom type, we need to either:

   • Ensure they have the same parameters, or
   • Ensure that they will be written to disk as different atom types with different parameters if the engine requires it

This is already implemented and we can support it with no/minor customization

4. Check that the specification for e.g., mixing rules, 1-4 interactions, etc, are consistent

Same as 3.

[mattwthompson]

You may wish to consider residue templates instead of force fields for cases like this. For small, custom force fields that are only valid for one molecule - and likely to only be used with one molecule - you can bypass the typing step and most of the complexity associated with force fields. In GROMACS world this is an ITP file and I think there is a similar equivalent in AMBER world for amino acids. What you describe in your first step is similar (but not identical) to what people do for proteins.

[rsdefever]

You may wish to consider residue templates instead of force fields for cases like this. For small, custom force fields that are only valid for one molecule - and likely to only be used with one molecule - you can bypass the typing step and most of the complexity associated with force fields. In GROMACS world this is an ITP file and I think there is a similar equivalent in AMBER world for amino acids. What you describe in your first step is similar (but not identical) to what people do for proteins.

Yes, and this is what we have in our group now, in a custom format. I have updated my original comment to not use the term "force field", because what I was envisioning was more like a residue template.

What formats are the molecules represented currently?

Currently they are in a custom Maginn force field archive (MaFFA) format.

@umesh-timalsina what is the status of #462? It looks like there is not yet a method to serialize an entire Topology?

This would be great, but I am assuming would require an extra step for conversion to mbuild <-> gmso.

I don't mind the conversion. I don't think it would be too difficult to allow that function to take a gmso.Topology, if we wanted.

This part is not clear to me. Can we discuss this further?

What I mean is that after (1) and (2) we would have fully parameterized gmso.Topology objects for a single molecule, and then a unparameterized gmso.Topology object with N molecules (from the mbuild.fill_box). We need some way to apply, as @mattwthompson calls it, the residue template to the N molecules in the system so that we have a fully parameterized gmso.Topology object for all N molecules.