Neuroscience based machine intelligence models by Gary Gaulin

Wednesday, October 25, 2023

Create Vesicle using Lennard Jones Potential, and propagate Traveling Waves over surface

import sys
import pygame
import random
import numpy as np
import os

program_dir = os.getcwd() # Get the Folder path this Python Programm is running from.
image_dir = program_dir + '\Images' # Directory to PNG Image files
#if not os.path.exists(image_dir): # If no directory folder for PNG files yet
# os.mkdir('Images') # then make new directory.

signalcollapsepoint = 0 # Program start point of Wave, oint where wave collapsed on other side.
signalactivity = 0 # Total number of connections that are spiking/on/set/1 during this time cycle.
defaultpointcount = 1024 # Number of particles when no StartupXYZ.txt file. Symmetrical are 32, 360.
pressureRampStep = .0025
pressureMin = 1
pressureMax = 1
pressureRampDirection = 1
adjustR = 1.01 # Shaking pressure like in jar of water, can help jog into place.
forceCutOff = 2 # Amount of spacing distance at which forces are too small to influence.
groupCutOff = 3.5
connectCutOff = 1.3 # Amount of spacing between neighbors to connect to.
rotationstep = .0125 # Amount of rotation per timestep.
rotation = 0 # Initialize rotation angle, to draw all points on screen.
epsilon = 1 # Lennard Jones variable, remains 1, normalized.
sigma = 1 # Lennard Jones variable, remains 1, normalized.
dt = .025 # Timestep size, how (slower but) precise it moves.
maxf = 2 # Velocity limit for when enough force to fly off sphere.
popwaitcount = 400
popwait = 0
pops = 0

height = 520 # Size in pixels for each of 2 squares to draw inside.
screensize = (height*2, height) # Pygame variable used to initialize screen window.
screenradius = int(height*.43) # Radius of vesicle in pixels, when drawn to screen.
xc1 = int(screensize[0]*.25) # Center x screen location for left square.
xc2 = int(screensize[0]*.75) # Center x screen location for square to right.
yc = int(screensize[1]*.5) # Center y screen location point for both squares.
timestepnumber = 0 # Counter for number of times program was stepped/looped one moment in time.

pi = 3.14159265358979323846264338327950288412
radian = pi * 2
threeQtrRadian = radian * 3/4

Black = (0, 0, 0)
Brown = (55, 0, 0)
Red = (255, 0, 0)
Orange = (255, 165, 0)
Yellow = (220, 220, 0)
Green = (0, 200, 0)
Blue = (0, 0, 150)
Violet = (255, 0, 255)
Gray = (120, 120, 120)
White = (255, 255, 255)
Cyan = (0, 255, 255)

# Colors used in electronics to mark resistors and other components with bands, to make easy to read not tiny numbers.
# These are used to similarly identify number of connections each cell has by color, to make easy to see what they are.
# 0=Black, 1=Brown, 2=Red, 3=Orange, 4=Yellow, 5=Green, 6=Blue, 7=Violet, 8=Gray, 9=White
ElectronicColorCode = Black, Brown, Red, Orange, Yellow, Green, Blue, Violet, Gray, White
# Colors used to draw connections that are either 0, 1 = action potential.
InputSignalColor = Black, Cyan # Colors for connections that just propagated the bit to next point/cell.
OutputSignalColor = Black, Red # Full on for entire timecycle is Red.


def csvwrite(sourcelist, filename): # Write Comma Separated Value list of lists [,,] or tuples (,,) to disk.
# For control over file behavior csv or other import module is not used for reading or writing.
with open(filename, 'w') as f: # Open file to write over, using 'f' as its name:
for line in sourcelist: # For each line of values in source list:
fs = str(line) # Convert line to a string containing commas.
f.write(fs[1:-1] + "\n") # Write all but [] or () chars, line feed at end.


# File to read must have no blank lines or comments added to it, comma separated values only.
def xyzread(filename): # Read disk file back into list of lists.
x = []
y = []
z = []
try:
with open(filename, 'r') as f: # Open file in read mode using 'f' as its name:
for line in f: # For each comma separated line of code in 'f':
fs = line.rstrip('\n') # Remove line x,y,z readings in line into s
if len(fs) > 0:
fs = fs.split(',') # Split line at commas into list named s
x.append(float(fs[0])) # Save first value in s as x reading.
y.append(float(fs[1])) # Save second value in s as y reading.
z.append(float(fs[2])) # Save third value in s as z reading.
except IOError: # XYZ File does not exist. Start new vesicle.
x = [0.0 for _ in range(defaultpointcount)]
y = [0.0 for _ in range(defaultpointcount)]
z = [0.0 for _ in range(defaultpointcount)]
# Randomly distribute points inside sphere by filling a cube but only save points inside.
for i in range(defaultpointcount): # For each x,y,z reading to save in list.
ra = screenradius + 1
while ra > screenradius: # Loop until radius is inside sphere.
x[i] = random.uniform(-screenradius, screenradius)
y[i] = random.uniform(-screenradius, screenradius)
z[i] = random.uniform(-screenradius, screenradius)
az, el, ra = cart2sphere(x[i], y[i], z[i]) # Generate random Cartesian point inside a cube.
return x, y, z


