OpenMM and Plumed RMSD Restraint

[jwan...@stanford.edu]

Hi Plumed Community!

My ultimate goal is to conduct targeted MD using RMSD as the collective variable and the RESTRAINT function. But I've been running into issues in which the particle coordinates turn to NA whenever I attempt this. I had the chance to conduct a series of experiments to test my system using the Plumed OpenMM plugin. 

1) I used RESTRAINT on a number of torsional angles and then examined the results in VMD. The results show that the external force was applied on the expected atoms, verifying that my atom numbering is correct in the reference PDB. This works fine!

2) I tried the RESTRAINT function on a dozen bond distances in my system. This also yielded the expected results, with the correct atoms being biased. This also works fine!

3) I computed RMSD without biasing the system and got an expected result. And this also works fine!

However, when I attempt to bias the system using RMSD as the collective variable, I inevitably end up with an error: particle coordinates turn NA. I'm confident that my atoms are numbered correctly with the correct OCCUPANCY and BETA factors set in the reference PDB. 

Given this situation, I have 2 question: 1) have other Plumed users had success biasing RMSD as a collective variable on an OpenMM platform?, 2) if so, what common mistakes could lead to a NA atom coordinate issue?

[Giovanni Bussi]

See this post: https://groups.google.com/d/msgid/plumed-users/97478daf-5a18-40e4-ab01-c646d8d46fcb%40googlegroups.com?utm_medium=email&utm_source=footer

The problem is with your pdb file. You should convert it to unix newlines with dos2unix.

I fixed the parser yesterday (https://github.com/plumed/plumed2/issues/223), so the problem should not appear in the future hopefully.

Giovanni

--
You received this message because you are subscribed to the Google Groups "PLUMED users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to plumed-users+unsubscribe@googlegroups.com.
To post to this group, send email to plumed...@googlegroups.com.
Visit this group at https://groups.google.com/group/plumed-users.
To view this discussion on the web visit https://groups.google.com/d/msgid/plumed-users/9d6b2bd8-9972-45a5-afdb-01fa923367ce%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Jason Ku Wang]

Thanks for the suggestion - I went ahead and made sure newlines were in unix form and the PDB reference file loads in fine. But after a few timesteps, I still receive the same particle NA error as before: 

"Traceback (most recent call last):

  File "targeted_md.py", line 100, in <module>

    sim.step(50000)

  File "PATH/simulation.py", line 132, in step

    self._simulate(endStep=self.currentStep+steps)

  File "PATH/simulation.py", line 194, in _simulate

    self.integrator.step(10) # Only take 10 steps at a time, to give Python more chances to respond to a control-c.

  File "PATH/openmm.py", line 4731, in step

    return _openmm.LangevinIntegrator_step(self, steps)

Exception: Particle coordinate is nan"

For further reference, here are the first few lines of my COLVAR output file: 

#! FIELDS time restraint.bias 

#! SET min_t 

#! SET max_t

 0.004000 25.691915

 0.008000 25.691915

 0.012000 nan

 0.016000 nan

 0.020000 nan

The remainder of the COLVAR prints nan. 

Thanks for the help once again - really looking forward to getting this portion of Plumed to work :)!

ᐧ

--
You received this message because you are subscribed to a topic in the Google Groups "PLUMED users" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/plumed-users/S0qbsPgGHjY/unsubscribe.
To unsubscribe from this group and all its topics, send an email to plumed-users+unsubscribe@googlegroups.com.

To post to this group, send email to plumed...@googlegroups.com.
Visit this group at https://groups.google.com/group/plumed-users.

To view this discussion on the web visit https://groups.google.com/d/msgid/plumed-users/CAPLm8ZKVunHFnn2z1Ycr%3DeJV2d%2BZFPMFd1verpCzfUq4Uqp84A%40mail.gmail.com.

For more options, visit https://groups.google.com/d/optout.

--

Jason Wang

Stanford University, Class of 2018

jwan...@stanford.edu

[Giovanni Bussi]

Perhaps the KAPPA is too large. Try to increase gradually to see if at some point the simulation has problems.

Giovanni

To view this discussion on the web visit https://groups.google.com/d/msgid/plumed-users/CAO%3DdQNg58W5ZKX-JfmcMdKU%3DiBggQjdthPGBKw0LLHAMDB%2B8PA%40mail.gmail.com.

[Jason Ku Wang]

Hi Professor Bussi, 

Thanks for that suggestion - I just tried running the targeted MD with KAPPAs ranging from 0.001 to 2092. All KAPPAs still result in the same NA atom coordinate issue as well as "ValueError: Energy is NaN" error.

