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:[PAR]",
131 "[GRK]theta[grk][SUB]1[sub]: the angle with the vector from the first atom of the solvent",
132 "molecule to the midpoint between atoms 2 and 3.[BR]",
133 "[GRK]theta[grk][SUB]2[sub]: 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.[PAR]",
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 [MATH][COS][GRK]theta[grk][SUB]1[sub][cos][math] for rmin<=r<=rmax.[PAR]",
142 "[TT]-no[tt]: distribution of [MATH][COS][GRK]theta[grk][SUB]2[sub][cos][math] for rmin<=r<=rmax.[PAR]",
143 "[TT]-ro[tt]: [MATH][CHEVRON][COS][GRK]theta[grk][SUB]1[sub][cos][chevron][math] and [MATH][CHEVRON]3[COS]^2[GRK]theta[grk][SUB]2[sub][cos]-1[chevron][math] as a function of the",
145 "[TT]-co[tt]: the sum over all solvent molecules within distance r",
146 "of [MATH][COS][GRK]theta[grk][SUB]1[sub][cos][math] and [MATH]3[COS]^2([GRK]theta[grk][SUB]2[sub])-1[cos][math] 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");