# Write Comma Separated Value list of lists [,,] or tuples (,,) to disk.
def xyzwrite(filename):
# For control over file behavior csv or other import module is not used for reading or writing.
with open(filename, 'w') as f: # Open file to write over, using 'f' as its name:
for i in range(points):
sxyz = str([xp[i], yp[i], zp[i]]) # For each value in the 3 source lists
f.write(sxyz[1:-1] + "\n") # Write all but [] or () chars, line feed at end.
screen.fill(White) # Make screen background flash white to indicate file save.


# File to read must have no blank lines or comments added to it, comma separated values only.
def csvread(filename): # Read disk file back into list of lists.
destination = [] # Clear destination list lists are appended to.
with open(filename, 'r') as f: # Open file in read mode using 'f' as its name:
for line in f: # For each comma separated line of code in 'f':
fs = line.rstrip('\n')
d = [] # Empty list named 'd' to fill with elements.
if len(fs) > 0:
fs = fs.split(',') # Split line at commas into list named 's'
for e in fs: # For each string element listed in 's' list:
d.append(float(e)) # Append floating point element value to list.
destination.append(d) # Append list of elements to destination list.
return destination # Returns list of lists only = [[,,],[,,],[,,]]


# Calculate spherical angle force vectors of a point are together most pointing towards.
# Magnitude becomes 0 when neighboring points are in opposite or opposing each other.
# icenter is the point to calculate.
def vectorsum(icenter, listname):
xsum = 0.0
ysum = 0.0
zsum = 0.0
for i in listname:
xctr = xp[icenter]
yctr = yp[icenter]
zctr = zp[icenter]
xsum += xp[i] - xctr
ysum += yp[i] - yctr
zsum += zp[i] - zctr
az, el, ra = cart2sphere(xsum, ysum, zsum)
return az, el, ra


# Cartesian to Spherical coordinates.
def cart2sphere(x, y, z):
hxy = np.hypot(x, y)
return np.arctan2(y, x), np.arctan2(z, hxy), np.hypot(hxy, z)


# Spherical to Cartesian coordinates.
def sphere2cart(az, el, ra):
theta = ra * np.cos(el)
x = theta * np.cos(az)
y = theta * np.sin(az)
z = ra * np.sin(el)
return x, y, z


# To make program faster: Return a list of nearbby points less than (too far to attract) CutOff distance away.
def update_grouplist():
print(timestepnumber, 'Compiling Group list')
groupsum = 0
for i in range(points): # For all the points in array
grouplist[i] = []
for j in range(points): # For all the other points in array
if i != j: # Except itself
dx = xp[j] - xp[i] # Calculate distance between points along x.
if abs(dx) <= groupR: # If (neg made positive) x distance OK continue.
dy = yp[j] - yp[i] # Calculate distance between points along y.
if abs(dy) <= groupR: # If (neg made positive) y distance OK continue.
dz = zp[j] - zp[i] # Calculate distance between points along z.
if abs(dz) <= groupR: # If (neg=pos.) z OK then time consuming sqrt precision.
dist = np.sqrt((dx * dx) + (dy * dy) + (dz * dz)) # Calculate 3D distance.
if dist < groupR: # If distance is within group range then.
grouplist[i].append(j) # Append i and j points to the end of list.
groupsum += len(grouplist[i])
print("Updated Group list, ", len(grouplist), "total distances, average", groupsum / points, "per point.")


