--- gmx_order.c 2011-03-04 06:10:44.000000000 -0500 +++ gmx_order.c 2011-06-09 17:57:11.586591000 -0400 @@ -379,6 +379,13 @@ real arcdist; gmx_rmpbc_t gpbc=NULL; +// BEGIN CN MOD ------------------------------------------------------------------------------------ +// variables for modification by CN +rvec ab,ac,bc,e,ce,eg,ef,g,f,cg,cf; +float dab,def,deg,gdot,fdot,edot; +// END CN MOD ------------------------------------------------------------------------------------ + + /* PBC added for center-of-mass vector*/ /* Initiate the pbc structure */ memset(&pbc,0,sizeof(pbc)); @@ -501,98 +508,66 @@ direction[0],direction[1],direction[2]);*/ } - if (bUnsat) { - /* Using convention for unsaturated carbons */ - /* first get Sz, the vector from Cn to Cn+1 */ - rvec_sub(x1[a[index[i+1]+j]], x1[a[index[i]+j]], dist); - length = norm(dist); - check_length(length, a[index[i]+j], a[index[i+1]+j]); - svmul(1/length, dist, Sz); - - /* this is actually the cosine of the angle between the double bond - and axis, because Sz is normalized and the two other components of - the axis on the bilayer are zero */ - if (use_unitvector) - { - sdbangle += acos(iprod(direction,Sz)); /*this can probably be optimized*/ - } - else - sdbangle += acos(Sz[axis]); - } else { - /* get vector dist(Cn-1,Cn+1) for tail atoms */ - rvec_sub(x1[a[index[i+1]+j]], x1[a[index[i-1]+j]], dist); - length = norm(dist); /* determine distance between two atoms */ - check_length(length, a[index[i-1]+j], a[index[i+1]+j]); - - svmul(1/length, dist, Sz); - /* Sz is now the molecular axis Sz, normalized and all that */ - } - - /* now get Sx. Sx is normal to the plane of Cn-1, Cn and Cn+1 so - we can use the outer product of Cn-1->Cn and Cn+1->Cn, I hope */ - rvec_sub(x1[a[index[i+1]+j]], x1[a[index[i]+j]], tmp1); - rvec_sub(x1[a[index[i-1]+j]], x1[a[index[i]+j]], tmp2); - cprod(tmp1, tmp2, Sx); - svmul(1/norm(Sx), Sx, Sx); - - /* now we can get Sy from the outer product of Sx and Sz */ - cprod(Sz, Sx, Sy); - svmul(1/norm(Sy), Sy, Sy); - - /* the square of cosine of the angle between dist and the axis. - Using the innerproduct, but two of the three elements are zero - Determine the sum of the orderparameter of all atoms in group - */ - if (use_unitvector) - { - cossum[XX] = sqr(iprod(Sx,direction)); /* this is allowed, since Sa is normalized */ - cossum[YY] = sqr(iprod(Sy,direction)); - cossum[ZZ] = sqr(iprod(Sz,direction)); - } - else - { - cossum[XX] = sqr(Sx[axis]); /* this is allowed, since Sa is normalized */ - cossum[YY] = sqr(Sy[axis]); - cossum[ZZ] = sqr(Sz[axis]); +// BEGIN CN MOD ------------------------------------------------------------------------------------ + if (bUnsat) { + //exactly as for the case without unsat, but here we select the ce vector to represent the hydrogen + // place the 2 imaginary hydrogen atoms by enscribing a tetrahedron in a cube + // a,b,c are the 3 backbone atoms, where a and b are in two cube corners and c is at the center + rvec_sub(x1[a[index[i-1]+j]], x1[a[index[i+1]+j]], ab); + rvec_sub(x1[a[index[i]+j]], x1[a[index[i+1]+j]], ac); + rvec_sub(x1[a[index[i]+j]], x1[a[index[i-1]+j]], bc); + // e is a point in the middle of the face opposite from the ab face + for(m=0;m + edot=ce[2]/norm(ce); + + // using frameorder[ZZ] to store the result + frameorder[ZZ] += 0.5 * (3 * edot*edot - 1); + + } else { + // place the 2 imaginary hydrogen atoms by enscribing a tetrahedron in a cube + // a,b,c are the 3 backbone atoms, where a and b are in two cube corners and c is at the center + rvec_sub(x1[a[index[i-1]+j]], x1[a[index[i+1]+j]], ab); + rvec_sub(x1[a[index[i]+j]], x1[a[index[i+1]+j]], ac); + rvec_sub(x1[a[index[i]+j]], x1[a[index[i-1]+j]], bc); + // e is a point in the middle of the face opposite from the ab face + for(m=0;m + gdot=cg[2]/norm(cg); + fdot=cf[2]/norm(cf); + + // using frameorder[ZZ] to store the result -- mult by 0.25 instead of 0.5 because we are adding two values + frameorder[ZZ] += 0.25 * (3 * gdot*gdot - 1); + frameorder[ZZ] += 0.25 * (3 * fdot*fdot - 1); + } - for (m = 0; m < DIM; m++) - frameorder[m] += 0.5 * (3 * cossum[m] - 1); - - if (bSliced) { - /* get average coordinate in box length for slicing, - determine which slice atom is in, increase count for that - slice. slFrameorder and slOrder are reals, not - rvecs. Only the component [axis] of the order tensor is - kept, until I find it necessary to know the others too - */ - - z1 = x1[a[index[i-1]+j]][axis]; - z2 = x1[a[index[i+1]+j]][axis]; - z_ave = 0.5 * (z1 + z2); - if (z_ave < 0) - z_ave += box[axis][axis]; - if (z_ave > box[axis][axis]) - z_ave -= box[axis][axis]; - - slice = (int)(0.5 + (z_ave / (*slWidth))) - 1; - slCount[slice]++; /* determine slice, increase count */ - - slFrameorder[slice] += 0.5 * (3 * cossum[axis] - 1); - } - else if (permolecule) - { - /* store per-molecule order parameter - * To just track single-axis order: (*slOrder)[j][i] += 0.5 * (3 * iprod(cossum,direction) - 1); - * following is for Scd order: */ - (*slOrder)[j][i] += -1* (0.3333 * (3 * cossum[XX] - 1) + 0.3333 * 0.5 * (3 * cossum[YY] - 1)); - } - if (distcalc) - { - /* bin order parameter by arc distance from reference group*/ - arcdist=acos(iprod(dvec,direction)); - (*distvals)[j][i]+=arcdist; - } +// END CN MOD ------------------------------------------------------------------------------------ } /* end loop j, over all atoms in group */ for (m = 0; m < DIM; m++) @@ -702,8 +677,11 @@ for (atom = 1; atom < ngrps - 1; atom++) { fprintf(ord,"%12d %12g %12g %12g\n", atom, order[atom][XX], order[atom][YY], order[atom][ZZ]); - fprintf(slOrd,"%12d %12g\n", atom, -1 * (0.6667 * order[atom][XX] + - 0.333 * order[atom][YY])); +// BEGIN CN MOD ------------------------------------------------------------------------------------ + // will place the order parameter value directly in the ZZ index. + // not sure why I need to take the negative value.... + fprintf(slOrd,"%12d %12g\n", atom, -1 * order[atom][ZZ]); +// END CN MOD ------------------------------------------------------------------------------------ } ffclose(ord);