Untitled

#1/9/25 EC Part A Stefan Botzolakis
#Bonding Morse Potential by Tom Collins

scene.height = 250 #sets the size of the animation window
gd = gdisplay(x=0, y=0, width=600, height=250, xtitle='distance (A)', ytitle='PE (kJ/mol)', foreground=color.black, background=color.white, xmax=2 * BL, xmin= BL/2, ymax=BE, ymin=- (BE+20))
BM = 1 #bond multiplicity, 1= single etc.

#---------Bond 1 Config---------------

BL = 0.74 #eqm bond length in Angstroms
BE = 432 #bond enthalpy in kjmol-1
k = sqrt(BM/2*BE) #sets the spring constant of the bond
M1 = 1 #mass of atom 1 in amu
M2 = 1 #mass of atom 2 in amu

#-------Bond 2 Config---------------------

M3 = 2
M4 = 2
BL2 = 1
BE2 = 300
k2 = sqrt(BM/2*BE2) 

max_PE = 0
max_PE2 = 0

#-------------Bond 2 Logic ------------------
nuc3 = sphere(pos=vector(-(BL2 / 2 + 0.5), 0.5, 0), radius=0.2 * M3 ** (1 / 3), color=color.blue)
nuc4 = sphere(pos=vector((BL2 / 2 + 0.5), 0.5, 0), radius=0.2 * M4 ** (1 / 3), color=color.yellow)
nuc3.eqm = vector(-BL2 / 2, 0.5, 0)
nuc4.eqm = vector(BL2 / 2, 0.5, 0)
#print("test")
bond2=helix(pos=nuc3.pos, axis=nuc4.pos - nuc3.pos, radius=0.1) 
nuc3.velocity = vector(0, 0, 0)
nuc4.velocity = vector(0, 0, 0)
nuc3.momentum = M3 * nuc3.velocity
nuc4.momentum = M4 * nuc4.velocity
f2 = gcurve(color=color.red, dot=True) #graph
#print("test")

#-----------------Bond 1 Origional Logic
nuc1=sphere(pos = vector(-(BL/2 + 0.2),0,0), radius = 0.2 * M1**(1/3), color = color.red)#sets nucleus 1 conditions, starting at 0.2 A further than eqm
nuc2=sphere(pos = vector((BL/2 + 0.2),0,0), radius = 0.2 * M2**(1/3), color = color.green)#sets nucleus 2 conditions, starting at 0.2 A further than eqm
nuc1.eqm=vector(-BL/2,0,0) #defines eqm position for nuc1
nuc2.eqm=vector(BL/2,0,0) #defines eqm position for nuc2
r_nuc_nuc= nuc2.pos-nuc1.pos #calculates the vector between the two nuclei
bond=helix (pos = nuc1.pos, axis = r_nuc_nuc, radius = 0.1) # shows the bond as a spring between the two nuclei
nuc1.velocity=vector(0,0,0) #gives nuc 1 an initial velocity of 0
nuc2.velocity=vector(0,0,0) #gives nuc 2 an initial velocity of 0
nuc1.montumme=M1*nuc1.velocity #defines the momentum of nuc 1
nuc2.momentum=M2*nuc2.velocity #defines the momentum of nuc 2
f1 = gcurve(color=color.cyan, dot = True) #graph

#--------------Simulation Loop (both bonds)---------

dt = 0.015 #sets the steps in time that the code recalculates
t = 0 #starts the "clock" at t = 0
while t<=2: 
  rate(300) #slows the animation to 100 calculations per real second
  nuc1.pos = nuc1.pos + nuc1.velocity*dt #tells the code to recalculate the position of the nucleus after every step
  nuc2.pos = nuc2.pos + nuc2.velocity*dt #tells the code to recalculate the position of the nucleus after every step
  r_nuc_nuc = nuc2.pos-nuc1.pos #recalculates the vector between nuclei
  bond.pos=nuc1.pos #moves the spring for each step
  bond.axis=r_nuc_nuc #changes the length of the spring for each step
  F1 = -k * (nuc1.pos-nuc1.eqm) #calculates the force on atom 1
  F2 = -k * (nuc2.pos-nuc2.eqm) #calculates the force on atom 2
  nuc1.momentum = nuc1.momentum + F1*dt #calculates the momentum on atom 1 after the force has acted
  nuc2.momentum = nuc2.momentum + F2*dt #calculates the momentum on atom 2 after the force has acted
  nuc1.velocity = nuc1.momentum/M1 #calculates the new velocity of atom 1
  nuc2.velocity = nuc2.momentum/M2 #calculates the new velocity of atom 2
  scene.center=nuc1.pos + 0.5*bond.axis # positions the center of the scene in the center of the bond
  scene.autoscale=False #stops the animation from rescaling after starting
  PE = BE * (exp(-2*k*(r_nuc_nuc.mag-BL)) - 2*exp(-k*(r_nuc_nuc.mag-BL))) #uses the Morse equation to calculate the PE of the system in kJmol-1
  f1.plot(pos=(r_nuc_nuc.mag,PE))
  
  
  nuc3.pos = nuc3.pos + nuc3.velocity * dt
  nuc4.pos = nuc4.pos + nuc4.velocity * dt
  r_nuc3_nuc4 = nuc4.pos - nuc3.pos
  bond2.pos = nuc3.pos
  bond2.axis = r_nuc3_nuc4
  F3 = -k2 * (nuc3.pos - nuc3.eqm)
  F4 = -k2 * (nuc4.pos - nuc4.eqm)
  nuc3.momentum += F3 * dt
  nuc4.momentum += F4 * dt
  nuc3.velocity = nuc3.momentum / M3
  nuc4.velocity = nuc4.momentum / M4
  PE2 = BE2 * (exp(-2 * k2 * (r_nuc3_nuc4.mag - BL2)) - 2 * exp(-k2 * (r_nuc3_nuc4.mag - BL2)))
  f2.plot(pos=(r_nuc3_nuc4.mag, PE2))

  if PE > max_PE:
    max_PE = PE
    #print("test")

  if PE2 > max_PE2:
    max_PE2 = PE2
    
  t = t + dt #updates the time after each step
  
#print("test")
print(f"The maximum PE for bond 1 was: {max_PE:.2f} kJ/mol")
print(f"The maximum PE for bond 2 was: {max_PE2:.2f} kJ/mol")