Search code examples
c#simulationphysics

Molecular MC simulation is not equilibrating

Suppose, a polymer has N monomers in its chain. I want to simulate its movement using the bead-spring model. However, there was no periodic boundary condition applied. Coz, points were generated such that they never cross the boundary.

So, I wrote the following program.

I am using 1 million steps. The energy is not fluctuating as expected. After several thousand steps the curve goes totally flat.

enter image description here

enter image description here

The X-axis is steps. Y-axis is total energy.

Can anyone check the source code and tell me what I should change?

N.B. I am especially concerned with the function that calculates the total energy of the polymer.

Probably, the algorithm is incorrect.

    public double GetTotalPotential()
    {
        double totalBeadPotential = 0.0;
        double totalSpringPotential = 0.0;

        // calculate total bead-energy
        for (int i = 0; i < beadsList.Count; i++)
        {
            Bead item_i = beadsList[i];
            Bead item_i_plus_1 = null;

            try
            {
                item_i_plus_1 = beadsList[i + 1];

                if (i != beadsList.Count - 1)
                {
                    // calculate total spring energy.
                    totalSpringPotential += item_i.GetHarmonicPotential(item_i_plus_1);
                }
            }
            catch { }

            for (int j = 0; j < beadsList.Count; j++)
            {
                if (i != j)
                {
                    Bead item_j = beadsList[j];
                    totalBeadPotential += item_i.GetPairPotential(item_j);
                    //Console.Write(totalBeadPotential + "\n");
                    //Thread.Sleep(100);
                }
            }
        }

        return totalBeadPotential + totalSpringPotential;
    }
Solution

  • Problem of this application is that simulations (Simulation.SimulateMotion) are run in separate thread in parallel to the draw timer (SimulationGuiForm.timer1_Tick) and share the same state (polymerChain) without any sync/signaling, so some mutations of polymerChain are skipped completely (not drawn) and when the simulation is finished (far before the finish of the drawing) the timer1_Tick will redraw the same polymerChain. You can easily check that by adding counter to Simulation and increasing it in the SimulateMotion:

    public class Simulation
{
    public static int Simulations = 0; // counter
    public static void SimulateMotion(PolymerChain polymerChain, int iterations)
    {
        Random random = new Random();
        
        for (int i = 0; i < iterations; i++)
        {
            Simulations++; // bump the counter
            // rest of the code
            // ...

    And checking it in timer1_Tick:

    private void timer1_Tick(object sender, EventArgs e)
{
    // ...
    // previous code
    if (Simulation.Simulations == totalIterations)
    {
        // breakpoint or Console.Writeline() ...
        // will be hit as soon as "the curve goes totally flat" 
    }

    DrawZGraph();            
}

    You need to rewrite your application in such way that SimulateMotion either stores iterations in some collection which is consumed by timer1_Tick (basically implementing producer-consumer pattern, for example you can try using BlockingCollection, like I do in the pull request) or performs it's actions only when the current state is rendered.