[AMBER-Developers] inconsistency in sander manual and code for restart velocities

From: Carlos Simmerling <carlos.simmerling.gmail.com>
Date: Tue, 6 Jan 2009 16:58:08 -0500

hi all- I've been tracking down a bug in sander for amber9/10 and know
that one exists, but it's hard to say if it's in the manual or the
I just want some quick feedback since there could be several ways to "fix" it.
it's related to all of the overly complicated combinations of irest,
ntx, tempi and init.

my understanding, along with page 26 of the amber10 manual, is that
for a PBC system, if one wants to restart with box info from inpcrd
(say from a constant pressure run), then one should use ntx=5.

from the manual:
ntx = 5 X and V are read formatted; box information will be read if
ntb>0. The velocity information
will only be used if irest=1.

it's the last part that's wrong. if I set ntx=5 to read the box, but
irest=0 so that my velocities are assigned by sander, setvel() is
still not called in sander.f (it is only called if ntx <=3). this
means the run is the same no matter what the random seed is. also,
init is passed to setvel but never used.

irest is indeed used, in mdread2(), to set the value of init:

   init = 3
   if (irest > 0) init = 4

the problem is that init isn't used for anything related to setting
velocities, but it is used in runmd to do the setup if it's not a

I got burned on this by starting a large # of runs from the same
structure with different ig values, thinking I would get different
runs. all were identical (with ntt=1, with ntt=3 they would not be
because the ig value would change the thermostat, but I still don't
think that's what we want).

SO, I think the problem is that we use ntx for deciding on calling
setvel, but we really should use ntx for reading them, and irest for
calling setvel (or init- this is really confusing since they map to
different values). alternately, the setvel call could be moved into
the init=3 part of runmd, though I haven't thought much about whether
there is a reason it isn't there already.

either that or we change the manual to say that any time we read
velocities we will use them (but then can we still read the box? it
appears that getcor() doesn't ever read the box for formatted files,
but does so for unformatted! so it's unclear to me right now if we
reliably read the box info for ntx=1 PBC restarts).

I could fix it but want to make sure I know our "intent" before making
changes or a bugfix. my suggestion is to fix the code, but it might
break old inputs. still, given the issues lately with needing to set
different ig values, we should allow users to actually do it the way
that the manual has been describing!


AMBER-Developers mailing list
Received on Tue Jan 06 2009 - 15:33:47 PST
Custom Search