# Return a list of points within cutoff range.
def update_cutofflist():
print(timestepnumber, 'Compiling CutOff list')
cutoffsum = 0
for i in range(len(grouplist)):
cutofflist[i] = []
for j in grouplist[i]:
if i != j: # Except itself
dx = xp[j] - xp[i] # Calculate distance between points along x.
if abs(dx) <= cutoffR: # If (neg made positive) x distance OK continue.
dy = yp[j] - yp[i] # Calculate distance between points along y.
if abs(dy) <= cutoffR: # If (neg made positive) y distance OK continue.
dz = zp[j] - zp[i] # Calculate distance between points along z.
if abs(dz) <= cutoffR: # If (neg=pos.) z OK then time consuming sqrt precision.
dist = np.sqrt((dx * dx) + (dy * dy) + (dz * dz)) # Calculate distance between.
if dist < cutoffR: # If distance is within cutoff range then
cutofflist[i].append(j) # Append i and j points to the end of list.
cutoffsum += len(cutofflist[i])
print("Updated CutOff list, ", len(cutofflist), "total distances, average", cutoffsum/points, "per point.")
return cutoffsum/points


# Return a list of particles within group range, and another list for close enough to form in/out connections to.
def update_connections():
global connections
for i in range(points):
connectlist[i] = []
for j in cutofflist[i]:
if i != j:
dx = xp[j] - xp[i] # Calculate x,y,z distances.
dy = yp[j] - yp[i]
dz = zp[j] - zp[i]
dist = np.sqrt((dx * dx) + (dy * dy) + (dz * dz)) # Calculate 3D distance, a time consuming command.
if dist < connectR: # If close enough to connect together then
connectlist[i].append(j) # save j point at this i, makes list of what connects to where.
connections[i] = len(connectlist[i])


# Update (azimuth) longitude and (elevation) latitude coordinates. Radius for all points is that of sphere, not needed.
def update_angles():
global latlong
latlong = []
for i in range(points):
az, el, ra = cart2sphere(xp[i], yp[i], zp[i]) # Spherical azimuth,elevation = longitude,latitude.
latlong.append([-el, az]) # Elevation is negative, north pole is positive latitude.
# print(latlong)


def update_positions():
global cutofflist
if holdradius:
for i in range(points):
# Hold to sphere. Convert points to spherical coordinates then back using sphere radius as center distance.
xp[i], yp[i], zp[i] = sphere2cart(np.arctan2(yp[i], xp[i]), # Calc azimuth.
np.arctan2(zp[i], np.hypot(xp[i], yp[i])), # Calc elevation.
vesicleR) # Desired sphere radius.
# Initialize velocity vector arrays for summing forces, combined angles and magnitudes.
xi = [0 for _ in range(points)]
yi = [0 for _ in range(points)]
zi = [0 for _ in range(points)]
xj = [0 for _ in range(points)]
yj = [0 for _ in range(points)]
zj = [0 for _ in range(points)]
# Clear sum, minimum and maximum
ljt = 0 # Initialize variable for totaling Lennard-Jones force.
ljmaximum = -10000 # Initialize for finding max, will be made greatest reading.
ljmaximum_i = 0 # For also saving max 'i' number of greatest reading.
ljminimum = 100000000000 # Initialize for finding min, will be made greatest reading.
ljminimum_i = 0 # For also saving min 'i' number of greatest reading.
ljminimum_j = 0 # For also saving min 'j' number of greatest reading.
atmax = 0
for i in range(len(cutofflist)): # For all the points in array.
ljsum[i] = 0
ljt1 = 0
for j in cutofflist[i]:
# Calculate x,y,z distances.
dx = xp[j] - xp[i]
dy = yp[j] - yp[i]
dz = zp[j] - zp[i]
# Use the three distances for calculating an always positive 3D distance.
dist = np.sqrt((dx * dx) + (dy * dy) + (dz * dz))
# Calculate Lennard Jones force.
ljr = dist/spacing # Lennard-Jones radius. At r=1 force is 0, spacing=r*2
lj = 4 * epsilon * ((sigma/ljr)**12 - (sigma/ljr)**6)
# Reduce timestep motion for smooth movement.
f = lj * dt
if f > maxf: # Limit amount of displacement to not fly off screen.
f = maxf
atmax = 1
# f = np.sqrt((x * x) + (y * y) + (z * z))
ljt1 += lj # Add to sum for average of all points.
ljsum[i] += f # Used for brightness to draw on screen.
if lj > ljmaximum: # If force is greater than last greatest then:
ljmaximum = lj # Save maximum reading.
ljmaximum_i = i # Also save max 'i' number of greatest reading.
if lj < ljminimum: # If force is lower than last lowest then:
ljminimum = lj # Save minimum reading.
ljminimum_i = i # Also save min 'i' number of lowest reading.
ljminimum_j = j # Also save min 'j' number of lowest reading.
# if f < -maxf: # Limit amount of displacement to not fly off screen.
# f = -maxf
# Use distance for calculating spherical azimuth and elevation angle r, as in cart2sphere.
x, y, z = sphere2cart(np.arctan2(dy, dx), np.arctan2(dz, dist), f)
# Update x,y,z force acting on i point, from this j point.
xi[i] -= x
yi[i] -= y
zi[i] -= z
# Update other point, j. For every action there is equal and opposite reaction.
xj[j] += x
yj[j] += y
zj[j] += z
# ljt += ljt1 / len(cutofflist[i]) # Add to sum for average of all points.
listLen = len(cutofflist[i])
if listLen > 0: # List length must be greater than zero or division by zero error.
ljt += ljt1 / len(cutofflist[i]) # Add to sum for average of all points.
# Add force displacement amounts, to points.
for i in range(points):
xp[i] += xi[i] + xj[i]
yp[i] += yi[i] + yj[i]
zp[i] += zi[i] + zj[i]
listLen = len(cutofflist[i])
if listLen > 0: # List length must be greater than zero or division by zero error.
ljsum[i] = ljsum[i] / len(cutofflist[i])
return ljt/points, ljminimum, ljminimum_i, ljminimum_j, ljmaximum, ljmaximum_i, atmax


