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. dawei
-- 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