CHARMM Development Project Forums User Discussion & Questions Molecular Dynamics Size of Water Box w/ PBCs

Recently, as I attempted to perform simulations analogous to previous ones I've performed but with a slightly different set up, I encountered an unexpected error. I was attempting to "restart" dynamics using the spatial coordinates of a point along my first directory and rescaled velocities, appropriate for the temperature of the simulations (i.e., reading in those spatial coordinates using the "dyna strt" command as opposed to "dyna restart").


As I had not had trouble doing this with previous setups, I was surprised by this error.

...especially because this error occurred seemingly before any dynamics had actually begun (and definitely before anything was written to either the restart or the trajectory files). There were any number of errors like this (see attached output file--er, no, it's too large) and then charmm crashes.

From my reading in these forums, I understand that a SHAKE error is often the symptom of a more fundamental problem. Looking at the atoms listed in the errors, they all seem to be waters, which is curious because I've chosen to do the water images by residue (now you can see the attached input file--so no water should be split between sides of the boundary. Also, when I look at the actual coordinates in the .crd file these atoms are an appropriate distance apart!)

So, the source of the error would seem to be elsewhere. Upon examining the original trajectory, I noticed that the protein drifts out of the box. At first, I thought, aha, despite what SHAKE is saying, this must be the problem, but then I read other forum threads (including this one), which suggests to me that this is not in and of itself problematical...

However, this unexpected error has raised questions in my mind about the size of the water box. The buffer I've given is smaller than the ctofnb distance (8 vs 12). I've checked (as suggested here) to make sure that the protein never 'sees' its own image, and it does not. However, there are clearly solvent molecules that 'see' both the protein and its image, and I wonder if this is what gives the protein its velocity (with ntrfrq set to 500--or is that an indication that this should be a lower number??). The discussion of how much solvent is needed on the charmmtutorial.org website makes it sound not a completely cut and dry issue. I am concerned, though, that my setup might detract from the applicability of my results, though consultations with more experienced charmm users suggested that this was not the case.

So, I guess that's two questions/issues:

1. Are my results invalidated by my choice of water buffer (smaller than the cutoff distance)?

2. What might be causing this unexpected error?

Thanks very much!

My guess is that the unit cell edge length in the script doesn't match the actual box size of the coordinate set. If you get the step 0 energy printout, check the VDW terms, esp. IMNBvdw.

That's certainly a logical suspect, but the unit cell edge length 92.4165 is exactly consistent with that specified in the original dynamics run (confirmed from both input and output files) and, naturally, the original solvation step (confirmed from output)...

There is no step 0 energy printout, alas. However, I will insert an "energy" command directly before the dynamics command...

There is another wrinkle that I forgot to mention in my first post. If, instead of starting the dynamics with velocities at 300K, I start them at 200K and heat to 300 (see attached input file), the simulation does not crash and I do not get these errors.

Here's the dynamics command that works:

Is that the actual box size for that specific dynamics frame? The initial values don't matter as much as having the exact value for that specific coordinate set. The ENERGY command should give the exact same values for the potential terms as observed in the previous dynamics, if the box size is correct.

I see. I wasn't aware that there were different box sizes associated with specific dynamics time frames...

Indeed, when I inserted the 'energy' command before the 'dyna strt' command, the energy terms were different...

How would I go about determining the box size for that specific coordinate set?

The simulation was run using constant P, which means the box edge length changes with each integration step, and fluctuates around a length consistent with the reference pressure.

The edge lengths (and cell angles) are printed in the dynamics output, so it should be straightforward to find the size in the output log. (I've assumed that the time point of the coord set is known, and matches one of those printed in the log.)

The unit cell data is also stored with each coordinate set in the trajectory files, and can be extracted in a number of ways, as long as the CRYSTAL setup is done before reading the trajectory files(s).

Ah, thanks, Rick. I'll find or extract the edge lengths...

So, I had been chosen to run CPT with the PMASS set to zero for my [NVT] dynamics because of the recommendation at the very bottom of this page. (And those above that about having pressure info printed out.) So, I take it that's indeed a legit way of running NVT... I hadn't realized that the box size would change, though, as this would seem to change the V in the NVT, no?

It wasn't obvious from the initial post that @PMASS was zero; it should be obvious from the output log file whether or not the edge length (and volume) is constant or not.

Still, you should be able to read a coordinate set from a previous MD run and get the same energy values, if everything is the same (esp. PBC and non-bond setup).

Thanks!