def render():
# Precalculate all the rotated X,Y screen locations and their colors.
for i in range(points):
az, el, r = cart2sphere(xp[i], yp[i], zp[i]) # Cartesian to Spherical (az, el, r)
xr, yr, zr = sphere2cart(az + rotation, el, r) # Spherical to Cartesian.
xyzi[i] = [xr, yr, zr, i]
unsortedxyzi[i] = xyzi[i]
# To properly draw on top of each other sort the list according to what is furthest away, [1] = by Y axis.
xyzi.sort(key=lambda x: x[1])
if drawcolors != 0:
for i in range(points): # Draw all points to screen.
drawvesicle(xc1, 1, i) # Draw left view i point to screen in proper color.
drawvesicle(xc2, -1, points-1-i) # Draw right view i point to screen in proper color.
if drawwaves is True:
for i in range(points):
drawsignallines(xc1, 1, i) # Draw left view signal lines from given i point.
drawsignallines(xc2, -1, points-1-i) # Draw right view signal lines from given i point.
if axeslines: # Draw x,y,z axis lines to screen.
ra = vesiclestartR * 1.1 # Length of axis line to draw.
drawaxisline(0, 0, ra, Yellow)
if rotation > threeQtrRadian: # Green x behind blue y?
drawaxisline(ra, 0, 0, Green)
drawaxisline(0, ra, 0, Blue)
if rotation < threeQtrRadian: # Green x on top of blue y?
drawaxisline(ra, 0, 0, Green)
if pops > 0:
drawpop(xc1, 1, White, xpop, ypop, zpop)
drawpop(xc2, -1, White, xpop, ypop, zpop)
drawpop(xc1, 1, Cyan, xnew, ynew, znew)
drawpop(xc2, -1, Cyan, xnew, ynew, znew)
pygame.display.update() # Display everything drawn into screen buffer.


# Draw one i point filled circle in sorted xyz list. For next loop in render draws all vesicle points.
def drawpop(xc, direction, c, x, y, z):
az, el, r = cart2sphere(x, y, z) # Cartesian to Spherical (az, el, r)
xr, yr, zr = sphere2cart(az + rotation, el, vesicleR) # Spherical to Cartesian.
if yr*direction < 0:
return
xs = int(xr * direction + xc) # Screen pixel x location to draw at.
ys = int(zr + yc) # Screen pixel y location to draw at.
cr = int(spacing/2) + 12 # Circle radius after distance away is factored in.
pygame.draw.circle(screen, c, [xs, ys], cr, 5) # Draw circle to screen.


# Draw one i point filled circle in sorted xyz list. For next loop in render draws all vesicle points.
def drawvesicle(xc, direction, i):
x, y, z, ui = xyzi[i] # Get x,y,z and unsorted 'i' from sorted list.
if y*direction < -spacing:
return
c = Black # Default fill color is black, when drawcolors=0
if drawcolors == 1: # Keyboard 'c' changes, 1=LennardJones colors.
c = ljcolor(ui) # Use Lennard-Jones amount, warm to hot colors.
if drawcolors == 2: # Keyboard 'c' changes, 2=Connections based colors.
c = connectcolor(ui) # Electronic color code, number of connections each.
xs = int((x*direction) + xc) # Screen pixel x location to draw at.
ys = int(z + yc) # Screen pixel y location to draw at.
d = (((y*direction)+vesiclestartR)/spherewidth*.4) + .6 # Fraction according to how far away, depth.
if d > 1: # If fractional depth is greater than 1 then
d = 1 # reduce to 1, close as possible to observer.
cr = int(d * spacing * .5) + 1 # Circle radius after distance away is factored in.
pygame.draw.circle(screen, c, [xs, ys], cr, 0) # Draw filled circle to screen.
if outlines:
pygame.draw.circle(screen, White, [xs, ys], cr, 1) # Draw outline circle around it.
if drawnumbers:
fnt = pygame.font.SysFont('arial', int(halfspacing))
txtsurf = fnt.render(str(ui), False, White)
textwidth = fnt.size(str(ui))[0]
screen.blit(txtsurf, (xs - textwidth/2, ys - quarterspacing))


