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

From: Brian Radak (bradak_at_anl.gov)
Date: Tue Sep 06 2016 - 16:43:30 CDT

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 <mailto: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
>> <mailto: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
>>> <mailto: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
>>>> <mailto: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 <mailto: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 <mailto: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
>>>>>> <mailto: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-dJEEx0SBlU2Fuekma?dl=0
>>>>>>> <https://www.dropbox.com/sh/mcnkxt9pefnqgrz/AADRRw4-dJEEx0SBlU2Fuekma?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
>>>>>>> <mailto: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
>>>>>>> <mailto: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
>>>>>> <mailto: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
>>>>> <mailto: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 <mailto: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 <mailto: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 <mailto: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 <mailto: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