Chain dynamics implementation

Strent (strain+entanglement) and chain dynamics are calculated using CUDA GPU kernels.

We have 2 types of kernels:

1. strent_kernel

Here strent-related quantities are calculted. We're using 2D (strents x chains) arrays and kernels to do so.

a. Apply strent creation or destruction calculated before and shift other strent respectively

float4 QN;
if (fetch_new_strent(i, oft))
    QN = d_new_strent[j]; //if new strent should go here
else
    QN = tex2D(t_a_QN, make_offset(i, oft), j);

float tcd;
    if (dCD_flag) { //If constraint dynamics is enabled
        if (fetch_new_strent(i, oft))
            tcd = d_new_tau_CD[j];
    else
    	tcd = tex2D(t_a_tCD, make_offset(i, oft), j);
} else
	tcd = 0;

surf2Dwrite(QN, s_b_QN, 16 * i, j);
surf2Dwrite(tcd, s_b_tCD, 4 * i, j);

b. Calculate probabilities for one Kuhn step to jump to or from next strent

float2 wsh = make_float2(0.0f, 0.0f); //Variable for shifting probability

float sig1 = __fdividef(0.75f, QN.w * (QN.w + 1));
float sig2 = __fdividef(0.75f, QN2.w * (QN2.w - 1));
float prefact1 = (Q == 0.0f) ? 1.0f : __fdividef(QN.w, (QN.w + 1));
float prefact2 = (Q2 == 0.0f) ? 1.0f : __fdividef(QN2.w, (QN2.w - 1));
float f1 = (Q == 0.0f) ? 2.0f * QN.w + 0.5f : QN.w;
float f2 = (Q2 == 0.0f) ? 2.0f * QN2.w - 0.5f : QN2.w;
float friction = __fdividef(2.0f, f1 + f2);
wsh.x = friction * __powf(prefact1 * prefact2, 0.75f)* __expf(Q * sig1 - Q2 * sig2); //P

c. Write sum of Wsh+ (SD jump to next strent), Wsh- (SD jump from next strent), Wcdd (CD destruction of entanglement) and Wcdc(CD creation of entanglement) for current strent

surf2Dwrite(wsh.x + wsh.y + dCD_flag * (tcd + d_CD_create_prefact * (QN.w - 1.0f)), s_sum_W, 4 * i, j);

2. chain_kernel

a. Sum up Wsh+/- and Wcdd/Wcdc probabilities for every strent (except of the last) previously calculated in strent_kernel

b. Calculate Wcdc (CD creation of entaglement) for the last strent

float W_CD_c_z = dCD_flag * d_CD_create_prefact * (QNtail.w - 1.0f);

c. Calculate Wsdc/Wsdd (SD creation/destruction probabilities) for first and last strents

d. Sum up probabilities of all possible processes in the chain

float sumW = sum_wshpm + W_SD_c_1 + W_SD_c_z + W_SD_d_1 + W_SD_d_z + W_CD_c_z;

e. Compare with random number assigned for current chain and pick a strent where something will happen

float pr = (sumW) * tex2D(t_uniformrand, rand_used[i], i);
rand_used[i]++;
int j = 0;
float tpr = 0.0f;
if (tz != 1)
    surf2Dread(&tpr, s_sum_W, 4 * j, i);
while ((pr >= tpr) && (j < tz - 2)) {
    pr -= tpr;
    j++;
    surf2Dread(&tpr, s_sum_W, 4 * j, i);
}

f1. Choose one of the following processes via comparing with the same random number successively

-CD destruction

-SD shift -

-SD shift +

-CD creation

f2. If we reached j=tz-2 in e. we skip f1 and compare the following probabilities:

-CD creation z (last strent)

-SD creation 1 (first strent)

-SD creation z (last strent)

-SD destruction 1 (first strent)

-SD destruction z (last strent)