KickJava   Java API By Example, From Geeks To Geeks.

Java > Open Source Codes > MonteCarlo


1 import java.io.*;
2 import JSci.maths.*;
3 import JSci.maths.vectors.AbstractDoubleVector;
4 import JSci.maths.vectors.Double3Vector;
5 import JSci.maths.statistics.*;
6 import JSci.io.*;
7
8 /**
9 * Monte Carlo calculation of Helium ground state energy.
10 * @author Mark Hale
11 * @version 1.0
12 */
13 public final class MonteCarlo implements Mapping {
14         private int N;
15         private double energy[];
16     private AbstractDoubleVector r1;
17     private AbstractDoubleVector r2;
18     /**
19     * Instantiate class.
20     */
21     public static void main(String JavaDoc arg[]) {
22             if(arg.length==1) {
23             int n=Integer.valueOf(arg[0]).intValue();
24                 new MonteCarlo(n);
25         } else {
26             System.out.println("Need to specify command line arguments:");
27             System.out.println("<number of iterations>");
28                 }
29     }
30         /**
31         * Constructor.
32         * @param n number of iterations
33         */
34     public MonteCarlo(int n) {
35                 N=n;
36         r1=new Double3Vector(Math.random(),Math.random(),Math.random());
37         r2=new Double3Vector(-Math.random(),-Math.random(),-Math.random());
38                 energy=new double[N];
39         compute();
40                 saveResults();
41     }
42         /**
43         * Compute the ground state energy.
44         */
45     private void compute() {
46                 AbstractDoubleVector tmpr1, tmpr2;
47         double prob,tmpprob;
48
49 // Stabilising
50 for(int i=0;i<1000;i++) {
51             tmpr1=r1.mapComponents(this);
52             tmpr2=r2.mapComponents(this);
53             tmpprob=trialWF(tmpr1,tmpr2);
54                         tmpprob*=tmpprob;
55             prob=trialWF(r1,r2);
56                         prob*=prob;
57             if(tmpprob/prob>Math.random()) {
58                 r1=tmpr1;
59                 r2=tmpr2;
60             }
61         }
62
63                 System.out.println("Integrating...");
64                 for(int i=0;i<N;i++) {
65             tmpr1=r1.mapComponents(this);
66             tmpr2=r2.mapComponents(this);
67             tmpprob=trialWF(tmpr1,tmpr2);
68                         tmpprob*=tmpprob;
69             prob=trialWF(r1,r2);
70                         prob*=prob;
71             if(tmpprob/prob>Math.random()) {
72                 r1=tmpr1;
73                 r2=tmpr2;
74             }
75             energy[i]=localEnergy(r1,r2);
76         }
77         }
78     /**
79     * Trial wavefunction.
80         * @param r1 position vector of electron 1
81         * @param r2 position vector of electron 2
82     */
83     private double trialWF(AbstractDoubleVector r1, AbstractDoubleVector r2) {
84         double modR1=r1.norm();
85         double modR2=r2.norm();
86         double modR12=r1.subtract(r2).norm();
87         return Math.exp(-2*modR1)*Math.exp(-2*modR2)*Math.exp(modR12/2);
88     }
89     /**
90     * Local energy calculation.
91         * @param r1 position vector of electron 1
92         * @param r2 position vector of electron 2
93     */
94     private double localEnergy(AbstractDoubleVector r1, AbstractDoubleVector r2) {
95         AbstractDoubleVector r12=r2.subtract(r1);
96         double termR112=r1.scalarProduct(r12)/(r1.norm()*r12.norm());
97         double termR212=r2.scalarProduct(r12)/(r2.norm()*r12.norm());
98         return -17/4-termR112+termR212;
99     }
100     /**
101     * Update electron co-ordinates.
102     */
103     public double map(double x) {
104         return x+(Math.random()-0.5);
105     }
106     /**
107     * Not used, dummy implementation for Mapping interface.
108     */
109     public Complex map(Complex z) {
110         return null;
111     }
112     /**
113     * Log results to disk.
114     */
115     private void saveResults() {
116                 File file=new File("results.dat");
117                 char buf[]=null;
118                 NormalDistribution norm=new NormalDistribution(energy);
119         double data[][]=new double[1][3];
120                 double mean=norm.getMean();
121                 double var=norm.getVariance();
122                 System.out.println("Energy: "+mean+" Var: "+var);
123
124                 data[0][0]=N;
125                 data[0][1]=mean;
126                 data[0][2]=var;
127 // Read in existing data
128 try {
129                         FileReader in=new FileReader(file);
130                         buf=new char[(int)file.length()];
131                         in.read(buf);
132                         in.close();
133                 } catch(Exception JavaDoc e) {
134                         System.out.println("No previous data - new file.");
135                 }
136 // Save all to file
137 try {
138                         TextWriter out=new TextWriter(file,',');
139                         if(buf!=null)
140                                 out.write(buf);
141                         out.write(data);
142                         out.close();
143                 } catch(Exception JavaDoc e) {
144                         System.out.println("Failed to save.");
145                 }
146     }
147 }
148
Popular Tags