GaMD crashing for CHARMM36 force field files

[mdpoleto]

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]

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]

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]

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]

@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]

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]

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

[mdpoleto]

@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]

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 <#22 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ABALA2LYQRP22XYFKVMTIMTXEQLOHANCNFSM6AAAAAAQCT32DA> . You are receiving this because you were mentioned.Message ID: ***@***.***>

[dyelar]

@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.