3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
10 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
11 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
12 * Copyright (c) 2001-2004, The GROMACS development team,
13 * check out http://www.gromacs.org for more information.
15 * This program is free software; you can redistribute it and/or
16 * modify it under the terms of the GNU General Public License
17 * as published by the Free Software Foundation; either version 2
18 * of the License, or (at your option) any later version.
20 * If you want to redistribute modifications, please consider that
21 * scientific software is very special. Version control is crucial -
22 * bugs must be traceable. We will be happy to consider code for
23 * inclusion in the official distribution, but derived work must not
24 * be called official GROMACS. Details are found in the README & COPYING
25 * files - if they are missing, get the official version at www.gromacs.org.
27 * To help us fund GROMACS development, we humbly ask that you cite
28 * the papers on the package - you can find them in the top README file.
30 * For more info, check our website at http://www.gromacs.org
33 * Green Red Orange Magenta Azure Cyan Skyblue
67 static int ocomp(const void *a
,const void *b
)
80 int gmx_trjorder(int argc
,char *argv
[])
82 const char *desc
[] = {
83 "trjorder orders molecules according to the smallest distance",
84 "to atoms in a reference group. It will ask for a group of reference",
85 "atoms and a group of molecules. For each frame of the trajectory",
86 "the selected molecules will be reordered according to the shortest",
87 "distance between atom number [TT]-da[tt] in the molecule and all the",
88 "atoms in the reference group. All atoms in the trajectory are written",
89 "to the output trajectory.[PAR]",
90 "trjorder can be useful for e.g. analyzing the n waters closest to a",
92 "In that case the reference group would be the protein and the group",
93 "of molecules would consist of all the water atoms. When an index group",
94 "of the first n waters is made, the ordered trajectory can be used",
95 "with any Gromacs program to analyze the n closest waters.",
97 "If the output file is a pdb file, the distance to the reference target",
98 "will be stored in the B-factor field in order to color with e.g. rasmol.",
100 "With option [TT]-nshell[tt] the number of molecules within a shell",
101 "of radius [TT]-r[tt] around the refernce group are printed."
103 static int na
=3,ref_a
=1;
105 static bool bCOM
= FALSE
;
107 { "-na", FALSE
, etINT
, {&na
},
108 "Number of atoms in a molecule" },
109 { "-da", FALSE
, etINT
, {&ref_a
},
110 "Atom used for the distance calculation" },
111 { "-com", FALSE
, etBOOL
, {&bCOM
},
112 "Use the distance to the center of mass of the reference group" },
113 { "-r", FALSE
, etREAL
, {&rcut
},
114 "Cutoff used for the distance calculation when computing the number of molecules in a shell around e.g. a protein" },
118 bool bNShell
,bPDBout
;
124 real t
,totmass
,mass
,rcut2
=0,n2
;
125 int natoms
,nwat
,ncut
;
126 char **grpname
,title
[256];
128 atom_id sa
,sr
,*swi
,**index
;
133 { efTRX
, "-f", NULL
, ffREAD
},
134 { efTPS
, NULL
, NULL
, ffREAD
},
135 { efNDX
, NULL
, NULL
, ffOPTRD
},
136 { efTRO
, "-o", "ordered", ffOPTWR
},
137 { efXVG
, "-nshell", "nshell", ffOPTWR
}
139 #define NFILE asize(fnm)
141 CopyRight(stderr
,argv
[0]);
142 parse_common_args(&argc
,argv
,PCA_CAN_TIME
| PCA_BE_NICE
,
143 NFILE
,fnm
,asize(pa
),pa
,asize(desc
),desc
,0,NULL
,&oenv
);
145 read_tps_conf(ftp2fn(efTPS
,NFILE
,fnm
),title
,&top
,&ePBC
,&x
,NULL
,box
,TRUE
);
148 /* get index groups */
149 printf("Select a group of reference atoms and a group of molecules to be ordered:\n");
153 get_index(&top
.atoms
,ftp2fn_null(efNDX
,NFILE
,fnm
),2,isize
,index
,grpname
);
155 natoms
=read_first_x(oenv
,&status
,ftp2fn(efTRX
,NFILE
,fnm
),&t
,&x
,box
);
156 if (natoms
> top
.atoms
.nr
)
157 gmx_fatal(FARGS
,"Number of atoms in the run input file is larger than in the trjactory");
159 for(j
=0; (j
<isize
[i
]); j
++)
160 if (index
[i
][j
] > natoms
)
161 gmx_fatal(FARGS
,"An atom number in group %s is larger than the number of atoms in the trajectory");
163 if ((isize
[SOL
] % na
) != 0)
164 gmx_fatal(FARGS
,"Number of atoms in the molecule group (%d) is not a multiple of na (%d)",
167 nwat
= isize
[SOL
]/na
;
169 gmx_fatal(FARGS
,"The reference atom can not be larger than the number of atoms in a molecule");
173 for(i
=0; (i
<natoms
); i
++)
178 bNShell
= ((opt2bSet("-nshell",NFILE
,fnm
)) ||
179 (opt2parg_bSet("-r",asize(pa
),pa
)));
183 fp
= xvgropen(opt2fn("-nshell",NFILE
,fnm
),"Number of molecules",
184 "Time (ps)","N",oenv
);
185 printf("Will compute the number of molecules within a radius of %g\n",
188 if (!bNShell
|| opt2bSet("-o",NFILE
,fnm
)) {
189 bPDBout
= (fn2ftp(opt2fn("-o",NFILE
,fnm
)) == efPDB
);
190 if (bPDBout
&& !top
.atoms
.pdbinfo
) {
191 fprintf(stderr
,"Creating pdbfino records\n");
192 snew(top
.atoms
.pdbinfo
,top
.atoms
.nr
);
194 out
= open_trx(opt2fn("-o",NFILE
,fnm
),"w");
197 rm_pbc(&top
.idef
,ePBC
,natoms
,box
,x
,x
);
198 set_pbc(&pbc
,ePBC
,box
);
203 for(i
=0; i
<isize
[REF
]; i
++) {
204 mass
= top
.atoms
.atom
[index
[REF
][i
]].m
;
207 xcom
[j
] += mass
*x
[index
[REF
][i
]][j
];
209 svmul(1/totmass
,xcom
,xcom
);
210 for(i
=0; (i
<nwat
); i
++) {
211 sa
= index
[SOL
][na
*i
];
212 pbc_dx(&pbc
,xcom
,x
[sa
+ref_a
],dx
);
214 order
[i
].d2
= norm2(dx
);
217 /* Set distance to first atom */
218 for(i
=0; (i
<nwat
); i
++) {
219 sa
= index
[SOL
][na
*i
];
220 pbc_dx(&pbc
,x
[index
[REF
][0]],x
[sa
+ref_a
],dx
);
222 order
[i
].d2
= norm2(dx
);
224 for(j
=1; (j
<isize
[REF
]); j
++) {
226 for(i
=0; (i
<nwat
); i
++) {
227 sa
= index
[SOL
][na
*i
];
228 pbc_dx(&pbc
,x
[sr
],x
[sa
+ref_a
],dx
);
230 if (n2
< order
[i
].d2
)
238 for(i
=0; (i
<nwat
); i
++)
239 if (order
[i
].d2
<= rcut2
)
241 fprintf(fp
,"%10.3f %8d\n",t
,ncut
);
244 qsort(order
,nwat
,sizeof(*order
),ocomp
);
245 for(i
=0; (i
<nwat
); i
++)
246 for(j
=0; (j
<na
); j
++)
247 swi
[index
[SOL
][na
*i
]+j
] = order
[i
].i
+j
;
249 /* Store the distance as the B-factor */
251 for(i
=0; (i
<nwat
); i
++) {
252 for(j
=0; (j
<na
); j
++) {
253 top
.atoms
.pdbinfo
[order
[i
].i
+j
].bfac
= sqrt(order
[i
].d2
);
257 write_trx(out
,natoms
,swi
,&top
.atoms
,0,t
,box
,x
,NULL
,NULL
);
259 } while(read_next_x(oenv
,status
,&t
,natoms
,x
,box
));