src/gmxlib/shift_util.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  * GROningen Mixture of Alchemy and Childrens' Stories
  34  */
  35 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  38
  39 #include <stdio.h>
  40 #include <math.h>
  41 #include "typedefs.h"
  42 #include "vec.h"
  43 #include "coulomb.h"
  44 #include "smalloc.h"
  45 #include "physics.h"
  46 #include "txtdump.h"
  47 #include "futil.h"
  48 #include "names.h"
  49 #include "writeps.h"
  50 #include "macros.h"
  51 #include "xvgr.h"
  52 #include "pppm.h"
  53 #include "gmxfio.h"
  54
  55 #ifdef GMX_THREADS
  56 #include "thread_mpi.h"
  57 #endif
  58
  59 #define p2(x) ((x)*(x))
  60 #define p3(x) ((x)*(x)*(x))
  61 #define p4(x) ((x)*(x)*(x)*(x))
  62
  63 static real A,A_3,B,B_4,C,c1,c2,c3,c4,c5,c6,One_4pi,FourPi_V,Vol,N0;
  64 #ifdef GMX_THREADS
  65 static tMPI_Thread_mutex_t shift_mutex=TMPI_THREAD_MUTEX_INITIALIZER;
  66 #endif
  67
  68
  69 void set_shift_consts(FILE *log,real r1,real rc,rvec box,t_forcerec *fr)
  70 {
  71 #ifdef GMX_THREADS
  72   /* at the very least we shouldn't allow multiple threads to set these
  73      simulataneously */
  74   tMPI_Thread_mutex_lock(&shift_mutex);
  75 #endif
  76   /* A, B and C are recalculated in tables.c */
  77   if (r1 < rc) {
  78     A   = (2*r1-5*rc)/(p3(rc)*p2(rc-r1));
  79     B   = (4*rc-2*r1)/(p3(rc)*p3(rc-r1));
  80     /*C   = (10*rc*rc-5*rc*r1+r1*r1)/(6*rc*rc); Hermans Eq. not correct */
  81   }
  82   else
  83     gmx_fatal(FARGS,"r1 (%f) >= rc (%f) in %s, line %d",
  84               r1,rc,__FILE__,__LINE__);
  85
  86   A_3 = A/3.0;
  87   B_4 = B/4.0;
  88   C   = 1/rc-A_3*p3(rc-r1)-B_4*p4(rc-r1);
  89   N0  = 2.0*M_PI*p3(rc)*p3(rc-r1);
  90
  91   Vol     =(box[XX]*box[YY]*box[ZZ]);
  92   FourPi_V=4.0*M_PI/Vol;
  93
  94   if (debug) {
  95       fprintf(debug,"Constants for short-range and fourier stuff:\n"
  96               "r1 = %10.3f,  rc = %10.3f\n"
  97               "A  = %10.3e,  B  = %10.3e,  C  = %10.3e, FourPi_V = %10.3e\n",
  98               r1,rc,A,B,C,FourPi_V);
  99   }
 100
 101   /* Constants derived by Mathematica */
 102   c1 = -40*rc*rc    + 50*rc*r1    - 16*r1*r1;
 103   c2 =  60*rc       - 30*r1;
 104   c3 = -10*rc*rc*rc + 20*rc*rc*r1 - 13*rc*r1*r1 + 3*r1*r1*r1;
 105   c4 = -20*rc*rc    + 40*rc*r1    - 14*r1*r1;
 106   c5 = -c2;
 107   c6 = -5*rc*rc*r1  +  7*rc*r1*r1 - 2*r1*r1*r1;
 108
 109   if (debug) {
 110       fprintf(debug,"c1 = %10.3e,  c2 = %10.3e,  c3 = %10.3e\n"
 111               "c4 = %10.3e,  c5 = %10.3e,  c6 = %10.3e,  N0 = %10.3e\n",
 112               c1,c2,c3,c4,c5,c6,N0);
 113   }
 114
 115   One_4pi = 1.0/(4.0*M_PI);
 116 #ifdef GMX_THREADS
 117   tMPI_Thread_mutex_unlock(&shift_mutex);
 118 #endif
 119 }
 120
 121 real gk(real k,real rc,real r1)
 122 /* Spread function in Fourier space. Eqs. 56-64 */
 123 {
 124   real N,gg;
 125
 126   N  = N0*p4(k);
 127
 128   /* c1 thru c6 consts are global variables! */
 129   gg = (FourPi_V / (N*k)) *
 130     ( c1*k*cos(k*rc) + (c2+c3*k*k)*sin(k*rc) +
 131       c4*k*cos(k*r1) + (c5+c6*k*k)*sin(k*r1) );
 132
 133   return gg;
 134 }
 135
 136 real gknew(real k,real rc,real r1)
 137 {
 138   real rck,rck2;
 139
 140   rck = rc*k;
 141   rck2= rck*rck;
 142
 143   return -15.0*((rck2-3.0)*sin(rck) + 3*rck*cos(rck))/(Vol*rck2*rck2*rck);
 144 }
 145
 146 real calc_dx2dx(rvec xi,rvec xj,rvec box,rvec dx)
 147 {
 148   int  m;
 149   real ddx,dx2;
 150
 151   dx2=0;
 152   for(m=0; (m<DIM); m++) {
 153     ddx=xj[m]-xi[m];
 154     if (ddx < -0.5*box[m])
 155       ddx+=box[m];
 156     else if (ddx >= 0.5*box[m])
 157       ddx-=box[m];
 158     dx[m]=ddx;
 159     dx2 += ddx*ddx;
 160   }
 161   return dx2;
 162 }
 163
 164 real calc_dx2(rvec xi,rvec xj,rvec box)
 165 {
 166   rvec dx;
 167
 168   return calc_dx2dx(xi,xj,box,dx);
 169 }
 170
 171 real shiftfunction(real r1,real rc,real R)
 172 {
 173   real dr;
 174
 175   if (R <= r1)
 176     return 0.0;
 177   else if (R >= rc)
 178     return -1.0/(R*R);
 179
 180   dr=R-r1;
 181
 182   return A*dr*dr+B*dr*dr*dr;
 183 }
 184
 185 real new_f(real r,real rc)
 186 {
 187   real rrc,rrc2,rrc3;
 188
 189   rrc  = r/rc;
 190   rrc2 = rrc*rrc;
 191   rrc3 = rrc2*rrc;
 192   return 1.5*rrc2*rrc3 - 2.5*rrc3 + 1.0;
 193 }
 194
 195 real new_phi(real r,real rc)
 196 {
 197   real rrr;
 198
 199   rrr = sqr(r/rc);
 200
 201   return 1/r-(0.125/rc)*(15 + 3*rrr*rrr - 10*rrr);
 202 }
 203
 204 real old_f(real r,real rc,real r1)
 205 {
 206   real dr,r2;
 207
 208   if (r <= r1)
 209     return 1.0;
 210   else if (r >= rc)
 211     return 0;
 212
 213   dr = r-r1;
 214   r2 = r*r;
 215   return 1+A*r2*dr*dr+B*r2*dr*dr*dr;
 216 }
 217
 218 real old_phi(real r,real rc,real r1)
 219 {
 220   real dr;
 221
 222   if (r <= r1)
 223     return 1/r-C;
 224   else if (r >= rc)
 225     return 0.0;
 226
 227   dr = r-r1;
 228
 229   return 1/r-A_3*dr*dr*dr-B_4*dr*dr*dr*dr - C;
 230 }
 231
 232 real spreadfunction(real r1,real rc,real R)
 233 {
 234   real dr;
 235
 236   if (R <= r1)
 237     return 0.0;
 238   else if (R >= rc)
 239     return 0.0;
 240
 241   dr=R-r1;
 242
 243   return -One_4pi*(dr/R)*(2*A*(2*R-r1)+B*dr*(5*R-2*r1));
 244 }
 245
 246 real potential(real r1,real rc,real R)
 247 {
 248   if (R < r1)
 249     return (1.0/R-C);
 250   else if (R <= rc)
 251     return (1.0/R - A_3*p3(R-r1) - B_4*p4(R-r1) - C);
 252   else
 253     return 0.0;
 254 }
 255
 256
 257
 258 real shift_LRcorrection(FILE *fp,int start,int natoms,
 259                         t_commrec *cr,t_forcerec *fr,
 260                         real charge[],t_blocka *excl,rvec x[],
 261                         gmx_bool bOld,matrix box,matrix lr_vir)
 262 {
 263   static gmx_bool bFirst=TRUE;
 264   static real Vself;
 265   int    i,i1,i2,j,k,m,iv,jv;
 266   int *AA;
 267   double qq; /* Necessary for precision */
 268   double isp=0.564189583547756;
 269   real   qi,dr,dr2,dr_1,dr_3,fscal,Vexcl,qtot=0;
 270   rvec   df,dx;
 271   real   r1=fr->rcoulomb_switch;
 272   real   rc=fr->rcoulomb;
 273   ivec   shift;
 274
 275   if (bFirst) {
 276     qq =0;
 277     for(i=start; (i<start+natoms); i++)
 278       qq  += charge[i]*charge[i];
 279
 280     /* Obsolete, until we implement dipole and charge corrections.
 281        qtot=0;
 282        for(i=0;i<nsb->natoms;i++)
 283        qtot+=charge[i];
 284     */
 285
 286     Vself = 0.5*C*ONE_4PI_EPS0*qq;
 287     if(debug) {
 288         fprintf(fp,"Long Range corrections for shifted interactions:\n");
 289         fprintf(fp,"r1 = %g, rc=%g\n",r1,rc);
 290         fprintf(fp,"start=%d,natoms=%d\n",start,natoms);
 291         fprintf(fp,"qq = %g, Vself=%g\n",qq,Vself);
 292     }
 293
 294   }
 295   AA = excl->a;
 296   Vexcl = 0;
 297
 298   for(i=start; (i<start+natoms); i++) {
 299     /* Initiate local variables (for this i-particle) to 0 */
 300     i1  = excl->index[i];
 301     i2  = excl->index[i+1];
 302     qi  = charge[i]*ONE_4PI_EPS0;
 303
 304     /* Loop over excluded neighbours */
 305     for(j=i1; (j<i2); j++) {
 306       k = AA[j];
 307       /*
 308        * First we must test whether k <> i, and then, because the
 309        * exclusions are all listed twice i->k and k->i we must select
 310        * just one of the two.
 311        * As a minor optimization we only compute forces when the charges
 312        * are non-zero.
 313        */
 314       if (k > i) {
 315         qq = qi*charge[k];
 316         if (qq != 0.0) {
 317           dr2 = 0;
 318           rvec_sub(x[i],x[k],dx);
 319           for(m=DIM-1; m>=0; m--) {
 320             if (dx[m] > 0.5*box[m][m])
 321               rvec_dec(dx,box[m]);
 322             else if (dx[m] < -0.5*box[m][m])
 323               rvec_inc(dx,box[m]);
 324
 325             dr2  += dx[m]*dx[m];
 326           }
 327           dr_1    = gmx_invsqrt(dr2);
 328           dr      = 1.0/dr_1;
 329           dr_3    = dr_1*dr_1*dr_1;
 330           /* Compute exclusion energy and scalar force */
 331
 332           Vexcl  += qq*(dr_1-potential(r1,rc,dr));
 333           fscal   = qq*(-shiftfunction(r1,rc,dr))*dr_3;
 334
 335           if ((fscal != 0) && debug)
 336             fprintf(debug,"i: %d, k: %d, dr: %.3f fscal: %.3f\n",i,k,dr,fscal);
 337
 338           /* The force vector is obtained by multiplication with the
 339            * distance vector
 340            */
 341           svmul(fscal,dx,df);
 342           rvec_inc(fr->f_novirsum[k],df);
 343           rvec_dec(fr->f_novirsum[i],df);
 344           for(iv=0;iv<DIM;iv++)
 345               for(jv=0;jv<DIM;jv++)
 346                   lr_vir[iv][jv]+=0.5*dx[iv]*df[jv];
 347         }
 348       }
 349     }
 350   }
 351   if (bFirst && debug)
 352     fprintf(fp,"Long Range correction: Vexcl=%g\n",Vexcl);
 353
 354   bFirst = FALSE;
 355   /* Return the correction to the energy */
 356   return (-(Vself+Vexcl));
 357 }
 358
 359 real phi_aver(int natoms,real phi[])
 360 {
 361   real phitot;
 362   int  i;
 363
 364   phitot=0;
 365   for(i=0; (i<natoms); i++)
 366     phitot=phitot+phi[i];
 367
 368   return (phitot/natoms);
 369 }
 370
 371 real symmetrize_phi(FILE *log,int natoms,real phi[],gmx_bool bVerbose)
 372 {
 373   real phitot;
 374   int  i;
 375
 376   phitot=phi_aver(natoms,phi);
 377   if (bVerbose)
 378     fprintf(log,"phi_aver = %10.3e\n",phitot);
 379
 380   for(i=0; (i<natoms); i++)
 381     phi[i]-=phitot;
 382
 383   return phitot;
 384 }
 385
 386 static real rgbset(real col)
 387 {
 388   int icol32;
 389
 390   icol32=32.0*col;
 391   return icol32/32.0;
 392 }
 393
 394
 395
 396 real analyse_diff(FILE *log,char *label,const output_env_t oenv,
 397                   int natom,rvec ffour[],rvec fpppm[],
 398                   real phi_f[],real phi_p[],real phi_sr[],
 399                   char *fcorr,char *pcorr,char *ftotcorr,char *ptotcorr)
 400 /* Analyse difference between forces from fourier (_f) and other (_p)
 401  * LR solvers (and potential also).
 402  * If the filenames are given, xvgr files are written.
 403  * returns the root mean square error in the force.
 404  */
 405 {
 406   int  i,m;
 407   FILE *fp=NULL,*gp=NULL;
 408   real f2sum=0,p2sum=0;
 409   real df,fmax,dp,pmax,rmsf;
 410
 411   fmax = fabs(ffour[0][0]-fpppm[0][0]);
 412   pmax = fabs(phi_f[0] - phi_p[0]);
 413
 414   for(i=0; (i<natom); i++) {
 415     dp     = fabs(phi_f[i] - phi_p[i]);
 416     pmax   = max(dp,pmax);
 417     p2sum += dp*dp;
 418     for(m=0; (m<DIM); m++) {
 419       df     = fabs(ffour[i][m] - fpppm[i][m]);
 420       fmax   = max(df,fmax);
 421       f2sum += df*df;
 422     }
 423   }
 424
 425   rmsf = sqrt(f2sum/(3.0*natom));
 426   fprintf(log,"\n********************************\nERROR ANALYSIS for %s\n",
 427           label);
 428   fprintf(log,"%-10s%12s%12s\n","Error:","Max Abs","RMS");
 429   fprintf(log,"%-10s  %10.3f  %10.3f\n","Force",fmax,rmsf);
 430   fprintf(log,"%-10s  %10.3f  %10.3f\n","Potential",pmax,sqrt(p2sum/(natom)));
 431
 432   if (fcorr) {
 433     fp = xvgropen(fcorr,"LR Force Correlation","Four-Force","PPPM-Force",oenv);
 434     for(i=0; (i<natom); i++) {
 435       for(m=0; (m<DIM); m++) {
 436         fprintf(fp,"%10.3f  %10.3f\n",ffour[i][m],fpppm[i][m]);
 437       }
 438     }
 439     gmx_fio_fclose(fp);
 440     do_view(oenv,fcorr,NULL);
 441   }
 442   if (pcorr)
 443     fp = xvgropen(pcorr,"LR Potential Correlation","Four-Pot","PPPM-Pot",oenv);
 444   if (ptotcorr)
 445     gp = xvgropen(ptotcorr,"Total Potential Correlation","Four-Pot","PPPM-Pot",
 446                   oenv);
 447   for(i=0; (i<natom); i++) {
 448     if (pcorr)
 449       fprintf(fp,"%10.3f  %10.3f\n",phi_f[i],phi_p[i]);
 450     if (ptotcorr)
 451       fprintf(gp,"%10.3f  %10.3f\n",phi_f[i]+phi_sr[i],phi_p[i]+phi_sr[i]);
 452   }
 453   if (pcorr) {
 454     gmx_fio_fclose(fp);
 455     do_view(oenv,pcorr,NULL);
 456   }
 457   if (ptotcorr) {
 458     gmx_fio_fclose(gp);
 459     do_view(oenv,ptotcorr,NULL);
 460   }
 461
 462   return rmsf;
 463 }
 464