michellab / BioSimSpace_examples

GNU General Public License v2.0
0 stars 0 forks source link

dummy coordinates in mapping code #11

Open jmichel80 opened 6 years ago

jmichel80 commented 6 years ago

Documenting a difference in behaviour between the MCS BSS and FESetup implementations.

scripts/test_mapping_methane_toluene.py

generates a merged molecule wth a topology comparable to FESetup (with a local copy of findmcs.cpp modified to match hydrogens to heavy atoms). However the coordinates of the dummy atoms in merged0.pdb (the ring atoms, apart from the hydrogen/carbon atom) are those from the aligned toluene.pdb molecule. This means that the bonded geometry around the dummy/non-dummy atoms is distorted. See picture below, showing long bond lengths and low angle for the side of the benzyl ring bonded to the methyl.

image

By contrast here is the morph input generated by FESetuo for methane->toluene

image

The difference is that in FESetup the Cartesian coordinates of the dummy atoms are worked out by mixing the coordinates of the mapped atoms in topology0 (methane), and the internal coordinates of the unmapped atoms in topology1 (toluene). In this case here in effect the benzyl ring is translated towards the methyl group so the methane molecule geometry is not distorted.

More testing is needed to determine whether one implementation is better than the other. To summarise:

The BSS implementation does not update coordinates, but generates initial conformations that are potentially highly strained The FESetup implementation generates initial conformations with low intramolecular energetics, but potentially changes intermolecular energetics.

lohedges commented 6 years ago

Thanks, this is good to know. It would be very easy to give the user the ability to specify how the dummy coordinates are set as a keyword option when calling the merge function. Do you know where in FESetup the coordinates are mixed so I can take a look at how its implemented.

As a related question, do you know how the various codes, e.g. Gromacs, handle the dummy atoms coordinates? That is, do they do any mixing for you, or do they just take the bare coordinates from the input file and run with those? I guess we will generally be doing some equilibration at the start of each lambda stage, so it's possible it doesn't matter too much anyway. (Apart from very strained configurations, perhaps.)

jmichel80 commented 6 years ago

You can read the source code of the FESetup implementation here

https://github.com/CCPBioSim/fesetup/blob/dev/mutate/util.py

specifically def dummy_coords line 882

I suspect all MD codes take the input coordinates and run with those. Equilibration may fail with distorted input structures, but I suspect energy minimisation followed by equilibration should work in many cases. So yes it may not matter too much and I suggest you don't work on more options to compute dummy coordinates unless we realise this is needed.