Da-Wei Li wrote:
Dear Justin
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.
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.
Right, forgot about that. I still think the cutoffs are a problem, though.
-Justin
best,
dawei
On Thu, Aug 11, 2011 at 8:47 AM, Justin A. Lemkul <jalem...@vt.edu
<mailto:jalem...@vt.edu>> wrote:
Da-Wei Li wrote:
Dear Mark and others
I did more tests and thought that it might come from numerical
error. The reasons are
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.
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.
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.
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.
To Mark only
Thanks. Here it is my script for
rerun: mdrun -v -s pbsa.tpr -rerun coor.xtc -e rerun
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)
rerun .mdp file:
******************************__****************************
; Run parameters
integrator = md ; leap-frog integrator
nsteps = 50000000 ; 100 ns
dt = 0.002 ; 2 fs
; Output control
nstxout = 500000 ; save coordinates every 1000 ps
nstvout = 500000 ; save velocities every 1000 ps
nstxtcout = 5000 ; xtc compressed trajectory output
every 1 ps
nstenergy = 5000 ; save energies every 1 ps
nstlog = 5000 ; update log file every 1 ps
xtc_grps = Protein ; save protein part only
; Bond parameters
continuation = yes ; Restarting after NPT
constraint_algorithm = lincs ; holonomic constraints
constraints = hbonds ; all bonds (even heavy atom-H bonds)
constrained
lincs_iter = 1 ; accuracy of LINCS
lincs_order = 4 ; also related to accuracy
; Neighborsearching
ns_type = grid ; search neighboring grid cels
nstlist = 10 ; 20 fs
rlist = 0.8 ; short-range neighborlist cutoff (in nm)
rcoulomb = 0.8 ; short-range electrostatic cutoff (in nm)
rvdw = 1.0 ; short-range van der Waals cutoff (in nm)
; Electrostatics
coulombtype = cut-off ; Particle Mesh Ewald for long-range
electrostatics
pme_order = 4 ; cubic interpolation
fourierspacing = 0.12 ; grid spacing for FFT
; Temperature coupling is on
tcoupl = no ; modified Berendsen thermostat
tc-grps = System ; two coupling groups - more accurate
tau_t = 0.1 ; time constant, in ps
ref_t = 300 ; reference temperature, one for each
group, in K
; Pressure coupling is on
pcoupl = no ; Pressure coupling on in NPT
pcoupltype = isotropic ; uniform scaling of box vectors
tau_p = 2.0 ; time constant, in ps
ref_p = 1.0 ; reference pressure, in bar
compressibility = 4.5e-5 ; isothermal compressibility of
water, bar^-1
; Periodic boundary conditions
pbc = no ; 3-D PBC
; Dispersion correction
;DispCorr = EnerPres ; account for cut-off vdW scheme
DispCorr = no
; Velocity generation
gen_vel = no ; Velocity generation is off
; IMPLICIT SOLVENT ALGORITHM
implicit_solvent = GBSA
gb_algorithm = OBC
nstgbradii = 1
rgbradii = 0.8
gb_epsilon_solvent = 80
gb_saltconc = 0
gb_obc_alpha = 1
gb_obc_beta = 0.8
gb_obc_gamma = 4.85
gb_dielectric_offset = 0.009
sa_algorithm = Ace-approximation
sa_surface_tension = 2.25936
******************************__******************************__***************************
Thanks all.
Using cutoffs this small may be the source of your problem. Proper
implicit solvent calculations require longer cutoffs than would
normally be used in explicit solvent MD. Try with longer (2.0 nm)
or infinite cutoffs and a fixed neighbor list (nstlist = 0) and see
if that smooths out the problem. What's likely happening now is
that you've got interactions moving very quickly in and out of the
very short cutoff, causing spikes in energy in between neighbor list
updates.
-Justin
--
==============================__==========
Justin A. Lemkul
Ph.D. Candidate
ICTAS Doctoral Scholar
MILES-IGERT Trainee
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu <http://vt.edu> | (540) 231-9080
<tel:%28540%29%20231-9080>
http://www.bevanlab.biochem.__vt.edu/Pages/Personal/justin
<http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin>
==============================__==========
--
gmx-users mailing list gmx-users@gromacs.org
<mailto:gmx-users@gromacs.org>
http://lists.gromacs.org/__mailman/listinfo/gmx-users
<http://lists.gromacs.org/mailman/listinfo/gmx-users>
Please search the archive at
http://www.gromacs.org/__Support/Mailing_Lists/Search
<http://www.gromacs.org/Support/Mailing_Lists/Search> before posting!
Please don't post (un)subscribe requests to the list. Use the www
interface or send it to gmx-users-requ...@gromacs.org
<mailto:gmx-users-requ...@gromacs.org>.
Can't post? Read http://www.gromacs.org/__Support/Mailing_Lists
<http://www.gromacs.org/Support/Mailing_Lists>
--
========================================
Justin A. Lemkul
Ph.D. Candidate
ICTAS Doctoral Scholar
MILES-IGERT Trainee
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin
========================================
--
gmx-users mailing list gmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
