Membrane elongation for large thread number

[a2maytin]

Hi Weria,

I am simulating a closed membrane with 1252 vertices. Running Metropolis Monte Carlo with thread number >15 causes the membrane to elongate from an originally spherical shape to an ellipsoid.

The elongation is greater for greater thread number, for example for -nt 30 the membrane forms a long, pill shape (shown below). The membrane does not noticeably elongate for -nt under 15.

Interestingly, the elongation always seems to occur in the [x,y,z]=[1,1,1] direction.

I suspect that when the number of threads is too high relative to the system size, there becomes a sort of bias with the locks on the vertices (which vertices are getting locked and for how long).

I can mitigate the issue by always using a small number of threads. However, I am working on simulating cell division with DNA-membrane coupling where DNA segregation breaks the rotational symmetry to bias elongation to lie along the direction between chromosomes. I am worried that this problem, even if small, will bias membrane elongation to lie in a direction other than that set up by the chromosomes.

As a temporary fix I am rotating the entire system by a random angle every 1000 steps, but is there a better way to address this issue?

I am using FreeDTS version 2.1. My CPU is an AMD Ryzen Threadripper PRO 3955WX 16-Cores (32 threads).

My input.dts is as follows:

Integrator_Type = MC_Simulation
Set_Steps = 1 5000000
Temperature = 1 0
Min_Max_Lenghts = 1 3
MinfaceAngle = -0.5

; Energy method
EnergyMethod = FreeDTS1.0_FF

; Membrane properties
Kappa = 60 0 0
Edge_Parameters = 5 0 0
VertexArea = 0 0.7 0 0

; Move parameters
VertexPositionIntegrator = MetropolisAlgorithmOpenMP 1 1 0.05
AlexanderMove = MetropolisAlgorithmOpenMP 1

; Output settings
TimeSeriesData_Period = 10000
VisualizationFormat = VTUFileFormat VTU_F 10000
NonbinaryTrajectory = TSI TrajTSI 10000
Restart_Period = 10000

[weria-pezeshkian]

Hi Andrew,

I know why this happens, and it was also expected when the ratio of vertices to threads decreases. There are several possible solutions that you can try (perhaps we should discuss them in a meeting). There is an easy approach, but it could be inefficient when using 30 threads for such a small system.

First, everything related to this issue happens in EvolveVerticesByMetropolisAlgorithmWithOpenMPType1.cpp, specifically within a small block of code (shown below). As you can see, I have two implementation approaches. In the first one, when a thread encounters a locked vertex (or neighbors of a locked one), it moves on to a new vertex. You can comment out this part and uncomment the second implementation, in which the thread waits until the vertex becomes free before proceeding. I have not tested this second method, but I would expect that it does not suffer from this problem (due to the physics of the system). However, algorithmically, it may get stuck in an infinite loop. It is dependent on Backoff_Factor parameter which is defined in SimDef.h.

In principle, you could develop a new method to parallelize the EvolveOneStep(int step) function in the EvolveVerticesByMetropolisAlgorithmWithOpenMPType1 object and find a better approach. For example, one could first define the exact sequence of vertices to be moved and then distribute them among the threads (which is perhaps the best solution). We could then create an EvolveVerticesByMetropolisAlgorithmWithOpenMPType2.cpp implementation (with right credit to you) and add it to the FreeDTS source code, provided you can come up with a solid algorithm.

Anyway, improvement of multithread features and tests is the priority and before end of 2026, I hope I have some strong scheme that we can easily use 30 and above thread numbers.

comment out below

            while (true) {
                int r_vid = m_pState->GetRandomNumberGenerator()->IntRNG(no_surf_v);  // Thread-local RNG
                pvertex = m_pSurfV[r_vid];

                // Try to lock the vertex and its neighbors
                if (pvertex->CheckLockVertex()) {
                    if (!pvertex->CheckLockNeighbourVertex()) {
                        pvertex->UnlockVertex();  // Unlock if neighbors can't be locked
                    } else {
                        break;  // Successfully locked vertex and neighbors, break the loop
                    }
                }
            }

uncomment this part

            /*int r_vid = m_pState->GetRandomNumberGenerator()->IntRNG(no_surf_v); // Thread-local RNG
            pvertex = m_pSurfV[r_vid];
        while (true) {
            if (pvertex->CheckLockVertex()) {
                if (!pvertex->CheckLockNeighbourVertex()) {
                    pvertex->UnlockVertex();  // Unlock the vertex if neighbors can't be locked
                } else {
                    break;  // Successfully locked vertex and neighbors, break the loop
                }
            }
        std::this_thread::sleep_for(std::chrono::nanoseconds(Backoff_Factor));
        }*/