Hi gromacs list, I'm about to start some REMD simulations using generalized Born solvent on a protein of about 5000 atoms. I have two questions, the first of which is about gromacs, the second more about REMD in general.
(1) I'm getting some pretty ugly energy drift (300K->500K in 1 ns) for an NVE MD test simulation using a 2 fs time step. It goes away if I use 1 fs, however I was under the impression that 2 fs is normally OK. I was wondering whether that could be caused by the use of the cut-off method which is required with the coloumb and VdW interactions when using GBSA? Or perhaps I'm doing something else wrong. I'll include the mdp file I'm using at the bottom of the email in case anyone feels like pointing out my foolishness to me. (2) I've been analyzing some data from an REMD simulation by my predecessor and see very slow replica flow rates. They are about two orders of magnitude smaller than the idealized rate of the exchange attempt frequency multiplied by the acceptance fraction (exchanges are attempted every 2 ps with a 0.4 acceptance fraction). If I look at the energy distribution for a given replica/temperature combination over a time scale of around 1 ns, it is clearly shifted from the average energy distribution for that temperature. The timescale for changes in this energy shift is around 10 ns. My current theory for the slow rate of replica flow is that the slow fluctuations in the energy of the protein are limiting replica flow, since a replica with lower than average energy will tend to remain at the bottom of the temperature range, while those with higher than average energies will tend to remain at the top. Has anyone else observed this kind of behavior? Is my reasoning wrong in any obvious way? Cheers, Ben ; An energy drift simulation of TaHSP. ; 7.3.2 Preprocssing ; ------------------ ; defines pass to the preprocessor define = ; 7.3.3 Run Control ; ----------------- ; group(s) for center of mass motion removal comm_grps = System ; Use the MD integrator (as opposed to minimization). integrator = md ; maximum number of steps to integrate nsteps = 10000 ; remove center of mass translation and rotation around centre of mass comm_mode = Angular ; [ps] time step for integration dt = 0.002 ; [steps] frequency of mass motion removal nstcomm = 10 ; [ps] starting time for run tinit = 0 ; 7.3.8 Output Control ; -------------------- ; [steps] freq to write velocities to trajectory nstvout = 0 ; [steps] freq to write energies to log file nstlog = 100 ; Write to energy file frequently. nstenergy = 100 ; group(s) to write to xtc trajectory xtc_grps = System ; [real] precision to write xtc trajectory xtc_precision = 1000 ; [steps] freq to write coordinates to xtc trajectory nstxtcout = 1000 ; [steps] freq to write coordinates to trajectory nstxout = 0 ; group(s) to write to energy file energygrps = System ; [steps] freq to write forces to trajectory nstfout = 0 ; 7.3.9 Neighbour Searching ; ------------------------- ; [nm] cut-off distance for the short-range neighbor list ; For Generalized Born this must be equal to the cut-off length for ; the born radius calculation. rlist = 1.6 ; method of updating neighbor list ns_type = grid ; [steps] freq to update neighbor list nstlist = 1 ; no periodic boundary conditions pbc = no ; 7.3.10 Electrostatics ; --------------------- ; apply a cut-off to electostatic coulombtype = cut-off ; coloumb cutoff radius rcoulomb = 1.6 ; 7.3.11 VdW ; ---------- ; Dispersion correction makes no sense without box size. DispCorr = no ; twin-range cut-off with rlist where rvdw >rlist vdwtype = cut-off ; Increasing VdW cutoff to same as everything else. rvdw = 1.6 ; 7.3.13 Ewald ; ------------ ; [nm] grid spacing for FFT grid when using PME fourierspacing = 0.12 ; relative strength of Ewald-shifted potential at rcoulomb ewald_rtol = 1e-05 ; interpolation order for PME, 4 cubic pme_order = 4 ; 7.3.14 Temperature Coupling ; --------------------------- ; Temperature. ref_t = 300 ; temperature coupling frequency nsttcouple = 1 ; Doing NVE simulation so no thermostat. tcoupl = no ; couple everything to same bath tc_grps = System ; [ps] time constant for coupling tau_t = 0.1 ; 7.3.15 Pressure Coupling ; ------------------------ ; [bar] reference pressure for coupling ref_p = 1.0 ; pressure coupling in x-y-z directions pcoupltype = isotropic ; [ps] time constant for coupling tau_p = 2.0 ; no pressure coupling if using generalized born pcoupl = no ; [bar^-1] compressibility compressibility = 4.5e-05 ; 7.3.17 Velocity Generation ; -------------------------- ; velocity generation turned off gen_vel = no ; 7.3.18 Bonds ; ------------ ; apply constraints to the start configuration continuation = yes ; [degrees] maximum angle that a bond can rotate before LINCS will complain lincs_warnangle = 30 ; highest order in the expansion of the contraint coupling matrix lincs_order = 4 ; number of iterations to correct for rotational lengthening lincs_iter = 1 ; LINear Constraint Solver constraint_algorithm = LINCS ; Bonds to H atoms are constraints constraints = hbonds ; 7.3.22 NMR refinement (dihedral restraints not in manual) ; --------------------------------------------------------- ; scaling factor for force constant for dihedral constraints dihre-fc = 10 ; No dihedral constraints are set. dihre = no ; 7.3.27 Implicit Solvent ; ----------------------- ; The cutoff radius for calculating the Born radii (must be equal to rlist). rgbradii = 1.6 ; Use a Generalized Born Solvent implicit_solvent = GBSA ; Use the Onufriev-Bashford-Case method to calculate the Born radii gb_algorithm = OBC ; Dielectric constant for implicit solvent. ; Default is 80 but we use 78.5 for consistency with amber. gb_epsilon_solvent = 78.5 -- 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