package etomica.potential.ewald;

import etomica.potential.Potential2Soft;
import etomica.potential.Potential2SoftSpherical;
import etomica.potential.TruncationFactory;
import etomica.space3d.Space3D;

import static etomica.math.SpecialFunctions.erfc;

public class P2Ewald1Real extends Potential2SoftSpherical {

    public static Potential2Soft makeTruncated(double qiqj, double alpha, TruncationFactory tf) {
        return tf.make(new P2Ewald1Real(qiqj, alpha));
    }

    private double qiqj;
    private final double alpha;

    private static final double SQRT_PI = Math.sqrt(Math.PI);

    public P2Ewald1Real(double qiqj, double alpha) {
        super(Space3D.getInstance());
        this.qiqj = qiqj;
        this.alpha = alpha;
    }

    public void setQiQj(double newQQ) {
        qiqj = newQQ;
    }

    public double getQiQj() {
        return qiqj;
    }

    @Override
    public double du(double r2) {
        double r = Math.sqrt(r2);
        return -qiqj * (2 * SQRT_PI * Math.exp(-alpha * alpha * r2) * alpha + erfc(alpha * r) / r);
    }

    @Override
    public void u012add(double r2, double[] u012) {
        double r = Math.sqrt(r2);
        double u = qiqj * erfc(alpha * r) / r;
        double derfc = -qiqj * 2 / SQRT_PI * Math.exp(-alpha * alpha * r2) * alpha;
        u012[0] += u;
        u012[1] += derfc - u;
        u012[2] += -derfc * (2 + alpha * alpha * 2 * r) + 2 * u;
    }

    @Override
    public double d2u(double r2) {
        double r = Math.sqrt(r2);
        return -qiqj * (2 * SQRT_PI * Math.exp(-alpha*alpha*r2) *(alpha*(1 - alpha*alpha*2*r)) + erfc(alpha*r)/r);
    }

    @Override
    public double integral(double rC) {
        return 0;
    }

    @Override
    public double u(double r2) {
        double r = Math.sqrt(r2);
        return qiqj*erfc(alpha*r)/r;
    }
}

@@ -1393,16 +1393,17 @@
  * 2 * r
+2
 ) + 2 * 
@@ -1401,24 +1401,24 @@
 2) + 2 * u;%0A
-
     %7D%0A%0A    @
@@ -1495,38 +1495,93 @@
 rt(r2);%0A        
-return
+double u = qiqj * erfc(alpha * r) / r;%0A        double derfc =
  -qiqj * (2 * SQ
@@ -1565,36 +1565,35 @@
 derfc = -qiqj * 
-(
 2 
-*
+/
  SQRT_PI * Math.
@@ -1606,68 +1606,92 @@
 lpha
-*
+ * 
 alpha
-*
+ * 
 r2) *
-(
+ 
 alpha
-*(1 -
+;%0A        return -derfc * (2 +
  alpha
-*
+ * 
 alpha
-*2*r)) + erfc(alpha*r)/r)
+ * 2 * r2) + 2 * u
 ;%0A  

Fix second derivative