Re: Constraint failures and TI electrostatic calcs

From: Courtney Taylor (courtney.b.taylor_at_gmail.com)
Date: Mon Dec 06 2010 - 14:18:40 CST

I would like to note that I just did a test where I added 10 steps of
minimization and it will run for 1000 steps. I don't really want to add
minimization though, since I am trying to start from the same coordinates
each lambda value, and minimization would change this (or at least
minimization DURING TI is not recommended in any of the good practices and
tutorials I have read)

Any thoughts?

Courtney Taylor
PhD Candidate
Department of Chemical and Biomolecular Engineering
Vanderbilt University
(225)-771-8554
(615)-343-3257

"Though the course may change sometimes the river always meets the
sea"--Zeppelin

On Mon, Dec 6, 2010 at 1:26 PM, Courtney Taylor <courtney.b.taylor_at_gmail.com
> wrote:

> All,
>
> I am working on TI calculations that mutate tyrosine to tryptophan on a
> larger protein backbone in a 50A x 50A x 50A box of water. I minimize the
> system for 10,000 steps, and then equilibrate for 2 ns with a 1 fs timestep.
> This portion of the calculation is fine. However, when I try to perform
> electrostatics calculations where I turn on and off the charges of the
> reactant and product, I get the immediate error:
>
> ERROR: Constraint failure in RATTLE algorithm for atom 513!
> ERROR: Constraint failure; simulation has become unstable.
>
> The constraint failures always occur within my hybrid atoms where TI is
> being implemented. I have looked over the mailing list for solutions to
> this, but the most I found was to increase the minimization steps, which I
> did. I also increased my PME parameters from 60x60x60 to 75x75x75. I've
> tried even lower timesteps and other adjustments to no avail. If I turn off
> SHAKE, I get the error that the atoms are moving too fast. The vdw portion
> of my calcs (where alchElecstart = 1 and alchvdwend = 1) appear so far to
> work correctly with no crashes in the first 20,000 steps
>
> My configuration file for lambda = 0 is below. I have run two other full
> calculations, tyrosine to alanine and tyrosine to tryptophan, with no issues
> at all. These are smaller mutations though, so I'm assuming that is what
> helped me there.
>
> #############################################################
> ## ADJUSTABLE PARAMETERS ##
> #############################################################
>
> set temperature 300
> set lambda 0.0
> set lambdaend 0.1
> set outputname y5wnc_elec_$lambda
> set dLambda 0.1
>
> firsttimestep 0
>
> #############################################################
> ## SIMULATION PARAMETERS ##
> #############################################################
>
> # Input
>
> paratypecharmm on
> parameters input/par_all27_prot_lipid.prm
> structure input/cbm_nc_y5w_namd.psf
> coordinates input/cbm_nc_y5w_namd.pdb
> Bincoordinates input/y5wnc_equil2.coor
> Binvelocities input/y5wnc_equil2.vel
> extendedsystem input/y5wnc_equil2.xsc
>
> #temperature $temperature
>
> # Force-Field Parameters
> exclude scaled1-4
> 1-4scaling 1.0
> cutoff 10.0
> switching on
> switchdist 9
> pairlistdist 13
> rigidBonds all
> rigidTolerance 1e-6
>
> # Integrator Parameters
> timestep 1.0
> stepspercycle 1
>
> # Periodic Boundary Conditions declare if xsc file unavailable
> #cellBasisVector1 50 0 0
> #cellBasisVector2 0 50 0
> #cellBasisVector3 0 0 50
> #cellOrigin 0 0 0
> # PME (for full-system periodic electrostatics)
> PME yes
> PMEInterpOrder 6
> PMETolerance 1e-6
> PMEGridSizeX 75
> PMEGridSizeY 75
> PMEGridSizeZ 75
>
> wrapall on
>
> # Constant Temperature Control
> langevin on ;# do langevin dynamics
> langevinDamping 5 ;# damping coefficient (gamma) of 5/ps
> langevinTemp $temperature
> langevinHydrogen off
>
> # Output
> outputName $outputname
> restartfreq 20000
> dcdfreq 10000
> xstFreq 20000
> outputEnergies 10000
>
>
> #Thermo Integration commands
> source fep.tcl
> alch on #Enable alch mode
> alchType ti #Use TI method
> alchfile input/y5wnc_alch.pdb #PDB file with alch flag
> alchCol B #Alch flag column
> alchOutfile $outputname.out #TI output file
> alchOutFreq 10000 #Output frequency
> alchEquilSteps 500000 #Equilibration before data collection
>
> #For elec only ElecLamStart = 0, vdwend = 0
> #For vdw only Elecstart = 1, vdwend = 1
> alchVdWShiftCoeff 5 #VDW soft core potential coeff (A^2)
> alchElecLambdaStart 0 #Introduce electrostatics
> alchDecouple on # Nonb interactions w environment scaled,
> within subset not scaled
> alchVdwLambdaEnd 0 #value to cancel vdw
>
> alchLambda $lambda # start lambda
> alchLambda2 $lambdaend # end lambda
> run 10000000
>
> Thanks for any advice you can give! I'm a bit stumped on what to do.
> Electrostatics calcs in the past have been the easiest to complete, so I'm
> surprised to run into this error now.
>
> Courtney Taylor
> PhD Candidate
> Department of Chemical and Biomolecular Engineering
> Vanderbilt University
>
>
>
>

This archive was generated by hypermail 2.1.6 : Wed Feb 29 2012 - 15:56:28 CST