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").
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.
** 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!