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 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);
215 read_tpx_top(ftp2fn(efTPX, NFILE, fnm),
216 ir, box, &natoms, NULL, NULL, NULL, top);
218 /* get index groups */
219 printf("Select a group of reference particles and a solvent group:\n");
223 get_index(&top->atoms, ftp2fn_null(efNDX, NFILE, fnm), 2, isize, index, grpname);
236 spol_atom2molindex(&(isize[1]), index[1], &(top->mols));
239 /* initialize reading trajectory: */
240 natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
242 rcut = 0.99*sqrt(max_cutoff2(ir->ePBC, box));
249 nbin = (int)(rcut*invbw)+2;
262 molindex = top->mols.index;
263 atom = top->atoms.atom;
265 gpbc = gmx_rmpbc_init(&top->idef, ir->ePBC, natoms);
267 /* start analysis of trajectory */
270 /* make molecules whole again */
271 gmx_rmpbc(gpbc, natoms, box, x);
273 set_pbc(&pbc, ir->ePBC, box);
276 calc_com_pbc(nrefat, top, x, &pbc, index[0], xref, ir->ePBC);
279 for (m = 0; m < isize[1]; m++)
283 a1 = molindex[mol+1];
284 for (i = 0; i < nrefgrp; i++)
286 pbc_dx(&pbc, x[a0+srefat], bCom ? xref : x[index[0][i]], trial);
287 rtry2 = norm2(trial);
288 if (i == 0 || rtry2 < rdx2)
290 copy_rvec(trial, dx);
296 hist[(int)(sqrt(rdx2)*invbw)+1]++;
298 if (rdx2 >= rmin2 && rdx2 < rmax2)
303 for (a = a0; a < a1; a++)
308 for (a = a0; a < a1; a++)
311 for (d = 0; d < DIM; d++)
316 for (d = 0; d < DIM; d++)
320 for (a = a0+1; a < a0+3; a++)
322 for (d = 0; d < DIM; d++)
324 dir[d] += 0.5*x[a][d];
329 svmul(ENM2DEBYE, dip, dip);
333 for (d = 0; d < DIM; d++)
335 sinp += dx[d]*dip[d];
336 sdinp += dx[d]*(dip[d] - refdip*dir[d]);
345 while (read_next_x(oenv, status, &t, x, box));
347 gmx_rmpbc_done(gpbc);
353 fprintf(stderr, "Average number of molecules within %g nm is %.1f\n",
354 rmax, (real)ntot/(real)nf);
361 fprintf(stderr, "Average dipole: %f (D), std.dev. %f\n",
362 sdip, sqrt(sdip2-sqr(sdip)));
363 fprintf(stderr, "Average radial component of the dipole: %f (D)\n",
365 fprintf(stderr, "Average radial component of the polarization: %f (D)\n",
369 fp = xvgropen(opt2fn("-o", NFILE, fnm),
370 "Cumulative solvent distribution", "r (nm)", "molecules", oenv);
372 for (i = 0; i <= nbin; i++)
375 fprintf(fp, "%g %g\n", i*bw, nmol/nf);
379 do_view(oenv, opt2fn("-o", NFILE, fnm), NULL);