Previous Thread
Next Thread
Print Thread
#27784 06/28/11 09:32 PM
Joined: Feb 2010
Posts: 29
Ascaris Offline OP
Forum Member
OP Offline
Forum Member
Joined: Feb 2010
Posts: 29
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").
Quote:
dyna leap cpt strt -
timestep 0.001 nstep 1000000 nprint 100 -
iunrea 51 iunwri 41 iuncrd 31 nsavcc 200 -
hoover reft 300.0 tmass @tmass -
pcons pint pref 1.0 pmass @pmass pgamma 0.0 -
ihtfrq 0 ieqfrq 0 ntrfrq 500 -
tstruc 300.0 firstt 300.0 -
iasors 1 iasvel 1 ichecw 0 -
iseed 3789047890 echeck 100.



As I had not had trouble doing this with previous setups, I was surprised by this error.
Quote:
** ERROR IN SHAKEA ** DEVIATION IN SHAKE TOO LARGE
NITER= 0 LL=57706 I=61143 J=61145
TOLER= 0.9162318400 DIFF=****************** RRPR=******************
X(I) -33.3024932200 -4.5760079566 -46.7711242069
XREF(I) -35.0977593259 -2.0397219126 -45.2191160940
X(J) -34.1208428564 -2.1318231170 -45.1587726655
XREF(J) -34.1456931937 -2.1288405805 -45.1759908908
*** LEVEL 1 WARNING *** BOMLEV IS 0

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

Attached Images
charmm.inp.txt (3.08 KB, 359 downloads)

Intestinal nematode
Bronx, NY
Joined: Sep 2003
Posts: 8,606
Likes: 24
rmv Online Content
Forum Member
Online Content
Forum Member
Joined: Sep 2003
Posts: 8,606
Likes: 24
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.


Rick Venable
computational chemist

Joined: Feb 2010
Posts: 29
Ascaris Offline OP
Forum Member
OP Offline
Forum Member
Joined: Feb 2010
Posts: 29
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:
Quote:
dyna leap verlet start -
timestep 0.001 nstep 50000 nprint 100 -
iunwri 41 iuncrd 31 nsavc 500 -
firstt 200.0 finalt 300.0 tbath 300.0 -
ihtfrq 2500 teminc 5 ieqfrq 0 -
iasors 1 iasvel 1 iscvel 0 ichecw 0 -
ntrfrq 500 -
iseed 40867920 -
echeck 100.0

Attached Images
charmm_heat.inp.txt (2.85 KB, 326 downloads)
Last edited by Ascaris; 06/29/11 06:33 PM.

Intestinal nematode
Bronx, NY
Joined: Sep 2003
Posts: 8,606
Likes: 24
rmv Online Content
Forum Member
Online Content
Forum Member
Joined: Sep 2003
Posts: 8,606
Likes: 24
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.


Rick Venable
computational chemist

Joined: Feb 2010
Posts: 29
Ascaris Offline OP
Forum Member
OP Offline
Forum Member
Joined: Feb 2010
Posts: 29
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?


Intestinal nematode
Bronx, NY
Joined: Sep 2003
Posts: 8,606
Likes: 24
rmv Online Content
Forum Member
Online Content
Forum Member
Joined: Sep 2003
Posts: 8,606
Likes: 24
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).


Rick Venable
computational chemist

Joined: Feb 2010
Posts: 29
Ascaris Offline OP
Forum Member
OP Offline
Forum Member
Joined: Feb 2010
Posts: 29
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?


Intestinal nematode
Bronx, NY
Joined: Sep 2003
Posts: 8,606
Likes: 24
rmv Online Content
Forum Member
Online Content
Forum Member
Joined: Sep 2003
Posts: 8,606
Likes: 24
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).


Rick Venable
computational chemist

Joined: Feb 2010
Posts: 29
Ascaris Offline OP
Forum Member
OP Offline
Forum Member
Joined: Feb 2010
Posts: 29
Thanks!

Last edited by Ascaris; 06/29/11 07:53 PM. Reason: Error

Intestinal nematode
Bronx, NY

Moderated by  BRBrooks, lennart, rmv 

Link Copied to Clipboard
Powered by UBB.threads™ PHP Forum Software 7.7.5
(Release build 20201027)
Responsive Width:

PHP: 7.3.31-1~deb10u1 Page Time: 0.010s Queries: 34 (0.006s) Memory: 0.7807 MB (Peak: 0.8776 MB) Data Comp: Off Server Time: 2022-10-01 01:23:23 UTC
Valid HTML 5 and Valid CSS