i, j, m, k, l, teller = 0,
slice, /* current slice number */
nr_frames = 0,
- *slCount; /* nr. of atoms in one slice */
+ *slCount; /* nr. of atoms in one slice */
real dbangle = 0, /* angle between double bond and axis */
sdbangle = 0; /* sum of these angles */
gmx_bool use_unitvector = FALSE; /* use a specified unit vector instead of axis to specify unit normal*/
atom_id *comidx = NULL, *distidx = NULL;
char *grpname = NULL;
t_pbc pbc;
- real arcdist;
+ real arcdist, tmpdist;
gmx_rmpbc_t gpbc = NULL;
/* PBC added for center-of-mass vector*/
use_unitvector = TRUE;
fprintf(stderr, "Select an index group to calculate the radial membrane normal\n");
get_index(&top->atoms, radfn, 1, &comsize, &comidx, &grpname);
- if (distcalc)
+ }
+ if (distcalc)
+ {
+ if (grpname != NULL)
{
- if (grpname != NULL)
- {
- sfree(grpname);
- }
- fprintf(stderr, "Select an index group to use as distance reference\n");
- get_index(&top->atoms, radfn, 1, &distsize, &distidx, &grpname);
- bSliced = FALSE; /*force slices off*/
+ sfree(grpname);
}
+ fprintf(stderr, "Select an index group to use as distance reference\n");
+ get_index(&top->atoms, radfn, 1, &distsize, &distidx, &grpname);
+ bSliced = FALSE; /*force slices off*/
}
if (use_unitvector && bSliced)
nslices, *slWidth);
}
+
#if 0
nr_tails = index[1] - index[0];
fprintf(stderr, "Number of elements in first group: %d\n", nr_tails);
so for DPPC ngrps = 16 and i runs from 1 to 14, including 14
*/
+
if (radial)
{
/*center-of-mass determination*/
rvec_inc(dist, x1[distidx[j]]);
}
svmul(1.0/distsize, dref, dref);
- pbc_dx(&pbc, dref, com, dvec);
- unitv(dvec, dvec);
+ if (radial)
+ {
+ pbc_dx(&pbc, dref, com, dvec);
+ unitv(dvec, dvec);
+ }
}
for (i = 1; i < ngrps - 1; i++)
the axis on the bilayer are zero */
if (use_unitvector)
{
- sdbangle += acos(iprod(direction, Sz)); /*this can probably be optimized*/
+ sdbangle += gmx_angle(direction, Sz); /*this can probably be optimized*/
}
else
{
}
if (distcalc)
{
- /* bin order parameter by arc distance from reference group*/
- arcdist = acos(iprod(dvec, direction));
- (*distvals)[j][i] += arcdist;
+ if (radial)
+ {
+ /* bin order parameter by arc distance from reference group*/
+ arcdist = gmx_angle(dvec,direction);
+ (*distvals)[j][i] += arcdist;
+ }
+ else if (i == 1)
+ {
+ /* Want minimum lateral distance to first group calculated */
+ tmpdist = trace(box); /* should be max value */
+ for (k=0;k<distsize;k++)
+ {
+ pbc_dx(&pbc, x1[distidx[k]], x1[a[index[i]+j]], dvec);
+ /* at the moment, just remove dvec[axis] */
+ dvec[axis] = 0;
+ tmpdist = min(tmpdist, norm2(dvec));
+ }
+ //fprintf(stderr, "Min dist %f; trace %f\n", tmpdist, trace(box));
+ (*distvals)[j][i]+=sqrt(tmpdist);
+ }
}
} /* end loop j, over all atoms in group */
{ "-radial", FALSE, etBOOL, {&radial},
"Compute a radial membrane normal" },
{ "-calcdist", FALSE, etBOOL, {&distcalc},
- "Compute distance from a reference (currently defined only for radial and permolecule)" },
+ "Compute distance from a reference" },
};
rvec *order; /* order par. for each atom */
fprintf(stderr, "Calculating Scd order parameters for each of %d molecules\n", nslices);
}
- if (distcalc)
+ if (radial)
{
- radial = TRUE;
fprintf(stderr, "Calculating radial distances\n");
if (!permolecule)
{