[gmx-users] micelles and trjconv -pbc cluster
Tsjerk Wassenaar
tsjerkw at gmail.com
Fri Apr 15 06:39:56 CEST 2011
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,
> but none which has neighbours in the rest group:
> 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
> */
On Fri, Apr 15, 2011 at 3:16 AM, Jianguo Li <ljggmx at yahoo.com.sg> wrote:
> 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
>
> On Thu, Apr 14, 2011 at 4:48 PM, Erik Marklund <erikm at xray.bmc.uu.se>wrote:
>
>> jim jack skrev 2011-04-14 16.42:
>>
>> Dear 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>
More information about the gromacs.org_gmx-users
mailing list