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 <jmda...@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 <jmda...@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 <u...@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 <jmda...@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 <u...@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@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-requ...@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@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-requ...@gromacs.org. >>> > * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists >>> >>> >>> -- >>> gmx-users mailing list gmx-users@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-requ...@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 -- gmx-users mailing list gmx-users@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-requ...@gromacs.org. * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
