Hi
I did not even remember that I had this problem, or how I solved it.
But I can show you how I do the analysis.
For iRED:
trajin ../mdcrd_prod_gas
vector v2 :2.N ired :2.H
vector v3 :3.N ired :3.H
vector v4 :4.N ired :4.H
vector v6 :6.N ired :6.H
vector v7 :7.N ired :7.H
vector v8 :8.N ired :8.H
... (and so on for the other vectors, I have 128 of them)
matrix ired name matired order 2
analyze matrix matired vecs 128 out ired.vec
Then I get the eigenvectors and eigenvalues in the ired.vec file.
These you will need you in order the calculate the order parameters.
As you see, you do not have to remove the overall tumbling with RMS
fit, this will be taken care of when you not using the 5 largest
eigenvalues in the calculation of the order parameters. I have only
got this to work with Amber10, with Amber9 the analysis do not give me
all of the eigenvectors. Be aware of this, that you really obtain all
the vectors you are specifying.
I attach a python script that could be used to calculate the order
parameters from the ired.vec file. Just type: python mat2s2.py
ired.vec . No warranties, though.
For auto-correlation:
trajin ../mdcrd_prod_gas
rms first out rmsd.txt ".C,CA,N"
vector v2 :2.N corr :2.H
analyze timecorr vec1 v2 vec2 v2 tstep 1.0 tcorr 3500.0 out Ct/v2.out
vector v3 :3.N corr :3.H
analyze timecorr vec1 v3 vec2 v3 tstep 1.0 tcorr 3500.0 out Ct/v3.out
vector v4 :4.N corr :4.H
analyze timecorr vec1 v4 vec2 v4 tstep 1.0 tcorr 3500.0 out Ct/v4.out
vector v6 :6.N corr :6.H
analyze timecorr vec1 v6 vec2 v6 tstep 1.0 tcorr 3500.0 out Ct/v6.out
… (and so on for the other vectors)
Here you will get an auto-correlation file for each of the bond
vectors v2, v3 etc. I have also fitted the structures to the first
snapshot to remove overall tumbling. The trajectory is 35 ns long,
which is why have I choose to calculate the correlation function up to
3.5 ns (3500 ps). Due to statistical uncertainties in the tail of the
correlation function I never feel safe to calculate it more than up to
1/10 of the simulation length. From the auto-correlation functions
(ACF) I calculate the order parameters by fitting the ACF to either a
single or double exponential, and the order parameter is taken to be
the plateau value. This is where the trouble of convergence come in.
If the ACF is not converged there is no plateau and you cannot use fit
a proper exponential to your ACF. I have seen some people to take the
order parameters as an average over the last X values of the ACF. I
think that is a good check for comparing with the exponential fit. If
the two do not match, the ACF have not converged.
I seems to me that you have mixed up the commands a little bit, or I'm
not sure what you are trying to do. I think new practitioners (and I
know I would have) benefit from a tutorial on these matters, because
as David Case pointed out, the manual is not clear at all.
I have not tried the new relaxation option that David Case mentioned.
It would be interesting to look at, but again to calculate any
correlation function is a non-trivial, and you must make sure that you
have a sufficiently long simulation.
Best regards,
Samuel
I think you
2010/5/3 <luzhenw1.msu.edu>:
> Thanks a lot for all the advice from Paul, carols and Samuel. It seems the maximum correlation time is very important for the calculation.
> Now I tried to calculate the order parameter using IRED approach and had the exact same problem as Samuel.
> http://archive.ambermd.org/200806/0316.html
> I got the same autocorrelation function for all the NH vectors. Samuel, can you help out?
> Here is my input file: (only use three residue for simplicity)
>
> trajin test.mdcrd
> vector r2 :2.N corr :2.H
> vector r3 :3.N corr :3.H
> vector r9 :9.N corr :9.H
> matrix ired order 2 name nhma
> analyze matrix nhma out nhmr vecs 7
> vector r2 :2.N corrired :2.H modes nhmr begin 1 end 3 npair 1
> analyze timecorr vec1 r2 vec2 r2 out r2.nh tcorr 500. tstep 2.0 norm out r2.nh
> vector r3 :3.N corrired :3.H modes nhmr begin 1 end 3 npair 2
> analyze timecorr vec1 r3 vec2 r3 out r3.nh tcorr 500. tstep 2.0 norm out r3.nh
> vector r9 :9.N corrired :9.H modes nhmr begin 1 end 3 npair 3
> analyze timecorr vec1 r9 vec2 r9 out r9.nh tcorr 500. tstep 2.0 norm out r9.nh
>
> all the best
>
> thQuoting Samuel Genheden <samuel.genheden.gmail.com>:
>
>> Hi
>>
>> We've started to abandon the autocorrelation approach, just because
>> the problems with convergence. You have to look carefully on each
>> residue to see if they have really converged or not. Therefore, we
>> recommend to use iRED instead. On our test cases we haven't observed
>> any large overall difference between iRED and
>> autocorrelation-approach, but the latter requires more hands-on job.
>>
>> / Samuel
>>
>> 2010/5/1 Paul S. Nerenberg <psn.berkeley.edu>:
>>> Hi,
>>>
>>> Your e-mail is unclear in terms of what you mean by "negative value of the
>>> auto-correlation function". It is not unexpected that the auto-correlation
>>> function becomes negative at some correlation time (tau), although it should
>>> usually be fluctuating around zero when this happens (i.e., not reaching a
>>> large negative value). Moreover, if you were to average your C(tau)
>>> calculation over several trajectories, then you should see that these
>>> fluctuations about zero generally become smaller (as this "noise" is
>>> averaged out).
>>>
>>> So long as your simulations are adequate in length with regards to the max.
>>> correlation time and you are reasonably certain that your system can be
>>> analyzed with the model-free approach (which is nothing more than a *model*
>>> that makes assumptions about the dynamics of your molecule), then you should
>>> not have a problem fitting your auto-correlation functions to that
>>> exponential form.
>>>
>>> Best,
>>>
>>> Paul
>>>
>>> P.S. The max. correlation time used in your analysis (and correspondingly
>>> your simulation length) should always be chosen with regard to the
>>> characteristic relaxation times of your molecule, as is explicitly
>>> demonstrated in the article that Carlos referenced. I'm sure that the more
>>> NMR-minded AMBER users out there probably have some other pointers as
>>> well...
>>>
>>>
>>> On Apr 30, 2010, at 2:39 PM, luzhenw1.msu.edu wrote:
>>>
>>>> Thanks for the information. I read the Bruschweiler paper for the IRED
>>>> method and think it is OK to not use IRED approach. The question still
>>>> remains.
>>>>
>>>> Quoting Carlos Simmerling <carlos.simmerling.gmail.com>:
>>>>
>>>>> this is a useful article for such calculations
>>>>> http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2426661/
>>>>>
>>>>>
>>>>> On Fri, Apr 30, 2010 at 4:43 PM, <luzhenw1.msu.edu> wrote:
>>>>>
>>>>>> According to model-free formalism, the internal correlation function of
>>>>>> NH
>>>>>> vector, CI(t)=S^2+(1-S^2)e^(-t/tau)>=0. I tried to calculate the order
>>>>>> parameter from MD calculation. However, I got negative value of the
>>>>>> auto-correlation function for some residues. Is the defination of
>>>>>> autocorrelation different between NMR and MD? Thanks.
>>>>>> _______________________________________________
>>>>>> AMBER mailing list
>>>>>> AMBER.ambermd.org
>>>>>> http://lists.ambermd.org/mailman/listinfo/amber
>>>>>>
>>>>> _______________________________________________
>>>>> AMBER mailing list
>>>>> AMBER.ambermd.org
>>>>> http://lists.ambermd.org/mailman/listinfo/amber
>>>>>
>>>>>
>>>> _______________________________________________
>>>> AMBER mailing list
>>>> AMBER.ambermd.org
>>>> http://lists.ambermd.org/mailman/listinfo/amber
>>>
>>>
>>> _______________________________________________
>>> AMBER mailing list
>>> AMBER.ambermd.org
>>> http://lists.ambermd.org/mailman/listinfo/amber
>>>
>>
>> _______________________________________________
>> AMBER mailing list
>> AMBER.ambermd.org
>> http://lists.ambermd.org/mailman/listinfo/amber
>>
>>
> _______________________________________________
> AMBER mailing list
> AMBER.ambermd.org
> http://lists.ambermd.org/mailman/listinfo/amber
>
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
- application/octet-stream attachment: mat2s2.py
Received on Mon May 03 2010 - 22:30:04 PDT
