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

From: faride badalkhani (farideh.khamseh_at_gmail.com)
Date: Mon Aug 22 2016 - 07:34:17 CDT

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
>

This archive was generated by hypermail 2.1.6 : Sun Dec 31 2017 - 23:20:39 CST