From: Jiyong Park (jiyong_at_mbi.ucla.edu)
Date: Mon Dec 10 2012 - 01:41:28 CST
Dear Users,
I came up with a suspicious dG values as I'm trying
to use Amber force field (FF99SB) in order to examine
dG of point mutation (Thr to Ala) of a protein.
I'm wondering if NAMD can handle AmberFF for FEP (or TI)
simulation. According to messages in namd-l archive,
one may handle this problem by setting up a proper PDB
file with beta-column indicating which atoms will be
annihilated and exnihilated. Details of my simulation
are explained below.
An alchemified (M2A) residue was prepared (see below)
then created a solvated system using tleap which produced
prmtop and inpcrd files. Also a pdb file was generated
to indicate which atoms will be annihilated (exnihilated).
Then initiated FEP simulation using NAMD v2.9
(NAMD_2.9_Linux-x86_64-ibverbs), after 1 ns equilibration.
The problem recognized was the magnitude of dG values
printed out looks too large to take. I reconize VdW energy
difference dominated overall non-covalent energy difference.
Below is a portion of FEP output file.
# STEP Elec vdW
dE dE_avg Temp dG
# l l+dl l
l+dl E(l+dl)-E(l)
FepEnergy: 50140 -158610.0634 -158608.2400 17985.9864
17916.7191 -67.4438 -67.8857 298.1082 -67.9218
FepEnergy: 50150 -158629.7685 -158627.9711 17876.9226
17807.6339 -67.4913 -67.8856 299.1902 -67.9217
FepEnergy: 50160 -159082.2776 -159080.4780 18454.3617
18384.8846 -67.6775 -67.8855 297.1587 -67.9216
FepEnergy: 50170 -158889.4018 -158887.5642 18208.2813
18138.9362 -67.5076 -67.8854 297.7808 -67.9215
FepEnergy: 50180 -158636.7925 -158635.2613 18066.9853
17997.6663 -67.7877 -67.8853 296.5642 -67.9214
FepEnergy: 50190 -158713.0381 -158711.4967 17915.1132
17845.7462 -67.8256 -67.8853 299.2589 -67.9214
Here lamda was 0.0 and delta lambda was 0.02. Corresponding test
run with param22+CMAP resulted dG ~ -0.3 kcal/mol. I guess dG ~ -68
kcal/mol cannot be a sound number.
ALCH section of NAMD input script and alchemified residue library
is included for your comprehension.
1) NAMD input
alchVdwLambdaEnd 0.7
alchElecLambdaStart 0.5
alchVdWShiftCoeff 5.0
alchDecouple on
2) alchemified library file (M2A.lib)
!!index array str
"T2A"
!entry.T2A.unit.atoms table str name str type int typex int resx
int flags int seq int elmnt dbl chg
"N" "N" 0 1 131072 1 7 -0.415700
"H" "H" 0 1 131072 2 1 0.271900
"CAA" "CT" 0 1 131072 3 6 -0.038900 (XXA : outgoing atoms)
"HAA" "H1" 0 1 131072 4 1 0.100700
"CBA" "CT" 0 1 131072 5 6 0.365400
"HBA" "H1" 0 1 131072 6 1 0.004300
"CG2A" "CT" 0 1 131072 7 6 -0.243800
"HG21A" "HC" 0 1 131072 8 1 0.064200
"HG22A" "HC" 0 1 131072 9 1 0.064200
"HG23A" "HC" 0 1 131072 10 1 0.064200
"OG1A" "OH" 0 1 131072 11 8 -0.676100
"HG1A" "HO" 0 1 131072 12 1 0.410200
"C" "C" 0 1 131072 13 6 0.597300
"O" "O" 0 1 131072 14 8 -0.567900
"CAB" "CT" 0 1 131072 15 6 0.033700 (XXB : incoming atoms)
"HAB" "H1" 0 1 131072 16 1 0.082300
"CBB" "CT" 0 1 131072 17 6 -0.182500
"HB1B" "HC" 0 1 131072 18 1 0.060300
"HB2B" "HC" 0 1 131072 19 1 0.060300
"HB3B" "HC" 0 1 131072 20 1 0.060300
!entry.T2A.unit.atomspertinfo table str pname str ptype int ptypex
int pelmnt dbl pchg
"N" "N" 0 -1 0.0
"H" "H" 0 -1 0.0
"CAA" "CT" 0 -1 0.0
"HAA" "H1" 0 -1 0.0
"CBA" "CT" 0 -1 0.0
"HBA" "H1" 0 -1 0.0
"CG2A" "CT" 0 -1 0.0
"HG21A" "HC" 0 -1 0.0
"HG22A" "HC" 0 -1 0.0
"HG23A" "HC" 0 -1 0.0
"OG1A" "OH" 0 -1 0.0
"HG1A" "HO" 0 -1 0.0
"C" "C" 0 -1 0.0
"O" "O" 0 -1 0.0
"CAB" "CT" 0 -1 0.0
"HAB" "H1" 0 -1 0.0
"CBB" "CT" 0 -1 0.0
"HB1B" "HC" 0 -1 0.0
"HB2B" "HC" 0 -1 0.0
"HB3B" "HC" 0 -1 0.0
!entry.T2A.unit.boundbox array dbl
-1.000000
0.0
0.0
0.0
0.0
!entry.T2A.unit.childsequence single int
2
!entry.T2A.unit.connect array int
1
13
!entry.T2A.unit.connectivity table int atom1x int atom2x int flags
1 2 1
1 3 1
3 4 1
3 5 1
3 13 1
5 6 1
5 7 1
5 11 1
7 8 1
7 9 1
7 10 1
11 12 1
13 14 1
1 15 1
15 16 1
15 17 1
15 13 1
17 18 1
17 19 1
17 20 1
!entry.T2A.unit.hierarchy table str abovetype int abovex str
belowtype int belowx
"U" 0 "R" 1
"R" 1 "A" 1
"R" 1 "A" 2
"R" 1 "A" 3
"R" 1 "A" 4
"R" 1 "A" 5
"R" 1 "A" 6
"R" 1 "A" 7
"R" 1 "A" 8
"R" 1 "A" 9
"R" 1 "A" 10
"R" 1 "A" 11
"R" 1 "A" 12
"R" 1 "A" 13
"R" 1 "A" 14
"R" 1 "A" 15
"R" 1 "A" 16
"R" 1 "A" 17
"R" 1 "A" 18
"R" 1 "A" 19
"R" 1 "A" 20
!entry.T2A.unit.name single str
"T2A"
This archive was generated by hypermail 2.1.6 : Tue Dec 31 2013 - 23:22:49 CST