Dear Mike, It's hard to know what's going on from your input. Did you check the thermodynamic components with g_energy (especially the pressure bond, LJ and Coulomb energies)? This may give you a hint. Also, I would do a test run with PME - you're anyway not following the parametrisation of Jorgensen and co-workers exactly. As a final note - the correct phase is something that's very difficult to reproduce exactly with classical MD, and while a single molecule interacting with a protein may work just fine, a bulk of such molecules does not always reproduce the correct phase at a given temperature. This is one reason why lipid models are still being developed.
Hope that helps a bit, Ran Mike Wykes wrote: > Dear All > > I would like to perform MD simulations with benzene as a solvent and > am observing some strange behaviour when I use Lincs to constrain all > the bonds. When I run a fully flexible NPT MD of a box of 320 benzene > molecules simulation at 298K and 1bar, the density comes out at 841 > g/l, not too far away from the experimental value of 876 g/l. However > when I constrain all bonds using Lincs, the system expands rapidly, > stabilising at a density of 2.62 g/l ! Both simulations started from > NVT equilibrated simulations fixed to the experimental density. In the > papers describing OPLS parametrisation, the MC simulations were indeed > performed with fully flexible molecules, but it surprises me that the > bond constraints would affect the density so strongly. Does anyone > have any ideas why this is occurring? I am using a cutoff of 1.5 nm > for VDW and Coulomb interactions without EWALD/PME but this is > consistent with how OPLS was parametrised. There are no Lincs warnings > in the log file of the constrained simulation. > > Please find my mdp and Benzene itp files below, the only difference > between the flexible and constrained runs being dt = 0.001/0.002 and > constraints = none/all-bonds respectively. > > Many thanks for your ideas/explanations as to what could be going on, > > Mike > > ; > ; File 'mdout.mdp' was generated > ; By user: mwykes (7017) > ; On host: node168 > ; At date: Thu Feb 4 20:05:46 2010 > ; > > ; VARIOUS PREPROCESSING OPTIONS > ; Preprocessor information: use cpp syntax. > ; e.g.: -I/home/joe/doe -I/home/mary/hoe > include = > ; e.g.: -DI_Want_Cookies -DMe_Too > define = -DFLEX_SPC > > ; RUN CONTROL PARAMETERS > integrator = md > ; Start time and timestep in ps > tinit = 0 > dt = 0.002 > nsteps = 5000000 > ; For exact run continuation or redoing part of a run > ; Part index is updated automatically on checkpointing (keeps files separate) > simulation_part = 1 > 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 = 1.0 > emstep = 0.1 > ; Max number of iterations in relax_shells > niter = 20 > ; 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 = 0 > nstvout = 0 > nstfout = 0 > ; Output frequency for energies to log file and energy file > nstlog = 500 > nstenergy = 500 > ; 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 = BNZ > ; Selection of energy groups > energygrps = BNZ > > ; NEIGHBORSEARCHING PARAMETERS > ; nblist update frequency > nstlist = 10 > ; ns algorithm (simple or grid) > nstype = grid > ; Periodic boundary conditions: xyz, no, xy > pbc = xyz > periodic_molecules = no > ; nblist cut-off > rlist = 1.5 > > ; OPTIONS FOR ELECTROSTATICS AND VDW > ; Method for doing electrostatics > coulombtype = Cut-off > rcoulomb-switch = 0 > rcoulomb = 1.5 > ; 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 ALGORITHM > implicit_solvent = 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 > ; Dielectric coefficient of the implicit solvent > gb_epsilon_solvent = 80 > ; Salt concentration in M for Generalized Born models > gb_saltconc = 0 > ; Scaling factors used in the OBC GB model. Default values are OBC(II) > gb_obc_alpha = 1 > gb_obc_beta = 0.8 > gb_obc_gamma = 4.85 > ; Surface tension (kJ/mol/nm^2) for the SA (nonpolar surface) part of GBSA > ; The default value (2.092) corresponds to 0.005 kcal/mol/Angstrom^2. > sa_surface_tension = 2.092 > > ; OPTIONS FOR WEAK COUPLING ALGORITHMS > ; Temperature coupling > tcoupl = berendsen > ; Groups to couple separately > tc-grps = BNZ > ; Time constant (ps) and reference temperature (K) > tau_t = 0.1 > ref_t = 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 = > ; 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 = no > gen-temp = 300 > gen-seed = -1 > > ; OPTIONS FOR BONDS > constraints = all-bonds > ; Type of constraint algorithm > constraint-algorithm = Lincs > ; Do not constrain the start configuration > continuation = 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 = 1 > ; 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_r_linpot = -1 > wall_atomtype = > wall_density = > wall_ewald_zfac = 3 > > ; COM PULLING > ; Pull type: no, umbrella, constraint or 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 = 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) and S to energy file > nstorireout = 100 > ; Dihedral angle restraints: No or Yes > dihre = no > dihre_fc = 9999.0 > > ; Free energy control stuff > free-energy = no > init-lambda = 0 > delta-lambda = 0 > sc-alpha = 0 > sc-power = 0 > sc-sigma = 0.3 > couple-moltype = > couple-lambda0 = vdw-q > couple-lambda1 = vdw-q > couple-intramol = no > > ; 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 > > benzene.itp: > > [ defaults ] > ; nbfunc comb-rule gen-pairs fudgeLJ fudgeQQ > 1 3 yes 0.5 0.5 > > [ atomtypes ] > ; full atom descriptions are available in ffoplsaa.atp > ; name bond_type mass charge ptype sigma epsilon > opls_145 CA 6 12.01100 -0.115 A 3.55000e-01 2.92880e-01 > opls_146 HA 1 1.00800 0.115 A 2.42000e-01 1.25520e-01 > > [ bondtypes ] > ; i j func b0 kb > CA HA 1 0.10800 307105.6 ; PHE, etc. > CA CA 1 0.14000 392459.2 ; TRP,TYR,PHE > > [ angletypes ] > ; i j k func th0 cth > CA CA HA 1 120.000 292.880 ; > CA CA CA 1 120.000 527.184 ; PHE(OL) > > [ dihedraltypes ] > ; i j k l func coefficients > ; OPLS Fourier dihedraltypes translated to Gromacs Ryckaert-Bellemans form > ; according to the formula in the Gromacs manual. > X CA CA X 3 30.33400 0.00000 -30.33400 > 0.00000 0.00000 0.00000 ; aromatic ring > > [ dihedraltypes ] > ; Improper OPLS dihedrals to keep groups planar. > ; (OPLS doesnt use impropers for chiral atoms). > ; Since these functions are periodic of the form 1-cos(2*x), the are actually > ; implemented as proper dihedrals [1+cos(2*x+180)] for the moment, > ; to keep things compatible. > ; The defines are used in ffoplsaa.rtp or directly in your .top ile. > #define improper_Z_CA_X_Y 180.0 4.60240 2 > > > [ moleculetype ] > ; Name nrexcl > BNZ 3 > > [ atoms ] > ; nr type resnr residue atom cgnr charge mass > typeB chargeB massB > 1 opls_145 1 BNZ CA1 1 -0.115 12.011 ; -0.1150 > 2 opls_145 1 BNZ CA2 2 -0.115 12.011 ; -0.2300 > 3 opls_145 1 BNZ CA3 3 -0.115 12.011 ; -0.3450 > 4 opls_145 1 BNZ CA4 4 -0.115 12.011 ; -0.4600 > 5 opls_145 1 BNZ CA5 5 -0.115 12.011 ; -0.5750 > 6 opls_145 1 BNZ CA6 6 -0.115 12.011 ; -0.6900 > 7 opls_146 1 BNZ HA7 1 0.115 1.008 ; -0.5750 > 8 opls_146 1 BNZ HA8 2 0.115 1.008 ; -0.4600 > 9 opls_146 1 BNZ HA9 3 0.115 1.008 ; -0.3450 > 10 opls_146 1 BNZ HA10 4 0.115 1.008 ; -0.2300 > 11 opls_146 1 BNZ HA11 5 0.115 1.008 ; -0.1150 > 12 opls_146 1 BNZ HA12 6 0.115 1.008 ; -0.0000 > > [ bonds ] > ; ai aj funct c0 c1 c2 c3 > 1 2 1 > 1 6 1 > 1 7 1 > 2 3 1 > 2 8 1 > 3 4 1 > 3 9 1 > 4 5 1 > 4 10 1 > 5 6 1 > 5 11 1 > 6 12 1 > > [ angles ] > ; ai aj ak funct c0 c1 c2 > c3 > 2 1 6 1 > 2 1 7 1 > 6 1 7 1 > 1 2 3 1 > 1 2 8 1 > 3 2 8 1 > 2 3 4 1 > 2 3 9 1 > 4 3 9 1 > 3 4 5 1 > 3 4 10 1 > 5 4 10 1 > 4 5 6 1 > 4 5 11 1 > 6 5 11 1 > 1 6 5 1 > 1 6 12 1 > 5 6 12 1 > > [ dihedrals ] > ; ai aj ak al funct c0 c1 > c2 c3 c4 c5 > 6 1 2 3 3 > 6 1 2 8 3 > 7 1 2 3 3 > 7 1 2 8 3 > 2 1 6 5 3 > 2 1 6 12 3 > 7 1 6 5 3 > 7 1 6 12 3 > 1 2 3 4 3 > 1 2 3 9 3 > 8 2 3 4 3 > 8 2 3 9 3 > 2 3 4 5 3 > 2 3 4 10 3 > 9 3 4 5 3 > 9 3 4 10 3 > 3 4 5 6 3 > 3 4 5 11 3 > 10 4 5 6 3 > 10 4 5 11 3 > 4 5 6 1 3 > 4 5 6 12 3 > 11 5 6 1 3 > 11 5 6 12 3 > > [ dihedrals ] > ; ai aj ak al funct c0 c1 > c2 c3 > 2 6 1 7 1 improper_Z_CA_X_Y > 6 2 1 7 1 improper_Z_CA_X_Y > 1 3 2 8 1 improper_Z_CA_X_Y > 3 1 2 8 1 improper_Z_CA_X_Y > 2 4 3 9 1 improper_Z_CA_X_Y > 4 2 3 9 1 improper_Z_CA_X_Y > 3 5 4 10 1 improper_Z_CA_X_Y > 5 3 4 10 1 improper_Z_CA_X_Y > 4 6 5 11 1 improper_Z_CA_X_Y > 6 4 5 11 1 improper_Z_CA_X_Y > 1 5 6 12 1 improper_Z_CA_X_Y > 5 1 6 12 1 improper_Z_CA_X_Y > > [ pairs ] > 6 3 1 > 6 8 1 > 7 3 1 > 7 8 1 > 2 5 1 > 2 12 1 > 7 5 1 > 7 12 1 > 1 4 1 > 1 9 1 > 8 4 1 > 8 9 1 > 2 5 1 > 2 10 1 > 9 5 1 > 9 10 1 > 3 6 1 > 3 11 1 > 10 6 1 > 10 11 1 > 4 1 1 > 4 12 1 > 11 1 1 > 11 12 1 > -- 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