package etomica.potential.ewald;

import etomica.potential.IPotential2;
import etomica.potential.TruncationFactory;

import static etomica.math.SpecialFunctions.factorial;

public class P2Ewald6FourierCancel implements IPotential2 {

    public static IPotential2 makeTruncated(double sigmaI, double epsilonI, double sigmaJ, double epsilonJ, double alpha6, TruncationFactory tf) {
        return tf.make(new P2Ewald6FourierCancel(sigmaI, epsilonI, sigmaJ, epsilonJ, alpha6));
    }

    private final double alpha6Sq;
    private final double alpha66;
    private final double Bij, f6;

    public P2Ewald6FourierCancel(double sigmaI, double epsilonI, double sigmaJ, double epsilonJ, double alpha6) {
        super();
        this.alpha6Sq = alpha6 * alpha6;
        this.alpha66 = alpha6Sq * alpha6Sq * alpha6Sq;
        double localBij = 0;
        for (int k = 0; k <= 6; k++) {
            long ck = factorial(6) / (factorial(6 - k) * factorial(k));
            double bik = 0.25 * Math.pow(sigmaI, k) * Math.sqrt(ck * epsilonI);
            double bjk = 0.25 * Math.pow(sigmaJ, 6-k) * Math.sqrt(ck * epsilonJ);
            localBij += bik * bjk;
        }
        this.Bij = localBij;
        f6 = 4*0.5*(sigmaI+sigmaJ)*Math.sqrt(epsilonI*epsilonJ);
    }

    @Override
    public double u(double r2) {
        double a2 = r2 * alpha6Sq;
        double a4 = a2*a2;
        return (-Bij*alpha66*(1+a2+a4/2)*Math.exp(-a2) + f6)/(a4*a2);
    }

    @Override
    public double du(double r2) {
        double a2 = r2 * alpha6Sq;
        double a4 = a2*a2;
        double a6 = a4*a2;
        double e = Math.exp(-a2);
        // (2a + 2a3) e/a6 - 2(1 + a2 + a4/2) e/a5 - 6(1 + a2 + a4/2) e/a7
        return (Bij * alpha66 *(6 + 6*a2 + 3*a4 + a6)*e - 6*f6)/a6;
    }

    @Override
    public void u012add(double r2, double[] u012) {
        double a2 = r2 * alpha6Sq;
        double a4 = a2 * a2;
        double a6 = a4 * a2;
        double e = Math.exp(-a2);
        u012[0] += (-Bij * alpha66 * (1 + a2 + a4 / 2) * e + f6)/ a6;
        u012[1] += (Bij * alpha66 * (6 + 6 * a2 + 3 * a4 + a6) * e - 6*f6) / a6;
        u012[2] += (-Bij * alpha66 * (42 + 42 * a2 + 21 * a4 + (7 + 2 * a2) * a6) * e + 42*f6) / a6;
    }

    @Override
    public double d2u(double r2) {
        double a2 = r2 * alpha6Sq;
        double a4 = a2*a2;
        double a6 = a4*a2;
        double e = Math.exp(-a2);
        return (-Bij * alpha66 *(42 + 42*a2 + 21*a4 + (7+2*a2)*a6)*e + 42*f6)/a6;
    }

}

@@ -698,25 +698,8 @@
 ) %7B%0A
-        super();%0A
     
@@ -731,24 +731,24 @@
 6 * alpha6;%0A
+
         this
@@ -1181,16 +1181,25 @@
  f6 = 4*
+Math.pow(
 0.5*(sig
@@ -1209,16 +1209,19 @@
 +sigmaJ)
+,6)
 *Math.sq
@@ -1378,14 +1378,8 @@
 urn 
-(-Bij*
 alph
@@ -1383,16 +1383,22 @@
 lpha66*(
+-Bij*(
 1+a2+a4/
@@ -1697,30 +1697,30 @@
  return 
-(Bij * 
 alpha66
+ * (Bij
  *(6 + 6
@@ -1965,39 +1965,39 @@
  u012%5B0%5D += 
-(-Bij * alpha66
+alpha66 * (-Bij
  * (1 + a2 +
@@ -2038,23 +2038,16 @@
 12%5B1%5D +=
- (Bij *
  alpha66
@@ -2048,16 +2048,23 @@
 lpha66 *
+ (Bij *
  (6 + 6 
@@ -2120,31 +2120,31 @@
 2%5B2%5D += 
-(-Bij * alpha66
+alpha66 * (-Bij
  * (42 +
@@ -2400,23 +2400,23 @@
 urn 
-(-Bij * alpha66
+alpha66 * (-Bij
  *(4

fix f6 and usage of alpha66