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

[Mark Abraham]

On 13/02/2011 6:07 PM, TJ Mustard wrote:
>
>
> On February 12, 2011 at 10:14 PM Mark Abraham 
> <Mark.Abraham at anu.edu.au> wrote:
>
>> 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> <mailto: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?
>
> I have ran these jobs without the restraints, and the FEP energy is 
> ~130kJ/mol. I can run the md (actually in sd) but I am not sure what 
> you mean.
>

When trouble-shooting, eliminating sources of problems is critical. If 
you can run a non-FEP simulation with your position restraints, then 
they are probably not the source of the problem. If FEP without 
restraints works, and non-FEP with restraints works, and FEP with 
restraints doesn't, then someone can work out what to try next. Until 
then, the simplest hypothesis is that your restraints don't work (either 
by design or implementation flaws). :-)

Mark
>
> Thank you,
>
> TJ Mustard
>
>> 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 <mailto:mustardt at onid.orst.edu>
>>>
>>
> 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/3fedeedf/attachment.html>