[gmx-users] problem with position restraints: X0 set to zero

[Julian Garrec]

Justin A. Lemkul wrote:
> Julian Garrec wrote:
>   
>> Justin A. Lemkul wrote:
>>     
>>> Julian Garrec wrote:
>>>  
>>>       
>>>> Dear GROMACS users,
>>>>
>>>> I am trying to equilibrate my system (monomeric protein in water) and 
>>>> I want to use position restraint on heavy atoms of the solute using 
>>>> the posre.itp file. For some reason, grompp applies correctly the 
>>>> force constant I want, but sets all the reference positions to zero. 
>>>> Of course I would like this ref. positions to be taken from my .gro 
>>>> file:
>>>>
>>>> ${GROMPP} -debug 10 -f inp.mdp -po out.mdp -c first.gro -r first.gro 
>>>> -p prm.top -o allprm.tpr -maxwarn 10
>>>>
>>>>
>>>> Here are the first line of my posre.itp file:
>>>>
>>>> [ position_restraints ]
>>>> ; atom  type      fx      fy      fz
>>>>     1     1 567 567 567
>>>>
>>>>
>>>> ... and the first line containing the "POSRES" flag in my tpr file 
>>>> (dumped with gmxdump):
>>>>
>>>> functype[441]=POSRES, pos0A=( 0.00000000e+00, 0.00000000e+00, 
>>>> 0.00000000e+00), fcA=( 5.67000000e+02, 5.67000000e+02, 
>>>> 5.67000000e+02), pos0B=( 0.00000000e+00, 0.00000000e+00, 
>>>> 0.00000000e+00), fcB=( 5.67000000e+02, 5.67000000e+02, 5.67000000e+02)
>>>>
>>>>
>>>> Using the -debug option with grompp didn't provide me further 
>>>> information.
>>>>
>>>> Does anybody know how to cure this problem ?
>>>>
>>>>     
>>>>         
>>> I don't know what pos0A is doing, but much earlier in the gmxdump 
>>> output should've been something along the lines of:
>>>
>>> posres_xA[    0]={ 3.19200e+00,  4.72500e+00,  2.66900e+00}
>>>
>>> This line holds the coordinates of each of the restrained atoms.
>>>
>>> -Justin
>>>   
>>>       
>> Hello Justin,
>>
>> Thanks for your answer. I wanted to perform further tests before 
>> answering you.
>>
>> Actually, I do have such lines like:
>>
>> posres_xA[    0]={ 3.19200e+00,  4.72500e+00,  2.66900e+00}
>>
>> As far as I know, the lines containing the flag "functype" in the dumped 
>> .tpr file define all the functions that are used to compute the 
>> potential energy. Since position restraints involve additive terms of 
>> the form k(x-x_0)^2, I expected that both k and x0 are defined in the 
>> same "functype" line.
>>
>> What I didn't say in my previous mail is that minimizing my system with 
>> such restraints still induces strong deviation from the experimental 
>> structure. Infact if I apply a huge force constant, I can see that the 
>> restraints are properly applied. It seems that their is a strong 
>> mismatching between the experimental structure and the FF I use...
>>
>>     
>
> I can't replicate your problem.  I've done minimizations of different systems 
> with different force constants for the position restraints; they always come out 
> fine.
>
> It might help to know what's in your system, how you set it up (the sequence of 
> commands), what force field you're using, your .mdp file, and any other 
> information you think is necessary.  Also be sure to note any error messages or 
> warnings that the Gromacs tools print out for you.  They are usually quite helpful.
>
> -Justin
>   
Sorry, I realize my previous mail was not so clear.

Actually, I have done several extra tests and typing the command:
${GROMPP} -f inp.mdp -po out.mdp -c first.gro -r first.gro -p prm.top -o 
allprm.tpr -maxwarn 10
gmxdump -s allprm.tpr > TPRDUMP.txt
I got different results with gromacs 3.3.3 and 4.0.7:

GROMACS 3.3.3 gives:
functype[ ]=POSRES, pos0A=( 3.19000000e+00, 5.21700000e+00, 
4.43800000e+00), fcA=( 1.00000000e+06, 1.00000000e+06, 1.00000000e+06), 
pos0B=( 3.19000000e+00, 5.21700000e+00, 4.43800000e+00), fcB=( 
1.00000000e+06, 1.00000000e+06, 1.00000000e+06)

while GROMACS 4.0.7 gives:
functype[ ]=POSRES, pos0A=( 0.00000000e+00, 0.00000000e+00, 
0.00000000e+00), fcA=( 1.00000000e+06, 1.00000000e+06, 1.00000000e+06), 
pos0B=( 0.00000000e+00, 0.00000000e+00, 0.00000000e+00), fcB=( 
1.00000000e+06, 1.00000000e+06, 1.00000000e+06)

Since i) I got distorted structures after minimization and ii) I was 
expecting that the 3 first numbers (pos0A) were the values of x0 in the 
additional harmonic potentials, I concluded that something was wrong in 
the way I was trying to apply harmonix restraint with GROMAC 4.0.7.

Doing other tests, I found that the distortions I got come from a 
mismatching between the experimental structure and the FF I use. 
Basically, starting from this structure and trying to minimize it with 
this FF (+ the protonation pattern + starting water box and counter 
ions) the system is very far from the closest minimum and experiences 
huuuuuuuuuuuuuuuudge forces. Only a force constant of 100000 kJ/mol.A^2 
is able to maintain the protein in a reasonable way. So clearly, my 
problem does not come from GROMACS but from my setup or the experimental 
structure.

Sorry again and thank for your reply.

Best regards,

Julian


>   
>> Thanks again,
>>
>> Julian
>>
>>     
>>>  
>>>       
>>>> Thanks,
>>>>
>>>> Julian
>>>>     
>>>>         
>>>   
>>>       
>>     
>
>