package etomica.potential;

import etomica.api.IAtomList;
import etomica.api.IBoundary;
import etomica.api.IBox;
import etomica.api.IMolecule;
import etomica.api.IMoleculeList;
import etomica.api.IPotentialAtomic;
import etomica.api.IPotentialMolecular;
import etomica.api.IVector;
import etomica.api.IVectorMutable;
import etomica.atom.DipoleSource;
import etomica.atom.IAtomOriented;
import etomica.space.ISpace;
import etomica.space.Space;
import etomica.space.Tensor;
import etomica.util.numerical.BesselFunction;

public class PotentialCalculationPhiSumHeisenberg implements PotentialCalculationMolecular {
	 protected final IVectorMutable ei,ej;
	 protected IVectorMutable dr;
	 protected double secondDerivativeSum= 0;
	 protected DipoleSource dipoleSource;
	 

	public PotentialCalculationPhiSumHeisenberg(ISpace space) {
	    dr = space.makeVector();
	    ei = space.makeVector();
	    ej = space.makeVector();
	    
	}

	public void doCalculation(IAtomList atoms, IPotentialAtomic potential) {
		if(!(potential instanceof IPotentialAtomicSecondDerivative)){
			return;
		}

		//don't need this for now
//		IPotentialAtomicSecondDerivative potentialSecondDerivative = (IPotentialAtomicSecondDerivative) potential;
//
//		Tensor[] t = potentialSecondDerivative.secondDerivative(atoms);
		
		IAtomOriented atom1 = (IAtomOriented)atoms.getAtom(0);
    	IAtomOriented atom2 = (IAtomOriented)atoms.getAtom(1);
    	ei.E(atom1.getOrientation().getDirection());
    	ej.E(atom2.getOrientation().getDirection());

		double c1 = ei.getX(0);//cost1
		double c2 = ej.getX(0);//cost2
		double s1 = ei.getX(1);//sint1
		double s2 = ej.getX(1);//sint2

		//ij phi_ij = -J*(1+c1*c2/s1/s2); J would be multiplied in meter
		//-J Sin[t1] Sin[t2] (Cos[t1] Cos[t2] + Sin[t1] Sin[t2])
//		secondDerivativeSum += -2.0*(s1*s2+c1*c2)*s1*s2;

		//ii phi_ii = J*s2^3/s1; J would be multiplied in meter
		//J Sin[t1] Sin[t2]
		//jj phi_jj = J*s1^3/s2; J would be multiplied in meter
		//J Sin[t1] Sin[t2]
		//ii and jj is the same
//		secondDerivativeSum += 2.0*s1*s2;

		//or you could combine ij ii and jj
//		secondDerivativeSum += -2.0*(s1*s2+c1*c2-1)*s1*s2;

		double Cos = ei.dot(ej);
		secondDerivativeSum += -2*Cos*Cos+2*Cos;
	}

	public void doCalculation(IMoleculeList molecules, IPotentialMolecular potential) {
		if(!(potential instanceof IPotentialMolecularSecondDerivative)){
			return;
		}
	}
	

	public void zeroSum() {
		secondDerivativeSum = 0.0;
	}

	/**
	 * Returns the current value of the energy sum.
	 */
	public double getSum() {
        return secondDerivativeSum;
    }
	

}

@@ -1081,220 +1081,8 @@
 %09%7D%0A%0A
-%09%09//don't need this for now%0A//%09%09IPotentialAtomicSecondDerivative potentialSecondDerivative = (IPotentialAtomicSecondDerivative) potential;%0A//%0A//%09%09Tensor%5B%5D t = potentialSecondDerivative.secondDerivative(atoms);%0A%09%09
 %0A%09%09I
@@ -1135,20 +1135,17 @@
 tom(0);%0A
-    
+%09
 %09IAtomOr
@@ -1192,20 +1192,17 @@
 tom(1);%0A
-    
+%09
 %09ei.E(at
@@ -1239,20 +1239,17 @@
 ion());%0A
-    
+%09
 %09ej.E(at
@@ -1293,707 +1293,788 @@
 %0A%0A%09%09
-double c1 = ei.getX(0);//cost1%0A%09%09double c2 = ej.getX(0);//cost2%0A%09%09double s1 = ei.getX(1);//sint1%0A%09%09double s2 = ej.getX(1);//sint2%0A%0A%09%09//ij phi_ij = -J*(1+c1*c2/s1/s2); J would be multiplied in meter%0A%09%09//-J Sin%5Bt1%5D Sin%5Bt2%5D (Cos%5Bt1%5D Cos%5Bt2%5D + Sin%5Bt1%5D Sin%5Bt2%5D)%0A//%09%09secondDerivativeSum += -2.0*(s1*s2+c1*c2)*s1*s2;%0A%0A%09%09//ii phi_ii = J*s2%5E3/s1; J would be multiplied in meter%0A%09%09//J Sin%5Bt1%5D Sin%5Bt2%5D%0A%09%09//jj phi_jj = J*s1%5E3/s2; J would be multiplied in meter%0A%09%09//J Sin%5Bt1%5D Sin%5Bt2%5D%0A%09%09//ii and jj is the same%0A//%09%09secondDerivativeSum += 2.0*s1*s2;%0A%0A%09%09//or you could combine ij ii and jj%0A//%09%09secondDerivativeSum += -2.0*(s1*s2+c1*c2-1)*s1*s2;%0A%0A%09%09double Cos = ei.dot(ej);%0A%09%09secondDerivativeSum += -2*Cos*Cos+2*Cos;
+//complicated way to get secondDerivativeSum, don't need the secondDerivative from p2Spin%0A%09%09IPotentialAtomicSecondDerivative potentialSecondDerivative = (IPotentialAtomicSecondDerivative) potential;%0A%09%09Tensor%5B%5D t = potentialSecondDerivative.secondDerivative(atoms);%0A%09%09secondDerivativeSum += 2*t%5B0%5D.component(0,0)*(ei.dot(ej));%0A%09%09secondDerivativeSum += t%5B1%5D.component(0,0);%0A%09%09secondDerivativeSum += t%5B2%5D.component(0,0);%0A%0A%09%09//much easier way but need to time back coupling J in the meter and here!!!!!!!!%0A//%09%09double Cos = ei.dot(ej);%0A//%09%09secondDerivativeSum += -2*Cos*Cos+2*Cos;%0A//%09%09double diff = 2*t%5B0%5D.component(0,0)*(ei.dot(ej))+t%5B1%5D.component(0,0)+t%5B2%5D.component(0,0)-1.5*( -2*Cos*Cos+2*Cos);%0A//%09%09System.out.println(%22diff = %22 + diff);//check these two approach is the same or not
 %0A%09%7D%0A

simplify the code and add a second method to compute secondDerivative