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, int ePBC)
55 const real tol = 1e-4;
61 /* First simple calculation */
64 for (m = 0; (m < nrefat); m++)
67 mass = top->atoms.atom[ai].m;
68 for (j = 0; (j < DIM); j++)
70 xref[j] += mass*x[ai][j];
74 svmul(1/mtot, xref, xref);
75 /* Now check if any atom is more than half the box from the COM */
82 for (m = 0; (m < nrefat); m++)
85 mass = top->atoms.atom[ai].m/mtot;
86 pbc_dx(pbc, x[ai], xref, dx);
87 rvec_add(xref, dx, xtest);
88 for (j = 0; (j < DIM); j++)
90 if (fabs(xtest[j]-x[ai][j]) > tol)
92 /* Here we have used the wrong image for contributing to the COM */
93 xref[j] += mass*(xtest[j]-x[ai][j]);
101 printf("COM: %8.3f %8.3f %8.3f iter = %d\n", xref[XX], xref[YY], xref[ZZ], iter);
109 void spol_atom2molindex(int *n, int *index, t_block *mols)
118 while (m < mols->nr && index[i] != mols->index[m])
124 gmx_fatal(FARGS, "index[%d]=%d does not correspond to the first atom of a molecule", i+1, index[i]+1);
126 for (j = mols->index[m]; j < mols->index[m+1]; j++)
128 if (i >= *n || index[i] != j)
130 gmx_fatal(FARGS, "The index group is not a set of whole molecules");
134 /* Modify the index in place */
137 printf("There are %d molecules in the selection\n", nmol);
142 int gmx_spol(int argc, char *argv[])
149 int nrefat, natoms, nf, ntot;
151 rvec *xtop, *x, xref, trial, dx = {0}, dip, dir;
156 atom_id **index, *molindex;
158 real rmin2, rmax2, rcut, rcut2, rdx2 = 0, rtry2, qav, q, dip2, invbw;
159 int nbin, i, m, mol, a0, a1, a, d;
160 double sdip, sdip2, sinp, sdinp, nmol;
163 gmx_rmpbc_t gpbc = NULL;
166 const char *desc[] = {
167 "[TT]g_spol[tt] analyzes dipoles around a solute; it is especially useful",
168 "for polarizable water. A group of reference atoms, or a center",
169 "of mass reference (option [TT]-com[tt]) and a group of solvent",
170 "atoms is required. The program splits the group of solvent atoms",
171 "into molecules. For each solvent molecule the distance to the",
172 "closest atom in reference group or to the COM is determined.",
173 "A cumulative distribution of these distances is plotted.",
174 "For each distance between [TT]-rmin[tt] and [TT]-rmax[tt]",
175 "the inner product of the distance vector",
176 "and the dipole of the solvent molecule is determined.",
177 "For solvent molecules with net charge (ions), the net charge of the ion",
178 "is subtracted evenly from all atoms in the selection of each ion.",
179 "The average of these dipole components is printed.",
180 "The same is done for the polarization, where the average dipole is",
181 "subtracted from the instantaneous dipole. The magnitude of the average",
182 "dipole is set with the option [TT]-dip[tt], the direction is defined",
183 "by the vector from the first atom in the selected solvent group",
184 "to the midpoint between the second and the third atom."
188 static gmx_bool bCom = FALSE, bPBC = FALSE;
189 static int srefat = 1;
190 static real rmin = 0.0, rmax = 0.32, refdip = 0, bw = 0.01;
192 { "-com", FALSE, etBOOL, {&bCom},
193 "Use the center of mass as the reference postion" },
194 { "-refat", FALSE, etINT, {&srefat},
195 "The reference atom of the solvent molecule" },
196 { "-rmin", FALSE, etREAL, {&rmin}, "Maximum distance (nm)" },
197 { "-rmax", FALSE, etREAL, {&rmax}, "Maximum distance (nm)" },
198 { "-dip", FALSE, etREAL, {&refdip}, "The average dipole (D)" },
199 { "-bw", FALSE, etREAL, {&bw}, "The bin width" }
203 { efTRX, NULL, NULL, ffREAD },
204 { efTPX, NULL, NULL, ffREAD },
205 { efNDX, NULL, NULL, ffOPTRD },
206 { efXVG, NULL, "scdist.xvg", ffWRITE }
208 #define NFILE asize(fnm)
210 if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_CAN_VIEW | PCA_BE_NICE,
211 NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &oenv))
218 read_tpx_top(ftp2fn(efTPX, NFILE, fnm),
219 ir, box, &natoms, NULL, NULL, NULL, top);
221 /* get index groups */
222 printf("Select a group of reference particles and a solvent group:\n");
226 get_index(&top->atoms, ftp2fn_null(efNDX, NFILE, fnm), 2, isize, index, grpname);
239 spol_atom2molindex(&(isize[1]), index[1], &(top->mols));
242 /* initialize reading trajectory: */
243 natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
245 rcut = 0.99*sqrt(max_cutoff2(ir->ePBC, box));
252 nbin = (int)(rcut*invbw)+2;
265 molindex = top->mols.index;
266 atom = top->atoms.atom;
268 gpbc = gmx_rmpbc_init(&top->idef, ir->ePBC, natoms);
270 /* start analysis of trajectory */
273 /* make molecules whole again */
274 gmx_rmpbc(gpbc, natoms, box, x);
276 set_pbc(&pbc, ir->ePBC, box);
279 calc_com_pbc(nrefat, top, x, &pbc, index[0], xref, ir->ePBC);
282 for (m = 0; m < isize[1]; m++)
286 a1 = molindex[mol+1];
287 for (i = 0; i < nrefgrp; i++)
289 pbc_dx(&pbc, x[a0+srefat], bCom ? xref : x[index[0][i]], trial);
290 rtry2 = norm2(trial);
291 if (i == 0 || rtry2 < rdx2)
293 copy_rvec(trial, dx);
299 hist[(int)(sqrt(rdx2)*invbw)+1]++;
301 if (rdx2 >= rmin2 && rdx2 < rmax2)
306 for (a = a0; a < a1; a++)
311 for (a = a0; a < a1; a++)
314 for (d = 0; d < DIM; d++)
319 for (d = 0; d < DIM; d++)
323 for (a = a0+1; a < a0+3; a++)
325 for (d = 0; d < DIM; d++)
327 dir[d] += 0.5*x[a][d];
332 svmul(ENM2DEBYE, dip, dip);
336 for (d = 0; d < DIM; d++)
338 sinp += dx[d]*dip[d];
339 sdinp += dx[d]*(dip[d] - refdip*dir[d]);
348 while (read_next_x(oenv, status, &t, x, box));
350 gmx_rmpbc_done(gpbc);
356 fprintf(stderr, "Average number of molecules within %g nm is %.1f\n",
357 rmax, (real)ntot/(real)nf);
364 fprintf(stderr, "Average dipole: %f (D), std.dev. %f\n",
365 sdip, sqrt(sdip2-sqr(sdip)));
366 fprintf(stderr, "Average radial component of the dipole: %f (D)\n",
368 fprintf(stderr, "Average radial component of the polarization: %f (D)\n",
372 fp = xvgropen(opt2fn("-o", NFILE, fnm),
373 "Cumulative solvent distribution", "r (nm)", "molecules", oenv);
375 for (i = 0; i <= nbin; i++)
378 fprintf(fp, "%g %g\n", i*bw, nmol/nf);
382 do_view(oenv, opt2fn("-o", NFILE, fnm), NULL);