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
Stay up to dateembedded code automagically updates, each snippet and article has a feed 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.