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 "[TT]trjorder[tt] orders molecules according to the smallest distance",
84 "to atoms in a reference group",
85 "or on z-coordinate (with option [TT]-z[tt]).",
86 "With distance ordering, it will ask for a group of reference",
87 "atoms and a group of molecules. For each frame of the trajectory",
88 "the selected molecules will be reordered according to the shortest",
89 "distance between atom number [TT]-da[tt] in the molecule and all the",
90 "atoms in the reference group. The center of mass of the molecules can",
91 "be used instead of a reference atom by setting [TT]-da[tt] to 0.",
92 "All atoms in the trajectory are written",
93 "to the output trajectory.[PAR]",
94 "[TT]trjorder[tt] can be useful for e.g. analyzing the n waters closest to a",
96 "In that case the reference group would be the protein and the group",
97 "of molecules would consist of all the water atoms. When an index group",
98 "of the first n waters is made, the ordered trajectory can be used",
99 "with any Gromacs program to analyze the n closest waters.",
101 "If the output file is a [TT].pdb[tt] file, the distance to the reference target",
102 "will be stored in the B-factor field in order to color with e.g. Rasmol.",
104 "With option [TT]-nshell[tt] the number of molecules within a shell",
105 "of radius [TT]-r[tt] around the reference group are printed."
107 static int na
=3,ref_a
=1;
109 static gmx_bool bCOM
=FALSE
,bZ
=FALSE
;
111 { "-na", FALSE
, etINT
, {&na
},
112 "Number of atoms in a molecule" },
113 { "-da", FALSE
, etINT
, {&ref_a
},
114 "Atom used for the distance calculation, 0 is COM" },
115 { "-com", FALSE
, etBOOL
, {&bCOM
},
116 "Use the distance to the center of mass of the reference group" },
117 { "-r", FALSE
, etREAL
, {&rcut
},
118 "Cutoff used for the distance calculation when computing the number of molecules in a shell around e.g. a protein" },
119 { "-z", FALSE
, etBOOL
, {&bZ
},
120 "Order molecules on z-coordinate" }
125 gmx_bool bNShell
,bPDBout
;
128 rvec
*x
,*xsol
,xcom
,dx
;
132 real t
,totmass
,mass
,rcut2
=0,n2
;
133 int natoms
,nwat
,ncut
;
134 char **grpname
,title
[256];
135 int i
,j
,d
,*isize
,isize_ref
=0,isize_sol
;
136 atom_id sa
,sr
,*swi
,**index
,*ind_ref
=NULL
,*ind_sol
;
139 { efTRX
, "-f", NULL
, ffREAD
},
140 { efTPS
, NULL
, NULL
, ffREAD
},
141 { efNDX
, NULL
, NULL
, ffOPTRD
},
142 { efTRO
, "-o", "ordered", ffOPTWR
},
143 { efXVG
, "-nshell", "nshell", ffOPTWR
}
145 #define NFILE asize(fnm)
147 CopyRight(stderr
,argv
[0]);
148 parse_common_args(&argc
,argv
,PCA_CAN_TIME
| PCA_BE_NICE
,
149 NFILE
,fnm
,asize(pa
),pa
,asize(desc
),desc
,0,NULL
,&oenv
);
151 read_tps_conf(ftp2fn(efTPS
,NFILE
,fnm
),title
,&top
,&ePBC
,&x
,NULL
,box
,TRUE
);
154 /* get index groups */
155 printf("Select %sa group of molecules to be ordered:\n",
156 bZ
? "" : "a group of reference atoms and ");
160 get_index(&top
.atoms
,ftp2fn_null(efNDX
,NFILE
,fnm
),bZ
? 1 : 2,
161 isize
,index
,grpname
);
164 isize_ref
= isize
[0];
165 isize_sol
= isize
[1];
169 isize_sol
= isize
[0];
173 natoms
=read_first_x(oenv
,&status
,ftp2fn(efTRX
,NFILE
,fnm
),&t
,&x
,box
);
174 if (natoms
> top
.atoms
.nr
)
175 gmx_fatal(FARGS
,"Number of atoms in the run input file is larger than in the trjactory");
177 for(j
=0; (j
<isize
[i
]); j
++)
178 if (index
[i
][j
] > natoms
)
179 gmx_fatal(FARGS
,"An atom number in group %s is larger than the number of atoms in the trajectory");
181 if ((isize_sol
% na
) != 0)
182 gmx_fatal(FARGS
,"Number of atoms in the molecule group (%d) is not a multiple of na (%d)",
187 gmx_fatal(FARGS
,"The reference atom can not be larger than the number of atoms in a molecule");
192 for(i
=0; (i
<natoms
); i
++)
197 bNShell
= ((opt2bSet("-nshell",NFILE
,fnm
)) ||
198 (opt2parg_bSet("-r",asize(pa
),pa
)));
202 fp
= xvgropen(opt2fn("-nshell",NFILE
,fnm
),"Number of molecules",
203 "Time (ps)","N",oenv
);
204 printf("Will compute the number of molecules within a radius of %g\n",
207 if (!bNShell
|| opt2bSet("-o",NFILE
,fnm
)) {
208 bPDBout
= (fn2ftp(opt2fn("-o",NFILE
,fnm
)) == efPDB
);
209 if (bPDBout
&& !top
.atoms
.pdbinfo
) {
210 fprintf(stderr
,"Creating pdbfino records\n");
211 snew(top
.atoms
.pdbinfo
,top
.atoms
.nr
);
213 out
= open_trx(opt2fn("-o",NFILE
,fnm
),"w");
215 gpbc
= gmx_rmpbc_init(&top
.idef
,ePBC
,natoms
,box
);
217 gmx_rmpbc(gpbc
,natoms
,box
,x
);
218 set_pbc(&pbc
,ePBC
,box
);
221 /* Calculate the COM of all solvent molecules */
222 for(i
=0; i
<nwat
; i
++) {
225 for(j
=0; j
<na
; j
++) {
226 sa
= ind_sol
[i
*na
+j
];
227 mass
= top
.atoms
.atom
[sa
].m
;
229 for(d
=0; d
<DIM
; d
++) {
230 xsol
[i
][d
] += mass
*x
[sa
][d
];
233 svmul(1/totmass
,xsol
[i
],xsol
[i
]);
236 /* Copy the reference atom of all solvent molecules */
237 for(i
=0; i
<nwat
; i
++) {
238 copy_rvec(x
[ind_sol
[i
*na
+ref_a
]],xsol
[i
]);
243 for(i
=0; (i
<nwat
); i
++) {
246 order
[i
].d2
= xsol
[i
][ZZ
];
251 for(i
=0; i
<isize_ref
; i
++) {
252 mass
= top
.atoms
.atom
[ind_ref
[i
]].m
;
255 xcom
[j
] += mass
*x
[ind_ref
[i
]][j
];
257 svmul(1/totmass
,xcom
,xcom
);
258 for(i
=0; (i
<nwat
); i
++) {
260 pbc_dx(&pbc
,xcom
,xsol
[i
],dx
);
262 order
[i
].d2
= norm2(dx
);
265 /* Set distance to first atom */
266 for(i
=0; (i
<nwat
); i
++) {
268 pbc_dx(&pbc
,x
[ind_ref
[0]],xsol
[i
],dx
);
270 order
[i
].d2
= norm2(dx
);
272 for(j
=1; (j
<isize_ref
); j
++) {
274 for(i
=0; (i
<nwat
); i
++) {
276 pbc_dx(&pbc
,x
[sr
],xsol
[i
],dx
);
278 if (n2
< order
[i
].d2
)
286 for(i
=0; (i
<nwat
); i
++)
287 if (order
[i
].d2
<= rcut2
)
289 fprintf(fp
,"%10.3f %8d\n",t
,ncut
);
292 qsort(order
,nwat
,sizeof(*order
),ocomp
);
293 for(i
=0; (i
<nwat
); i
++)
294 for(j
=0; (j
<na
); j
++)
295 swi
[ind_sol
[na
*i
]+j
] = order
[i
].i
+j
;
297 /* Store the distance as the B-factor */
299 for(i
=0; (i
<nwat
); i
++) {
300 for(j
=0; (j
<na
); j
++) {
301 top
.atoms
.pdbinfo
[order
[i
].i
+j
].bfac
= sqrt(order
[i
].d2
);
305 write_trx(out
,natoms
,swi
,&top
.atoms
,0,t
,box
,x
,NULL
,NULL
);
307 } while(read_next_x(oenv
,status
,&t
,natoms
,x
,box
));
313 gmx_rmpbc_done(gpbc
);