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
Writing an article is easy - try our reStructured Text demo

# Lennard-Jones potential Molecular dynamics monte-carlo 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.
 Language C
# '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}125126void 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 data149	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	}157158	x[0] *= 0.8f;159//	x[1] *= 1.1f;160	161	// Initialize the accelerations162	for(int j = 0; j < particleCount; ++j)163	{164		if((j > 0) && (j < particleCount - 1)) // We are in the middle somewhere165		{166			F[j] = ljForce(x[j], x[j + 1]) + ljForce(x[j], x[j - 1]);167		}168		if(j == 0) // We are at the beginning169		{170			F[j] = ljForce(x[0], x[1])  + ljForce(x[0], x[particleCount - 1] + length);171		}172		if(j == (particleCount - 1)) //At the end173		{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;184185	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		}195196		for(int j = 0; j < particleCount; ++j)197		{198			if((j > 0) && (j < particleCount - 1)) // We are in the middle somewhere199			{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 beginning204			{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 end209			{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);					//ok227		u_dist->Fill(U[0]);		//ok228		229		k_data[i] = getKineticE(v, particleCount);230		u_data[i] = getPotentialE(x, particleCount);231		t_data[i] = i;			// ok232		x_data[i] = x[0];		// ok233		e_data[i] = k_data[i] + u_data[i];	// Total energy234	}235	236	//for(int i = 0; i < eventCount; ++i)237		//cout<<u_data[i]<<endl;238}239240double 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}249250double 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 somewhere259		{260			U += ljPotential(x[j], x[j + 1]); // - ljPotential(x[j], x[j - 1]);261		}262		if(j == 0) // We are at the beginning263		{264			U += ljPotential(x[0], x[1]); // - ljPotential(x[0], x[particleCount - 1] - length);265		}266		if(j == (particleCount - 1)) //At the end267		{268			U += ljPotential(x[j], x[0] + length); // - ljPotential(x[j], x[j - 1]);269		}270	}271272	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.