There are two distinct approaches: either treat the rigid body as a set of point masses, and alter the forces to fit the constraints or use the true degrees of freedom and write the Lagrangian in terms of centre of mass + moment of inertia, force + torque. The latter is extremely difficult, as it involves event detection across the entire surface area of the body. Thus in practice point masses with constraints are usual.
The simplest case to start with, which will enable you to use much of your code from checkpoint 1, is the diatomic Lennard Jones molecule.
Conceptually, for the point mass approach we calculate forces on the two particles as usual F1', F2' . Then we evaluate the component of the force along the direction of the bond (F1'-F2').r/2r and subtract this from each atom to give new the forces F1, F2 which are actually implemented in the code.
In practice, this is better done by integrating the equations of motion as usual, then correcting the positions and velocities of the particles. If a diatomic molecule has length l , but after a verlet timestep this grows to |r2-r1| = r12 , we can adjust the bondlength which conserving the centre of mass. We also reset the velocities such that |v2-v1|.r12=0 .
A quick and dirty method based on Verlet is to correct the new positions explicitly, shrinking bonds which end up oversized and zeroing the current velocities. However this does not properly respect time-reversal symmetry (it's like doing an Euler integration) and loses all the benefits of Verlet.
The SHAKE method, originally by Ryckaert, Ciccotti, and Berendsen 1977 does it properly. Its Verlet-specific application is explained in appendix A of: "Constrained dynamics and extraction of normal modes from ab initio molecular dynamics" M.M. Siddick, G.J. Ackland, and C.A. Morrison Journal of Chemical Physics, 125 064707 (2006)
While the "quick and dirty" method can be implemented as an independent method after a force update of all particles, to implement SHAKE you need to modify the Verlet integrator itself, and you need to store both new and previous values of the positions (or the constrained vectors).
Try updating your MD code to include rigid diatomic molecules. Note that this does not work well with fixed boundaries. You can implement hardish boundaries by adding a repulsive potential between boundary and particle: the easiest way is to write a component-wise potential where a0 and b0 are the lengths of the box in x and y direction s respectively.
More complex rigid and semi-rigid shapes can always be described by a set of point particles and constrained lengths. In this case each particle is involved in more than one constraint. In practice, one can repeatedly apply the constraints cycling round until all are satisfied. Again the conservation of angular and linear momentum follows from the pairwise nature of each constraint. The number of time each constraint is applied is the usual molecular dynamics trade off between accuracy and speed.
Try setting up a chain of connected LJ particles (a "polymer") with the bondlength between each and the next constrained. Investigate how many times you need to loop through the constraints to get good energy conservation.
Finally, try setting up a set of squares. To do this you need
constrain only two edges and two diagonals, provided you set up the
starting conditions correctly: a parallelogram is a valid solution to
the constraints problem, but cannot be accessed from the square in 2D
dynamics. This should give you pause for thought about how exactly to
define the necessary constraints.