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