[gmx-users] FEP energy errors with positional restraints?
Mark.Abraham at anu.edu.au
Sun Feb 13 06:31:04 CET 2011
On 13/02/2011 3:57 PM, TJ Mustard wrote:
> On February 12, 2011 at 5:35 PM Mark Abraham <Mark.Abraham at anu.edu.au>
>> On 13/02/2011 11:49 AM, TJ Mustard wrote:
>>> Hi all,
>>> I have been testing the ability of taking a sphere of a protein
>>> around a ligand, and positionally restrain the specified alpha
>>> carbons. I was hoping to keep non connected protein chains from
>>> drifting apart. I have been able to run these md/fep jobs, but I get
>>> huge interaction energies for the ligand, which has no positional
>>> restraints on it. I also don't restrain any atoms within the rvdw,
>>> rcoulomb and rlist radii. Am I thinking this is a possibility when
>>> it is physically impossible to simulate?
>> I doubt it.
>>> Currently I am selecting all residues around the ligand that have an
>>> atom within 20 Angstroms. I then save this as a pdb file and then
>>> run it through pdb2gmx, manually create a posres.itp file for each
>>> "chain" with their first and last residue's alpha carbon. Once I
>>> turn off these positional restraints the FEP energies drop down to
>>> "normal" levels.
>>> Does anyone have an idea what is happening?
>>> And if you do, can you please give a recommendation?
>> I'd guess you're not restraining how you think you are :-) Bear in
>> mind that position restraints are indexed relative to a
>> [moleculetype] (and must be #included there), and not the whole system.
> First I tried setting this globally it the .top file with the
> numbering corresponding to the atoms in question there, but found that
> the .top file only runs the solvent molecules and all of my reference
> atoms were outside the parameters (atoms 1-3 for water).
By default -DPOSRES will only restrain solvent because a) that's all
it's meant to do, because b) that [position_restraints] block is local
to the SOL [moleculetype], like I said.
> Then I setup the restraints argument in the Protein_chain_A.itp.itp
> file to reference the posre_Protein_ends_chain_A.itp file I made.
> ; Include Position restraint file
> #ifdef POSRES_PROTEIN
> #include "posre_Protein_ends_chain_A.itp"
> In the posre_Protein_ends_chain_A.itp file I entered the alpha-C atoms
> (2 total both the first residue and last) for the protein.
> [ position_restraints ]
> ; atom type fx fy fz
> 5 1 1000 1000 1000
> 262 1 1000 1000 1000
Well, that should work, if you put this in the right [moleculetype]
block. Whether that's enough of a restraint can't be said.
>> The combination of g_select and genrestr is probably a more reliable
>> and documentable way to generate your position restraints.
> I tried using genrestr directly and found that the selection was
> limited to the presets. I then tried to run a make_ndx to select the
> residues that were the starting and terminating ends, but I don't know
> how and if this program can do this task. I could manually find the
> residue numbers and input them directly but knowing that I was looking
> to restrain 6 atoms I decided to do it manually.
I thought you were trying to get all the alpha carbons outside a sphere,
sorry. For just 6 atoms, sure do it by hand.
> I have yet to look at g_select, as I am just seeing if this idea would
> work before I start making a automated script for this.
The point here is that it is straightforward to use g_select to "make an
index group of alpha carbons further than a given distance from some
location" and now you can use genrestr to make position restraints for
that whole group.
> Would there be a better way to hold these atoms in the respective
> locations. I could hold there distances constant.
Also possible, but (IIRC) the atoms have to be part of the same
[moleculetype], which is awkward for inter-chain restraints. Either way,
you have to address whether your restraints are perturbing the dynamics.
-------------- next part --------------
An HTML attachment was scrubbed...
More information about the gromacs.org_gmx-users