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

From: faride badalkhani (farideh.khamseh_at_gmail.com)
Date: Mon Sep 05 2016 - 06:03:13 CDT

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
>

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