search

mosdef-hub / mbuild

A hierarchical, component based molecule builder
https://mbuild.mosdef.org
Other
174 stars 81 forks source link

energy_minimize not minimizing ethane bonds to correct length using foyer TraPPE-UA FF #940

Closed bc118 closed 3 years ago

bc118 commented 3 years ago

Bug summary

The energy_minimize function is not minimizing ethane bonds to correct length using foyer TraPPE-UA FF, even when starting from a mol2 file at nearly a perfect bond distance. It actually moves them to a bond distance quite far away from the equilibrium bond distance. I'm assuming this is going to be an issue for larger molecules as well.

This is especially an issue for MC engines, as the only time they sample bonds is during a regrowth move. Therefore, it may or may not be fixed in the simulation, if the simulation will even start with these high energy bonds. Although, this is likely an issue for some MD engines, at least at the start of the simulation.

Code to reproduce the behavior

minimization_error.zip

See the attached code (minimazation_error_trappe.py) and the attached excel spreadsheet (Minimzation_excel.xlsx). To compare bond distances manually, you will need to calculate the between the 2 carbon atoms in the system:

Software versions

justinGilmer commented 3 years ago

I think we might be able to leverage the energy_minimize(scale_bonds) function argument to scale the bond force constant to 0. That should keep the bonds at their starting distance, although I have not used this argument before.

justinGilmer commented 3 years ago

I realized I used the incorrect scaling factor above. We should scale the bond constants to large values to hold the bond lengths in place.

bc118 commented 3 years ago

Sorry, I do not really understand what the above comment means. Does it mean that is the source of the error?

justinGilmer commented 3 years ago

I’m not entirely sure why the bonds are being stretched from their equilibrium positions that much. But I think we can massively scale the bond constants to incur a massive energy penalty of they deviate from their equilibrium positions (that they essentially start in). That should still let the angles and dihedrals settle into a decently minimized orientation, barring no major metastability issues.

justinGilmer commented 3 years ago

Actually, i just found a working fix, we can add additional support for the other options in parmed.structure.createSystem to pass in simtk.openmm.app.AllBonds

I tested it locally and it does appear to work as expected.