To dig deeper, I tried printing out some metrics. In particular, we see that once the Plumed RMSD Restraint is enacted during the first timestep, the kinetic energy explodes. I've included the first printout of a control simulation of the same system without a Plumed restraint as a reference. 

"Step","Time (ps)","Potential Energy (kcal/mol)","Kinetic Energy (kcal/mol)","Total Energy (kcal/mol)","Temperature (K)","Box Volume (angstrom**3)"

With Plumed Force: 1,0.004,-402081.389505,1.05313646317e+12,1.05313606109e+12,3795867991.35,1366599.66129

Without Plumed Force:

1,0.004,-455816.488053,85885.7311361,-369930.756917,309.561874586,1365753.36761

Given these results, I wonder if this issue is a result of a large system (my system has 138108 atoms), Plumed, or the OpenMM Plugin to Plumed that I am using.

Best,

Jason

ᐧ

To view this discussion on the web visit https://groups.google.com/d/msgid/plumed-users/CAPLm8ZKPDjjKCj-pP3Am20Tf5YFkA3s9M%3DyrH_ALLdjKLbEDYA%40mail.gmail.com.

For more options, visit https://groups.google.com/d/optout.

[Jason Ku Wang]

It looks like I've pinpointed an exact issue but not the cause:

The energy inevitably explodes after a few steps no matter how small I set the KAPPA, and how long I equilibrate. Has Plumed's RMSD collective variable been shown to work on large systems (e.g. 138108 atoms)?

Thanks once again,

Jason

ᐧ

[Giovanni Bussi]

PLUMED is routinely used on large systems, and using RMSD as a variable for metadynamics or restrained MD.

I opened your PDB file again and it looks not correctly aligned. Please check if it follows PDB specification for all the fields

http://www.wwpdb.org/documentation/file-format-content/format33/sect9.html#ATOM

Giovanni

To view this discussion on the web visit https://groups.google.com/d/msgid/plumed-users/CAO%3DdQNiEt%2B-RRMOje121hbQDLvpUtxKX9CspUJkz05HEoDS3_A%40mail.gmail.com.

[jwan...@stanford.edu]

That's likely it! Thanks so much!

Additionally, it seems that Plumed takes in PDB formatted reference files. In the case of systems that are >99,999 atoms (PDB format only allows numbering up to 99,999), how would we go about loading in such a reference?

[Giovanni Bussi]

There's no way currently. Notice that computing an RMSD with 100000 atoms probably does not make very much sense. Anyway, we should tackle this problem sooner or later, so I will open a github issue not to

forget it

Giovanni

--

You received this message because you are subscribed to the Google Groups "PLUMED users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to plumed-users+unsubscribe@googlegroups.com.
To post to this group, send email to plumed...@googlegroups.com.
Visit this group at https://groups.google.com/group/plumed-users.

To view this discussion on the web visit https://groups.google.com/d/msgid/plumed-users/e34899e0-76e8-4961-b093-dd4eab446d75%40googlegroups.com.

[jwan...@stanford.edu]

Hi Giovanni!

Thanks to your help, the targeted MD portion works great with RMSD restraints and moving restraints. But we'd still love to be able to read in atoms indices that are numbered >99,999. For the case of a large system in which water and ions are numbered first, the protein atoms are numbered >99,999. Perhaps a fix could be, being able to read atom indices in hexadecimal?

Would such a fix be simple to do/would it be possible in the near future?

Best,

Jason

On Friday, April 7, 2017 at 2:14:23 AM UTC-7, jwan...@stanford.edu wrote:

[Giovanni Bussi]

There's an issue open on github: https://github.com/plumed/plumed2/issues/226

I am personally not sure about the hex format for the reasons explained on github. Anyway, you can easily implement it in your local version. The conversion is made in file src/tools/PDB.cpp, around line 127.

Just replace the following line:

      Tools::convert(serial,a);

with something like

  unsigned tmp;

  sscanf(serial.c_str(),"%whateverformatyoulike",&tmp);

  a.setSerial(tmp);

serial is a string made of 5 characters.

Giovanni

--

You received this message because you are subscribed to the Google Groups "PLUMED users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to plumed-users+unsubscribe@googlegroups.com.
To post to this group, send email to plumed...@googlegroups.com.
Visit this group at https://groups.google.com/group/plumed-users.

To view this discussion on the web visit https://groups.google.com/d/msgid/plumed-users/e3bf1fb8-4982-4676-9c39-8fa0b2dd884b%40googlegroups.com.