[gmx-users] TPI Results differ in v4.5.7 and v4.6.1

João M. Damas jmdamas at itqb.unl.pt
Mon Sep 16 19:27:50 CEST 2013


I am sorry for the late follow-up on this subject.

My results using PME electrostatics do not show any differences between
versions 4.5.4 and 4.6.3 (I also checked with 4.6.1 and it gave the same
results):

https://www.dropbox.com/s/o2kswrw5eq8fcsp/plot_batch-1_pme.png

I hope that you, Niels, have found a solution to your problem, which
probably has a different origin from what I initially thought...

Best,
João


On Mon, Jul 1, 2013 at 1:29 PM, João M. Damas <jmdamas at itqb.unl.pt> wrote:

> I have run TPI using three versions (4.0.4, 4.5.4 and 4.6.1) and three
> different insertions particles: CH4 (uncharged), Cl- (negative) and Na+
> (positive). All TPI were run on the same simulations of SPC water and the
> three particles taken from GROMOS force-field. Reaction-field was used for
> electrostatics.
>
> https://www.dropbox.com/s/in66zx8t2mrprgt/plot_batch-1.png
>
> As you can see, for an uncharged particle all the versions provide the
> same result. The same does not happen when charged particles are inserted,
> which hint on problems with the electrostatics between 4.0.4 and 4.[56].X.
> Like when Niels reported for the Cut-off scheme, the reaction-field
> provides the same results for 4.5.4 and 4.6.1. Is it indeed a problem with
> PME?
>
> I still have not tested PME, but I think Niels' test foretells the results
> I'm going to get. Niels, can you confirm this issue is only happening with
> charged particles? And are you going to file an issue like Szilárd
> suggested?
>
> Best,
> João
>
>
>
>
>
>
> On Sat, Jun 29, 2013 at 5:21 PM, João M. Damas <jmdamas at itqb.unl.pt>wrote:
>
>> Niels,
>>
>> Which force-field did you use? I guess an uncharged CH4 shouldn't be
>> giving different results for TPI when changing coulomb... Actually, coulomb
>> is turned off if there's no charge in the particles to insert, if I
>> remember the code correctly.
>>
>> João
>>
>>
>> On Mon, Jun 24, 2013 at 3:40 PM, Niels Müller <uni at nielsm.de> wrote:
>>
>>> Hi João,
>>>
>>> Indeed your instinct seems to be good! When switching the Coulomb-Type
>>> to Cut-Off, there doesn't seem to be a difference between 4.6 and 4.5.
>>> Apparently its an issue with the PME sum. We will investigate further.
>>>
>>>
>>> Am 24.06.2013 um 14:42 schrieb João M. Damas <jmdamas at itqb.unl.pt>:
>>>
>>> > Niels,
>>> >
>>> > This is very interesting. At our group, a colleague of mine and I have
>>> also
>>> > identified differences in the TPI integrator between 4.0.X and 4.5.X,
>>> but
>>> > we still haven't had the time to report it properly, since we are
>>> using a
>>> > slightly modified version of the TPI algorithm.
>>> >
>>> > Instinctively, we were attributing it to some different behaviours in
>>> the
>>> > RF that are observed between those versions. We also know that the TPI
>>> > algorithm began allowing PME treatment from 4.5.X onwards, so maybe
>>> there
>>> > are some differences going on the electrostatics level? But, IIRC, no
>>> > modifications to the TPI code were on the release notes from 4.5.X to
>>> > 4.6.X...
>>> >
>>> > We'll try to find some time to report our findings as soon as possible.
>>> > Maybe they are related.
>>> >
>>> > Best,
>>> > João
>>> >
>>> >
>>> > On Mon, Jun 24, 2013 at 10:19 AM, Niels Müller <uni at nielsm.de> wrote:
>>> >
>>> >> Hi GMX Users,
>>> >>
>>> >> We are computing the chemical potential of different gas molecules in
>>> a
>>> >> polymer melt with the tpi integrator.
>>> >> The computations are done for CO2 and CH4.
>>> >> The previous computations were done with v4.5.5 or 4.5.7 and gave
>>> equal
>>> >> results.
>>> >>
>>> >> I recently switched to gromacs version 4.6.1, and the chemical
>>> potential
>>> >> computed by this version is shifted by a nearly constant factor,
>>> which is
>>> >> different for the two gas molecules.
>>> >> We are perplexed what causes this shift. Was there any change in the
>>> new
>>> >> version that affects the tpi integration? I will provide the mdp file
>>> we
>>> >> used below.
>>> >>
>>> >> The tpi integration is run on basis of the last 10 ns of a 30 ns NVT
>>> >> simulation with 'mdrun -rerun'.
>>> >>
>>> >> Best regards,
>>> >> Niels.
>>> >>
>>> >> #########################
>>> >> The mdp file:
>>> >> #########################
>>> >>
>>> >> ; VARIOUS PREPROCESSING OPTIONS
>>> >> cpp                      = cpp
>>> >> include                =
>>> >> define                  =
>>> >>
>>> >> ; RUN CONTROL PARAMETERS
>>> >> integrator               = tpi
>>> >> ; Start time and timestep in ps
>>> >> tinit                    = 0
>>> >> dt                       = 0.001
>>> >> nsteps                   = 1000000
>>> >> ; 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
>>> >> ; Temperature, friction coefficient (amu/ps) and random seed
>>> >> bd-fric                  = 0.5
>>> >> ld-seed                  = 1993
>>> >>
>>> >> ; ENERGY MINIMIZATION OPTIONS
>>> >> ; Force tolerance and initial step-size
>>> >> emtol                    = 100
>>> >> emstep                   = 0.01
>>> >> ; Max number of iterations in relax_shells
>>> >> niter                    = 20
>>> >> ; Step size (1/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                  = 100
>>> >> nstvout                  = 0
>>> >> nstfout                  = 0
>>> >> ; Checkpointing helps you continue after crashes
>>> >> nstcheckpoint            = 100
>>> >> ; Output frequency for energies to log file and energy file
>>> >> nstlog                   = 100
>>> >> nstenergy                = 100
>>> >> ; Output frequency and precision for xtc file
>>> >> nstxtcout                = 0
>>> >> 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                 =
>>> >> ; Selection of energy groups
>>> >> energygrps               =
>>> >>
>>> >> ; NEIGHBORSEARCHING PARAMETERS
>>> >> ; nblist update frequency
>>> >> nstlist                  = 5
>>> >> ; 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
>>> >> ; Dielectric constant (DC) for cut-off or DC of reaction field
>>> >> epsilon-r                = 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                 = EnerPres
>>> >> ; Extension of the potential lookup tables beyond the cut-off
>>> >> table-extension          = 1
>>> >> ; 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                   = V-rescale
>>> >> ; Groups to couple separately
>>> >> tc-grps                  = System
>>> >> ; Time constant (ps) and reference temperature (K)
>>> >> tau_t                    = 0.1
>>> >> ref_t                    = 318
>>> >> ; Pressure coupling
>>> >> Pcoupl                 = Parrinello-Rahman
>>> >> Pcoupltype               = isotropic
>>> >> ; Time constant (ps), compressibility (1/bar) and reference P (bar)
>>> >> tau_p                    = 5.0
>>> >> compressibility          = 4.5e-5
>>> >> ref_p                    = 1.0
>>> >> ; Random seed for Andersen thermostat
>>> >> andersen_seed            = 815131
>>> >>
>>> >> ; SIMULATED ANNEALING
>>> >> ; Type of annealing for each temperature group (no/single/periodic)
>>> >> annealing                = 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                 = 400
>>> >> gen_seed                 = 1993
>>> >>
>>> >> ; OPTIONS FOR BONDS
>>> >> ;constraints              = none
>>> >> 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          = Conservative
>>> >> ; Use sqrt of the time averaged times the instantaneous violation
>>> >> disre-mixed              = no
>>> >> disre-fc                 = 1000
>>> >> disre-tau                = 0
>>> >> ; 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) 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-sigma                 = 0.3
>>> >>
>>> >> ; Non-equilibrium MD stuff
>>> >> acc-grps                 =
>>> >> accelerate               =
>>> >> freezegrps               =
>>> >> freezedim                =
>>> >> cos-acceleration         = 0
>>> >>
>>> >>
>>> >> --
>>> >> gmx-users mailing list    gmx-users at gromacs.org
>>> >> http://lists.gromacs.org/mailman/listinfo/gmx-users
>>> >> * Please search the archive at
>>> >> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
>>> >> * Please don't post (un)subscribe requests to the list. Use the
>>> >> www interface or send it to gmx-users-request at gromacs.org.
>>> >> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>>> >>
>>> >
>>> >
>>> >
>>> > --
>>> > João M. Damas
>>> > PhD Student
>>> > Protein Modelling Group
>>> > ITQB-UNL, Oeiras, Portugal
>>> > Tel:+351-214469613
>>> > --
>>> > gmx-users mailing list    gmx-users at gromacs.org
>>> > http://lists.gromacs.org/mailman/listinfo/gmx-users
>>> > * Please search the archive at
>>> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
>>> > * Please don't post (un)subscribe requests to the list. Use the
>>> > www interface or send it to gmx-users-request at gromacs.org.
>>> > * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>>>
>>>
>>> --
>>> gmx-users mailing list    gmx-users at gromacs.org
>>> http://lists.gromacs.org/mailman/listinfo/gmx-users
>>> * Please search the archive at
>>> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
>>> * Please don't post (un)subscribe requests to the list. Use the
>>> www interface or send it to gmx-users-request at gromacs.org.
>>> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>>>
>>
>>
>>
>> --
>> João M. Damas
>> PhD Student
>> Protein Modelling Group
>> ITQB-UNL, Oeiras, Portugal
>> Tel:+351-214469613
>>
>
>
>
> --
> João M. Damas
> PhD Student
> Protein Modelling Group
> ITQB-UNL, Oeiras, Portugal
> Tel:+351-214469613
>



-- 
João M. Damas
PhD Student
Protein Modelling Group
ITQB-UNL, Oeiras, Portugal
Tel:+351-214469613



More information about the gromacs.org_gmx-users mailing list