src/tools/gmx_spol.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   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, 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.
  11  *
  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.
  16  *
  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.
  21  *
  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.
  26  *
  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.
  34  *
  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.
  37  */
  38 #ifdef HAVE_CONFIG_H
  39 #include <config.h>
  40 #endif
  41
  42 #include "macros.h"
  43 #include "statutil.h"
  44 #include "smalloc.h"
  45 #include "copyrite.h"
  46 #include "gstat.h"
  47 #include "vec.h"
  48 #include "xvgr.h"
  49 #include "pbc.h"
  50 #include "index.h"
  51 #include "tpxio.h"
  52 #include "physics.h"
  53 #include "gmx_ana.h"
  54
  55
  56 static void calc_com_pbc(int nrefat,t_topology *top,rvec x[],t_pbc *pbc,
  57                          atom_id index[],rvec xref,int ePBC,matrix box)
  58 {
  59   const real tol=1e-4;
  60   gmx_bool  bChanged;
  61   int   m,j,ai,iter;
  62   real  mass,mtot;
  63   rvec  dx,xtest;
  64
  65   /* First simple calculation */
  66   clear_rvec(xref);
  67   mtot = 0;
  68   for(m=0; (m<nrefat); m++) {
  69     ai = index[m];
  70     mass = top->atoms.atom[ai].m;
  71     for(j=0; (j<DIM); j++)
  72       xref[j] += mass*x[ai][j];
  73     mtot += mass;
  74   }
  75   svmul(1/mtot,xref,xref);
  76   /* Now check if any atom is more than half the box from the COM */
  77   if (ePBC != epbcNONE) {
  78     iter = 0;
  79     do {
  80       bChanged = FALSE;
  81       for(m=0; (m<nrefat); m++) {
  82         ai   = index[m];
  83         mass = top->atoms.atom[ai].m/mtot;
  84         pbc_dx(pbc,x[ai],xref,dx);
  85         rvec_add(xref,dx,xtest);
  86         for(j=0; (j<DIM); j++)
  87           if (fabs(xtest[j]-x[ai][j]) > tol) {
  88             /* Here we have used the wrong image for contributing to the COM */
  89             xref[j] += mass*(xtest[j]-x[ai][j]);
  90             x[ai][j] = xtest[j];
  91             bChanged = TRUE;
  92           }
  93       }
  94       if (bChanged)
  95         printf("COM: %8.3f  %8.3f  %8.3f  iter = %d\n",xref[XX],xref[YY],xref[ZZ],iter);
  96       iter++;
  97     } while (bChanged);
  98   }
  99 }
 100
 101 void spol_atom2molindex(int *n,int *index,t_block *mols)
 102 {
 103   int nmol,i,j,m;
 104
 105   nmol = 0;
 106   i=0;
 107   while (i < *n) {
 108     m=0;
 109     while(m < mols->nr && index[i] != mols->index[m])
 110       m++;
 111     if (m == mols->nr)
 112       gmx_fatal(FARGS,"index[%d]=%d does not correspond to the first atom of a molecule",i+1,index[i]+1);
 113     for(j=mols->index[m]; j<mols->index[m+1]; j++) {
 114       if (i >= *n || index[i] != j)
 115         gmx_fatal(FARGS,"The index group is not a set of whole molecules");
 116       i++;
 117     }
 118     /* Modify the index in place */
 119     index[nmol++] = m;
 120   }
 121   printf("There are %d molecules in the selection\n",nmol);
 122
 123   *n = nmol;
 124 }
 125
 126 int gmx_spol(int argc,char *argv[])
 127 {
 128   t_topology *top;
 129   t_inputrec *ir;
 130   t_atom     *atom;
 131   char     title[STRLEN];
 132   t_trxstatus *status;
 133   int      nrefat,natoms,nf,ntot;
 134   real     t;
 135   rvec     *xtop,*x,xref,trial,dx={0},dip,dir;
 136   matrix   box;
 137
 138   FILE    *fp;
 139   int     *isize,nrefgrp;
 140   atom_id **index,*molindex;
 141   char    **grpname;
 142   real    rmin2,rmax2,rcut,rcut2,rdx2=0,rtry2,qav,q,dip2,invbw;
 143   int     nbin,i,m,mol,a0,a1,a,d;
 144   double  sdip,sdip2,sinp,sdinp,nmol;
 145   int     *hist;
 146   t_pbc   pbc;
 147   gmx_rmpbc_t  gpbc=NULL;
 148
 149
 150   const char *desc[] = {
 151     "[TT]g_spol[tt] analyzes dipoles around a solute; it is especially useful",
 152     "for polarizable water. A group of reference atoms, or a center",
 153     "of mass reference (option [TT]-com[tt]) and a group of solvent",
 154     "atoms is required. The program splits the group of solvent atoms",
 155     "into molecules. For each solvent molecule the distance to the",
 156     "closest atom in reference group or to the COM is determined.",
 157     "A cumulative distribution of these distances is plotted.",
 158     "For each distance between [TT]-rmin[tt] and [TT]-rmax[tt]",
 159     "the inner product of the distance vector",
 160     "and the dipole of the solvent molecule is determined.",
 161     "For solvent molecules with net charge (ions), the net charge of the ion",
 162     "is subtracted evenly from all atoms in the selection of each ion.",
 163     "The average of these dipole components is printed.",
 164     "The same is done for the polarization, where the average dipole is",
 165     "subtracted from the instantaneous dipole. The magnitude of the average",
 166     "dipole is set with the option [TT]-dip[tt], the direction is defined",
 167     "by the vector from the first atom in the selected solvent group",
 168     "to the midpoint between the second and the third atom."
 169   };
 170
 171   output_env_t oenv;
 172   static gmx_bool bCom = FALSE,bPBC = FALSE;
 173   static int  srefat=1;
 174   static real rmin=0.0,rmax=0.32,refdip=0,bw=0.01;
 175   t_pargs pa[] = {
 176     { "-com",  FALSE, etBOOL,  {&bCom},
 177       "Use the center of mass as the reference postion" },
 178     { "-refat",  FALSE, etINT, {&srefat},
 179       "The reference atom of the solvent molecule" },
 180     { "-rmin",  FALSE, etREAL, {&rmin}, "Maximum distance (nm)" },
 181     { "-rmax",  FALSE, etREAL, {&rmax}, "Maximum distance (nm)" },
 182     { "-dip",   FALSE, etREAL, {&refdip}, "The average dipole (D)" },
 183     { "-bw",    FALSE, etREAL, {&bw}, "The bin width" }
 184   };
 185
 186   t_filenm fnm[] = {
 187     { efTRX, NULL,  NULL,  ffREAD },
 188     { efTPX, NULL,  NULL,  ffREAD },
 189     { efNDX, NULL,  NULL,  ffOPTRD },
 190     { efXVG, NULL,  "scdist.xvg",  ffWRITE }
 191   };
 192 #define NFILE asize(fnm)
 193
 194   CopyRight(stderr,argv[0]);
 195   parse_common_args(&argc,argv,PCA_CAN_TIME | PCA_CAN_VIEW | PCA_BE_NICE,
 196                     NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL,&oenv);
 197
 198   snew(top,1);
 199   snew(ir,1);
 200   read_tpx_top(ftp2fn(efTPX,NFILE,fnm),
 201                ir,box,&natoms,NULL,NULL,NULL,top);
 202
 203   /* get index groups */
 204   printf("Select a group of reference particles and a solvent group:\n");
 205   snew(grpname,2);
 206   snew(index,2);
 207   snew(isize,2);
 208   get_index(&top->atoms,ftp2fn_null(efNDX,NFILE,fnm),2,isize,index,grpname);
 209
 210   if (bCom) {
 211     nrefgrp = 1;
 212     nrefat  = isize[0];
 213   } else {
 214     nrefgrp = isize[0];
 215     nrefat  = 1;
 216   }
 217
 218   spol_atom2molindex(&(isize[1]),index[1],&(top->mols));
 219   srefat--;
 220
 221   /* initialize reading trajectory:                         */
 222   natoms=read_first_x(oenv,&status,ftp2fn(efTRX,NFILE,fnm),&t,&x,box);
 223
 224   rcut  = 0.99*sqrt(max_cutoff2(ir->ePBC,box));
 225   if (rcut == 0)
 226     rcut = 10*rmax;
 227   rcut2 = sqr(rcut);
 228   invbw = 1/bw;
 229   nbin = (int)(rcut*invbw)+2;
 230   snew(hist,nbin);
 231
 232   rmin2 = sqr(rmin);
 233   rmax2 = sqr(rmax);
 234
 235   nf = 0;
 236   ntot = 0;
 237   sdip  = 0;
 238   sdip2 = 0;
 239   sinp  = 0;
 240   sdinp = 0;
 241
 242   molindex = top->mols.index;
 243   atom     = top->atoms.atom;
 244
 245   gpbc = gmx_rmpbc_init(&top->idef,ir->ePBC,natoms,box);
 246
 247   /* start analysis of trajectory */
 248   do {
 249     /* make molecules whole again */
 250     gmx_rmpbc(gpbc,natoms,box,x);
 251
 252     set_pbc(&pbc,ir->ePBC,box);
 253     if (bCom)
 254       calc_com_pbc(nrefat,top,x,&pbc,index[0],xref,ir->ePBC,box);
 255
 256     for(m=0; m<isize[1]; m++) {
 257       mol = index[1][m];
 258       a0 = molindex[mol];
 259       a1 = molindex[mol+1];
 260       for(i=0; i<nrefgrp; i++) {
 261         pbc_dx(&pbc,x[a0+srefat],bCom ? xref : x[index[0][i]],trial);
 262         rtry2 = norm2(trial);
 263         if (i==0 || rtry2 < rdx2) {
 264           copy_rvec(trial,dx);
 265           rdx2 = rtry2;
 266         }
 267       }
 268       if (rdx2 < rcut2)
 269         hist[(int)(sqrt(rdx2)*invbw)+1]++;
 270       if (rdx2 >= rmin2 && rdx2 < rmax2) {
 271         unitv(dx,dx);
 272         clear_rvec(dip);
 273         qav = 0;
 274         for(a=a0; a<a1; a++) {
 275           qav += atom[a].q;
 276         }
 277         qav /= (a1 - a0);
 278         for(a=a0; a<a1; a++) {
 279           q = atom[a].q - qav;
 280           for(d=0; d<DIM; d++)
 281             dip[d] += q*x[a][d];
 282         }
 283         for(d=0; d<DIM; d++)
 284           dir[d] = -x[a0][d];
 285         for(a=a0+1; a<a0+3; a++) {
 286           for(d=0; d<DIM; d++)
 287             dir[d] += 0.5*x[a][d];
 288         }
 289         unitv(dir,dir);
 290
 291         svmul(ENM2DEBYE,dip,dip);
 292         dip2 = norm2(dip);
 293         sdip  += sqrt(dip2);
 294         sdip2 += dip2;
 295         for(d=0; d<DIM; d++) {
 296           sinp  += dx[d]*dip[d];
 297           sdinp += dx[d]*(dip[d] - refdip*dir[d]);
 298         }
 299
 300         ntot++;
 301       }
 302     }
 303     nf++;
 304
 305   }  while (read_next_x(oenv,status,&t,natoms,x,box));
 306
 307   gmx_rmpbc_done(gpbc);
 308
 309   /* clean up */
 310   sfree(x);
 311   close_trj(status);
 312
 313   fprintf(stderr,"Average number of molecules within %g nm is %.1f\n",
 314           rmax,(real)ntot/(real)nf);
 315   if (ntot > 0) {
 316     sdip  /= ntot;
 317     sdip2 /= ntot;
 318     sinp  /= ntot;
 319     sdinp /= ntot;
 320     fprintf(stderr,"Average dipole:                               %f (D), std.dev. %f\n",
 321             sdip,sqrt(sdip2-sqr(sdip)));
 322     fprintf(stderr,"Average radial component of the dipole:       %f (D)\n",
 323             sinp);
 324     fprintf(stderr,"Average radial component of the polarization: %f (D)\n",
 325             sdinp);
 326   }
 327
 328   fp = xvgropen(opt2fn("-o",NFILE,fnm),
 329                 "Cumulative solvent distribution","r (nm)","molecules",oenv);
 330   nmol = 0;
 331   for(i=0; i<=nbin; i++) {
 332     nmol += hist[i];
 333     fprintf(fp,"%g %g\n",i*bw,nmol/nf);
 334   }
 335   ffclose(fp);
 336
 337   do_view(oenv,opt2fn("-o",NFILE,fnm),NULL);
 338
 339   thanx(stderr);
 340
 341   return 0;
 342 }