notebook.community

Edit and run

TRaPPE FF Parameter Conversion

TraPPE Force Field Parameters & Functional Form

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*}

Nonbonded Potential

\begin{equation*} \ u_{NB} = 4\epsilon_{ij}[(\frac{\sigma_{ij}}r_{ij})^{12} - (\frac{\sigma_{ij}}r_{ij})^6] \end{equation*}

Combination rules

\begin{equation*} \sigma_{ij} = \frac{1}2 (\sigma_{ii} + \sigma_{jj}) \end{equation*}\begin{equation*} \epsilon_{ij} = (\epsilon_{ii}\epsilon_{jj})^{1/2} \end{equation*}

Nonbonded Parameters

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

Bonded Potentials

Bond Stretching Parameters

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

Angles Bending Parameters

\begin{equation*} u_{bend} = \frac{k_{\theta}}2 (\theta - \theta_{eq} )^{2} \end{equation*}
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

Dihedral Potential

\begin{equation*} \ u_{torsion}(\phi) = c_{0} + c_{1}[1 + cos(\phi)] + c_{2}[1 - cos(2\phi)] + c_{3}[1 + cos(3\phi)] \end{equation*}\begin{equation*} \ u_{torsion_{cyclic}}(\phi) = c_{0} + c_{1}[cos(\phi)] + c_{2}[cos(2\phi)] + c_{3}[cos(3\phi)] \end{equation*}

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

FF Units

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

Amber FF Conversion

\begin{equation*} \ U_{total} = \sum_{bonds}{k_{b}(r-r_{0})^2} + \sum_{angles}{k_{\theta}{(\theta - \theta_o}^2)} + \sum_{dihedrals}{(V_{n}[1 + cos(n\phi - \gamma)]}) + \sum_{i=1}^{N-1}{\sum_{j=i+1}^{N}{ \frac{A_{ij}}{R_{ij}^{12}} - \frac{B_{ij}}{R_{ij}^{6}} }} \end{equation*}
  • Note - electrostatics omitted here for TRaPPE alkanes

A good resource : http://alma.karlov.mff.cuni.cz/bio/99_Studenti/00_Dalsi/ParamFit/2013_ParamFit_AmberTools13.pdf

Nonbonded Parameter Considerations

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)

Dihedral Considerations

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*}

Converting from TraPPE Dihedral form to Amber dihedral form

\begin{equation*} u_{dihedralTrappe}(\phi) = c_{0} + c_{1}[1 + cos(\phi)] + c_{2}[1 - cos(2\phi)] + c_{3}[1 + cos(3\phi)] \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*}

FF Units

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()



    


















    






tetrahydrofuran
cyclopentane
cyclopentane
cyclohexane
cyclohexane
ethanol
methanol
1-2-dimethoxyethane
propane
1-3-5-trioxane
1-3-5-trioxane
butane
ethane

















In [ ]:

Content source: dgasmith/EEX_scratch

Similar notebooks:

notebook.community | gallery | about