// Java code for WCA dimer

import java.lang.Math;

public class WCADimer implements Potential {

  private WCAPotential system; // WCA fluid solvent

  private int natoms; // number of atoms
  private double r0; // compact state separation
  private double h; // barrier height
  private double w; // second minimum is at r0 + 2*w

  private double box_edge_length;
  private double half_box_edge;
  
  /*
   * Create a WCA dimer system.
   *
   * @param natoms the number of atoms in the system
   */
  public WCADimer(int natoms, double sigma, double epsilon, double box_edge_length, double barrier_height) {
    // Create fluid.
    this.system = new WCAPotential(natoms, sigma, epsilon, box_edge_length);
    
    this.natoms = natoms;
    this.box_edge_length = box_edge_length;
    this.half_box_edge = this.box_edge_length / 2.0;

    // 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 imagedSeparationVector(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];

      // Image into box.
      while (dr < - this.half_box_edge)
        dr += this.box_edge_length;
      while (dr > this.half_box_edge)
        dr -= this.box_edge_length;

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

      // Accumulate square distance.
      r2 += dr*dr;
    }
    
    return r2;
  }
    
  public double computePotential(double [] positions) {
    double potential = this.system.computePotential(positions); // WCA fluid energy

    // Compute rij = r_j - r_i
    double [] rij = new double[3];
    double r2 = imagedSeparationVector(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);
    potential += this.h * y*y;
      
    return potential;
  }

  public double computeGradient(double [] gradient, double [] positions) {
    double potential = this.system.computeGradient(gradient, positions); // WCA fluid gradient

    // Compute rij = r_j - r_i
    double [] rij = new double[3];
    double r2 = imagedSeparationVector(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);
    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;
  }

}