michellab / BioSimSpace

Code and resources for the EPSRC BioSimSpace project.
https://biosimspace.org
GNU General Public License v3.0
77 stars 19 forks source link

Check atom names in BSS generated pert files #24

Closed jmichel80 closed 5 years ago

jmichel80 commented 6 years ago

SOMD may require that each atom in a pert file has a unique Sire AtomName. As BSS does not rename input atoms, it may be necessary to generate unique atom names when writing prm7 and pert files for a SOMD free energy calculation. Thus the SOMD prm7 may also require SOMD specific atom names, together with heuristics for settings the atomic mass of perturbed elements.

lohedges commented 6 years ago

This is what will be done by BioSimSpace when writing the pert file and AMBER prm/rst files as part of the SOMD free energy protocol. I plan on making a copy of the user input molecule, then applying my pert file writer (which has to rename the atoms in order to make them unique). You can look at the _toPertFile method in BioSimSpace/python/_SireWrappers/_molecue.py to see how the renaming works.

I'll leave this issue as a reminder and close when the pert file / free energy protocol is fully implemented.

jmichel80 commented 6 years ago

On a related note, check also def make_pert_file in https://github.com/CCPBioSim/fesetup/blob/dev/mutate/topol/pertfile.py

I don't know how BSS merge molecules handle this, but there was some complicated logic in FESetup to make sure that dihedrals and impropers were correctly written in the pert file. The issue is that the same dihedral or improper may be defined with different atom indices in the respective end-states, and care must be taken to avoid accidentally duplicating forcefield terms.

There are also some rules around the values of bond/angles or dihedral/impropers for terms involving dummy atoms that could be controlled with keyword arguments.

lohedges commented 6 years ago

I believe I handle this correctly in both the Gromacs and pert file writer. Since all of the terms are stored as Sire objects, with corresponding BondIDs, AngleIDs, and DihedralDs, it's easy to check for matching terms in the end-states, e.g. by checking the the ID and its mirror.

I'll take a look through the FESetup logic to see if I'm missing anything. For the Gromacs writer, I believe that my procedure is consistent with everything from the latest manual.