Hey :)
As there seems to be interest :)
Below is the diff between native and my version of gmx_trjconv.c 4.5, the
complete source file is attached. I'd fancy comments and suggestions.
Hope it helps,
Tsjerk
###
82a83,208
static void taw_pbc_cluster(int nrefat,t_topology *top,rvec *x, atom_id
*ndx,
matrix box,real cutoff)
{
int i,j,id,na,nc,cid,nr,m,*cluster;
atom_id *active,*done,*rest;
matrix S,T,A;
rvec *y,rd;
real ax2,ay2,az2,axy,axz,ayz,d,cut2;
gmx_bool bClustered;
cut2 = cutoff*cutoff;
transpose(box,S);
m_inv(S,T);
ax2 = iprod( S[0], S[0] );
ay2 = iprod( S[1], S[1] );
az2 = iprod( S[2], S[2] );
axy = 2*iprod( S[0], S[1] );
axz = 2*iprod( S[0], S[2] );
ayz = 2*iprod( S[1], S[2] );
snew(rest, nrefat);
snew(active, nrefat);
snew(done, nrefat);
snew(cluster,nrefat);
snew(y, nrefat);
nr = nrefat;
nc = m = 0;
cid = 0;
for (i=0; i<nr; i++)
{
rest[i] = i;
/* Transform to box coordinates */
mvmul(T,x[ndx[i]],y[i]);
}
while (nr)
{
if (nr<0)
exit(1);
/*
If we get here we have _nc_ clustered atoms,
Time for a new cluster.
Pick one atom from the rest group and add it to the clustered ones.
Rest counter nr is decremented.
The new active atom is not counted as clustered yet.
The cluster id was already set.
When entering the following while loop we have exactly one active
atom.
*/
done[nc] = rest[--nr];
cluster[nc] = cid;
na = 1;
m = 0;
while (na)
{
/*
Run over active atoms.
The active atoms run from nc to nc+na
*/
/* Iterate over atoms in rest group */
for (i=0; i<nr; i++)
{
bClustered = FALSE;
/* For each atom in rest group check distance from active group */
for (j=nc; j<nc+na; j++)
{
/* Distance in periodic system in box coordinates */
rvec_sub( y[rest[i]], y[done[j]], rd );
/* Truncate coordinates for minimial distances */
rd[0] -= floor( rd[0] + 0.5 );
rd[1] -= floor( rd[1] + 0.5 );
rd[2] -= floor( rd[2] + 0.5 );
/* Distance follows from d = r'S'Sr = r'Ar = */
d =
rd[0]*(rd[0]*ax2+rd[1]*axy+rd[2]*axz)+rd[1]*(rd[1]*ay2+rd[2]*ayz)+rd[2]*rd[2]*az2;
/* If atom i is within the cutoff distance of j, it is part of the
same cluster */
if (d<cut2)
{
/* Add the atom to the active cluster group */
id = nc+na+m;
done[id] = rest[i];
cluster[id] = cid;
m++;
/* Relocate the atoms such that it is positioned properly */
rvec_add(rd,y[done[j]],y[done[id]]);
/* Now pass on to the next atom in the rest group */
break;
}
/*
We only get here if the atom could not be assigned to a cluster.
Then we shift the atom down in the rest array.
*/
rest[i-m] = rest[i];
}
}
/* Subtract the number of atoms taken from the rest group */
nr -= m;
/* Now mark all previously active atoms clustered */
nc = nc + na;
/* And mark all atoms added active */
na = m;
m = 0;
} /* while (na) */
cid++;
} /* while (nr) */
fprintf(stderr,"Found %d clusters",cid);
sfree(rest);
sfree(active);
sfree(done);
sfree(cluster);
for (i=0; i<nrefat; i++)
{
mvmul(S,y[i],x[ndx[i]]);
}
sfree(y);
}
638c764
< static real dropunder=0,dropover=0;
---
static real dropunder=0,dropover=0,dclus=0.5;
722c848,851
< "coarse grained ones" } };
---
"coarse grained ones" },
{ "-dclus", FALSE, etREAL,
{ &dclus },
"Cutoff distance for clustering" } };
932c1061,1062
< if (0 == top.mols.nr && (bCluster || bPBCcomMol))
---
/* if (0 == top.mols.nr && (bCluster || bPBCcomMol))*/
if (0 == top.mols.nr && bPBCcomMol)
1197c1327,1328
<
---
taw_pbc_cluster(ifit,&top,fr.x,ind_fit,fr.box,dclus);
/*
1198a1330
*/
Hi Tsjerk,
I am doing protein-micelle simultions and could you also send me your
modified trjconv code?
Thank you very much!
best regards,
Jianguo
------------------------------
*From:* Tsjerk Wassenaar <tsjerkw at gmail.com>
*To:* Discussion list for GROMACS users <gmx-users at gromacs.org>
*Sent:* Thursday, 14 April 2011 23:44:07
*Subject:* Re: [gmx-users] micelles and trjconv -pbc cluster
Hi George,
Recently I wrote an alternative, non-iterative clustering routine, that
does not suffer from convergence failures. If you want, I can send you the
modified trjconv source code. Note that it does not bother about the center
of mass of the cluster, but just builds a network of neighbours, until there
are no more. If you're clusters are not periodic, it won't matter.
Let me know...
Cheers,
Tsjerk
Post by jim jackDear GROMACS users,
I am trying to simulate an SDS micelle in water. As simulation time
goes by, the micelle approaches the edge of the box and consequently some of
these molecules get in from the other side. This leads to incorrect radius
of gyration, eccentricity, etc. A solution to this problem is the option trjconv
-pbc cluster as described in the page
http://www.gromacs.org/Documentation/How-tos/Micelle_Clustering. In this
case, the problem is that it takes a lot of time and a huge file (several
GB) is created due to this procedure. Is there any other alternative?
Thanks in advance
George Koros
I don't think that the cluster option always converges. You could, if
your micelle is intact at frame 0, first do trjconv -pbc nojump, then
optionally trjconv -center. That should give a trajectory from which you
could calculate the radius of gyration. If SDS molecules occationally leave
the micelle and recombine with a priodic, then you might have a problem with
the suggested approach.
--
-----------------------------------------------
Erik Marklund, PhD student
Dept. of Cell and Molecular Biology, Uppsala University.
Husargatan 3, Box 596, 75124 Uppsala, Sweden
phone: +46 18 471 4537 fax: +46 18 511 755erikm at xray.bmc.uu.se http://folding.bmc.uu.se/
--
gmx-users mailing list gmx-users at gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-request at gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
--
Tsjerk A. Wassenaar, Ph.D.
post-doctoral researcher
Molecular Dynamics Group
* Groningen Institute for Biomolecular Research and Biotechnology
* Zernike Institute for Advanced Materials
University of Groningen
The Netherlands
--
Tsjerk A. Wassenaar, Ph.D.
post-doctoral researcher
Molecular Dynamics Group
* Groningen Institute for Biomolecular Research and Biotechnology
* Zernike Institute for Advanced Materials
University of Groningen
The Netherlands
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20110415/2f4b3043/attachment.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: gmx_trjconv.c.gz
Type: application/x-gzip
Size: 17278 bytes
Desc: not available
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20110415/2f4b3043/attachment.gz>