Hi,<br><br>I'm using GROMACS 4.0.5 to do a MARTINI simulation of a CG protein in a CG water box, looking at the diffusion constant of the protein among other things. I'm using the Nose-Hoover thermostat and Parrinello-Rahman barostat with tc-grps = System. In CHARMM, a protein-in-water simulation with an extended ensemble will show equal temperatures for the protein and the water. However, I'm finding that, with ref_t = 323, I get a SOL (both W and WF particles) temperature of 323 and a Protein temperature of 295. Is this what I should expect?<br>
<br>The system has 62.5k particles, the box is 20 x 20 x 20, and the relevant .mdp parameters I've used are:<br><br>tcoupl = Nose-Hoover<br>tc-grps = System<br>tau_t = 2<br>
ref_t = 323<br>Pcoupl = Parrinello-Rahman<br>Pcoupltype = isotropic<br>tau_p = 5<br>compressibility = 5e-5<br>ref_p = 1.0<br>
<br>These parameters seem to lead to stable systems for bilayer simulations. I've included the whole .mdp file at the end in case I've forgotten to mention something relevant.<br><br>I calculated temperatures with<br>
<br>g_traj -f md2.trr -s md2.tpr -n index.ndx -ot -ng 5<br><br>Thanks,<br><br>-Michael<br><br>------ begin md2.mdp ------<br><br>; RUN CONTROL PARAMETERS =<br>integrator = md<br>; start time and timestep in ps =<br>
tinit = 0.0<br>dt = 0.020<br>nsteps = 25500000<br>; mode for center of mass motion removal<br>comm-mode = linear<br>; number of steps for center of mass motion removal<br>
nstcomm = 1<br>; group(s) for center of mass motion removal<br>comm-grps = System<br><br><br>; OUTPUT CONTROL OPTIONS =<br>; Output frequency for coords (x), velocities (v) and forces (f) =<br>
nstxout = 50000<br>nstvout = 50000<br>nstfout = 0<br>; Output frequency for energies to log file and energy file =<br>nstlog = 2500<br>nstenergy = 2500<br>
; Output frequency and precision for xtc file =<br>nstxtcout = 2500<br>xtc_precision = 100<br>; This selects the subset of atoms for the xtc file. You can =<br>; select multiple groups. By default all atoms will be written. =<br>
xtc-grps =<br>; Selection of energy groups =<br>energygrps = System Protein W WF<br><br>; NEIGHBORSEARCHING PARAMETERS =<br>; nblist update frequency =<br>nstlist = 5<br>; ns algorithm (simple or grid) =<br>
ns_type = grid<br>; Periodic boundary conditions: xyz or none =<br>pbc = xyz<br>; nblist cut-off =<br>rlist = 1.2<br><br>; OPTIONS FOR ELECTROSTATICS AND VDW =<br>
; Method for doing electrostatics =<br>coulombtype = Shift<br>rcoulomb_switch = 0.0<br>rcoulomb = 1.2<br>; Dielectric constant (DC) for cut-off or DC of reaction field =<br>epsilon_r = 15<br>
; Method for doing Van der Waals =<br>vdw_type = Shift<br>; cut-off lengths =<br>rvdw_switch = 0.9<br>rvdw = 1.2<br>; Apply long range dispersion corrections for Energy and Pressure =<br>
DispCorr = No<br>; Spacing for the PME/PPPM FFT grid =<br>fourierspacing = 0.12<br>; FFT grid size, when a value is 0 fourierspacing will be used =<br>fourier_nx = 10<br>fourier_ny = 10<br>
fourier_nz = 10<br>; EWALD/PME/PPPM parameters =<br>pme_order = 4<br>ewald_rtol = 1e-05<br>epsilon_surface = 0<br>optimize_fft = no<br><br>; OPTIONS FOR WEAK COUPLING ALGORITHMS =<br>
; Temperature coupling =<br>tcoupl = Nose-Hoover<br>; Groups to couple separately =<br>tc-grps = System<br>; Time constant (ps) and reference temperature (K) =<br>tau_t = 2<br>
ref_t = 323<br>; Pressure coupling =<br>Pcoupl = Parrinello-Rahman<br>Pcoupltype = isotropic<br>; Time constant (ps), compressibility (1/bar) and reference P (bar) =<br>
tau_p = 5<br>compressibility = 5e-5<br>ref_p = 1.0<br><br>; SIMULATED ANNEALING CONTROL =<br>annealing = no<br><br>; GENERATE VELOCITIES FOR STARTUP RUN =<br>continuation = yes<br>
gen_vel = no<br>;gen_temp = 323<br>;gen_seed = 473529<br><br>; OPTIONS FOR BONDS =<br>constraints = none<br>; Type of constraint algorithm =<br>constraint_algorithm = Lincs<br>
; Relative tolerance of shake =<br>shake_tol = 0.0001<br>; Highest order in the expansion of the constraint coupling matrix =<br>lincs_order = 4<br>; Lincs will write a warning to the stderr if in one step a bond =<br>
; rotates over more degrees than =<br>lincs_warnangle = 30<br>; Convert harmonic bonds to morse potentials =<br>morse = no<br><br>; NMR refinement stuff =<br>; Distance restraints type: No, Simple or Ensemble =<br>
disre = No<br>; Force weighting of pairs in one distance restraint: Equal or Conservative =<br>disre_weighting = Equal<br>; Use sqrt of the time averaged times the instantaneous violation =<br>
disre_mixed = no<br>disre_fc = 1000<br>disre_tau = 1.25<br>; Output frequency for pair distances to energy file =<br>nstdisreout = 100<br>------- end md2.mdp -------<br>
<br><br clear="all"><br>-- <br>Michael Lerner, Ph.D.<br>IRTA Postdoctoral Fellow<br>Laboratory of Computational Biology NIH/NHLBI<br>5635 Fishers Lane, Room T909, MSC 9314<br>Rockville, MD 20852 (UPS/FedEx/Reality)<br>Bethesda MD 20892-9314 (USPS)<br>