#include "viewit.h"
#include "xvgr.h"
#include "gmx_ana.h"
+#include "gmx_sort.h"
enum { euSel,euRect, euTric, euCompact, euNR};
-static real calc_isquared(int nmol,rvec m_com[],rvec m_shift[],rvec clust_com)
+
+static int
+sort_comdist2(void *thunk, const void *a, const void *b)
{
- real Isq=0;
- int i;
- rvec m0,dx;
-
- for(i=0; (i<nmol); i++) {
- rvec_add(m_com[i],m_shift[i],m0);
- rvec_sub(m0,clust_com,dx);
- Isq += iprod(dx,dx);
+ /* Thunk should point to a real array with the distance to the cluster COM for each molecule,
+ * a/b point to integers that refer to the molecule number.
+ */
+ real *pcomdist2 = thunk;
+ int ia = * (int *)a;
+ int ib = * (int *)b;
+ int rc;
+
+ if(pcomdist2[ia]<pcomdist2[ib])
+ {
+ rc=-1;
+ }
+ else if (pcomdist2[ia]>pcomdist2[ib])
+ {
+ rc=1;
}
- return Isq;
+ else
+ {
+ rc=0;
+ }
+ return rc;
}
+
static void calc_pbc_cluster(int ecenter,int nrefat,t_topology *top,int ePBC,
rvec x[],atom_id index[],
- rvec clust_com,matrix box)
+ rvec clust_com,matrix box, rvec clustercenter)
{
- const real tol=1e-3;
- gmx_bool bChanged;
- int m,i,j,j0,j1,jj,ai,iter,is;
- real fac,Isq,min_dist2;
- rvec dx,ddx,xtest,xrm,box_center;
- int nmol,nmol_cl,imol_center;
+ int m,i,j,j0,j1,jj,ai,aj;
+ int imin,jmin;
+ real fac,min_dist2;
+ rvec dx,xtest,box_center;
+ int nmol,imol_center;
atom_id *molind;
- gmx_bool *bMol,*bTmp;
- rvec *m_com,*m_shift,m0;
+ gmx_bool *bMol,*bTmp;
+ rvec *m_com,*m_shift;
t_pbc pbc;
-
+ real *com_dist2;
+ int *cluster;
+ int *added;
+ int ncluster,nadded;
+ real tmp_r2;
+
calc_box_center(ecenter,box,box_center);
/* Initiate the pbc structure */
snew(bMol,nmol);
snew(m_com,nmol);
snew(m_shift,nmol);
+ snew(cluster,nmol);
+ snew(added,nmol);
snew(bTmp,top->atoms.nr);
- nmol_cl = 0;
+
for(i=0; (i<nrefat); i++) {
/* Mark all molecules in the index */
ai = index[i];
*/
min_dist2 = 10*sqr(trace(box));
imol_center = -1;
- for(i=0; (i<nmol); i++) {
- for(j=molind[i]; (j<molind[i+1]); j++) {
+ ncluster = 0;
+ for(i=0; i<nmol; i++)
+ {
+ for(j=molind[i]; j<molind[i+1]; j++)
+ {
if (bMol[i] && !bTmp[j])
- gmx_fatal(FARGS,"Molecule %d marked for clustering but not atom %d",
- i+1,j+1);
+ {
+ gmx_fatal(FARGS,"Molecule %d marked for clustering but not atom %d in it - check your index!",i+1,j+1);
+ }
else if (!bMol[i] && bTmp[j])
- gmx_fatal(FARGS,"Atom %d marked for clustering but not molecule %d",
- j+1,i+1);
- else if (bMol[i]) {
- /* Compute center of geometry of molecule */
+ {
+ gmx_fatal(FARGS,"Atom %d marked for clustering but not molecule %d - this is an internal error...",j+1,i+1);
+ }
+ else if (bMol[i])
+ {
+ /* Make molecule whole, move 2nd and higher atom to same periodicity as 1st atom in molecule */
+ if(j>molind[i])
+ {
+ pbc_dx(&pbc,x[j],x[molind[i]],dx);
+ rvec_add(x[molind[i]],dx,x[j]);
+ }
+ /* Compute center of geometry of molecule - m_com[i] was zeroed when we did snew() on it! */
rvec_inc(m_com[i],x[j]);
}
}
- if (bMol[i]) {
+ if (bMol[i])
+ {
/* Normalize center of geometry */
fac = 1.0/(molind[i+1]-molind[i]);
for(m=0; (m<DIM); m++)
m_com[i][m] *= fac;
/* Determine which molecule is closest to the center of the box */
pbc_dx(&pbc,box_center,m_com[i],dx);
- if (iprod(dx,dx) < min_dist2) {
- min_dist2 = iprod(dx,dx);
+ tmp_r2=iprod(dx,dx);
+
+ if (tmp_r2 < min_dist2)
+ {
+ min_dist2 = tmp_r2;
imol_center = i;
}
- nmol_cl++;
+ cluster[ncluster++]=i;
}
}
sfree(bTmp);
- if (nmol_cl <= 0) {
+ if (ncluster <= 0)
+ {
fprintf(stderr,"No molecules selected in the cluster\n");
return;
- } else if (imol_center == -1) {
+ }
+ else if (imol_center == -1)
+ {
fprintf(stderr,"No central molecules could be found\n");
return;
}
-
-
- /* First calculation is incremental */
- clear_rvec(clust_com);
- Isq = 0;
- for(i=m=0; (i<nmol); i++) {
- /* Check whether this molecule is part of the cluster */
- if (bMol[i]) {
- if ((i > 0) && (m > 0)) {
- /* Compute center of cluster by dividing by number of molecules */
- svmul(1.0/m,clust_com,xrm);
- /* Distance vector between molecular COM and cluster */
- pbc_dx(&pbc,m_com[i],xrm,dx);
- rvec_add(xrm,dx,xtest);
- /* xtest is now the image of m_com[i] that is closest to clust_com */
- rvec_inc(clust_com,xtest);
- rvec_sub(xtest,m_com[i],m_shift[i]);
- }
- else {
- rvec_inc(clust_com,m_com[i]);
- }
- m++;
- }
- }
- assert(m == nmol_cl);
- svmul(1/nmol_cl,clust_com,clust_com);
- put_atom_in_box(box,clust_com);
-
- /* Now check if any molecule is more than half the box from the COM */
- if (box) {
- iter = 0;
- do {
- bChanged = FALSE;
- for(i=0; (i<nmol) && !bChanged; i++) {
- if (bMol[i]) {
- /* Sum com and shift from com */
- rvec_add(m_com[i],m_shift[i],m0);
- pbc_dx(&pbc,m0,clust_com,dx);
- rvec_add(clust_com,dx,xtest);
- rvec_sub(xtest,m0,ddx);
- if (iprod(ddx,ddx) > tol) {
- /* Here we have used the wrong image for contributing to the COM */
- rvec_sub(xtest,m_com[i],m_shift[i]);
- for(j=0; (j<DIM); j++)
- clust_com[j] += (xtest[j]-m0[j])/nmol_cl;
- bChanged = TRUE;
- }
+
+ nadded = 0;
+ added[nadded++] = imol_center;
+ bMol[imol_center] = FALSE;
+
+ while (nadded<ncluster)
+ {
+ /* Find min distance between cluster molecules and those remaining to be added */
+ min_dist2 = 10*sqr(trace(box));
+ imin = -1;
+ jmin = -1;
+ /* Loop over added mols */
+ for(i=0;i<nadded;i++)
+ {
+ ai = added[i];
+ /* Loop over all mols */
+ for(j=0;j<ncluster;j++)
+ {
+ aj = cluster[j];
+ /* check those remaining to be added */
+ if(bMol[aj])
+ {
+ pbc_dx(&pbc,m_com[aj],m_com[ai],dx);
+ tmp_r2=iprod(dx,dx);
+ if (tmp_r2 < min_dist2)
+ {
+ min_dist2 = tmp_r2;
+ imin = ai;
+ jmin = aj;
+ }
}
}
- Isq = calc_isquared(nmol,m_com,m_shift,clust_com);
- put_atom_in_box(box,clust_com);
-
- if (bChanged && (iter > 0))
- printf("COM: %8.3f %8.3f %8.3f iter = %d Isq = %8.3f\n",
- clust_com[XX],clust_com[YY],clust_com[ZZ],iter,Isq);
- iter++;
- } while (bChanged);
- }
- /* Now transfer the shift to all atoms in the molecule */
- for(i=0; (i<nmol); i++) {
- if (bMol[i]) {
- for(j=molind[i]; (j<molind[i+1]); j++)
- rvec_inc(x[j],m_shift[i]);
}
+
+ /* Add the best molecule */
+ added[nadded++] = jmin;
+ bMol[jmin] = FALSE;
+ /* Calculate the shift from the ai molecule */
+ pbc_dx(&pbc,m_com[jmin],m_com[imin],dx);
+ rvec_add(m_com[imin],dx,xtest);
+ rvec_sub(xtest,m_com[jmin],m_shift[jmin]);
+ rvec_inc(m_com[jmin],m_shift[jmin]);
+
+ for(j=molind[jmin]; j<molind[jmin+1]; j++)
+ {
+ rvec_inc(x[j],m_shift[jmin]);
+ }
+ fprintf(stdout,"\rClustering iteration %d of %d...",nadded,ncluster);
+ fflush(stdout);
}
+
+ sfree(added);
+ sfree(cluster);
sfree(bMol);
sfree(m_com);
sfree(m_shift);
+
+ fprintf(stdout,"\n");
}
static void put_molecule_com_in_box(int unitcell_enum,int ecenter,
"results if you in fact have a cluster. Luckily that can be checked",
"afterwards using a trajectory viewer. Note also that if your molecules",
"are broken this will not work either.[BR]",
+ "The separate option [TT]-clustercenter[tt] can be used to specify an",
+ "approximate center for the cluster. This is useful e.g. if you have",
+ "two big vesicles, and you want to maintain their relative positions.[BR]",
"[TT]* whole[tt] only makes broken molecules whole.[PAR]",
"Option [TT]-ur[tt] sets the unit cell representation for options",
"[TT]mol[tt], [TT]res[tt] and [TT]atom[tt] of [TT]-pbc[tt].",
static char *exec_command=NULL;
static real dropunder=0,dropover=0;
static gmx_bool bRound=FALSE;
-
+ static rvec clustercenter = {0,0,0};
+
t_pargs
pa[] =
{
{ "-center", FALSE, etBOOL,
{ &bCenter }, "Center atoms in box" },
{ "-boxcenter", FALSE, etENUM,
- { center_opt }, "Center for -pbc and -center" },
+ { center_opt }, "Center for -pbc and -center" },
{ "-box", FALSE, etRVEC,
{ newbox },
"Size for new cubic box (default: read from input)" },
+ { "-clustercenter", FALSE, etRVEC,
+ { clustercenter },
+ "Optional starting point for pbc cluster option" },
{ "-trans", FALSE, etRVEC,
{ trans },
"All coordinates will be translated by trans. This "
else if (bCluster) {
rvec com;
- calc_pbc_cluster(ecenter,ifit,&top,ePBC,fr.x,ind_fit,com,fr.box);
+ calc_pbc_cluster(ecenter,ifit,&top,ePBC,fr.x,ind_fit,com,fr.box,clustercenter);
}
if (bPFit) {
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff