I have written the program attached bellow that is supposed to calculate and display the electrostatic potential of Silicon atoms arranged in a CFC unit cell.

The program is arranged in a way that define the electrostatic potential function first, followed by the atom position array, that is followed by the mesh definition and electrostatic potential initialization, finally the calculation for the potential is performed defining the potential inside the atom radius as positive infinity, for the graph plot as a volume type in plotly.

The problem is that the code won’t compile nor give an error that I can detect and potentially fix.

# -*- coding: utf-8 -*-
"""
Created on Thu Apr 25 11:02:14 2024

@author: Luiz Felipe Pompermaier

Silicon Lattice with periodic eletrostatic potencial.
"""

import plotly.graph_objects as go
import numpy as np

def V(q,r):
    '''
    Parameters
    ----------
    q : Integer
        A integer multiple of the fundamental charge e.
    r : float
        Radial distance from the point charge.

    Returns
    -------
    V : Float
        Value for the electrostatic potential at a given point.
    '''
    e0 = 8.8541878128*10**(-12) # Vacuum Permittivity (C / V m)
    e = 1.602176634*10**(-19) # Fundamental charge (C)
    
    V = (1/(4*np.pi*e0))*((q*e)/r)
    
    return V

print("Runing ...")

# List of attoms in the first unit cell in fractions of the unit cell vectors.
Atom = np.array(
    [[0,0,0],[1,0,0],[0,1,0],[0,0,1],[1,1,0],[1,0,1],[0,1,1],[1,1,1], # 8 corners
    [.5,.5,0],[.5,.5,1],[.5,0,.5],[.5,1,.5],[0,.5,.5],[1,.5,.5], # 6 center faces
    [.25,.25,.25], [.75,.75,.25],[.75,.25,.75],[.25,.75,.75]] # 4 central atoms
    )
print(Atom)

# Value for the lattice constant of a cubic silicon cell
a = 5.430986 # in angstrom

# Radius of a silicon atom in the lattice
r = 1.11 # in angstrom

# Deffine the mesh of the axis (1000 points/axis)
mesh = 1000
x,y,z = np.linspace(0, (a*10**10),mesh,True)
print(x)

# Inicialize the eletrostatic potential
value = np.empty((mesh+1,mesh+1,mesh+1), float())
print(value)

# For every atom in Atom calculate the eletrostatic potential generated
# by that atom at a given position in the lattice
# If inside a certain atom radius in the lattice then define V = infinity
for ai in range(len(Atom)):
    for xi in range(len(x)):
        for yi in range(len(y)):
            for zi in range(len(z)):
                radius = np.sqrt((Atom[ai][0]*a-xi)^2+(Atom[ai][1]*a-yi)**2+(Atom[ai][2]*a-zi)**2)
                if radius <= r:
                    value[xi][yi][zi] = float('inf')
                else:
                    value[xi][yi][zi] += V(4,radius*10**(-10))
print(value)
                    
# Get the minimum and maximum values of potential
Vmin = value.min()
Vmax = value.max()
print(Vmin,Vmax)
# Plot the graph as a volume plot
fig = go.Figure(data=go.Volume(
    x=x.flatten(),
    y=y.flatten(),
    z=z.flatten(),
    valeus= value.flatten(),
    isomin = Vmin,
    isomax = Vmax,
    opacity=0.1, # needs to be small to see through all surfaces
    surface_count=25, # needs to be a large number for good volume rendering
    ))

fig.show()

I have already read the code in detail to search for possible mistakes and improvements, although I have found and made some improvements from the initial code, I haven’t found any mistakes.