How can I minimize a potential energy function in 2D with periodic boundary conditions using python?

Right now, I am trying to minimize an electrostatic potential energy function in 2D with periodic boundary conditions and two charges q1 and q2. Because of the conditions, the two charges should end up on the diagonal at some point that is L/sqrt(2) away from each other where L is the length of the box. As of now, I have the following code:

import matplotlib.pyplot as plt
from scipy.optimize import minimize

def PE_func(x):
    L,q1,q2,x1,y1=13,1,1,0,0
    #this starts the periodic boundary conditions
    if x[0]>L:
        x[0]=x[0]-L
    if x[0]<0:
        x[0]=x[0]+L
    if x[1]>L:
        x[1]=x[1]-L
    if x[1]<0:
        x[1]=x[1]+L
    delta_X=np.abs(x1-x[0])
    delta_Y=np.abs(y1-x[1])
    if (delta_X > L/np.sqrt(2)):
        delta_X=L-delta_X
    if (delta_Y > L/np.sqrt(2)):
        delta_Y=L-delta_Y
    #this ends the periodic boundary conditions
    r=np.sqrt((delta_X**2)+(delta_Y**2))
    PE=(9e9*q1*q2)/r
    print(f'The distance between the particles is {r} m.')
    return PE

#coords=[X1,X2]
coords=[np.random.randint(0, 13),0]
print(f'Initial coordinates are {coords}.')

result=minimize(PE_func,coords,tol=1e-6)
print(f'Final coordinates are {result.x}.')

This code minimizes the PE function in 1D where one charge is kept still at the origin and the other charge is allowed to move along a "wire" of length L. It still has periodic boundary conditions. It works 90% of the time and gives an end result of a separated distance roughly equal to L/sqrt(2). However, it still doesn't give me the coordinates that correspond to that distance. Also, if I try to turn it into 2D by using this code:

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize

def PE_func(x):
    L,q1,q2=13,1,1
    #this starts the periodic boundary conditions
    if x[0]>L:
        x[0]=x[0]-L
    if x[0]<0:
        x[0]=x[0]+L
    if x[1]>L:
        x[1]=x[1]-L
    if x[1]<0:
        x[1]=x[1]+L
    delta_X=np.abs(x[1]-x[0])
    delta_Y=np.abs(x[3]-x[2])
    if (delta_X > L/np.sqrt(2)):
        delta_X=L-delta_X
    if (delta_Y > L/np.sqrt(2)):
        delta_Y=L-delta_Y
    #this ends the periodic boundary conditions
    r=np.sqrt((delta_X**2)+(delta_Y**2))
    PE=(9e9*q1*q2)/r
    print(f'The distance between the particles is {r} m.')
    return PE

#coords=[X1,X2,Y1,Y2]
coords=[np.random.randint(0, 13,4)]
print(f'Initial coordinates are {coords}.')

result=minimize(PE_func,coords,tol=1e-6)
print(f'Final coordinates are {result.x}.')

It doesn't work and gives different results everytime. Any help would be appreciated