From: Rajan Vatassery (rajan_at_umn.edu)
Date: Fri Dec 14 2012 - 13:38:02 CST
Hi Eddie,
The NpT simulation provides a log file which I extract the volume
information from. I created a python code, below, but I think VMD has
tools for this extraction as well in the extensions tab. I use xmgrace
to delete the first ~1000 points from the resulting file to make sure
that I am in a region where the volume is changing only due to random
fluctuations. Then I take the avg of the remaining points. I suppose
there would be some way to get NAMD to do all this automatically, but I
think most people would prefer to see the graph themselves and decide
where to start and stop dropping points.
As for the temperature, you're right about the fast heating. In my case
the system does not have complicated coordinated motions such as a
protein or some very large molecule would. My molecules are all under 60
atoms, and so I believe rapid heating will not contribute any
non-physical interactions. The non-physical interactions that I have
detected and fixed in my system are all the result of changing the box
size from the NpT to NVT simulations. In other, more complicated
systems, I would strongly suggest people to be careful with heating as
you have mentioned.
rajan
-----------------------------------------------------------------------------------
#!/usr/bin/python
import os,sys,math
file1 = open(sys.argv[1],'r')
strin = file1.readlines()
TS = []
BOND = []
ANGLE = []
DIHED = []
IMPRP = []
ELECT = []
VDW = []
BOUNDARY = []
MISC = []
KINETIC = []
TOTAL = []
TEMP = []
POTENTIAL = []
TOTAL3 = []
TEMPAVG = []
PRESSURE = []
GPRESSURE = []
VOLUME = []
PRESSAVG = []
GPRESSAVG = []
for X in range(len(strin)):
if strin[X][0:6] == 'ENERGY':
str1 = strin[X].split()
TS.append(str1[1]),BOND.append(str1[2]),ANGLE.append(str1[3]),DIHED.append(str1[4]),IMPRP.append(str1[5]),ELECT.append(str1[6]),VDW.append(str1[7]),BOUNDARY.append(str1[8]),MISC.append(str1[9]),KINETIC.append(str1[10]),TOTAL.append(str1[11]),TEMP.append(str1[12]),POTENTIAL.append(str1[13]),TOTAL3.append(str1[14]),TEMPAVG.append(str1[15]),PRESSURE.append(str1[16]),GPRESSURE.append(str1[17]),VOLUME.append(str1[18]),PRESSAVG.append(str1[19]),GPRESSAVG.append(str1[20])
print str1
print
TS[len(TS)-1],BOND[len(BOND)-1],ANGLE[len(ANGLE)-1],DIHED[len(DIHED)-1],IMPRP[len(IMPRP)-1],ELECT[len(ELECT)-1],VDW[len(VDW)-1],BOUNDARY[len(BOUNDARY)-1],MISC[len(MISC)-1],KINETIC[len(KINETIC)-1],TOTAL[len(TOTAL)-1],TEMP[len(TEMP)-1],POTENTIAL[len(POTENTIAL)-1],TOTAL3[len(TOTAL3)-1],TEMPAVG[len(TEMPAVG)-1],PRESSURE[len(PRESSURE)-1],GPRESSURE[len(GPRESSURE)-1],VOLUME[len(VOLUME)-1],PRESSAVG[len(PRESSAVG)-1],GPRESSAVG[len(GPRESSAVG)-1]
newfile = sys.argv[1].split('.')
if sys.argv[3] == 'w':
file2 = open(newfile[0]+'.'+sys.argv[2],'w')
elif sys.argv[3] == 'a':
file2 = open(newfile[0]+'.'+sys.argv[2],'a')
else:
print 'Error with last command line argument'
arr1 = eval(sys.argv[2].upper())
for X in range(len(arr1)):
#print arr1[X]
file2.write('%9d %12.4f\n' % (int(TS[X]),float(arr1[X])))
file1.close()
file2.close()
-----------------------------------------------------------------------------------
USAGE:
FILENAME.py LOGFILE.log volume w
"volume" can be replaced with any of the other entries on the ENERGY
line, such as "angle", "dihed", "potential", etc. The "w" is an artifact
and just needs to be appended to the end of each command.
On Fri, 2012-12-14 at 10:33 -0600, Dr. Eddie wrote:
> Rajan,
> Thank you so much for this! It is very helpful!
>
>
> I do have a few questions. The cell size in SU was put in manually
> from the xsc file? Or is there some automatic way of making namd read
> that cell size?
>
>
> You set the temperature to 298 but later on you minimize the structure
> before running. Thus your temperature will be set to 0. The langevin
> temperature gauge will bring it back up to 298 though. Is this way of
> doing heating smooth in your experience compared with a loop over
> temperatures?
> Thanks!
> Eddie
>
>
>
>
>
>
>
>
> On Mon, Dec 10, 2012 at 3:55 PM, Rajan Vatassery <rajan_at_umn.edu>
> wrote:
> Hi Eddie,
> Here's the gzipped tarball. I have included the
> following:
>
> I8i conf file: Initial NpT simulation conf file that starts
> from the
> initial inputs of I8i14_6426.psf/pdb.
>
> this conf file dictates that the output should be placed in
> the
> folder /lustre/vatasser/ClassicalMech/, and the files
> I8i14_6426.coor
> and I8i14_6426.vel will be your input for the NVT simulation.
>
> SU conf file: NVT simulation that sets up the NVE simulation
> after a ns.
> Note the difference between the cellBasisVectors of the NpT
> and the NVT
> simulation.
>
> NVE conf file: NVE simulation that runs for 20 ns based on
> inputs from
> SU run. Note the cellBasisVectors are commented out.
>
> The comments in the conf files are generally irrelevant so
> probably best
> to ignore them.
>
> Good luck and let me know if there's anything else you need.
>
> rajan
>
> On Mon, 2012-12-10 at 15:00 -0600, Dr. Eddie wrote:
> > Yes, if you don't mind that would be great. I have some idea
> of how to
> > do what you describe but seeing it would really help!
> >
> >
> > On Mon, Dec 10, 2012 at 1:33 PM, Rajan Vatassery
> <rajan_at_umn.edu>
> > wrote:
> > Eddie,
> > I usually start with a PDB/PSF that comes
> from VMD,
> > and using the
> > default setting:
> >
> > binaryoutput yes
> >
> > I receive a .coor, .vel, and a .xsc file from the
> first (NpT)
> > simulation. This simulation minimizes and
> equilibrates the
> > system at a
> > variable volume.
> > The next simulation, NVT, receives the
> following input
> > from the NpT
> > simluation:
> >
> > binCoordinates NpT.coor
> > binVelocities NpT.vel
> >
> > I find the average volume that was produced from the
> NpT
> > simulation
> > after the volume stops changing drastically, and
> this number
> > is set as
> > the fixed volume with the lines:
> >
> > cellBasisVector1 xvector 0 0
> > cellBasisVector2 0 yvector 0
> > cellBasisVector3 0 0 zvector
> > cellOrigin 0 0 0
> >
> > where xvector*yvector*zvector is your total volume.
> NVT runs
> > and outputs
> > the same as the NpT simulation, and this time I use
> all three
> > (.vel, .coor, and .xsc) files as inputs to the NVE
> simulation.
> >
> > binCoordinates NVT.coor
> > binVelocities NVT.vel
> > extendedSystem NVT.xsc
> >
> > Because in the NVE simulation I am using
> the .xsc file
> > to determine the
> > PBC boundaries, I comment out the cellBasisVector
> lines. Note
> > that going
> > from NpT to NVT sometimes produces erroneously high
> energies.
> > Basically
> > the source of this problem is that you are taking a
> volume of
> > your
> > simulation (the average over the NpT) that doesn't
> necessarily
> > correspond to the .vel and .coor files (which are
> just the
> > snapshots of
> > the last point in the NpT). If the volume of your
> last point
> > is much
> > different than the average, some strange things can
> happen as
> > a result
> > of the compression/expansion of your system. The way
> to detect
> > if your
> > system has suffered from this is to watch the
> energies of your
> > system
> > from NpT to NVT to NVE. If the system energy has
> changed
> > drastically
> > from NpT to NVE, and that change isn't just a result
> of a
> > small T
> > difference from removing the thermostat, then you'll
> have to
> > find the
> > coordinates and velocities from the .dcd file of a
> > configuration where
> > the volume is closer to your average volume. This is
> a
> > complicated
> > process, and even more complicated to describe
> without
> > figures/equations, but the solution is simple.
> > I can also send over my conf files if you
> would like.
> >
> > rajan
> >
> >
> >
> > On Mon, 2012-12-10 at 12:20 -0600, Dr. Eddie wrote:
> > > Rajan,
> > > May I ask how you did the start and restart
> exactly? Did you
> > build a
> > > new pdb file in vmd and start with that after the
> Langevin
> > dynamics
> > > were done?
> > > Thanks,
> > > Eddie
> > >
> > >
> > > On Mon, Dec 10, 2012 at 11:00 AM, Rajan Vatassery
> > <rajan_at_umn.edu>
> > > wrote:
> > > Hi Eddie,
> > > Unfortunately, as Norman Geist
> suggested,
> > the best way
> > > to accomplish
> > > what I was trying to do was to set up
> multiple conf
> > files. I
> > > looked into
> > > the "checkpoint" keyword as well, but I
> couldn't
> > find enough
> > > details to
> > > figure out how to use it properly. It
> sounds like
> > you're
> > > encountering a
> > > similar situation.
> > > My obsessive-compulsive side
> really wanted
> > to keep my
> > > directories
> > > simple by using only one conf file and one
> set of
> > output
> > > files, but
> > > separating things was just easier. If you
> find out
> > how to use
> > > the
> > > checkpoint/revert keywords, I'm sure there
> are a lot
> > of people
> > > out there
> > > that could use a few pointers.
> > >
> > > Good luck,
> > >
> > > Rajan
> > >
> > > On Mon, 2012-12-10 at 09:08 -0600, Dr.
> Eddie wrote:
> > > > Thanks Rajan. So how did you "solve"
> turning the
> > langevin
> > > dynamics
> > > > off? Restart the simulation somehow?
> Were you able
> > to
> > > restart the
> > > > simulation in 1 config file?
> > > >
> > > >
> > > > I know there are commands 'checkpoint'
> and
> > 'revert', from
> > > here, but
> > > > the documentation on these is scarce at
> best! It
> > sounds like
> > > I should
> > > > be able to make a checkpoint and then
> revert back
> > changing
> > > some system
> > > > parameters, but I have no idea how!
> > > > Thanks!
> > > > Eddie
> > > >
> > > >
> > > > On Sun, Dec 9, 2012 at 5:02 PM, Rajan
> Vatassery
> > > <rajan_at_umn.edu> wrote:
> > > > Eddie,
> > > > I ran into a similar
> problem as
> > your last
> > > error
> > > > message a while back:
> > > >
> > > >
> > >
> >
> http://www.ks.uiuc.edu/Research/namd/mailing_list/namd-l/14609.html
> > > >
> > > > It's not possible to change all
> of the
> > conf file
> > > > keywords/variables
> > > > after the first "run" command in
> the conf
> > file.
> > > Here's the
> > > > page to read:
> > > >
> > > >
> > >
> >
> http://www.ks.uiuc.edu/Research/namd/2.8/ug/node10.html
> > > >
> > > > I haven't seen your first error
> before so
> > I can't
> > > comment on
> > > > that.
> > > >
> > > >
> > > > Good luck,
> > > >
> > > >
> > > > rajan
> > > >
> > > > On Sun, 2012-12-09 at 10:16
> -0600, Dr.
> > Eddie wrote:
> > > > > Hi all,
> > > > > I've been trying to setup a
> single
> > config file
> > > that would do
> > > > the
> > > > > minimization, then slow
> heating with
> > pressure
> > > equilibration
> > > > and then
> > > > > turn off the Langevin dynamics
> once the
> > > structure/solvent
> > > > has reached
> > > > > the desired
> temperature/pressure.
> > > > >
> > > > > I have tried many things and
> would like
> > some
> > > advice to make
> > > > this work
> > > > > (if it is possible) and ensure
> I'm not
> > making some
> > > gross
> > > > error.
> > > > > The relevant part of the
> config file is
> > (my number
> > > of
> > > > timesteps is
> > > > > small so I can see the crashes
> quickly)
> > > > >
> > > > > langevin on
> > > > > langevinDamping 1.
> > > > > langevinTemp 100
> # <----
> > 1st problem
> > > is if this
> > > > (and all
> > > > > the temperatures below) is
> lower than
> > 100
> > > > >
> > > #REINITIALIZING
> > > > VELOCITIES AT
> > > > > STEP 16 TO 25 KELVIN.
> > > > >
> > #TCL:
> > > Running for 48
> > > > steps
> > > > >
> > #FATAL
> > > ERROR:
> > > > Periodic cell has
> > > > > become too small for original
> patch
> > grid!
> > > > > langevinHydrogen no
> > > > > useGroupPressure yes
> > > > > useFlexibleCell no
> > > > > useConstantArea no
> > > > > langevinPiston on
> > > > > langevinPistonTarget 1.01325
> > > > > langevinPistonPeriod 100.
> > > > > langevinPistonDecay 50.
> > > > > langevinPistonTemp 100
> > > > > temperature 100
> > > > > for { set MIN 0 } { $MIN < 3 }
> { incr
> > MIN 1 } {
> > > > > minimize 16
> > > > > reinitvels 100
> > > > > run 48
> > > > > }
> > > > >
> > > > > # I need to start my heating
> at 100 k
> > (set TEMP
> > > 100 below)
> > > > otherwise I
> > > > > also get the error:
> > > > > #TCL: Setting parameter
> > langevinPistonTemp to 25
> > > > > #TCL: Setting parameter
> langevinTemp to
> > 25
> > > > > #REINITIALIZING VELOCITIES AT
> STEP 192
> > TO 25
> > > KELVIN.
> > > > > #TCL: Running for 48 steps
> > > > > #FATAL ERROR: Periodic cell
> has become
> > too small
> > > for
> > > > original patch
> > > > > grid!
> > > > > for { set TEMP 100} { $TEMP <
> 285 }
> > { incr TEMP
> > > 25 } {
> > > > > langevinPistonTemp $TEMP
> > > > > langevinTemp $TEMP
> > > > > reinitvels $TEMP
> > > > > run 48
> > > > > }
> > > > > temperature 310
> #<------- If
> > I use 100
> > > K for
> > > > minimization
> > > > > etc This command fails with:
> > > > >
> > #TCL: Setting
> > > parameter
> > > > > temperature to 310
> > > > >
> > #FATAL ERROR:
> > > Setting
> > > > parameter
> > > > > temperature from script
> failed!
> > > > >
> # DO
> > I need
> > > it?
> > > > >
> > > > > langevinPistonTemp 310
> > > > > langevinTemp 310
> > > > > reinitvels 310
> > > > > run 50000
> > > > >
> > > > > langevinPiston off
> #<---- doing
> > all of the
> > > above
> > > > (which seems
> > > > > dubious) this fails with:
> > > > >
> #TCL:
> > Setting
> > > parameter
> > > > > langevinPiston to off
> > > > >
> #TCL:
> > Setting
> > > parameter
> > > > langevin to
> > > > > off
> > > > >
> #FATAL
> > ERROR:
> > > Setting
> > > > parameter
> > > > > langevin from script failed!
> > > > > langevin off
> > > > > run 20000000
> > > > >
> > > > >
> > > > > I'd like to end up with a
> minimized,
> > slowly heated
> > > system at
> > > > the right
> > > > > temperature/pressure. Then I'd
> like to
> > run that
> > > system
> > > > without the
> > > > > Langevin forces on to see how
> much for
> > the
> > > protein's
> > > > sidechains alone
> > > > > move (rmsd, distances to
> other CA's
> > etc) without
> > > any random
> > > > forces
> > > > > increasing the distance.
> > > > >
> > > > > Thank you very much for any
> help!
> > > > > Eddie
> > > > >
> > > > >
> > > > > --
> > > > > Eddie
> > > > >
> > > >
> > > >
> > > >
> > > >
> > > >
> > > >
> > > > --
> > > > Eddie
> > > >
> > > >
> > >
> > >
> > >
> > >
> > >
> > >
> > > --
> > > Eddie
> > >
> > >
> >
> >
> >
> >
> >
> >
> > --
> > Eddie
> >
> >
>
>
>
>
>
> --
> Eddie
>
>
This archive was generated by hypermail 2.1.6 : Mon Dec 31 2012 - 23:22:22 CST