import java.util.Arrays; import java.util.Scanner; public class Main { public static void main(String[] args) throws CloneNotSupportedException { new Main().run(); } void run() throws CloneNotSupportedException { solver(); } void solver() throws CloneNotSupportedException { Scanner sc = new Scanner(System.in); double potential = 0; double balancer = 1; double volume = 1; // 単位構造は(1,0,0),(0,1,0),(0,0,1) // 逆格子ベクトルは(1,0,0),(0,1,0),(0,0,1) double[][] rt = new double[16][]; for (int i = 0; i < 8; i++) { rt[i] = new double[] { (i >> 2) % 2 / 2.0, (i >> 1) % 2 / 2.0, i % 2 / 2.0 }; } for (int i = 0; i < 8; i++) { rt[i + 8] = new double[] { (i >> 2) % 2 / 2.0 + 1.0 / 4.0, (i >> 1) % 2 /2.0 + 1.0 / 4.0, i % 2 / 2.0 + 1.0 / 4.0 }; } double[] qt = new double[16]; for (int i = 0; i < 16; i++) { qt[i] = sc.nextDouble(); } potential += shortRangeInteraction(rt, qt, balancer); potential += longRangeInteracion(rt, qt, balancer, 1).real; System.out.println(potential); } double shortRangeInteraction(double[][] rt, double[] qt, double balancer) { double ret = 0; int L = 30; for (int x = -L; x <= L; x++) { for (int y = -L; y <= L; y++) { for (int z = -L; z <= L; z++) { for (int i = 0; i < qt.length; i++) { double[] nrt = new double[] { rt[i][0] + x, rt[i][1] + y, rt[i][2] + z }; if (nrt[0] == 0 && nrt[1] == 0 && nrt[2] == 0) continue; ret += qt[i] / norm(nrt) * erfc(Math.sqrt(balancer) * norm(nrt)); } } } } return ret; } Complex longRangeInteracion(double[][] rt, double[] qt, double balancer, double volume) throws CloneNotSupportedException { Complex ret = new Complex(0, 0); int L = 30; for (int x = 0; x < L; x++) { for (int y = 0; y < L; y++) { for (int z = 0; z < L; z++) { if (x == 0 && y == 0 && z == 0) continue; // Gは逆格子ベクトル double[] G = new double[] { 2 * Math.PI * x, 2 * Math.PI * y, 2 * Math.PI * z }; // r2=G^2 double r2 = 4 * Math.PI * Math.PI * (x * x + y * y + z * z); ret.add(structuralFactor(qt, G, rt).multiply(1 / r2 * Math.exp(-r2 / (4 * balancer)))); } } } ret = ret.multiply(4 * Math.PI / volume); // consider the effect of the self potential ret.add(new Complex(-2 * qt[0] * Math.sqrt(balancer / Math.PI), 0)); return ret; } double norm(double[] vector) { double ret = 0; for (int i = 0; i < vector.length; i++) { ret += vector[i] * vector[i]; } ret = Math.sqrt(ret); return ret; } Complex structuralFactor(double[] qt, double[] G, double[][] rt) throws CloneNotSupportedException { Complex ret = new Complex(0, 0); for (int i = 0; i < qt.length; i++) { ret.add(exp(dotProduct(G, rt[i])).multiply(qt[i])); } return ret; } // =exp(-i(theta)) Complex exp(double theta) { return new Complex(-Math.cos(theta), -Math.sin(theta)); } double dotProduct(double[] v1, double[] v2) { double ret = 0; for (int i = 0; i < v1.length; i++) { ret = v1[i] * v2[i]; } return ret; } // fractional error in math formula less than 1.2 * 10 ^ -7. // although subject to catastrophic cancellation when z in very close to 0 // from Chebyshev fitting formula for erf(z) from Numerical Recipes, 6.2 public static double erf(double z) { double t = 1.0 / (1.0 + 0.5 * Math.abs(z)); // use Horner's method double ans = 1 - t * Math.exp(-z * z - 1.26551223 + t * (1.00002368 + t * (0.37409196 + t * (0.09678418 + t * (-0.18628806 + t * (0.27886807 + t * (-1.13520398 + t * (1.48851587 + t * (-0.82215223 + t * (0.17087277)))))))))); if (z >= 0) return ans; else return -ans; } public static double erfc(double z) { return 1 - erf(z); } class Complex implements Cloneable { double real; double imaginal; public Complex(double real, double imaginal) { this.real = real; this.imaginal = imaginal; } void add(Complex c) { this.real += c.real; this.imaginal += c.imaginal; } void show() { System.out.println("real:" + real + " comples:" + imaginal); } Complex multiply(double a) throws CloneNotSupportedException { this.real *= a; this.imaginal *= a; return (Complex) this.clone(); } } Complex add(Complex... c) { Complex ret = new Complex(0, 0); for (int i = 0; i < c.length; i++) { ret.add(c[i]); } return ret; } void tr(Object... o) { System.out.println(Arrays.deepToString(o)); } }