src/tools/gmx_spol.c

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