// Java code for WCA dimer in vacuum

import java.lang.Math;

public class WCADimerVacuum implements Potential {

  private double r0; // compact state separation
  private double h; // barrier height
  private double w; // second minimum is at r0 + 2*w
  
  /*
   * Create a WCA dimer system.
   *
   * @param natoms the number of atoms in the system
   */
  public WCADimerVacuum(double sigma, double epsilon, double barrier_height) {
    // Store parameters for potential.
    this.r0 = Math.pow(2.0, (1.0/6.0)) * sigma;
    this.w = 0.5 * this.r0;
    this.h = barrier_height;
  }

  private double separationVector(double [] rij, double [] positions, int i, int j) {
    // Compute squared distance.
    double r2 = 0.0;

    for(int k = 0; k < 3; k++) {
      // Compute interparticle separation.
      double dr = positions[3*j+k] - positions[3*i+k];

      // Store interparticle separation.
      rij[k] = dr;

      // Accumulate square distance.
      r2 += dr*dr;
    }
    
    return r2;
  }
    
  public double computePotential(double [] positions) {
    // Compute rij = r_j - r_i
    double [] rij = new double[3];
    double r2 = separationVector(rij, positions, 0, 1);
    double r = Math.sqrt(r2);

    // Dimer potential
    // h*(1-((r-r0-w)/w)^2)^2
    double x = ((r - this.r0 - this.w)/this.w);
    double y = (1.0 - x*x);
    double potential = this.h * y*y;
      
    return potential;
  }

  public double computeGradient(double [] gradient, double [] positions) {
    // Compute rij = r_j - r_i
    double [] rij = new double[3];
    double r2 = separationVector(rij, positions, 0, 1);
    double r = Math.sqrt(r2);

    // Compute unit vector separation.
    double [] nij = new double[3];
    for(int k = 0; k < 3; k++)
      nij[k] = rij[k] / r;               

    // Dimer potential
    // h*(1-((r-r0-w)/w)^2)^2
    double x = ((r - this.r0 - this.w)/this.w);
    double y = (1.0 - x*x);
    double potential = this.h * y*y;
    
    double dxdr = 1.0 / this.w;
    double dydr = -2.0*x*dxdr;
    double dUdr = this.h * 2.0 * y * dydr;

    for (int k = 0; k < 3; k++) {
      gradient[3*0+k] = -dUdr * nij[k];
      gradient[3*1+k] = +dUdr * nij[k];
    }
      
    return potential;
  }

}