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