Dear Justin<br><br>You are right about the cut-off. The vdw energy spike disappeared after I increased the cut-off in the rerun. But I still don't understand why? Can the cutoff error build up to several thousand kj/mol for 100AA protein.<br>
<br>Again I would like to emphasize NOT to use a xtc file in MM/PBSA type calculation. The error is just too large.<br><br>best,<br><br>dawei<br><br><div class="gmail_quote">On Thu, Aug 11, 2011 at 9:03 AM, Justin A. Lemkul <span dir="ltr"><<a href="mailto:jalemkul@vt.edu">jalemkul@vt.edu</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin: 0pt 0pt 0pt 0.8ex; border-left: 1px solid rgb(204, 204, 204); padding-left: 1ex;"><div class="im"><br>
<br>
Da-Wei Li wrote:<br>
<blockquote class="gmail_quote" style="margin: 0pt 0pt 0pt 0.8ex; border-left: 1px solid rgb(204, 204, 204); padding-left: 1ex;">
Dear Justin<br>
<br>
An implicit water simulaiton with this short cutoff is problematic but I think it is fine for rerun. I want to exactly repeat the original energies in the explicit water MD.<br>
<br>
The manu say "neighbor list searching will be performed for every frame" with rerun option. So that I don't think this is the cause.<br>
<br>
</blockquote>
<br></div>
Right, forgot about that. I still think the cutoffs are a problem, though.<br>
<br>
-Justin<br>
<br>
<blockquote class="gmail_quote" style="margin: 0pt 0pt 0pt 0.8ex; border-left: 1px solid rgb(204, 204, 204); padding-left: 1ex;">
best,<br>
<br>
dawei<div><div></div><div class="h5"><br>
<br>
On Thu, Aug 11, 2011 at 8:47 AM, Justin A. Lemkul <<a href="mailto:jalemkul@vt.edu" target="_blank">jalemkul@vt.edu</a> <mailto:<a href="mailto:jalemkul@vt.edu" target="_blank">jalemkul@vt.edu</a>>> wrote:<br>
<br>
<br>
<br>
Da-Wei Li wrote:<br>
<br>
Dear Mark and others<br>
<br>
I did more tests and thought that it might come from numerical<br>
error. The reasons are<br>
<br>
1. If I use .trr file instead of the low precision xtc file,<br>
things become better, ie, I get much less snapshots that has<br>
high energy.<br>
<br>
2. I supplied -pforce in my mdrun -rerun and found that the high<br>
vdw energy was usually caused by one pair of atoms, whose<br>
distance was just very near the clash zone, so that small error<br>
on the coordinates would cause large energy error. The force is<br>
always around 10000.<br>
<br>
3. Actually bond length and bond angle energies are also<br>
affected. I can fully reproduce these two energies if I use .trr<br>
file in my rerun but will get several tens of kj/mol error if I<br>
use .xtc file, for a protein of size of 100 AA.<br>
<br>
<br>
Now the question I still have is whether numerical error can be<br>
so large? The xtc file has a precision of 0.001 nm. Anyway, I<br>
will test more by using double precision Gromacs and define<br>
energy groups so that I can compare energy of protein directly<br>
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<br>
nojump (em.tpr is generated for energy minimization, protein is<br>
in the middle of the box)<br>
<br>
rerun .mdp file:<br>
<br>
******************************<u></u>__****************************<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<br>
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)<br>
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<br>
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<br>
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<br>
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>
******************************<u></u>__****************************<u></u>**__**************************<u></u>*<br>
<br>
Thanks all.<br>
<br>
<br>
Using cutoffs this small may be the source of your problem. Proper<br>
implicit solvent calculations require longer cutoffs than would<br>
normally be used in explicit solvent MD. Try with longer (2.0 nm)<br>
or infinite cutoffs and a fixed neighbor list (nstlist = 0) and see<br>
if that smooths out the problem. What's likely happening now is<br>
that you've got interactions moving very quickly in and out of the<br>
very short cutoff, causing spikes in energy in between neighbor list<br>
updates.<br>
<br>
-Justin<br>
<br>
-- ==============================<u></u>__==========<br>
<br>
Justin A. Lemkul<br>
Ph.D. Candidate<br>
ICTAS Doctoral Scholar<br>
MILES-IGERT Trainee<br>
Department of Biochemistry<br>
Virginia Tech<br>
Blacksburg, VA<br></div></div>
jalemkul[at]<a href="http://vt.edu" target="_blank">vt.edu</a> <<a href="http://vt.edu" target="_blank">http://vt.edu</a>> | <a href="tel:%28540%29%20231-9080" value="+15402319080" target="_blank">(540) 231-9080</a><br>
<tel:%28540%29%20231-9080><div class="im"><br>
<a href="http://www.bevanlab.biochem." target="_blank">http://www.bevanlab.biochem.</a>__<a href="http://vt.edu/Pages/Personal/justin" target="_blank"><u></u>vt.edu/Pages/Personal/justin</a><br>
<<a href="http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin" target="_blank">http://www.bevanlab.biochem.<u></u>vt.edu/Pages/Personal/justin</a>><br>
<br>
==============================<u></u>__==========<br>
<br>
-- gmx-users mailing list <a href="mailto:gmx-users@gromacs.org" target="_blank">gmx-users@gromacs.org</a><br></div>
<mailto:<a href="mailto:gmx-users@gromacs.org" target="_blank">gmx-users@gromacs.org</a>><div class="im"><br>
<a href="http://lists.gromacs.org/__mailman/listinfo/gmx-users" target="_blank">http://lists.gromacs.org/__<u></u>mailman/listinfo/gmx-users</a><br>
<<a href="http://lists.gromacs.org/mailman/listinfo/gmx-users" target="_blank">http://lists.gromacs.org/<u></u>mailman/listinfo/gmx-users</a>><br>
Please search the archive at<br>
<a href="http://www.gromacs.org/__Support/Mailing_Lists/Search" target="_blank">http://www.gromacs.org/__<u></u>Support/Mailing_Lists/Search</a><br>
<<a href="http://www.gromacs.org/Support/Mailing_Lists/Search" target="_blank">http://www.gromacs.org/<u></u>Support/Mailing_Lists/Search</a>> before posting!<br>
Please don't post (un)subscribe requests to the list. Use the www<br>
interface or send it to <a href="mailto:gmx-users-request@gromacs.org" target="_blank">gmx-users-request@gromacs.org</a><br></div>
<mailto:<a href="mailto:gmx-users-request@gromacs.org" target="_blank">gmx-users-request@<u></u>gromacs.org</a>>.<div class="im"><br>
Can't post? Read <a href="http://www.gromacs.org/__Support/Mailing_Lists" target="_blank">http://www.gromacs.org/__<u></u>Support/Mailing_Lists</a><br></div><div class="im">
<<a href="http://www.gromacs.org/Support/Mailing_Lists" target="_blank">http://www.gromacs.org/<u></u>Support/Mailing_Lists</a>><br>
<br>
<br>
</div></blockquote><div><div></div><div class="h5">
<br>
-- <br>
==============================<u></u>==========<br>
<br>
Justin A. Lemkul<br>
Ph.D. Candidate<br>
ICTAS Doctoral Scholar<br>
MILES-IGERT Trainee<br>
Department of Biochemistry<br>
Virginia Tech<br>
Blacksburg, VA<br>
jalemkul[at]<a href="http://vt.edu" target="_blank">vt.edu</a> | <a href="tel:%28540%29%20231-9080" value="+15402319080" target="_blank">(540) 231-9080</a><br>
<a href="http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin" target="_blank">http://www.bevanlab.biochem.<u></u>vt.edu/Pages/Personal/justin</a><br>
<br>
==============================<u></u>==========<br>
-- <br>
gmx-users mailing list <a href="mailto:gmx-users@gromacs.org" target="_blank">gmx-users@gromacs.org</a><br>
<a href="http://lists.gromacs.org/mailman/listinfo/gmx-users" target="_blank">http://lists.gromacs.org/<u></u>mailman/listinfo/gmx-users</a><br>
Please search the archive at <a href="http://www.gromacs.org/Support/Mailing_Lists/Search" target="_blank">http://www.gromacs.org/<u></u>Support/Mailing_Lists/Search</a> before posting!<br>
Please don't post (un)subscribe requests to the list. Use the www interface or send it to <a href="mailto:gmx-users-request@gromacs.org" target="_blank">gmx-users-request@gromacs.org</a>.<br>
Can't post? Read <a href="http://www.gromacs.org/Support/Mailing_Lists" target="_blank">http://www.gromacs.org/<u></u>Support/Mailing_Lists</a><br>
</div></div></blockquote></div><br>