Oh yeah, someone finally asked something about my field.

So, maybe COMSOL is what you’re looking for, but if you want the cool, smooth and amazing stuff of the pros:

Time to learn numerical calculus, like verlet integration, maybe some calculus for turning potential fields into force (maybe some forward differentiation with dual numbers, just like sympy does hahaha, you can use that too). If it’s a hairy equation, that takes some hypergeometric, bessel, etc, scipy has all the special equations.

I think verlet + sympy for turning the potential fields into the force fields at the particle spots will do amazing, if you don’t want paper ready simulations. Like, pos[n, 3], v[n, 3], a[n, 3] (maybe some torque if you want) would be almost all you need.

For visualization, matplotlib has a lot of stuff, you can make animations and videos.

If you want to simulate the fields itself, you can use finite differences or finite element methods.

If you need a lot of performance, I recommend learning Julia or the GODLY Fortran 2008 (fprettify, fortls, gfortran will help to build a modern Fortran program that you can call from python using cdll and ctypes with the ieee c libs for fortran).

For advanced stuff, you have Geant4. Fully open source software developed by the CERN and collaborators. It has a steep learning curve, but the feeling of simulating advanced physics (and somewhat underestading them) is amazing.

At a past job a physicist on my team used a free code that was developed at CERN for radiation simulations. Can’t recall what it was called but look into open stuff they’ve developed

Could you elaborate on what you want to make? If it’s lots of particles (thousands) then look into DEM software like ansys rocky or LIGGGHTS

A lot of people already mentioned it but geant4 is the gold standard for that kind of simulation. It could also be worth looking into pythia if you’re interested in simulating high energy events with geant4

Edit: Corrected the integration

This is a simulation I made for a ball falling (the principles can be used everywhere):

import numpy as np
import matplotlib.pyplot as plt


def calc_acceleration(pos, v):
    return np.array([0, 0, -9.8])


def integrate(pos, v, a_old, delta_t):
    pos_new = pos + v * delta_t + 0.5 * a_old * delta_t * delta_t

    a_new = calc_acceleration(pos_new, v)

    v_new = v + 0.5 * (a_old + a_new) * delta_t

    if pos_new[2] < 0.0:
        v_new[2] = -v_new[2]

    return pos_new, v_new, a_new


# Constants
delta_t = 0.01
t_max = 10

# Time evolution
times: np.ndarray = np.arange(0, t_max, delta_t)

# Initial state
v = np.array([15.0, -3.0, 15.0])
pos = np.array([0.0, 0.0, 7.0])
positions = np.zeros([times.size, 3])
a_old = calc_acceleration(pos, v)

# Simulating
for i in range(times.size):
    pos, v, a_old = integrate(pos, v, a_old, delta_t)
    positions[i, :] = pos

# Plotting
fig = plt.figure()
ax = fig.add_subplot(projection="3d")
ax.scatter(positions[:, 0], positions[:, 1], positions[:, 2])
plt.show()

Very simple, you can see a ball falling and getting up again from that starting point, so you just see at least some starting point on how to do this stuff, the basics are not that difficult. But every simulation is like this, just more complicated hahaha.