forcefields_example

This notebook imports the fundamental objects of the streamm.forcefields module and goes through the functionality of each

from pprint import pprint

import logging
logging.basicConfig(filename='forcefield_example.log',level=logging.DEBUG)

from pathlib2 import Path
import os

import streamm.forcefields.particletype as particletype

Let’s start with an ethane molecule from the buildingblocks_example.ipynb example

We have an sp3 carbon bonded to hydrogens

Let’s create the force field parameters object for a ‘CT’ carbon and an ‘HC’ hydrogen

CT = particletype.Particletype('CT')
HC = particletype.Particletype('HC')

Set some parameters from J. Am. Chem. Soc., 1996, 118 (45), pp 11225–11236

https://pubs.acs.org/doi/abs/10.1021/ja9621760

In general you should pick a force field that has been shown to work well for your system and set up the parameters

Check that we have our units set right

print(CT.unit_conf['energy'],CT.unit_conf['length'])

Our potential is in kCal/mol (kCalmol) so let’s get the unit dictionary and create our own default

import copy

oplsaa_unit_conf = copy.deepcopy(CT.unit_conf)

oplsaa_unit_conf['energy'] = 'kCalmol'

pprint(oplsaa_unit_conf)

CT.update_units(oplsaa_unit_conf)

HC.update_units(oplsaa_unit_conf)

CT.epsilon = 0.066 # kcal/mol
CT.sigma = 3.5 # Angstroms

HC.epsilon = 0.03 # kcal/mol
HC.sigma = 2.5 # Angstroms

Set mass using periodic table

import pymatgen_core.core.periodic_table as periodic_table

CT.mass =  periodic_table.Element['C'].atomic_mass.real
HC.mass =  periodic_table.Element['H'].atomic_mass.real

Set the bond stretching parameters

import streamm.forcefields.bondtype as bondtype

C_H = bondtype.Bondtype('CT','HC',unit_conf=oplsaa_unit_conf)

C_H.setharmonic(1.080,367.0)

print(C_H)

C_C = bondtype.Bondtype('CT','CT',unit_conf=oplsaa_unit_conf)
C_C.setharmonic(1.53,268.0)

import streamm.forcefields.angletype as angletype

H_C_H = angletype.Angletype('HC','CT','HC',unit_conf=oplsaa_unit_conf)

H_C_H.setharmonic(110.7,37.50)

print(H_C_H)

H_C_C = angletype.Angletype('HC','CT','CT',unit_conf=oplsaa_unit_conf)
H_C_C.setharmonic(110.7,37.50)

Now we need a dihedral potential for the HC-CT-CT-HC dihedral

import streamm.forcefields.dihtype as dihtype

H_C_C_H = dihtype.Dihtype('HC','CT','CT','HC',unit_conf=oplsaa_unit_conf)

H_C_C_H.type ='opls'

H_C_C_H.setopls(0.0,0.0,0.3,0.0)

Let’s create a parameter container to keep track of our parameters

import streamm.forcefields.parameters as parameters

paramC = parameters.Parameters('oplsaa',unit_conf=oplsaa_unit_conf)

Add parameters to the container

paramC.add_particletype(CT)

paramC.add_particletype(HC)

paramC.add_bondtype(C_H)
paramC.add_bondtype(C_C)

paramC.add_angletype(H_C_H)
paramC.add_angletype(H_C_C)

paramC.add_dihtype(H_C_C_H)

print(paramC)

for ptkey,pt in paramC.particletypes.items():
    print(ptkey,pt,pt.unit_conf['energy'],pt.unit_conf['length'])

for btkey,bt in paramC.bondtypes.items():
    print(btkey,bt,bt.unit_conf['harm_bond_coeff'],pt.unit_conf['length'])

for atkey,at in paramC.angletypes.items():
    print(atkey,at,at.unit_conf['energy'],at.unit_conf['length'])

print(paramC.tag)

paramC.unit_conf

We can dump a pickle file

paramC.dump_pickle()

Or we can export a json object

paramC_json = paramC.export_json()

Read in ethane .json file from the structures example

import streamm.structures.buildingblock as bb

need_files = ['ethane_struc.json']
for f in need_files:
    path = Path(f)
    if not path.is_file():
        print("Need to run buildingblocks_example.ipynb")
        os.system("jupyter nbconvert --to python  buildingblocks_example.ipynb")
        os.system("python buildingblocks_example.py")

mol_i = bb.Buildingblock('ethane')

mol_i.import_json()

print(mol_i.print_properties())

Let’s set the paramkey for each particle based on the symbol.

for pk,p in mol_i.particles.items():
    print(p.symbol)
    if( p.symbol == 'C' ):
        p.paramkey = 'CA'
    elif( p.symbol == 'H' ):
        p.paramkey = 'HA'
    print(p.paramkey ,mol_i.bonded_nblist.calc_nnab(pk))

This is a bit redundant, but we can think of a more complex molecule where we could use the number of neighbors to write a more complex routine

print(mol_i.n_particles)

Now we can set the particles, bonds, bond angles and dihedrals of the molecule to have parameters

First lets set the particle types

for pk,p in mol_i.particles.items():
    if( p.paramkey == 'CA' ):
        p.param = CT
        p.param_index = 0
    elif( p.paramkey == 'HA' ):
        p.param = HC
        p.param_index = 1

Now we can set the bond types

for bk,b in mol_i.bonds.items():
    b.param = C_H
    b.param_index = 0

for ak,a in mol_i.angles.items():
    a.param = H_C_H
    b.param_index = 0

for dk,d in mol_i.dihedrals.items():
    d.param = H_C_C_H
    d.param_index = 0

print("Particles ")
for pk,p in mol_i.particles.items():
    print(p,p.param, p.param_index)
print("\n Bonds ")
for bk,b in mol_i.bonds.items():
    print(b,b.param, b.param_index)
print("\n Bond angles ")
for ak,a in mol_i.angles.items():
    print(a,a.param, a.param_index)
print("\n Dihedrals ")
for ak,a in mol_i.dihedrals.items():
    print(a,a.param, a.param_index)

Now our molecule has forcefield parameters for all the interactions

Now let’s say we want to use a software package like GROMACS that used kJ/mol instead of kCal/mol

gromacs_unit_conf = copy.deepcopy(oplsaa_unit_conf)

gromacs_unit_conf['energy'] = 'kJmol'
gromacs_unit_conf['length'] = 'nm'

gromacs_unit_conf['harm_bond_coeff'] = 'kJmolsqnm' #*

• The harmonic bond coefficient harm_bond_coeff has to be changed as well since it has special units of energy/length^2

pprint(gromacs_unit_conf)

mol_i.update_units(gromacs_unit_conf)

print("Particles ")
for pk,p in mol_i.particles.items():
    print(p,p.param, p.param_index)
print("\n Bonds ")
for bk,b in mol_i.bonds.items():
    print(b,b.param, b.param_index)
print("\n Bond angles ")
for ak,a in mol_i.angles.items():
    print(a,a.param, a.param_index)
print("\n Dihedrals ")
for ak,a in mol_i.dihedrals.items():
    print(a,a.param, a.param_index)

mol_i.update_units(oplsaa_unit_conf)

print("Particles ")
for pk,p in mol_i.particles.items():
    print(p,p.param, p.param_index)
print("\n Bonds ")
for bk,b in mol_i.bonds.items():
    print(b,b.param, b.param_index)
print("\n Bond angles ")
for ak,a in mol_i.angles.items():
    print(a,a.param, a.param_index)
print("\n Dihedrals ")
for ak,a in mol_i.dihedrals.items():
    print(a,a.param, a.param_index)