# Draw signal lines outward from one given i in sorted xyz list. For next loop in render draws all vesicle points.
def drawsignallines(xc, direction, i):
global connectlist
global outs
x, y, z, ui = xyzi[i] # Get variables from sorted list of x,y,z points.
if y*direction < -spacing:
return
x1 = int((x*direction) + xc) # Screen pixel x location to draw at.
y1 = int(z + yc) # Screen pixel y location to draw at.
d = (((y*direction)+vesicleR) / spherewidth) # Fraction according to how far away, depth.
if d > 1: # If fractional depth is greater than 1 then
d = 1 # reduce to 1, close as possible to observer.
outbits = outs[ui]
for j in range(len(connectlist[ui])):
x2, y2, z2, ui2 = unsortedxyzi[connectlist[ui][j]]
b = outbits[j]
if b == 1:
thick = int((b * d * spacing / 3) + 1)
c = OutputSignalColor[b]
if c != Black:
pygame.draw.line(screen, c, [x1, y1], [(x2*direction)+xc, z2+yc], thick)
inbits = ins[ui]
for j in range(len(connectlist[ui])):
x2, y2, z2, ui2 = unsortedxyzi[connectlist[ui][j]]
b = inbits[j]
thick = int((b * d * spacing / 9) + 1)
c = InputSignalColor[b]
if c != Black:
pygame.draw.line(screen, c, [x1, y1], [(x2*direction)+xc, z2+yc], thick)


def save():
os.system("ffmpeg -r 1 -i img%01d.png -vcodec mpeg4 -y movie.mp4")


def drawaxisline(x, y, z, c):
az, el, ra = cart2sphere(x, y, z) # Cartesian to Spherical (az, el, r)
x1, y1, z1 = sphere2cart(az + rotation, el, ra) # Convert from Spherical to Cartesian.
pygame.draw.line(screen, c, (xc1, yc), (x1+xc1, z1+yc), 1)


# Vary brightness according to amount of force. Positive is orange to white, negative (attraction) blue.
def ljcolor(t):
if ljsum[t] >= 0:
poscolor = int(ljsum[t] * (256/maxf) * 10)
if poscolor > 255:
poscolor = int(poscolor/200)
if poscolor > 255:
poscolor = 255
return 255, 240, poscolor
else:
return poscolor, poscolor * .75, 0
if ljsum[t] < 0:
negcolor = int(-ljsum[t] * 2000)
if negcolor > 255:
negcolor = 255
return 0, negcolor*.6, negcolor

def connectcolor(n):
cns = len(connectlist[n])
if cns > 9:
cns = 9
return ElectronicColorCode[cns]


def timestepsignals():
global ins
global outs
global signalactivity
global signalcollapsepoint
signalactivity = 0
ins = [[] for _ in range(points)] # Clear Input list output bits from neighbors are gathered into.

# For all points in array, gather bits from neighbors to a list for each. Example: ins[i,[1,0,0,0,1,0]]
for i in range(points): # For each point in vesicle, each with connections to 'j's
for j in range(len(connectlist[i])): # For each of the (usually 6) 'j' connections to other cells,
if j < len(outs[i]): # If output connection still exists in connection list then
b = outs[i][j] # bit to save is from output list.
else: # Else
b = 0 # to prevent program error append with 0 to fill void later.
ins[connectlist[i][j]].append(b) # Append bit from 'j' output to corresponding input line at 'i'

# If any bits are 1 then propagate wave using NOT to output inverted/opposite bit pattern. 1=0 and 0=1
for i in range(points): # For each 'i' point in vesicle, each with connections to 'j's
outs[i] = []
setbits = sum(ins[i])
signalactivity += setbits # Signal activity equals total number of bits set, entire surface.
if setbits > 0: # If there are any bits set then
outs[i] = [int(not b) for b in ins[i]] # all 1's in output list become 0, and all 0 becomes 1, negates.
else: # Else
outs[i] = ins[i] # Output bits remain all 0.

