Generally speaking, most work in molecular dynamics tries to simulate how actual molecules behave (i.e. quantum mechanics) and it doesn't sound like you want to go down that rabbit hole. I completely agree, but I'll begin with a disclaimer that looking up "molecular dynamics" probably won't turn up the kind of results you want.
Since your comments describe wanting a simple simulation based on kinetic energy, forces and collisions (i.e. Newtonian mechanics), it will probably be more helpful to look into more general numerical simulations of Newtonian mechanics.
The type of model you are thinking of is very similar to how Newtonian gravity is usually modeled. They are basically identical from a programming perspective, the only difference is that charges are signed (so you can have repulsion as well as attraction).
The simplest way to approach this kind of problem is with numerical integration - set up your initial conditions $x(0)$ and $v(0)$, then update both little bit in each time step $\Delta t$. The most basic form is:
$$ v_x(t + \Delta t) = v_x(t) + a_x(t) \Delta t $$
$$ x(t + \Delta t) = x(t) + v_x(t) \Delta t $$
I used the x-axis here, but you can update the y and z components the same way. This is a really crude numerical integrator (Euler's method) but it's a good start to understand the concept. Once you feel comfortable with it, you can use more advanced, accurate numerical integration schemes.
This acceleration function is really general and you can include as many relevant terms as you see fit. To capture the force of the electric field, you can calculate the acceleration using Coulomb's law, but you can expand it with additional components (e.g. the force generated by the Lennard-Jones potential).
The process to calculate acceleration is typically to calculate the net force on a particle, then divide by its mass to get its acceleration. Again, this is essentially how gravitational acceleration is calculated (so feel free to refer to that material for guidance).
If you want a (really, really) simple model of molecules, you can treat them using rigid body dynamics. You start off by calculating the force on each atom in exactly the same way, then combine those into a net force and torque on each molecule. This can give you a more sophisticated charge distribution than simple spheres (though it will complicate accurate collision detection).
Collisions are more complex but can be simplified if you treat them as a collision of two spheres. At that point it's really a 2D ball collision detection problem - I'll refer to other SE questions for how to handle that.
This can be a really tricky problem when you are dealing with discrete jumps in time (what if the particles would have collided in between those steps?). Funnily enough, it helps that Euler's method is so bad - it will force you to use a small $\Delta t$ so you'll catch more collisions, but you'll still miss some unless you implement some clever methods to handle intermediate collisions.
Once you've implemented these - congratulations! You've coded a simple model of classical particles. There are plenty of ways to improve it:
- This numerical integrator is terrible (see: Runge-Kutta 4 or symplectic integration)
- It won't support too many particles (see: Barnes-Hut or fast multipole methods)
- The collision detection could be dramatically improved (see: particle-particle collisions)
- Include simplified approximations of bonding potential (see: Max's answer, Lennard-Jones potential)
To name a few, but you've taken your first step towards a functional model of classical charged particles. Over time you can refine this simulation, better approximate their interactions and perhaps even go down that quantum rabbit hole.