<html>
<head>
<meta content="text/html; charset=ISO-8859-1"
http-equiv="Content-Type">
</head>
<body bgcolor="#FFFFFF" text="#000000">
On 11/08/2011 10:22 PM, Da-Wei Li wrote:
<blockquote
cite="mid:CAFs1QTc9JUt2bdraYDWXCvp=2yDjsHda8vdtpe7oQwA1N_-BHw@mail.gmail.com"
type="cite">Dear Mark and others<br>
<br>
I did more tests and thought that it might come from numerical
error. The reasons are<br>
<br>
1. If I use .trr file instead of the low precision xtc file,
things become better, ie, I get much less snapshots that has high
energy.<br>
</blockquote>
<br>
That doesn't make sense if the underlying frames are the same in
your .xtc and .trr files. However, if you're talking about different
numbers of frames, then you probably have transient periodicity
artefacts.<br>
<br>
<blockquote
cite="mid:CAFs1QTc9JUt2bdraYDWXCvp=2yDjsHda8vdtpe7oQwA1N_-BHw@mail.gmail.com"
type="cite">
<br>
2. I supplied -pforce in my mdrun -rerun and found that the high
vdw energy was usually caused by one pair of atoms, whose distance
was just very near the clash zone, so that small error on the
coordinates would cause large energy error. The force is always
around 10000.<br>
</blockquote>
<br>
Sounds like you have a non-equilibrium trajectory (from your EM?),
in which case small deviations can have large effects.<br>
<br>
<blockquote
cite="mid:CAFs1QTc9JUt2bdraYDWXCvp=2yDjsHda8vdtpe7oQwA1N_-BHw@mail.gmail.com"
type="cite">
<br>
3. Actually bond length and bond angle energies are also affected.
I can fully reproduce these two energies if I use .trr file in my
rerun but will get several tens of kj/mol error if I use .xtc
file, for a protein of size of 100 AA.<br>
</blockquote>
<br>
Several tens of kJ/mol error in the non-bonded energy looks OK
compared to the magnitude of the non-bonded energy. Anyway, if your
purpose is to compare energies computed from the same
configurations, then you should use the same configurations for
computing the energies, not approximations to those configurations.<br>
<br>
<blockquote
cite="mid:CAFs1QTc9JUt2bdraYDWXCvp=2yDjsHda8vdtpe7oQwA1N_-BHw@mail.gmail.com"
type="cite">
<br>
<br>
Now the question I still have is whether numerical error can be so
large? The xtc file has a precision of 0.001 nm. Anyway, I will
test more by using double precision Gromacs and define energy
groups so that I can compare energy of protein directly between
original MD and rerun.<br>
<br>
<br>
<br>
To Mark only<br>
<br>
Thanks. Here it is my script for <br>
<br>
rerun: mdrun -v -s pbsa.tpr -rerun coor.xtc -e rerun<br>
superpose: trjconv -s em.tpr -f coor.xtc -o nojump.xtc -pbc
nojump (em.tpr is generated for energy minimization, protein is
in the middle of the box)<br>
</blockquote>
<br>
Here you are generating a nojump trajectory file, but you are not
doing your rerun on it. That could explain some symptoms - you might
not have whole molecules.<br>
<br>
Mark<br>
<br>
<blockquote
cite="mid:CAFs1QTc9JUt2bdraYDWXCvp=2yDjsHda8vdtpe7oQwA1N_-BHw@mail.gmail.com"
type="cite">
<br>
rerun .mdp file:<br>
<br>
**********************************************************<br>
<br>
; Run parameters<br>
integrator = md ; leap-frog integrator<br>
nsteps = 50000000 ; 100 ns<br>
dt = 0.002 ; 2 fs<br>
; Output control<br>
nstxout = 500000 ; save coordinates every 1000 ps<br>
nstvout = 500000 ; save velocities every 1000 ps<br>
nstxtcout = 5000 ; xtc compressed trajectory output
every 1 ps <br>
nstenergy = 5000 ; save energies every 1 ps<br>
nstlog = 5000 ; update log file every 1 ps<br>
xtc_grps = Protein ; save protein part only<br>
; Bond parameters<br>
continuation = yes ; Restarting after NPT <br>
constraint_algorithm = lincs ; holonomic constraints <br>
constraints = hbonds ; all bonds (even heavy atom-H bonds)
constrained<br>
lincs_iter = 1 ; accuracy of LINCS<br>
lincs_order = 4 ; also related to accuracy<br>
; Neighborsearching<br>
ns_type = grid ; search neighboring grid cels<br>
nstlist = 10 ; 20 fs<br>
rlist = 0.8 ; short-range neighborlist cutoff (in
nm)<br>
rcoulomb = 0.8 ; short-range electrostatic cutoff (in
nm)<br>
rvdw = 1.0 ; short-range van der Waals cutoff (in
nm)<br>
; Electrostatics<br>
coulombtype = cut-off ; Particle Mesh Ewald for long-range
electrostatics<br>
pme_order = 4 ; cubic interpolation<br>
fourierspacing = 0.12 ; grid spacing for FFT<br>
; Temperature coupling is on<br>
tcoupl = no ; modified Berendsen thermostat<br>
tc-grps = System ; two coupling groups - more accurate<br>
tau_t = 0.1 ; time constant, in ps<br>
ref_t = 300 ; reference temperature, one for each
group, in K<br>
; Pressure coupling is on<br>
pcoupl = no ; Pressure coupling on in NPT<br>
pcoupltype = isotropic ; uniform scaling of box vectors<br>
tau_p = 2.0 ; time constant, in ps<br>
ref_p = 1.0 ; reference pressure, in bar<br>
compressibility = 4.5e-5 ; isothermal compressibility of water,
bar^-1<br>
; Periodic boundary conditions<br>
pbc = no ; 3-D PBC<br>
; Dispersion correction<br>
;DispCorr = EnerPres ; account for cut-off vdW scheme<br>
DispCorr = no<br>
; Velocity generation<br>
gen_vel = no ; Velocity generation is off <br>
<br>
<br>
<br>
; IMPLICIT SOLVENT ALGORITHM<br>
implicit_solvent = GBSA <br>
gb_algorithm = OBC <br>
nstgbradii = 1<br>
rgbradii = 0.8 <br>
gb_epsilon_solvent = 80<br>
gb_saltconc = 0<br>
gb_obc_alpha = 1<br>
gb_obc_beta = 0.8<br>
gb_obc_gamma = 4.85<br>
gb_dielectric_offset = 0.009<br>
sa_algorithm = Ace-approximation<br>
sa_surface_tension = 2.25936<br>
<br>
***************************************************************************************<br>
<br>
Thanks all.<br>
<br>
dawei<br>
<br>
<br>
<br>
<br>
<br>
<br>
<br>
<br>
<fieldset class="mimeAttachmentHeader"></fieldset>
<br>
</blockquote>
<br>
</body>
</html>