structures_example

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

from pprint import pprint
import copy

Set up a log file for this example so we can read what exactly streamm is doing, if we feel like it.

import logging
logging.basicConfig(filename='structures_example.log',level=logging.INFO)

Let’s start with the Particle object

from streamm.structures.particle import Particle

Create a particle object with label ‘C1’

p_i = Particle(label='C1')

print(p_i)

Assign the carbon element to the particle

p_i.set_element('C')

Let’s oxidize the carbon just to make the charge non-zero

p_i.charge = -1.0

Check that the element properties were set to the particle

print(p_i.show_attributes())

Say we want to change the units to SI

Let’s look at the current units of the particle instance

default_unit_conf = copy.deepcopy(p_i.unit_conf)
pprint(default_unit_conf)

Create a dictionary with new units

new_unit_conf = {'length':'m','mass':'kg','charge':'C'}

p_i.update_units(new_unit_conf)

print(p_i.show_attributes())

That’s cool, but we should stick with the default units, so let’s change them back

p_i.update_units(default_unit_conf)

print(p_i.show_attributes())

Let’s create another particle and set the element to hydrogen

p_j = Particle(symbol='H')

print(p_j.show_attributes())

Let’s make an empty structure container

from streamm.structures.structure import Structure

mol_i = Structure('methane')

Now let’s construct a molecule

We can add the carbon at the origin using the add_partpos() function.

pos_i = [0.0,0.0,0.0]
mol_i.add_partpos(p_i,pos_i)

for p_index,particle_i in mol_i.particles.items():
    if( particle_i.symbol == 'H' ):
        particle_i.residue = 1

        h_cnt += 1

for p_index,particle_i in mol_i.particles.items():
    print(p_index,particle_i)

print("Now the structure container has {} particle ".format(mol_i.n_particles))

Find the positions of the hydrogens to give a tetrahedral molecular geometry

import numpy as np
import decimal

bond_length = p_i.bonded_radius + p_j.bonded_radius

print(bond_length,mol_i.unit_conf['length'])

tet_a = bond_length/np.sqrt(3)

print(tet_a)

Add hydrogens

pos_j = [tet_a,tet_a,tet_a]
mol_i.add_partpos(p_j,pos_j)

for p_index,particle_i in mol_i.particles.items():
    print(p_index,particle_i)

We can add the subsequent hydrogens using the same particle object since add_partpos makes a deepcopy of the object when adding to the structure container

pos_j = [-tet_a,-tet_a,tet_a]
mol_i.add_partpos(p_j,pos_j)

pos_j = [-tet_a,tet_a,-tet_a]
mol_i.add_partpos(p_j,pos_j)

pos_j = [tet_a,-tet_a,-tet_a]
mol_i.add_partpos(p_j,pos_j)

Check the position array

print(mol_i.positions)

The particles instance variable of the structure container is a dictionary, so we can just loop over that using the iteritems() function.

for p_index,particle_i in mol_i.particles.items():
    print(p_index,particle_i)

Hum, let’s fix the labels of the hydrogens…

h_cnt = 1
for p_index,particle_i in mol_i.particles.items():
    if( particle_i.symbol == 'H' ):
        particle_i.label = 'H{}'.format(h_cnt)

        h_cnt += 1

for p_index,particle_i in mol_i.particles.items():
    print(p_index,particle_i)

Okay, that looks better

Print .xyz file and check geometry with a molecular viewer such as Avogadro (https://avogadro.cc/)

mol_i.write_xyz()

Looks good, you should have the geometry of a methane molecule with a C-H bond length of 1.2 Angstroms

However, we have not told streamm about the bonds. There are a few ways to do this, let’s do it explicitly with the Bond object fist.

from streamm.structures.bond import Bond

based on the particle index values

b_ij = Bond(0,1)

Now add the bond to the bonds dictionary in the structure container

mol_i.add_bond(b_ij)

print("Now the structure container has {} particle/s and {} bond/s".format(mol_i.n_particles,mol_i.n_bonds))

Neat, but adding all the bonds, bond angles and dihedrals explicitly would be pretty tedious, so let’s use some functions to do that.

First, let’s guess the bonded_nblist of the molecule based on the bonded_radius of each particle (atom)

mol_i.bonded_nblist = mol_i.guess_nblist(0,radii_buffer=1.35)

print(mol_i.bonded_nblist)

Let’s take a look at the neighbor lists list and index instance variables

print(mol_i.bonded_nblist.list)
print(mol_i.bonded_nblist.index)

Looking at the index for particle 0, we get that it has neighbors in the list from 0:3 (index[0]:index[0+1]-1). Therefore, we know particle 0 has [1, 2, 3, 4] for neighbors.

print(mol_i.bonded_nblist.calc_nnab(0))

Now we can use the bonded neighbor list to construct the bonds, bond angles and dihedrals

mol_i.bonded_bonds()
mol_i.bonded_angles()
mol_i.bonded_dih()

print(mol_i.print_properties())

A little easier than adding everything by hand

We can dump a json file methane_struc.json for checkpointing the entire structure

mol_json = mol_i.export_json()

print(mol_i.tag)

We can reload the molecule as another instance of the Structure object

mol_test = Structure('methane')

mol_test.import_json()

Then we can use the structure’s eq and ne functions to make sure they are the same

if( mol_test != mol_i ):
    print("Molecules not equal")
elif( mol_test == mol_i ):
    print("Molecules equal")
else:
    print("Uh Oh")

Now let’s set some groups. This is a little unnecessary for methane, but it will come in super handy if you have a large simulation of thousands of molecules.

To do this we will set the residue variable for each particle.

mol_i.particles[0].residue = 0
for p_index,particle_i in mol_i.particles.items():
    if( particle_i.symbol == 'H' ):
        particle_i.residue = 1
    print(particle_i, particle_i.residue)

import streamm.structures.group as group

groups_i = group.Groups('methane_residues',mol_i)

Find groups based on residue variable

groups_i.group_prop('residue',groups_i.tag)

for g_index,group_i in groups_i.groups.items():
    print(group_i.pkeys)
    print(group_i.write_coord())

Looks good. We have two groups in the group container, the first with the carbon particle index 0 and the rest are the hydrogens.

res_json = groups_i.export_json()

Now let’s change the units

mol_i.update_units({'length':'pm'})

Check the positions

print(mol_i.positions)

Check the particle bond radii

for p_index,particle_i in mol_i.particles.items():
    print(particle_i,particle_i.bonded_radius)