I wrote a Python code to do MD for a Lennard-Jones fluid but the VACF is wrong: What might be the problem?

I am trying to write a molecular dynamics simulation for a Lennard-Jones fluid in a box with periodic boundary conditions. The box has no net momentum. I am writing this in python.

I have written a library of functions to set up my box of particles. Then I am implementing it in a separate script.
Here is dynamics.py:

import itertools
import numpy as np
import random
import time
import math

random.seed(time.time())

#create a box of particles
#make object Box which will hold all the particles
class Box:
    def __init__(self, numberOfParticles, boxLength, dimension, sigma, epsilon, temperature, dt):
        self.numberOfParticles = numberOfParticles
        self.boxLength = boxLength
        self.dimension = dimension
        self.sigma = sigma
        self.epsilon = epsilon
        self.temperature = temperature
        self.dt = dt #time step
        ##### non given quantities
        self.nrho = numberOfParticles/(boxLength**(dimension)) #number density
        self.particlePositions = np.zeros((numberOfParticles, dimension)) #do a cubic lattice
        self.particleVelocities = self.boxLength*(np.random.rand(numberOfParticles, dimension)-0.5) #assign randomly
        self.particleForces = np.zeros((numberOfParticles, dimension))
    #
    #
    #now to evaluate energy of configuration
    #evaluating kinetic energy of the system

    def latticePositions(self):

        pointsInLattice = math.ceil(self.numberOfParticles**(1/3))

        spots = np.linspace(0, self.boxLength, num=pointsInLattice, endpoint=False)
        count = 0

        for p in itertools.product(spots, repeat=3):
            p = np.asarray(list(p))
            self.particlePositions[count, :] = p
            count += 1
            if count>self.numberOfParticles-1:
                break
        #
        return self
#

    def evaluateKineticEnergy(self):
        #square every element, add up the elements of each row
        kineticEnergy = 0.5*np.sum(np.square(self.particleVelocities))
        return kineticEnergy
#
    #I will be selecting a particle, and summing up all the potential energy arising
    #due to interactions with every other particle

    def evaluatePotentialEnergy(self):
        energy = 0
        for i in range(self.numberOfParticles):
            for j in range(i+1, self.numberOfParticles):
                displacement = self.particlePositions[i,:]-self.particlePositions[j,:]
                for k in range(self.dimension):
                    if abs(displacement[k])>self.boxLength/2:
                        displacement[k] -= self.boxLength*np.sign(displacement[k])
                r = np.linalg.norm(displacement,2) #finding euclidean distance between two particles
                energy += (4*self.epsilon*((self.sigma/r)**12-(self.sigma/r)**6)) #evaluating potential energy, multiply by 2?
        return energy
    #sum of potential and kinetic energy is equal to the total energy

    def evaluateTotalEnergy(self):
        totalEnergy = self.evaluatePotentialEnergy()+self.evaluateKineticEnergy()
        return totalEnergy
    #end of energy calculations
    #
    #
    #find the force each particle is experiencing due to the other particles
    #force = - gradient of potential
    def evaluateForce(self):
        self.particleForces = np.zeros((self.numberOfParticles, self.dimension))
        def LJForce(displacement):
            r = np.linalg.norm(displacement, 2)
            force = 48/(r**2)*(1/(r**12)-0.5*1/(r**6))*displacement
            return force
        for i in range(self.numberOfParticles):
            for j in range(i+1, self.numberOfParticles):
                rij = self.particlePositions[i,:]-self.particlePositions[j,:]
                for k in range(self.dimension):
                    if abs(rij[k])>self.boxLength/2:
                        rij[k] -= self.boxLength*np.sign(rij[k])
                rji = -rij
                self.particleForces[i,:] += LJForce(rij)
                self.particleForces[j,:] += -self.particleForces[i,:]
        return self
    #end of force evaluations
    #make sure total momentum of box is zero
    def stationaryCenterOfMass(self):
        #ensure center of mass is stationary
        v_cm = np.mean(self.particleVelocities, axis=0)
        self.particleVelocities = self.particleVelocities - v_cm
        return self
#
#
#
def VelocityVerletTimeStepping(currentBox):
    previousParticleForces = currentBox.particleForces
    currentBox.particlePositions = (currentBox.particlePositions + currentBox.particleVelocities*currentBox.dt + 0.5*currentBox.particleForces*(currentBox.dt)**2)%(currentBox.boxLength)
    currentBox = currentBox.evaluateForce()
    currentBox.particleVelocities = currentBox.particleVelocities + 0.5*(previousParticleForces + currentBox.particleForces)*currentBox.dt
    return currentBox
#

Now I am calling these functions to perform time-evolution. I am evaluating energy and ACF and plotting them against time to see if I have done this right.

import dynamics
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
import math

#parameters of simulation
numberOfParticles = 100
dimension = 3
sigma = 1
epsilon = 1
temperature = 1
#max_iterations = 100
boxLength = 10*sigma
dt = 0.001 #time step
kb = 1 #boltzmann
##
#
#

nmax = 100 #number of time steps to take

#set up the box

currentBox = dynamics.Box(numberOfParticles, boxLength, dimension, sigma, epsilon, temperature, dt)

#placing particles in a lattice
currentBox = currentBox.latticePositions()

#ensuring box has no net momentum
currentBox = currentBox.stationaryCenterOfMass()

#calculating forces on particles in the box
currentBox = currentBox.evaluateForce()


#making a list of particle positions and velocities for acf and what not
particlePositionsList = [currentBox.particlePositions]
particleVelocityList = [currentBox.particleVelocities]

#making a list of energies at various time steps to plot later
energy = np.zeros(nmax+1,)
energy[0] = currentBox.evaluateTotalEnergy()

timepoints = np.arange(nmax+1)*currentBox.dt

#start time stepping routine
#time points = 0, 1, 2, ..., nmax
for i in range(nmax):

    #do the time step, knowing that currentBox already knows the particle forces at the moment
    currentBox = dynamics.VelocityVerletTimeStepping(currentBox) #evaluates forces on particles, updates particle positions and velocities
    energy[i+1] = currentBox.evaluateTotalEnergy()

    particlePositionsList.append(currentBox.particlePositions)
    particleVelocityList.append(currentBox.particleVelocities)

#
#print(energy)

ACF = np.zeros(nmax+1,)

for i in range(nmax+1):
    for j in range(nmax+1-i):
        ACF[i] = ACF[i] + np.sum(particleVelocityList[j]*particleVelocityList[j+i])
        #ACF[j] = ACF[j] + np.sum(particleVelocityList[i]*particleVelocityList[j+i]) #this one works
    #
    ACF[i] = ACF[i]/(nmax+1-i)
#


plt.plot(timepoints, energy)
plt.title("Energy over time")
plt.xlabel("time")
plt.ylabel("energy")
plt.ylim(np.amin(energy)*0.999, np.amax(energy)*1.001)
plt.show()

plt.plot(timepoints, ACF)
plt.title("Normalized VACF plot")
plt.xlabel("time")
plt.ylabel("VACF")
plt.show()

Here are the results:

Energy looks good, but ACF looks really wrong. This is how it is supposed to look like:

I am unsure where I am going wrong here. I am a lone software engineer thrust into the world of physics and molecular modelling, so any advice you have would be appreciated!!

Edit 1: After initializing velocity according to a Gaussian distribution,

self.particleVelocities = self.boxLength*(np.random.normal(0, 1, (numberOfParticles, dimension))) #assign velocity as per normal distr

These are the figures I get: