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
53 static void calc_com_pbc(int nrefat, t_topology *top, rvec x[], t_pbc *pbc,
54 atom_id index[], rvec xref, int ePBC, matrix box)
56 const real tol = 1e-4;
62 /* First simple calculation */
65 for (m = 0; (m < nrefat); m++)
68 mass = top->atoms.atom[ai].m;
69 for (j = 0; (j < DIM); j++)
71 xref[j] += mass*x[ai][j];
75 svmul(1/mtot, xref, xref);
76 /* Now check if any atom is more than half the box from the COM */
83 for (m = 0; (m < nrefat); m++)
86 mass = top->atoms.atom[ai].m/mtot;
87 pbc_dx(pbc, x[ai], xref, dx);
88 rvec_add(xref, dx, xtest);
89 for (j = 0; (j < DIM); j++)
91 if (fabs(xtest[j]-x[ai][j]) > tol)
93 /* Here we have used the wrong image for contributing to the COM */
94 xref[j] += mass*(xtest[j]-x[ai][j]);
102 printf("COM: %8.3f %8.3f %8.3f iter = %d\n", xref[XX], xref[YY], xref[ZZ], iter);
110 void spol_atom2molindex(int *n, int *index, t_block *mols)
119 while (m < mols->nr && index[i] != mols->index[m])
125 gmx_fatal(FARGS, "index[%d]=%d does not correspond to the first atom of a molecule", i+1, index[i]+1);
127 for (j = mols->index[m]; j < mols->index[m+1]; j++)
129 if (i >= *n || index[i] != j)
131 gmx_fatal(FARGS, "The index group is not a set of whole molecules");
135 /* Modify the index in place */
138 printf("There are %d molecules in the selection\n", nmol);
143 int gmx_spol(int argc, char *argv[])
150 int nrefat, natoms, nf, ntot;
152 rvec *xtop, *x, xref, trial, dx = {0}, dip, dir;
157 atom_id **index, *molindex;
159 real rmin2, rmax2, rcut, rcut2, rdx2 = 0, rtry2, qav, q, dip2, invbw;
160 int nbin, i, m, mol, a0, a1, a, d;
161 double sdip, sdip2, sinp, sdinp, nmol;
164 gmx_rmpbc_t gpbc = NULL;
167 const char *desc[] = {
168 "[TT]g_spol[tt] analyzes dipoles around a solute; it is especially useful",
169 "for polarizable water. A group of reference atoms, or a center",
170 "of mass reference (option [TT]-com[tt]) and a group of solvent",
171 "atoms is required. The program splits the group of solvent atoms",
172 "into molecules. For each solvent molecule the distance to the",
173 "closest atom in reference group or to the COM is determined.",
174 "A cumulative distribution of these distances is plotted.",
175 "For each distance between [TT]-rmin[tt] and [TT]-rmax[tt]",
176 "the inner product of the distance vector",
177 "and the dipole of the solvent molecule is determined.",
178 "For solvent molecules with net charge (ions), the net charge of the ion",
179 "is subtracted evenly from all atoms in the selection of each ion.",
180 "The average of these dipole components is printed.",
181 "The same is done for the polarization, where the average dipole is",
182 "subtracted from the instantaneous dipole. The magnitude of the average",
183 "dipole is set with the option [TT]-dip[tt], the direction is defined",
184 "by the vector from the first atom in the selected solvent group",
185 "to the midpoint between the second and the third atom."
189 static gmx_bool bCom = FALSE, bPBC = FALSE;
190 static int srefat = 1;
191 static real rmin = 0.0, rmax = 0.32, refdip = 0, bw = 0.01;
193 { "-com", FALSE, etBOOL, {&bCom},
194 "Use the center of mass as the reference postion" },
195 { "-refat", FALSE, etINT, {&srefat},
196 "The reference atom of the solvent molecule" },
197 { "-rmin", FALSE, etREAL, {&rmin}, "Maximum distance (nm)" },
198 { "-rmax", FALSE, etREAL, {&rmax}, "Maximum distance (nm)" },
199 { "-dip", FALSE, etREAL, {&refdip}, "The average dipole (D)" },
200 { "-bw", FALSE, etREAL, {&bw}, "The bin width" }
204 { efTRX, NULL, NULL, ffREAD },
205 { efTPX, NULL, NULL, ffREAD },
206 { efNDX, NULL, NULL, ffOPTRD },
207 { efXVG, NULL, "scdist.xvg", ffWRITE }
209 #define NFILE asize(fnm)
211 parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_CAN_VIEW | PCA_BE_NICE,
212 NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &oenv);
216 read_tpx_top(ftp2fn(efTPX, NFILE, fnm),
217 ir, box, &natoms, NULL, NULL, NULL, top);
219 /* get index groups */
220 printf("Select a group of reference particles and a solvent group:\n");
224 get_index(&top->atoms, ftp2fn_null(efNDX, NFILE, fnm), 2, isize, index, grpname);
237 spol_atom2molindex(&(isize[1]), index[1], &(top->mols));
240 /* initialize reading trajectory: */
241 natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
243 rcut = 0.99*sqrt(max_cutoff2(ir->ePBC, box));
250 nbin = (int)(rcut*invbw)+2;
263 molindex = top->mols.index;
264 atom = top->atoms.atom;
266 gpbc = gmx_rmpbc_init(&top->idef, ir->ePBC, natoms, box);
268 /* start analysis of trajectory */
271 /* make molecules whole again */
272 gmx_rmpbc(gpbc, natoms, box, x);
274 set_pbc(&pbc, ir->ePBC, box);
277 calc_com_pbc(nrefat, top, x, &pbc, index[0], xref, ir->ePBC, box);
280 for (m = 0; m < isize[1]; m++)
284 a1 = molindex[mol+1];
285 for (i = 0; i < nrefgrp; i++)
287 pbc_dx(&pbc, x[a0+srefat], bCom ? xref : x[index[0][i]], trial);
288 rtry2 = norm2(trial);
289 if (i == 0 || rtry2 < rdx2)
291 copy_rvec(trial, dx);
297 hist[(int)(sqrt(rdx2)*invbw)+1]++;
299 if (rdx2 >= rmin2 && rdx2 < rmax2)
304 for (a = a0; a < a1; a++)
309 for (a = a0; a < a1; a++)
312 for (d = 0; d < DIM; d++)
317 for (d = 0; d < DIM; d++)
321 for (a = a0+1; a < a0+3; a++)
323 for (d = 0; d < DIM; d++)
325 dir[d] += 0.5*x[a][d];
330 svmul(ENM2DEBYE, dip, dip);
334 for (d = 0; d < DIM; d++)
336 sinp += dx[d]*dip[d];
337 sdinp += dx[d]*(dip[d] - refdip*dir[d]);
346 while (read_next_x(oenv, status, &t, natoms, x, box));
348 gmx_rmpbc_done(gpbc);
354 fprintf(stderr, "Average number of molecules within %g nm is %.1f\n",
355 rmax, (real)ntot/(real)nf);
362 fprintf(stderr, "Average dipole: %f (D), std.dev. %f\n",
363 sdip, sqrt(sdip2-sqr(sdip)));
364 fprintf(stderr, "Average radial component of the dipole: %f (D)\n",
366 fprintf(stderr, "Average radial component of the polarization: %f (D)\n",
370 fp = xvgropen(opt2fn("-o", NFILE, fnm),
371 "Cumulative solvent distribution", "r (nm)", "molecules", oenv);
373 for (i = 0; i <= nbin; i++)
376 fprintf(fp, "%g %g\n", i*bw, nmol/nf);
380 do_view(oenv, opt2fn("-o", NFILE, fnm), NULL);