src/contrib/ehanal.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-2008, 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 <stdlib.h>
  43 #include <math.h>
  44 #include <string.h>
  45 #include "typedefs.h"
  46 #include "smalloc.h"
  47 #include "macros.h"
  48 #include "copyrite.h"
  49 #include "statutil.h"
  50 #include "gmx_fatal.h"
  51 #include "random.h"
  52 #include "pdbio.h"
  53 #include "futil.h"
  54 #include "physics.h"
  55 #include "xvgr.h"
  56 #include "vec.h"
  57 #include "names.h"
  58 #include "ehdata.h"
  59 #include "pdbio.h"
  60
  61 t_histo *init_histo(int np,real minx,real maxx)
  62 {
  63   t_histo *h;
  64
  65   snew(h,1);
  66   snew(h->y,np+1);
  67   snew(h->nh,np+1);
  68   h->np   = np;
  69   if (maxx <= minx)
  70     gmx_fatal(FARGS,"minx (%f) should be less than maxx (%f) in init_histo",minx,maxx);
  71   h->minx = minx;
  72   h->maxx = maxx;
  73   h->dx_1 = np/(maxx-minx);
  74   h->dx   = (maxx-minx)/np;
  75
  76   return h;
  77 }
  78
  79 void done_histo(t_histo *h)
  80 {
  81   sfree(h->y);
  82   sfree(h->nh);
  83   h->np = 0;
  84 }
  85
  86 void add_histo(t_histo *h,real x,real y)
  87 {
  88   int n;
  89
  90   n = (x-h->minx)*h->dx_1;
  91   if ((n < 0) || (n > h->np))
  92     gmx_fatal(FARGS,"Invalid x (%f) in add_histo. Should be in %f - %f",x,h->minx,h->maxx);
  93   h->y[n] += y;
  94   h->nh[n]++;
  95 }
  96
  97 void dump_histo(t_histo *h,char *fn,char *title,char *xaxis,char *yaxis,
  98                 int enorm,real norm_fac)
  99 {
 100   FILE *fp;
 101   int  i,nn;
 102
 103   for(nn=h->np; (nn > 0); nn--)
 104     if (h->nh[nn] != 0)
 105       break;
 106   for(i=0; (i<nn); i++)
 107     if (h->nh[i] > 0)
 108       break;
 109   fp = xvgropen(fn,title,xaxis,yaxis);
 110   for(  ; (i<nn); i++) {
 111     switch (enorm) {
 112     case enormNO:
 113       fprintf(fp,"%12f  %12f  %d\n",h->minx+h->dx*i,h->y[i],h->nh[i]);
 114       break;
 115     case enormFAC:
 116       fprintf(fp,"%12f  %12f  %d\n",h->minx+h->dx*i,h->y[i]*norm_fac,h->nh[i]);
 117       break;
 118     case enormNP:
 119       if (h->nh[i] > 0)
 120         fprintf(fp,"%12f  %12f  %d\n",
 121                 h->minx+h->dx*i,h->y[i]*norm_fac/h->nh[i],h->nh[i]);
 122       break;
 123     default:
 124       gmx_fatal(FARGS,"Wrong value for enorm (%d)",enorm);
 125     }
 126   }
 127   ffclose(fp);
 128 }
 129
 130 /*******************************************************************
 131  *
 132  * Functions to analyse and monitor scattering
 133  *
 134  *******************************************************************/
 135
 136 void add_scatter_event(t_ana_scat *scatter,rvec pos,gmx_bool bInel,
 137                        real t,real ekin)
 138 {
 139   int np = scatter->np;
 140
 141   if (np == scatter->maxp) {
 142     scatter->maxp += 32;
 143     srenew(scatter->time,scatter->maxp);
 144     srenew(scatter->ekin,scatter->maxp);
 145     srenew(scatter->bInel,scatter->maxp);
 146     srenew(scatter->pos,scatter->maxp);
 147   }
 148   scatter->time[np]  = t;
 149   scatter->bInel[np] = np;
 150   scatter->ekin[np]  = ekin;
 151   copy_rvec(pos,scatter->pos[np]);
 152   scatter->np++;
 153 }
 154
 155 void reset_ana_scat(t_ana_scat *scatter)
 156 {
 157   scatter->np = 0;
 158 }
 159
 160 void done_scatter(t_ana_scat *scatter)
 161 {
 162   scatter->np = 0;
 163   sfree(scatter->time);
 164   sfree(scatter->ekin);
 165   sfree(scatter->bInel);
 166   sfree(scatter->pos);
 167 }
 168
 169 void analyse_scatter(t_ana_scat *scatter,t_histo *hmfp)
 170 {
 171   int   i,n;
 172   rvec  dx;
 173
 174   if (scatter->np > 1) {
 175     for(i=1; (i<scatter->np); i++) {
 176       rvec_sub(scatter->pos[i],scatter->pos[i-1],dx);
 177       add_histo(hmfp,scatter->ekin[i],norm(dx));
 178     }
 179   }
 180 }
 181
 182 /*******************************************************************
 183  *
 184  * Functions to analyse structural changes
 185  *
 186  *******************************************************************/
 187
 188 t_ana_struct *init_ana_struct(int nstep,int nsave,real timestep,
 189                               int maxparticle)
 190 {
 191   t_ana_struct *anal;
 192
 193   snew(anal,1);
 194   anal->nanal = 1.2*((nstep / nsave)+1);
 195   anal->index = 0;
 196   anal->dt    = nsave*timestep;
 197   snew(anal->t,anal->nanal);
 198   snew(anal->maxdist,anal->nanal);
 199   snew(anal->d2_com,anal->nanal);
 200   snew(anal->d2_origin,anal->nanal);
 201   snew(anal->nion,anal->nanal);
 202   anal->nstruct   = 1;
 203   anal->nparticle = 1;
 204   anal->maxparticle = maxparticle;
 205   snew(anal->q,1);
 206   snew(anal->x,1);
 207   snew(anal->x[0],maxparticle);
 208
 209   return anal;
 210 }
 211
 212 void done_ana_struct(t_ana_struct *anal)
 213 {
 214   int i;
 215
 216   sfree(anal->t);
 217   sfree(anal->maxdist);
 218   sfree(anal->d2_com);
 219   sfree(anal->d2_origin);
 220   sfree(anal->nion);
 221   sfree(anal->q);
 222   for(i=0; (i<anal->nstruct); i++)
 223     sfree(anal->x[i]);
 224   sfree(anal->x);
 225 }
 226
 227 void reset_ana_struct(t_ana_struct *anal)
 228 {
 229   int i;
 230
 231   for(i=0; (i<anal->nanal); i++) {
 232     anal->t[i] = 0;
 233     anal->maxdist[i] = 0;
 234     clear_rvec(anal->d2_com[i]);
 235     clear_rvec(anal->d2_origin[i]);
 236     anal->nion[i] = 0;
 237   }
 238   anal->index = 0;
 239 }
 240
 241 void add_ana_struct(t_ana_struct *total,t_ana_struct *add)
 242 {
 243   int i,m;
 244
 245   if (total->index == 0)
 246     total->index = add->index;
 247   else if (total->index != add->index)
 248     gmx_fatal(FARGS,"Analysis incompatible (total: %d, add: %d) %s, %d",
 249                 total->index,add->index,__FILE__,__LINE__);
 250   for(i=0; (i<total->index); i++) {
 251     if (total->t[i] == 0)
 252       total->t[i] = add->t[i];
 253     else if (total->t[i] != add->t[i])
 254       gmx_fatal(FARGS,"Inconsistent times in analysis (%f-%f) %s, %d",
 255                   total->t[i],add->t[i],__FILE__,__LINE__);
 256     if (add->maxdist[i] > total->maxdist[i])
 257       total->maxdist[i]  = add->maxdist[i];
 258     for(m=0; (m<DIM); m++) {
 259       total->d2_com[i][m]    += add->d2_com[i][m]/add->nion[i];
 260       total->d2_origin[i][m] += add->d2_origin[i][m]/add->nion[i];
 261     }
 262     total->nion[i]     += add->nion[i];
 263   }
 264 }
 265
 266 static void do_add_struct(t_ana_struct *anal,int nparticle,rvec x[])
 267 {
 268   int i,j;
 269
 270   if (nparticle > anal->nparticle) {
 271     for(i=0; (i<anal->nstruct); i++) {
 272       for(j=anal->nparticle; (j<nparticle); j++)
 273         copy_rvec(x[j],anal->x[i][j]);
 274     }
 275   }
 276   anal->nparticle=nparticle;
 277   srenew(anal->x,anal->nstruct+1);
 278   snew(anal->x[anal->nstruct],anal->maxparticle);
 279   for(j=0; (j<nparticle); j++)
 280     copy_rvec(x[j],anal->x[anal->nstruct][j]);
 281   anal->nstruct++;
 282 }
 283
 284 void analyse_structure(t_ana_struct *anal,real t,rvec center,
 285                        rvec x[],int nparticle,real charge[])
 286 {
 287   int  i,j,m,nel,n=0;
 288   rvec dx,com;
 289   real dx2,dx1;
 290
 291   j = anal->index;
 292   if (j >= anal->nanal)
 293     gmx_fatal(FARGS,"Too many points in analyse_structure");
 294   anal->t[j]       = t;
 295   anal->maxdist[j] = 0;
 296   clear_rvec(com);
 297   nel = 0;
 298   for(i=0; (i<nparticle); i++) {
 299     if (charge[i] < 0) {
 300       rvec_inc(com,x[i]);
 301       nel++;
 302     }
 303   }
 304   if (nel > 0)
 305     for(m=0; (m<3); m++)
 306       com[m] /= nel;
 307   for(i=0; (i<nparticle); i++) {
 308     if (charge[i] < 0) {
 309       rvec_sub(x[i],com,dx);
 310       for(m=0; (m<DIM); m++) {
 311         anal->d2_com[j][m]    += sqr(dx[m]);
 312         anal->d2_origin[j][m] += sqr(x[i][m]);
 313       }
 314       dx2 = iprod(x[i],x[i]);
 315       dx1 = sqrt(dx2);
 316       if (dx1 > anal->maxdist[j])
 317         anal->maxdist[j] = dx1;
 318       n++;
 319     }
 320   }
 321   do_add_struct(anal,nparticle,x);
 322   anal->nion[j] = n;
 323   anal->index++;
 324 }
 325
 326 void dump_ana_struct(char *rmax,char *nion,char *gyr_com,char *gyr_origin,
 327                      t_ana_struct *anal,int nsim)
 328 {
 329   FILE *fp,*gp,*hp,*kp;
 330   int  i,j;
 331   real t,d2;
 332   char *legend[] = { "Rg", "RgX", "RgY", "RgZ" };
 333
 334   fp = xvgropen(rmax,"rmax","Time (fs)","r (nm)");
 335   gp = xvgropen(nion,"N ion","Time (fs)","N ions");
 336   hp = xvgropen(gyr_com,"Radius of gyration wrt. C.O.M.",
 337                 "Time (fs)","Rg (nm)");
 338   xvgr_legend(hp,asize(legend),legend);
 339   kp = xvgropen(gyr_origin,"Radius of gyration wrt. Origin",
 340                 "Time (fs)","Rg (nm)");
 341   xvgr_legend(kp,asize(legend),legend);
 342   for(i=0; (i<anal->index); i++) {
 343     t = 1000*anal->t[i];
 344     fprintf(fp,"%12g  %10.3f\n",t,anal->maxdist[i]);
 345     fprintf(gp,"%12g  %10.3f\n",t,(1.0*anal->nion[i])/nsim-1);
 346     d2 = anal->d2_com[i][XX] + anal->d2_com[i][YY] + anal->d2_com[i][ZZ];
 347     fprintf(hp,"%12g  %10.3f  %10.3f  %10.3f  %10.3f\n",
 348             t,sqrt(d2/nsim),
 349             sqrt(anal->d2_com[i][XX]/nsim),
 350             sqrt(anal->d2_com[i][YY]/nsim),
 351             sqrt(anal->d2_com[i][ZZ]/nsim));
 352     d2 = anal->d2_origin[i][XX] + anal->d2_origin[i][YY] + anal->d2_origin[i][ZZ];
 353     fprintf(kp,"%12g  %10.3f  %10.3f  %10.3f  %10.3f\n",
 354             t,sqrt(d2/nsim),
 355             sqrt(anal->d2_origin[i][XX]/nsim),
 356             sqrt(anal->d2_origin[i][YY]/nsim),
 357             sqrt(anal->d2_origin[i][ZZ]/nsim));
 358   }
 359   ffclose(hp);
 360   ffclose(gp);
 361   ffclose(fp);
 362   ffclose(kp);
 363 }
 364
 365 void dump_as_pdb(char *pdb,t_ana_struct *anal)
 366 {
 367   FILE *kp;
 368   int  i,j;
 369   real t;
 370
 371   kp = ffopen(pdb,"w");
 372   for(i=0; (i<anal->nstruct); i++) {
 373     t = 1000*anal->t[i];
 374     fprintf(kp,"MODEL  %d  time %g fs\n",i+1,t);
 375     for(j=0; (j<anal->nparticle); j++) {
 376       fprintf(kp,pdbformat,"ATOM",i+1,(j < anal->nion[i]) ? "O" : "N",
 377               "PLS",' ',1,
 378               anal->x[i][j][XX]/100,
 379               anal->x[i][j][YY]/100,
 380               anal->x[i][j][ZZ]/100);
 381       fprintf(kp,"\n");
 382     }
 383     fprintf(kp,"ENDMDL\n");
 384   }
 385   ffclose(kp);
 386 }
 387
 388 char *enms[eNR] = {
 389   "Coulomb", "Repulsion", "Potential",
 390   "EkHole",  "EkElectron", "EkLattice", "Kinetic",
 391   "Total"
 392 };
 393
 394 void add_ana_ener(t_ana_ener *ae,int nn,real e[])
 395 {
 396   int i;
 397
 398   /* First time around we are constantly increasing the array size */
 399   if (nn >= ae->nx) {
 400     if (ae->nx == ae->maxx) {
 401       ae->maxx += 1024;
 402       srenew(ae->e,ae->maxx);
 403     }
 404     for(i=0; (i<eNR); i++)
 405       ae->e[ae->nx][i] = e[i];
 406     ae->nx++;
 407   }
 408   else {
 409     for(i=0; (i<eNR); i++)
 410       ae->e[nn][i] += e[i];
 411   }
 412 }
 413
 414 void dump_ana_ener(t_ana_ener *ae,int nsim,real dt,char *edump,
 415                    t_ana_struct *total)
 416 {
 417   FILE *fp;
 418   int  i,j;
 419   real fac;
 420
 421   fac = 1.0/(nsim*ELECTRONVOLT);
 422   fp=xvgropen(edump,"Energies","Time (fs)","E (eV)");
 423   xvgr_legend(fp,eNR,enms);
 424   fprintf(fp,"@ s%d legend \"Ek/Nelec\"\n",eNR);
 425   fprintf(fp,"@ type nxy\n");
 426   for(i=0; (i<ae->nx); i++) {
 427     fprintf(fp,"%10f",1000.0*dt*i);
 428     for(j=0; (j<eNR); j++)
 429       fprintf(fp,"  %8.3f",ae->e[i][j]*fac);
 430     fprintf(fp,"  %8.3f\n",ae->e[i][eELECTRON]/(ELECTRONVOLT*total->nion[i]));
 431   }
 432   fprintf(fp,"&\n");
 433   ffclose(fp);
 434 }
 435