http://chem-siepmann.oit.umn.edu/siepmann/trappe/index.html
Overall, the potential energy function takes the following form:
\begin{equation*} \ U_{total} = \sum_{angles}{\frac{k_{\theta}}2 (\theta - \theta_{eq} )^{2}} + \sum_{dihedrals}{c_{0} + c_{1}[1 + cos(\phi)] + c_{2}[1 - cos(2\phi)] + c_{3}[1 + cos(3\phi)]} + \sum_{i=1}^{N-1}{\sum_{j=i+1}^{N}{ 4\epsilon_{ij}[(\frac{\sigma_{ij}}r_{ij})^{12} - (\frac{\sigma_{ij}}r_{ij})^6] }} \end{equation*}Here, CHx indicates linear alkane carbons while Ccx indicates carbons in cyclic rings of size x.
(psuedo)atom | type | $ \frac{\epsilon}k_B $ [K] | $ \sigma $ [Angstrom] | q [e] |
---|---|---|---|---|
CH4 | CH4 | 148.0 | 3.730 | 0.000 |
CH3 | [CH3]-CHx | 98.0 | 3.750 | 0.000 |
CH2 | CHx-[CH2]-CHx | 46.0 | 3.950 | 0.000 |
Cc5 | CH2-[CH2]-CH2 | 56.3 | 3.880 | 0.000 |
Cc6 | CH2-[CH2]-CH2 | 52.5 | 3.910 | 0.000 |
TraPPE uses fixed bond lengths. We will be adding flexible bonds in our implementation, which are based on the parameters for sp3 carbons in the general amber force field (gaff)
Type | Length [Angstrom] | $k_{b}$ (kcal/mol) |
---|---|---|
CHx-CHy | 1.540 | 300.9 |
Type | $ \theta $ | $ \frac{k_{\theta}}{k_{B}} [K/rad^2] $ |
---|---|---|
CHx-(CH2)-CHx | 114.0 | 62500 |
Cc5-(Cc5)-Cc5 | 105.5 | 62500 |
Cc6-(Cc6)-Cc6 | 114.0 | 62500 |
Because of the different form of the cyclic alkanes, we must transformt the equation...
\begin{equation*} \ u_{torsion_{cyclic}}(\phi) = c_{0} + c_{1}[cos(\phi)] + c_{2}[cos(2\phi)] + c_{3}[cos(3\phi)] = (c_{0} - c_{1} - c_{2} - c_{3}) + c_{1}[1 + cos(\phi)] + c_{2}[1 - cos(2\phi)] + c_{3}[1 + cos(3\phi)] \end{equation*}*The table below gives the parameters for the first form above.
cc5 $c_{0}$ originial 31394
cc6 $c_{0}$ original 5073
Type | $ \frac{c_{0}}{k_{B}} [K] $ | $ \frac{c_{1}}{k_{B}} [K] $ | $ \frac{c_{2}}{k_{B}} [K] $ | $ \frac{c_{3}}{k_{B}} [K] $ |
---|---|---|---|---|
CHx-(CH2)-(CH2)-CHy | 0 | 355.03 | -68.19 | 791.32 |
Cc5-(Cc5)-(Cc5)-Cc5 | -32534 | 45914 | 16518 | 1496 |
Cc6-(Cc6)-(Cc6)-Cc6 | -5339 | 6840 | 3509 | 63 |
Parameter Symbol | Name | Units | Special notes |
---|---|---|---|
$ Length $ | Bond Length | Angstrom | |
$ k_{\theta} $ | Harmonic Angle Constant | kcal/mol/radian$^2$ | Appears as $k_{\theta}/2$ in functional |
$ \theta_{eq} $ | Equilibrium Angle | degrees | |
$ c_{n} $ | Torsion Barrier | K | Given as $c_{n}/k_{B}$ |
$ {\epsilon} $ | well depth | K | Given as ${\epsilon}/k_{B}$ above |
$ {\sigma} $ | K | Given as ${\sigma}/k_{B}$ above |
A good resource : http://alma.karlov.mff.cuni.cz/bio/99_Studenti/00_Dalsi/ParamFit/2013_ParamFit_AmberTools13.pdf
The functional form of the Amber force field considers $ A_{ij} $ and $B_{ij} $ as the functional form, these are put into tleap using an frcmod with parameters $\sigma$ and $R_{min}$. A and B parameters are then calculated using the following combination rules:
\begin{equation*} \ A = {\epsilon}{R_{min}}^{12} \ B = 2{\epsilon}{R_{min}}^{6} \end{equation*}where
\begin{equation*} \ R_{min} = 2^{1/6}{\sigma} \end{equation*}and mixing rules are the same as outlined above for the TraPPE forcefield (Lorentz/Berthelot mixing rules)
Amber is able to handle dihedrals with more than one torsional term, as is needed for TraPPE.
For Trappe, there are multiple $c_{n}$'s corresponding to $V_{n}$'s in Amber. All TraPPE terms for straight alkanes have a phase $({\gamma}$) of zero.
For the case of $c_{0}$, n = 0, $cos(1)=0$, so insert phase shift of 90 degrees
\begin{equation*} \ U_{dihedral} = V_{n}[1 + cos(n{\phi}-{\gamma})] \end{equation*}The minus sign in the second term can be accounted for by incorporating a phase shift, ${\gamma}$ of 180 degrees into a multi- term Amber equation, I.E. -
\begin{equation*} c_{0} + c_{1}[1 + cos(\phi)] + c_{2}[1 - cos(2\phi)] + c_{3}[1 + cos(3\phi)] = V_{0}[1 + cos(0)] + V_{1}[1+cos({\phi})] + V_{2}[1+cos(2{\phi} - 90)] + V_{3}[1+cos(3{\phi})] \end{equation*}Parameter Symbol | Name | Units | Special notes |
---|---|---|---|
$ k_{b} $ | Harmonic Bond Constant | kcal/mol/Angstrom$^2 $ | |
$ r_{0} $ | Equilibrium Bond Length | Angstrom | |
$ k_{\theta} $ | Harmonic Angle Constant | kcal/mol/radian$^2$ | |
$ \theta_{0} $ | Equilibrium Angle | degrees | |
$ V_{n} $ | Torsion Barrier | kcal/mol | |
$ \gamma $ | Torision Phase | degrees | |
$ n $ | Torsion Periodicity | ||
$ A_{ij} $ | Nonbonded A parameter | Used in prmtop file (generated by tleap) | |
$ B_{ij} $ | Nonbonded B parameter | Used in prmtop file (generated by tleap) | |
$ R_{min}$ | vdW Radius | Angstrom | Input by user (frcmod file) |
$ {\epsilon} $ | well depth | kcal/mol | Input by user (frcmd file) |
In [1]:
import scipy.constants as const
import pandas as pd
import os
import glob
import re
import math
In [2]:
# Functions to convert units
def K_to_kcal_mol(joule_value):
'''
Convert value which is in units of K/kB to units of kcal/mol
'''
return const.Avogadro*(joule_value*const.k/(const.calorie*1000))
def convert_nb(trappe_epsilon, trappe_sigma):
'''
Convert nb parameters (epsilon and sigma) from Trappe FF to Amber FF.
Input -
Epsilon units - K
Sigma units - Angstrom
Output -
Rmin units Angstrom
epsilon units kcal/mol
'''
r_min = 2.**(1./6.)*trappe_sigma/2
epsilon = K_to_kcal_mol(trappe_epsilon)
return r_min, epsilon
def calculate_bead_masses(bead_name):
masses = {
"O" : 15.999,
"CH2" : 14.02658,
"CH3" : 15.0452,
"CH" : 13.01864,
"CH4" : 16.0425,
"H" : 1.00794,
}
try:
return [masses[x] for x in bead_name]
except KeyError:
print("Entry %s not found in mass dictionary" %(bead_name))
In [3]:
cyclic = ["cyclopentane" , "cyclohexane", "1-3-5-trioxane"]
def process_name_string(types, length):
"""
Processes names for atoms in bond and angle interactions to fit AMBER specifications
"""
return_names = []
name_list = [re.sub("\(|\)", "", x ) for x in types.values]
name_split = [[re.split("-| ", x) for x in name_list][n] for n in range(0,len(name_list))]
for n in range(0, len(name_split)):
name_split[n] = ["X " if "x" in x else x for x in name_split[n][0:length]]
name_split[n] = ["X " if "y" in x else x for x in name_split[n]]
name_split[n] = [s[0:3:2] if len(s) > 2 else s for s in name_split[n]]
return_names.append('-'.join(name_split[n]))
return return_names
def read_trappe_csv(f, name):
# Read data to pandas dataframe
data = pd.read_csv(f, names=range(0,20), skiprows=1)
data.dropna(inplace=True, axis=1, how='all')
# Find section starts
atom_section_start = data[data[1]=="(pseudo)atom"].index.tolist()
bond_section_start = data[data[1]=="stretch"].index.tolist()
angle_section_start = data[data[1]=="bend"].index.tolist()
torsion_section_start = data[data[1]=="torsion"].index.tolist()
# Start dictionary of relevant sections
sections = {}
## Create dataframes from sections
# Read in atom information
atoms = data.iloc[atom_section_start[0]+1:bond_section_start[0]].copy()
atoms.columns = data.iloc[atom_section_start[0]]
atoms.dropna(inplace=True, axis=1)
sections['atoms'] = atoms
# Get bond data from file. All bonds are set to the same k value for simplicity
test_val = float(data.iloc[bond_section_start[0]+1][0])
if not math.isnan(test_val):
bonds = data.iloc[bond_section_start[0]+1:angle_section_start[0]].copy()
bonds.columns = data.iloc[bond_section_start[0]]
bonds['k_b'] = 300.9
bonds.dropna(inplace=True, axis=1)
sections['bonds'] = bonds
# Get bond information
test_val = float(data.iloc[angle_section_start[0]+1][0])
if not math.isnan(test_val):
if len(torsion_section_start) > 0:
angles = data.iloc[angle_section_start[0]+1:torsion_section_start[0]].copy()
else:
angles = data.iloc[angle_section_start[0]+1:].copy()
angles.dropna(inplace=True, axis=1)
angles.columns = data.iloc[angle_section_start[0]][:5].values
sections['angles'] = angles
# Handle torsions
if len(torsion_section_start) > 0:
torsions = data.iloc[torsion_section_start[0]+1:].copy()
torsions.dropna(inplace=True, axis=1)
torsions.columns = data.iloc[torsion_section_start[0]].values
torsion_names = process_name_string(torsions['type'], 4)
torsions['type'] = torsion_names
if name in cyclic:
hold = torsions['c0/kB [K]'].astype(float) - torsions['c1/kB [K]'].astype(float) - \
torsions['c2/kB [K]'].astype(float) - torsions['c3/kB [K]'].astype(float)
torsions['c0/kB [K]'] = hold
print(name)
sections['torsions'] = torsions
return sections
def calc_amber_parameters(section_data, name):
## Process stored data frames
# Make write list dictionary
write_list = {}
for k,v in section_data.items():
if k == 'atoms':
# Calculate atom masses
v['masses'] = calculate_bead_masses(v['(pseudo)atom'].values)
# Manipulate atom names for atoms frame
v['(pseudo)atom'] = (v['(pseudo)atom'].str[0:3:2])
# Calculate epsilon in kcal/mol and rmin
v['r_min'], v['epsilon [kcal/mol]'] = convert_nb(v['epsilon/kB [K]'].astype(float),
v['sigma [Ang.]'].astype(float))
if k == 'angles':
# Process angles
v["k_theta [kcal/mol]"] = K_to_kcal_mol(v["k_theta/kB [K/rad^2]"].astype(float)/2.)
v['type'] = v['bend']
# Process angle type labels
v.drop(["bend"], axis=1, inplace=True)
if k == 'bonds':
v['type'] = v['stretch']
if k == 'torsions':
# Process torsions
v['ee_scaling'] = "SCEE=0.0"
v['nb_scaling'] = "SCNB=0.0"
v['type'] = v['torsion']
if name in cyclic:
v['phase_c0'] = 90
for x in range(3,-1,-1):
v['c%s/kB [kcal/mol]' %(x)] = K_to_kcal_mol(v['c%s/kB [K]' %(x)].astype(float))
v['%s' %(x)] = -x
if name not in cyclic:
if x == 2:
v['phase_c%s' %(x)] = -180
else:
v['phase_c%s' %(x)] = 0
elif x is not 0:
v['phase_c%s' %(x)] = 0
v['divider'] = 1
# Change type to atom names
types = v['type']
types = [re.sub("'|'| ","", x) for x in types]
proc_type = [re.split("-", x) for x in types]
if proc_type[0][0].isdigit():
name_list = []
for x in proc_type:
name_string = ""
atom_list = [section_data['atoms']['(pseudo)atom'][(section_data['atoms']["#"] == y)].values[0] for y in x]
name_string = "-".join(atom_list)
name_list.append(name_string)
v['type'] = name_list
# Remove unnecessary colums (unconverted units)
drop_list = [col for col in v.columns if '[K' in col]
v.drop(drop_list, inplace=True, axis=1)
if k is not 'atoms' and k is not 'torsions':
v.drop_duplicates(subset=["type"], inplace=True, keep='first')
elif k is 'torsions':
v.drop(['torsion',"#"], inplace=True, axis=1)
v.drop_duplicates(inplace=True, keep='first')
section_data['atoms'].drop_duplicates(subset=['type'], inplace=True, keep='first')
return section_data
def write_frcmod(name, amber):
# Get filename for frcmod & open
frcmod = open("frcmod_trappe.%s" %(name), "w")
frcmod.write("# Implementation of TraPPE FF for %s in Amber\n\n" %(name))
for k,v in amber.items():
if k == 'atoms':
# Write atom and nonbond data to frcmod
frcmod.write('MASS\n')
v.to_csv(frcmod, header=False, index=False, columns=['(pseudo)atom', 'masses'], sep="\t")
frcmod.write('\nNONBOND\n')
v.to_csv(frcmod, header=False, index=False, columns=["(pseudo)atom", "r_min",
"epsilon [kcal/mol]"], sep="\t")
if k == 'bonds':
# Write bond data to frcmod
frcmod.write('\nBONDS\n')
v.to_csv(frcmod, header=False,index=False, columns=["type", "k_b", "length [Ang.]"], sep="\t")
if k == 'angles':
frcmod.write('\nANGLE\n')
v.to_csv(frcmod, header=False,index=False, columns=["type", "k_theta [kcal/mol]",
"theta [degrees]"], sep="\t")
if k == 'torsions':
frcmod.write('\nDIHED\n')
for x in range(3,-1,-1):
v.to_csv(frcmod, header=False,index=False, columns=["type", "divider", 'c%s/kB [kcal/mol]' %(x), \
"phase_c%s" %(x),"%s" %(x), "ee_scaling", "nb_scaling"], sep="\t")
frcmod.close()
return 0
def write_leap_in(name, amber):
file_name = name + "_leap.in"
f = open(file_name, "w")
f.write("loadamberparams frcmod_trappe.%s\n" %(name))
f.write("sys = loadmol2 %s.mol2\n" %(name))
f.write('setbox sys "centers" 50\n')
f.write("saveamberparm sys trappe_%s_single_molecule.prmtop trappe_%s_single_molecule.inpcrd\n" %(name, name))
f.write("quit\n")
f.close()
In [4]:
files = glob.glob("trappe_*.csv")
run_file = "make_systems.sh"
rf = open(run_file, "w")
for f in files:
# Get molecule name
name = re.split("_|\.",f)[1]
print(name)
# Add to sh file
rf.write("tleap -f %s_leap.in\n" %(name))
# Process 'csv' from trappe
output = read_trappe_csv(f, name)
# Transform to amber format and units
amber = calc_amber_parameters(output, name)
# Print frcmod
write_frcmod(name, amber)
# Print accompanying leap input file
write_leap_in(name, amber)
rf.close()
In [ ]: