Dear Gromacs Users, I am trying to simulate a system consisting of a vacuum/condensed phase interface in which a 6x6x12nm condensed phase region is flanked on both ends (in the z-dimension) by a 6x6x12nm vacuum region to form overall box dimensions of 6x6x36 nm. The system is a binary liquid mixture of methanol (0.125 mole fraction methanol) in water using a polarizable (charge on a spring) force field (COS/M methanol and COS/G2 water) at 300K and 1bar. The system is stable in Gromacs 4.5.3; however, mdrun gives a segmentation fault in Gromacs 4.6 when attempting to do dynamics (energy minimization completes with no apparent problems). If I remove the vacuum region, mdrun works. If I incrementally add 2 Angstroms to the z-dimension until I reached a vacuum region of 34 Angstroms total (17 Angstroms on both sides of the condensed phase region) and try to simulate these systems, mdrun works every time. When I reach 36 Angstroms, the segmentation fault re-appears. Although not the system I am actually interested in, I did some simulations using Gromacs 4.6 with the 34 Angstrom vacuum region system and observed an undulating and very turbulant "vacuum"/condensed phase interface in which a column of water/methanol mixture came out of the condensed phase region to connect the two interfaces. Also, the center of mass motion of the system appeared to not have been removed. In contrast, using Gromacs 4.5.3, the interface is not undulating, but is "calm" and qualitatively planar, no column forms to connect the interfaces, and there is no problem with the center of mass motion removal. Some of the molecules do enter the "vacuum" region when running with 4.5.3, but this appears to be due to the movement of individual molecules, not a collective motion of many molecules. This system runs fine with a non-polarizable force field in Gromacs 4.6.
I have also compared several properties (using g_energy) of the bulk system (no vacuum region) using Gromacs 4.6 and 4.5.3 and they are not the same for the polarizable force field, but they are the same for the non-polarizable force field. Specifically, with the polarizable force field, the LJ(SR) energy is more positive with 4.6, the LJ(LR) energy is more negative with 4.6, the Coulomb(SR) energy is more negative with 4.6, the Could. recip. energy is more negative with 4.6, the polarization energy is more positive with 4.6, the potential energy is more negative in 4.6, the average kinetic energy (and temperature) is the same but the fluctuations are greater in 4.6, the total energy is more negative in 4.6, the pressure looks fine, and the volume looks fine. Where I've noted differences, these are all statistically significant differences. I would like to know if I can just use 4.5.3 and assume the differences between the results of 4.5.3 and 4.6 are due to some problem in 4.6. I am using all the same input files and commands with both versions, only different executables. I ran the regression tests for 4.6 when I installed it, and passed them all. Sincerely, Elizabeth ----- ADDITIONAL DETAILS: Below is where the problem appears for the interface system when z-dimension of the vacuum region is >= 36 Angstroms total. eq4 is my first attempt at dynamics, after three successful energy minimizations (1st: charges screened and no bond constraints, 2nd: charges felt but no bond constraints, 3rd: charges felt and bond constraints on) [ploetz@cluster AddRemainingVacuumBack]$ grompp -f eq4.mdp -c em3.pdb -o eq4.tpr -n index.ndx -p sys.top -nice 0 :-) G R O M A C S (-: Green Red Orange Magenta Azure Cyan Skyblue :-) VERSION 4.6 (-: Contributions from Mark Abraham, Emile Apol, Rossen Apostolov, Herman J.C. Berendsen, Aldert van Buuren, Pär Bjelkmar, Rudi van Drunen, Anton Feenstra, Gerrit Groenhof, Christoph Junghans, Peter Kasson, Carsten Kutzner, Per Larsson, Pieter Meulenhoff, Teemu Murtola, Szilard Pall, Sander Pronk, Roland Schulz, Michael Shirts, Alfons Sijbers, Peter Tieleman, Berk Hess, David van der Spoel, and Erik Lindahl. Copyright (c) 1991-2000, University of Groningen, The Netherlands. Copyright (c) 2001-2012,2013, The GROMACS development team at Uppsala University & The Royal Institute of Technology, Sweden. check out http://www.gromacs.org for more information. This program is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation; either version 2.1 of the License, or (at your option) any later version. :-) grompp (-: Option Filename Type Description ------------------------------------------------------------ -f eq4.mdp Input grompp input file with MD parameters -po mdout.mdp Output grompp input file with MD parameters -c em3.pdb Input Structure file: gro g96 pdb tpr etc. -r conf.gro Input, Opt. Structure file: gro g96 pdb tpr etc. -rb conf.gro Input, Opt. Structure file: gro g96 pdb tpr etc. -n index.ndx Input, Opt! Index file -p sys.top Input Topology file -pp processed.top Output, Opt. Topology file -o eq4.tpr Output Run input file: tpr tpb tpa -t traj.trr Input, Opt. Full precision trajectory: trr trj cpt -e ener.edr Input, Opt. Energy file -ref rotref.trr In/Out, Opt. Full precision trajectory: trr trj cpt Option Type Value Description ------------------------------------------------------ -[no]h bool no Print help info and quit -[no]version bool no Print version info and quit -nice int 0 Set the nicelevel -[no]v bool no Be loud and noisy -time real -1 Take frame at or first after this time. -[no]rmvsbds bool yes Remove constant bonded interactions with virtual sites -maxwarn int 0 Number of allowed warnings during input processing. Not for normal use and may generate unstable systems -[no]zero bool no Set parameters for bonded interactions without defaults to zero instead of generating an error -[no]renum bool yes Renumber atomtypes and minimize number of atomtypes Ignoring obsolete mdp entry 'cpp' Ignoring obsolete mdp entry 'domain-decomposition' Ignoring obsolete mdp entry 'andersen_seed' Ignoring obsolete mdp entry 'nstcheckpoint' Replacing old mdp entry 'unconstrained-start' by 'continuation' NOTE 1 [file eq4.mdp]: The Berendsen thermostat does not generate the correct kinetic energy distribution. You might want to consider using the V-rescale thermostat. Generated 80 of the 91 non-bonded parameter combinations Excluding 2 bonded neighbours molecule type 'COSM' turning all bonds into constraints... Excluding 2 bonded neighbours molecule type 'COSG2' turning all bonds into constraints... Velocities were taken from a Maxwell distribution at 300.15 K Cleaning up constraints and constant bonded interactions with virtual sites Number of degrees of freedom in T-Coupling group MOH is 8243.62 Number of degrees of freedom in T-Coupling group SOL is 57717.38 Largest charge group radii for Van der Waals: 0.118, 0.118 nm Largest charge group radii for Coulomb: 0.118, 0.118 nm Calculating fourier grid dimensions for X Y Z Using a fourier grid of 48x48x288, spacing 0.120 0.120 0.120 Estimate for the relative computational load of the PME mesh part: 0.37 This run will generate roughly 622 Mb of data There was 1 note gcq#228: "Bailed Out Of Edge Synchronization After 10,000 Iterations" (X/Motif) [ploetz@cluster AddRemainingVacuumBack]$ more eq4.mdp ; 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 = 500000 ; 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 = 500 ; group(s) for center of mass motion removal comm-grps = ; OUTPUT CONTROL OPTIONS ; Output frequency for coords (x), velocities (v) and forces (f) nstxout = 500 nstvout = 0 nstfout = 0 ; Checkpointing helps you continue after crashes nstcheckpoint = 5000 ; Output frequency for energies to log file and energy file nstlog = 500 nstenergy = 500 ; 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 = MOH 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 = 1.0 domain-decomposition = no ; OPTIONS FOR ELECTROSTATICS AND VDW ; Method for doing electrostatics coulombtype = PME rcoulomb-switch = 0 rcoulomb = 1.0 ; 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 = 1.5 ; 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 ; 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 = MOH SOL ; Time constant (ps) and reference temperature (K) tau-t = 0.1 0.1 ref-t = 300.15 300.15 ; Pressure coupling Pcoupl = no Pcoupltype = isotropic ; Time constant (ps), compressibility (1/bar) and reference P (bar) tau-p = 0.5 compressibility = 4.5e-5 ref-p = 1 ; Random seed for Andersen thermostat andersen_seed = 815131 ; GENERATE VELOCITIES FOR STARTUP RUN gen-vel = yes gen-temp = 300.15 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 = 0.0001 ; 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 = 4 ; 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 = [ploetz@cluster AddRemainingVacuumBack]$ [ploetz@cluster AddRemainingVacuumBack]$ mdrun -v -s eq4.tpr -c eq4.pdb -g eq4.log -o eq4.trr -nice 0 :-) G R O M A C S (-: GRoups of Organic Molecules in ACtion for Science :-) VERSION 4.6 (-: Contributions from Mark Abraham, Emile Apol, Rossen Apostolov, Herman J.C. Berendsen, Aldert van Buuren, Pär Bjelkmar, Rudi van Drunen, Anton Feenstra, Gerrit Groenhof, Christoph Junghans, Peter Kasson, Carsten Kutzner, Per Larsson, Pieter Meulenhoff, Teemu Murtola, Szilard Pall, Sander Pronk, Roland Schulz, Michael Shirts, Alfons Sijbers, Peter Tieleman, Berk Hess, David van der Spoel, and Erik Lindahl. Copyright (c) 1991-2000, University of Groningen, The Netherlands. Copyright (c) 2001-2012,2013, The GROMACS development team at Uppsala University & The Royal Institute of Technology, Sweden. check out http://www.gromacs.org for more information. This program is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation; either version 2.1 of the License, or (at your option) any later version. :-) mdrun (-: Option Filename Type Description ------------------------------------------------------------ -s eq4.tpr Input Run input file: tpr tpb tpa -o eq4.trr Output Full precision trajectory: trr trj cpt -x traj.xtc Output, Opt. Compressed trajectory (portable xdr format) -cpi state.cpt Input, Opt. Checkpoint file -cpo state.cpt Output, Opt. Checkpoint file -c eq4.pdb Output Structure file: gro g96 pdb etc. -e ener.edr Output Energy file -g eq4.log Output Log file -dhdl dhdl.xvg Output, Opt. xvgr/xmgr file -field field.xvg Output, Opt. xvgr/xmgr file -table table.xvg Input, Opt. xvgr/xmgr file -tabletf tabletf.xvg Input, Opt. xvgr/xmgr file -tablep tablep.xvg Input, Opt. xvgr/xmgr file -tableb table.xvg Input, Opt. xvgr/xmgr file -rerun rerun.xtc Input, Opt. Trajectory: xtc trr trj gro g96 pdb cpt -tpi tpi.xvg Output, Opt. xvgr/xmgr file -tpid tpidist.xvg Output, Opt. xvgr/xmgr file -ei sam.edi Input, Opt. ED sampling input -eo edsam.xvg Output, Opt. xvgr/xmgr file -j wham.gct Input, Opt. General coupling stuff -jo bam.gct Output, Opt. General coupling stuff -ffout gct.xvg Output, Opt. xvgr/xmgr file -devout deviatie.xvg Output, Opt. xvgr/xmgr file -runav runaver.xvg Output, Opt. xvgr/xmgr file -px pullx.xvg Output, Opt. xvgr/xmgr file -pf pullf.xvg Output, Opt. xvgr/xmgr file -ro rotation.xvg Output, Opt. xvgr/xmgr file -ra rotangles.log Output, Opt. Log file -rs rotslabs.log Output, Opt. Log file -rt rottorque.log Output, Opt. Log file -mtx nm.mtx Output, Opt. Hessian matrix -dn dipole.ndx Output, Opt. Index file -multidir rundir Input, Opt., Mult. Run directory -membed membed.dat Input, Opt. Generic data file -mp membed.top Input, Opt. Topology file -mn membed.ndx Input, Opt. Index file Option Type Value Description ------------------------------------------------------ -[no]h bool no Print help info and quit -[no]version bool no Print version info and quit -nice int 0 Set the nicelevel -deffnm string Set the default filename for all file options -xvg enum xmgrace xvg plot formatting: xmgrace, xmgr or none -[no]pd bool no Use particle decompostion -dd vector 0 0 0 Domain decomposition grid, 0 is optimize -ddorder enum interleave DD node order: interleave, pp_pme or cartesian -npme int -1 Number of separate nodes to be used for PME, -1 is guess -nt int 0 Total number of threads to start (0 is guess) -ntmpi int 0 Number of thread-MPI threads to start (0 is guess) -ntomp int 0 Number of OpenMP threads per MPI process/thread to start (0 is guess) -ntomp_pme int 0 Number of OpenMP threads per MPI process/thread to start (0 is -ntomp) -pin enum auto Fix threads (or processes) to specific cores: auto, on or off -pinoffset int 0 The starting logical core number for pinning to cores; used to avoid pinning threads from different mdrun instances to the same core -pinstride int 0 Pinning distance in logical cores for threads, use 0 to minimize the number of threads per physical core -gpu_id string List of GPU id's to use -[no]ddcheck bool yes Check for all bonded interactions with DD -rdd real 0 The maximum distance for bonded interactions with DD (nm), 0 is determine from initial coordinates -rcon real 0 Maximum distance for P-LINCS (nm), 0 is estimate -dlb enum auto Dynamic load balancing (with DD): auto, no or yes -dds real 0.8 Minimum allowed dlb scaling of the DD cell size -gcom int -1 Global communication frequency -nb enum auto Calculate non-bonded interactions on: auto, cpu, gpu or gpu_cpu -[no]tunepme bool yes Optimize PME load between PP/PME nodes or GPU/CPU -[no]testverlet bool no Test the Verlet non-bonded scheme -[no]v bool yes Be loud and noisy -[no]compact bool yes Write a compact log file -[no]seppot bool no Write separate V and dVdl terms for each interaction type and node to the log file(s) -pforce real -1 Print all forces larger than this (kJ/mol nm) -[no]reprod bool no Try to avoid optimizations that affect binary reproducibility -cpt real 15 Checkpoint interval (minutes) -[no]cpnum bool no Keep and number checkpoint files -[no]append bool yes Append to previous output files when continuing from checkpoint instead of adding the simulation part number to all file names -nsteps int -2 Run this number of steps, overrides .mdp file option -maxh real -1 Terminate after 0.99 times this time (hours) -multi int 0 Do multiple simulations in parallel -replex int 0 Attempt replica exchange periodically with this period (steps) -nex int 0 Number of random exchanges to carry out each exchange interval (N^3 is one suggestion). -nex zero or not specified gives neighbor replica exchange. -reseed int -1 Seed for replica exchange, -1 is generate a seed -[no]ionize bool no Do a simulation including the effect of an X-Ray bombardment on your system Reading file eq4.tpr, VERSION 4.6 (single precision) Using 12 MPI threads starting mdrun 'COSM in water x_COSM 0.125' 500000 steps, 1000.0 ps. Segmentation fault [ploetz@cluster AddRemainingVacuumBack]$ [ploetz@cluster AddRemainingVacuumBack]$ more sys.top #include "polff.itp" #include "cosm.itp" #include "cosg2.itp" [ system ] ; Name COSM in water x_COSM 0.125 [ molecules ] ; Compound #mols COSM 1374 COSG2 9620 [ploetz@cluster AddRemainingVacuumBack]$ more polff.itp ;Polarizable force field atom types ;Includes types for COS/G2 (water) COS/M (MeOH 1-site) and CPC (MeOH 2-site) [ defaults ] LJ Geometric [ atomtypes ] ;name mass charge ptype c6 c12 WO 15.99940 0.0 A 0.0 0.0 WH 1.00800 0.0 A 0.0 0.0 WD 0.00000 0.0 D 0.0 0.0 CSMO 15.99940 0.0 A 0.0 0.0 CSMH 1.00800 0.0 A 0.0 0.0 CSMC 15.03500 0.0 A 0.0 0.0 CPCO 15.99940 0.0 A 0.0 0.0 CPCH 1.00800 0.0 A 0.0 0.0 CPCC 15.03500 0.0 A 0.0 0.0 WS 0.00000 0.0 S 0.0 0.0 CSMS 0.00000 0.0 S 0.0 0.0 CPCSC 0.00000 0.0 S 0.0 0.0 CPCSO 0.00000 0.0 S 0.0 0.0 [ nonbond_params ] ;NB: mix of COSM and CPC not included ;NB: includes special COSMC - COS/G2 C12 parameters WO WO 1 3.24434e-03 3.45765e-06 WO CSMO 1 2.71125e-03 3.00305e-06 WO CSMC 1 5.36612e-03 7.71936e-06 WO CPCO 1 2.68847e-03 2.74273e-06 WO CPCC 1 6.35721e-03 10.57112e-06 CSMO CSMO 1 2.26576e-03 2.60823e-06 CSMO CSMC 1 4.48440e-03 7.10600e-06 CSMC CSMC 1 8.87552e-03 19.36000e-06 CPCO CPCO 1 2.22784e-03 2.17563e-06 CPCO CPCC 1 5.26799e-03 8.38538e-06 CPCC CPCC 1 12.45679e-03 32.31923e-06 [ploetz@cluster AddRemainingVacuumBack]$ [ploetz@cluster AddRemainingVacuumBack]$ more cosm.itp ; Topology for COS/M methanol model ; Yu et al. J., Comput. Chem. 27 (1494-1504), 2006 [ moleculetype ] ; molname nrexcl COSM 2 [ atoms ] ; id at type res nr residu name at name cg nr charge 1 CSMO 1 MOH OM 1 7.470 2 CSMH 1 MOH HM 1 0.360 3 CSMC 1 MOH CM 1 0.170 4 CSMS 1 MOH SOM 1 -8.000 [ polarization ] ;Center Shell funct alpha (nm^3) 1 4 1 0.001320 [ constraints ] ; ai aj funct length 1 2 2 0.1000 1 3 2 0.1430 2 3 2 0.1988 [ exclusions ] ; iatom excluded from interaction with i 1 2 3 4 2 1 3 4 3 1 2 4 4 1 2 3 #ifdef POSRES [ position_restraints ] ; iatom type fx fy fz 1 1 100 100 100 3 1 100 100 100 #endif [ploetz@cluster AddRemainingVacuumBack]$ [ploetz@cluster AddRemainingVacuumBack]$ more cosg2.itp ; Topology for COS/G2 water model ; Yu and Van Gunsteren, J. Chem. Phys. 121 (9549-9564), 2004 [ moleculetype ] ; molname nrexcl COSG2 2 [ atoms ] ; id at type res nr residu name at name cg nr charge 1 WO 1 SOL OW 1 0.0000 2 WH 1 SOL HW1 1 0.5265 3 WH 1 SOL HW2 1 0.5265 4 WD 1 SOL MW 1 6.9470 5 WS 1 SOL SW 1 -8.0000 [ polarization ] ;Center Shell funct alpha (nm^3) 4 5 1 0.001255 [ settles ] ; i funct dOH dHH 1 1 0.09572 0.15139 [ dummies3 ] ; The position of the dummies is computed as follows: ; ; O ; ; D ; ; H H ; ; 2 * b = distance (OD) / [ cos (angle(DOH)) * distance (OH) ] ; 0.022 nm / [ cos (104.52 / 2 deg) * 0.09572 nm ] ; 0.022 nm / 0.058588228 nm ; Dummy pos x4 = x1 + a*(x2-x1) + b*(x3-X1) ; Dummy from funct a b 4 1 2 3 1 0.187751028 0.187751028 [ exclusions ] ; iatom excluded from interaction with i 1 2 3 4 5 2 1 3 4 5 3 1 2 4 5 4 1 2 3 5 5 1 2 3 4 #ifdef POSRES [ position_restraints ] ; iatom type fx fy fz 1 1 100 100 100 #endif [ploetz@cluster AddRemainingVacuumBack]$ -- View this message in context: http://gromacs.5086.x6.nabble.com/Gromacs-4-6-4-5-3-qualitative-differences-4-6-instability-in-polarizable-force-field-vacuum-liquid-ms-tp5012174.html Sent from the GROMACS Users Forum mailing list archive at Nabble.com. -- 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