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   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   int      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
 145   const char *desc[] = {
 146     "g_spol analyzes dipoles around a solute; it is especially useful",
 147     "for polarizable water. A group of reference atoms, or a center",
 148     "of mass reference (option [TT]-com[tt]) and a group of solvent",
 149     "atoms is required. The program splits the group of solvent atoms",
 150     "into molecules. For each solvent molecule the distance to the",
 151     "closest atom in reference group or to the COM is determined.",
 152     "A cumulative distribution of these distances is plotted.",
 153     "For each distance between [TT]-rmin[tt] and [TT]-rmax[tt]",
 154     "the inner product of the distance vector",
 155     "and the dipole of the solvent molecule is determined.",
 156     "For solvent molecules with net charge (ions), the net charge of the ion",
 157     "is subtracted evenly at all atoms in the selection of each ion.",
 158     "The average of these dipole components is printed.",
 159     "The same is done for the polarization, where the average dipole is",
 160     "subtracted from the instantaneous dipole. The magnitude of the average",
 161     "dipole is set with the option [TT]-dip[tt], the direction is defined",
 162     "by the vector from the first atom in the selected solvent group",
 163     "to the midpoint between the second and the third atom."
 164   };
 165
 166   output_env_t oenv;
 167   static bool bCom = FALSE,bPBC = FALSE;
 168   static int  srefat=1;
 169   static real rmin=0.0,rmax=0.32,refdip=0,bw=0.01;
 170   t_pargs pa[] = {
 171     { "-com",  FALSE, etBOOL,  {&bCom},
 172       "Use the center of mass as the reference postion" },
 173     { "-refat",  FALSE, etINT, {&srefat},
 174       "The reference atom of the solvent molecule" },
 175     { "-rmin",  FALSE, etREAL, {&rmin}, "Maximum distance (nm)" },
 176     { "-rmax",  FALSE, etREAL, {&rmax}, "Maximum distance (nm)" },
 177     { "-dip",   FALSE, etREAL, {&refdip}, "The average dipole (D)" },
 178     { "-bw",    FALSE, etREAL, {&bw}, "The bin width" }
 179   };
 180
 181   t_filenm fnm[] = {
 182     { efTRX, NULL,  NULL,  ffREAD },
 183     { efTPX, NULL,  NULL,  ffREAD },
 184     { efNDX, NULL,  NULL,  ffOPTRD },
 185     { efXVG, NULL,  "scdist.xvg",  ffWRITE }
 186   };
 187 #define NFILE asize(fnm)
 188
 189   CopyRight(stderr,argv[0]);
 190   parse_common_args(&argc,argv,PCA_CAN_TIME | PCA_CAN_VIEW | PCA_BE_NICE,
 191                     NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL,&oenv);
 192
 193   snew(top,1);
 194   snew(ir,1);
 195   read_tpx_top(ftp2fn(efTPX,NFILE,fnm),
 196                ir,box,&natoms,NULL,NULL,NULL,top);
 197
 198   /* get index groups */
 199   printf("Select a group of reference particles and a solvent group:\n");
 200   snew(grpname,2);
 201   snew(index,2);
 202   snew(isize,2);
 203   get_index(&top->atoms,ftp2fn_null(efNDX,NFILE,fnm),2,isize,index,grpname);
 204
 205   if (bCom) {
 206     nrefgrp = 1;
 207     nrefat  = isize[0];
 208   } else {
 209     nrefgrp = isize[0];
 210     nrefat  = 1;
 211   }
 212
 213   spol_atom2molindex(&(isize[1]),index[1],&(top->mols));
 214   srefat--;
 215
 216   /* initialize reading trajectory:                         */
 217   natoms=read_first_x(oenv,&status,ftp2fn(efTRX,NFILE,fnm),&t,&x,box);
 218
 219   rcut  = 0.99*sqrt(max_cutoff2(ir->ePBC,box));
 220   if (rcut == 0)
 221     rcut = 10*rmax;
 222   rcut2 = sqr(rcut);
 223   invbw = 1/bw;
 224   nbin = (int)(rcut*invbw)+2;
 225   snew(hist,nbin);
 226
 227   rmin2 = sqr(rmin);
 228   rmax2 = sqr(rmax);
 229
 230   nf = 0;
 231   ntot = 0;
 232   sdip  = 0;
 233   sdip2 = 0;
 234   sinp  = 0;
 235   sdinp = 0;
 236
 237   molindex = top->mols.index;
 238   atom     = top->atoms.atom;
 239
 240   /* start analysis of trajectory */
 241   do {
 242     /* make molecules whole again */
 243     rm_pbc(&top->idef,ir->ePBC,natoms,box,x,x);
 244
 245     set_pbc(&pbc,ir->ePBC,box);
 246     if (bCom)
 247       calc_com_pbc(nrefat,top,x,&pbc,index[0],xref,ir->ePBC,box);
 248
 249     for(m=0; m<isize[1]; m++) {
 250       mol = index[1][m];
 251       a0 = molindex[mol];
 252       a1 = molindex[mol+1];
 253       for(i=0; i<nrefgrp; i++) {
 254         pbc_dx(&pbc,x[a0+srefat],bCom ? xref : x[index[0][i]],trial);
 255         rtry2 = norm2(trial);
 256         if (i==0 || rtry2 < rdx2) {
 257           copy_rvec(trial,dx);
 258           rdx2 = rtry2;
 259         }
 260       }
 261       if (rdx2 < rcut2)
 262         hist[(int)(sqrt(rdx2)*invbw)+1]++;
 263       if (rdx2 >= rmin2 && rdx2 < rmax2) {
 264         unitv(dx,dx);
 265         clear_rvec(dip);
 266         qav = 0;
 267         for(a=a0; a<a1; a++) {
 268           qav += atom[a].q;
 269         }
 270         qav /= (a1 - a0);
 271         for(a=a0; a<a1; a++) {
 272           q = atom[a].q - qav;
 273           for(d=0; d<DIM; d++)
 274             dip[d] += q*x[a][d];
 275         }
 276         for(d=0; d<DIM; d++)
 277           dir[d] = -x[a0][d];
 278         for(a=a0+1; a<a0+3; a++) {
 279           for(d=0; d<DIM; d++)
 280             dir[d] += 0.5*x[a][d];
 281         }
 282         unitv(dir,dir);
 283
 284         svmul(ENM2DEBYE,dip,dip);
 285         dip2 = norm2(dip);
 286         sdip  += sqrt(dip2);
 287         sdip2 += dip2;
 288         for(d=0; d<DIM; d++) {
 289           sinp  += dx[d]*dip[d];
 290           sdinp += dx[d]*(dip[d] - refdip*dir[d]);
 291         }
 292
 293         ntot++;
 294       }
 295     }
 296     nf++;
 297
 298   }  while (read_next_x(oenv,status,&t,natoms,x,box));
 299
 300   /* clean up */
 301   sfree(x);
 302   close_trj(status);
 303
 304   fprintf(stderr,"Average number of molecules within %g nm is %.1f\n",
 305           rmax,(real)ntot/(real)nf);
 306   if (ntot > 0) {
 307     sdip  /= ntot;
 308     sdip2 /= ntot;
 309     sinp  /= ntot;
 310     sdinp /= ntot;
 311     fprintf(stderr,"Average dipole:                               %f (D), std.dev. %f\n",
 312             sdip,sqrt(sdip2-sqr(sdip)));
 313     fprintf(stderr,"Average radial component of the dipole:       %f (D)\n",
 314             sinp);
 315     fprintf(stderr,"Average radial component of the polarization: %f (D)\n",
 316             sdinp);
 317   }
 318
 319   fp = xvgropen(opt2fn("-o",NFILE,fnm),
 320                 "Cumulative solvent distribution","r (nm)","molecules",oenv);
 321   nmol = 0;
 322   for(i=0; i<=nbin; i++) {
 323     nmol += hist[i];
 324     fprintf(fp,"%g %g\n",i*bw,nmol/nf);
 325   }
 326   ffclose(fp);
 327
 328   do_view(oenv,opt2fn("-o",NFILE,fnm),NULL);
 329
 330   thanx(stderr);
 331
 332   return 0;
 333 }