Molecular Dynamics simulations of Silicon-Energy Problems

[danielakkerma]

Hello everyone!
I have constructed a model of a silicon crystal, consisting of 64 atoms, for my various needs of testing physical phenomena in that substance(Silicon).
At present, rather annoyingly, I am still stuck at the basic operational level, that is, calibrating the system to an expected equilibrium, then releasing it, using a specified potential, and (fingers-crossed), hoping that it would maintain that structure, but most of all, would at least conserve fundamental physical properties, most notably, the total Energy; however, it doesn't work!
My problem is that the energy starts to spike up for no apparent reason, and that the whole system therefore starts to melt unexpectedly, without any external input of energy of any kind!
***
I have taken as my potential the Stillinger-Weber one, personally derived the formulae for the force(including the 3-body-part), made sure all quantities were physically sound, and then applied it.
I proceeded then to allocate velocities using the Maxwell distribution, and to remove, thereafter, the center-of-mass drift.
I assigned coordinates using what every novice in Material/Solid State/Nano physics knows:
A similar rule can be found here.
Finally, I tucked the matter neatly into a Java program I had written for it, that integrates the equations of Newtonian motion via the Predictor-Corrector methods that are well-known, particularly Gear's approach; Then hoped for the best.
What am I missing? What have I done wrong?
Most sincerely thankful,
Beholden for your attention,
And reliant on your aid,
Daniel
P.S.
I've attached a text file denoting the evolution of time, Kinetic Energy, Potential Energy in that order. My time step is extremely small, 1 fs, a Femtosecond, but that file only includes data from every 10th iteration, to save space.;
Another is the set of coordinates of the atoms, in Angstroms(10^-10 meters).
The temperature was set to 300° Kelvin at the beginning of the simulation. You can see how rapidly the system swings out of control!
All other constants can be found in the links I've provided.
If there's anything else I should, or need to post, please let me know!

 

[danielakkerma]

Folks, I really hate to bump this thread, but perhaps I am posting in the wrong board? Maybe the moderators could direct me to which niche this belongs, and move it there?
Thanks, and sorry for boosting what I fear is a very boring question,
:D!
Daniel

 

[NemoReally]

I'm afraid it's difficult to say what the problem could be without looking at the code - there are a myriad of ways in which a programming error could create difficulties. Swinging out of control, however, is sometimes a symptom of an error in the integrator, either in its implementation or with the step size versus rate of change of the equations - you might be better of trying an adaptive DE solver.

 

[danielakkerma]

Hey, thanks for a very prompt and instructive reply!
Could posting the integration scheme help? (this is the third-order Gear Predictor-Corrector algorithm) Anyway, it goes as follows:

while(i < Mesh.getCount())
{
c = Mesh.get(i);
	c.r.Add(Vec.Dot(dt, c.v));
	c.r.Add(Vec.Dot(dt*dt/(2.0d), c.a));
	c.r.Add(Vec.Dot(dt*dt*dt/(6.0d), c.ap));
	c.v.Add(Vec.Dot(dt, c.a));
	c.v.Add(Vec.Dot(dt*dt/2.0d, c.ap));
	c.a.Add(Vec.Dot(dt, c.ap));
	c.apg = c.a;
i++;
}

Here, we have Mesh, as a linked list containing Particles(of which type is the variable "c").
Each particle contains vectors r, v, a, ap, apg, that correspond to the various vectors needed for this operation.
Static functions, belonging to the class(Vec), i.e., Add, Dot, carry out the vectorial addition and dot product, respectively.
After "predicting" the expected positions and accelerations, we compute the actual value for each c.a(that is, each particle's acceleration, by employing a function- "Forces"-, that collects and sums the particle's interactions using the selected potential, Stillinger-Weber).

i = 0; 
while(i < Mesh.getCount())
		{
			Forces(Mesh.get(i), i);
                        i++;
		}

Then we "correct" our erstwhile hunch re the future positions, by comparing the conjectured acceleration with the actual one:

		i = 0; 
		while(i < Mesh.getCount())
		{
		c = Mesh.get(i);
			da = Vec.Add(c.a, Vec.Dot(-1.0d, c.apg));
			c.r.Add(Vec.Dot(3.0d*dt*dt/(16.0d*d_c), da));
			c.v.Add(Vec.Dot(251.0d*dt/360.d, da));
                        c.ap.Add(Vec.Dot(11.0d/(dt*18.0d), da));
               i++;
}

All units are fully SI, meaning, dt = 10^-15 seconds, the distances likewise in meters.
Where could I have gone wrong?
Yours most thankful!
Daniel

 

[LudusRex]

Dear Daniel,

Thank you for sharing your work and concerns with us. As a fellow scientist, I understand the frustration of encountering unexpected results in simulations.

Firstly, I commend you on your approach and the level of detail you have provided. It is clear that you have put a lot of effort into your model and simulations.

From what you have described, it appears that your model may not be accurately representing the physical properties of silicon. The Stillinger-Weber potential is a well-established potential for silicon, but it is possible that your implementation or calculations may have errors. I would recommend double-checking your calculations and equations, and perhaps comparing your results with those of other studies that have used the same potential.

Another possibility is that your initial conditions may not be appropriate for the system you are simulating. The Maxwell distribution is commonly used for assigning velocities, but it may not be suitable for your specific system. You may want to try using a different initial temperature or velocity distribution and see if that affects the stability of your simulation.

Additionally, it is important to consider the size of your system and the time scale of your simulation. A 64-atom system may be too small to accurately represent the behavior of a bulk material like silicon. You may want to consider increasing the size of your system or running the simulation for a longer period of time to observe any long-term trends.

Overall, it is difficult to pinpoint the exact issue without more information or access to your code. I would suggest seeking the help of colleagues or mentors who have experience with molecular dynamics simulations and can provide more specific guidance.

Best of luck with your research and I hope you are able to resolve the issues you are encountering. Keep up the good work!

Sincerely,