[gmx-users] FEP energy errors with positional restraints?

[Mark Abraham]

On 13/02/2011 5:03 PM, TJ Mustard wrote:
>
>
> On February 12, 2011 at 9:31 PM Mark Abraham <Mark.Abraham at anu.edu.au> 
> wrote:
>
>> 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> <mailto:Mark.Abraham at anu.edu.au> wrote:
>>>
>>>> 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"
>>> #endif
>>>
>>> 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.
>
> Well my major problem is the energy being transmitted about 12 
> angstroms away into my FEP molecule. I have consistently gotten an 
> interaction energy of around ~130 kJ/mol. With these restraints on the 
> energy spikes to ~125000 kJ/mol. Since I have my rvdw, rcoulomb, rlist 
> being 1.0 (nm), do I need to move these constraints father away from 
> the ligand?
>

Hence my thinking that you're somehow applying restraints to the wrong 
[moleculetype] and the FEP is amplifying the problem. Can you run normal 
MD normally, with and without restraints?

Mark
>>>
>>>> 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.
>
> I have just looked into this and I will hopefully be able to use this 
> for the other projects I have.
>
>>> 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.
>
> Well I will stay away from this for now then.
>
> Thank you,
>
> TJ Mustard
>
>> Mark
>
> TJ Mustard
> Email: mustardt at onid.orst.edu
>

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20110213/785bb1e8/attachment.html>