2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team,
6 * check out http://www.gromacs.org for more information.
7 * Copyright (c) 2012,2013, by the GROMACS development team, led by
8 * David van der Spoel, Berk Hess, Erik Lindahl, and including many
9 * others, as listed in the AUTHORS file in the top-level source
10 * directory and at http://www.gromacs.org.
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
61 /****************************************************************************/
62 /* This program calculates the ordering of water molecules across a box, as */
63 /* function of the z-coordinate. This implies averaging over slices and over*/
64 /* time. Output is the average cos of the angle of the dipole with the */
65 /* normal, per slice. */
66 /* In addition, it calculates the average dipole moment itself in three */
68 /****************************************************************************/
70 void calc_h2order(const char *fn, atom_id index[], int ngx, rvec **slDipole,
71 real **slOrder, real *slWidth, int *nslices,
72 t_topology *top, int ePBC,
73 int axis, gmx_bool bMicel, atom_id micel[], int nmic,
74 const output_env_t oenv)
76 rvec *x0, /* coordinates with pbc */
77 dipole, /* dipole moment due to one molecules */
79 com; /* center of mass of micel, with bMicel */
80 rvec *dip; /* sum of dipoles, unnormalized */
81 matrix box; /* box (3x3) */
83 real t, /* time from trajectory */
84 *sum, /* sum of all cosines of dipoles, per slice */
85 *frame; /* order over one frame */
86 int natoms, /* nr. atoms in trj */
88 slice=0, /* current slice number */
89 *count; /* nr. of atoms in one slice */
90 gmx_rmpbc_t gpbc=NULL;
92 if ((natoms = read_first_x(oenv,&status,fn,&t,&x0,box)) == 0)
93 gmx_fatal(FARGS,"Could not read coordinates from statusfile\n");
96 *nslices = (int)(box[axis][axis] * 10); /* default value */
102 normal[0] = 1; normal[1] = 0; normal[2] = 0;
105 normal[0] = 0; normal[1] = 1; normal[2] = 0;
108 normal[0] = 0; normal[1] = 0; normal[2] = 1;
111 gmx_fatal(FARGS,"No valid value for -axis-. Exiting.\n");
112 /* make compiler happy */
113 normal[0] = 1; normal[1] = 0; normal[2] = 0;
117 snew(count, *nslices);
120 snew(frame, *nslices);
122 *slWidth = box[axis][axis]/(*nslices);
123 fprintf(stderr,"Box divided in %d slices. Initial width of slice: %f\n",
128 gpbc = gmx_rmpbc_init(&top->idef,ePBC,natoms,box);
129 /*********** Start processing trajectory ***********/
132 *slWidth = box[axis][axis]/(*nslices);
135 gmx_rmpbc(gpbc,natoms,box,x0);
138 calc_xcm(x0, nmic, micel, top->atoms.atom, com, FALSE);
140 for (i = 0; i < ngx/3; i++)
142 /* put all waters in box */
143 for (j = 0; j < DIM; j++)
145 if (x0[index[3*i]][j] < 0)
147 x0[index[3*i]][j] += box[j][j];
148 x0[index[3*i+1]][j] += box[j][j];
149 x0[index[3*i+2]][j] += box[j][j];
151 if (x0[index[3*i]][j] > box[j][j])
153 x0[index[3*i]][j] -= box[j][j];
154 x0[index[3*i+1]][j] -= box[j][j];
155 x0[index[3*i+2]][j] -= box[j][j];
159 for (j = 0; j < DIM; j++)
161 x0[index[3*i]][j] * top->atoms.atom[index[3*i]].q +
162 x0[index[3*i+1]][j] * top->atoms.atom[index[3*i+1]].q +
163 x0[index[3*i+2]][j] * top->atoms.atom[index[3*i+2]].q;
165 /* now we have a dipole vector. Might as well safe it. Then the
166 rest depends on whether we're dealing with
167 a flat or a spherical interface.
171 { /* this is for spherical interfaces */
172 rvec_sub(com, x0[index[3*i]], normal); /* vector from Oxygen to COM */
173 slice = norm(normal)/(*slWidth); /* spherical slice */
175 sum[slice] += iprod(dipole, normal) / (norm(dipole) * norm(normal));
176 frame[slice] += iprod(dipole, normal) / (norm(dipole) * norm(normal));
180 /* this is for flat interfaces */
182 /* determine which slice atom is in */
183 slice = (x0[index[3*i]][axis] / (*slWidth));
184 if (slice < 0 || slice >= *nslices)
186 fprintf(stderr,"Coordinate: %f ",x0[index[3*i]][axis]);
187 fprintf(stderr,"HELP PANIC! slice = %d, OUT OF RANGE!\n",slice);
191 rvec_add(dipole, dip[slice], dip[slice]);
192 /* Add dipole to total. mag[slice] is total dipole in axis direction */
193 sum[slice] += iprod(dipole,normal)/norm(dipole);
194 frame[slice] += iprod(dipole,normal)/norm(dipole);
195 /* increase count for that slice */
201 } while (read_next_x(oenv,status,&t,natoms,x0,box));
202 /*********** done with status file **********/
204 fprintf(stderr,"\nRead trajectory. Printing parameters to file\n");
205 gmx_rmpbc_done(gpbc);
207 for (i = 0; i < *nslices; i++) /* average over frames */
209 fprintf(stderr,"%d waters in slice %d\n",count[i],i);
210 if (count[i] > 0) /* divide by number of molecules in each slice */
212 sum[i] = sum[i] / count[i];
213 dip[i][XX] = dip[i][XX] / count[i];
214 dip[i][YY] = dip[i][YY] / count[i];
215 dip[i][ZZ] = dip[i][ZZ] / count[i];
218 fprintf(stderr,"No water in slice %d\n",i);
221 *slOrder = sum; /* copy a pointer, I hope */
223 sfree(x0); /* free memory used by coordinate arrays */
226 void h2order_plot(rvec dipole[], real order[], const char *afile,
227 int nslices, real slWidth, const output_env_t oenv)
229 FILE *ord; /* xvgr files with order parameters */
230 int slice; /* loop index */
231 char buf[256]; /* for xvgr title */
232 real factor; /* conversion to Debye from electron*nm */
234 /* factor = 1e-9*1.60217733e-19/3.336e-30 */
235 factor = 1.60217733/3.336e-2;
236 fprintf(stderr,"%d slices\n",nslices);
237 sprintf(buf,"Water orientation with respect to normal");
238 ord = xvgropen(afile,buf,
239 "box (nm)","mu_x, mu_y, mu_z (D), cosine with normal",oenv);
241 for (slice = 0; slice < nslices; slice++)
242 fprintf(ord,"%8.3f %8.3f %8.3f %8.3f %e\n", slWidth*slice,
243 factor*dipole[slice][XX], factor*dipole[slice][YY],
244 factor*dipole[slice][ZZ], order[slice]);
249 int gmx_h2order(int argc,char *argv[])
251 const char *desc[] = {
252 "[TT]g_h2order[tt] computes the orientation of water molecules with respect to the normal",
253 "of the box. The program determines the average cosine of the angle",
254 "between the dipole moment of water and an axis of the box. The box is",
255 "divided in slices and the average orientation per slice is printed.",
256 "Each water molecule is assigned to a slice, per time frame, based on the",
257 "position of the oxygen. When [TT]-nm[tt] is used, the angle between the water",
258 "dipole and the axis from the center of mass to the oxygen is calculated",
259 "instead of the angle between the dipole and a box axis."
261 static int axis = 2; /* normal to memb. default z */
262 static const char *axtitle="Z";
263 static int nslices = 0; /* nr of slices defined */
265 { "-d", FALSE, etSTR, {&axtitle},
266 "Take the normal on the membrane in direction X, Y or Z." },
267 { "-sl", FALSE, etINT, {&nslices},
268 "Calculate order parameter as function of boxlength, dividing the box"
269 " in this number of slices."}
271 const char *bugs[] = {
272 "The program assigns whole water molecules to a slice, based on the first "
273 "atom of three in the index file group. It assumes an order O,H,H. "
274 "Name is not important, but the order is. If this demand is not met, "
275 "assigning molecules to slices is different."
279 real *slOrder, /* av. cosine, per slice */
280 slWidth = 0.0; /* width of a slice */
282 char *grpname, /* groupnames */
284 int ngx, /* nr. of atomsin sol group */
285 nmic=0; /* nr. of atoms in micelle */
286 t_topology *top; /* topology */
288 atom_id *index, /* indices for solvent group */
290 gmx_bool bMicel = FALSE; /* think we're a micel */
291 t_filenm fnm[] = { /* files for g_order */
292 { efTRX, "-f", NULL, ffREAD }, /* trajectory file */
293 { efNDX, NULL, NULL, ffREAD }, /* index file */
294 { efNDX, "-nm", NULL, ffOPTRD }, /* index with micelle atoms */
295 { efTPX, NULL, NULL, ffREAD }, /* topology file */
296 { efXVG,"-o", "order", ffWRITE }, /* xvgr output file */
299 #define NFILE asize(fnm)
301 CopyRight(stderr,argv[0]);
302 parse_common_args(&argc, argv,
303 PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE, NFILE,
304 fnm, asize(pa),pa,asize(desc),desc,asize(bugs),bugs,&oenv);
305 bMicel = opt2bSet("-nm",NFILE,fnm);
307 top = read_top(ftp2fn(efTPX,NFILE,fnm),&ePBC); /* read topology file */
309 rd_index(ftp2fn(efNDX,NFILE,fnm),1,&ngx,&index,&grpname);
312 rd_index(opt2fn("-nm",NFILE,fnm), 1, &nmic, &micelle, &micname);
314 calc_h2order(ftp2fn(efTRX,NFILE,fnm), index, ngx, &slDipole, &slOrder,
315 &slWidth, &nslices, top, ePBC, axis, bMicel, micelle, nmic,
318 h2order_plot(slDipole, slOrder, opt2fn("-o",NFILE,fnm), nslices,
321 do_view(oenv,opt2fn("-o",NFILE,fnm),"-nxy"); /* view xvgr file */