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. 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? Any help would be thankfully apreciated. Best regards, Matthias Ernst -- 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
