Hey :)<br><br>As there seems to be interest :)<br>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.<br><br>Hope it helps,<br>
<br>Tsjerk<br><br>###<br><br>82a83,208<br>> static void taw_pbc_cluster(int nrefat,t_topology *top,rvec *x, atom_id *ndx,<br>> matrix box,real cutoff)<br>> {<br>> int i,j,id,na,nc,cid,nr,m,*cluster;<br>
> atom_id *active,*done,*rest;<br>> matrix S,T,A;<br>> rvec *y,rd;<br>> real ax2,ay2,az2,axy,axz,ayz,d,cut2;<br>> gmx_bool bClustered;<br>> <br>> cut2 = cutoff*cutoff;<br>> <br>> transpose(box,S);<br>
> m_inv(S,T);<br>> ax2 = iprod( S[0], S[0] );<br>> ay2 = iprod( S[1], S[1] );<br>> az2 = iprod( S[2], S[2] );<br>> axy = 2*iprod( S[0], S[1] );<br>> axz = 2*iprod( S[0], S[2] );<br>> ayz = 2*iprod( S[1], S[2] );<br>
> <br>> snew(rest, nrefat);<br>> snew(active, nrefat);<br>> snew(done, nrefat);<br>> snew(cluster,nrefat);<br>> snew(y, nrefat);<br>> <br>> nr = nrefat;<br>> nc = m = 0;<br>
> cid = 0;<br>> for (i=0; i<nr; i++)<br>> {<br>> rest[i] = i;<br>> /* Transform to box coordinates */<br>> mvmul(T,x[ndx[i]],y[i]);<br>> }<br>> while (nr)<br>> {<br>> if (nr<0)<br>
> exit(1);<br>> <br>> /* <br>> If we get here we have _nc_ clustered atoms, <br>> but none which has neighbours in the rest group:<br>> Time for a new cluster.<br>> Pick one atom from the rest group and add it to the clustered ones.<br>
> Rest counter nr is decremented.<br>> The new active atom is not counted as clustered yet.<br>> The cluster id was already set.<br>> When entering the following while loop we have exactly one active atom.<br>
> */<br>> done[nc] = rest[--nr];<br>> cluster[nc] = cid;<br>> na = 1;<br>> m = 0;<br>> while (na)<br>> {<br>> /*<br>> Run over active atoms.<br>> The active atoms run from nc to nc+na<br>
> */<br>> <br>> /* Iterate over atoms in rest group */<br>> for (i=0; i<nr; i++)<br>> {<br>> bClustered = FALSE;<br>> <br>> /* For each atom in rest group check distance from active group */<br>
> for (j=nc; j<nc+na; j++)<br>> {<br>> /* Distance in periodic system in box coordinates */<br>> rvec_sub( y[rest[i]], y[done[j]], rd );<br>> /* Truncate coordinates for minimial distances */<br>
> rd[0] -= floor( rd[0] + 0.5 );<br>> rd[1] -= floor( rd[1] + 0.5 );<br>> rd[2] -= floor( rd[2] + 0.5 );<br>> /* Distance follows from d = r'S'Sr = r'Ar = */<br>> 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;<br>
> <br>> /* If atom i is within the cutoff distance of j, it is part of the same cluster */ <br>> if (d<cut2)<br>> {<br>> /* Add the atom to the active cluster group */ <br>
> id = nc+na+m;<br>> done[id] = rest[i];<br>> cluster[id] = cid;<br>> m++;<br>> /* Relocate the atoms such that it is positioned properly */<br>> rvec_add(rd,y[done[j]],y[done[id]]);<br>
> /* Now pass on to the next atom in the rest group */<br>> break;<br>> }<br>> <br>> /* <br>> We only get here if the atom could not be assigned to a cluster.<br>> Then we shift the atom down in the rest array. <br>
> */<br>> rest[i-m] = rest[i];<br>> }<br>> }<br>> /* Subtract the number of atoms taken from the rest group */<br>> nr -= m;<br>> /* Now mark all previously active atoms clustered */<br>
> nc = nc + na;<br>> /* And mark all atoms added active */<br>> na = m;<br>> m = 0;<br>> } /* while (na) */<br>> cid++;<br>> } /* while (nr) */<br>> fprintf(stderr,"Found %d clusters",cid);<br>
> <br>> <br>> sfree(rest);<br>> sfree(active);<br>> sfree(done);<br>> sfree(cluster);<br>> <br>> for (i=0; i<nrefat; i++)<br>> {<br>> mvmul(S,y[i],x[ndx[i]]);<br>> }<br>
> <br>> sfree(y);<br>> }<br>> <br>638c764<br>< static real dropunder=0,dropover=0;<br>---<br>> static real dropunder=0,dropover=0,dclus=0.5;<br>722c848,851<br>< "coarse grained ones" } };<br>
---<br>> "coarse grained ones" },<br>> { "-dclus", FALSE, etREAL,<br>> { &dclus },<br>> "Cutoff distance for clustering" } };<br>
932c1061,1062<br>< if (0 == <a href="http://top.mols.nr">top.mols.nr</a> && (bCluster || bPBCcomMol))<br>---<br>> /* if (0 == <a href="http://top.mols.nr">top.mols.nr</a> && (bCluster || bPBCcomMol))*/<br>
> if (0 == <a href="http://top.mols.nr">top.mols.nr</a> && bPBCcomMol)<br>1197c1327,1328<br>< <br>---<br>> taw_pbc_cluster(ifit,&top,fr.x,ind_fit,fr.box,dclus);<br>> /*<br>
1198a1330<br>> */<br><br><br><br><br><div class="gmail_quote">On Fri, Apr 15, 2011 at 3:16 AM, Jianguo Li <span dir="ltr"><<a href="mailto:ljggmx@yahoo.com.sg">ljggmx@yahoo.com.sg</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin: 0pt 0pt 0pt 0.8ex; border-left: 1px solid rgb(204, 204, 204); padding-left: 1ex;"><div><div style="font-family: times new roman,new york,times,serif; font-size: 12pt;"><div>Hi Tsjerk,<br>
<br>I am doing protein-micelle simultions and could you also send me your modified trjconv code?<br>Thank you very much!<br><br>best regards,<br>Jianguo<br></div><div style="font-family: times new roman,new york,times,serif; font-size: 12pt;">
<br><div style="font-family: times new roman,new york,times,serif; font-size: 12pt;"><font face="Tahoma" size="2"><hr size="1"><div class="im"><b><span style="font-weight: bold;">From:</span></b> Tsjerk Wassenaar <<a href="mailto:tsjerkw@gmail.com" target="_blank">tsjerkw@gmail.com</a>><br>
</div><div class="im"><b><span style="font-weight: bold;">To:</span></b> Discussion list for GROMACS users <<a href="mailto:gmx-users@gromacs.org" target="_blank">gmx-users@gromacs.org</a>><br></div><b><span style="font-weight: bold;">Sent:</span></b> Thursday, 14 April 2011 23:44:07<br>
<b><span style="font-weight: bold;">Subject:</span></b> Re: [gmx-users] micelles and trjconv -pbc
cluster<br></font><div><div></div><div class="h5"><br>Hi George,<br><br>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.<br>
<br>Let me know...<br><br>Cheers,<br><br>Tsjerk<br><br><div class="gmail_quote">On Thu, Apr 14, 2011 at 4:48 PM, Erik Marklund <span dir="ltr"><<a rel="nofollow" href="mailto:erikm@xray.bmc.uu.se" target="_blank">erikm@xray.bmc.uu.se</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin: 0pt 0pt 0pt 0.8ex; border-left: 1px solid rgb(204, 204, 204); padding-left: 1ex;">
<div>
jim jack skrev 2011-04-14 16.42:
<div><div></div><div><blockquote type="cite">
<table border="0" cellpadding="0" cellspacing="0">
<tbody>
<tr>
<td style="font: inherit;" valign="top">Dear GROMACS users,
<br>
<br>
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 <span style="font-style: italic;">trjconv -pbc
cluste</span>r as described in the page <a rel="nofollow">http://www.gromacs.org/Documentation/How-tos/Micelle_Clustering</a>.
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?<br>
<br>
Thanks in advance<br>
<br>
George Koros<br>
<br>
</td>
</tr>
</tbody>
</table>
</blockquote></div></div>
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.<br>
<br>
<pre>-- <br>-----------------------------------------------<br>Erik Marklund, PhD student<br>Dept. of Cell and Molecular Biology, Uppsala University.<br>Husargatan 3, Box 596, 75124 Uppsala, Sweden<br>phone: <a rel="nofollow">+46 18 471 4537</a> fax: <a rel="nofollow">+46 18 511 755</a>
<a rel="nofollow" href="mailto:erikm@xray.bmc.uu.se" target="_blank">erikm@xray.bmc.uu.se</a> <a rel="nofollow" href="http://folding.bmc.uu.se/" target="_blank">http://folding.bmc.uu.se/</a>
</pre>
</div>
<br>--<br>
gmx-users mailing list <a rel="nofollow" href="mailto:gmx-users@gromacs.org" target="_blank">gmx-users@gromacs.org</a><br>
<a rel="nofollow" href="http://lists.gromacs.org/mailman/listinfo/gmx-users" target="_blank">http://lists.gromacs.org/mailman/listinfo/gmx-users</a><br>
Please search the archive at <a rel="nofollow" href="http://www.gromacs.org/Support/Mailing_Lists/Search" target="_blank">http://www.gromacs.org/Support/Mailing_Lists/Search</a> before posting!<br>
Please don't post (un)subscribe requests to the list. Use the<br>
www interface or send it to <a rel="nofollow" href="mailto:gmx-users-request@gromacs.org" target="_blank">gmx-users-request@gromacs.org</a>.<br>
Can't post? Read <a rel="nofollow" href="http://www.gromacs.org/Support/Mailing_Lists" target="_blank">http://www.gromacs.org/Support/Mailing_Lists</a><br></blockquote></div><br><br clear="all"><br>-- <br>Tsjerk A. Wassenaar, Ph.D.<br>
<br>post-doctoral researcher<br>Molecular Dynamics Group<br>* Groningen Institute for Biomolecular Research and Biotechnology<br>* Zernike Institute for Advanced Materials<br>University of Groningen<br>The Netherlands<br>
</div></div></div></div>
</div></div></blockquote></div><br><br clear="all"><br>-- <br>Tsjerk A. Wassenaar, Ph.D.<br><br>post-doctoral researcher<br>Molecular Dynamics Group<br>* Groningen Institute for Biomolecular Research and Biotechnology<br>
* Zernike Institute for Advanced Materials<br>University of Groningen<br>The Netherlands<br>