Error in thermodynamic integration for calculating ligand absolute binding free energy

From: Marawan Hussien (marawanhussain_at_yahoo.com)
Date: Mon Jul 17 2017 - 11:45:22 CDT

Hi,I am having trouble running absolute binding free energy calculation for a solvated protein ligand complex built with AMBER.Here is my NAMD configuration file:
######parmfile          com_solvated.topambercoor  eq_5_npt.rst7

bincoordinates    com_2.restart.coorbinVelocities     com_2.restart.velextendedsystem    com_2.restart.xsc
set temperature           300;
set outputname com_4;

outputName         $outputname; # base name for output from this run                                       # NAMD writes two files at the end, final coord and vel                                       # in the format of first-dyn.coor and first-dyn.velrestartfreq        1000;                # 500 steps = every 1psdcdfreq           1000;

xstFreq           1000;                # XSTFreq: control how often the extended systen configuration                                       # will be appended to the XST fileoutputEnergies     1000;                # 125 steps = every 0.25ps                                       # The number of timesteps between each energy output of NAMDoutputTiming      1000;                # The number of timesteps between each timing output shows                                       # time per step and time to completion
# Force-Field Parametersamber onuseSettle      on rigidTolerance 1.0e-8 cutoff         13.5 pairlistdist   16.0 switching      off exclude        scaled1-4 readexclusions yes 1-4scaling     0.83333333 scnb           2.0 zeromomentum   on ljcorrection   onwatermodel     tip3
stepspercycle       20;                # 20 redo pairlists every ten stepspairlistsPerCycle    2;                # 2 is the default                                        # cycle represents the number of steps between atom reassignments                                       # this means every 20/2=10 steps the pairlist will be updated                       # Integrator Parameterstimestep            2.0;               # fs/steprigidBonds          all;               # Bound constraint all bonds involving H are fixed in lengthnonbondedFreq       1;                 # nonbonded forces every stepfullElectFrequency  1;                 # PME every step

wrapWater            on;               # wrap water to central cellwrapAll              on;               # wrap other molecules too
# PME (for full-system periodic electrostatics)PME                yes;PMEGridSpacing     1;margin 6
# Pressure and volume controluseGroupPressure       yes;            # use a hydrogen-group based pseudo-molecular viral to calcualte pressure and                                       # has less fluctuation, is needed for rigid bonds (rigidBonds/SHAKE)

langevin                onlangevinDamping         1
langevinHydrogen        off
# constant pressurelangevinPiston          onlangevinPistonTarget    1.01325langevinPistonPeriod    50.0langevinPistonDecay     25.0langevinTemp 310langevinPistonTemp      300

minimize 500run 1000
source fep.tcl#alch                  onalchFile              bound.fepalchCol               BalchOutFreq           10alchOutFile           system.fepoutalchEquilSteps        500#set nSteps      2000#
runFEP 0.0 1.0 0.1 $numSteps #####

Whenever I run the file, I got this error:
#####FINISHED WRITING RESTART VELOCITIESThe last velocity output (seq=1000) takes 0.014 seconds, 202.688 MB of memory in useInfo: Benchmark time: 1024 CPUs 0.0103227 s/step 0.0597378 days/ns 202.688 MB memoryTCL: Setting parameter alch to onFATAL ERROR: Setting parameter alch from script failed!FATAL ERROR: Setting parameter alch from script failed!#####

Nothing helped, I tried some suggestions, such as including this keyword:global at k r0

Or adding explicit parameters for softcore potentials, but in vain. I also tried to use the same input while imposing constrains on the ligand (5 kcal/mol) cartesian constrains, but this also did not help. I am not sure if there is any incompatibility between the tcl forces from the configuraion file and the fep.tcl script. Any suggestion?
Regards,Marawan

This archive was generated by hypermail 2.1.6 : Sun Dec 31 2017 - 23:21:28 CST