search

MiaoLab20 / gamd-openmm

Gaussian Accelerated Molecular Dynamics (GaMD) is a computational method for both unconstrained enhanced sampling and free energy calculations of biomolecules. It works by adding a harmonic boost potential to smooth the potential energy surface and reduce system energy barriers. http://miao.compbio.ku.edu/GaMD/
MIT License
65 stars 18 forks source link

GaMD crashing for CHARMM36 force field files #22

Open mdpoleto opened 2 years ago

mdpoleto commented 2 years ago

Hi all,

Our lab has tried to use the GaMD-OpenMM code to simulate CHARMM36 systems but the stock code seems to be not working for them. @hmmichel has been working on this as well.

The first issue we had was how the code was setup to read, parse and build C36 systems using OpenMM functions. Those were generally easy fixes, which I am making a PR for. Basically, psf.createSystem() requires the box vectors to be assigned beforehand. To that end, I created a new tag in the config.py to make possible for the user to define the box vector within the XML file. The last modification was to make the XML parser more compatible with CHARMM-GUI inputs, in which the toppar files are listed in a toppar.str file and a small python function reads the file and grep the list of topology files to be read. I am including a tar ball to help reproducing the errors mentioned here.

c36_gamd_testfiles.zip

With the system being read and built correctly (we checked the forces of each force group), we attempted to simulate an alanine dipeptide in order to reproduce the data from the original paper. The conventional MD steps work just fine, but the system crashes in the first steps of GaMD equilibration. I could not pin-point the problem yet, but I believe it has something to do to the topology parsing and how the force groups as being assigned. When we break down the force groups of each system, we obtain:

FF14SB: 0 HarmonicBondForce 1 HarmonicAngleForce 2 PeriodicTorsionForce 3 NonbondedForce 4 CMMotionRemover 5 MonteCarloBarostat

FF19SB: 0 HarmonicBondForce 1 HarmonicAngleForce 2 PeriodicTorsionForce 3 CMAPTorsionForce 4 NonbondedForce 5 CMMotionRemover 6 MonteCarloBarostat

CHARMM36: 0 HarmonicBondForce 1 HarmonicAngleForce 2 HarmonicBondForce (Urey-Bradley) 3 PeriodicTorsionForce 4 CustomTorsionForce (Impropers) 5 CMAPTorsionForce 6 NonbondedForce 7 CustomNonbondedForce (NBFIX) 8 CMMotionRemover 9 MonteCarloBarostat

Maybe the group assignment is not handling those differences? I am happy to help/test any solutions.

mdpoleto commented 2 years ago

I made a PR (#23) with the necessary fixes to be able to read in and build a C36 system. Any suggestions are very welcomed!

dyelar commented 1 year ago

Right now, we set the PeriodicTorsionForce to force group 2 and the NonbondedForce to force group 1, when we set those values. The values are set by name and it appears the names are the same. In your example, you are doing a lower-dihedral boost, so it would only set the PeriodicTorsionForce to force group 2 for acting upon it and the NonbondedForce would not be set.

Were you expecting it to boost on the CustomTorsionForce instead?

mdpoleto commented 1 year ago

No, CustomTorsionForce is used for Impropers in CHARMM. My goal was to boost torsions with PeriodicTorsionForce. Since we are on the topic, shouldn't CMAPTorsionForce be included (if available) in the boosted group along with PeriodicTorsioForce? I mean, if a torsional energy is calculated by both dihedral and CMAP terms, I would imagine these 2 groups would be assigned to a specific group to be boosted. I think the same rationale could be applied to NonbondedForce and CustomNonbondedForce (NBFIXes).

mdpoleto commented 1 year ago

@lvotapka mentioned something via email that I believe solves the issue.

In OpenMM, "amber_file_parser.py" parses the topology and assigns all ForceGroups to group0. However, for CHARMM, "charmmpsffile.py" assigns each ForceGroup to their own group (bonds to 0, angles to 1, Urey-Bradley to 2, etc). Correct me if I am wrong here, but as far as I understood, gamd-openmm requires all ForceGroups to be in group0 and only the boosted groups to be separated from the rest.

Using the inputs in my example, I manually set all ForceGroups to group0 and the simulation does not blow up anymore on stage3 (gamd equilibration). I will test this fix further and I let you know how well it works.

GMdSilva commented 1 year ago

Hi everyone,

Were there any developments on this issue? I'd like to use GaMD-OpenMM with a system built in the CHARMM36 forcefield, but it's not clear to me how to proceed.

Thanks!

MiaoLab20 commented 1 year ago

I believe we have a solution for this. @dyelar and @mdpoleto, could you probably advise on this? Thanks!

mdpoleto commented 1 year ago

@MiaoLab20 and @GMdSilva, the current code does not crash when using CHARMM36 input files anymore. However, it also does not consider CMAP and NBFIXes terms when boosting torsional and nonbonded potentials, so simulations using CHARMM36 or AmberFF19SB input files are most likely to be incorrect.

I believe @dyelar was testing a new version (that can be found in the devel branch) but I do not know if it is ready to be merged into the main branch. He might be able to provide more details on this.

dyelar commented 1 year ago

I have code for this in the test branch. I was updating the necessary testing code for running validation. Everything got put on hold, when my son broke his neck, needed to have surgery, and is needing care. I've had to be very specific about the order in which I get things done across my different projects, as I either had no work time or only a few hours of work time per day. I already had on my schedule for tomorrow to be back working on this. My apologies for the delay.

On Thu, May 4, 2023 at 4:11 PM Marcelo D. Polêto @.***> wrote:

@MiaoLab20 https://github.com/MiaoLab20 and @GMdSilva https://github.com/GMdSilva, the current code does not crash when using CHARMM36 input files anymore. However, it also does not consider CMAP and NBFIXes terms when boosting torsional and nonbonded potentials, so simulations using CHARMM36 or AmberFF19SB input files are most likely to be incorrect.

I believe @dyelar https://github.com/dyelar was testing a new version (that can be found in the devel branch) but I do not know if it is ready to be merged into the main branch. He might be able to provide more details on this.

— Reply to this email directly, view it on GitHub https://github.com/MiaoLab20/gamd-openmm/issues/22#issuecomment-1535417474, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABALA2LYQRP22XYFKVMTIMTXEQLOHANCNFSM6AAAAAAQCT32DA . You are receiving this because you were mentioned.Message ID: @.***>

dyelar commented 1 year ago

@mdpoleto @GMdSilva @MiaoLab20 @lvotapka

I split out the code from the other code I was working with into a new branch. The tests looked good for the changes that were made to include the CustomNonbondedForce and the CMAPTorsionForce. The code has been committed and merged into the main branch. I've tagged with 0.9.2.

I'll leave this issue open for the next week, in case something comes up when you attempt to utilize. After that, I'll close this issue.