Moeed wrote:
Hello Justin,
[Actually, the reason I posted some parts of my mdp file was because you
said I am sending redundant information so I just deleted the first
lines of mdp file :) ]
What I was referring to was the meaningless grompp outputs, EM results, etc. An
.mdp file should always be posted in full if it is relevant.
To get the xtc file initially I used nrexcl=3 in the top file. Then I
used mdrun -rerun to breakdown energies on this xtc file with new top
file nrexcl=7 to exclude all intramolecular nonbonded interactions (you
definitely recall all these details :)
1- Does it matter in the first run (to get xtc) what option I use for
nrexcl in top file. Actually, I thought maybe in the normal run I have
to set this option to zero. I did that and I got lincs warning (the same
applies for nrexcl=1). So I tried nrexcl=2 and I noticed the results are
not the same as nrexcl 3. However, this means I am still having
nonbonded intramolecular interaction between atoms that are less than 2
bonds away. Am I right?
Haphazardly setting different values for nrexcl will break down the physics!
You must understand the intrinsic assumptions in most MD force fields before
playing with these values. If you set nrexcl = 1 or 2 then nonbonded
interactions are being calculated for atom pairs that are normally governed only
by bonded interactions. Anything less than 3 (or greater than three in most
cases, for that matter) is going to lead to a breakdown and unreasonable
results, if you manage to get the simulation running at all. That's why your
nrexcl = 7 has to be done as a rerun. If you tried this kind of simulation on
its own, the results would be gibberish.
2- I am wondering why in this breakdown I still see LJ and coloumb 1-4
interactions while I have nrexcl=2.
Because 1-4 interactions are separated by three bonds. Using nrexcl = 2 has no
effect here.
3- In the second set of results (nrexcl=3) potential energy is 2766 >
2156 (the value for nrexcl==2). Logically, when 3 bonds away are
excluded potential energy should be less I think!..
Hard to say. You're completely changing the balance of power between bonded and
nonbonded interactions.
4- Regarding the units, in chapter 2 of manual the units for energy is
in KJ/mol. Here is this table I see heat capacity is given in J/mol K.
Could you please kindly tell me the units what is the unit for energy
temrs LJ or coulomb?
The heat capacity you get is not useful. I believe the method for calculating
it has been updated in the development code, but basically you'll get the same
result (or nearly so) for just about every simulation that you run, regardless
of what's in it :)
5- Very important question: If I want to compute the interaction
energies BETWEEN molecuels, do I have to sum up LJ (SR) =(-3026 e.g) and
coloumb (SR)=74 ? and can I conclude that there is a net attraction
energy in the system?
That depends. Are there long range interactions? PME? Dispersion correction?
6- If I want to assure that equilibrium is reached, can I make the
simulation time long enough such that T reaches 300 K? now it is 297 K.
You can make run for as much time as you want. Sounds to me like you're
relatively close to where you need to be, but you've got to reach the target
temperature, ensure that it's stable, and then and only then should you be
collecting data.
-Justin
I appreciate your help.
Thanks
Moeed
_mdrun -rerun on trajectory for nrexcl=2 in top file_
Energy Average RMSD Fluct. Drift
Tot-Drift
-------------------------------------------------------------------------------
Angle 4286.19 258.765 231.378 40.1269
401.349
Ryckaert-Bell. 722.966 69.0616 61.1677 11.1046
111.068
*LJ-14 498.432 31.2891 30.9138 1.67344
16.7378*
* *
*Coulomb-14 -322.264 11.5056 10.1507
1.87612 18.765*
* *
*LJ (SR) -3026.03 45.5198 38.6106 8.35011
83.5178*
* *
*Coulomb (SR) 74.4564 2.53605 2.53537 0.02034
0.203441*
Coul. recip. -77.5276 2.3709 1.31777 -0.682623
-6.8276
Potential 2156.22 335.429 282.807 62.4689
624.814
Kinetic En. 0 0 0
0 0
Total Energy 2156.22 335.429 282.807 62.4689
624.814
Temperature 0 0 0
0 0
Pressure (bar) 1563.28 161.924 161.475 4.17097
41.718
Cons. rmsd () 0 0 0
0 0
Mu-X -0.00385475 0.893851 0.893682 -0.00602335
-0.0602456
Heat Capacity Cv: nan J/mol K (factor = nan)
_
_
_mdrun -rerun on trajectory for nrexcl=3 in top file_
Statistics over 5001 steps [ 0.0000 thru 10.0000 ps ], 14 data sets
All averages are over 501 frames
Energy Average RMSD Fluct. Drift
Tot-Drift
-------------------------------------------------------------------------------
Angle 4400.09 284.259 237.268 54.2184
542.293
Ryckaert-Bell. 1009.62 124.468 94.0379 28.2417
282.474
*LJ-14 648.411 40.6844 39.926 2.70793
27.0848*
**
*Coulomb-14 -278.654 34.4227 11.8071
11.1987 112.01*
**
*LJ (SR) -3010.16 47.968 42.5388
7.67736 76.789*
**
*Coulomb (SR) 74.7846 2.5186 2.47007 0.170417
1.70451*
Coul. recip. -78.0436 2.58606 1.54725 -0.717662
-7.17805
Potential 2766.05 454.268 342.141 103.497
1035.18
Kinetic En. 0 0 0
0 0
Total Energy 2766.05 454.268 342.141 103.497
1035.18
Temperature 0 0 0
0 0
Pressure (bar) 1973.01 177.068 172.803 13.3792
133.819
Cons. rmsd () 0 0 0
0 0
Mu-X -0.0523974 0.826473 0.826318 -0.00553723
-0.0553834
Heat Capacity Cv: nan J/mol K (factor = nan)
; Run control
integrator = md
dt = 0.002 ; ps !
nsteps = 5000 ; total 1.0 ps.
nstcomm = 1 ; frequency for center of mass motion
removal
; Output control
nstenergy = 10 ; frequency to write energies to energy
file. i.e., energies and other statistical data are stored every 10 steps
nstxout = 10 ; frequency to write
coordinates/velocity/force to output trajectory file
nstvout = 0
nstfout = 10
nstlog = 10 ; frequency to write energies to log file
nstxtcout = 10
; Neighbor searching
nstlist = 10 ; neighborlist will be updated at least
every 10 steps
;ns_type = grid
; Electrostatics/VdW
coulombtype = PME
vdw-type = cut-off
; Cut-offs
rlist = 1.0
rcoulomb = 1.0
rvdw = 1.0
; Temperature coupling Berendsen temperature coupling is on in
two groups
Tcoupl = berendsen
tc-grps = HEX ;sol
tau_t = 0.1 ;0.1
ref_t = 300 ;300
; Pressure coupling: Pressure coupling is not on
Pcoupl = no
tau_p = 0.5
compressibility = 4.5e-5
ref_p = 1.0
; Velocity generation Generate velocites is on at 300 K.
Manual p155
gen_vel = yes
gen_temp = 300.0
gen_seed = 173529
; Bonds
constraints = all-bonds
constraint-algorithm = lincs
pbc=xyz
--
========================================
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 [email protected]
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at http://www.gromacs.org/search before posting!
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to [email protected].
Can't post? Read http://www.gromacs.org/mailing_lists/users.php
