[gmx-developers] pressure calculation in GROMACS CVS
mzink at gwdg.de
mzink at gwdg.de
Wed May 14 19:30:21 CEST 2008
Hi Berk,
here are the results from my complete simulation system. The simulations
were performed with CVS version from 2007-04-07. When preprocessing with
this version, the pressure is much too low and the volume changes
significantly during pressure scaling, resulting in a compression of my
protein (gyration radius declines by 0.2nm in the first steps, berendsen
tau_p=1ps):
Time(ps) Volume
0.000000 33352.519531
1.000000 32016.798828
2.000000 31509.839844
3.000000 31247.882812
4.000000 31129.593750
5.000000 31069.384766
Time(ps) Pressure
0.000000 -1977.399414
1.000000 -507.769196
2.000000 -282.439423
3.000000 -90.428642
4.000000 -78.975014
5.000000 -46.453594
When I did grompp with the cvs version from 2007-06-21 or gromacs-331, the
pressure was about less then 100bar at step 0 and the box volume did not
change.
Although the starting structure are identical, mdrun is always the same
and I checked tpr files with gmxcheck and gmxdump, I cannot find the
mistake. At the moment I have to perform the simulations with constant
volume or with cvs versions BEFORE july 2007, otherwise the box
compression destroys my system. Otherwise I have to run grompp with older
versions, than it is working with current cvs-mdrun. For me it seems so
that grompp is the problem and not mdrun but I don't have a clue what's
going on. Maybe you have an idea?
Regards,
Mareike
> Hi,
>
> The results you posted show, according to me, no significant difference
> in pressure and volume,
> except for the pressure at step 0.
>
> I assume, but you did not state this clearly, that both runs were done
> with exactly the same
> mdrun binary. I changed something with the constraining of the initial
> velocities in mdrun
> some time ago, which could explain the pressure difference at step 0.
> If it is really only a grompp version difference and not mdrun, I do not
> understand how this
> is possible.
>
> You can compare tpr files with:
> gmxcheck -s1 ... -s2 ...
>
> Berk.
>
>
> Mareike Zink wrote:
>> Dear GROMACS developers,
>>
>> I may found a bug in the pressure calculation of the current GROMACS CVS
>> version:
>>
>> My system (protein in water, approx. 320,000 particles) is equilibrated
>> in an
>> NPT ensemble and I want to run it with CVS version from 2007-07-04 or
>> the
>> current version from 2008-04-15 with a constant pressure of 1bar. The
>> following pressure and box volumes are given as an output:
>>
>> Time(ps) Pressure(bar)
>> 0.000000 1519.933594
>> 0.200000 1.404368
>> 0.400000 53.713745
>> 0.600000 -16.645607
>> 0.800000 -40.373356
>> 1.000000 64.932564
>>
>> Time(ps) Volume(nm³)
>> 0.000000 2458.407227
>> 0.200000 2458.920166
>> 0.400000 2458.617188
>> 0.600000 2459.025146
>> 0.800000 2459.828613
>> 1.000000 2459.592285
>>
>> When I did the preprocessing with a CVS version from 2007-06-21 or
>> GROMACS-3.3.1 but run the simulation with the CVS version from
>> 2007-07-04 or
>> 2008-04-15, the pressure and volume look much different:
>>
>> Time(ps) Pressure(bar)
>> 0.000000 95.616234
>> 0.200000 4.654475
>> 0.400000 59.303020
>> 0.600000 -7.656351
>> 0.800000 -1.787814
>> 1.000000 -9.527036
>>
>> Time(ps) Volume
>> 0.000000 2458.407227
>> 0.200000 2458.871582
>> 0.400000 2458.587158
>> 0.600000 2458.832031
>> 0.800000 2459.606445
>> 1.000000 2459.485107
>>
>>
>> Maybe there is a problem with the preprocessor. The problem with the
>> much too
>> high pressure only occures if grompp was done with a CVS version from
>> 2007-07-04 or later.
>>
>> Furthermore, the system I was using here, is just a subsystem of my
>> complete
>> simulation system which contains 4,5 Million atoms. Here I have the same
>> problem but the pressure is much too small (-2000bar). Due to the wrong
>> pressure, the pressure coupling rescales the box size within 1ps (my
>> tau_p)
>> and compresses the box to a much too small volume. Maybe the pressure
>> calculation is buggy?
>> Do you kno anything about this problem?
>>
>> My mdp-files for grompp with gromacs331 and cvs are below.
>>
>> If you need more information or files, please contact me.
>>
>> Thanx,
>>
>> Mareike
>>
>>
>> GROMACS 331:
>>
>> ;
>> ; File 'mdout.mdp' was generated
>> ; By user: mzink (37571)
>> ; On host: beany
>> ; At date: Mon Dec 18 12:45:54 2006
>> ;
>>
>> ; VARIOUS PREPROCESSING OPTIONS
>> title =
>> ; Preprocessor - specify a full path if necessary.
>> cpp = /lib/cpp
>> include =
>> define =
>>
>> ; RUN CONTROL PARAMETERS
>> integrator = md
>> ; Start time and timestep in ps
>> tinit = 0
>> dt = 0.002
>> nsteps = 5000
>> ; For exact run continuation or redoing part of a run
>> init_step = 0
>> ; mode for center of mass motion removal
>> comm-mode = Linear
>> ; number of steps for center of mass motion removal
>> nstcomm = 1
>> ; group(s) for center of mass motion removal
>> comm-grps =
>>
>> ; LANGEVIN DYNAMICS OPTIONS
>> ; Friction coefficient (amu/ps) and random seed
>> bd-fric = 0
>> ld_seed = 1993
>>
>> ; ENERGY MINIMIZATION OPTIONS
>> ; Force tolerance and initial step-size
>> emtol = 0.000001
>> emstep = 0.01
>> ; Max number of iterations in relax_shells
>> niter = 100
>> ; Step size (ps^2) for minimization of flexible constraints
>> fcstep = 0
>> ; Frequency of steepest descents steps when doing CG
>> nstcgsteep = 1000
>> nbfgscorr = 10
>>
>> ; OUTPUT CONTROL OPTIONS
>> ; Output frequency for coords (x), velocities (v) and forces (f)
>> nstxout = 5000
>> nstvout = 5000
>> nstfout = 5000
>> ; Checkpointing helps you continue after crashes
>> nstcheckpoint = 5000
>> ; Output frequency for energies to log file and energy file
>> nstlog = 100
>> nstenergy = 100
>> ; Output frequency and precision for xtc file
>> nstxtcout = 500
>> xtc_precision = 1000
>> ; This selects the subset of atoms for the xtc file. You can
>> ; select multiple groups. By default all atoms will be written.
>> xtc-grps = Protein CA
>> ; Selection of energy groups
>> energygrps = Protein CA SOL
>>
>> ; NEIGHBORSEARCHING PARAMETERS
>> ; nblist update frequency
>> nstlist = 10
>> ; ns algorithm (simple or grid)
>> ns_type = grid
>> ; Periodic boundary conditions: xyz (default), no (vacuum)
>> ; or full (infinite systems only)
>> pbc = xyz
>> ; nblist cut-off
>> rlist = 0.9
>> domain-decomposition = no
>>
>> ; OPTIONS FOR ELECTROSTATICS AND VDW
>> ; Method for doing electrostatics
>> coulombtype = PME
>> rcoulomb_switch = 0
>> rcoulomb = 0.9
>> ; Relative dielectric constant for the medium and the reaction field
>> epsilon_r = 1
>> epsilon_rf = 1
>> ; Method for doing Van der Waals
>> vdw-type = Cut-off
>> ; cut-off lengths
>> rvdw_switch = 0
>> rvdw = 0.9
>> ; Apply long range dispersion corrections for Energy and Pressure
>> DispCorr = No
>> ; Extension of the potential lookup tables beyond the cut-off
>> table-extension = 1
>> ; Seperate tables between energy group pairs
>> energygrp_table =
>> ; Spacing for the PME/PPPM FFT grid
>> fourierspacing = 0.12
>> ; FFT grid size, when a value is 0 fourierspacing will be used
>> fourier_nx = 0
>> fourier_ny = 0
>> fourier_nz = 0
>> ; EWALD/PME/PPPM parameters
>> pme_order = 4
>> ewald_rtol = 1e-05
>> ewald_geometry = 3d
>> epsilon_surface = 0
>> optimize_fft = no
>>
>> ; GENERALIZED BORN ELECTROSTATICS
>> ; Algorithm for calculating Born radii
>> gb_algorithm = Still
>> ; Frequency of calculating the Born radii inside rlist
>> nstgbradii = 1
>> ; Cutoff for Born radii calculation; the contribution from atoms
>> ; between rlist and rgbradii is updated every nstlist steps
>> rgbradii = 2
>> ; Salt concentration in M for Generalized Born models
>> gb_saltconc = 0
>>
>> ; IMPLICIT SOLVENT (for use with Generalized Born electrostatics)
>> implicit_solvent = No
>>
>> ; OPTIONS FOR WEAK COUPLING ALGORITHMS
>> ; Temperature coupling
>> Tcoupl = berendsen
>> ; Groups to couple separately
>> tc_grps = Protein Other
>> ; Time constant (ps) and reference temperature (K)
>> tau_t = 0.1 0.1
>> ref_t = 300 300
>> ; Pressure coupling
>> Pcoupl = berendsen
>> Pcoupltype = isotropic
>> ; Time constant (ps), compressibility (1/bar) and reference P (bar)
>> tau_p = 1.0
>> compressibility = 4.5e-5
>> ref_p = 1.0
>> ; Random seed for Andersen thermostat
>> andersen_seed = 815131
>>
>> ; OPTIONS FOR QMMM calculations
>> QMMM = no
>> ; Groups treated Quantum Mechanically
>> QMMM-grps =
>> ; QM method
>> QMmethod =
>> ; QMMM scheme
>> QMMMscheme = normal
>> ; QM basisset
>> QMbasis =
>> ; QM charge
>> QMcharge =
>> ; QM multiplicity
>> QMmult =
>> ; Surface Hopping
>> SH =
>> ; CAS space options
>> CASorbitals =
>> CASelectrons =
>> SAon =
>> SAoff =
>> SAsteps =
>> ; Scale factor for MM charges
>> MMChargeScaleFactor = 1
>> ; Optimization of QM subsystem
>> bOPT =
>> bTS =
>>
>> ; SIMULATED ANNEALING
>> ; Type of annealing for each temperature group (no/single/periodic)
>> annealing = no no
>> ; Number of time points to use for specifying annealing in each group
>> annealing_npoints =
>> ; List of times at the annealing points for each group
>> annealing_time =
>> ; Temp. at each annealing point, for each group.
>> annealing_temp =
>>
>> ; GENERATE VELOCITIES FOR STARTUP RUN
>> gen_vel = yes
>> gen_temp = 300
>> gen_seed = 173529
>>
>> ; OPTIONS FOR BONDS
>> constraints = all-bonds
>> ; Type of constraint algorithm
>> constraint-algorithm = Lincs
>> ; Do not constrain the start configuration
>> unconstrained_start = no
>> ; Use successive overrelaxation to reduce the number of shake iterations
>> Shake-SOR = no
>> ; Relative tolerance of shake
>> shake_tol = 1e-04
>> ; Highest order in the expansion of the constraint coupling matrix
>> lincs-order = 4
>> ; Number of iterations in the final step of LINCS. 1 is fine for
>> ; normal simulations, but use 2 to conserve energy in NVE runs.
>> ; For energy minimization with constraints it should be 4 to 8.
>> lincs-iter = 1
>> ; Lincs will write a warning to the stderr if in one step a bond
>> ; rotates over more degrees than
>> lincs-warnangle = 30
>> ; Convert harmonic bonds to morse potentials
>> morse = no
>>
>> ; ENERGY GROUP EXCLUSIONS
>> ; Pairs of energy groups for which all non-bonded interactions are
>> excluded
>> energygrp_excl =
>>
>> ; NMR refinement stuff
>> ; Distance restraints type: No, Simple or Ensemble
>> disre = No
>> ; Force weighting of pairs in one distance restraint: Conservative or
>> Equal
>> disre_weighting = Equal
>> ; Use sqrt of the time averaged times the instantaneous violation
>> disre_mixed = no
>> disre_fc = 1000
>> disre_tau = 1.25
>> ; Output frequency for pair distances to energy file
>> nstdisreout = 100
>> ; Orientation restraints: No or Yes
>> orire = no
>> ; Orientation restraints force constant and tau for time averaging
>> orire-fc = 0
>> orire-tau = 0
>> orire-fitgrp =
>> ; Output frequency for trace(SD) and S to energy file
>> nstorireout = 100
>> ; Dihedral angle restraints: No, Simple or Ensemble
>> dihre = No
>> dihre-fc = 1000
>> dihre-tau = 0
>> ; Output frequency for dihedral values to energy file
>> nstdihreout = 100
>>
>> ; Free energy control stuff
>> free_energy = no
>> init_lambda = 0
>> delta_lambda = 0
>> sc-alpha = 0
>> sc-power = 0
>> sc-sigma = 0.3
>>
>> ; Non-equilibrium MD stuff
>> acc-grps =
>> accelerate =
>> freezegrps =
>> freezedim =
>> cos-acceleration = 0
>> deform =
>>
>> ; Electric fields
>> ; Format is number of terms (int) and for all terms an amplitude (real)
>> ; and a phase angle (real)
>> E-x =
>> E-xt =
>> E-y =
>> E-yt =
>> E-z =
>> E-zt =
>>
>> ; User defined thingies
>> user1-grps =
>> user2-grps =
>> userint1 = 0
>> userint2 = 0
>> userint3 = 0
>> userint4 = 0
>> userreal1 = 0
>> userreal2 = 0
>> userreal3 = 0
>> userreal4 = 0
>>
>>
>> GROMACS CVS from 2007-04-07
>> ;
>> ; File 'mdout.mdp' was generated
>> ; By user: mzink (37571)
>> ; On host: beany
>> ; At date: Mon Mar 17 20:14:43 2008
>> ;
>>
>> ; VARIOUS PREPROCESSING OPTIONS
>> title =
>> ; Preprocessor - specify a full path if necessary.
>> cpp = /lib/cpp
>> include =
>> define =
>>
>> ; RUN CONTROL PARAMETERS
>> integrator = md
>> ; Start time and timestep in ps
>> tinit = 0
>> dt = 0.002
>> nsteps = 5000
>> ; For exact run continuation or redoing part of a run
>> init_step = 0
>> ; mode for center of mass motion removal
>> comm-mode = Linear
>> ; number of steps for center of mass motion removal
>> nstcomm = 1
>> ; group(s) for center of mass motion removal
>> comm-grps =
>>
>> ; LANGEVIN DYNAMICS OPTIONS
>> ; Friction coefficient (amu/ps) and random seed
>> bd-fric = 0
>> ld_seed = 1993
>>
>> ; ENERGY MINIMIZATION OPTIONS
>> ; Force tolerance and initial step-size
>> emtol = 0.000001
>> emstep = 0.01
>> ; Max number of iterations in relax_shells
>> niter = 100
>> ; Step size (ps^2) for minimization of flexible constraints
>> fcstep = 0
>> ; Frequency of steepest descents steps when doing CG
>> nstcgsteep = 1000
>> nbfgscorr = 10
>>
>> ; TEST PARTICLE INSERTION OPTIONS
>> rtpi = 0.05
>>
>> ; OUTPUT CONTROL OPTIONS
>> ; Output frequency for coords (x), velocities (v) and forces (f)
>> nstxout = 5000
>> nstvout = 5000
>> nstfout = 5000
>> ; Checkpointing helps you continue after crashes
>> nstcheckpoint = 5000
>> ; Output frequency for energies to log file and energy file
>> nstlog = 100
>> nstenergy = 100
>> ; Output frequency and precision for xtc file
>> nstxtcout = 500
>> xtc_precision = 1000
>> ; This selects the subset of atoms for the xtc file. You can
>> ; select multiple groups. By default all atoms will be written.
>> xtc-grps = Protein CA
>> ; Selection of energy groups
>> energygrps = Protein CA SOL
>>
>> ; NEIGHBORSEARCHING PARAMETERS
>> ; nblist update frequency
>> nstlist = 10
>> ; ns algorithm (simple or grid)
>> ns_type = grid
>> ; Periodic boundary conditions: xyz, no, xy
>> pbc = xyz
>> periodic_molecules = no
>> ; nblist cut-off
>> rlist = 0.9
>>
>> ; OPTIONS FOR ELECTROSTATICS AND VDW
>> ; Method for doing electrostatics
>> coulombtype = PME
>> rcoulomb_switch = 0
>> rcoulomb = 0.9
>> ; Relative dielectric constant for the medium and the reaction field
>> epsilon_r = 1
>> epsilon_rf = 1
>> ; Method for doing Van der Waals
>> vdw-type = Cut-off
>> ; cut-off lengths
>> rvdw_switch = 0
>> rvdw = 0.9
>> ; Apply long range dispersion corrections for Energy and Pressure
>> DispCorr = No
>> ; Extension of the potential lookup tables beyond the cut-off
>> table-extension = 1
>> ; Seperate tables between energy group pairs
>> energygrp_table =
>> ; Spacing for the PME/PPPM FFT grid
>> ;fourierspacing = 0.12
>> ; FFT grid size, when a value is 0 fourierspacing will be used
>> fourier_nx = 120
>> fourier_ny = 120
>> fourier_nz = 120
>> ; EWALD/PME/PPPM parameters
>> pme_order = 4
>> ewald_rtol = 1e-05
>> ewald_geometry = 3d
>> epsilon_surface = 0
>> optimize_fft = no
>>
>> ; GENERALIZED BORN ELECTROSTATICS
>> ; Algorithm for calculating Born radii
>> gb_algorithm = Still
>> ; Frequency of calculating the Born radii inside rlist
>> nstgbradii = 1
>> ; Cutoff for Born radii calculation; the contribution from atoms
>> ; between rlist and rgbradii is updated every nstlist steps
>> rgbradii = 2
>> ; Salt concentration in M for Generalized Born models
>> gb_saltconc = 0
>>
>> ; IMPLICIT SOLVENT (for use with Generalized Born electrostatics)
>> implicit_solvent = No
>>
>> ; OPTIONS FOR WEAK COUPLING ALGORITHMS
>> ; Temperature coupling
>> Tcoupl = berendsen
>> ; Groups to couple separately
>> tc_grps = Protein Other
>> ; Time constant (ps) and reference temperature (K)
>> tau_t = 0.1 0.1
>> ref_t = 300 300
>> ; Pressure coupling
>> Pcoupl = berendsen
>> Pcoupltype = isotropic
>> ; Time constant (ps), compressibility (1/bar) and reference P (bar)
>> tau_p = 1.0
>> compressibility = 4.5e-5
>> ref_p = 1.0
>> ; Scaling of reference coordinates, No, All or COM
>> refcoord_scaling = No
>> ; Random seed for Andersen thermostat
>> andersen_seed = 815131
>>
>> ; OPTIONS FOR QMMM calculations
>> QMMM = no
>> ; Groups treated Quantum Mechanically
>> QMMM-grps =
>> ; QM method
>> QMmethod =
>> ; QMMM scheme
>> QMMMscheme = normal
>> ; QM basisset
>> QMbasis =
>> ; QM charge
>> QMcharge =
>> ; QM multiplicity
>> QMmult =
>> ; Surface Hopping
>> SH =
>> ; CAS space options
>> CASorbitals =
>> CASelectrons =
>> SAon =
>> SAoff =
>> SAsteps =
>> ; Scale factor for MM charges
>> MMChargeScaleFactor = 1
>> ; Optimization of QM subsystem
>> bOPT =
>> bTS =
>>
>> ; SIMULATED ANNEALING
>> ; Type of annealing for each temperature group (no/single/periodic)
>> annealing = no no
>> ; Number of time points to use for specifying annealing in each group
>> annealing_npoints =
>> ; List of times at the annealing points for each group
>> annealing_time =
>> ; Temp. at each annealing point, for each group.
>> annealing_temp =
>>
>> ; GENERATE VELOCITIES FOR STARTUP RUN
>> gen_vel = yes
>> gen_temp = 300
>> gen_seed = 173529
>>
>> ; OPTIONS FOR BONDS
>> constraints = all-bonds
>> ; Type of constraint algorithm
>> constraint-algorithm = Lincs
>> ; Do not constrain the start configuration
>> continuation = yes
>> ; Use successive overrelaxation to reduce the number of shake iterations
>> Shake-SOR = no
>> ; Relative tolerance of shake
>> shake_tol = 1e-04
>> ; Highest order in the expansion of the constraint coupling matrix
>> lincs-order = 4
>> ; Number of iterations in the final step of LINCS. 1 is fine for
>> ; normal simulations, but use 2 to conserve energy in NVE runs.
>> ; For energy minimization with constraints it should be 4 to 8.
>> lincs-iter = 1
>> ; Lincs will write a warning to the stderr if in one step a bond
>> ; rotates over more degrees than
>> lincs-warnangle = 30
>> ; Convert harmonic bonds to morse potentials
>> morse = no
>>
>> ; ENERGY GROUP EXCLUSIONS
>> ; Pairs of energy groups for which all non-bonded interactions are
>> excluded
>> energygrp_excl =
>>
>> ; WALLS
>> ; Number of walls, type, atom types, densities and box-z scale factor
>> for
>> Ewald
>> nwall = 0
>> wall_type = 9-3
>> wall_atomtype =
>> wall_density =
>> wall_ewald_zfac = 3
>>
>> ; COM PULLING
>> ; Pull type: no, umbrella, constraint, constant_force
>> pull = no
>>
>> ; NMR refinement stuff
>> ; Distance restraints type: No, Simple or Ensemble
>> disre = No
>> ; Force weighting of pairs in one distance restraint: Conservative or
>> Equal
>> disre_weighting = Equal
>> ; Use sqrt of the time averaged times the instantaneous violation
>> disre_mixed = no
>> disre_fc = 1000
>> disre_tau = 1.25
>> ; Output frequency for pair distances to energy file
>> nstdisreout = 100
>> ; Orientation restraints: No or Yes
>> orire = no
>> ; Orientation restraints force constant and tau for time averaging
>> orire-fc = 0
>> orire-tau = 0
>> orire-fitgrp =
>> ; Output frequency for trace(SD) and S to energy file
>> nstorireout = 100
>> ; Dihedral angle restraints: No, Simple or Ensemble
>> dihre = No
>> dihre-fc = 1000
>> dihre-tau = 0
>> ; Output frequency for dihedral values to energy file
>> nstdihreout = 100
>>
>> ; Free energy control stuff
>> free_energy = no
>> init_lambda = 0
>> delta_lambda = 0
>> sc-alpha = 0
>> sc-power = 0
>> sc-sigma = 0.3
>>
>> ; Non-equilibrium MD stuff
>> acc-grps =
>> accelerate =
>> freezegrps =
>> freezedim =
>> cos-acceleration = 0
>> deform =
>>
>> ; Electric fields
>> ; Format is number of terms (int) and for all terms an amplitude (real)
>> ; and a phase angle (real)
>> E-x =
>> E-xt =
>> E-y =
>> E-yt =
>> E-z =
>> E-zt =
>>
>> ; User defined thingies
>> user1-grps =
>> user2-grps =
>> userint1 = 0
>> userint2 = 0
>> userint3 = 0
>> userint4 = 0
>> userreal1 = 0
>> userreal2 = 0
>> userreal3 = 0
>> userreal4 = 0
>>
>>
>>
>> _______________________________________________
>> gmx-developers mailing list
>> gmx-developers at gromacs.org
>> http://www.gromacs.org/mailman/listinfo/gmx-developers
>> Please don't post (un)subscribe requests to the list. Use the
>> www interface or send it to gmx-developers-request at gromacs.org.
>>
>> This email was Anti Virus checked by Astaro Security Gateway.
>> http://www.astaro.com
>>
>>
>
>
More information about the gromacs.org_gmx-developers
mailing list