Hello, In umbrella sampling simulations I have observed unusual unfolding events, which appear after restarting simulations.
I have simulated a homo-dimer with about 150 residues for each monomer. The system has been set up with TIP3P explicit solvent and the CHARMM22* force field. In the umbrellas I use the pull code (see mdp-file below) to confine the monomers. I also use the enforced rotation module keep the orientation of each monomer separately fixed. The problem arises when I restart simulations. I restart them by using the check point file and either the original tpr file or a new tpr file for extending the simulation time (see below). In about 10 of my 200 trajectories the protein had unfolded right after restarting. The unfolding looks artificial and it happens within a few ps. Before the protein had been simulated for 20ns without any problems and the unfolding does not seem to have a physical reason. In some case almost the whole contour length of the backbone is unfolded. Also the unfolded configurations sometimes seem to be squeezed into a plane and happen strangly symmetric with respect to the two monomers. In other cases the unfolding is not as severe and only local, but still seems artificial from visual observations of the protein structure. For restarting I use: mdrun -deffnm umbrella -pf pullf-umbrella.xvg -px pullx-umbrella.xvg -cpi umbrella.cpt (The tpr file is called umbrella.tpr) or: mdrun -s umbrella_ext.tpr -deffnm umbrella -pf pullf-umbrella.xvg -px pullx-umbrella.xvg -cpi umbrella.cpt The umbrella_ext.tpr is created by: tpbconv -s umbrella.tpr -extend 10000 -o umbrella_ext.tpr I was using gromacs-4.6.2 and 4.6.5. My system has about 130k atoms and was using between 96 and 144 CPU's. I was wondering if the issue might be related to the enforced rotation module, since I have only experienced the unfolding events in the case where I was using this module. Please find the content of my mdp-file below. I would be happy to know, if anyone help me to find out what is the issue. Cheers, Daniel MDP file: title = Umbrella sampling explicit solvent ; Run parameters integrator = md dt = 0.002 tinit = 0 nsteps = 10000000 ; 20 ns comm-mode = None ; Output parameters nstxout = 50000 ; every 100 ps nstvout = 50000 nstfout = 50000 nstxtcout = 5000 ; every 10 ps nstenergy = 5000 ; Bond parameters constraint_algorithm = lincs constraints = all-bonds continuation = yes ; Single-range cutoff scheme nstlist = 5 ns_type = grid rlist = 1.4 rcoulomb = 1.4 rvdw = 1.4 ; PME electrostatics parameters coulombtype = PME fourierspacing = 0.12 fourier_nx = 0 fourier_ny = 0 fourier_nz = 0 pme_order = 4 ewald_rtol = 1e-5 optimize_fft = yes ; Berendsen temperature coupling is on in two groups Tcoupl = Nose-Hoover tc_grps = Protein Non-Protein tau_t = 0.5 0.5 ref_t = 300 300 ; Pressure coupling is on Pcoupl = Parrinello-Rahman pcoupltype = isotropic tau_p = 1.0 compressibility = 4.5e-5 ref_p = 1.0 refcoord_scaling = com ; Generate velocities is off gen_vel = no ; Periodic boundary conditions are on in all directions pbc = xyz ; Long-range dispersion correction DispCorr = EnerPres ; Pull code pull = umbrella pull_geometry = position pull_dim = Y Y Y ; umbrella potential in all three dimensions pull_start = yes ; define initial COM distance > 0 pull_ngroups = 1 pull_group0 = freeze pull_group1 = pull pull_vec1 = 0 0 1 ; pull along z pull_rate1 = 0.0 ; nm per ps pull_k1 = 10000 ; kJ mol^-1 nm^-2 pull-nstxout = 10 pull-nstfout = 10 ; ENFORCED ROTATION ; Enforced rotation: No or Yes rotation = Yes ; Output frequency for angle, torque and rotation potential energy for the whole group rot-nstrout = 10 ; Output frequency for per-slab data (angles, torques and slab centers) rot-nstsout = 10 ; Number of rotation groups rot-ngroups = 2 ; Rotation group name rot-group0 = freeze_Ca rot-group1 = pull_Ca ; Rotation potential. Can be iso, iso-pf, pm, pm-pf, rm, rm-pf, rm2, rm2-pf, flex, flex-t, flex2, flex2-t rot-type0 = rm2-pf rot-type1 = rm2-pf ; Use mass-weighting of the rotation group positions rot-massw0 = no ; only Ca, all masses equal rot-massw1 = no ; Rotation vector, will get normalized rot-vec0 = 1 0 0 ; specify arbitrary vector, rot mat is id regardless, because rot-rate = 0 rot-vec1 = 1 0 0 ; Rotation rate [degree/ps] and force constant [kJ/(mol*nm^2)] rot-rate0 = 0.0 rot-k0 = 1000.0 rot-rate1 = 0.0 rot-k1 = 1000.0 ; Fitting method to determine angle of rotation group (rmsd, norm, or potential) rot-fit-method0 = rmsd rot-fit-method1 = rmsd ; Value of additive constant epsilon' [nm^2] for rm2* and flex2* potentials rot-eps0 = 0.01 rot-eps1 = 0.01 -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.