Dear Justin, Thanks and I tried your suggestion, that is minimizing the system without restraints and increasing the Fmax to 1000, the mdp file used is as follows:
; minim.mdp - used as input into grompp to generate em.tpr ; Parameters describing what to do, when to stop and what to save integrator = steep ; Algorithm (steep = steepest descent minimization) emtol = 1000.0 ; Stop minimization when the maximum force < 1000.0 kJ/mol/nm emstep = 0.01 ; Energy step size nsteps = 50000 ; Maximum number of (minimization) steps to perform ; Parameters describing how to find the neighbors of each atom and how to calculate the interactions nstlist = 1 ; Frequency to update the neighbor list and long range forces ns_type = grid ; Method to determine neighbor list (simple, grid) rlist = 1.2 ; Cut-off for making neighbor list (short range forces) coulombtype = PME ; Treatment of long range electrostatic interactions rcoulomb = 1.2 ; Short-range electrostatic cut-off rvdw = 1.2 ; Short-range Van der Waals cut-off pbc = xyz ; Periodic Boundary Conditions (yes/no) and the ouput is : Steepest Descents converged to Fmax < 1000 in 1192 steps Potential Energy = -4.5402134e+05 Maximum force = 9.6024463e+02 on atom 3573 Norm of force = 3.8811920e+01 Later is used this as input for the NVT equillibration with the following nvt.mdp file: title = NVT equilibration define = -DPOSRES ; position restrain the protein ; Run parameters integrator = md ; leap-frog integrator nsteps = 50000 ; 2 * 50000 = 100 ps dt = 0.002 ; 2 fs ; Output control nstxout = 100 ; save coordinates every 0.2 ps nstvout = 100 ; save velocities every 0.2 ps nstenergy = 100 ; save energies every 0.2 ps nstlog = 100 ; update log file every 0.2 ps ; Bond parameters continuation = no ; first dynamics run constraint_algorithm = lincs ; holonomic constraints constraints = all-bonds ; all bonds (even heavy atom-H bonds) constrained lincs_iter = 1 ; accuracy of LINCS lincs_order = 4 ; also related to accuracy ; Neighborsearching ns_type = grid ; search neighboring grid cels nstlist = 5 ; 10 fs rlist = 1.2 ; short-range neighborlist cutoff (in nm) rcoulomb = 1.2 ; short-range electrostatic cutoff (in nm) rvdw = 1.2 ; short-range van der Waals cutoff (in nm) ; Electrostatics coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics pme_order = 4 ; cubic interpolation fourierspacing = 0.16 ; grid spacing for FFT ; Temperature coupling is on tcoupl = V-rescale ; modified Berendsen thermostat tc-grps = Protein DPPC SOL_CL- ; three coupling groups - more accurate tau_t = 0.1 0.1 0.1 ; time constant, in ps ref_t = 323 323 323 ; reference temperature, one for each group, in K ; Pressure coupling is off pcoupl = no ; no pressure coupling in NVT ; Periodic boundary conditions pbc = xyz ; 3-D PBC ; Dispersion correction DispCorr = EnerPres ; account for cut-off vdW scheme ; Velocity generation gen_vel = yes ; assign velocities from Maxwell distribution gen_temp = 323 ; temperature for Maxwell distribution gen_seed = -1 ; generate a random seed ; COM motion removal ; These options remove motion of the protein/bilayer relative to the solvent/ions nstcomm = 1 comm-mode = Linear comm-grps = Protein_DPPC SOL_CL- and the ouput nvt.log file is as under: Input Parameters: integrator = md nsteps = 50000 init_step = 0 ns_type = Grid nstlist = 5 ndelta = 2 nstcomm = 1 comm_mode = Linear nstlog = 100 nstxout = 100 nstvout = 100 nstfout = 0 nstenergy = 100 nstxtcout = 0 init_t = 0 delta_t = 0.002 xtcprec = 1000 nkx = 48 nky = 48 nkz = 60 pme_order = 4 ewald_rtol = 1e-05 ewald_geometry = 0 epsilon_surface = 0 optimize_fft = FALSE ePBC = xyz bPeriodicMols = FALSE bContinuation = FALSE bShakeSOR = FALSE etc = V-rescale epc = No epctype = Isotropic tau_p = 1 ref_p (3x3): ref_p[ 0]={ 0.00000e+00, 0.00000e+00, 0.00000e+00} ref_p[ 1]={ 0.00000e+00, 0.00000e+00, 0.00000e+00} ref_p[ 2]={ 0.00000e+00, 0.00000e+00, 0.00000e+00} compress (3x3): compress[ 0]={ 0.00000e+00, 0.00000e+00, 0.00000e+00} compress[ 1]={ 0.00000e+00, 0.00000e+00, 0.00000e+00} compress[ 2]={ 0.00000e+00, 0.00000e+00, 0.00000e+00} refcoord_scaling = No posres_com (3): posres_com[0]= 0.00000e+00 posres_com[1]= 0.00000e+00 posres_com[2]= 0.00000e+00 posres_comB (3): posres_comB[0]= 0.00000e+00 posres_comB[1]= 0.00000e+00 posres_comB[2]= 0.00000e+00 andersen_seed = 815131 rlist = 1.2 rtpi = 0.05 coulombtype = PME rcoulomb_switch = 0 rcoulomb = 1.2 vdwtype = Cut-off rvdw_switch = 0 rvdw = 1.2 epsilon_r = 1 epsilon_rf = 1 tabext = 1 implicit_solvent = No gb_algorithm = Still gb_epsilon_solvent = 80 nstgbradii = 1 rgbradii = 2 gb_saltconc = 0 gb_obc_alpha = 1 gb_obc_beta = 0.8 gb_obc_gamma = 4.85 sa_surface_tension = 2.092 DispCorr = EnerPres free_energy = no init_lambda = 0 sc_alpha = 0 sc_power = 0 sc_sigma = 0.3 delta_lambda = 0 nwall = 0 wall_type = 9-3 wall_atomtype[0] = -1 wall_atomtype[1] = -1 wall_density[0] = 0 wall_density[1] = 0 wall_ewald_zfac = 3 pull = no disre = No disre_weighting = Conservative disre_mixed = FALSE dr_fc = 1000 dr_tau = 0 nstdisreout = 100 orires_fc = 0 orires_tau = 0 nstorireout = 100 dihre-fc = 1000 em_stepsize = 0.01 em_tol = 10 niter = 20 fc_stepsize = 0 nstcgsteep = 1000 nbfgscorr = 10 ConstAlg = Lincs shake_tol = 1e-04 lincs_order = 4 lincs_warnangle = 30 lincs_iter = 1 bd_fric = 0 ld_seed = 1993 cos_accel = 0 deform (3x3): deform[ 0]={ 0.00000e+00, 0.00000e+00, 0.00000e+00} deform[ 1]={ 0.00000e+00, 0.00000e+00, 0.00000e+00} deform[ 2]={ 0.00000e+00, 0.00000e+00, 0.00000e+00} userint1 = 0 userint2 = 0 userint3 = 0 userint4 = 0 userreal1 = 0 userreal2 = 0 userreal3 = 0 userreal4 = 0 grpopts: nrdf: 6048.01 12219 45018 ref_t: 323 323 323 tau_t: 0.1 0.1 0.1 anneal: No No No ann_npoints: 0 0 0 acc: 0 0 0 nfreeze: N N N energygrp_flags[ 0]: 0 efield-x: n = 0 efield-xt: n = 0 efield-y: n = 0 efield-yt: n = 0 efield-z: n = 0 efield-zt: n = 0 bQMMM = FALSE QMconstraints = 0 QMMMscheme = 0 scalefactor = 1 qm_opts: ngQM = 0 Table routines are used for coulomb: TRUE Table routines are used for vdw: FALSE Will do PME sum in reciprocal space. ++++ PLEASE READ AND CITE THE FOLLOWING REFERENCE ++++ U. Essman, L. Perela, M. L. Berkowitz, T. Darden, H. Lee and L. G. Pedersen A smooth particle mesh Ewald method J. Chem. Phys. 103 (1995) pp. 8577-8592 -------- -------- --- Thank You --- -------- -------- Using a Gaussian width (1/beta) of 0.384195 nm for Ewald Cut-off's: NS: 1.2 Coulomb: 1.2 LJ: 1.2 System total charge: 0.000 Generated table with 1100 data points for Ewald. Tabscale = 500 points/nm Generated table with 1100 data points for LJ6. Tabscale = 500 points/nm Generated table with 1100 data points for LJ12. Tabscale = 500 points/nm Generated table with 1100 data points for 1-4 COUL. Tabscale = 500 points/nm Generated table with 1100 data points for 1-4 LJ6. Tabscale = 500 points/nm Generated table with 1100 data points for 1-4 LJ12. Tabscale = 500 points/nm Enabling SPC water optimization for 7497 molecules. Configuring nonbonded kernels... Testing AMD 3DNow support... not present. Testing ia32 SSE support... present. Removing pbc first time Initializing LINear Constraint Solver ++++ PLEASE READ AND CITE THE FOLLOWING REFERENCE ++++ B. Hess and H. Bekker and H. J. C. Berendsen and J. G. E. M. Fraaije LINCS: A Linear Constraint Solver for molecular simulations J. Comp. Chem. 18 (1997) pp. 1463-1472 -------- -------- --- Thank You --- -------- -------- The number of constraints is 9045 ++++ PLEASE READ AND CITE THE FOLLOWING REFERENCE ++++ S. Miyamoto and P. A. Kollman SETTLE: An Analytical Version of the SHAKE and RATTLE Algorithms for Rigid Water Models J. Comp. Chem. 13 (1992) pp. 952-962 -------- -------- --- Thank You --- -------- -------- Center of mass motion removal mode is Linear We have the following groups for center of mass motion removal: 0: Protein_DPPC 1: SOL_CL- ++++ PLEASE READ AND CITE THE FOLLOWING REFERENCE ++++ G. Bussi, D. Donadio and M. Parrinello Canonical sampling through velocity rescaling J. Chem. Phys. 126 (2007) pp. 014101 -------- -------- --- Thank You --- -------- -------- There are: 31609 Atoms Max number of connections per atom is 28 Total number of connections is 138849 Max number of graph edges per atom is 4 Total number of graph edges is 48078 Constraining the starting coordinates (step 0) Constraining the coordinates at t0-dt (step 0) RMS relative constraint deviation after constraining: 4.83e-04 Initial temperature: 333.793 K Started mdrun on node 0 Sun Jul 5 18:04:45 2009 Step Time Lambda 0 0.00000 0.00000 Grid: 10 x 10 x 13 cells Long Range LJ corr.: <C6> 9.7327e-04 Long Range LJ corr.: Epot -2200.19, Pres: -136.406, Vir: 2200.19 Energies (kJ/mol) Angle G96Angle Proper Dih. Ryckaert-Bell. Improper Dih. 1.19670e+04 8.38636e+03 8.15057e+03 4.20390e+03 3.08430e+03 LJ-14 Coulomb-14 LJ (SR) Disper. corr. Coulomb (SR) 4.45274e+03 5.37063e+04 3.47509e+06 -2.20019e+03 -4.22052e+05 Coul. recip. Position Rest. Potential Kinetic En. Total Energy -1.63954e+05 1.03877e+01 2.98084e+06 5.03020e+09 5.03318e+09 Conserved En. Temperature Pressure (bar) Cons. rmsd () 5.03318e+09 1.91196e+07 1.03974e+08 9.31478e+00 and if i remove the the constraints and lowering the time step using the following nvt.mdp file it is running: title = NVT equilibration define = -DPOSRES ; position restrain the protein ; Run parameters integrator = md ; leap-frog integrator nsteps = 100000 ; 1 * 100000 = 100 ps dt = 0.001 ; 1 fs ; Output control nstxout = 100 ; save coordinates every 0.2 ps nstvout = 100 ; save velocities every 0.2 ps nstenergy = 100 ; save energies every 0.2 ps nstlog = 100 ; update log file every 0.2 ps ; Bond parameters continuation = no ; first dynamics run constraint_algorithm = lincs ; holonomic constraints constraints = none ; all bonds (even heavy atom-H bonds) constrained lincs_iter = 1 ; accuracy of LINCS lincs_order = 4 ; also related to accuracy ; Neighborsearching ns_type = grid ; search neighboring grid cels nstlist = 5 ; 10 fs rlist = 1.2 ; short-range neighborlist cutoff (in nm) rcoulomb = 1.2 ; short-range electrostatic cutoff (in nm) rvdw = 1.2 ; short-range van der Waals cutoff (in nm) ; Electrostatics coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics pme_order = 4 ; cubic interpolation fourierspacing = 0.16 ; grid spacing for FFT ; Temperature coupling is on tcoupl = V-rescale ; modified Berendsen thermostat tc-grps = Protein DPPC SOL_CL- ; three coupling groups - more accurate tau_t = 0.1 0.1 0.1 ; time constant, in ps ref_t = 323 323 323 ; reference temperature, one for each group, in K ; Pressure coupling is off pcoupl = no ; no pressure coupling in NVT ; Periodic boundary conditions pbc = xyz ; 3-D PBC ; Dispersion correction DispCorr = EnerPres ; account for cut-off vdW scheme ; Velocity generation gen_vel = yes ; assign velocities from Maxwell distribution gen_temp = 323 ; temperature for Maxwell distribution gen_seed = -1 ; generate a random seed ; COM motion removal ; These options remove motion of the protein/bilayer relative to the solvent/ions nstcomm = 1 comm-mode = Linear comm-grps = Protein_DPPC SOL_CL- please suggest me is it ok to remove the constraints and run the NVT equillibration. Thanks Ram On Fri, Oct 16, 2009 at 9:27 PM, Justin A. Lemkul <jalem...@vt.edu> wrote: > > > ram bio wrote: >> >> Dear Gromacs users, >> >> I am doing protein in lipid-bilayer simulation and i am following the >> procedure as per justin tutorial. I am able to insert the protein in >> lipid bilayer and minimize the system as per Inflategro >> procedure,during the total procedure the system was minimized in every >> step.Then, I solvated and ionized sytem and minimized using the >> following mdp file: >> >> ; minim.mdp - used as input into grompp to generate em.tpr >> ; Parameters describing what to do, when to stop and what to save >> define = -DSTRONG_POSRES >> integrator = steep ; Algorithm (steep = steepest descent >> minimization) >> emtol = 500.0 ; Stop minimization when the maximum >> force < 1000.0 kJ/mol/nm >> emstep = 0.01 ; Energy step size >> nsteps = 50000 ; Maximum number of (minimization) >> steps to perform >> >> ; Parameters describing how to find the neighbors of each atom and how >> to calculate the interactions >> nstlist = 1 ; Frequency to update the neighbor >> list and long range forces >> ns_type = grid ; Method to determine neighbor list >> (simple, grid) >> rlist = 1.2 ; Cut-off for making neighbor list >> (short range forces) >> coulombtype = PME ; Treatment of long range >> electrostatic interactions >> rcoulomb = 1.2 ; Short-range electrostatic cut-off >> rvdw = 1.2 ; Short-range Van der Waals cut-off >> pbc = xyz ; Periodic Boundary Conditions (yes/no) >> >> the output was as follows: >> >> Steepest Descents converged to Fmax < 500 in 4770 steps >> Potential Energy = -3.8820288e+05 >> Maximum force = 4.4803549e+02 on atom 3573 >> Norm of force = 1.7854408e+01 >> >> As the potential energy and Fmax values were agreeable , I proceeded >> to equillibrate the system using NVT. >> > > Did you minimize the structure without restraints, prior to NVT? > > <snip> > >> Angle G96Angle Proper Dih. Ryckaert-Bell. Improper Dih. >> 2.20077e+04 8.54042e+03 6.78950e+03 4.34650e+03 3.03266e+03 >> LJ-14 Coulomb-14 LJ (SR) Disper. corr. Coulomb (SR) >> 4.76527e+03 5.50236e+04 8.36617e+06 -2.20019e+03 -4.46844e+05 >> Coul. recip. Position Rest. Potential Kinetic En. Total Energy >> -1.65524e+05 1.07769e+01 7.85612e+06 5.88813e+10 5.88892e+10 >> Conserved En. Temperature Pressure (bar) Cons. rmsd () >> 5.88892e+10 2.23805e+08 1.21713e+09 3.10245e+01 >> > > In my experience, the combination of an astronomically high temperature and > a repulsive temperature is indicative of restraining an unrestrainable > starting structure. Try the EM I suggested above. Other than that, as I've > suggested before, see the Advanced Troubleshooting page I created in the > tutorial. > > -Justin > >> Please help me to proceed further and let me know where are the >> mistakes lying and how to overcome them. >> >> Thanks in advance, >> >> Ram >> _______________________________________________ >> 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 >> > > -- > ======================================== > > Justin A. Lemkul > Ph.D. Candidate > ICTAS Doctoral Scholar > Department of Biochemistry > Virginia Tech > Blacksburg, VA > jalemkul[at]vt.edu | (540) 231-9080 > http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin > > ======================================== > _______________________________________________ > 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 > _______________________________________________ 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