Constraint failures and TI electrostatic calcs

From: Courtney Taylor (courtney.b.taylor_at_gmail.com)
Date: Mon Dec 06 2010 - 13:26:52 CST

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 - 05:23:28 CST