# If no signal activity then start wave by setting all output connections of a randomly selected place to 1.
if signalactivity == 0: # If no signal activity at all then
# startwave(random.randint(0, points-1)) # Randomly choose the point to start wave from.
startwave(4) # Start wave at same given point each time.

# If uncontrolled signal activity then clear all outputs.
if signalactivity > points * 2: # If more than average two connections per cell then
clearwave()


# Start wave from cell at point number given by value of i
def startwave(i):
global signalactivity
outs[i] = [1] * len(connectlist[i]) # Output 1, action potential, to all neighbors.
signalactivity += 6
print('~~~~~~~~~~ Wave started from point number', i, '~~~~~~~~~~')


# For when extreme hexagonal misalignment causes uncontrollable seixure-like activity.
def clearwave():
for i in range(points): # for all points in network.
outs[i] = [0] * len(outs[i]) # fill list with 0's, same length as there are out connections.
screen.fill(White)
myfont = pygame.font.SysFont('arial', 68)
textsurface = myfont.render("STOP", False, Black)
screen.blit(textsurface, (height * .82, height * .8))
print('~~~~~~~~~~ All connection outputs were set to 0 ~~~~~~~~~~')


# To Squeeze and Release pressure from all around.
def pressureRamp(press, direction):
newpressure = press + (pressureRampStep * direction)
if newpressure > pressureMax:
newpressure = pressureMax
direction = -1
if newpressure < pressureMin:
newpressure = pressureMin
direction = 1
return newpressure, direction


def listlennardjones():
lj = 0
for i in range(32):
d = i / 16
# ljr = d / sp # Lennard-Jones radius. At ljr=1 force is 0.
if d > 0: # Prevent divide by zero error when LJ radius is zero.
lj = 4 * epsilon * ((sigma / d) ** 12 - (sigma / d) ** 6)
print(lj, d)


#######################################################################################################################
xp, yp, zp = xyzread("StartupXYZ.txt") # Load either startup file, or random positions for new vesicle.
points = len(xp) # Number of points total in vesicle.
ljavr = 0.0
ljmax = 0.0
ljsum = [0 for _ in range(points)] # Initialize Lennard-Jones sum array.
xyzi = [[] for _ in range(points)] # List to sort, to draw according to their distance from observer.
unsortedxyzi = [[] for _ in range(points)] # Above list before sorting.
grouplist = [[] for _ in range(points)] # Nearby points in larger group, stay with for very long time.
cutofflist = [[] for _ in range(points)] # Points within cutoff range.
connectlist = [[] for _ in range(points)] # Connections from each i point, to neighbors j.
connections = [[] for _ in range(points)] # Connections from each i point, to neighbors j.
latlong = [[] for _ in range(points)] # Spherical coordinates for each point.
cutoffconnects = [[] for _ in range(points)] # Number of connections, usually 6.
ins = [[] for _ in range(points)] # Input (action potential) bits, from each 1 bit connection.
outs = [[] for _ in range(points)] # Output (action potentials) bits, to each 1 bit connection.

# PyGame Initialization
pygame.init()
screen = pygame.display.set_mode(screensize)
pygame.display.set_caption("Lennard Jones Potential Vescicle. A=Axis, C=ColorCode, F=Forces, H=HoldToSphere, N=indexNumber, P=Pressure, W=Waves, Spacebar=Pause, S=Save-XYZ.txt")
fps = pygame.time.Clock()
pygame.font.init()
# Display 'Initializing nnnnn points' message.
myfont = pygame.font.SysFont('arial', 50)
s = 'Initializing ' + str(points) + ' points'
textsurface = myfont.render(s, False, White)
screen.blit(textsurface, (10, 10))
pygame.display.update() # Display everything drawn into screen buffer.
fps.tick(2)

# Initialize vesicle size parameters.
sas = 4*pi*(screenradius*screenradius) # Area of sphere at this given screen radius.
sar = sas/points # Divide for area of one Lennard-Jones radius, between points.
rbp = np.sqrt(sar/pi) * 2 # Radius between points, given area.
spacing = rbp-((1/points)*rbp) # Radius between points, ideal spacing.
halfspacing = spacing / 2
quarterspacing = spacing / 4

listlennardjones()

avcutoffs = 0
cutoffcalctime = 10 # Time steps between updates of i,j points within cut-off range.
cutoffupdate = cutoffcalctime # Time cycles per update of cutoff list. +1 to upda
groupcalctime = 10 # Time steps between updates of i,j points within larger group.
groupupdate = groupcalctime # Time cycles per update of group list. +1 to update on startup.

