Re: Autocorrelation function of Rg goes down without fluctuation around zero

From: faride badalkhani (farideh.khamseh_at_gmail.com)
Date: Wed Sep 07 2016 - 01:14:08 CDT

Hi,

With "very good" mean it fluctuate around zero that is a sign of getting
equilibrated, as I have read in publications. But, the autocorrelation
function goes down at 20-30 ns. You could find the the Rg, autocorrelation
of Rg, and RMSD in the following link:
https://www.dropbox.com/sh/ol53056nvgplpws/AAD_6s60gbntQ8rXK99CBhl_a?dl=0

Regards,
Farideh

On Wed, Sep 7, 2016 at 2:13 AM, Brian Radak <bradak_at_anl.gov> wrote:

> I'm not sure I understand the question. How do you mean "very good" here?
> What "goes down"? Are you referring to the correlation function at 20-30
> ns? Remember, the observable that goes into a correlation function is the
> time lag, so that you only have one observation at the full length of your
> simulation - the function should get progressively noisier as time goes on.
> I would usually only compute the function out to half the simulation
> length, although that is quite arbitrary a time.
>
> A rule of thumb that I have heard (although I could not give the exact
> justification from memory) is that you should simulate out to ~1000 times
> the characteristic time (i.e. twice the integral of the correlation
> function). That is no guarantee on accuracy though.
>
> HTH,
>
> Brian
>
> On 09/05/2016 06:03 AM, faride badalkhani wrote:
> Hi,
>
> I have run a simulation for 2 ns in NVT ensemble and for 30 ns in NPT
> ensemble. When I check the autocorrelation function of Rg for this system
> it is very good for the first 20 ns but after that (20-30 ns) it goes down.
> I want to know is there any optimum time for a simulation run?
> Should I consider the data for the first 20 ns or continue the simulation
> for more time?
>
> Regards
> Farideh
>
> On Tue, Aug 23, 2016 at 11:58 PM, Brian Radak <brian.radak.accts_at_gmail.com
> > wrote:
>
>> On 08/23/2016 02:16 PM, faride badalkhani wrote:
>>
>> Hi,
>> I mean did I include all the necessary parameters in the NVT
>> configuration file?
>> As you said I performed a NVT (for 300 ps) then NPT (for 300 ps) and then
>> NPT for 20 ns (it is running now) and I have some questions:
>> 1) when I performed NVT simulation water box shape changed to a semi
>> spherical one but after performing NPT for 300 ps it got the cubic shape
>> again. Could you tell me if what happened in NVT should not have happened?
>>
>> I could not say. This all depends on your boundary conditions and the
>> Boltzmann factor of the initial configuration (is it a reasonable structure
>> for these parameters?).
>>
>> 2) what is the difference between config files we execute for first NPT
>> (the short one) and the final NPT run? For example, in GROMACS for NVT and
>> NPT we impose some forces on atoms and then for production run we remove
>> all the forces and atoms can move. Is there something similar in NAMD, too?
>>
>> That all depends. If you want to include a restraint or constraint
>> protocol, that is separate from NpT. Generally constraints can be confusing
>> for NpT because the virial must be modified - in my opinion restraints are
>> almost always better.
>>
>> Your help is really appreciated.
>>
>> Regards,
>> Farideh
>>
>> On Aug 23, 2016 11:29 PM, "Brian Radak" <brian.radak.accts_at_gmail.com>
>> wrote:
>>
>>> I'm not sure I understand the question. How do you mean "very good"
>>> here? What "goes down"? Are you referring to the correlation function at
>>> 20-30 ns? Remember, the observable that goes into a correlation function is
>>> the time lag, so that you only have one observation at the full length of
>>> your simulation - the function should get progressively noisier as time
>>> goes on. I would usually only compute the function out to half the
>>> simulation length, although that is quite arbitrary a time.
>>>
>>> A rule of thumb that I have heard (although I could not give the exact
>>> justification from memory) is that you should simulate out to ~1000 times
>>> the characteristic time (i.e. twice the integral of the correlation
>>> function). That is no guarantee on accuracy though.
>>>
>>> HTH,
>>>
>>> Brian
>>>
>>> On 09/05/2016 06:03 AM, faride badalkhani wrote:
>>>
>>> Hi,
>>>
>>> I have run a simulation for 2 ns in NVT ensemble and for 30 ns in NPT
>>> ensemble. When I check the autocorrelation function of Rg for this system
>>> it is very good for the first 20 ns but after that (20-30 ns) it goes down.
>>> I want to know is there any optimum time for a simulation run?
>>> Should I consider the data for the first 20 ns or continue the
>>> simulation for more time?
>>>
>>> Regards
>>> Farideh
>>>
>>> On Tue, Aug 23, 2016 at 11:58 PM, Brian Radak <
>>> brian.radak.accts_at_gmail.com> wrote:
>>>
>>>> On 08/23/2016 02:16 PM, faride badalkhani wrote:
>>>>
>>>> Hi,
>>>> I mean did I include all the necessary parameters in the NVT
>>>> configuration file?
>>>> As you said I performed a NVT (for 300 ps) then NPT (for 300 ps) and
>>>> then NPT for 20 ns (it is running now) and I have some questions:
>>>> 1) when I performed NVT simulation water box shape changed to a semi
>>>> spherical one but after performing NPT for 300 ps it got the cubic shape
>>>> again. Could you tell me if what happened in NVT should not have happened?
>>>>
>>>> I could not say. This all depends on your boundary conditions and the
>>>> Boltzmann factor of the initial configuration (is it a reasonable structure
>>>> for these parameters?).
>>>>
>>>> 2) what is the difference between config files we execute for first NPT
>>>> (the short one) and the final NPT run? For example, in GROMACS for NVT and
>>>> NPT we impose some forces on atoms and then for production run we remove
>>>> all the forces and atoms can move. Is there something similar in NAMD, too?
>>>>
>>>> That all depends. If you want to include a restraint or constraint
>>>> protocol, that is separate from NpT. Generally constraints can be confusing
>>>> for NpT because the virial must be modified - in my opinion restraints are
>>>> almost always better.
>>>>
>>>> Your help is really appreciated.
>>>>
>>>> Regards,
>>>> Farideh
>>>>
>>>> On Aug 23, 2016 11:29 PM, "Brian Radak" <brian.radak.accts_at_gmail.com>
>>>> wrote:
>>>>
>>>>> How do you mean "correct"? Does the file execute successfully? It can
>>>>> be very useful to try this and then read the INFO header lines.
>>>>>
>>>>> NpT generally requires *more* options than NVT. So long as you set
>>>>> "langevin on" (or some other proper thermostat) and do not set
>>>>> "langevinPiston on", you should get NVT results.
>>>>>
>>>>> On 08/22/2016 07:34 AM, faride badalkhani wrote:
>>>>>
>>>>> Dear NAMD users,
>>>>>
>>>>> I want to perform a sequence of equilibration at NVT, then NPT, and
>>>>> then production at NPT. I have some questions:
>>>>> 1- I have never performed NVT simulations before, so could you tell me
>>>>> if the following configuration file is correct?
>>>>>
>>>>> #############################################################
>>>>> ## JOB DESCRIPTION ##
>>>>> #############################################################
>>>>>
>>>>> # Minimization and Equilibration of
>>>>> # G5 OH terminated dendrimer in a Water Box
>>>>> # NVT ensemble
>>>>>
>>>>>
>>>>> #############################################################
>>>>> ## ADJUSTABLE PARAMETERS ##
>>>>> #############################################################
>>>>>
>>>>> structure ./OH-OH/input/OH_wb.psf
>>>>> coordinates ./OH-OH/input/OH_wb.pdb
>>>>>
>>>>> set temperature 300
>>>>> set outputname ./OH-OH/output/OH_eq
>>>>>
>>>>> firsttimestep 0
>>>>>
>>>>>
>>>>> #############################################################
>>>>> ## SIMULATION PARAMETERS ##
>>>>> #############################################################
>>>>>
>>>>> # Input
>>>>> paraTypeCharmm on
>>>>> parameters ./OH-OH/input/par_all36_OH.prm
>>>>> temperature $temperature
>>>>>
>>>>>
>>>>> # Force-Field Parameters
>>>>> exclude scaled1-4
>>>>> 1-4scaling 1.0
>>>>> cutoff 12.0
>>>>> switching on
>>>>> switchdist 10.0
>>>>> pairlistdist 14.0
>>>>>
>>>>>
>>>>> # Integrator Parameters
>>>>> timestep 1.0 ;# 1 fs/step
>>>>> rigidBonds water ;# needed for 1 fs steps
>>>>> nonbondedFreq 1
>>>>> fullElectFrequency 2
>>>>> stepspercycle 10
>>>>>
>>>>>
>>>>> # Constant Temperature Control
>>>>> langevin on ;# do langevin dynamics
>>>>> langevinDamping 5 ;# damping coefficient (gamma) of 5/ps
>>>>> langevinTemp $temperature
>>>>> langevinHydrogen off ;# don't couple langevin bath to hydrogens
>>>>>
>>>>>
>>>>> # Periodic Boundary Conditions
>>>>> cellBasisVector1 73.5 0.0 0.0
>>>>> cellBasisVector2 0.0 86.5 0.0
>>>>> cellBasisVector3 0.0 0.0 88.4
>>>>> cellOrigin -12.7 -12.5 -26.3
>>>>>
>>>>> wrapAll on
>>>>>
>>>>>
>>>>> # PME (for full-system periodic electrostatics)
>>>>> PME yes
>>>>>
>>>>>
>>>>> # let NAMD determine grid
>>>>> PMEGridSpacing 1.0
>>>>>
>>>>>
>>>>> useGroupPressure yes ;# needed for rigidBonds
>>>>> useFlexibleCell no
>>>>> useConstantArea no
>>>>>
>>>>>
>>>>> # Output
>>>>> outputName $outputname
>>>>>
>>>>>
>>>>> restartfreq 5000 ;# 5000 steps = every 5.0 ps
>>>>> dcdfreq 5000
>>>>> xstFreq 5000
>>>>> outputEnergies 5000
>>>>> outputPressure 5000
>>>>>
>>>>>
>>>>> #############################################################
>>>>> ## EXTRA PARAMETERS ##
>>>>> #############################################################
>>>>>
>>>>>
>>>>> #############################################################
>>>>> ## EXECUTION SCRIPT ##
>>>>> #############################################################
>>>>>
>>>>> # Minimization
>>>>> minimize 20000
>>>>> reinitvels $temperature
>>>>>
>>>>> run 300000 ;# 300 ps
>>>>>
>>>>>
>>>>> 2- It would be necessary to set "rigidBonds water " when time
>>>>> step is 1 fs?
>>>>>
>>>>> 3- How many time steps are standard values for the first NVT and NPT
>>>>> simulations? (for example in Gromacs it is recommended to run NVT and NPT
>>>>> for 300 ps if someone wants to run the priduction run for 30 ns).
>>>>>
>>>>> Regards,
>>>>> Farideh
>>>>>
>>>>> On Thu, Aug 18, 2016 at 9:07 PM, Brian Radak <bradak_at_anl.gov> wrote:
>>>>>
>>>>>> I would think the standard sequence is equilibrate at NpT, then NVT,
>>>>>> and then production at NVT. Some people do things differently, and that's
>>>>>> fine. There's nothing wrong with that - it just might be a bit slower (you
>>>>>> could even check this during equilibration). It's up to you.
>>>>>>
>>>>>> Restarting probably shouldn't affect results. It might if you were
>>>>>> trying to get very precise Newtonian dynamics, but there are plenty of
>>>>>> other concerns in that case also.
>>>>>>
>>>>>> On 08/18/2016 11:31 AM, faride badalkhani wrote:
>>>>>>
>>>>>> And I have one more questions! Could you tell me if the restarting
>>>>>> the simulations every 5 ns affects the results or not?
>>>>>>
>>>>>> Regards,
>>>>>> Farideh
>>>>>>
>>>>>> On Aug 18, 2016 8:57 PM, "faride badalkhani" <
>>>>>> farideh.khamseh_at_gmail.com> wrote:
>>>>>>
>>>>>>> So, do you think it is more reasonable to perform minimization and
>>>>>>> then NVT simulation for a few time steps before production run?
>>>>>>>
>>>>>>> Regards,
>>>>>>> Farideh
>>>>>>>
>>>>>>> On Aug 18, 2016 8:50 PM, "Brian Radak" <bradak_at_anl.gov> wrote:
>>>>>>>
>>>>>>>> After equilibration to a desired density, pressure control is
>>>>>>>> usually just an added cost (maybe 5-10%?) with not much value. I usually do
>>>>>>>> everything at NVT. For equilibration at the force field density, these are
>>>>>>>> usually indistinguishable (unless you want temperature or pressure
>>>>>>>> dependent properties).
>>>>>>>> On 08/18/2016 11:16 AM, faride badalkhani wrote:
>>>>>>>>
>>>>>>>> Usually, the autocorrelation function of Rg is used to find if the
>>>>>>>> structure is well equilibrated or not, and to find relaxation time for
>>>>>>>> them.
>>>>>>>> What do you mean with "It's also not common to use pressure control
>>>>>>>> for these kinds of calculations"? Which method do you suggest for
>>>>>>>> perfotming simulations?
>>>>>>>>
>>>>>>>> Regards,
>>>>>>>> Farideh
>>>>>>>>
>>>>>>>> On Aug 18, 2016 8:30 PM, "Brian Radak" <bradak_at_anl.gov> wrote:
>>>>>>>>
>>>>>>>>> I've not used either of "measure rgyr" or "g_analyze", so I can
>>>>>>>>> only assume they are correctly implemented.
>>>>>>>>>
>>>>>>>>> I have pretty significant reservations about including your
>>>>>>>>> equilibration components in the calculation. It's also not common to use
>>>>>>>>> pressure control for these kinds of calculations, but I don't have a
>>>>>>>>> specific reason as to why that's bad other than that it changes the
>>>>>>>>> physical meaning of the result (I don't know what you want C(t) for, I
>>>>>>>>> assume this is not a problem).
>>>>>>>>>
>>>>>>>>> My guess is that you just have very poor statistics due to a long
>>>>>>>>> characteristic time. This is not unexpected for a (presumably large) system
>>>>>>>>> like a dendrimer.
>>>>>>>>>
>>>>>>>>> Cheers,
>>>>>>>>>
>>>>>>>>> Brian
>>>>>>>>>
>>>>>>>>> On 08/18/2016 10:31 AM, faride badalkhani wrote:
>>>>>>>>>
>>>>>>>>> This system is a Dendrimer, which is a hyperbranched polymer. I
>>>>>>>>> used the following script to calculate the radius of gyration:
>>>>>>>>>
>>>>>>>>> set outfile [open rg.dat w]
>>>>>>>>> set mol [molinfo top]
>>>>>>>>> set sel [atomselect top "all not waters"]
>>>>>>>>> set frames [molinfo $mol get numframes]
>>>>>>>>> for {set i 0} {$i < $frames} {incr i} {
>>>>>>>>> $sel frame $i
>>>>>>>>> $sel update
>>>>>>>>> puts $outfile "$i [measure rgyr $sel]"
>>>>>>>>> }
>>>>>>>>> close $outfile
>>>>>>>>> $sel delete
>>>>>>>>>
>>>>>>>>> For calculating the C(t) I used the g_analyze in gromacs and the
>>>>>>>>> Rg computed using the above script. for plotting the data I used the whole
>>>>>>>>> simulation time (40 ns). I loaded all the trajectory files in VMD and then
>>>>>>>>> calculated the Rg.
>>>>>>>>>
>>>>>>>>> Regards,
>>>>>>>>> Farideh
>>>>>>>>>
>>>>>>>>> On Thu, Aug 18, 2016 at 7:51 PM, Brian Radak <bradak_at_anl.gov>
>>>>>>>>> wrote:
>>>>>>>>>
>>>>>>>>>> What system is this? How did you compute Rg and C(t)? Which parts
>>>>>>>>>> of the trajectory were used?
>>>>>>>>>>
>>>>>>>>>> Correlation functions can be very slow to converge, 5 ns might
>>>>>>>>>> not be enough.
>>>>>>>>>>
>>>>>>>>>> On 08/18/2016 10:16 AM, faride badalkhani wrote:
>>>>>>>>>>
>>>>>>>>>> Thank you for your answer and time. At first I performed
>>>>>>>>>> minimization for 20 ps and then NPT simulation for 10 ns. After that, I
>>>>>>>>>> restarted the simulation for 5 more ns and continued simulation in 5 ns
>>>>>>>>>> steps till 40 ns. I did not get any error during the simulations, but the
>>>>>>>>>> autocorrelation function does not fluctuate around zero. I have shared the
>>>>>>>>>> Rg, RMSD, and C(t) ans also .top, .par, .conf, .pdb, and the .log file for
>>>>>>>>>> the first 10 ns in the link below:
>>>>>>>>>>
>>>>>>>>>> https://www.dropbox.com/sh/mcnkxt9pefnqgrz/AADRRw4-dJEEx0SBl
>>>>>>>>>> U2Fuekma?dl=0
>>>>>>>>>>
>>>>>>>>>> p.s. Bond distances, angles, dihedrals and impropers were taken
>>>>>>>>>> form CHARMM GENERAL FORCE FIELD and the charges were assigned using CGen FF
>>>>>>>>>> on paramchem.
>>>>>>>>>>
>>>>>>>>>> Thank you so much for your help!
>>>>>>>>>> Regards,
>>>>>>>>>> Farideh
>>>>>>>>>>
>>>>>>>>>> On Thu, Aug 18, 2016 at 7:17 PM, Brian Radak <bradak_at_anl.gov>
>>>>>>>>>> wrote:
>>>>>>>>>>
>>>>>>>>>>> You either made a mistake in your calculation or else your
>>>>>>>>>>> statistics are bad. Some details might help distinguish which of those is
>>>>>>>>>>> the case.
>>>>>>>>>>>
>>>>>>>>>>> HTH,
>>>>>>>>>>> Brian
>>>>>>>>>>>
>>>>>>>>>>>
>>>>>>>>>>> On 08/18/2016 04:21 AM, faride badalkhani wrote:
>>>>>>>>>>>
>>>>>>>>>>>> Dear NAMD users,
>>>>>>>>>>>>
>>>>>>>>>>>> I have plotted the autocorrelation function of radius of
>>>>>>>>>>>> gyration as a function of time to investigate the equilibration of a system
>>>>>>>>>>>> of hyperbranched polymer. However, the plot goes down and does not
>>>>>>>>>>>> fluctuate around zero.
>>>>>>>>>>>>
>>>>>>>>>>>> Could anybody tell me what is the reason and what should I do?
>>>>>>>>>>>> Any help will be appreciated.
>>>>>>>>>>>>
>>>>>>>>>>>> Regards,
>>>>>>>>>>>> Farideh
>>>>>>>>>>>>
>>>>>>>>>>>
>>>>>>>>>>> --
>>>>>>>>>>> Brian Radak
>>>>>>>>>>> Postdoctoral Appointee
>>>>>>>>>>> Leadership Computing Facility
>>>>>>>>>>> Argonne National Laboratory
>>>>>>>>>>>
>>>>>>>>>>> 9700 South Cass Avenue, Bldg. 240
>>>>>>>>>>> Argonne, IL 60439-4854
>>>>>>>>>>> (630) 252-8643
>>>>>>>>>>> brian.radak_at_anl.gov
>>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>> --
>>>>>>>>>> Brian Radak
>>>>>>>>>> Postdoctoral Appointee
>>>>>>>>>> Leadership Computing Facility
>>>>>>>>>> Argonne National Laboratory
>>>>>>>>>>
>>>>>>>>>> 9700 South Cass Avenue, Bldg. 240
>>>>>>>>>> Argonne, IL 60439-4854
>>>>>>>>>> (630) 252-8643
>>>>>>>>>> brian.radak_at_anl.gov
>>>>>>>>>>
>>>>>>>>>
>>>>>>>>>
>>>>>>>>> --
>>>>>>>>> Brian Radak
>>>>>>>>> Postdoctoral Appointee
>>>>>>>>> Leadership Computing Facility
>>>>>>>>> Argonne National Laboratory
>>>>>>>>>
>>>>>>>>> 9700 South Cass Avenue, Bldg. 240
>>>>>>>>> Argonne, IL 60439-4854
>>>>>>>>> (630) 252-8643
>>>>>>>>> brian.radak_at_anl.gov
>>>>>>>>>
>>>>>>>>
>>>>>>>> --
>>>>>>>> Brian Radak
>>>>>>>> Postdoctoral Appointee
>>>>>>>> Leadership Computing Facility
>>>>>>>> Argonne National Laboratory
>>>>>>>>
>>>>>>>> 9700 South Cass Avenue, Bldg. 240
>>>>>>>> Argonne, IL 60439-4854
>>>>>>>> (630) 252-8643
>>>>>>>> brian.radak_at_anl.gov
>>>>>>>>
>>>>>>>
>>>>>> --
>>>>>> Brian Radak
>>>>>> Postdoctoral Appointee
>>>>>> Leadership Computing Facility
>>>>>> Argonne National Laboratory
>>>>>>
>>>>>> 9700 South Cass Avenue, Bldg. 240
>>>>>> Argonne, IL 60439-4854
>>>>>> (630) 252-8643
>>>>>> brian.radak_at_anl.gov
>>>>>>
>>>>>
>>>>>
>>>>> --
>>>>> Brian Radak
>>>>> Postdoctoral Scholar
>>>>> Gordon Center for Integrative Science, W323A
>>>>> Department of Biochemistry & Molecular Biology
>>>>> University of Chicago
>>>>> 929 E. 57th St.
>>>>> Chicago, IL 60637-1454
>>>>> Tel: 773/834-2812
>>>>> email: radak_at_uchicago.edu
>>>>>
>>>>
>>>> --
>>>> Brian Radak
>>>> Postdoctoral Scholar
>>>> Gordon Center for Integrative Science, W323A
>>>> Department of Biochemistry & Molecular Biology
>>>> University of Chicago
>>>> 929 E. 57th St.
>>>> Chicago, IL 60637-1454
>>>> Tel: 773/834-2812
>>>> email: radak_at_uchicago.edu
>>>>
>>>
>>>
>>> --
>>> Brian Radak
>>> Postdoctoral Appointee
>>> Leadership Computing Facility
>>> Argonne National Laboratory
>>>
>>> 9700 South Cass Avenue, Bldg. 240
>>> Argonne, IL 60439-4854
>>> (630) 252-8643
>>> brian.radak_at_anl.gov
>>>
>>

This archive was generated by hypermail 2.1.6 : Tue Dec 27 2016 - 23:22:26 CST