[gmx-users] Calculation of unsaturated deuterium order parameters for POPC
chris.neale at utoronto.ca
chris.neale at utoronto.ca
Fri Jun 10 18:35:46 CEST 2011
I have patched g_order version 4.5.4 to compute the order parameters
in a different way. Despite trying, I don't understand the original
g_order implementation. In the version that I am providing, the
hydrogens are explicitly built and the coordinates of these hydrogens
are output to stdout.
Compared to the original g_order: there are large differences for the
2 carbons in the double bond and all single bonded carbons are the
same within what appears to be rounding error.
Compared to the VMD script that Thomas Piggot mentioned, the results
differ more, but only because that VMD script is using the real
hydrogen positions and not some with idealized geometry. There is a
larger difference for the double-bonded carbons, but if I use the
idealized hydrogen positions (output by my new version of g_order)
with the VMD .tcl script then I get the same values.
notes on the patch:
1. it is not complete for all options. Use it only with g_order -od
and with or without -unsat.
2. There are extraneous print statement that will provide you with the
coordinates of the atoms that are built. This is for debugging. Either
run with g_order >/dev/null or comment out the 3 fprintf lines near
the comments containing the word LOOK in capitals.
3. I did my best, but it could always be in error. Please use with caution.
To apply the patch (version 4.5.4)
cd src/tools
patch < g_order.patch
If it will not compile, then the spacing/hard-returns of the patch
must have been mangled in the process of posting. Contact me off-list
and I will provide you with a copy.
Here is the patch:
--- 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<DIM;m++){
+ ce[m]=(ac[m]+bc[m])/2;
}
+ rvec_add(x1[a[index[i]+j]],ce,e);
+ //UNCOMMENT THE NEXT LINE TO PRINT OUT THE POSITION OF THE
PSEUDO-ATOM FOR TESTING <<<<<<<<<<-------------- LOOK LOOK LOOK LOOK
+ fprintf(stdout,"ATOM E %f %f %f\n",e[0],e[1],e[2]);
+ //now take the angle between ce and z and this gets plugged into
S_CD=1/2<3cos^2(ang)-1>
+ 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<DIM;m++){
+ ce[m]=(ac[m]+bc[m])/2;
+ }
+ rvec_add(x1[a[index[i]+j]],ce,e);
+ dab=norm(ab);
+ //the eg vector is the cross product of the ac and ce vectors
(and ef is cross prod of bc and ce)
+ cprod(ac,ce,eg);
+ cprod(bc,ce,ef);
+ deg=norm(eg);
+ def=norm(ef);
+ //progress along eg by the correct distance from e to find g (and
along ef from e to find f)
+ for(m=0;m<DIM;m++){
+ g[m]=e[m]-eg[m]*(dab/deg)/2;
+ f[m]=e[m]-ef[m]*(dab/def)/2;
+ }
+ //now find the cg and cf vectors, which are the CD vectors for S_CD
+ rvec_sub(g,x1[a[index[i]+j]],cg);
+ rvec_sub(f,x1[a[index[i]+j]],cf);
+ //UNCOMMENT THE NEXT TWO LINES TO PRINT OUT THE POSITION OF THE
PSEUDO-ATOM FOR TESTING <<<<<<<<<<-------------- LOOK LOOK LOOK LOOK
+ fprintf(stdout,"ATOM G %f %f %f\n",g[0],g[1],g[2]);
+ fprintf(stdout,"ATOM F %f %f %f\n",f[0],f[1],f[2]);
+
+ //now take the angle between cg and z --- and also for cf and z
-- and this gets plugged into S_CD=1/2<3cos^2(ang)-1>
+ 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);
-------------- next part --------------
A non-text attachment was scrubbed...
Name: g_order.patch
Type: text/x-patch
Size: 7849 bytes
Desc: not available
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20110610/0aef9792/attachment.bin>
More information about the gromacs.org_gmx-users
mailing list