License Apache License, 2.0
Lines 253
Keywords
Lennard-Jones (1) molecular dynamics (1) monte carlo (1) ROOT (1) simulation (1)
Permissions
Owner: Stou S.
Viewable by Everyone
Editable by All Siafoo Users
Hide
Siafoo is here to make coding less frustrating and to save you time. Join Siafoo Now or Learn More

Lennard-Jones potential Molecular dynamics monte-carlo Atom Feed 0

In Brief This is some sort of monte-carlo simulation of 10 particles in a Lennard-Jones potential that I wrote for what appears to be thermodynamics. It uses the ROOT framework so good luck getting it to compile.
# 's
  1//////////////////////////
2// Stou Sandalski
3// Physics 112: Fall 2005
4// Computer Problem 1:
5// Molecular simulator
6
7
8#include <cmath>
9
10// Some constants to make life easier
11const double epsilon = 1.0e-20; // J // d
12const double sigma = 3.0e-10; // m // d
13const double mass = 1.0e-26; // kg // d
14const double delta = 1.0e-9;
15
16double *x_data = NULL; // X value for one particle per time
17double *u_data = NULL; // Total potential energy per time
18double *k_data = NULL; // Total potential energy per time
19double *p_data = NULL; // Momentum per time
20double *e_data = NULL; // Total Energy
21double *t_data = NULL; // Time data
22
23double keTotal(double v)
24{
25 return (0.5f)*mass*(v*v);
26}
27
28double ljPotential(double x1, double x2)
29{
30 double ra = x1 - x2;
31
32// if(ra < 0)
33// ra *= -1;
34
35 return ljPotential(ra);
36}
37
38double ljPotential(double r)
39{
40 if(r == 0)
41 {
42 cout<<"Division by zero in potential calculation: r= "<<r<<endl;
43 return 0;
44 }
45
46 double rt = epsilon*(pow((sigma/r), 12) - pow((sigma/r),6));
47 return rt;
48}
49
50double ljForce(double x1, double x2)
51{
52 double ra = x1 - x2;
53
54 if(ra < 0)
55 ra *= -1;
56
57 return ljForce(ra);
58}
59
60double ljForce(double r)
61{
62 if(r == 0)
63 {
64 cout<<"Division by zero in force calculation: r= "<<r<<endl;
65 return 0;
66 }
67
68 double rt = (epsilon/sigma)*(pow((sigma/r), 13) - pow((sigma/r),7));
69 return rt;
70}
71
72
73void runSim(int eventCount, int particleCount)
74{
75 x_data = new double[eventCount]; // X value for one particle
76 u_data = new double[eventCount]; // Total potential energy for a given time
77 k_data = new double[eventCount]; // Total Kinetic energy for a given time
78 p_data = new double[eventCount]; // Momentum
79 e_data = new double[eventCount]; // Total Energy
80 t_data = new double[eventCount]; // Time data
81
82 TH1D *p_dist = new TH1D("p_dist", "Momentum Distribution", 100, -1.0e-36, 1.0e-36);
83 TH1D *u_dist = new TH1D("u_dist", "Potential Energy Distribution", 100, -1.0e-36, 1.0e-36);
84
85 TFile homework1("problem1.root", "RECREATE");
86
87 runSim(eventCount, particleCount, p_dist, u_dist);
88
89 TGraph *u_graph = new TGraph(eventCount, t_data, u_data);
90 TGraph *k_graph = new TGraph(eventCount, t_data, k_data);
91 TGraph *e_graph = new TGraph(eventCount, t_data, e_data);
92 TGraph *p_graph = new TGraph(eventCount, t_data, p_data);
93 TGraph *x_graph = new TGraph(eventCount, t_data, x_data);
94 TGraph *xp_graph = new TGraph(eventCount, p_data, x_data);
95
96 u_graph->Write("u_Graph");
97 k_graph->Write("k_Graph");
98 e_graph->Write("e_Graph");
99 p_graph->Write("p_Graph");
100 x_graph->Write("x_Graph");
101 xp_graph->Write("xp_Graph");
102
103 p_dist->Write();
104 u_dist->Write();
105
106 delete p_dist;
107 delete u_dist;
108
109 delete u_graph;
110 delete k_graph;
111 delete e_graph;
112 delete p_graph;
113 delete x_graph;
114 delete xp_graph;
115
116 delete[] x_data;
117 delete[] u_data;
118 delete[] k_data;
119 delete[] p_data;
120 delete[] e_data;
121 delete[] t_data;
122
123 homework1.Close();
124}
125
126void runSim(int eventCount, int particleCount, TH1D *p_dist, TH1D *u_dist)
127{
128 particleCount = 2;
129 double length = particleCount*pow(2, (1.0f/6.0f))*3.0e-10;
130 double x[2];
131 double v_half[2];
132 double v[2];
133 double a[2];
134 double F[2];
135 double M[2];
136 double U[2];
137
138 for(int i = 0; i < eventCount; ++i)
139 {
140 x_data[i] = 0;
141 u_data[i] = 0;
142 p_data[i] = 0;
143 k_data[i] = 0;
144 e_data[i] = 0;
145 t_data[i] = 0;
146 }
147
148 // Initialize the data
149 for(int j = 0; j < particleCount; ++j)
150 {
151 x[j] = j*pow(2, (1.0f/6.0f))*3.0e-10;
152 v[j] = 0;
153 a[j] = 0;
154 F[j] = 0;
155 M[j] = 0;
156 }
157
158 x[0] *= 0.8f;
159// x[1] *= 1.1f;
160
161 // Initialize the accelerations
162 for(int j = 0; j < particleCount; ++j)
163 {
164 if((j > 0) && (j < particleCount - 1)) // We are in the middle somewhere
165 {
166 F[j] = ljForce(x[j], x[j + 1]) + ljForce(x[j], x[j - 1]);
167 }
168 if(j == 0) // We are at the beginning
169 {
170 F[j] = ljForce(x[0], x[1]) + ljForce(x[0], x[particleCount - 1] + length);
171 }
172 if(j == (particleCount - 1)) //At the end
173 {
174 F[j] = ljForce(x[j], length) + ljForce(x[j], x[j - 1]);
175 }
176 }
177
178 a[0] = F[0]/mass;
179
180 cout<<"Initial Fs "<<F[0]<<" , "<<F[1]<<endl;
181 cout<<"Initial Xs "<<x[0]<<" , "<<x[1]<<endl;
182 cout<<"Initial Vs "<<v[0]<<" , "<<v[1]<<endl;
183 cout<<"Initial As "<<a[0]<<" , "<<a[1]<<endl;
184
185 double temp = 0;
186
187 for(int i = 0; i < eventCount; ++i)
188 {
189
190 for(int j = 0; j < particleCount; ++j)
191 {
192 x[j] = x[j] + v[j]*delta*((double)i) + (0.5f)*a[j]*(pow(delta*((double)i), 2));
193 v_half[j] = v[j] + ((0.5f)*a[j]*delta*((double)i));
194 }
195
196 for(int j = 0; j < particleCount; ++j)
197 {
198 if((j > 0) && (j < particleCount - 1)) // We are in the middle somewhere
199 {
200 F[j] = ljForce(x[j], x[j + 1]) + ljForce(x[j], x[j - 1]);
201 U[j] = ljPotential(x[j], x[j + 1]); // + ljPotential(x[j], x[j - 1]);
202 }
203 if(j == 0) // We are at the beginning
204 {
205 F[j] = ljForce(x[0], x[1]) + ljForce(x[0], x[particleCount - 1] + length);
206 U[j] = ljPotential(x[0], x[1]); // + ljPotential(x[0], x_image[1]);
207 }
208 if(j == (particleCount - 1)) //At the end
209 {
210 F[j] = ljForce(x[j], x[0] + length) + ljForce(x[j], x[j - 1]);
211 U[j] = ljPotential(x[j], x[0] + length); // + ljPotential(x[j], x[j - 1]);
212 }
213
214 }
215
216 for(int j = 0; j < particleCount; ++j)
217 {
218 a[j] = F[j]/mass;
219 v[j] = v_half[j] + (0.5f)*a[j]*delta*((double)i);
220 }
221
222 cout<<v[1]<<endl;
223 // cout<<a[0]<<endl;
224// cout<<v_half<<endl;
225
226 p_dist->Fill(v[0]*mass); //ok
227 u_dist->Fill(U[0]); //ok
228
229 k_data[i] = getKineticE(v, particleCount);
230 u_data[i] = getPotentialE(x, particleCount);
231 t_data[i] = i; // ok
232 x_data[i] = x[0]; // ok
233 e_data[i] = k_data[i] + u_data[i]; // Total energy
234 }
235
236 //for(int i = 0; i < eventCount; ++i)
237 //cout<<u_data[i]<<endl;
238}
239
240double getKineticE(double *v, int particleCount)
241{
242 double T = 0;
243
244 for(int i = 0; i < particleCount; ++i)
245 T += (0.5f)*mass*(v[i]*v[i]);
246
247 return T;
248}
249
250double getPotentialE(double *x, int particleCount)
251{
252 double length = particleCount*pow(2, (1.0f/6.0f))*3.0e-10;
253
254 double U = 0;
255
256 for(int j = 0; j < particleCount; ++j)
257 {
258 if((j > 0) && (j < particleCount - 1)) // We are in the middle somewhere
259 {
260 U += ljPotential(x[j], x[j + 1]); // - ljPotential(x[j], x[j - 1]);
261 }
262 if(j == 0) // We are at the beginning
263 {
264 U += ljPotential(x[0], x[1]); // - ljPotential(x[0], x[particleCount - 1] - length);
265 }
266 if(j == (particleCount - 1)) //At the end
267 {
268 U += ljPotential(x[j], x[0] + length); // - ljPotential(x[j], x[j - 1]);
269 }
270 }
271
272 return U;
273}

This is some sort of monte-carlo simulation of 10 particles in a Lennard-Jones potential that I wrote for what appears to be thermodynamics. It uses the ROOT framework so good luck getting it to compile.