simulation.py

import random
CC_reqm = 4 #In angstrom units
CH_reqm = 3

CC_C12 = 2516582.4
CC_C6 = 1228.8

CH_C12 = 29108.222
CH_C6 = 79.857949

init_C = 259.679732
init_H = 313.5755
    

#Ground state energies for C and H atoms

class Atom(object):

    def __init__(self, name, valence):
        self.value = name #Name of Atom
        self.valence = valence #No. of valence electrons
        self.bonds = 0 #No. of bonds by default is 0

    
    def add(self, Atom1):
        #Atom1 to be added to current
        if Atom1.bonds+1<=Atom1.valence and 1+self.bonds<=self.valence:
            Atom1.bonds += 1
            self.bonds += 1
            return 1
        else:
            return 0 #Don't have enough valence electrons

def init(atom):
    if atom.value == 'C':
        return init_C
    else:
        return init_H

def energy(atom1, atom2):
    if (atom1.value=='C' and atom2.value=='C'):
        return (CC_C12/(CC_reqm**12) - CC_C6/(CC_reqm**6))
    else:
        return (CH_C12/(CH_reqm**12) - CH_C6/(CH_reqm**6))

##atoms = [Atom('C', 4), Atom('C', 4), Atom('H', 1), Atom('H', 1)]
atoms = list()
for i in range(0, 2):
    item = [Atom('C', 4), len(atoms)]
    atoms.append(item)

for i in range(0, 2):
    item = [Atom('H', 1), len(atoms)]
    atoms.append(item)

bonds = [[0, 0], [0, 0]]
mark = list() #List of atoms that are currently bonded.
energy_system = 0
energy_system = init_C #Start with a single C atom
mark.append(atoms[0])
del atoms[0]


while len(atoms)>0:
    idx = random.randint(0, len(atoms)-1)
    elem = atoms[idx] #The atom that is unmarked i.e. non-bonded
    while 1:
        idx_marked = random.randint(0, len(mark) - 1)
        print energy_system + energy(elem[0], mark[idx_marked][0])
        print init(mark[idx_marked][0])
        if elem[0].add(mark[idx_marked][0]) and (energy_system + energy(elem[0], mark[idx_marked][0])<init(mark[idx_marked][0])):
            mark.append(elem) #Add the atom to the marked pool
            print str(elem[1])+" "+str(mark[idx_marked][1])
            del atoms[idx]
            break
        else:
            print "ok"
            continue