Cleaned up

#include "forces.h"
#include <assert.h>
#include <math.h>
#include "ndmalloc.h"

void interaction_pairs_alloc(interaction_pairs_t* p, long long Npairsmax)
{
    p->npairs = Npairsmax;
    p->pairi  = ndmalloc(sizeof(int), 1, Npairsmax);
    p->pairj  = ndmalloc(sizeof(int), 1, Npairsmax);
    p->dj     = ndmalloc(sizeof(packed_int), 1, Npairsmax);
}

void interaction_pairs_free(interaction_pairs_t* p)
{
    ndfree(p->pairi);
    ndfree(p->pairj);
    ndfree(p->dj);
}

double computeForces(int N, atom_t atoms[], interaction_pairs_t* p, double L, bool accum)
/* accum: true if the forces need to be added too, false if the forces should be replaced
   returns the potential enenty.*/
{
    /* initialize energy and forces to zero */
    double Usum = 0;
    if (!accum) {
        for (int i = 0; i < N; ++i)
            atoms[i].fx = atoms[i].fy = atoms[i].fz = 0.0;
    }
    /* determine interactions */
    for (size_t k = 0; k < p->npairs; ++k) {
        int i = p->pairi[k];
        int j = p->pairj[k];
        int d = p->dj[k];
        /* determine distance in periodic geometry */
        double dx = atoms[i].rx - atoms[j].rx - L*pack2x(d);
        double dy = atoms[i].ry - atoms[j].ry - L*pack2y(d);
        double dz = atoms[i].rz - atoms[j].rz - L*pack2z(d);
        double r2 = dx*dx + dy*dy + dz*dz;
        if (r2 < rc*rc) {
            double r2i = 1/r2;
            double r6i = r2i*r2i*r2i;
            double fij = 48*r2i*r6i*(r6i-0.5); /* force multiplier */
            double eij = 4*r6i*(r6i-1); /* potential energy between i and j */
            /* within smooth cutoff region? */
            if (r2 > rcp*rcp) {
                /* in out part of potential, between rcp and rc,      */
                /* replace phi by alpha * phi for smoother potential  */
                /* based on a slight rewriting of the alpha factor to */
                /*  alpha=1/2-1/4x(x^2-3)                             */
                /* where                                              */
                /*  x=(2r-rcp-rc)/(rcp-rc)                            */
                double r = sqrt(r2);
                double x = (2*r-rcp-rc)/(rcp-rc);
                double alpha  = 0.5-0.25*x*(x*x-3);
                double dalpha = 1.5*(x*x-1)/(r*(rcp-rc));
                fij = alpha*fij+dalpha*eij;
                eij = alpha*eij;
            }
            Usum += eij;
            atoms[i].fx += fij*dx;
            atoms[i].fy += fij*dy;
            atoms[i].fz += fij*dz;
            atoms[j].fx -= fij*dx;
            atoms[j].fy -= fij*dy;
            atoms[j].fz -= fij*dz;
        }
    }
    return Usum;
}