Looks like new XML file is properly recognizing the ACE-ARG-NME patch residue without any problem. Thanks for your work. Really appreciate it!

I see in the new xml file the lines for the hydrogen. I am still getting a particle position = nan error, but the stack trace is different than reported in #3828

Traceback (most recent call last):
  File "/data/B/bandrews/OpenMM/GGG_CHARMM_Drude_2019_SWM4-NDP_water/openmm_setup.py", line 75, in <module>
    simulation.step(600000000)
  File "/home/bandrews/miniconda3/envs/openMM/lib/python3.9/site-packages/openmm/app/simulation.py", line 134, in step
    self._simulate(endStep=self.currentStep+steps)
  File "/home/bandrews/miniconda3/envs/openMM/lib/python3.9/site-packages/openmm/app/simulation.py", line 236, in _simulate
    self._generate_reports(wrapped, True)
  File "/home/bandrews/miniconda3/envs/openMM/lib/python3.9/site-packages/openmm/app/simulation.py", line 258, in _generate_reports
    reporter.report(self, state)
  File "/home/bandrews/miniconda3/envs/openMM/lib/python3.9/site-packages/openmm/app/dcdreporter.py", line 107, in report
    self._dcd.writeModel(state.getPositions(), periodicBoxVectors=state.getPeriodicBoxVectors())
  File "/home/bandrews/miniconda3/envs/openMM/lib/python3.9/site-packages/openmm/app/dcdfile.py", line 125, in writeModel
    raise ValueError('Particle position is NaN')
ValueError: Particle position is NaN

I could look into which particle is producing the NaN value. But that will take some time. Thanks for making the change to the forcefield! I think the stack trace changing is at least encouraging.

Brian

Does it become nan within the first few steps, or does it take a while? If the latter, it probably is something different. The FAQ has advice on debugging problems of this sort.

It has varied between 50K and 250K time steps.

I am finding my system with a negatively charged C terminus also results in a particle coordinate is nan but after 1 million time steps, so I have not realized this system was also breaking with the same stack trace. Reducing time steps or further energy minimization has not resolved this issue so far. I will continue trouble shooting. Thanks again for your help.

@andrewsb8 Hi? I also had a similar issue while running a single amino acid residue (but not GLY) with negatively charged C residue. From here, this is a bit my guess about this issue. (I did not try it yet) For charmm_polar_2019.xml file, <AllowPatch name="CNEG"/> and <AllowPatch name="CTEG"/> are used for C-termius residue. If you check the original .str file for such force field in here (file: drude_toppar_2019.str), below is the definition of CNEG and CTEG residue.

PRES CNEG          0.000 ! protonated C-terminus for polypeptide with terminal glycine
                         ! Yes, the patch name is not ideal; the total charge is zero
!
!    OT2-HT
!    /
! --C
!    \\
!    OT1
!
DELETE ATOM O
DELETE ATOM DO
DELETE ATOM LPOA
DELETE ATOM LPOB

PRES CTEG         -1.000 ! C-terminus for polypeptide terminating with a glycine
!
!     OT2(-)
!    /
! --C
!    \\
!     OT1
!
DELETE ATOM O
DELETE ATOM DO
DELETE ATOM LPOA
DELETE ATOM LPOB

I don't know the specific detail about force field parameterization but at least it shows that these residues suitable for polypeptide terminating with a glycine. Well for your single amino acid case, there are other C-terminus residues used for single amino acid.

PRES CNES          0.000 ! protonated C-terminus for single amino acid
! Note: with ACE use CNEU
!
!    OT2-HT
!    /
! --C
!    \\
!    OT1
!
DELETE ATOM O
DELETE ATOM DO
DELETE ATOM LPOA
DELETE ATOM LPOB

PRES CTES         -1.000 ! C-terminus for single amino acid
! Note: with ACE use CTER
!
!     OT2(-)
!    /
! --C
!    \\
!     OT1
!
DELETE ATOM O
DELETE ATOM DO
DELETE ATOM LPOA
DELETE ATOM LPOB

In current charmm_polar_2019.xml file, forcefield.createsystem does not automatically check whether C-terminus patch is for single amino acid or polypeptide, and since both [CTES,CNES] and [CTEG,CNEG] have same topology information, only one of them is defined in xml file. I'm not sure if it would work or not but you may add parameters for [CTES,CNES] patch, change the <AllowPatch> from [CTEG,CNEG] to [CTES,CNES], and run the simulation for trobleshooting.

This is kind of an impedance mismatch between CHARMM and OpenMM. In CHARMM, patches are a tool for modelling. The user tells it, "Apply this patch to this residue," and it modifies the residue as specified. It's up to the user to read the comments and figure out what patch makes sense. It doesn't specify programmatically that one patch is only intended for single residues and a different one for polypeptides.

In OpenMM, patches are a tool for parametrization. The force field receives a model that as already been modified in some way, and it has to try to figure out what patch was used to do it so it can pick correct parameters. For all it knows, you could have applied CTES to a polypeptide or CTEG to a single amino acid. That isn't how they were intended to be used, but that intent isn't recorded in any machine readable way.

I doubt this would cause the simulation to blow up though. The two patches are very similar. As far as the simulation code is concerned, both versions should work equally well. One may give a more realistic description of how the charge is distributed in a particular real molecule, but the simulation code doesn't care how physically accurate your model is, just whether it's numerically well behaved.

@smlee99 Thank you for the suggestion! I am simulating triglycine so the issue between patches for amino acids and polypeptides shouldn't apply here since I am simulating the latter. I also agree with Dr. Eastman that either patch should not cause the system to blow up.