// Java code for WCA dimer

from java.lang.Math import pow;

double compute_energy(double [] positions, double sigma, double epsilon) {
  // Determine number of atoms from size of vector.
  double natoms = (int)(x.size / 3);

  // Pre-compute squared cutoff.
  double rmin = pow(2., (1./6.)) * sigma;  // cutoff at minimum of potential
  double cutoff2 = rmin*rmin;
  double sigma2 = sigma*sigma;

  // Compute interaction potential.
  double potential = 0.0;
  for(int i = 0; i < natoms; i++) 
    for(int j = 0; j < i; j++) {
      // Compute rij = r_j - r_i
      double rij[3];
      for(int k = 0; k < 3; k++) 
        rij[k] = positions[3*j+k] - positions[3*i+k];
      // Compute squared distance.
      double r2 = 0.0;
      for (int k = 0; k < 3; k++)
        r2 += rij[k] * rij[k];
      // Compute interaction potential.
      if (r2 < cutoff2) {
        // 4 epsilon [ (r/sigma)^12 - (r/sigma)^6 ]
        double x2 = r2 / sigma2;
        double x6 = x2*x2*x2;
        double x12 = x6*x6;
        potential += 4.0 * epsilon * (x12 - x6);
      }
    }
      
  return potential;
}

double compute_force(double [] force, double [] positions, double sigma, double epsilon) {
  // Determine number of atoms from size of vector.
  double natoms = (int)(x.size / 3);

  // Pre-compute squared cutoff.
  double rmin = pow(2., (1./6.)) * sigma;  // cutoff at minimum of potential
  double cutoff2 = rmin*rmin;
  double sigma2 = sigma*sigma;

  // Compute force.
  double potential = 0.0;
  for(int i = 0; i < natoms; i++) 
    for(int j = 0; j < i; j++) {
      // Compute rij = r_j - r_i
      double rij[3];
      for(int k = 0; k < 3; k++) 
        rij[k] = positions[3*j+k] - positions[3*i+k];
      // Compute squared distance.
      double r2 = 0.0;
      for (int k = 0; k < 3; k++)
        r2 += rij[k] * rij[k];
      // Compute interaction potential.
      if (r2 < cutoff2) {
        // 4 epsilon [ (r/sigma)^12 - (r/sigma)^6 ]
        double x2 = r2 / sigma2;
        double x6 = x2*x2*x2;
        double x12 = x6*x6;
        potential += 4.0 * epsilon * (x12 - x6);

        // Force contribution.
        double r = sqrt(r2);
        double x5 = x2*x2*(r/sigma);
        double x11 = x6 * x5;
        double nij[3];
        for(int k = 0; k < 3; k++)
          nij[k] = rij[k] / r;               
        for(int k = 0; k < 3; k++) {
          force[3*i+k] = - 4.0 * epsilon * (12.0 * x11 - 6.0 * x5) * nij[k];
          force[3*j+k] = + 4.0 * epsilon * (12.0 * x11 - 6.0 * x5) * nij[k];
        }
      }
    }
      
  return potential;
}