src/gmxlib/shift_util.c

   1 /*
   2  * $Id$
   3  *
   4  *       This source code is part of
   5  *
   6  *        G   R   O   M   A   C   S
   7  *
   8  * GROningen MAchine for Chemical Simulations
   9  *
  10  *               VERSION 2.0
  11  *
  12  * Copyright (c) 1991-1999
  13  * BIOSON Research Institute, Dept. of Biophysical Chemistry
  14  * University of Groningen, The Netherlands
  15  *
  16  * Please refer to:
  17  * GROMACS: A message-passing parallel molecular dynamics implementation
  18  * H.J.C. Berendsen, D. van der Spoel and R. van Drunen
  19  * Comp. Phys. Comm. 91, 43-56 (1995)
  20  *
  21  * Also check out our WWW page:
  22  * http://md.chem.rug.nl/~gmx
  23  * or e-mail to:
  24  * gromacs@chem.rug.nl
  25  *
  26  * And Hey:
  27  * Green Red Orange Magenta Azure Cyan Skyblue
  28  */
  29 static char *SRCID_shift_util_c = "$Id$";
  30
  31 #include <stdio.h>
  32 #include <math.h>
  33 #include "typedefs.h"
  34 #include "vec.h"
  35 #include "shift_util.h"
  36 #include "smalloc.h"
  37 #include "physics.h"
  38 #include "txtdump.h"
  39 #include "futil.h"
  40 #include "names.h"
  41 #include "fftgrid.h"
  42 #include "writeps.h"
  43 #include "macros.h"
  44 #include "xvgr.h"
  45
  46 #define p2(x) ((x)*(x))
  47 #define p3(x) ((x)*(x)*(x))
  48 #define p4(x) ((x)*(x)*(x)*(x))
  49
  50 static real A,A_3,B,B_4,C,c1,c2,c3,c4,c5,c6,One_4pi,FourPi_V,Vol,N0;
  51
  52 void set_shift_consts(FILE *log,real r1,real rc,rvec box,t_forcerec *fr)
  53 {
  54   /* A, B and C are recalculated in tables.c */
  55   if (r1 < rc) {
  56     A   = (2*r1-5*rc)/(p3(rc)*p2(rc-r1));
  57     B   = (4*rc-2*r1)/(p3(rc)*p3(rc-r1));
  58     /*C   = (10*rc*rc-5*rc*r1+r1*r1)/(6*rc*rc); Hermans Eq. not correct */
  59   }
  60   else
  61     fatal_error(0,"r1 (%f) >= rc (%f) in %s, line %d",
  62                 r1,rc,__FILE__,__LINE__);
  63
  64   A_3 = A/3.0;
  65   B_4 = B/4.0;
  66   C   = 1/rc-A_3*p3(rc-r1)-B_4*p4(rc-r1);
  67   N0  = 2.0*M_PI*p3(rc)*p3(rc-r1);
  68
  69   Vol     =(box[XX]*box[YY]*box[ZZ]);
  70   FourPi_V=4.0*M_PI/Vol;
  71
  72   if(debug)
  73       fprintf(log,"Constants for short-range and fourier stuff:\n"
  74               "r1 = %10.3f,  rc = %10.3f\n"
  75               "A  = %10.3e,  B  = %10.3e,  C  = %10.3e, FourPi_V = %10.3e\n",
  76               r1,rc,A,B,C,FourPi_V);
  77
  78   /* Constants derived by Mathematica */
  79   c1 = -40*rc*rc    + 50*rc*r1    - 16*r1*r1;
  80   c2 =  60*rc       - 30*r1;
  81   c3 = -10*rc*rc*rc + 20*rc*rc*r1 - 13*rc*r1*r1 + 3*r1*r1*r1;
  82   c4 = -20*rc*rc    + 40*rc*r1    - 14*r1*r1;
  83   c5 = -c2;
  84   c6 = -5*rc*rc*r1  +  7*rc*r1*r1 - 2*r1*r1*r1;
  85
  86   if(debug)
  87       fprintf(log,"c1 = %10.3e,  c2 = %10.3e,  c3 = %10.3e\n"
  88               "c4 = %10.3e,  c5 = %10.3e,  c6 = %10.3e,  N0 = %10.3e\n",
  89               c1,c2,c3,c4,c5,c6,N0);
  90
  91   One_4pi = 1.0/(4.0*M_PI);
  92 }
  93
  94 real gk(real k,real rc,real r1)
  95 /* Spread function in Fourier space. Eqs. 56-64 */
  96 {
  97   real N,gg;
  98
  99   N  = N0*p4(k);
 100
 101   /* c1 thru c6 consts are global variables! */
 102   gg = (FourPi_V / (N*k)) *
 103     ( c1*k*cos(k*rc) + (c2+c3*k*k)*sin(k*rc) +
 104       c4*k*cos(k*r1) + (c5+c6*k*k)*sin(k*r1) );
 105
 106   return gg;
 107 }
 108
 109 real gknew(real k,real rc,real r1)
 110 {
 111   real rck,rck2;
 112
 113   rck = rc*k;
 114   rck2= rck*rck;
 115
 116   return -15.0*((rck2-3.0)*sin(rck) + 3*rck*cos(rck))/(Vol*rck2*rck2*rck);
 117 }
 118
 119 real calc_dx2dx(rvec xi,rvec xj,rvec box,rvec dx)
 120 {
 121   int  m;
 122   real ddx,dx2;
 123
 124   dx2=0;
 125   for(m=0; (m<DIM); m++) {
 126     ddx=xj[m]-xi[m];
 127     if (ddx < -0.5*box[m])
 128       ddx+=box[m];
 129     else if (ddx >= 0.5*box[m])
 130       ddx-=box[m];
 131     dx[m]=ddx;
 132     dx2 += ddx*ddx;
 133   }
 134   return dx2;
 135 }
 136
 137 real calc_dx2(rvec xi,rvec xj,rvec box)
 138 {
 139   rvec dx;
 140
 141   return calc_dx2dx(xi,xj,box,dx);
 142 }
 143
 144 real shiftfunction(real r1,real rc,real R)
 145 {
 146   real dr;
 147
 148   if (R <= r1)
 149     return 0.0;
 150   else if (R >= rc)
 151     return -1.0/(R*R);
 152
 153   dr=R-r1;
 154
 155   return A*dr*dr+B*dr*dr*dr;
 156 }
 157
 158 real new_f(real r,real rc)
 159 {
 160   real rrc,rrc2,rrc3;
 161
 162   rrc  = r/rc;
 163   rrc2 = rrc*rrc;
 164   rrc3 = rrc2*rrc;
 165   return 1.5*rrc2*rrc3 - 2.5*rrc3 + 1.0;
 166 }
 167
 168 real new_phi(real r,real rc)
 169 {
 170   real rrr;
 171
 172   rrr = sqr(r/rc);
 173
 174   return 1/r-(0.125/rc)*(15 + 3*rrr*rrr - 10*rrr);
 175 }
 176
 177 real old_f(real r,real rc,real r1)
 178 {
 179   real dr,r2;
 180
 181   if (r <= r1)
 182     return 1.0;
 183   else if (r >= rc)
 184     return 0;
 185
 186   dr = r-r1;
 187   r2 = r*r;
 188   return 1+A*r2*dr*dr+B*r2*dr*dr*dr;
 189 }
 190
 191 real old_phi(real r,real rc,real r1)
 192 {
 193   real dr;
 194
 195   if (r <= r1)
 196     return 1/r-C;
 197   else if (r >= rc)
 198     return 0.0;
 199
 200   dr = r-r1;
 201
 202   return 1/r-A_3*dr*dr*dr-B_4*dr*dr*dr*dr - C;
 203 }
 204
 205 real phi_sr(FILE *log,int nj,rvec x[],real charge[],real rc,real r1,rvec box,
 206             real phi[],t_block *excl,rvec f_sr[],bool bOld)
 207 {
 208   int  i,j,k,m,ni,i1,i2;
 209   real pp,r2,R,R_1,R_2,rc2;
 210   real qi,qj,vsr,eps,fscal;
 211   rvec dx;
 212
 213   vsr = 0.0;
 214   eps = ONE_4PI_EPS0;
 215   rc2 = rc*rc;
 216   ni=0;
 217   for(i=0; (i<nj-1); i++) {
 218     qi=charge[i];
 219     for(j=i+1; (j<nj); j++) {
 220       i1=excl->index[i];
 221       i2=excl->index[i+1];
 222       for(k=i1; (k<i2); k++)
 223         if (excl->a[k] == j)
 224           break;
 225       if (k == i2) {
 226         r2=calc_dx2dx(x[i],x[j],box,dx);
 227         if (r2 < rc2) {
 228           qj    = charge[j];
 229           R_1   = invsqrt(r2);
 230           R_2   = R_1*R_1;
 231           R     = invsqrt(R_2);
 232           if (bOld) {
 233             fscal = old_f(R,rc,r1)*R_2;
 234             pp    = old_phi(R,rc,r1);
 235           }
 236           else {
 237             fscal = new_f(R,rc)*R_2;
 238             pp    = new_phi(R,rc);
 239           }
 240           phi[i] += eps*qj*pp;
 241           phi[j] += eps*qi*pp;
 242           vsr    += eps*qj*qi*pp;
 243           for(m=0; (m<DIM); m++) {
 244             f_sr[i][m] += dx[m]*fscal;
 245             f_sr[j][m] -= dx[m]*fscal;
 246           }
 247           ni++;
 248         }
 249       }
 250     }
 251   }
 252   fprintf(log,"There were %d short range interactions, vsr=%g\n",ni,vsr);
 253
 254   return vsr;
 255 }
 256
 257 real spreadfunction(real r1,real rc,real R)
 258 {
 259   real dr;
 260
 261   if (R <= r1)
 262     return 0.0;
 263   else if (R >= rc)
 264     return 0.0;
 265
 266   dr=R-r1;
 267
 268   return -One_4pi*(dr/R)*(2*A*(2*R-r1)+B*dr*(5*R-2*r1));
 269 }
 270
 271 real potential(real r1,real rc,real R)
 272 {
 273   if (R < r1)
 274     return (1.0/R-C);
 275   else if (R <= rc)
 276     return (1.0/R - A_3*p3(R-r1) - B_4*p4(R-r1) - C);
 277   else
 278     return 0.0;
 279 }
 280
 281
 282
 283 real shift_LRcorrection(FILE *fp,t_nsborder *nsb,t_commrec *cr,t_forcerec *fr,
 284                         real charge[],t_block *excl,rvec x[],
 285                         bool bOld,rvec box_size,matrix lr_vir)
 286 {
 287   static bool bFirst=TRUE;
 288   static real Vself;
 289   int    i,i1,i2,j,k,m,iv,jv;
 290   unsigned int *AA;
 291   double qq; /* Necessary for precision */
 292   double isp=0.564189583547756;
 293   real   qi,dr,ddd,dr2,dr_1,dr_3,fscal,Vexcl,qtot=0;
 294   rvec   df,dx;
 295   real   r1=fr->rcoulomb_switch;
 296   real   rc=fr->rcoulomb;
 297   ivec   shift;
 298   int    start=START(nsb);
 299   int    natoms=HOMENR(nsb);
 300
 301   if (bFirst) {
 302     qq =0;
 303     for(i=start; (i<start+natoms); i++)
 304       qq  += charge[i]*charge[i];
 305
 306     /* Obsolete, until we implement dipole and charge corrections.
 307        qtot=0;
 308        for(i=0;i<nsb->natoms;i++)
 309        qtot+=charge[i];
 310     */
 311
 312     Vself = 0.5*C*ONE_4PI_EPS0*qq;
 313     if(debug) {
 314         fprintf(fp,"Long Range corrections for shifted interactions:\n");
 315         fprintf(fp,"r1 = %g, rc=%g\n",r1,rc);
 316         fprintf(fp,"start=%d,natoms=%d\n",start,natoms);
 317         fprintf(fp,"qq = %g, Vself=%g\n",qq,Vself);
 318     }
 319
 320   }
 321   AA = excl->a;
 322   Vexcl = 0;
 323
 324   for(i=start; (i<start+natoms); i++) {
 325     /* Initiate local variables (for this i-particle) to 0 */
 326     i1  = excl->index[i];
 327     i2  = excl->index[i+1];
 328     qi  = charge[i]*ONE_4PI_EPS0;
 329
 330     /* Loop over excluded neighbours */
 331     for(j=i1; (j<i2); j++) {
 332       k = AA[j];
 333       /*
 334        * First we must test whether k <> i, and then, because the
 335        * exclusions are all listed twice i->k and k->i we must select
 336        * just one of the two.
 337        * As a minor optimization we only compute forces when the charges
 338        * are non-zero.
 339        */
 340       if (k > i) {
 341         qq = qi*charge[k];
 342         if (qq != 0.0) {
 343           /* Compute distance vector, no PBC check! */
 344           dr2 = 0;
 345           for(m=0; (m<DIM); m++) {
 346             ddd = x[i][m] - x[k][m];
 347               if(ddd>box_size[m]/2) {  /* ugly hack,   */
 348                 ddd-=box_size[m];      /* to fix pbc.. */
 349               } else if (ddd<-box_size[m]/2)
 350                   ddd+=box_size[m];
 351
 352             dx[m] = ddd;
 353             dr2  += ddd*ddd;
 354           }
 355           dr_1    = invsqrt(dr2);
 356           dr      = 1.0/dr_1;
 357           dr_3    = dr_1*dr_1*dr_1;
 358           /* Compute exclusion energy and scalar force */
 359
 360           Vexcl  += qq*(dr_1-potential(r1,rc,dr));
 361           fscal   = qq*(-shiftfunction(r1,rc,dr))*dr_3;
 362
 363           if ((fscal != 0) && debug)
 364             fprintf(debug,"i: %d, k: %d, dr: %.3f fscal: %.3f\n",i,k,dr,fscal);
 365
 366           /* The force vector is obtained by multiplication with the
 367            * distance vector
 368            */
 369           svmul(fscal,dx,df);
 370           rvec_inc(fr->flr[k],df);
 371           rvec_dec(fr->flr[i],df);
 372           for(iv=0;iv<DIM;iv++)
 373               for(jv=0;jv<DIM;jv++)
 374                   lr_vir[iv][jv]+=0.5*dx[iv]*df[jv];
 375         }
 376       }
 377     }
 378   }
 379   if (bFirst && debug)
 380     fprintf(fp,"Long Range correction: Vexcl=%g\n",Vexcl);
 381
 382   bFirst = FALSE;
 383
 384   return (Vself+Vexcl);
 385 }
 386
 387
 388 void calc_ener(FILE *fp,char *title,bool bHeader,int nmol,
 389                int natoms,real phi[],real charge[],t_block *excl)
 390 {
 391   int  i,i1,i2,j,k;
 392   real qq,qi,vv,V,Vex,Vc,Vt;
 393
 394   qq   = 0;
 395   V    = 0;
 396   for(i=0; (i<natoms); i++) {
 397     vv   = charge[i]*phi[i];
 398     V   += vv;
 399     qq  += charge[i]*charge[i];
 400   }
 401   V  = 0.5*V;
 402   Vc = 0.5*C*ONE_4PI_EPS0*qq;
 403
 404   qq = 0;
 405   for(i=0; (i<excl->nr); i++) {
 406     i1 = excl->index[i];
 407     i2 = excl->index[i+1];
 408     qi = charge[i];
 409     for(j=i1; (j<i2); j++) {
 410       k = excl->a[j];
 411       if (k != i)
 412         qq+=qi*charge[k];
 413     }
 414   }
 415   Vex = qq*0.5*C*ONE_4PI_EPS0;
 416
 417   Vt = V - Vc - Vex;
 418
 419   if (bHeader) {
 420     fprintf(fp,"%12s  %12s  %12s  %12s  %12s  %12s\n",
 421             "","Vphi","Vself","Vexcl","Vtot","1Mol");
 422
 423   }
 424   fprintf(fp,"%12s  %12.5e  %12.5e  %12.5e  %12.5e  %12.5e\n",
 425           title,V,Vc,Vex,Vt,Vt/nmol);
 426 }
 427
 428 real phi_aver(int natoms,real phi[])
 429 {
 430   real phitot;
 431   int  i;
 432
 433   phitot=0;
 434   for(i=0; (i<natoms); i++)
 435     phitot=phitot+phi[i];
 436
 437   return (phitot/natoms);
 438 }
 439
 440 real symmetrize_phi(FILE *log,int natoms,real phi[],bool bVerbose)
 441 {
 442   real phitot;
 443   int  i;
 444
 445   phitot=phi_aver(natoms,phi);
 446   if (bVerbose)
 447     fprintf(log,"phi_aver = %10.3e\n",phitot);
 448
 449   for(i=0; (i<natoms); i++)
 450     phi[i]-=phitot;
 451
 452   return phitot;
 453 }
 454
 455 static real rgbset(real col)
 456 {
 457   int icol32;
 458
 459   icol32=32.0*col;
 460   return icol32/32.0;
 461 }
 462
 463 void plot_phi(char *fn,rvec box,int natoms,rvec x[],real phi[])
 464 {
 465   FILE *eps;
 466   real phi_max,rr,gg,bb,fac,dx,x0,y0;
 467   real offset;
 468   int  i;
 469
 470   phi_max=phi[0];
 471   rr=gg=bb=0.0;
 472   for(i=0; (i<natoms); i++)
 473     phi_max=max(phi_max,fabs(phi[i]));
 474
 475   if (phi_max==0.0) {
 476     fprintf(stderr,"All values zero, see .out file\n");
 477     return;
 478   }
 479   offset=20.0;
 480   fac=15.0;
 481 #ifdef DEBUG
 482   fprintf(stderr,"Scaling box by %g\n",fac);
 483 #endif
 484   eps=ps_open(fn,0,0,
 485               (real)(fac*box[XX]+2*offset),(real)(fac*box[YY]+2*offset));
 486   ps_translate(eps,offset,offset);
 487   ps_color(eps,0,0,0);
 488   ps_box(eps,1,1,(real)(fac*box[XX]-1),(real)(fac*box[YY]-1));
 489   dx=0.15*fac;
 490   for(i=0; (i<natoms); i++) {
 491     rr=gg=bb=1.0;
 492     if (phi[i] < 0)
 493       gg=bb=(1.0+(phi[i]/phi_max));
 494     else
 495       rr=gg=(1.0-(phi[i]/phi_max));
 496     rr=rgbset(rr);
 497     gg=rgbset(gg);
 498     bb=rgbset(bb);
 499     ps_color(eps,rr,gg,bb);
 500     x0=fac*x[i][XX];
 501     y0=fac*x[i][YY];
 502     ps_fillbox(eps,(real)(x0-dx),(real)(y0-dx),(real)(x0+dx),(real)(y0+dx));
 503   }
 504   ps_close(eps);
 505 }
 506
 507 void plot_qtab(char *fn,int nx,int ny,int nz,real ***qtab)
 508 {
 509   rvec box;
 510   rvec *xx;
 511   real *phi;
 512   int  i,npt,ix,iy,iz;
 513
 514   box[XX]=nx;
 515   box[YY]=ny;
 516   box[ZZ]=nz;
 517
 518   npt=(box[XX]*box[YY]*box[ZZ]);
 519   snew(xx,npt);
 520   snew(phi,npt);
 521   nx/=2;
 522   ny/=2;
 523   nz/=2;
 524   i=0;
 525   for(ix=-nx; (ix<nx); ix++)
 526     for(iy=-ny; (iy<ny); iy++)
 527       for(iz=-nz; (iz<nz); iz++,i++) {
 528         xx[i][XX]=ix+nx+0.5;
 529         xx[i][YY]=iy+ny+0.5;
 530         xx[i][ZZ]=iz+nz+0.5; /* onzin */
 531         phi[i]=qtab[ix+nx][iy+ny][iz+nz];
 532       }
 533
 534   plot_phi(fn,box,npt,xx,phi);
 535
 536   sfree(xx);
 537   sfree(phi);
 538 }
 539
 540 void print_phi(char *fn,int natoms,rvec x[],real phi[])
 541 {
 542   FILE *fp;
 543   int  i;
 544
 545   fp=ffopen(fn,"w");
 546   for(i=0; (i<natoms); i++)
 547     fprintf(fp,"%10d  %12.5e\n",i,phi[i]);
 548   fclose(fp);
 549 }
 550
 551 void write_pqr(char *fn,t_atoms *atoms,rvec x[],real phi[],real dx)
 552 {
 553   FILE *fp;
 554   int  i,rnr;
 555
 556   fp=ffopen(fn,"w");
 557   for(i=0; (i<atoms->nr); i++) {
 558     rnr=atoms->atom[i].resnr;
 559     fprintf(fp,"%-6s%5d  %-4.4s%3.3s %c%4d    %8.3f%8.3f%8.3f%6.2f%6.2f\n",
 560             "ATOM",(i+1),*atoms->atomname[i],*atoms->resname[rnr],' ',
 561             (rnr+1) % 10000,
 562             10*(dx+x[i][XX]),10*x[i][YY],10*(x[i][ZZ]),0.0,phi[i]);
 563   }
 564   fclose(fp);
 565 }
 566
 567 void write_grid_pqr(char *fn,int nx,int ny,int nz,real ***phi)
 568 {
 569   FILE *fp;
 570   int  i,j,k,rnr=0;
 571   real fac=4.0;
 572
 573   fp=ffopen(fn,"w");
 574   for(i=0; (i<nx); i++)
 575     for(j=0; (j<ny); j++)
 576       for(k=0; (k<nz); k++,rnr++)
 577         fprintf(fp,"%-6s%5d  %-4.4s%3.3s %c%4d    %8.3f%8.3f%8.3f%6.2f%6.2f\n",
 578                 "ATOM",(i+1),"C","C",' ',
 579                 1+((rnr+1) % 10000),fac*i,fac*j,fac*k,0.0,phi[i][j][k]);
 580   fclose(fp);
 581 }
 582
 583
 584 real analyse_diff(FILE *log,char *label,
 585                   int natom,rvec ffour[],rvec fpppm[],
 586                   real phi_f[],real phi_p[],real phi_sr[],
 587                   char *fcorr,char *pcorr,char *ftotcorr,char *ptotcorr)
 588 {
 589   int  i,m;
 590   FILE *fp=NULL,*gp=NULL;
 591   real f2sum=0,p2sum=0;
 592   real df,fmax,dp,pmax,rmsf;
 593
 594   fmax = fabs(ffour[0][0]-fpppm[0][0]);
 595   pmax = fabs(phi_f[0] - phi_p[0]);
 596
 597   for(i=0; (i<natom); i++) {
 598     dp     = fabs(phi_f[i] - phi_p[i]);
 599     pmax   = max(dp,pmax);
 600     p2sum += dp*dp;
 601     for(m=0; (m<DIM); m++) {
 602       df     = fabs(ffour[i][m] - fpppm[i][m]);
 603       fmax   = max(df,fmax);
 604       f2sum += df*df;
 605     }
 606   }
 607
 608   rmsf = sqrt(f2sum/(3.0*natom));
 609   fprintf(log,"\n********************************\nERROR ANALYSIS for %s\n",
 610           label);
 611   fprintf(log,"%-10s%12s%12s\n","Error:","Max Abs","RMS");
 612   fprintf(log,"%-10s  %10.3f  %10.3f\n","Force",fmax,rmsf);
 613   fprintf(log,"%-10s  %10.3f  %10.3f\n","Potential",pmax,sqrt(p2sum/(natom)));
 614
 615   if (fcorr) {
 616     fp = xvgropen(fcorr,"LR Force Correlation","Four-Force","PPPM-Force");
 617     for(i=0; (i<natom); i++) {
 618       for(m=0; (m<DIM); m++) {
 619         fprintf(fp,"%10.3f  %10.3f\n",ffour[i][m],fpppm[i][m]);
 620       }
 621     }
 622     ffclose(fp);
 623     xvgr_file(fcorr,NULL);
 624   }
 625   if (pcorr)
 626     fp = xvgropen(pcorr,"LR Potential Correlation","Four-Pot","PPPM-Pot");
 627   if (ptotcorr)
 628     gp = xvgropen(ptotcorr,"Total Potential Correlation","Four-Pot","PPPM-Pot");
 629   for(i=0; (i<natom); i++) {
 630     if (pcorr)
 631       fprintf(fp,"%10.3f  %10.3f\n",phi_f[i],phi_p[i]);
 632     if (ptotcorr)
 633       fprintf(gp,"%10.3f  %10.3f\n",phi_f[i]+phi_sr[i],phi_p[i]+phi_sr[i]);
 634   }
 635   if (pcorr) {
 636     ffclose(fp);
 637     xvgr_file(pcorr,NULL);
 638   }
 639   if (ptotcorr) {
 640     ffclose(gp);
 641     xvgr_file(ptotcorr,NULL);
 642   }
 643
 644   return rmsf;
 645 }
 646