xpop = 0
ypop = 0
zpop = 0
xnew = 0
ynew = 0
znew = 0

# Initialize keyboard entered toggle controls.
paused = False
axeslines = False
pointschange = True
forces = True # Must start True or list of points within cutoff range is not generated.
drawcolors = 2 # 0=none, 1=LennardJones, 2=Connections
pressureChanges = False
drawwaves = True # To enable or disable drawing of signal propagation lines.
drawnumbers = False
holdradius = True
if points < 2000:
outlines = True
else:
outlines = False
screen.fill(Black) # Clear screen of "Intitializing nnnn Points" message.

while True:
for event in pygame.event.get():
if event.type == pygame.QUIT:
pygame.quit()
sys.exit()
# Space=pause A = Axes,x=green,y=bluez=yellow. F=force. H=hold Radius. O=outlines.
# P=Pressure. S=save. W=waves. C= color, LJ, # connections. none waves only.
if event.type == pygame.KEYUP:
if event.key == pygame.K_SPACE: # Space bar = pause.
paused = not paused
if event.key == pygame.K_a: # a = Axis lines in left draw: x=green, y=blue, z=yellow.
axeslines = not axeslines
if event.key == pygame.K_f: # f = Lennard-Jones Force on/off. None = fast stable waves.
forces = not forces
if event.key == pygame.K_h: # h = Hold to (an inner volume) sphere. None = particles drift
holdradius = not holdradius
if event.key == pygame.K_n: # n = Number the points drawn to screen with their index number.
if points < 5000:
drawnumbers = not drawnumbers
if event.key == pygame.K_o: # o = Outline around circles.
outlines = not outlines
if event.key == pygame.K_p: # p = Pressure changes like pounding surf waves, shaken in jar.
pressurechanges = not pressurechanges
if event.key == pygame.K_s: # s = Save new StartupXYZ.txt file, for later automatic restart.
xyzwrite('StartupXYZ.txt')
if event.key == pygame.K_w: # w = Wave signal draw to screen. Inputs turn violet, Outputs white.
drawwaves = not drawwaves
if not drawwaves and drawcolors == 0:
drawcolors = 1 # No blank screen.
if event.key == pygame.K_c: # c = change colors: None, Lennard-Jones, Number of connections.
drawcolors += 1
if drawcolors > 2:
if drawwaves: # Draw nothing while showing waves.
drawcolors = 0
else: # Else toggle to Lennard Jones so that a blank screen is not shown.
drawcolors = 1
if not paused:
myfont = pygame.font.SysFont('arial', 25)
# Pressure changes as in shaking in a jar of water, or caught inside pounding shoreline waves.
vesiclestartR = int(screenradius) # Start size of vesicle sphere, from calculated full size.
vesicleR = vesiclestartR * adjustR
cutoffR = forceCutOff * spacing * adjustR # Distance two points are too far away to exert force.
groupR = groupCutOff * spacing * adjustR
connectR = spacing * connectCutOff * adjustR # Distance two points are touching, bumping each other.
spherewidth = vesiclestartR * 2
ljdiff = 0

# Update positions of points, and rotation to drawn on screen.
if forces:
if pressureChanges:
adjustR, pressureRampDirection = pressureRamp(adjustR, pressureRampDirection)
# adjustR = 1 + random.uniform(-pressure, 0) # Pressure changes of churning water, shaking.
else:
adjustR = 1 # Stay full radius.

cutoffupdate += 1
if cutoffupdate > cutoffcalctime:
cutoffupdate = 0
groupupdate += 1
if groupupdate > groupcalctime:
groupupdate = 0
update_grouplist()
avcutoffs = update_cutofflist()
ljmaxwas = ljmax
ljavr, ljmin, ljmin_i, ljmin_j, ljmax, ljmax_i, atmaxforce = update_positions()
if atmaxforce != 0:
pressureRampDirection = 1
else:
pressureRampDirection = -1

if ljmax < 10: # If activity is low then increase Lennard Jones force.
dt *= 1.1
if dt >= .5:
dt = .5
else: # else decrease.
dt *= .95
if dt <= .001:
dt = .001
# if ljavr < .5: # If activity is low then increase Lennard Jones force.
# dt *= 1.1
# if dt >= .5:
# dt = .5
# else: # else decrease.
# dt *= .95
# if dt <= .001:
# dt = .001
update_connections()
update_angles()
cutoffcalctime = int(1000 / ljmax) * 5 # Timesteps before recalculating distances.
# cutoffcalctime = int(1 / (np.sqrt(abs(ljavr)) / 5) + 1) # Timesteps before recalculating distances.
# cutoffcalctime = int(2 / ljavr) * 5 # Timesteps before recalculating distances.
if cutoffcalctime > 200:
cutoffcalctime = 200
print(timestepnumber, " LJavr=", ljavr, " LJmin=", ljmin, " LJmax=", ljmax,
" dt=", dt, " r=", vesicleR, "of", screenradius)

