?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Fixed make_edi.c
[gromacs/rigid-bodies.git] / src / tools / gmx_sorient.c
blob9b97c56a13402c12e6aeb9ffeaaf21631391adc2
1 /*
2 *
3 * This source code is part of
4 *
5 * G R O M A C S
6 *
7 * GROningen MAchine for Chemical Simulations
8 *
9 * VERSION 3.2.0
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.
14
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.
19 *
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.
26 *
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.
29 *
30 * For more info, check our website at http://www.gromacs.org
31 *
32 * And Hey:
33 * Green Red Orange Magenta Azure Cyan Skyblue
34 */
35 #ifdef HAVE_CONFIG_H
36 #include <config.h>
37 #endif
38
39 #include "macros.h"
40 #include "statutil.h"
41 #include "smalloc.h"
42 #include "copyrite.h"
43 #include "gstat.h"
44 #include "vec.h"
45 #include "xvgr.h"
46 #include "pbc.h"
47 #include "index.h"
48 #include "tpxio.h"
49 #include "gmx_ana.h"
50
51
52 static void calc_com_pbc(int nrefat,t_topology *top,rvec x[],t_pbc *pbc,
53 atom_id index[],rvec xref,bool bPBC,matrix box)
54 {
55 const real tol=1e-4;
56 bool bChanged;
57 int m,j,ai,iter;
58 real mass,mtot;
59 rvec dx,xtest;
60
61 /* First simple calculation */
62 clear_rvec(xref);
63 mtot = 0;
64 for(m=0; (m<nrefat); m++) {
65 ai = index[m];
66 mass = top->atoms.atom[ai].m;
67 for(j=0; (j<DIM); j++)
68 xref[j] += mass*x[ai][j];
69 mtot += mass;
70 }
71 svmul(1/mtot,xref,xref);
72 /* Now check if any atom is more than half the box from the COM */
73 if (bPBC) {
74 iter = 0;
75 do {
76 bChanged = FALSE;
77 for(m=0; (m<nrefat); m++) {
78 ai = index[m];
79 mass = top->atoms.atom[ai].m/mtot;
80 pbc_dx(pbc,x[ai],xref,dx);
81 rvec_add(xref,dx,xtest);
82 for(j=0; (j<DIM); j++)
83 if (fabs(xtest[j]-x[ai][j]) > tol) {
84 /* Here we have used the wrong image for contributing to the COM */
85 xref[j] += mass*(xtest[j]-x[ai][j]);
86 x[ai][j] = xtest[j];
87 bChanged = TRUE;
88 }
89 }
90 if (bChanged)
91 printf("COM: %8.3f %8.3f %8.3f iter = %d\n",xref[XX],xref[YY],xref[ZZ],iter);
92 iter++;
93 } while (bChanged);
94 }
95 }
96
97 int gmx_sorient(int argc,char *argv[])
98 {
99 t_topology top;
100 int ePBC;
101 char title[STRLEN];
102 int status;
103 int natoms;
104 real t;
105 rvec *xtop,*x;
106 matrix box;
107
108 FILE *fp;
109 int i,j,p,sa0,sa1,sa2,n,ntot,nf,m,*hist1,*hist2,*histn,nbin1,nbin2,nrbin;
110 real *histi1,*histi2,invbw,invrbw;
111 double sum1,sum2;
112 int *isize,nrefgrp,nrefat;
113 atom_id **index;
114 char **grpname;
115 real inp,outp,two_pi,nav,normfac,rmin2,rmax2,rcut,rcut2,r2,r,mass,mtot;
116 real c1,c2;
117 char str[STRLEN];
118 bool bTPS;
119 rvec xref,dx,dxh1,dxh2,outer;
120 t_pbc pbc;
121 char *legr[] = { "<cos(\\8q\\4\\s1\\N)>",
122 "<3cos\\S2\\N(\\8q\\4\\s2\\N)-1>" };
123 char *legc[] = { "cos(\\8q\\4\\s1\\N)",
124 "3cos\\S2\\N(\\8q\\4\\s2\\N)-1" };
125
126 const char *desc[] = {
127 "g_sorient analyzes solvent orientation around solutes.",
128 "It calculates two angles between the vector from one or more",
129 "reference positions to the first atom of each solvent molecule:[BR]",
130 "theta1: the angle with the vector from the first atom of the solvent",
131 "molecule to the midpoint between atoms 2 and 3.[BR]",
132 "theta2: the angle with the normal of the solvent plane, defined by the",
133 "same three atoms, or when the option [TT]-v23[tt] is set",
134 "the angle with the vector between atoms 2 and 3.[BR]",
135 "The reference can be a set of atoms or",
136 "the center of mass of a set of atoms. The group of solvent atoms should",
137 "consist of 3 atoms per solvent molecule.",
138 "Only solvent molecules between [TT]-rmin[tt] and [TT]-rmax[tt] are",
139 "considered for [TT]-o[tt] and [TT]-no[tt] each frame.[PAR]",
140 "[TT]-o[tt]: distribtion of cos(theta1) for rmin<=r<=rmax.[PAR]",
141 "[TT]-no[tt]: distribution of cos(theta2) for rmin<=r<=rmax.[PAR]",
142 "[TT]-ro[tt]: <cos(theta1)> and <3cos^2(theta2)-1> as a function of the",
143 "distance.[PAR]",
144 "[TT]-co[tt]: the sum over all solvent molecules within distance r",
145 "of cos(theta1) and 3cos^2(theta2)-1 as a function of r.[PAR]",
146 "[TT]-rc[tt]: the distribution of the solvent molecules as a function of r"
147 };
148
149 output_env_t oenv;
150 static bool bCom = FALSE,bVec23=FALSE,bPBC = FALSE;
151 static real rmin=0.0,rmax=0.5,binwidth=0.02,rbinw=0.02;
152 t_pargs pa[] = {
153 { "-com", FALSE, etBOOL, {&bCom},
154 "Use the center of mass as the reference postion" },
155 { "-v23", FALSE, etBOOL, {&bVec23},
156 "Use the vector between atoms 2 and 3" },
157 { "-rmin", FALSE, etREAL, {&rmin}, "Minimum distance (nm)" },
158 { "-rmax", FALSE, etREAL, {&rmax}, "Maximum distance (nm)" },
159 { "-cbin", FALSE, etREAL, {&binwidth}, "Binwidth for the cosine" },
160 { "-rbin", FALSE, etREAL, {&rbinw}, "Binwidth for r (nm)" },
161 { "-pbc", FALSE, etBOOL, {&bPBC}, "Check PBC for the center of mass calculation. Only necessary when your reference group consists of several molecules." }
162 };
163
164 t_filenm fnm[] = {
165 { efTRX, NULL, NULL, ffREAD },
166 { efTPS, NULL, NULL, ffREAD },
167 { efNDX, NULL, NULL, ffOPTRD },
168 { efXVG, NULL, "sori.xvg", ffWRITE },
169 { efXVG, "-no", "snor.xvg", ffWRITE },
170 { efXVG, "-ro", "sord.xvg", ffWRITE },
171 { efXVG, "-co", "scum.xvg", ffWRITE },
172 { efXVG, "-rc", "scount.xvg", ffWRITE }
173 };
174 #define NFILE asize(fnm)
175
176 CopyRight(stderr,argv[0]);
177 parse_common_args(&argc,argv,PCA_CAN_TIME | PCA_CAN_VIEW | PCA_BE_NICE,
178 NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL,&oenv);
179
180 two_pi = 2/M_PI;
181
182 bTPS = (opt2bSet("-s",NFILE,fnm) || !opt2bSet("-n",NFILE,fnm) || bCom);
183 if (bTPS) {
184 read_tps_conf(ftp2fn(efTPS,NFILE,fnm),title,&top,&ePBC,&xtop,NULL,box,
185 bCom);
186 }
187
188 /* get index groups */
189 printf("Select a group of reference particles and a solvent group:\n");
190 snew(grpname,2);
191 snew(index,2);
192 snew(isize,2);
193 if (bTPS) {
194 get_index(&top.atoms,ftp2fn_null(efNDX,NFILE,fnm),2,isize,index,grpname);
195 } else {
196 get_index(NULL,ftp2fn(efNDX,NFILE,fnm),2,isize,index,grpname);
197 }
198
199 if (bCom) {
200 nrefgrp = 1;
201 nrefat = isize[0];
202 } else {
203 nrefgrp = isize[0];
204 nrefat = 1;
205 }
206
207 if (isize[1] % 3)
208 gmx_fatal(FARGS,"The number of solvent atoms (%d) is not a multiple of 3",
209 isize[1]);
210
211 /* initialize reading trajectory: */
212 natoms=read_first_x(oenv,&status,ftp2fn(efTRX,NFILE,fnm),&t,&x,box);
213
214 rmin2 = sqr(rmin);
215 rmax2 = sqr(rmax);
216 rcut = 0.99*sqrt(max_cutoff2(guess_ePBC(box),box));
217 if (rcut == 0)
218 rcut = 10*rmax;
219 rcut2 = sqr(rcut);
220
221 invbw = 1/binwidth;
222 nbin1 = (int)(2*invbw + 0.5);
223 nbin2 = (int)(invbw + 0.5);
224
225 invrbw = 1/rbinw;
226
227 snew(hist1,nbin1+1);
228 snew(hist2,nbin2+1);
229 nrbin = rcut/rbinw;
230 if (nrbin == 0)
231 nrbin = 1;
232 snew(histi1,nrbin);
233 snew(histi2,nrbin);
234 snew(histn,nrbin);
235
236 ntot = 0;
237 nf = 0;
238 sum1 = 0;
239 sum2 = 0;
240
241 /* start analysis of trajectory */
242 do {
243 if (bTPS) {
244 /* make molecules whole again */
245 rm_pbc(&top.idef,ePBC,natoms,box,x,x);
246 }
247
248 set_pbc(&pbc,ePBC,box);
249 n = 0;
250 inp = 0;
251 outp = 0;
252 for(p=0; (p<nrefgrp); p++) {
253 if (bCom)
254 calc_com_pbc(nrefat,&top,x,&pbc,index[0],xref,bPBC,box);
255 else
256 copy_rvec(x[index[0][p]],xref);
257
258 for(m=0; m<isize[1]; m+=3) {
259 sa0 = index[1][m];
260 sa1 = index[1][m+1];
261 sa2 = index[1][m+2];
262 pbc_dx(&pbc,x[sa0],xref,dx);
263 r2 = norm2(dx);
264 if (r2 < rcut2) {
265 r = sqrt(r2);
266 if (!bVec23) {
267 /* Determine the normal to the plain */
268 rvec_sub(x[sa1],x[sa0],dxh1);
269 rvec_sub(x[sa2],x[sa0],dxh2);
270 rvec_inc(dxh1,dxh2);
271 svmul(1/r,dx,dx);
272 unitv(dxh1,dxh1);
273 inp = iprod(dx,dxh1);
274 cprod(dxh1,dxh2,outer);
275 unitv(outer,outer);
276 outp = iprod(dx,outer);
277 } else {
278 /* Use the vector between the 2nd and 3rd atom */
279 rvec_sub(x[sa2],x[sa1],dxh2);
280 unitv(dxh2,dxh2);
281 outp = iprod(dx,dxh2)/r;
282 }
283 (histi1[(int)(invrbw*r)]) += inp;
284 (histi2[(int)(invrbw*r)]) += 3*sqr(outp) - 1;
285 (histn[(int)(invrbw*r)])++;
286 if (r2>=rmin2 && r2<rmax2) {
287 (hist1[(int)(invbw*(inp + 1))])++;
288 (hist2[(int)(invbw*fabs(outp))])++;
289 sum1 += inp;
290 sum2 += outp;
291 n++;
292 }
293 }
294 }
295 }
296 ntot += n;
297 nf++;
298
299 } while (read_next_x(oenv,status,&t,natoms,x,box));
300
301 /* clean up */
302 sfree(x);
303 close_trj(status);
304
305 /* Add the bin for the exact maximum to the previous bin */
306 hist1[nbin1-1] += hist1[nbin1];
307 hist2[nbin2-1] += hist2[nbin2];
308
309 nav = (real)ntot/(nrefgrp*nf);
310 normfac = invbw/ntot;
311
312 fprintf(stderr, "Average nr of molecules between %g and %g nm: %.1f\n",
313 rmin,rmax,nav);
314 if (ntot > 0) {
315 sum1 /= ntot;
316 sum2 /= ntot;
317 fprintf(stderr,"Average cos(theta1) between %g and %g nm: %6.3f\n",
318 rmin,rmax,sum1);
319 fprintf(stderr,"Average 3cos2(theta2)-1 between %g and %g nm: %6.3f\n",
320 rmin,rmax,sum2);
321 }
322
323 sprintf(str,"Solvent orientation between %g and %g nm",rmin,rmax);
324 fp=xvgropen(opt2fn("-o",NFILE,fnm), str,"cos(\\8q\\4\\s1\\N)","",oenv);
325 if (output_env_get_print_xvgr_codes(oenv))
326 fprintf(fp,"@ subtitle \"average shell size %.1f molecules\"\n",nav);
327 for(i=0; i<nbin1; i++) {
328 fprintf(fp,"%g %g\n",(i+0.5)*binwidth-1,2*normfac*hist1[i]);
329 }
330 ffclose(fp);
331
332 sprintf(str,"Solvent normal orientation between %g and %g nm",rmin,rmax);
333 fp=xvgropen(opt2fn("-no",NFILE,fnm), str,"cos(\\8q\\4\\s2\\N)","",oenv);
334 if (output_env_get_print_xvgr_codes(oenv))
335 fprintf(fp,"@ subtitle \"average shell size %.1f molecules\"\n",nav);
336 for(i=0; i<nbin2; i++) {
337 fprintf(fp,"%g %g\n",(i+0.5)*binwidth,normfac*hist2[i]);
338 }
339 ffclose(fp);
340
341
342 sprintf(str,"Solvent orientation");
343 fp=xvgropen(opt2fn("-ro",NFILE,fnm),str,"r (nm)","",oenv);
344 if (output_env_get_print_xvgr_codes(oenv))
345 fprintf(fp,"@ subtitle \"as a function of distance\"\n");
346 xvgr_legend(fp,2,legr,oenv);
347 for(i=0; i<nrbin; i++)
348 fprintf(fp,"%g %g %g\n",(i+0.5)*rbinw,
349 histn[i] ? histi1[i]/histn[i] : 0,
350 histn[i] ? histi2[i]/histn[i] : 0);
351 ffclose(fp);
352
353 sprintf(str,"Cumulative solvent orientation");
354 fp=xvgropen(opt2fn("-co",NFILE,fnm),str,"r (nm)","",oenv);
355 if (output_env_get_print_xvgr_codes(oenv))
356 fprintf(fp,"@ subtitle \"as a function of distance\"\n");
357 xvgr_legend(fp,2,legc,oenv);
358 normfac = 1.0/(nrefgrp*nf);
359 c1 = 0;
360 c2 = 0;
361 fprintf(fp,"%g %g %g\n",0.0,c1,c2);
362 for(i=0; i<nrbin; i++) {
363 c1 += histi1[i]*normfac;
364 c2 += histi2[i]*normfac;
365 fprintf(fp,"%g %g %g\n",(i+1)*rbinw,c1,c2);
366 }
367 ffclose(fp);
368
369 sprintf(str,"Solvent distribution");
370 fp=xvgropen(opt2fn("-rc",NFILE,fnm),str,"r (nm)","molecules/nm",oenv);
371 if (output_env_get_print_xvgr_codes(oenv))
372 fprintf(fp,"@ subtitle \"as a function of distance\"\n");
373 normfac = 1.0/(rbinw*nf);
374 for(i=0; i<nrbin; i++) {
375 fprintf(fp,"%g %g\n",(i+0.5)*rbinw,histn[i]*normfac);
376 }
377 ffclose(fp);
378
379 do_view(oenv, opt2fn("-o",NFILE,fnm),NULL);
380 do_view(oenv, opt2fn("-no",NFILE,fnm),NULL);
381 do_view(oenv, opt2fn("-ro",NFILE,fnm),"-nxy");
382 do_view(oenv, opt2fn("-co",NFILE,fnm),"-nxy");
383
384 thanx(stderr);
385
386 return 0;
387 }