Re: Pressure Discrepancy

From: Peter Freddolino (petefred_at_ks.uiuc.edu)
Date: Wed Mar 28 2007 - 13:40:36 CDT

Hi Morad,
I saw that you had a piston period of 200 fs, but decay time of 500 fs.
As noted at http://www.ks.uiuc.edu/Research/namd/2.6/ug/node32.html,
this will result in underdamped dynamics, and thus pressure control of
your system. You should try a decay time of 50 or 100 fs and see if
that stabilizes your pressure.
Peter

Morad Alawneh wrote:
> *No, I do not.
> I have attached the NAMD script file to those who want to check my the
> input parameters.
>
> Thanks
>
> Morad Alawneh
> *
> Marcos Sotomayor wrote:
>>
>> Are you using multiple time stepping?
>>
>> Marcos
>>
>> On Wed, 28 Mar 2007, Morad Alawneh wrote:
>>
>>> *There is one suggestion, which is to use NVT ensemble, by having
>>> different box sizes and choose the box size that gives you the right
>>> pressure. This suggestion is fine if you can use NVT ensemble, which
>>> varies the pressure and keeps the box size constant.
>>>
>>> Unfortunately, this suggestion is not suitable for my case, since I
>>> have
>>> to use NPAT ensemble.
>>>
>>> I think there are some bugs in the NAMD algorithms, since I found other
>>> two problems. I was not able to constrain the com of my protein around
>>> the origin for whatever the force constant I use in both SMD or Free
>>> Energy of Conformational Change Calculations.
>>>
>>> I sent a copy of my notes to James Phillips, NAMD developer, hopping to
>>> here back from him soon.
>>>
>>> Thanks
>>>
>>> Morad Alawneh
>>> **
>>> *
>>> tamal_at_iitk.ac.in wrote:
>>>> I am also finding the same problem.It never relates with our input
>>>> pressure.Have you got any reply?.
>>>> Thanks for your referred paper in your mail.It is really a
>>>> fundamental paper.
>>>>
>>>> thanks
>>>> Tamal Banerjee
>>>>
>>>>
>>>>> *Dear NAMD Developers and users,
>>>>>
>>>>> I am trying to evaluate the surface tension of my system (a protein
>>>>> surrounded by a lipid bilayer in water and salt) using this equation:
>>>>> 0.5*Lz*(Pzz - Pt), where
>>>>>
>>>>> Lz: the box length in z direction, which is the reaction coordinate.
>>>>> Pzz: the pressure sensor in z direction
>>>>> pt = 0.5*(Pxx + Pyy)
>>>>>
>>>>> I am using the NPAT ensemble with the following settings in NAMD2.6:
>>>>>
>>>>> # Langevin Dynamics Parameters
>>>>> langevin on
>>>>> langevinTemp 303.15
>>>>> langevinDamping 5.0
>>>>> langevinHydrogen on
>>>>> # Pressure Control
>>>>> useGroupPressure yes
>>>>> useFlexibleCell yes
>>>>> useConstantRatio no
>>>>> useConstantArea yes
>>>>> # Nose-Hoover Langevin Piston Pressure Control
>>>>> LangevinPiston on
>>>>> LangevinPistonTarget 1.01325 ;# 1 atm 1.01325 bar
>>>>> LangevinPistonPeriod 200.0
>>>>> LangevinPistonDecay 500.0
>>>>> LangevinPistonTemp 303.15
>>>>> # Periodic Boundary Conditions
>>>>> cellBasisVector1 40.0 00.0 00.0
>>>>> cellBasisVector2 00.0 40.0 00.0
>>>>> cellBasisVector3 00.0 00.0 90.0
>>>>> cellOrigin 0.0 0.0 0.0
>>>>> wrapWater on
>>>>> wrapAll on
>>>>> wrapNearest off
>>>>>
>>>>> According to the article
>>>>> ** ZHANG YH, FELLER SE, BROOKS BR, et al.
>>>>> COMPUTER-SIMULATION OF LIQUID/LIQUID INTERFACES .1. THEORY AND
>>>>> APPLICATION TO OCTANE/WATER
>>>>> <http://apps.isiknowledge.com/WoS/CIW.cgi?SID=2EPNc3mJfLCGDA6hnmm&Func=Abstract&doc=9/5>
>>>>>
>>>>> JOURNAL OF CHEMICAL PHYSICS 103 (23): 10252-10266 DEC 15 1995 *
>>>>> *Pzz should take the value of Pn, which is ~ 1 bar, but when I tested
>>>>> its values it had different values.
>>>>> If I take the average of those pressure values, they are far away
>>>>> from 1
>>>>> bar, here is an example:
>>>>>
>>>>> Pxx Pyy Pzz Pt
>>>>> P <P>
>>>>> 720.2230 415.1900 -427.7030 567.7070 235.9037
>>>>> 235.9037
>>>>> 196.3230 577.8960 266.0390 387.1090 346.7525
>>>>> 107.5640
>>>>> ............... ............... ...............
>>>>> ................ ............... .............
>>>>> **simulation average: -16.4359 -112.3742 -80.3948
>>>>> -54.0773
>>>>>
>>>>> I know there is a huge fluctuations in pressure, but the average
>>>>> values
>>>>> should take care of that and be close to 1 bar for Pzz, P, or <P>.
>>>>>
>>>>> Would you explain to me why do I have this kind of pressure
>>>>> discrepancy
>>>>> for Pzz, P, or <P>?
>>>>>
>>>>> Thanks
>>>>>
>>>>> Morad Alawneh
>>>>>
>>>>>
>>>>> **
>>>>>
>>>>>
>>>>>
>>>>> *
>>>>>
>>>>>
>>>>>
>>>>
>>>>
>>>> *****************************************
>>>> Tamal Banerjee
>>>> Research Scholar
>>>> Department of Chemical Engineering
>>>> I.I.T Kanpur
>>>> Phone: +91-512-2597533 (Lab)
>>>> +91-512-2595433 (Residence)
>>>> URL: home.iitk.ac.in/~tamal
>>>> *****************************************
>>>>
>>>>
>>>
>>
> ------------------------------------------------------------------------
>
> # Equilibrate the system under NPAT
>
> set InputName sys
> set Process0 equil10
> set Process prod
> set LastRun 0
> set OutputName [concat "$InputName\_$Process"]
> #set FixedAtoms fixed.pdb
> #set Restraints [concat "restraints_$Process0\.pdb"]
> set ESF [concat "../prep/$InputName\_$Process0\.xsc"]
> set Numer_of_Steps 5000000
> set FirstTimeStep [expr {$Numer_of_Steps*$LastRun}]period
> set Temp 303.15
> set Dist_Temp 303.15
> #set Gamma ??.0 ;# dyn/cm
> ####################################################################################
> # Input Files
>
> structure "../psfgen/$InputName.psf"
> coordinates "../psfgen/$InputName.pdb"
> bincoordinates [concat "../prep/$InputName\_$Process0\.coor"]
> binvelocities [concat "../prep/$InputName\_$Process0\.vel"];# use this or
> ;# temperature
> parameters "../psfgen/par_all27_prot_lipid.prm"
> paraTypeCharmm on
> ####################################################################################
> # Output Files
>
> outputname $OutputName
> binaryoutput yes ;# default
> #restartname res_$OutputName ;# produces pdb restart file
> #restartfreq 100 ;# generate restart file every N timesteps
> DCDfile $OutputName.dcd
> DCDfreq 500
> ####################################################################################
> # Standard Output
>
> outputEnergies 100 ;# timesteps between energy output
> outputMomenta 100 ;# timesteps between momenta output
> outputPressure 100 ;# timesteps between pressure output
> outputTiming 100 ;# timesteps between timing output
> ####################################################################################
> # Timestep Parameters
>
> numsteps $Numer_of_Steps ;# run stops when this step is reached
> timestep 2.0 ;# fs
> firsttimestep $FirstTimeStep ;# use this instead of restart from 0
> stepspercycle 20 ;# default
> ####################################################################################
> # Simulation Space Partitioning
>
> switching on
> switchDist 8.0
> cutoff 12.0
> pairlistdist 14.0
> margin 2.0
> pairlistsPerCycle 4
> ####################################################################################
> # Basic Dynamics
>
> #temperature $Temp ;# use this or binvelocities
> exclude scaled1-4 ;# specified by CHARMM
> 1-4scaling 1.0 ;# specified by CHARMM
> seed 2006
> rigidBonds all ;# activate ShakeH
> rigidTolerance 0.00000001 ;# default
> rigidIterations 100 ;# default
> useSettle on ;# default
> dielectric 1.0 ;# default
> COMmotion no ;# default
> zeroMomentum yes
> ####################################################################################
> # PME Parameters
>
> PME yes
> PMEGridSizeX 64 ;# Should be > 1.5 the box width
> PMEGridSizeY 64 ;# Should be > 1.5 the box height
> PMEGridSizeZ 256 ;# Should be > 1.5 the box length
> PMEInterpOrder 6
> ####################################################################################
> # Multiple Timestep Parameters
>
> fullElectFrequency 2
> nonbondedFreq 1
> MTSAlgorithm impulse ;# default
> longSplitting c1 ;# default
> ####################################################################################
> # Harmonic Constraint Parameters
>
> #constraints on
> #consexp 2 ;# default, use only even integers
> #consref $Restraints
> #conskfile $Restraints
> #conskcol B
> #selectConstraints on ;# Planer restraint
> #selectConstrZ on
> ####################################################################################
> # Fixed Atoms Parameters
>
> #fixedAtoms on
> #fixedAtomsForces on ;# keep it off at CP
> #fixedAtomsFile $FixedAtoms
> #fixedAtomsCol B
> ####################################################################################
> # Temperature Control
>
> # Langevin Dynamics Parameters
> langevin on
> langevinTemp $Temp
> langevinDamping 5.0 ;# /ps
> langevinHydrogen on ;# default
>
> # Temperature Coupling Parameters
> #tCouple on
> #tCoupleTemp $Temp
> #tCoupleFile
> #tCoupleCol B
>
> # Temperature Rescalling Parameters
> #rescaleFreq 50
> #rescaleTemp $Temp
>
> # Temperature Reassignment Parameters, usefull for annealing
> #reassignFreq 50
> #reassignTemp $Temp
> #reassignIncr 5.0
> #reassignHold $Dist_Temp
> ####################################################################################
> # Boundary Conditions
>
> # Spherical Harmonic Conditions
>
> # Cylindrical Harmonic Conditions
>
> # Periodic Boundary Conditions
> cellBasisVector1 40.0 00.0 00.0
> cellBasisVector2 00.0 40.0 00.0
> cellBasisVector3 00.0 00.0 90.0
> cellOrigin 0.0 0.0 0.0 ;# default
> wrapWater on
> wrapAll on
> wrapNearest off ;# use for non-rectangular cells
> extendedSystem $ESF
> ####################################################################################
> # Pressure Control
>
> useGroupPressure yes
> useFlexibleCell yes ;# anisotropic cell fluctuations
> useConstantRatio no ;# constant shape in x-y plane
> useConstantArea yes ;# constant area; constant dimentions in x-y
> ;#plane
>
> # Nose-Hoover Langevin Piston Pressure Control
> LangevinPiston on
> LangevinPistonTarget 1.01325 ;# 1 atm 1.01325 bar
> LangevinPistonPeriod 200.0 ;# fs
> LangevinPistonDecay 500.0 ;# fs
> LangevinPistonTemp $Temp
> #SurfaceTensionTarget $Gamma
> ####################################################################################
> # Interactive Molecular Dynamics
>
> IMDon on
> IMDport 2006
> IMDfreq 1
> IMDignore yes
> ####################################################################################
> # Energy Minimization: Conjugate Graidient Parameters
>
> #minimization on
> #minTinyStep 0.000001 ;# defalut
> #minBabyStep 0.0100 ;# default
> ####################################################################################
> # # Tcl Forces and Analysis
> #
> # set PDB [concat "../psfgen/$InputName\.pdb"]
> # ###----------------
> # set SelSegID1 {GA1}
> # set SelSegID2 {GA2}
> # set SelSegID3 {}
> # set SelSegID4 {}
> # ###----------------
> # set Protein_Center [list 0.0 0.0 0.0]
> # ###----------------
> # tclForces on
> # tclForcesScript "../psfgen/restraints.tcl"
> ####################################################################################
>
> #minimize 5000
>
> run $Numer_of_Steps

This archive was generated by hypermail 2.1.6 : Wed Feb 29 2012 - 15:44:31 CST