.. _structures_example:
  
structures_example
========================
 

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

.. code:: ipython3

    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.

.. code:: ipython3

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

Let’s start with the Particle object

.. code:: ipython3

    from streamm.structures.particle import Particle

Create a particle object with label ‘C1’

.. code:: ipython3

    p_i = Particle(label='C1')

.. code:: ipython3

    print(p_i)

Assign the carbon element to the particle

.. code:: ipython3

    p_i.set_element('C')

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

.. code:: ipython3

    p_i.charge = -1.0

Check that the element properties were set to the particle

.. code:: ipython3

    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

.. code:: ipython3

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

Create a dictionary with new units

.. code:: ipython3

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

.. code:: ipython3

    p_i.update_units(new_unit_conf)

.. code:: ipython3

    print(p_i.show_attributes())

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

.. code:: ipython3

    p_i.update_units(default_unit_conf)

.. code:: ipython3

    print(p_i.show_attributes())

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

.. code:: ipython3

    p_j = Particle(symbol='H')

.. code:: ipython3

    print(p_j.show_attributes())

Let’s make an empty structure container

.. code:: ipython3

    from streamm.structures.structure import Structure

.. code:: ipython3

    mol_i = Structure('methane')

Now let’s construct a molecule

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

.. code:: ipython3

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

.. code:: ipython3

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

.. code:: ipython3

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

.. code:: ipython3

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

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

.. code:: ipython3

    import numpy as np
    import decimal

.. code:: ipython3

    bond_length = p_i.bonded_radius + p_j.bonded_radius 

.. code:: ipython3

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

.. code:: ipython3

    tet_a = bond_length/np.sqrt(3)

.. code:: ipython3

    print(tet_a)

Add hydrogens

.. code:: ipython3

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

.. code:: ipython3

    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

.. code:: ipython3

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

.. code:: ipython3

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

.. code:: ipython3

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

Check the position array

.. code:: ipython3

    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.

.. code:: ipython3

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

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

.. code:: ipython3

    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
            

.. code:: ipython3

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

.. code:: ipython3

    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.

.. code:: ipython3

    from streamm.structures.bond import Bond

based on the particle index values

.. code:: ipython3

    b_ij = Bond(0,1)

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

.. code:: ipython3

    mol_i.add_bond(b_ij)

.. code:: ipython3

    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)

.. code:: ipython3

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

.. code:: ipython3

    print(mol_i.bonded_nblist)

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

.. code:: ipython3

    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.

.. code:: ipython3

    print(mol_i.bonded_nblist.calc_nnab(0))

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

.. code:: ipython3

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

.. code:: ipython3

    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

.. code:: ipython3

    mol_json = mol_i.export_json()

.. code:: ipython3

    print(mol_i.tag)

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

.. code:: ipython3

    mol_test = Structure('methane')

.. code:: ipython3

    mol_test.import_json()

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

.. code:: ipython3

    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.

.. code:: ipython3

    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)

.. code:: ipython3

    import streamm.structures.group as group

.. code:: ipython3

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

Find groups based on residue variable

.. code:: ipython3

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

.. code:: ipython3

    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.

.. code:: ipython3

    res_json = groups_i.export_json()

Now let’s change the units

.. code:: ipython3

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

Check the positions

.. code:: ipython3

    print(mol_i.positions)

Check the particle bond radii

.. code:: ipython3

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