# Popping out under pressure.
indexmin = ljsum.index(min(ljsum))
ljdiffwas = ljdiff
ljdiff = abs(ljmax - ljmaxwas)
if ljdiff < .00002:
print(ljdiff)
popwait += 1
if ljmax > 10: # If there is enough force on any one point.
pressureRampDirection = -1
if popwait > popwaitcount: # If there has been enough time since last swap
# Pop out the particle with most force around it
popwait = -1 # start new count for waiting again before swapping again
indexmax = ljsum.index(max(ljsum))
xpop = xp[indexmax]
ypop = yp[indexmax]
zpop = zp[indexmax]
azm, ele, rad = vectorsum(indexmin, cutofflist[indexmin])
xo, yo, zo = sphere2cart(azm, ele, rad)
xnew = xp[indexmin] + xo
ynew = yp[indexmin] + yo
znew = zp[indexmin] + zo
xp[indexmax] = xnew
yp[indexmax] = ynew
zp[indexmax] = znew
print(timestepnumber, '>>>>>> One popped out and was sent to one can pop in. <<<<<<')
myfont = pygame.font.SysFont('arial', 26)
textsurface = myfont.render('Pop-Out, 1 Moved', False, White)
screen.blit(textsurface, (height * .8, height * .84)) # height is half the width of screen
myfont = pygame.font.SysFont('arial', 22)
textsurface = myfont.render('From Blue to White Circle', False, White)
screen.blit(textsurface, (height * .76, height * .90)) # height is half the width of screen
pops =+ pops
cutoffupdate = cutoffcalctime # set to next update list of points in cutoff range.
groupupdate = groupcalctime # and set to next update list of points in larger groups.
# Propagate signals.
timestepsignals()

# Label screen:
myfont = pygame.font.SysFont('arial', 27)
textsurface = myfont.render('Front', False, White)
screen.blit(textsurface, (height * .83, height/10)) # height is half the width of screen
textsurface = myfont.render('Back', False, White)
screen.blit(textsurface, (height * 1.06, height/10))
myfont = pygame.font.SysFont('arial', 18)
textsurface = myfont.render(str(points) + ' mol/cell/column', False, White)
screen.blit(textsurface, (5, 0))
textsurface = myfont.render('LJmax: ' + str(round(ljmax, 4)), False, White)
screen.blit(textsurface, (5, height * .04))
textsurface = myfont.render('On: ' + str(signalactivity), False, White)
screen.blit(textsurface, (5, height * .08))
myfont = pygame.font.SysFont('arial', 20)
textsurface = myfont.render('TimeStep ' + str(timestepnumber), False, White)
screen.blit(textsurface, (height * .75, 0))
myfont = pygame.font.SysFont('arial', 18)
textsurface = myfont.render("LJ-Forces=" + str(forces) + " HoldRadius=" + str(holdradius) + " Pressure=" + str(pressureChanges), False, White)
screen.blit(textsurface, (5, height - 24))
if drawcolors == 2:
textsurface = myfont.render('Connections: Red=2 Orng=3 Yell=4 Grn=5 Blu=6 Viol=7 Gray=8', False, White)
if drawcolors == 1:
textsurface = myfont.render('Lennard Jones Potential, White = Extreme Pressure', False, White)
if drawcolors == 0:
textsurface = myfont.render('Shown as Points', False, White)
screen.blit(textsurface, (height, height - 24))
myfont = pygame.font.SysFont('arial', 18)
textsurface = myfont.render('Gary Gaulin, 2023', False, White)
screen.blit(textsurface, (height * 1.7, 4)) # Height is half the width of screen
render() # Draw to screen.

pygame.display.update # Display screen that was being drawn, to monitor display.
# pygame.image.save(screen, 'Images/' + str(timestepnumber) + '.png') # Save image file for making gif video.
fps.tick(50)
screen.fill(Black) # Screen can now be cleared for next draw.
# Rotate sphere a small amount for next draw.
rotation += rotationstep # Add screen rotation amount, for next rendering.
if rotation > radian: # Keep rotation amount in range from 0 to 1 Radian.
rotation = 0
timestepnumber += 1

No comments: