Implemented improved force-stalling extrusion protocol

Author: mtortora

examples/tutorial_stall.ipynb

@@ -13,7 +13,7 @@
         cuda_module = xp.RawModule(code=cuda_code, options=('--use_fast_math',))
 
     _single_leg_search = cuda_module.get_function('_single_leg_search')
-    _harmonic_boltzmann_filter = cuda_module.get_function('_harmonic_boltzmann_filter')
+    _harmonic_distance_filter = cuda_module.get_function('_harmonic_distance_filter')
 
 except ImportError:
     import numpy as xp
@@ -59,12 +59,17 @@ def update_topology(system, bond_list, local=True, thermalize=False):
         system.state.thermalize_particle_momenta(filter=hoomd.filter.All(), kT=1.0)
 
 
-def boltzmann_criterion(system, current_bond_list, trial_bond_list, step_dist=0.4, rest_dist=0.5, threads_per_block=256):
-    """Apply (3D) Boltzmann criterion to list of attempted (1D) extruder moves, based on harmonic bond potential"""
+def stall_criterion(system, current_bond_list, trial_bond_list, bonded_force_dict, threads_per_block=256):
+    """Apply (3D) tension-based stall criterion to list of attempted (1D) extruder moves, based on harmonic bond potential"""
 
-    mu = xp.float64(rest_dist)
-    sigma2 = xp.float64(step_dist**2) * 2.
+    assert bonded_force_dict["Backbone"]["Type"] == 'Harmonic'
 
+    rest_length = bonded_force_dict["Backbone"]["Rest length"]
+    wiggle_dist = bonded_force_dict["Backbone"]["Wiggle distance"]
+
+    rest_length = xp.float64(rest_length)
+    k_stretch = xp.float64(1./wiggle_dist**2)
+
     hbox = xp.asarray(system.state.box.L, dtype=xp.float64) / 2.
 
     old_bond_array = xp.asarray(current_bond_list, dtype=xp.int32)
@@ -76,17 +81,17 @@ def boltzmann_criterion(system, current_bond_list, trial_bond_list, step_dist=0.
     with system.state.gpu_local_snapshot as local_snap:
         N = int(new_bond_array.shape[0])
         rng = xp.random.random(N).astype(xp.float64)
-        
+
         rtags = local_snap.particles.rtag._coerce_to_ndarray()
         positions = local_snap.particles.position._coerce_to_ndarray()
         positions = positions.astype(xp.float64)
 
         num_blocks = (N+threads_per_block-1) // threads_per_block
 
-        _harmonic_boltzmann_filter(
+        _harmonic_distance_filter(
 			(num_blocks,),
 			(threads_per_block,),
-			(N, mu, sigma2,
+			(N, k_stretch, rest_length,
 			 rng, hbox, positions, rtags,
 			 old_bond_array, new_bond_array)
         )

examples/tutorial_stretch.ipynb

@@ -1,8 +1,8 @@
 extern "C" {
-__global__ void _harmonic_boltzmann_filter(
+__global__ void _harmonic_distance_filter(
         const unsigned int N,
-        const double mu,
-        const double sigma2,
+        const double k_stretch,
+        const double rest_length,
         const double* rng,
         const double* hdims,
         const double* positions,
@@ -21,45 +21,45 @@ __global__ void _harmonic_boltzmann_filter(
     if ( (old_bond[i].x == -1) || (old_bond[i].y == -1) )
         return;
 
-     if ( (new_bond[i].x == -1) || (new_bond[i].y == -1) )
-        return;   
+    if ( (new_bond[i].x == -1) || (new_bond[i].y == -1) )
+        return;
 
     double3* hdim = (double3*) hdims;
     double3* position = (double3*) positions;
 
     bool has_moved_left = (old_bond[i].x != new_bond[i].x);
     bool has_moved_right = (old_bond[i].y != new_bond[i].y);
+    
+    auto harmonic_work = [=](unsigned int& i, unsigned int& j) {
+    
+        double dx = position[j].x - position[i].x;
+        double dy = position[j].y - position[i].y;
+        double dz = position[j].z - position[i].z;
+		
+        dx -= (copysign(hdim->x, dx-hdim->x) + copysign(hdim->x, dx+hdim->x));
+        dy -= (copysign(hdim->y, dy-hdim->y) + copysign(hdim->y, dy+hdim->y));
+        dz -= (copysign(hdim->z, dz-hdim->z) + copysign(hdim->z, dz+hdim->z));
+        
+        double step_size = sqrt(dx*dx + dy*dy + dz*dz);
+        double harmonic_force = k_stretch * (step_size - rest_length);
+        
+        return harmonic_force * step_size;
+    };
 
-    if ( has_moved_left || has_moved_right ) {
-        double delta_r2 = 0.;
+    if ( has_moved_left ) {
+        unsigned int rtag1 = rtags[old_bond[i].x];
+        unsigned int rtag2 = rtags[new_bond[i].x];
 
-        for ( unsigned int j = 0; j < 2; ++j ) {
-            int2 tags = ((j==0) ? old_bond[i] : new_bond[i]);
+        if ( rng[i] > exp(-harmonic_work(rtag1, rtag2)) )
+            new_bond[i].x = old_bond[i].x;
+    }
 
-            unsigned int rtag1 = rtags[tags.x];
-            unsigned int rtag2 = rtags[tags.y];
-			
-            double dx = position[rtag2].x - position[rtag1].x;
-            double dy = position[rtag2].y - position[rtag1].y;
-            double dz = position[rtag2].z - position[rtag1].z;
-			
-            dx -= (copysign(hdim->x, dx-hdim->x) + copysign(hdim->x, dx+hdim->x));
-            dy -= (copysign(hdim->y, dy-hdim->y) + copysign(hdim->y, dy+hdim->y));
-            dz -= (copysign(hdim->z, dz-hdim->z) + copysign(hdim->z, dz+hdim->z));
-
-            double delta_r = sqrt(dx*dx + dy*dy + dz*dz) - mu;
-            delta_r2 += ((j==0) ? delta_r*delta_r : -delta_r*delta_r);
-        }
+    if ( has_moved_right ) {
+        unsigned int rtag1 = rtags[old_bond[i].y];
+        unsigned int rtag2 = rtags[new_bond[i].y];
 
-        double boltzmann_weight = exp(delta_r2/sigma2);
-
-        if ( rng[i] > boltzmann_weight ) {
-            if ( has_moved_left )
-                new_bond[i].x = old_bond[i].x;
-
-            if ( has_moved_right )
-                new_bond[i].y = old_bond[i].y;
-        }
+        if ( rng[i] > exp(-harmonic_work(rtag1, rtag2)) )
+            new_bond[i].y = old_bond[i].y;
     }
 }