[gmx-users] RE: LIE methodology check

[Tom Dupree]

Greetings all,
I now realise I asked too many questions with that one post.
I would like to check that mdrun -rerun can be used to do single point energy calculations on the co-ordinate frames from my trajectory using a different electrostatic model by changing the .mdp file used by grompp to build the .tpr file? 

All the best,
Tom


___________________________________________

Message: 2
Date: Mon, 21 May 2012 03:26:07 +0000
From: Tom Dupree <t.dupree at unsw.edu.au>
Subject: [gmx-users] LIE methodology check (mdrun -rerun cutoffs and
	box	vectors)
To: "gmx-users at gromacs.org" <gmx-users at gromacs.org>
Message-ID:
	<0A0E67195A48DC4A96321968E2B384B024BF646D at INFPWXM007.ad.unsw.edu.au>
Content-Type: text/plain; charset="windows-1252"

Greetings all,

Background:
I have a series of MD runs of a protein ligand complex where I have used PME, they represent far more calculations than I have time to repeat in full. I recently became aware that the coulombic interaction energy reported when PME is used is not valid for use with g_lie (or any LIE method). I have read some relevant LIE papers and it seems that reaction field is the electrostatic method used.
http://pubs.acs.org/doi/abs/10.1021/jm7012198

Current methodology:
I found that using mdrun ?rerun and a new .tpr file I can change reported energies. Are these energies corresponding to the settings in the new .mdp file? Can I change it from PME to Generalised-Reaction-Field? Are there any errors associated with changing it, as opposed to errors implicit in either method? I am understanding ?rerun to only be calculated for the co-ordinate frames that are outputted into the .xtc file irrespective of any settings in the new .mdp file and/or old .mdp file for the energy.

As I only outputted .xtc files I don?t have velocities (or did I set that in the .mdp file?)  but I assume that errors in kinetic energy, temperature and pressure are not important in the rerun as I am only concerned with interaction energy?

Right assuming all the above is valid I would like to improve my accuracy as much as possible.
My understanding of cut offs is that the bigger they are the more accurate and the more computationally expensive your calculations are until they approach half the shortest box vector where pbc issues cause errors. As I am only doing reruns now I should use the largest possible cut offs that don?t hit the pbc limit?
Is there an accuracy difference between the short range and long range calculations for the LJ potential? (is there a reason to use long range at all?)
For the reruns can I ignore the problem of a water molecule seeing both sides of my protein and use a cut off greater than the  editconf -d setting (in my case 1.6)? Or to be strictly accurate should I work out the distance to the other side of the mirror protein from the binding site via the pbc boundary and use half that?

Which brings me to my issue with the unit cell. I set -bt dodecahedron. But I don't see it as one (vmd/pymol display of water atoms/unit cell). How do I relate the vectors in the.gro file to the pbc arrangement?

------------------------------------
2584CL      CL 7792   5.230   2.213   3.017  0.0075  0.3566  0.5809
2585CL      CL 7793   1.448   2.710   1.489  0.5771  0.2054  0.0146
2586CL      CL 7794   3.070   1.270   3.024 -0.1756 -0.0490 -0.1351
2587CL      CL 7795   5.112   2.558   1.562 -0.1543 -0.0592 -0.0072
2588CL      CL 7796   2.906   1.133   1.094 -0.2265 -0.3952  0.1190
   4.80953   4.80953   3.40085   0.00000   0.00000   0.00000   0.00000   2.40476   2.40476
----------------------

In my trials I got the error
?ERROR 1 [file topol.top, line 50]:
  ERROR: The cut-off length is longer than half the shortest box vector or
  longer than the smallest box diagonal element. Increase the box size or
  decrease rlist.?
So if I understand the .gro file correctly (now) these vectors  are the long diagonals of three rhombuses. Two are at right angles to each other and the third connects two a lower edge on each of those. Trial and error lead me to the maximum cut off that grompp will accept is 2.40476, which is d/2, if I read this correctly. My question is, why don't my waters and/or unit cell appear to be a dodecahedron when viewed in VMD/pymol? For that matter why don't the waters and unit cell agree on a shape when viewed in VMD? (both are rhombic prisims).

All the best,
Tom