Dear all, by doing MD simulations of a protein complex embedded in 2 membranes (inner and outer membrane of bacteria, POPE), we observe a bilayer splitting in one of the membranes. This has the effect that the bilayer forms "bubbles" with vacuum inside.
It might have something to do with the PME handling of GroMACS 4.0.3 (GROMOS96-53a6 ff). - there is no splitting using a 1 nm smaller boxsize in exactly the same model - if we switch off the PME mode and use only Cut-off's we don't have any splitting effect - the behavior of the bilayer splitting depends on the number of used cluster/cpu nodes Changing the "fourierspacing" parameter to create different PME grids has no effect to avoid bilayer splitting. Changing: thermostat | barostat (semiisotropic): - berendsen | berendsen: "splitting" - v-rescale | parrinello-rahman: "splitting"/"keep together" (50:50) Here's the mdp file of our last try: title = Yep, sometimes I will cause a bilayer splitting; ; The following lines tell the program the standard locations where to find certain files cpp = cpp; Preprocessor include = -I../top; Directories to include in the topology format define = -DPOSRES; Apply position restraints. ; RUN CONTROL PARAMETERS integrator = MD ; Start time and timestep in ps tinit = 0 dt = 0.002 nsteps = 10000000 ; 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 = Protein POPE_center POPE_inner POPE_outer SOL_NA+_CL- ; ENERGY MINIMIZATION OPTIONS ; Force tolerance and initial step-size emtol = 500 emstep = 0.01 ; Max number of iterations in relax_shells niter = 0 ; Step size (1/ps^2) for minimization of flexible constraints fcstep = 0 ; OUTPUT CONTROL OPTIONS ; Output frequency for coords (x), velocities (v) and forces (f) nstxout = 500000 nstvout = 500000 nstfout = 0 ; Checkpointing helps you continue after crashes nstcheckpoint = 50000 ; Output frequency for energies to log file and energy file nstlog = 5000 nstenergy = 5000 ; Output frequency and precision for xtc file nstxtcout = 50000 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 POPE_center POPE_inner POPE_outer SOL_NA+_CL- ; Selection of energy groups energygrps = Protein POPE_center POPE_inner POPE_outer SOL_NA+_CL- ; 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.15 domain-decomposition = no ; OPTIONS FOR ELECTROSTATICS AND VDW ; Method for doing electrostatics coulombtype = PME rcoulomb-switch = 0 rcoulomb = 1.15 ; 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 = 1.4 ; Apply long range dispersion corrections for Energy and Pressure DispCorr = No ; 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 = yes ; OPTIONS FOR WEAK COUPLING ALGORITHMS ; Temperature coupling ; Tcoupl = berendsen Tcoupl = V-rescale ; Groups to couple separately tc-grps = Protein POPE_center POPE_inner POPE_outer SOL_NA+_CL- ; Time constant (ps) and reference temperature (K) tau-t = 0.1 0.1 0.1 0.1 0.1 ref-t = 310 310 310 310 310 ; Pressure coupling ; Pcoupl = berendsen Pcoupl = Parrinello-Rahman Pcoupltype = semiisotropic ; Time constant (ps), compressibility (1/bar) and reference P (bar) tau-p = 4 4 compressibility = 4.5e-5 4.5e-5 ref-p = 1.0 1.0 ; Random seed for Andersen thermostat andersen_seed = 815131 ; GENERATE VELOCITIES FOR STARTUP RUN gen_vel = yes gen-temp = 310 gen-seed = 24071998 ; OPTIONS FOR BONDS constraints = all-bonds ; Type of constraint algorithm constraint-algorithm = Lincs ; Do not constrain the start configuration unconstrained-start = 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 = 4 ; Lincs will write a warning to the stderr if in one step a bond ; rotates over more degrees than lincs-warnangle = 90 ; Convert harmonic bonds to morse potentials morse = no Does anyone have any idea? Many thanks for your time and any help. Cheers, Thomas -- Thomas H. Schmidt, PhD student Computational Structural Biology Life Science Informatics, B-IT LIMES-Institute, University of Bonn Dahlmannstrasse 2, D-53113 Bonn, Germany Phone: +49-(0)228-2699 323 Fax: +49-(0)228-2699 341 Email: schm...@bit.uni-bonn.de http://www.csb.bit.uni-bonn.de -- 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/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/mailing_lists/users.php