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
52 static void calc_com_pbc(int nrefat,t_topology *top,rvec x[],t_pbc *pbc,
53 atom_id index[],rvec xref,gmx_bool bPBC,matrix box)
61 /* First simple calculation */
64 for(m=0; (m<nrefat); m++) {
66 mass = top->atoms.atom[ai].m;
67 for(j=0; (j<DIM); j++)
68 xref[j] += mass*x[ai][j];
71 svmul(1/mtot,xref,xref);
72 /* Now check if any atom is more than half the box from the COM */
77 for(m=0; (m<nrefat); 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]);
91 printf("COM: %8.3f %8.3f %8.3f iter = %d\n",xref[XX],xref[YY],xref[ZZ],iter);
97 int gmx_sorient(int argc,char *argv[])
109 int i,j,p,sa0,sa1,sa2,n,ntot,nf,m,*hist1,*hist2,*histn,nbin1,nbin2,nrbin;
110 real *histi1,*histi2,invbw,invrbw;
112 int *isize,nrefgrp,nrefat;
115 real inp,outp,two_pi,nav,normfac,rmin2,rmax2,rcut,rcut2,r2,r,mass,mtot;
119 rvec xref,dx,dxh1,dxh2,outer;
120 gmx_rmpbc_t gpbc=NULL;
122 const char *legr[] = { "<cos(\\8q\\4\\s1\\N)>",
123 "<3cos\\S2\\N(\\8q\\4\\s2\\N)-1>" };
124 const char *legc[] = { "cos(\\8q\\4\\s1\\N)",
125 "3cos\\S2\\N(\\8q\\4\\s2\\N)-1" };
127 const char *desc[] = {
128 "[TT]g_sorient[tt] analyzes solvent orientation around solutes.",
129 "It calculates two angles between the vector from one or more",
130 "reference positions to the first atom of each solvent molecule:[BR]",
131 "theta1: the angle with the vector from the first atom of the solvent",
132 "molecule to the midpoint between atoms 2 and 3.[BR]",
133 "theta2: the angle with the normal of the solvent plane, defined by the",
134 "same three atoms, or when the option [TT]-v23[tt] is set",
135 "the angle with the vector between atoms 2 and 3.[BR]",
136 "The reference can be a set of atoms or",
137 "the center of mass of a set of atoms. The group of solvent atoms should",
138 "consist of 3 atoms per solvent molecule.",
139 "Only solvent molecules between [TT]-rmin[tt] and [TT]-rmax[tt] are",
140 "considered for [TT]-o[tt] and [TT]-no[tt] each frame.[PAR]",
141 "[TT]-o[tt]: distribtion of cos(theta1) for rmin<=r<=rmax.[PAR]",
142 "[TT]-no[tt]: distribution of cos(theta2) for rmin<=r<=rmax.[PAR]",
143 "[TT]-ro[tt]: <cos(theta1)> and <3cos^2(theta2)-1> as a function of the",
145 "[TT]-co[tt]: the sum over all solvent molecules within distance r",
146 "of cos(theta1) and 3cos^2(theta2)-1 as a function of r.[PAR]",
147 "[TT]-rc[tt]: the distribution of the solvent molecules as a function of r"
151 static gmx_bool bCom = FALSE,bVec23=FALSE,bPBC = FALSE;
152 static real rmin=0.0,rmax=0.5,binwidth=0.02,rbinw=0.02;
154 { "-com", FALSE, etBOOL, {&bCom},
155 "Use the center of mass as the reference postion" },
156 { "-v23", FALSE, etBOOL, {&bVec23},
157 "Use the vector between atoms 2 and 3" },
158 { "-rmin", FALSE, etREAL, {&rmin}, "Minimum distance (nm)" },
159 { "-rmax", FALSE, etREAL, {&rmax}, "Maximum distance (nm)" },
160 { "-cbin", FALSE, etREAL, {&binwidth}, "Binwidth for the cosine" },
161 { "-rbin", FALSE, etREAL, {&rbinw}, "Binwidth for r (nm)" },
162 { "-pbc", FALSE, etBOOL, {&bPBC}, "Check PBC for the center of mass calculation. Only necessary when your reference group consists of several molecules." }
166 { efTRX, NULL, NULL, ffREAD },
167 { efTPS, NULL, NULL, ffREAD },
168 { efNDX, NULL, NULL, ffOPTRD },
169 { efXVG, NULL, "sori.xvg", ffWRITE },
170 { efXVG, "-no", "snor.xvg", ffWRITE },
171 { efXVG, "-ro", "sord.xvg", ffWRITE },
172 { efXVG, "-co", "scum.xvg", ffWRITE },
173 { efXVG, "-rc", "scount.xvg", ffWRITE }
175 #define NFILE asize(fnm)
177 CopyRight(stderr,argv[0]);
178 parse_common_args(&argc,argv,PCA_CAN_TIME | PCA_CAN_VIEW | PCA_BE_NICE,
179 NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL,&oenv);
183 bTPS = (opt2bSet("-s",NFILE,fnm) || !opt2bSet("-n",NFILE,fnm) || bCom);
185 read_tps_conf(ftp2fn(efTPS,NFILE,fnm),title,&top,&ePBC,&xtop,NULL,box,
189 /* get index groups */
190 printf("Select a group of reference particles and a solvent group:\n");
195 get_index(&top.atoms,ftp2fn_null(efNDX,NFILE,fnm),2,isize,index,grpname);
197 get_index(NULL,ftp2fn(efNDX,NFILE,fnm),2,isize,index,grpname);
209 gmx_fatal(FARGS,"The number of solvent atoms (%d) is not a multiple of 3",
212 /* initialize reading trajectory: */
213 natoms=read_first_x(oenv,&status,ftp2fn(efTRX,NFILE,fnm),&t,&x,box);
217 rcut = 0.99*sqrt(max_cutoff2(guess_ePBC(box),box));
223 nbin1 = 1+(int)(2*invbw + 0.5);
224 nbin2 = 1+(int)(invbw + 0.5);
230 nrbin = 1+(int)(rcut/rbinw);
243 /* make molecules whole again */
244 gpbc = gmx_rmpbc_init(&top.idef,ePBC,natoms,box);
246 /* start analysis of trajectory */
249 /* make molecules whole again */
250 gmx_rmpbc(gpbc,natoms,box,x);
253 set_pbc(&pbc,ePBC,box);
257 for(p=0; (p<nrefgrp); p++) {
259 calc_com_pbc(nrefat,&top,x,&pbc,index[0],xref,bPBC,box);
261 copy_rvec(x[index[0][p]],xref);
263 for(m=0; m<isize[1]; m+=3) {
267 range_check(sa0,0,natoms);
268 range_check(sa1,0,natoms);
269 range_check(sa2,0,natoms);
270 pbc_dx(&pbc,x[sa0],xref,dx);
275 /* Determine the normal to the plain */
276 rvec_sub(x[sa1],x[sa0],dxh1);
277 rvec_sub(x[sa2],x[sa0],dxh2);
281 inp = iprod(dx,dxh1);
282 cprod(dxh1,dxh2,outer);
284 outp = iprod(dx,outer);
286 /* Use the vector between the 2nd and 3rd atom */
287 rvec_sub(x[sa2],x[sa1],dxh2);
289 outp = iprod(dx,dxh2)/r;
292 int ii = (int)(invrbw*r);
293 range_check(ii,0,nrbin);
295 histi2[ii] += 3*sqr(outp) - 1;
298 if ((r2>=rmin2) && (r2<rmax2)) {
299 int ii1 = (int)(invbw*(inp + 1));
300 int ii2 = (int)(invbw*fabs(outp));
302 range_check(ii1,0,nbin1);
303 range_check(ii2,0,nbin2);
316 } while (read_next_x(oenv,status,&t,natoms,x,box));
321 gmx_rmpbc_done(gpbc);
323 /* Add the bin for the exact maximum to the previous bin */
324 hist1[nbin1-1] += hist1[nbin1];
325 hist2[nbin2-1] += hist2[nbin2];
327 nav = (real)ntot/(nrefgrp*nf);
328 normfac = invbw/ntot;
330 fprintf(stderr, "Average nr of molecules between %g and %g nm: %.1f\n",
335 fprintf(stderr,"Average cos(theta1) between %g and %g nm: %6.3f\n",
337 fprintf(stderr,"Average 3cos2(theta2)-1 between %g and %g nm: %6.3f\n",
341 sprintf(str,"Solvent orientation between %g and %g nm",rmin,rmax);
342 fp=xvgropen(opt2fn("-o",NFILE,fnm), str,"cos(\\8q\\4\\s1\\N)","",oenv);
343 if (output_env_get_print_xvgr_codes(oenv))
344 fprintf(fp,"@ subtitle \"average shell size %.1f molecules\"\n",nav);
345 for(i=0; i<nbin1; i++) {
346 fprintf(fp,"%g %g\n",(i+0.5)*binwidth-1,2*normfac*hist1[i]);
350 sprintf(str,"Solvent normal orientation between %g and %g nm",rmin,rmax);
351 fp=xvgropen(opt2fn("-no",NFILE,fnm), str,"cos(\\8q\\4\\s2\\N)","",oenv);
352 if (output_env_get_print_xvgr_codes(oenv))
353 fprintf(fp,"@ subtitle \"average shell size %.1f molecules\"\n",nav);
354 for(i=0; i<nbin2; i++) {
355 fprintf(fp,"%g %g\n",(i+0.5)*binwidth,normfac*hist2[i]);
360 sprintf(str,"Solvent orientation");
361 fp=xvgropen(opt2fn("-ro",NFILE,fnm),str,"r (nm)","",oenv);
362 if (output_env_get_print_xvgr_codes(oenv))
363 fprintf(fp,"@ subtitle \"as a function of distance\"\n");
364 xvgr_legend(fp,2,legr,oenv);
365 for(i=0; i<nrbin; i++)
366 fprintf(fp,"%g %g %g\n",(i+0.5)*rbinw,
367 histn[i] ? histi1[i]/histn[i] : 0,
368 histn[i] ? histi2[i]/histn[i] : 0);
371 sprintf(str,"Cumulative solvent orientation");
372 fp=xvgropen(opt2fn("-co",NFILE,fnm),str,"r (nm)","",oenv);
373 if (output_env_get_print_xvgr_codes(oenv))
374 fprintf(fp,"@ subtitle \"as a function of distance\"\n");
375 xvgr_legend(fp,2,legc,oenv);
376 normfac = 1.0/(nrefgrp*nf);
379 fprintf(fp,"%g %g %g\n",0.0,c1,c2);
380 for(i=0; i<nrbin; i++) {
381 c1 += histi1[i]*normfac;
382 c2 += histi2[i]*normfac;
383 fprintf(fp,"%g %g %g\n",(i+1)*rbinw,c1,c2);
387 sprintf(str,"Solvent distribution");
388 fp=xvgropen(opt2fn("-rc",NFILE,fnm),str,"r (nm)","molecules/nm",oenv);
389 if (output_env_get_print_xvgr_codes(oenv))
390 fprintf(fp,"@ subtitle \"as a function of distance\"\n");
391 normfac = 1.0/(rbinw*nf);
392 for(i=0; i<nrbin; i++) {
393 fprintf(fp,"%g %g\n",(i+0.5)*rbinw,histn[i]*normfac);
397 do_view(oenv, opt2fn("-o",NFILE,fnm),NULL);
398 do_view(oenv, opt2fn("-no",NFILE,fnm),NULL);
399 do_view(oenv, opt2fn("-ro",NFILE,fnm),"-nxy");
400 do_view(oenv, opt2fn("-co",NFILE,fnm),"-nxy");