Re: [gmx-users] Problems in putting restraints during free energy perturbation calculations
Hi Mark, Thanks for reminding me about this. Yes, I've read up on this on a tutorial and I probably interchanged which interactions should be turned on first. Thank you again. Have a wonderful day! Searle On Mon, Mar 12, 2018 at 12:32 PM, Mark Abraham wrote: > Hi, > > Don't do that. You have charges on atoms that have no VDW, so there is > nothing to stop them flying all over the place and imparting crazy forces > on other atoms that have VDW, etc. Turn the VDW fully on before adding any > charge. This is a fairly well known issue but something that people > regularly get wrong, too. (Is there maybe a tutorial you did that missed > giving a pointer on this issue?) > > Mark > > On Mon, Mar 12, 2018 at 5:21 PM Searle Duay wrote: > > > Good day! > > > > I am doing some free energy perturbation to calculate the free energy of > > binding of an ion to a peptide. For one of my processes, I am slowly > > turning on the coulombic interactions by an increment of 0.05 lambda, > > followed by turning on the van der Waals interactions by an increment of > > 0.05 lambda. This is being done while a pull restraint is present to keep > > the ion bound to the peptide. Following is a copy of the MDP file of my > > equilibration run under NPT ensemble: > > > > ; Run control > > integrator = sd ; Langevin dynamics > > tinit= 0 > > dt = 0.0005 > > nsteps = 20; 100 ps > > nstcomm = 200 > > ; Output control > > nstxout = 1000 > > nstvout = 1000 > > nstfout = 0 > > nstlog = 1000 > > nstenergy= 1000 > > nstxout-compressed = 0 > > ; Neighborsearching and short-range nonbonded interactions > > cutoff-scheme= Verlet > > nstlist = 20 > > ns-type = grid > > pbc = xyz > > rlist= 1.2 > > ; Electrostatics > > coulombtype = pme > > rcoulomb = 1.2 > > ; van der Waals > > vdwtype = Cut-off > > vdw-modifier = Force-switch > > rvdw-switch = 1.0 > > rvdw = 1.2 > > ; Apply long range dispersion corrections for Energy and Pressure > > DispCorr = EnerPres > > ; Spacing for the PME/PPPM FFT grid > > fourierspacing = 0.12 > > ; EWALD/PME/PPPM parameters > > pme-order= 6 > > ewald-rtol = 1e-06 > > epsilon-surface = 0 > > ; Temperature coupling > > ; tcoupl is implicitly handled by the sd integrator > > tcoupl = berendsen > > tc-grps = PROT SOL_ION > > tau-t= 1.01.0 > > ref-t= 300 300 > > ; Pressure coupling is on for NPT > > pcoupl = berendsen > > pcoupltype = isotropic > > tau-p= 5.0 > > compressibility = 4.5e-5 4.5e-5 > > ref-p= 1.0 1.0 > > ; Free energy control stuff > > free-energy = yes > > init-lambda-state= 1 > > delta-lambda = 0 > > calc-lambda-neighbors= 1; only immediate neighboring windows > > ; Vectors of lambda specified here > > ; Each combination is an index that is retrieved from init_lambda_state > for > > each simulation > > ; init_lambda_state0123456789 > > 10 11 12 13 14 15 16 17 18 19 20 21 22 23 > 24 > > 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 > > 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 > > 55 56 57 58 59 60 > > vdw-lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 > > 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.05 0.10 > 0.15 > > 0.20 0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.60 0.65 0.70 0.75 0.80 0.85 > 0.90 > > 0.95 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 > 1.00 > > 1.00 1.00 1.00 1.00 1.00 1.00 1.00 > > coul-lambdas = 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 > > 0.45 0.50 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95 1.00 1.00 1.00 > 1.00 > > 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 > 1.00 > > 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 > 1.00 > > 1.00 1.00 1.00 1.00 1.00 1.00 1.00 > > ; We are not transforming any bonded or restrained interactions > > bonded-lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 > > 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 > 0.00 > > 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 > 0.00 > > 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 > 0.00 > > 0.00 0.00 0.00 0.00 0.00 0.00 0.00 > > restraint-lambdas= 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 > > 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.
Re: [gmx-users] Problems in putting restraints during free energy perturbation calculations
Hi, Don't do that. You have charges on atoms that have no VDW, so there is nothing to stop them flying all over the place and imparting crazy forces on other atoms that have VDW, etc. Turn the VDW fully on before adding any charge. This is a fairly well known issue but something that people regularly get wrong, too. (Is there maybe a tutorial you did that missed giving a pointer on this issue?) Mark On Mon, Mar 12, 2018 at 5:21 PM Searle Duay wrote: > Good day! > > I am doing some free energy perturbation to calculate the free energy of > binding of an ion to a peptide. For one of my processes, I am slowly > turning on the coulombic interactions by an increment of 0.05 lambda, > followed by turning on the van der Waals interactions by an increment of > 0.05 lambda. This is being done while a pull restraint is present to keep > the ion bound to the peptide. Following is a copy of the MDP file of my > equilibration run under NPT ensemble: > > ; Run control > integrator = sd ; Langevin dynamics > tinit= 0 > dt = 0.0005 > nsteps = 20; 100 ps > nstcomm = 200 > ; Output control > nstxout = 1000 > nstvout = 1000 > nstfout = 0 > nstlog = 1000 > nstenergy= 1000 > nstxout-compressed = 0 > ; Neighborsearching and short-range nonbonded interactions > cutoff-scheme= Verlet > nstlist = 20 > ns-type = grid > pbc = xyz > rlist= 1.2 > ; Electrostatics > coulombtype = pme > rcoulomb = 1.2 > ; van der Waals > vdwtype = Cut-off > vdw-modifier = Force-switch > rvdw-switch = 1.0 > rvdw = 1.2 > ; Apply long range dispersion corrections for Energy and Pressure > DispCorr = EnerPres > ; Spacing for the PME/PPPM FFT grid > fourierspacing = 0.12 > ; EWALD/PME/PPPM parameters > pme-order= 6 > ewald-rtol = 1e-06 > epsilon-surface = 0 > ; Temperature coupling > ; tcoupl is implicitly handled by the sd integrator > tcoupl = berendsen > tc-grps = PROT SOL_ION > tau-t= 1.01.0 > ref-t= 300 300 > ; Pressure coupling is on for NPT > pcoupl = berendsen > pcoupltype = isotropic > tau-p= 5.0 > compressibility = 4.5e-5 4.5e-5 > ref-p= 1.0 1.0 > ; Free energy control stuff > free-energy = yes > init-lambda-state= 1 > delta-lambda = 0 > calc-lambda-neighbors= 1; only immediate neighboring windows > ; Vectors of lambda specified here > ; Each combination is an index that is retrieved from init_lambda_state for > each simulation > ; init_lambda_state0123456789 > 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 > 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 > 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 > 55 56 57 58 59 60 > vdw-lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 > 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.05 0.10 0.15 > 0.20 0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 > 0.95 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 > 1.00 1.00 1.00 1.00 1.00 1.00 1.00 > coul-lambdas = 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 > 0.45 0.50 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95 1.00 1.00 1.00 1.00 > 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 > 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 > 1.00 1.00 1.00 1.00 1.00 1.00 1.00 > ; We are not transforming any bonded or restrained interactions > bonded-lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 > 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 > 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 > 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 > 0.00 0.00 0.00 0.00 0.00 0.00 0.00 > restraint-lambdas= 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 > 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 > 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 > 1.00 1.00 0.95 0.90 0.85 0.80 0.75 0.70 0.65 0.60 0.55 0.50 0.45 0.40 0.35 > 0.30 0.25 0.20 0.15 0.10 0.05 0.00 > ; Masses are not changing (particle identities are the same at lambda = 0 > and lambda = 1) > mass-lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 > 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 > 0.00 0.00 0.0
