Re: [gmx-users] Pulling ligand - Different Profiles (Force vs time)
Thank you Justin. On Wed, Jun 27, 2012 at 2:39 PM, Justin A. Lemkul wrote: > > > On 6/27/12 9:36 AM, Steven Neumann wrote: >> >> On Wed, Jun 27, 2012 at 1:51 PM, Justin A. Lemkul wrote: >>> >>> >>> >>> On 6/27/12 7:48 AM, Steven Neumann wrote: Dear Gmx Users, I obtained a protein-ligand complex from 100ns simulation. Now I am pulling my ligand away from the protein after the energy minimzation in water and equilibration of 100ps (two coupling baths: Protein, LIG_Water_and_ions). Then I proceed my pulling : grompp -f pull.mdp -c npt.gro -p topol.top -n index.ndx -t npt.cpt -o pull.tpr mdrun -s pull.tpr -deffnm pull title = Umbrella pulling simulation define = -DPOSRES ; Run parameters integrator = md dt = 0.002 tinit = 0 nsteps = 50 ; 1 ns nstcomm = 10 ; Output parameters nstxout = 0 nstvout = 0 nstfout = 500 nstxtcout = 1000 ; every 1 ps nstenergy = 500 ; Bond parameters constraint_algorithm = lincs constraints = all-bonds continuation = yes ; continuing from NPT ; Single-range cutoff scheme nstlist = 5 ns_type = grid rlist = 1.4 rcoulomb = 1.4 rvdw = 1.2 vdwtype = Switch rvdw-switch = 1.0 ; 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 ; Temperature coupling is on tcoupl = V-rescale ; modified Berendsen thermostat tc_grps = Protein LIG_Water_and_ions ; two coupling groups - more accurate tau_t = 0.1 0.1 ; time constant, in ps ref_t = 298 298 ; reference temperature, one for each group, in K ; Pressure coupling is on Pcoupl = Parrinello-Rahman pcoupltype = isotropic tau_p = 2.0 compressibility = 4.5e-5 ref_p = 1.0 ; 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 = distance ; simple distance increase pull_dim = N N Y pull_start = yes ; define initial COM distance > 0 pull_ngroups = 1 pull_group0 = Protein pull_group1 = LIG pull_rate1 = 0.004 ; 0.004 nm per ps = 4 nm per ns pull_k1 = 500 ; kJ mol^-1 nm^-2 I run 3 pulling simulations with the same mdp and I obtain 3 different profiles (Force vs time). Then I used 2xlonger pulling with the same pulling distance and I run 3 simulations again. Each time I obtain different profile. Can anyone explain me this? I am using velocities from npt simulation as above (gen_vel = no and continuation = yes) so I presume the output should be similar. Please, advice. >>> >>> I assume you're passing a checkpoint file to grompp? If you're relying >>> on >>> velocities from the .gro file, they are of insufficient precision to >>> guarantee proper continuation. >> >> >> Thank you Justin. I am using according to your tutorial: >> >> grompp -f pull.mdp -c npt.gro -p topol.top -n index.ndx -t npt.cpt -o >> pull.tpr >> mdrun -s pull.tpr -deffnm pull >> >> Would you suggest: >> >> grompp -f pull.mdp -c npt.gro -p topol.top -n index.ndx -t npt.cpt -o >> pull.tpr >> mdrun -s pull.tpr -cpi npt.cpt -deffnm pull ?? >> > > No, I would not, especially if the NPT run uses position restraints, in > which case the two phases are different. I missed the command line in the > earlier message. What you are doing makes sense. > > >> Profiles do not vary slightly - the maximum pulling force (breaking >> point) varies from 290 to 500 kJ/mol nm2 which is really a lot. >> > > Consult the points below and watch your trajectories. If you're getting > different forces, your ligands are experiencing different interactions. SMD > is a path-dependent, non-equilibrium process. Good sampling and a > justifiable path are key. > > -Justin > > >>> >>> Small variations are inherent in any simulation set, and in the case of >>> pulling, small changes (though intentional) are the basis for Jarzynski's >>> method. In any case, all MD simulations are chaotic and so it depends on >>> what your definition of "different" is in the context of whether or not >>> there are meaningful changes imparted through the course of each >>> simulation. >>> Also note that in the absence of the -reprod flag, the same .tpr fi
Re: [gmx-users] Pulling ligand - Different Profiles (Force vs time)
On 6/27/12 9:36 AM, Steven Neumann wrote: On Wed, Jun 27, 2012 at 1:51 PM, Justin A. Lemkul wrote: On 6/27/12 7:48 AM, Steven Neumann wrote: Dear Gmx Users, I obtained a protein-ligand complex from 100ns simulation. Now I am pulling my ligand away from the protein after the energy minimzation in water and equilibration of 100ps (two coupling baths: Protein, LIG_Water_and_ions). Then I proceed my pulling : grompp -f pull.mdp -c npt.gro -p topol.top -n index.ndx -t npt.cpt -o pull.tpr mdrun -s pull.tpr -deffnm pull title = Umbrella pulling simulation define = -DPOSRES ; Run parameters integrator = md dt = 0.002 tinit = 0 nsteps = 50; 1 ns nstcomm = 10 ; Output parameters nstxout = 0 nstvout = 0 nstfout = 500 nstxtcout = 1000 ; every 1 ps nstenergy = 500 ; Bond parameters constraint_algorithm= lincs constraints = all-bonds continuation= yes ; continuing from NPT ; Single-range cutoff scheme nstlist = 5 ns_type = grid rlist = 1.4 rcoulomb= 1.4 rvdw= 1.2 vdwtype = Switch rvdw-switch = 1.0 ; 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 ; Temperature coupling is on tcoupl = V-rescale ; modified Berendsen thermostat tc_grps = Protein LIG_Water_and_ions ; two coupling groups - more accurate tau_t = 0.1 0.1 ; time constant, in ps ref_t = 298 298 ; reference temperature, one for each group, in K ; Pressure coupling is on Pcoupl = Parrinello-Rahman pcoupltype = isotropic tau_p = 2.0 compressibility = 4.5e-5 ref_p = 1.0 ; 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 = distance ; simple distance increase pull_dim= N N Y pull_start = yes ; define initial COM distance > 0 pull_ngroups= 1 pull_group0 = Protein pull_group1 = LIG pull_rate1 = 0.004 ; 0.004 nm per ps = 4 nm per ns pull_k1 = 500 ; kJ mol^-1 nm^-2 I run 3 pulling simulations with the same mdp and I obtain 3 different profiles (Force vs time). Then I used 2xlonger pulling with the same pulling distance and I run 3 simulations again. Each time I obtain different profile. Can anyone explain me this? I am using velocities from npt simulation as above (gen_vel = no and continuation = yes) so I presume the output should be similar. Please, advice. I assume you're passing a checkpoint file to grompp? If you're relying on velocities from the .gro file, they are of insufficient precision to guarantee proper continuation. Thank you Justin. I am using according to your tutorial: grompp -f pull.mdp -c npt.gro -p topol.top -n index.ndx -t npt.cpt -o pull.tpr mdrun -s pull.tpr -deffnm pull Would you suggest: grompp -f pull.mdp -c npt.gro -p topol.top -n index.ndx -t npt.cpt -o pull.tpr mdrun -s pull.tpr -cpi npt.cpt -deffnm pull ?? No, I would not, especially if the NPT run uses position restraints, in which case the two phases are different. I missed the command line in the earlier message. What you are doing makes sense. Profiles do not vary slightly - the maximum pulling force (breaking point) varies from 290 to 500 kJ/mol nm2 which is really a lot. Consult the points below and watch your trajectories. If you're getting different forces, your ligands are experiencing different interactions. SMD is a path-dependent, non-equilibrium process. Good sampling and a justifiable path are key. -Justin Small variations are inherent in any simulation set, and in the case of pulling, small changes (though intentional) are the basis for Jarzynski's method. In any case, all MD simulations are chaotic and so it depends on what your definition of "different" is in the context of whether or not there are meaningful changes imparted through the course of each simulation. Also note that in the absence of the -reprod flag, the same .tpr file may result in a slightly different outcome. The implications of these outcomes are limited by sampling; the ensemble should converge with sufficient time and/or replicates. For non-equilibrium processes like pulling, convergence is probably harder, but again you have to ask whether the differences are meaningful. http://www.gromacs.org/Documentation/Terminology/Reproducibility -Justin -- Justin A. Lemkul, Ph.D. Research Scientist Department of Biochemistry Virginia Tech Blacksburg, VA jalemkul[at]vt.edu | (540) 231-9080 http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin ===
Re: [gmx-users] Pulling ligand - Different Profiles (Force vs time)
On Wed, Jun 27, 2012 at 1:51 PM, Justin A. Lemkul wrote: > > > On 6/27/12 7:48 AM, Steven Neumann wrote: >> >> Dear Gmx Users, >> >> I obtained a protein-ligand complex from 100ns simulation. Now I am >> pulling my ligand away from the protein after the energy minimzation >> in water and equilibration of 100ps (two coupling baths: Protein, >> LIG_Water_and_ions). >> Then I proceed my pulling : >> >> grompp -f pull.mdp -c npt.gro -p topol.top -n index.ndx -t npt.cpt -o >> pull.tpr >> >> mdrun -s pull.tpr -deffnm pull >> >> >> title = Umbrella pulling simulation >> define = -DPOSRES >> ; Run parameters >> integrator = md >> dt = 0.002 >> tinit = 0 >> nsteps = 50 ; 1 ns >> nstcomm = 10 >> ; Output parameters >> nstxout = 0 >> nstvout = 0 >> nstfout = 500 >> nstxtcout = 1000 ; every 1 ps >> nstenergy = 500 >> ; Bond parameters >> constraint_algorithm = lincs >> constraints = all-bonds >> continuation = yes ; continuing from NPT >> ; Single-range cutoff scheme >> nstlist = 5 >> ns_type = grid >> rlist = 1.4 >> rcoulomb = 1.4 >> rvdw = 1.2 >> vdwtype = Switch >> rvdw-switch = 1.0 >> ; 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 >> ; Temperature coupling is on >> tcoupl = V-rescale ; modified >> Berendsen thermostat >> tc_grps = Protein LIG_Water_and_ions ; two coupling groups - more >> accurate >> tau_t = 0.1 0.1 ; time constant, >> in ps >> ref_t = 298 298 ; reference >> temperature, one for each group, in K >> ; Pressure coupling is on >> Pcoupl = Parrinello-Rahman >> pcoupltype = isotropic >> tau_p = 2.0 >> compressibility = 4.5e-5 >> ref_p = 1.0 >> ; 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 = distance ; simple distance increase >> pull_dim = N N Y >> pull_start = yes ; define initial COM distance > 0 >> pull_ngroups = 1 >> pull_group0 = Protein >> pull_group1 = LIG >> pull_rate1 = 0.004 ; 0.004 nm per ps = 4 nm per ns >> pull_k1 = 500 ; kJ mol^-1 nm^-2 >> >> I run 3 pulling simulations with the same mdp and I obtain 3 >> different profiles (Force vs time). Then I used 2xlonger pulling with >> the same pulling distance and I run 3 simulations again. Each time I >> obtain different profile. Can anyone explain me this? I am using >> velocities from npt simulation as above (gen_vel = no and continuation >> = yes) so I presume the output should be similar. Please, advice. >> > > I assume you're passing a checkpoint file to grompp? If you're relying on > velocities from the .gro file, they are of insufficient precision to > guarantee proper continuation. Thank you Justin. I am using according to your tutorial: grompp -f pull.mdp -c npt.gro -p topol.top -n index.ndx -t npt.cpt -o pull.tpr mdrun -s pull.tpr -deffnm pull Would you suggest: grompp -f pull.mdp -c npt.gro -p topol.top -n index.ndx -t npt.cpt -o pull.tpr mdrun -s pull.tpr -cpi npt.cpt -deffnm pull ?? Profiles do not vary slightly - the maximum pulling force (breaking point) varies from 290 to 500 kJ/mol nm2 which is really a lot. Steven > > Small variations are inherent in any simulation set, and in the case of > pulling, small changes (though intentional) are the basis for Jarzynski's > method. In any case, all MD simulations are chaotic and so it depends on > what your definition of "different" is in the context of whether or not > there are meaningful changes imparted through the course of each simulation. > Also note that in the absence of the -reprod flag, the same .tpr file may > result in a slightly different outcome. The implications of these outcomes > are limited by sampling; the ensemble should converge with sufficient time > and/or replicates. For non-equilibrium processes like pulling, convergence > is probably harder, but again you have to ask whether the differences are > meaningful. > > http://www.gromacs.org/Documentation/Terminology/Reproducibility > > -Justin > > -- > > > Justin A. Lemkul, Ph.D. > Research Scientist > 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 > * Only plain text messages are allowed!
Re: [gmx-users] Pulling ligand - Different Profiles (Force vs time)
On 6/27/12 7:48 AM, Steven Neumann wrote: Dear Gmx Users, I obtained a protein-ligand complex from 100ns simulation. Now I am pulling my ligand away from the protein after the energy minimzation in water and equilibration of 100ps (two coupling baths: Protein, LIG_Water_and_ions). Then I proceed my pulling : grompp -f pull.mdp -c npt.gro -p topol.top -n index.ndx -t npt.cpt -o pull.tpr mdrun -s pull.tpr -deffnm pull title = Umbrella pulling simulation define = -DPOSRES ; Run parameters integrator = md dt = 0.002 tinit = 0 nsteps = 50; 1 ns nstcomm = 10 ; Output parameters nstxout = 0 nstvout = 0 nstfout = 500 nstxtcout = 1000 ; every 1 ps nstenergy = 500 ; Bond parameters constraint_algorithm= lincs constraints = all-bonds continuation= yes ; continuing from NPT ; Single-range cutoff scheme nstlist = 5 ns_type = grid rlist = 1.4 rcoulomb= 1.4 rvdw= 1.2 vdwtype = Switch rvdw-switch = 1.0 ; 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 ; Temperature coupling is on tcoupl = V-rescale ; modified Berendsen thermostat tc_grps = Protein LIG_Water_and_ions ; two coupling groups - more accurate tau_t = 0.1 0.1 ; time constant, in ps ref_t = 298 298 ; reference temperature, one for each group, in K ; Pressure coupling is on Pcoupl = Parrinello-Rahman pcoupltype = isotropic tau_p = 2.0 compressibility = 4.5e-5 ref_p = 1.0 ; 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 = distance ; simple distance increase pull_dim= N N Y pull_start = yes ; define initial COM distance > 0 pull_ngroups= 1 pull_group0 = Protein pull_group1 = LIG pull_rate1 = 0.004 ; 0.004 nm per ps = 4 nm per ns pull_k1 = 500 ; kJ mol^-1 nm^-2 I run 3 pulling simulations with the same mdp and I obtain 3 different profiles (Force vs time). Then I used 2xlonger pulling with the same pulling distance and I run 3 simulations again. Each time I obtain different profile. Can anyone explain me this? I am using velocities from npt simulation as above (gen_vel = no and continuation = yes) so I presume the output should be similar. Please, advice. I assume you're passing a checkpoint file to grompp? If you're relying on velocities from the .gro file, they are of insufficient precision to guarantee proper continuation. Small variations are inherent in any simulation set, and in the case of pulling, small changes (though intentional) are the basis for Jarzynski's method. In any case, all MD simulations are chaotic and so it depends on what your definition of "different" is in the context of whether or not there are meaningful changes imparted through the course of each simulation. Also note that in the absence of the -reprod flag, the same .tpr file may result in a slightly different outcome. The implications of these outcomes are limited by sampling; the ensemble should converge with sufficient time and/or replicates. For non-equilibrium processes like pulling, convergence is probably harder, but again you have to ask whether the differences are meaningful. http://www.gromacs.org/Documentation/Terminology/Reproducibility -Justin -- Justin A. Lemkul, Ph.D. Research Scientist 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 listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users * Only plain text messages are allowed! * 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
