search

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

"forcefield" property #147

Open SofiaBariami opened 4 years ago

SofiaBariami commented 4 years ago

Hello, I am not sure whether I should post this question here or at the BioSimSpace repository, but I am trying to understand how we can set the "forcefield" property that is required when we want to merge two molecules. The standard process to set properties is to create an editable molecule and atoms and then set all the properties separately, as we can see in amber.cpp, but I can't find anything on forcefield. At the BSS code, we can see at _molecule.py that this property is required and it should also have a specific Amber format. Does anyone have any thoughts about that? Thanks a lot!

lohedges commented 4 years ago

Hi @SofiaBariami, I agree that is kind of a dual Sire/BioimSpace question but I've transferred this issue over here since the issue (for now) is likely specific to the creation of merged molecules in BioSimSpace.

The forcefield property of a Sire molecule is an MMDetail object. For an AMBER-style system it might look something like this. (Here I'm querying the underlying Sire object of a molecule in BioSimSpace that was loaded from AMBER format files.)

molecule._sire_object.property("forcefield")
MM ForceField{ amber::ff,
               combining_rules = arithmetic,
               1-4 scaling = 0.833333, 0.5,
               nonbonded = coulomb, lj,
               bond = harmonic, angle = harmonic,
               dihedral = cosine }

You can see how this property is set in the new AmberPrm parser here, and similarly for the GroTop parser here. It looks like this property is not used or set in the old amber.cpp parser that you refer to.

In BioSimSpace, we currently only support merges between two molecules that have compatible AMBER-style forcefields. The MMDetail object provides a convenient way of checking for compatibility and could be extended to support additional forcefield types in future.

The following code snippet shows how the MMDetail object is used within the private _merge method of the BioSimSpace Molecule class to check whether it is possible to merge two molecules based on their forcefield property.

# Get the user name for the "forcefield" property.
ff0 = inv_property_map0.get("forcefield", "forcefield")
ff1 = inv_property_map1.get("forcefield", "forcefield")

# Force field information is missing.
if not molecule0.hasProperty(ff0):
    raise _IncompatibleError("Cannot determine 'forcefield' of 'molecule0'!")
if not molecule1.hasProperty(ff1):
    raise _IncompatibleError("Cannot determine 'forcefield' of 'molecule1'!")
# The force fields are incompatible.
if not molecule0.property(ff0).isCompatibleWith(molecule1.property(ff1)):
    raise _IncompatibleError("Cannot merge molecules with incompatible force fields!")

I hope this makes sense. I'm happy to discuss more if you need a clearer explanation. If you think that you might need to extend the MMDetail class to support any additional forcefields that you are implementing then I should be able to offer some guidance.

SofiaBariami commented 4 years ago

Hi Lester, thank you very much. This is very helpful. I am indeed trying to use BioSimSpace to work with molecules that are parameterised with another forcefield, and since the forcefield property was not set, python was complaining. Maybe extending the MMDetail class is the best way to proceed with this, especially if we are planning to use QuBe in the future as well. I have managed to overcome some of the errors in hacky ways. For example I substituted ff.electrostatic14ScaleFactor() and ff.vdw14ScaleFactor() of the _merge function with the actual values of the scale factors, but this practice is not the best, since the code is complaining in other steps of the process as well. For example when I am calling _toPertFile() I am getting an Exception 'SireBase::missing_property' because of the same thing: 

--> 709   if not self._sire_object.property("forcefield0").isAmberStyle():
     710       raise _IncompatibleError("Can only write perturbation files for AMBER style force fields.")

@jmichel80 what do you think? Should we extend MMDetail, or it is indeed too much of a detail, so find another way to overcome the issue?

jmichel80 commented 4 years ago

hi @SofiaBariami

Can you post the code you use to load a pair of QUBE molecules to be merged and sample input files (two small ligands in vacuum) ? Post also code for the AMBER ff use case. I think this happens because you are not using the same amber parser that @lohedges wrote for BioSimSpace. We haven't ported this to somd as it requires some additional code changes but we may have to consider this. It could also help solve the performance issues you have seen when loading large molecules parameterised with QUBE

BW

Julien

Dr. Julien Michel, Senior Lecturer Room 263, School of Chemistry University of Edinburgh David Brewster road Edinburgh, EH9 3FJ United Kingdom phone: +44 (0)131 650 4797 http://www.julienmichel.net/

On Mon, Mar 9, 2020 at 10:07 PM Sofia Bariami notifications@github.com<mailto:notifications@github.com> wrote:

Hi Lester, thank you very much. This is very helpful. I am indeed trying to use BioSimSpace to work with molecules that are parameterised with another forcefield, and since the forcefield property was not set, python was complaining. Maybe extending the MMDetail class is the best way to proceed with this, especially if we are planning to use QuBe in the future as well. I have managed to overcome some of the errors in hacky ways. For example I substituted ff.electrostatic14ScaleFactor() and ff.vdw14ScaleFactor() of the _merge functionhttps://github.com/michellab/BioSimSpace/blob/devel/python/BioSimSpace/_SireWrappers/_molecule.py#L1942 with the actual values of the scale factors, but this practice is not the best, since the code is complaining in other steps of the process as well. For example when I am calling _toPertFile() I am getting an Exception 'SireBase::missing_property' because of the same thing:

--> 709 if not self._sire_object.property("forcefield0").isAmberStyle(): 710 raise _IncompatibleError("Can only write perturbation files for AMBER style force fields.")

@jmichel80https://github.com/jmichel80 what do you think? Should we extend MMDetail, or it is indeed too much of a detail, so find another way to overcome the issue?

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHubhttps://github.com/michellab/BioSimSpace/issues/147?email_source=notifications&email_token=ACZN3ZD5LISH3DDECI4CLLDRGVSAZA5CNFSM4LEQDBXKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOEOJHX2I#issuecomment-596802537, or unsubscribehttps://github.com/notifications/unsubscribe-auth/ACZN3ZEGX7Y5UAZ3C6KC7LDRGVSAZANCNFSM4LEQDBXA.

The University of Edinburgh is a charitable body, registered in Scotland, with registration number SC005336.

lohedges commented 4 years ago

Much as I'd like to take credit for it, the updated AMBER parser is all @chryswoods work :-)

SofiaBariami commented 4 years ago

you can find the files and the script here: https://github.com/SofiaBariami/qube_project/tree/master/merge_molecules

@jmichel80 when you say to post code for the Amber ff case, do you mean to refer to the BSS code that we use? If so, we are using the following main functions:

jmichel80 commented 4 years ago

Hi Sofia,

Can you post a detailed description on github of the issue. Show what code from BioSimSpace you tried to use, what error messages you got. What modifications you think are needed to the various BioSimSpace functions. From quickly glancing at the code I get the impression you have copied/pasted several functions from BioSimSpace and then modified the code to make it work with your inputs. It would be better to identify what needs to be changed to your input molecules to make them compatible with the existing BioSimSpace code, and to make minimal modifications to BioSimSpace.

When I meant show the Amber ff use case, I meant to describe in the github issue exactly how you call various functions in BioSimSpace to produce a set of somd input files from a set of input prm7/rst7 files. This so we understand better what is it that makes the qube input incompatible whereas the amber input files work.

Dr. Julien Michel, Senior Lecturer Room 263, School of Chemistry University of Edinburgh David Brewster road Edinburgh, EH9 3FJ United Kingdom phone: +44 (0)131 650 4797 http://www.julienmichel.net/

On Tue, Mar 10, 2020 at 10:09 AM Sofia Bariami notifications@github.com<mailto:notifications@github.com> wrote:

you can find the files and the script here: https://github.com/SofiaBariami/qube_project/tree/master/merge_molecules

@jmichel80https://github.com/jmichel80 when you say to post code for the Amber ff case, do you mean to refer to the BSS code that we use? If so, we are using the following main functions:

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHubhttps://github.com/michellab/BioSimSpace/issues/147?email_source=notifications&email_token=ACZN3ZFDEX6FMZZSBBHKS7LRGYGT3A5CNFSM4LEQDBXKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOEOKZC6A#issuecomment-597004664, or unsubscribehttps://github.com/notifications/unsubscribe-auth/ACZN3ZA4LWXLZMOS6VHGWNLRGYGT3ANCNFSM4LEQDBXA.

The University of Edinburgh is a charitable body, registered in Scotland, with registration number SC005336.

SofiaBariami commented 4 years ago

Thanks for the clarifications Julien.

  1. When it comes to Amber, two sets of prm7/rst7 files are read with readMolecules(), to form two BioSimSpace._SireWrappers._system.System. From these system, we can extract the molecules that are of type BioSimSpace._SireWrappers._molecule.Molecule. Next we need to use matchAtoms() that finds mappings between atom indices in molecule 0 to those in molecule 1. The input of this function must be of type Molecule <BioSimSpace._SireWrappers.Molecule>. The first problem here is that when reading pdb/xml files, the output molecule will be a Sire.Mol._Mol.Molecule. However, with this is solved easily with BSS._SireWrappers.Molecule(), that turns a Sire molecule into a BSS molecule. This was the main problem that made me change the functions, but with this fix, we can use the BSS functions as they are.
  2. The second problem is the one with that we discussed with Lester above. During the merge of two molecules, BSS checks whether their forcefields are compatible. When we are using Amber files, the forcefield property is set, and looks like this: 
    In [20]: mol1._sire_object.property("forcefield")                                                           
Out[20]: 
MM ForceField{ amber::ff,
           combining_rules = arithmetic,
           1-4 scaling = 0.833333, 0.5,
           nonbonded = coulomb, lj,
           bond = harmonic, angle = harmonic,
           dihedral = cosine }

    However, the molecule loaded from xml/pdb, does not have this property at all, that results to an Exception 'SireBase::missing_property'. In order to solve this problem we should modify mmdetail and amberprm.cpp to accept input from xml/pdb with the qube forcefield. Since there are separate files for the file formats in Sire.IO (grotop, amberprm) I am not sure whether we should make new files for the xml/pdb files. The easy, temporary fix that I did, in order to proceed with my project quickly, is that I just commended out the lines where the forcefield property was checked. I know that this doesn't solve the issue, but I decided to ignore this bit until we come up with a solution, without letting it interfere with the progress of the project.

lohedges commented 4 years ago

I'm not completely sure how you are currently constructing a Sire/BioSimSpace system parameterised with the QuBe forcefield, but this is what I think would be the ideal solution...

Hopefully you are using the Sire PDB parser, then adding properties to the System object with the information that you parse from the XML file. (Although, perhaps you are doing this the other way around.) This is how the Sire parsers are intended to work: one parser constructs the initial system (using the molecular topology), then other parsers can add additional information to the existing System. You can read more details about this in the BioSimSpace paper here.

Ideally there would be a new Sire parser called QubeXML (or similar) that could parse the forcefield information from XML file. This would create its own MMDetail object which would be added as a property of each molecule in the System. If the forcefield can have terms that we don't currently support then we'd also need to update MMDetail to reflect this.

Obviously all of the above is quite involved and requires a lot of knowledge of Sire.

lohedges commented 4 years ago

Hi @SofiaBariami, I noticed that you recently pushed to the feature-vSites branch. Is the feature_virtual-sites branch now redundant? If so, I'll delete it since I'd like to avoid a proliferation of stale branches.

Cheers.

SofiaBariami commented 4 years ago

hi Lester, yes, can you delete the feature_virtual-sites branch please? I left it to delete = it myself later on, but I don't need it anymore