sim1.py

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

CC_C12 = 2516582.4
CC_C6 = 1228.8

CH_C12 = 29108.222 #kcal mol-1A12 http://www.csb.yale.edu/userguides/datamanip/autodock/html/Using_AutoDock_305.a.html
CH_C6 = 79.857949

init_C = 259.679732
init_H = 313.5755
    
R = 1.9872036*(1e-3) #https://en.wikipedia.org/wiki/Gas_constant
Na = 6.022140857*(1e23) #https://en.wikipedia.org/wiki/Avogadro_constant
T = 300 #Kelvin

acc_ratio = 0.5
def acceptance_ratio(delta):
    return math.exp((-delta*Na)/(R*T))
#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()

no_C = 2 #no of carbon
no_H = 2 #no_of_Hyderogen

for i in range(0, no_C):
    item = [Atom('C', 4), len(atoms)] #len(atoms) tracks index of the generated atom.
    atoms.append(item)

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

bonds = [[0 for x in range(0, no_C)] for y in range(0, no_H)]
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])
##        print acceptance_ratio(energy(elem[0], 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