/* ** This program calculates j-factor for a chain with arbitrary shape and ** flexibility using the algorithm described in ** ** Podtelezhnikov, A. A.; C. Mao; N. C. Seeman and A. V. Vologodskii. (2000) ** Multimerization-cyclization of DNA fragments as a method of conformational ** analysis. Biophys. J. 79(5) : 2692-2704. ** ** jfm2full.c, version 1.5 ** ** Copyright (c) 2000 - 2005 Alexei Podtelezhnikov */ #define VER "jfm2full 1.5, Copyright (c) 2000 - 2005 Alexei Podtelezhnikov\n" #include #include #include #include #define NOTHING 1.0e+37 /* a huge number */ #ifndef M_PI # define M_PI 3.14159265358979323846 /* pi */ #endif #ifndef M_SQRT2 # define M_SQRT2 1.41421356237309504880 /* sqrt(2) */ #endif /* intentionally dangerous internal macro */ #ifndef WITH_SINCOS # define sin_cos(si, co, x) si = sin(x); co = cos(x) #else # define sin_cos(si, co, x) asm ("fsincos" : "=t" (co), "=u" (si) : "0" (x)) #endif #define within(i, j, k) ((i <= j && j < k) || (k < i && (i <= j || j < k))) typedef double vector[3]; typedef vector triplet[3]; typedef double matrix[3][3]; double *add(vector, vector, vector); double *subtract(vector, vector, vector); double dotprod(vector, vector); double *crossprod(vector, vector, vector); double distance(vector, vector); void normalize(vector); void normalize_1(vector); /* other prototypes */ double totenergy(void); double tottwist(void); double twist(void); void chain4view(void); double writhe(void); double writhe1(void); double writhe2(void); /* storage locations */ struct Data { triplet xi; /* ground state */ double kac, kb; /* half-rigidities */ } *set; triplet *xx; vector *rr; double *erg; /* all defaults */ int NS = 0, Lk = 0; int FREQ = 100000; double bps = 3.4829; /* gives density */ double diam2 = 0.0; double rdata[3] = { 0.1, 0.1, 0.3 }; const double relaxer[3] = { 4.0, 4.0, 2.0 }; double bound[3]; double condit[3]; double current[3]; /* ********** Trivial vector, triplet, and matrix functions */ /* vector addition */ double *add(vector apb, vector a, vector b) { apb[0] = a[0] + b[0]; apb[1] = a[1] + b[1]; apb[2] = a[2] + b[2]; return apb; } /* vector subtraction */ double *subtract(vector amb, vector a, vector b) { amb[0] = a[0] - b[0]; amb[1] = a[1] - b[1]; amb[2] = a[2] - b[2]; return amb; } double dotprod(vector a, vector b) { return a[0] * b[0] + a[1] * b[1] + a[2] * b[2]; } double *crossprod(vector ab, vector a, vector b) { ab[0] = a[1] * b[2] - a[2] * b[1]; ab[1] = a[2] * b[0] - a[0] * b[2]; ab[2] = a[0] * b[1] - a[1] * b[0]; return ab; } void normalize(vector a) { double l; l = 1.0 / sqrt(dotprod(a, a)); a[0] *= l; a[1] *= l; a[2] *= l; } /* approximate normalization of a vector with the length close to 1.0 */ void normalize_1(vector a) { double l; l = 1.5 - 0.5 * dotprod(a, a); a[0] *= l; a[1] *= l; a[2] *= l; } void fixtriplet(triplet x) { double c1; normalize_1(x[2]); c1 = dotprod(x[1], x[2]); x[1][0] -= c1 * x[2][0]; x[1][1] -= c1 * x[2][1]; x[1][2] -= c1 * x[2][2]; normalize_1(x[1]); crossprod(x[0], x[1], x[2]); } void printout(triplet x) { printf("(%g %g %g) %g\n", x[0][0], x[0][1], x[0][2], dotprod(x[0], x[0])); printf("(%g %g %g) %g\n", x[1][0], x[1][1], x[1][2], dotprod(x[1], x[1])); printf("(%g %g %g) %g\n", x[2][0], x[2][1], x[2][2], dotprod(x[2], x[2])); printf(" %g %g %g\n", dotprod(x[0], x[1]), dotprod(x[1], x[2]), dotprod(x[2], x[0])); } void casttriplet(triplet x, triplet y) { x[0][0] = y[0][0]; x[0][1] = y[0][1]; x[0][2] = y[0][2]; x[1][0] = y[1][0]; x[1][1] = y[1][1]; x[1][2] = y[1][2]; x[2][0] = y[2][0]; x[2][1] = y[2][1]; x[2][2] = y[2][2]; } /* t=trans(x) */ void transset(matrix t, triplet x) { t[0][0] = x[0][0]; t[0][1] = x[1][0]; t[0][2] = x[2][0]; t[1][0] = x[0][1]; t[1][1] = x[1][1]; t[1][2] = x[2][1]; t[2][0] = x[0][2]; t[2][1] = x[1][2]; t[2][2] = x[2][2]; } void matrixvector(vector ta, matrix t, vector a) { ta[0] = dotprod(t[0], a); ta[1] = dotprod(t[1], a); ta[2] = dotprod(t[2], a); } void rotation(triplet tx, matrix t, triplet x) { matrixvector(tx[2], t, x[2]); matrixvector(tx[1], t, x[1]); crossprod(tx[0], tx[1], tx[2]); } void randvector(vector z) { double sb, cb, gamma, sc, cc; const double discrete = 2.0 / RAND_MAX; cb = discrete * rand() - 1.0; sb = sqrt(1.0 - cb * cb); gamma = M_PI * discrete * rand(); sin_cos(sc, cc, gamma); z[0] = sb * cc; z[1] = sb * sc; z[2] = cb; } /* set rotation matrix given the direction and the angle using the well-known quaternion-to-matrix transformation */ void rotmatrix(matrix t, vector a, double alpha) { double si, co, q0, q1, q2, q3, p; alpha *= 0.5; sin_cos(si, co, alpha); si *= M_SQRT2; co *= M_SQRT2; q0 = co; q1 = a[0] * si; q2 = a[1] * si; q3 = a[2] * si; p = q0 * q0 - 1.0; t[0][0] = p + q1 * q1; t[1][1] = p + q2 * q2; t[2][2] = p + q3 * q3; t[0][1] = t[1][0] = q1 * q2; t[1][2] = t[2][1] = q2 * q3; t[2][0] = t[0][2] = q3 * q1; p = q0 * q1; t[2][1] += p; t[1][2] -= p; p = q0 * q2; t[0][2] += p; t[2][0] -= p; p = q0 * q3; t[1][0] += p; t[0][1] -= p; } /* set triplet given the three Euler angles (floating x-convention) */ void eulerset(triplet x, double alpha, double beta, double gamma) { double sa, ca, sb, cb, sc, cc; sin_cos(sa, ca, alpha); sin_cos(sb, cb, beta); sin_cos(sc, cc, gamma); x[0][0] = ca * cc - sa * cb * sc; x[0][1] = sa * cc + ca * cb * sc; x[0][2] = sb * sc; x[1][0] = -ca * sc - sa * cb * cc; x[1][1] = -sa * sc + ca * cb * cc; x[1][2] = sb * cc; x[2][0] = sa * sb; x[2][1] = -ca * sb; x[2][2] = cb; } /* the Euler beta angle (fail-safe, no gimbal lock) */ double euler_bend(triplet x, triplet x0) { double cosb, sinbx, sinby; sinbx = dotprod(x[2], x0[0]); sinby = dotprod(x[2], x0[1]); cosb = dotprod(x[2], x0[2]); return atan2(sqrt(sinbx * sinbx + sinby * sinby), cosb); } /* the Euler alpha+gamma angle (fail-safe, no gimbal lock) */ double euler_twist(triplet x, triplet x0) { double pcosac, psinac; pcosac = dotprod(x[0], x0[0]) + dotprod(x[1], x0[1]); psinac = dotprod(x[0], x0[1]) - dotprod(x[1], x0[0]); return atan2(psinac, pcosac); } /* triplet orientation comparison in terms of convenient metrics */ void divergence(double *hsqbend, double *hsqtwist, triplet x, triplet y) { double cosb, pcosac; cosb = dotprod(x[2], y[2]); *hsqbend = 1.0 - cosb; pcosac = dotprod(x[0], y[0]) + dotprod(x[1], y[1]); *hsqtwist = 1.0 - pcosac / (1.0 + cosb); } /* squared distance between two points */ double distance(vector a, vector b) { vector h; subtract(h, a, b); return dotprod(h, h); } /* cylindrical clash check returns zero when the clash is found or suggests the distance to the next possible clash */ int clash(vector a, vector b, vector u, vector v, double d2) { vector ba; double a2, ddd, uv, b2, rna, rnb, ak, bk; subtract(ba, b, a); /* point-to-point */ ddd = a2 = dotprod(ba, ba); /* think of it as a chain of cylinders with beads in the joints */ /* let each segment have one bead and one cylinder */ /* checking collisions of beads on the other ends is just redundant */ /* quickly abandon segments that are too far apart to collide */ if (a2 > 4.0 + d2) return 1 + (int) (sqrt(a2) - sqrt(4.0 + d2)); if (ddd < d2) return 0; rna = dotprod(ba, u); if (0.0 < rna && rna < 1.0) { /* point-to-line */ ddd = a2 - rna * rna; if (ddd < d2) return 0; } rnb = dotprod(ba, v); if (-1.0 < rnb && rnb < 0.0) { /* point-to-line */ ddd = a2 - rnb * rnb; if (ddd < d2) return 0; } uv = dotprod(u, v); b2 = 1.0 - uv * uv; if (fabs(b2) > 1.0e-9) { ak = (rna - rnb * uv); bk = (-rnb + rna * uv); if (ak > 0.0 && ak < b2 && bk > 0.0 && bk < b2) { /* line-to-line */ ddd = a2 - (rna * ak - rnb * bk) / b2; if (ddd < d2) return 0; } } return 1; } /* ********** Metropolis Monte Carlo procedure */ double energy(triplet x, triplet x0, struct Data *p) { triplet xp; double ac, b; /* project */ rotation(xp, x0, x); b = euler_bend(xp, p->xi); ac = euler_twist(xp, p->xi); return p->kb * b * b + p->kac * ac * ac; } /* check for clashes between start-to-end fragment and the rest of the chain acceptable diameter-to-length ratio should be sufficiently less than two */ int excluded(int start, int end) { int i, j, jmp; /* virtual segment NS is not considered */ if (start == NS) start = 0; if (end == NS) end = 0; /* moving just a few segments would make this nested loop faster */ for (i = start; i != end;) { for (j = end; within(end, j, start);) { if (abs(i - j) > 1) { jmp = clash(rr[i], rr[j], xx[i][2], xx[j][2], diam2); if (jmp == 0) { /* fprintf(stderr, "%d %d\n", i, j); */ return 1; } } else jmp = 1; j += jmp; if (j >= NS) j = 0; /* wrap around */ } i++; if (i >= NS) i = 0; /* wrap around */ } return 0; } /* rotate a chain of segments wrapping around the end < start cases */ void rotcontour(int start, int end, matrix t) { int i; triplet x; /* update orientations */ for (i = start; i != end;) { rotation(x, t, xx[i]); casttriplet(xx[i], x); i++; if (i > NS) i = 0; /* wrap around */ } /* update positions */ if (start != 0) /* i.e., the start reference didn't move */ for (i = start; i != end && i < NS; i++) add(rr[i + 1], rr[i], xx[i][2]); if (start == 0 || end < start) /* also avoiding redundancy */ for (i = end; i != start && i > 0; i--) subtract(rr[i - 1], rr[i], xx[i - 1][2]); } int allowed(int k1, int k2, matrix rot) { double erg1, erg2, loss, hsqtwist, hsqbend, sqdist; triplet x; matrix tor; vector e1, e2; if (k1 == k2) /* ouch... */ return 0; sqdist = current[2]; if (k1 == 0) { subtract(e1, rr[0], rr[k2]); matrixvector(e2, rot, e1); subtract(e1, rr[NS], rr[k2]); sqdist = distance(e1, e2); } if (k2 == 0) { subtract(e1, rr[NS], rr[k1]); matrixvector(e2, rot, e1); subtract(e1, rr[0], rr[k1]); sqdist = distance(e1, e2); } if (sqdist > condit[2]) return 0; hsqbend = current[1]; hsqtwist = current[0]; if (condit[1] != NOTHING || condit[0] != NOTHING) { if (k1 == 0) { rotation(x, rot, xx[0]); divergence(&hsqbend, &hsqtwist, x, xx[NS]); } if (k2 == 0) { rotation(x, rot, xx[NS]); divergence(&hsqbend, &hsqtwist, xx[0], x); } if (hsqbend > condit[1] || hsqtwist > condit[0]) return 0; } erg1 = erg[k1]; erg2 = erg[k2]; if (k1 != 0) { rotation(x, rot, xx[k1]); erg1 = energy(x, xx[k1 - 1], set + k1); } if (k2 != 0) { rotation(x, rot, xx[k2 - 1]); erg2 = energy(xx[k2], x, set + k2); } /* Metropolis criteria */ loss = erg[k1] + erg[k2] - erg1 - erg2; if (loss < 0.0) if (exp(loss) * RAND_MAX < rand()) return 0; rotcontour(k1, k2, rot); /* update conformation */ /* lastly, very expensive self-avoidance check */ if (diam2 != 0.0 && excluded(k1, k2)) { transset(tor, rot); rotcontour(k1, k2, tor); /* undo move */ return 0; } current[0] = hsqtwist; current[1] = hsqbend; current[2] = sqdist; erg[k1] = erg1; erg[k2] = erg2; return 1; } void move() { int k1, k2, *cnt; double *ampl; vector a; matrix rot; static int cnt1 = 0, cnt2 = 0; static double ampl1 = 0.05, ampl2 = 0.2; const double discrete = 2.0 / RAND_MAX; k1 = rand(); if (k1 > RAND_MAX / 2) { /* short kink */ k1 = k1 % NS + 1; /* 1 <= k1 <= NS */ if (k1 > NS / 2) k2 = 0; else { k2 = k1; k1 = 0; } randvector(a); ampl = &l1; cnt = &cnt1; } else { /* crankshaft rotation */ k2 = rand() % (NS / 2) + 1; /* length, 1 <= k2 <= NS/2 */ k1 = k1 % (NS + 1); /* 0 <= k1 <= NS */ k2 = k1 + k2; if (k2 > NS) k2 -= NS + 1; /* wrap around */ subtract(a, rr[k2], rr[k1]); normalize(a); ampl = &l2; cnt = &cnt2; } rotmatrix(rot, a, (*ampl) * (discrete * rand() - 1.0)); if (allowed(k1, k2, rot)) { if (++*cnt > 32) { *cnt = 0; *ampl *= 1.1; /* printf("%.2g ", *ampl); */ } } else { if (--*cnt < -32) { *cnt = 0; *ampl *= 0.9; /* printf("%.2g ", *ampl); */ } } } /* ********** Simulation setup and j-factor calculations */ double jjfactor() { double v; char sym; static int count = 0, dt = 0, dtall = 0; static double prob = 1.; dtall++; if (current[count] < bound[count]) dt++; if (dtall != FREQ) return NOTHING; prob *= (double) dt / dtall; /* printf("%.2g ", (double) dt / dtall); */ /* initially condit[0] == bound[0] to relax the conformation */ if (dt == dtall && condit[count] != bound[count]) { condit[count] = NOTHING; /* remove it */ count++; } if (count < 3) { bound[count] = condit[count]; condit[count] *= relaxer[count]; /* relax it */ v = NOTHING; sym = '.'; } else { /* the calculation is over */ v = prob; sym = '\n'; prob = 1.; count = 0; } dtall = dt = 0; putchar(sym); fflush(stdout); return v; } void initiate() { int i; double angle; triplet x; matrix t; /* initialize conformation to a planar closed "circle" */ xx[0][0][0] = 1.; xx[0][0][1] = 0.; xx[0][0][2] = 0.; xx[0][1][0] = 0.; xx[0][1][1] = 1.; xx[0][1][2] = 0.; xx[0][2][0] = 0.; xx[0][2][1] = 0.; xx[0][2][2] = 1.; rr[0][0] = rr[0][1] = rr[0][2] = 0.; angle = 2 * M_PI / NS; rotmatrix(t, xx[0][0], angle); /* x-convention */ for (i = 0; i < NS; i++) { rotation(xx[i + 1], t, xx[i]); add(rr[i + 1], rr[i], xx[i][2]); } if (Lk != 0) { angle = 2 * M_PI * Lk / NS; for (i = 1; i <= NS; i++) { rotmatrix(t, xx[i][2], angle * i); rotation(x, t, xx[i]); casttriplet(xx[i], x); } } /* set confinements, current end separation and energies */ bound[0] = condit[0] = rdata[0]; bound[1] = condit[1] = rdata[1]; bound[2] = condit[2] = rdata[2]; divergence(current + 1, current, xx[0], xx[NS]); current[2] = distance(rr[0], rr[NS]); erg[0] = 0.; for (i = 1; i <= NS; i++) erg[i] = energy(xx[i], xx[i - 1], set + i); } double simulation() { int i, step = 0; double v = NOTHING; initiate(); /* the main simulation */ while (v == NOTHING) { move(); v = jjfactor(); /* printf("%g\n", totenergy()); */ step++; /* fix accumulating errors once in a while */ if (step % 1024 == 0) for (i = 0; i <= NS; i++) fixtriplet(xx[i]); } /* for (i = 0; i <= NS; i++) printout(xx[i]); */ return v; } void allocmemory() { xx = (triplet *) malloc(sizeof(triplet) * (NS + 1)); rr = (vector *) malloc(sizeof(vector) * (NS + 1)); erg = (double *) malloc(sizeof(double) * (NS + 1)); set = (struct Data *) malloc(sizeof(struct Data) * (NS + 1)); if (xx == NULL || rr == NULL || erg == NULL || set == NULL) { printf("Not enough memory\n"); exit(EXIT_FAILURE); } } /* The next function reads the data file. --- begin --- ## This is an example data file for jfm2full.c ## The following line specifies the total number of the model segments. ## It has to precede the rest. Number of segments 20 ## The next line defines the length of each segment in base pairs. ## All segments have the same length. This line is necessary ## for conversion of the result to moles per liter. Segment length 10.5 ## The next line specifies the segment diameter-to-length ratio. ## Any non-zero value turns on hardcore excluded-volume interactions ## and makes the chain self-avoiding. Use this feature only ## if you are absolutely sure that these effects are important ## because the program will run considerably slower. Diameter-to-length ratio 0.0 ## The next part contains detailed description of all the model ## chain segments. Each line contains 6 numbers: ## * number of sequential segments of a certain kind. ## * then THREE Euler angles in radians describing ## the equilibrium bend and twist in each segment of the kind. ## * then rms fluctuations of the twist. This is 0.272 ## for 10.5-bp segments if torsional rigidity is 2e-19 erg*cm. ## * then rms fluctuations of the bend. More precisely, this is the ## reciprocal square root rigidity. It equals to 0.268 ## for 10.5-bp segments if statistical length is 100 nm. Segments: number alpha beta gamma dac db 10 0.0 0.312 0.208 0.272 0.268 10 0.0 0.000 0.000 0.272 0.268 ## The next line specifies the desired linking number of the model ## chain in the circular topoisomer. Linking number 0 ## The next part is the restrictions on twist and bend ## angles and distance between the chain ends that provide ## an accurate calculation of the j-factor in most cases. ## They can also be used to remove certain restrictions altogether. ## First two numbers are in radians. The last is in base pairs. Conformation restriction: tw bd r 0.1 0.1 3.1 ## The next is the number of Metropolis movements for evaluation ## of each conditional probability. Increase it to improve the accuracy. Number of movements 10000 --- end --- */ void readdata() { int i = 1, numbr; char c[83]; double twist_eq = 0., arc_eq = 0.; double alpha, beta, gamma, sac, sb; while (fgets(c, sizeof(c), stdin) != NULL) { if (c[0] == '#' || c[0] == '\n') continue; if (NS == 0 && sscanf(c, "Number of segments %d", &NS) == 1) { allocmemory(); continue; } if (sscanf(c, "Segments: number alpha beta gamma dac db %*d") == EOF) { i = 1; twist_eq = 0.; arc_eq = 0.; continue; } if (sscanf(c, "%d %lf %lf %lf %lf %lf", &numbr, &alpha, &beta, &gamma, &sac, &sb) == 6) { for (; numbr > 0 && i <= NS; numbr--, i++) { eulerset(set[i].xi, alpha, beta, gamma); set[i].kac = 0.5 / (sac * sac); set[i].kb = 0.5 / (sb * sb); twist_eq += alpha + gamma; arc_eq += beta; } continue; } if (sscanf(c, "Diameter-to-length ratio %lf", &diam2) == 1) { diam2 *= diam2; continue; } if (sscanf(c, "Linking number %d", &Lk) == 1) continue; if (sscanf(c, "Segment length %lf", &bps) == 1) continue; if (sscanf(c, "Conformation restriction: tw bd r %*d") == EOF) { fscanf(stdin, "%lf %lf %lf", rdata, rdata + 1, rdata + 2); rdata[2] /= bps; /* to segment length */ continue; } if (sscanf(c, "Number of movements %d", &FREQ) == 1) continue; } NS = i - 1; /* report what was read */ printf("Segments = %d Arc = %g Tw = %g Lk = %d SqDiam = %g\n", NS, arc_eq / (2 * M_PI), twist_eq / (2 * M_PI), Lk, diam2); } int main(int argc, char *argv[]) { int i, exps = 9; double wolume, v, mv = 0.0, ev = 0.0; srand((unsigned) time(NULL)); if (argc > 1) if (freopen(argv[1], "rt", stdin) == NULL) { printf(VER "Can't open %s\nUsage: %s [infile]\n", argv[1], argv[0]); return EXIT_FAILURE; } readdata(); wolume = 1.5 / (rdata[0] * (1 - cos(rdata[1])) * rdata[2] * rdata[2] * rdata[2]); /* to moles per liter */ wolume /= 6.022e23 * 3.4e-9 * 3.4e-9 * 3.4e-9 * bps * bps * bps; /* reset parameters */ rdata[2] = rdata[2] * rdata[2]; rdata[1] = 1.0 - cos(rdata[1]); rdata[0] = 1.0 - cos(rdata[0]); printf("Simulations:\n"); for (i = 1; i <= exps; i++) { v = simulation(); printf("Probability = %e J-factor = %e M\n", v, v * wolume); /* printf("%f %f %f\n", writhe(), writhe1(), writhe2()); */ mv += (v - mv) / i; ev += (v * v - ev) / i; } ev = sqrt((ev - mv * mv) / exps); printf("Average J-factor = %e +/- %e M\n", mv * wolume, ev * wolume); return EXIT_SUCCESS; } /* ********** Useful functions */ double totenergy() { int i; double toten = 0.; for (i = 1; i <= NS; i++) toten += erg[i]; return toten; } double tottwist() { int i; double tottw = 0.; for (i = 0; i < NS; i++) tottw += euler_twist(xx[i + 1], xx[i]); return tottw; } /* Complex number multiplication with phase accumulation (2D-rotation) */ void phasiply(double *x0, double *y0, int *k0, double x, double y, int k) { double xo, yo, xn, yn; xo = *x0; yo = *y0; xn = xo * x - yo * y; yn = xo * y + yo * x; if (yo > 0.0 && y > 0.0 && yn < 0.0) ++k; if (yo < 0.0 && y < 0.0 && yn > 0.0) --k; *x0 = xn; *y0 = yn; *k0 += k; } /* avoid repeated calls to atan2 in the above tottwist implementation */ double twist() { int i, k0 = 0; double x0 = 1., y0 = 0., x, y; for (i = 0; i < NS; i++) { x = dotprod(xx[i + 1][0], xx[i][0]) + dotprod(xx[i + 1][1], xx[i][1]); y = dotprod(xx[i + 1][0], xx[i][1]) - dotprod(xx[i + 1][1], xx[i][0]); phasiply(&x0, &y0, &k0, x, y, 0); } return atan2(y0, x0) + 2 * M_PI * k0; } /* print out the chain coordinates for Vologodskii's view program */ void chain4view() { int i, j, l; static int k = 1; printf("\n%d %f %f", k, tottwist(), distance(rr[0], rr[NS])); for (l = 0, j = 0; j < 3; j++) for (i = 0; i <= NS; i++) { if (l % 5 == 0) putchar('\n'); l++; printf("%12.6f ", xx[i][2][j]); } k++; } /*** UNRELATED DEVELOPMENT OF OPTIMIZED WRITHE CALCULATION ALGORITHM ***/ /* writhe calculation algorithm that uses phasiply rather than calls atan2 */ double writhe() { int i, k, k0 = 0; vector r, ri, ri1, rk, rk1; double *zi, *zi1, *zk, *zk1; double x0 = 1., y0 = 0., dst, x, y; for (i = 0; i < NS - 1; i++) { zi = xx[i][2]; zi1 = (i) ? xx[i - 1][2] : xx[NS - 1][2]; for (k = i + 1; k < NS; k++) { zk = xx[k][2]; zk1 = xx[k - 1][2]; subtract(r, rr[k], rr[i]); normalize(r); crossprod(ri, zi, r); crossprod(ri1, zi1, r); crossprod(rk, r, zk); crossprod(rk1, r, zk1); x = dotprod(ri1, rk1); y = dotprod(ri1, zk1); phasiply(&x0, &y0, &k0, x, -y, 0); x = dotprod(ri, rk1); y = dotprod(ri, zk1); phasiply(&x0, &y0, &k0, x, y, 0); x = dotprod(ri1, rk); y = dotprod(ri1, zk); phasiply(&x0, &y0, &k0, x, y, 0); x = dotprod(ri, rk); y = dotprod(ri, zk); phasiply(&x0, &y0, &k0, x, -y, 0); //fprintf(stderr, "%g %g %g %d\n", dst, x0, y0, k0); } x = 1.0 / (fabs(x0) + fabs(y0)); x0 *= x; y0 *= x; } return atan2(y0, x0) / (2 * M_PI) + k0; } double writhe1() { int i, k; vector r, ri, ri1, rk, rk1; double *zi, *zi1, *zk, *zk1; double an, ad, dst, x1, x2, x3, x4, www = 0.; for (i = 0; i < NS - 1; i++) { zi = xx[i][2]; zi1 = (i) ? xx[i - 1][2] : xx[NS - 1][2]; for (k = i + 1; k < NS; k++) { zk = xx[k][2]; zk1 = xx[k - 1][2]; subtract(r, rr[k], rr[i]); dst = sqrt(dotprod(r, r)); crossprod(ri, r, zi); crossprod(ri1, r, zi1); crossprod(rk, r, zk); crossprod(rk1, r, zk1); an = 0.; ad = 1.; x1 = dotprod(ri1, rk1); x2 = dst * dotprod(ri1, zk1); x3 = ad * x1 - an * x2; x4 = ad * x2 + an * x1; ad = x3; an = x4; x1 = dotprod(ri, rk1); x2 = dst * dotprod(ri, zk1); x3 = ad * x1 - an * x2; x4 = ad * x2 + an * x1; ad = x3; an = x4; x1 = dotprod(ri1, rk); x2 = dst * dotprod(ri1, zk); x3 = ad * x1 - an * x2; x4 = ad * x2 + an * x1; ad = x3; an = x4; x1 = dotprod(ri, rk); x2 = dst * dotprod(ri, zk); x3 = ad * x1 - an * x2; x4 = ad * x2 + an * x1; ad = x3; an = x4; www += atan2(an, ad); } } return www / (2 * M_PI); } double writhe2() { int i, k; vector r, ri, ri1, rk, rk1; double *zi, *zi1, *zk, *zk1; double an, ad, dst, x, y, www = 0.; for (i = 0; i < NS - 1; i++) { zi = xx[i][2]; zi1 = (i) ? xx[i - 1][2] : xx[NS - 1][2]; for (k = i + 1; k < NS; k++) { zk = xx[k][2]; zk1 = xx[k - 1][2]; subtract(r, rr[k], rr[i]); normalize(r); crossprod(ri, zi, r); crossprod(ri1, zi1, r); crossprod(rk, r, zk); crossprod(rk1, r, zk1); x = dotprod(ri1, rk1); y = dotprod(ri1, zk1); www -= atan2(y, x); x = dotprod(ri, rk1); y = dotprod(ri, zk1); www += atan2(y, x); x = dotprod(ri1, rk); y = dotprod(ri1, zk); www += atan2(y, x); x = dotprod(ri, rk); y = dotprod(ri, zk); www -= atan2(y, x); } } return www / (2 * M_PI); }