[gmx-users] strange vdw energy from rerun with GBSA model.
lidawei at gmail.com
Thu Aug 11 16:31:37 CEST 2011
Please see my response below.
On Thu, Aug 11, 2011 at 10:23 AM, Mark Abraham <Mark.Abraham at anu.edu.au>wrote:
> On 11/08/2011 10:22 PM, 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.
> 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.
They are the same but different precision. I used a mdp file without PBC. I
believe this is reasonable because I can fully reproduce the bond length and
angle energy if I use trr file but cannot reproduce if I use .xtc file.
> 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.
> Sounds like you have a non-equilibrium trajectory (from your EM?), in which
> case small deviations can have large effects.
No. It is from a stable NPT trajectory with explicit water models. The
forces are around 10000 but are not larger.
> 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.
> 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.
That is why I say precesion in xtc file is not enough for bond length, angle
potential and vdw potential.
> 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
As Justin suggested, increase cutoff in rerun can remove the vdw energy
> 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)
> 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.
I am very sorry for this typo. I do use the nojump.xtc in my rerun.
> 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)
> 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
> 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
> ; 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.
> gmx-users mailing list gmx-users at gromacs.org
> 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-request at gromacs.org.
> Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
-------------- next part --------------
An HTML attachment was scrubbed...
More information about the gromacs.org_gmx-users