diff --git a/src/core/polymer.cpp b/src/core/polymer.cpp index 11f2f442a64..5b7a4817d9c 100755 --- a/src/core/polymer.cpp +++ b/src/core/polymer.cpp @@ -448,454 +448,10 @@ int counterionsC(PartCfg & partCfg, int N_CI, int part_id, int mode, double shie return (std::max(max_cnt, cnt1)); } -int saltC(PartCfg & partCfg, int N_pS, int N_nS, int part_id, int mode, double shield, int max_try, - double val_pS, double val_nS, int type_pS, int type_nS, double rad) { - int n, cnt1, max_cnt; - double pos[3], dis2; - - cnt1 = max_cnt = 0; - - /* Place positive salt ions */ - for (n = 0; n < N_pS; n++) { - for (cnt1 = 0; cnt1 < max_try; cnt1++) { - if (rad > 0.) { - pos[0] = rad * (2. * d_random() - 1.); - pos[1] = rad * (2. * d_random() - 1.); - pos[2] = rad * (2. * d_random() - 1.); - dis2 = pos[0] * pos[0] + pos[1] * pos[1] + pos[2] * pos[2]; - pos[0] += box_l[0] * 0.5; - pos[1] += box_l[1] * 0.5; - pos[2] += box_l[2] * 0.5; - if (((mode != 0) || (collision(partCfg, pos, shield, 0, nullptr) == 0)) && - (dis2 < (rad * rad))) - break; - } else { - pos[0] = box_l[0] * d_random(); - pos[1] = box_l[1] * d_random(); - pos[2] = box_l[2] * d_random(); - if ((mode != 0) || (collision(partCfg, pos, shield, 0, nullptr) == 0)) - break; - } - POLY_TRACE(printf("p"); fflush(nullptr)); - } - if (cnt1 >= max_try) - return (-1); - if (place_particle(part_id, pos) == ES_PART_ERROR) - return (-3); - if (set_particle_q(part_id, val_pS) == ES_ERROR) - return (-3); - if (set_particle_type(part_id, type_pS) == ES_ERROR) - return (-3); - part_id++; - - max_cnt = std::max(cnt1, max_cnt); - POLY_TRACE(printf("P"); fflush(nullptr)); - } - POLY_TRACE(printf(" %d->%d \n", cnt1, max_cnt)); - if (cnt1 >= max_try) - return (-1); - - /* Place negative salt ions */ - for (n = 0; n < N_nS; n++) { - for (cnt1 = 0; cnt1 < max_try; cnt1++) { - if (rad > 0.) { - pos[0] = rad * (2. * d_random() - 1.); - pos[1] = rad * (2. * d_random() - 1.); - pos[2] = rad * (2. * d_random() - 1.); - dis2 = pos[0] * pos[0] + pos[1] * pos[1] + pos[2] * pos[2]; - pos[0] += box_l[0] * 0.5; - pos[1] += box_l[1] * 0.5; - pos[2] += box_l[2] * 0.5; - if (((mode != 0) || (collision(partCfg, pos, shield, 0, nullptr) == 0)) && - (dis2 < (rad * rad))) - break; - } else { - pos[0] = box_l[0] * d_random(); - pos[1] = box_l[1] * d_random(); - pos[2] = box_l[2] * d_random(); - if ((mode != 0) || (collision(partCfg, pos, shield, 0, nullptr) == 0)) - break; - } - POLY_TRACE(printf("n"); fflush(nullptr)); - } - if (cnt1 >= max_try) - return (-1); - if (place_particle(part_id, pos) == ES_PART_ERROR) - return (-3); - if (set_particle_q(part_id, val_nS) == ES_ERROR) - return (-3); - if (set_particle_type(part_id, type_nS) == ES_ERROR) - return (-3); - part_id++; - max_cnt = std::max(cnt1, max_cnt); - POLY_TRACE(printf("N"); fflush(nullptr)); - } - POLY_TRACE(printf(" %d->%d \n", cnt1, max_cnt)); - if (cnt1 >= max_try) - return (-2); - return (std::max(max_cnt, cnt1)); -} - -double velocitiesC(double v_max, int part_id, int N_T) { - double v[3], v_av[3]; - int i; - v_av[0] = v_av[1] = v_av[2] = 0.0; - for (i = part_id; i < part_id + N_T; i++) { - do { - v[0] = v_max * 2. * (d_random() - .5) * time_step; - v[1] = v_max * 2. * (d_random() - .5) * time_step; - v[2] = v_max * 2. * (d_random() - .5) * time_step; - // note that time_step == -1, as long as it is not yet set - } while (sqrt(SQR(v[0]) + SQR(v[1]) + SQR(v[2])) > v_max * fabs(time_step)); - v_av[0] += v[0]; - v_av[1] += v[1]; - v_av[2] += v[2]; - if (set_particle_v(i, v) == ES_ERROR) { - fprintf(stderr, "INTERNAL ERROR: failed upon setting one of the " - "velocities in Espresso (current average: %f)!\n", - sqrt(SQR(v_av[0]) + SQR(v_av[1]) + SQR(v_av[2]))); - fprintf(stderr, "Aborting...\n"); - errexit(); - } - } - // note that time_step == -1, as long as it is not yet set - return (sqrt(SQR(v_av[0]) + SQR(v_av[1]) + SQR(v_av[2])) / fabs(time_step)); -} - -double maxwell_velocitiesC(int part_id, int N_T) { - double v[3], v_av[3], uniran[2]; - int i; - int flag = 1; - uniran[0] = d_random(); - uniran[1] = d_random(); - v_av[0] = v_av[1] = v_av[2] = 0.0; - for (i = part_id; i < part_id + N_T; i++) { - if (flag == 1) { - v[0] = pow((-2. * log(uniran[0])), 0.5) * cos(2. * PI * uniran[1]) * - time_step; - v[1] = pow((-2. * log(uniran[1])), 0.5) * sin(2. * PI * uniran[0]) * - time_step; - uniran[0] = d_random(); - uniran[1] = d_random(); - v[2] = pow((-2. * log(uniran[0])), 0.5) * cos(2. * PI * uniran[1]) * - time_step; - flag = 0; - } else { - v[0] = pow((-2. * log(uniran[1])), 0.5) * sin(2. * PI * uniran[0]) * - time_step; - uniran[0] = d_random(); - uniran[1] = d_random(); - v[1] = pow((-2. * log(uniran[0])), 0.5) * cos(2. * PI * uniran[1]) * - time_step; - v[2] = pow((-2. * log(uniran[1])), 0.5) * sin(2. * PI * uniran[0]) * - time_step; - flag = 1; - } - // printf("%f \n %f \n %f \n",v[0],v[1],v[2]); - v_av[0] += v[0]; - v_av[1] += v[1]; - v_av[2] += v[2]; - if (set_particle_v(i, v) == ES_ERROR) { - fprintf(stderr, "INTERNAL ERROR: failed upon setting one of the " - "velocities in Espresso (current average: %f)!\n", - sqrt(SQR(v_av[0]) + SQR(v_av[1]) + SQR(v_av[2]))); - fprintf(stderr, "Aborting...\n"); - errexit(); - } - } - // note that time_step == -1, as long as it is not yet set - return (sqrt(SQR(v_av[0]) + SQR(v_av[1]) + SQR(v_av[2])) / fabs(time_step)); -} - -int collectBonds(PartCfg & partCfg, int mode, int part_id, int N_P, int MPC, int type_bond, - std::vector bond, std::vector> bonds) { - int i, j, k, ii, size; - - partCfg.update_bonds(); - - if (mode == 1) { - /* Find all the bonds leading to and from the ending monomers of the chains. - */ - bond.resize(2 * N_P); - bonds.resize(2 * N_P); - for (i = 0; i < 2 * N_P; i++) { - bond[i] = 0; - bonds[i].resize(1); - } - for (k = part_id; k < N_P * MPC + part_id; k++) { - i = 0; - while (i < partCfg[k].bl.n) { - size = bonded_ia_params[partCfg[k].bl.e[i]].num; - if (partCfg[k].bl.e[i++] == type_bond) { - for (j = 0; j < size; j++) { - if ((partCfg[k].p.identity % MPC == 0) || - ((partCfg[k].p.identity + 1) % MPC == 0)) { - ii = partCfg[k].p.identity % MPC - ? 2 * (partCfg[k].p.identity + 1) / MPC - 1 - : 2 * partCfg[k].p.identity / MPC; - bonds[i].resize(bond[i] + 1) ; - bonds[ii][bond[ii]++] = partCfg[k].bl.e[i]; - } else if ((partCfg[k].bl.e[i] % MPC == 0) || - ((partCfg[k].bl.e[i] + 1) % MPC == 0)) { - ii = partCfg[k].bl.e[i] % MPC ? 2 * (partCfg[k].bl.e[i] + 1) / MPC - 1 - : 2 * partCfg[k].bl.e[i] / MPC; - bonds[i].resize(bond[i] + 1); - bonds[ii][bond[ii]++] = partCfg[k].p.identity; - } - i++; - } - } else - i += size; - } - } - POLY_TRACE(for (i = 0; i < 2 * N_P; i++) { - printf("(%d) %d:\t", i, i % 2 ? (i + 1) * MPC / 2 - 1 : i * MPC / 2); - if (bond[i] > 0) - for (j = 0; j < bond[i]; j++) - printf("%d ", bonds[i][j]); - printf("\t=%d\n", bond[i]); - }); - } else if (mode == 2) { - /* Find all the bonds leading to and from each monomer. */ - bond.resize(N_P * MPC); - bonds.resize(N_P * MPC); - for (i = 0; i < N_P * MPC + part_id; i++) { - bond[i] = 0; - bonds[i].resize(1); - } - for (k = part_id; k < N_P * MPC + part_id; k++) { - i = 0; - while (i < partCfg[k].bl.n) { - size = bonded_ia_params[partCfg[k].bl.e[i]].num; - if (partCfg[k].bl.e[i++] == type_bond) { - for (j = 0; j < size; j++) { - ii = partCfg[k].bl.e[i]; - bonds[k].resize(bond[k] + 1); - bonds[k][bond[k]++] = ii; - bonds[ii].resize(bond[ii] + 1); - bonds[ii][bond[ii]++] = k; - i++; - } - } else - i += size; - } - } - POLY_TRACE(for (i = 0; i < N_P * MPC + part_id; i++) { - printf("%d:\t", i); - if (bond[i] > 0) - for (j = 0; j < bond[i]; j++) - printf("%d ", bonds[i][j]); - printf("\t=%d\n", bond[i]); - }); - } else { - fprintf(stderr, "Unknown mode %d requested!\nAborting...\n", mode); - fflush(nullptr); - return (-2); - } - - return (0); -} - -int crosslinkC(PartCfg & partCfg, int N_P, int MPC, int part_id, double r_catch, int link_dist, - int chain_dist, int type_bond, int max_try) { - int i, j, k, ii, size, bondN[2], crossL; - std::vector bond; - std::vector> bonds; - - /* Find all the bonds leading to and from each monomer. */ - if (collectBonds(partCfg, 2, part_id, N_P, MPC, type_bond, bond, bonds)) - return (-2); - POLY_TRACE(for (i = 0; i < N_P * MPC + part_id; i++) { - printf("%d:\t", i); - if (bond[i] > 0) - for (j = 0; j < bond[i]; j++) - printf("%d ", bonds[i][j]); - printf("\t=%d\n", bond[i]); - }); - - /* Find all possible binding partners in the neighbourhood of the unconnected - * ending monomers. */ - std::vector link(2 * N_P); - std::vector> links(2 * N_P); - for (i = 0; i < N_P; i++) { - for (k = 0; k < 2; k++) { - if (bond[i * MPC + k * (MPC - 1)] == 1) { - links[2 * i + k].resize(n_part); - link[2 * i + k] = mindist3(partCfg,i * MPC + k * (MPC - 1) + part_id, r_catch, - links[2 * i + k].data()); - links[2 * i + k].resize(link[2 * i + k]); - } else if (bond[i * MPC + k * (MPC - 1)] == 2) - link[2 * i + k] = -1; /* Note that links[2*i+k] will not be malloc()ed - now (taken care of at end)!!! */ - else { - fprintf( - stderr, - "Runaway end-monomer %d detected (has %d bonds)!\nAborting...\n", - i * N_P + k * (MPC - 1) + part_id, bond[i * MPC + k * (MPC - 1)]); - fflush(nullptr); - return (-2); - } - POLY_TRACE(printf("%d: ", i * MPC + k * (MPC - 1) + part_id); - for (j = 0; j < link[2 * i + k]; j++) - printf("%d ", links[2 * i + k][j]); - printf("\t=%d\n", link[2 * i + k]); fflush(nullptr)); - } - } - - /* Throw out all the monomers which are ends, which are too close to the - * ending monomers on the same chain, or which are no monomers at all. */ - for (i = 0; i < N_P; i++) { - for (k = 0; k < 2; k++) { - size = 0; - ii = i * MPC + k * (MPC - 1) + part_id; - if (link[2 * i + k] >= 0) { - for (j = 0; j < link[2 * i + k]; j++) { /* only monomers && ((same - chain, but sufficiently far - away) || (different chain)) - */ - if ((links[2 * i + k][j] < N_P * MPC + part_id) && - (((abs(links[2 * i + k][j] - ii) > chain_dist) || - (abs(links[2 * i + k][j] - i * MPC) > (1. * MPC))))) - if ((links[2 * i + k][j] % MPC != 0) && - ((links[2 * i + k][j] + 1) % MPC != 0)) - links[2 * i + k][size++] = - links[2 * i + k][j]; /* no ends accepted */ - } - link[2 * i + k] = size; - links[2 * i + k].resize(link[2 * i + k]); - } - POLY_TRACE(printf("%d: ", ii); for (j = 0; j < link[2 * i + k]; j++) - printf("%d ", links[2 * i + k][j]); - printf("\t=%d\n", link[2 * i + k]); fflush(nullptr)); - } - } - - /* Randomly choose a partner (if not available -> '-1') for each polymer - * chain's end if it's not already been crosslinked (-> '-2'). */ - std::vector cross(2 * N_P); - crossL = 0; - for (i = 0; i < 2 * N_P; i++) - if (link[i] > 0) { - cross[i] = links[i][(int)dround(d_random() * (link[i] - 1))]; - crossL++; - } else { - cross[i] = -1 + link[i]; - crossL -= link[i]; - } - POLY_TRACE(for (i = 0; i < 2 * N_P; i++) - printf("%d -> %d \t", - i % 2 ? (i + 1) * MPC / 2 - 1 : i * MPC / 2, cross[i]); - printf("=> %d\n", crossL); fflush(nullptr)); - - /* Remove partners (-> '-3') if they are less than link_dist apart and retry. - */ - k = 0; - ii = 1; - while ((k < max_try) && (ii > 0)) { - POLY_TRACE(printf("Check #%d: ", k)); - for (i = 0; i < 2 * N_P; i++) { - if (cross[i] >= 0) { - for (j = 0; j < 2 * N_P; - j++) { /* In the neighbourhood of each partner shall be no future - crosslinks (preventing stiffness). */ - if ((j != i) && (cross[j] >= 0) && - (abs(cross[j] - cross[i]) < link_dist)) { - cross[i] = -3; - cross[j] = -3; - crossL -= 2; - POLY_TRACE(printf("%d->%d! ", i, j)); - break; - } - } - if (cross[i] == -3) - continue; /* Partners shall not be too close to the chain's ends - (because these will be crosslinked at some point). */ - if ((cross[i] % MPC < link_dist) || - (cross[i] % MPC >= MPC - link_dist)) { - cross[i] = -3; - crossL--; - POLY_TRACE(printf("%d->end! ", i)); - } else { /* In the neighbourhood of each partner there shall be no other - crosslinks (preventing stiffness). */ - for (j = cross[i] - link_dist + 1; j < cross[i] + link_dist - 1; - j++) { - if ((j % MPC == 0) || ((j + 1) % MPC == 0)) - size = 1; - else - size = 2; - if ((bond[j] > size) && (j - floor(i / 2.) * MPC < MPC)) { - cross[i] = -3; - crossL--; - POLY_TRACE(printf("%d->link! ", i)); - break; - } - } - } - } - } - POLY_TRACE(printf("complete => %d CL left; ", crossL)); - if (k == max_try - 1) - break; - else - ii = 0; /* Get out if max_try is about to be reached, preventing dangling - unchecked bond suggestions. */ - if (crossL < 2 * N_P) { - for (i = 0; i < 2 * N_P; - i++) { /* If crosslinks violated the rules & had to be removed, - create new ones now. */ - if (cross[i] == -3) { - ii++; - if (link[i] > 0) { - cross[i] = links[i][(int)dround(d_random() * (link[i] - 1))]; - crossL++; - } else { - return (-2); - } - } - } - } - POLY_TRACE(printf("+ %d new = %d CL.\n", ii, crossL)); - if (ii > 0) - k++; - } - POLY_TRACE(for (i = 0; i < 2 * N_P; i++) - printf("%d -> %d \t", - i % 2 ? (i + 1) * MPC / 2 - 1 : i * MPC / 2, cross[i]); - printf("=> %d\n", crossL); fflush(nullptr)); - - /* Submit all lawful partners as new bonds to Espresso (observing that bonds - * are stored with the higher-ID particle only). */ - if (k >= max_try) - return (-1); - - { - size = 0; - for (i = 0; i < N_P; i++) { - if (cross[2 * i] >= 0) { - bondN[0] = type_bond; - bondN[1] = i * MPC + part_id; - size++; - if (change_particle_bond(cross[2 * i], bondN, 0) == ES_ERROR) - return (-3); - } - if (cross[2 * i + 1] >= 0) { - bondN[0] = type_bond; - bondN[1] = cross[2 * i + 1]; - size++; - if (change_particle_bond(i * MPC + (MPC - 1) + part_id, bondN, 0) == - ES_ERROR) - return (-3); - } - } - POLY_TRACE(printf("Created %d new bonds; now %d ends are crosslinked!\n", - size, crossL)); - return (crossL); - } -} int diamondC(PartCfg & partCfg, double a, double bond_length, int MPC, int N_CI, double val_nodes, double val_cM, double val_CI, int cM_dist, int nonet) { diff --git a/src/core/polymer.hpp b/src/core/polymer.hpp index 21e122ebfbe..d90d296b93d 100755 --- a/src/core/polymer.hpp +++ b/src/core/polymer.hpp @@ -122,72 +122,6 @@ int polymerC(PartCfg &, int N_P, int MPC, double bond_length, int part_id, doubl int counterionsC(PartCfg &, int N_CI, int part_id, int mode, double shield, int max_try, double val_CI, int type_CI); -/** C implementation of 'salt \ \ [options]', - @param N_pS = number of positive salt ions to create - @param N_nS = number of negative salt ions to create - @param part_id = particle number of the first salt ion (defaults to - 'n_total_particles') - @param mode = selects setup mode: Self avoiding walk (SAW) or plain - random walk (RW) (defaults to 'SAW') - @param shield = shield around each particle another particle's - position may not enter if using SAW (defaults to '0') - @param max_try = how often a monomer should be reset if current - position collides with a previous particle (defaults to '30000') - @param val_pS = valency of the positive salt ions (defaults to '1') - @param val_nS = valencies of the negative salt ions (default to '-1') - @param type_pS = type numbers to be used for positive salt ions - (defaults to '3') - @param type_nS = type numbers to be used for negative salt ions - (defaults to '4') - @param rad = radius of the sphere for the cell model. If value is - >0 the ions are distributed in a sphere centered - in the simulation box. - @return Returns how often the attempt to place a particle failed in the - worst case. */ -int saltC(PartCfg &, int N_pS, int N_nS, int part_id, int mode, double shield, int max_try, - double val_pS, double val_nS, int type_pS, int type_nS, double rad); - -/** C implementation of 'velocities \ [options]', - @param v_max = maximum velocity to be used - @param part_id = particle number of the first of the \ particles - (defaults to '0') - @param N_T = number of particles of which the velocities should be - set (defaults to 'n_total_particles - part_id') - @return Returns the averaged velocity when done. */ -double velocitiesC(double v_max, int part_id, int N_T); - -/** C implementation of 'maxwell_velocities [options]', - @param part_id = particle number of the first of the \ particles - (defaults to '0') - @param N_T = number of particles of which the velocities should be - set (defaults to 'n_total_particles - part_id') - @return Returns the averaged velocity when done. */ -double maxwell_velocitiesC(int part_id, int N_T); - -/** Collects the bonds leading to and from the ending monomers of the chains - (mode == 1) or - all the bonds leading to and from each monomer (mode == 2). - @return Returns '0' upon success, '-2' otherwise. */ -int collectBonds(PartCfg &, int mode, int part_id, int N_P, int MPC, int type_bond, - std::vector bond, std::vector> bonds); - -/** C implementation of 'crosslink \ \ [options]', - @param N_P = number of polymer chains - @param MPC = monomers per chain - @param part_id = particle number of the start monomer (defaults to '0') - @param r_catch = maximum length of a crosslink (defaults to '1.9') - @param link_dist = minimum distance between the indices of two - crosslinked monomers with different binding partners (defaults to '2') - @param chain_dist = same as \, but for monomers of the same - bond (defaults to \ =\> no bonds to the same chain allowed) - @param type_FENE = type number of the FENE-typed bonded interaction bonds - to be set between the monomers (defaults to '0') - @param max_try = how often crosslinks should be removed if they are too - close to other links (defaults to '30000') - @return Returns how many ends are now successfully linked. */ -int crosslinkC(PartCfg &, int N_P, int MPC, int part_id, double r_catch, int link_dist, - int chain_dist, int type_FENE, int max_try); - /** C implementation of 'diamond \ \ \ [options]' */ int diamondC(PartCfg &, double a, double bond_length, int MPC, int N_CI, double val_nodes, double val_cM, double val_CI, int cM_dist, int nonet);