001: import java.io.*;
002: import JSci.maths.*;
003: import JSci.maths.vectors.AbstractDoubleVector;
004: import JSci.maths.vectors.Double3Vector;
005: import JSci.maths.statistics.*;
006: import JSci.io.*;
007:
008: /**
009: * Monte Carlo calculation of Helium ground state energy.
010: * @author Mark Hale
011: * @version 1.0
012: */
013: public final class MonteCarlo implements Mapping {
014: private int N;
015: private double energy[];
016: private AbstractDoubleVector r1;
017: private AbstractDoubleVector r2;
018:
019: /**
020: * Instantiate class.
021: */
022: public static void main(String arg[]) {
023: if (arg.length == 1) {
024: int n = Integer.valueOf(arg[0]).intValue();
025: new MonteCarlo(n);
026: } else {
027: System.out
028: .println("Need to specify command line arguments:");
029: System.out.println("<number of iterations>");
030: }
031: }
032:
033: /**
034: * Constructor.
035: * @param n number of iterations
036: */
037: public MonteCarlo(int n) {
038: N = n;
039: r1 = new Double3Vector(Math.random(), Math.random(), Math
040: .random());
041: r2 = new Double3Vector(-Math.random(), -Math.random(), -Math
042: .random());
043: energy = new double[N];
044: compute();
045: saveResults();
046: }
047:
048: /**
049: * Compute the ground state energy.
050: */
051: private void compute() {
052: AbstractDoubleVector tmpr1, tmpr2;
053: double prob, tmpprob;
054:
055: // Stabilising
056: for (int i = 0; i < 1000; i++) {
057: tmpr1 = r1.mapComponents(this );
058: tmpr2 = r2.mapComponents(this );
059: tmpprob = trialWF(tmpr1, tmpr2);
060: tmpprob *= tmpprob;
061: prob = trialWF(r1, r2);
062: prob *= prob;
063: if (tmpprob / prob > Math.random()) {
064: r1 = tmpr1;
065: r2 = tmpr2;
066: }
067: }
068:
069: System.out.println("Integrating...");
070: for (int i = 0; i < N; i++) {
071: tmpr1 = r1.mapComponents(this );
072: tmpr2 = r2.mapComponents(this );
073: tmpprob = trialWF(tmpr1, tmpr2);
074: tmpprob *= tmpprob;
075: prob = trialWF(r1, r2);
076: prob *= prob;
077: if (tmpprob / prob > Math.random()) {
078: r1 = tmpr1;
079: r2 = tmpr2;
080: }
081: energy[i] = localEnergy(r1, r2);
082: }
083: }
084:
085: /**
086: * Trial wavefunction.
087: * @param r1 position vector of electron 1
088: * @param r2 position vector of electron 2
089: */
090: private double trialWF(AbstractDoubleVector r1,
091: AbstractDoubleVector r2) {
092: double modR1 = r1.norm();
093: double modR2 = r2.norm();
094: double modR12 = r1.subtract(r2).norm();
095: return Math.exp(-2 * modR1) * Math.exp(-2 * modR2)
096: * Math.exp(modR12 / 2);
097: }
098:
099: /**
100: * Local energy calculation.
101: * @param r1 position vector of electron 1
102: * @param r2 position vector of electron 2
103: */
104: private double localEnergy(AbstractDoubleVector r1,
105: AbstractDoubleVector r2) {
106: AbstractDoubleVector r12 = r2.subtract(r1);
107: double termR112 = r1.scalarProduct(r12)
108: / (r1.norm() * r12.norm());
109: double termR212 = r2.scalarProduct(r12)
110: / (r2.norm() * r12.norm());
111: return -17 / 4 - termR112 + termR212;
112: }
113:
114: /**
115: * Update electron co-ordinates.
116: */
117: public double map(double x) {
118: return x + (Math.random() - 0.5);
119: }
120:
121: /**
122: * Not used, dummy implementation for Mapping interface.
123: */
124: public Complex map(Complex z) {
125: return null;
126: }
127:
128: /**
129: * Log results to disk.
130: */
131: private void saveResults() {
132: File file = new File("results.dat");
133: char buf[] = null;
134: NormalDistribution norm = new NormalDistribution(energy);
135: double data[][] = new double[1][3];
136: double mean = norm.getMean();
137: double var = norm.getVariance();
138: System.out.println("Energy: " + mean + " Var: " + var);
139:
140: data[0][0] = N;
141: data[0][1] = mean;
142: data[0][2] = var;
143: // Read in existing data
144: try {
145: FileReader in = new FileReader(file);
146: buf = new char[(int) file.length()];
147: in.read(buf);
148: in.close();
149: } catch (Exception e) {
150: System.out.println("No previous data - new file.");
151: }
152: // Save all to file
153: try {
154: TextWriter out = new TextWriter(file, ',');
155: if (buf != null)
156: out.write(buf);
157: out.write(data);
158: out.close();
159: } catch (Exception e) {
160: System.out.println("Failed to save.");
161: }
162: }
163: }
|