Thank you very much for your help, Justin. I tried it with those mdp options:
title = INITIAL_EM0 cpp = /lib/cpp constraints = all-bonds integrator = md rlist = 0.0 rcoulomb = 0.0 rvdw = 0.0 pbc = no This resulted in Angle Proper Dih. Improper Dih. LJ-14 Coulomb-14 1.88961e+01 6.04179e+01 5.60480e-04 3.74069e+01 -2.06394e+02 LJ (SR) Coulomb (SR) Potential Kinetic En. Total Energy 6.94305e-03 -5.66241e+01 -1.46290e+02 2.88457e-01 -1.46002e+02 Temperature Pressure (bar) Constr. rmsd 1.21730e+00 0.00000e+00 0.00000e+00 So, now the 1-4-values are slightly better (37.11 <-> 37.04 LJ, -207.7 <-> -206.3 Coulomb), while the SR-LJ and Coulomb energys are at two orders of magnitude too small (LJ: 0.7127 <-> 0.0069430, Coulomb: -4354.5 <-> -56.62). In fact, I considered looking through the code but wanted to avoid it because I fear that there are so much optimisations implemented in Gromacs that I probably will not understand whats going on and why. Maybe I will try it, but first I will play with the forcefield parameters to see what happens if I switch everything off except certain pairs, then I can compare them. Besides, I'm wondering, is there a way to disable all optimisations and just do a "plain" Lennard-Jones and Coulomb evaluation in vacuo without any box, solvent, boundary conditions and without cutoffs in Gromacs? Thanks and best regards, Matthias Am Montag, den 12.09.2011, 14:37 +0200 schrieb Justin A. Lemkul: > > Matthias Ernst wrote: > > Dear Gromacs users, > > > > I am currently working on porting forcefields from Gromacs to another > > program with the aim of simulating DNA. Therefore, I wrote a little > > script which reads in a PDB file, evaluates the Lennard-Jones and the > > Coulomb energy in vacuo with the forcefield data from Gromacs. I tried > > to implement the search for the neighbours by comparing inter-atomic > > distances to the van-der-Waals-radii and a flexible list-based factor > > weighting of the various contributions (where -1 is a default value, 0 > > is the atom itsself, 1 a bond to a next neighbour and so on, so my > > factor "3" should be equivalent to the "1-4-interaction" in Gromacs). > > > > To verify what I have done, I tried to do a calculation of a single "DA" > > nucleic acid with Gromacs to compare the LJ and Coulomb energy. > > Unfortunately, the results do not match and I do not understand why. > > > > Please find the PDB code I used here http://pastebin.com/zB3sQHdZ and > > the python code I used for evaluating the potentials here > > http://pastebin.com/2v4nVdg3 > > > > The results I got with my code are > > -> Total Lennard-Jones Energy [kJ/mol]: 0.712752493372 > > -> Total 1-4-Lennard-Jones Energy [kJ/mol]: 37.1121824156 > > -> Total Coulomb Energy [kJ/mol]: -4354.55160305 > > -> Total Coulomb 1-4 Energy [kJ/mol]: -207.712332451 > > > > Using Gromacs, I got > > Energies (kJ/mol) > > Angle Proper Dih. Improper Dih. LJ-14 Coulomb-14 > > 2.00975e+01 6.05485e+01 4.28869e-04 3.69509e+01 -2.06010e+02 > > LJ (SR) Coulomb (SR) Potential Pressure (bar) Constr. rmsd > > -1.45022e+00 -5.67518e+01 -1.46615e+02 0.00000e+00 1.80690e-04 > > > > The option in my mdp are > > title = INITIAL_EM0 > > cpp = /lib/cpp > > constraints = all-bonds > > integrator = steep > > nsteps = 50 > > nstcomm = 1 > > ns_type = grid > > rlist = 0.0 > > rcoulomb = 0.0 > > rvdw = 0.0 > > > > ; > > ; Energy minimizing stuff > > ; > > emtol = 1500.0 > > emstep = 0.01 > > > > Tcoupl = no > > Pcoupl = no > > gen_vel = no > > pbc = no > > > > I have not yet done a lot of calculations with Gromacs so maybe I did > > something wrong. If so, please correct me. > > > > You're better off with a single-point energy calculation, rather than a > 50-step > EM. Use the MD integrator with nsteps = 0 to compare configurations. If I > recall, steepest descents makes a move before step zero, so you're not > comparing > the same quantities. > > > What I wonder is that the 1-4-energies seem quite reasonable (LJ: 37.11 > > vs. 36.95, Coulomb -207.7 vs. -206.0) but the rest is rather different, > > especially the sign of the LJ term and the order of magnitude of the > > Coulomb term. Why is that? What does the "SR" in LJ and Coulomb in the > > Gromacs output mean? > > > > SR means short-range, i.e. within the cutoff. Since you've set all cutoffs > to > zero, this effectively means all interactions are calculated in this case. > > I can't offer an explanation for the different values observed, but I'll > suggest > you look into the Gromacs code to see how it calculates the energies and > compare > your program to see if you are, in fact, doing things the same way. > > -Justin > > -- > ======================================== > > 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/Support/Mailing_Lists/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/Support/Mailing_Lists
