src/kernel/mk_ghat.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  * GRowing Old MAkes el Chrono Sweat
  28  */
  29 static char *SRCID_mk_ghat_c = "$Id$";
  30
  31 #include <math.h>
  32 #include <stdio.h>
  33 #include "copyrite.h"
  34 #include "macros.h"
  35 #include "smalloc.h"
  36 #include "typedefs.h"
  37 #include "futil.h"
  38 #include "xvgr.h"
  39 #include "vec.h"
  40 #include "maths.h"
  41 #include "shift_util.h"
  42 #include "physics.h"
  43 #include "statutil.h"
  44 #include "tpxio.h"
  45 #include "fftgrid.h"
  46 #include "copyrite.h"
  47
  48 const real tol = 1e-8;
  49
  50 real sinx_x(real x)
  51 {
  52   if (x == 0)
  53     return 1.0;
  54   else
  55     return sin(x)/x;
  56 }
  57
  58 real  uhat(int porder,real k1,real k2,real k3,real h1,real h2,real h3)
  59 {
  60   real fac1,fac2,fac3;
  61   real f123,fff;
  62   int  i;
  63
  64   fac1 = sinx_x(k1*h1*0.5);
  65   fac2 = sinx_x(k2*h2*0.5);
  66   fac3 = sinx_x(k3*h3*0.5);
  67
  68   fff  = 1;
  69   f123 = fac1*fac2*fac3;
  70   for(i=1; (i<=porder+1); i++)
  71     fff *= f123;
  72
  73   if (fabs(fff) < tol)
  74     return 0.0;
  75   else
  76     return fff;
  77 }
  78
  79 real  uhat1D(int porder,real k1,real h1)
  80 {
  81   real fac1;
  82   real fff;
  83   int  i;
  84
  85   fac1 = sinx_x(k1*h1*0.5);
  86
  87   fff  = 1;
  88   for(i=1; (i<=porder+1); i++)
  89     fff *= fac1;
  90
  91   if (fabs(fff) < tol)
  92     return 0.0;
  93   else
  94     return fff;
  95 }
  96
  97 real shat(real acut,real kmag,real r1)
  98 {
  99   return gk(kmag,acut,r1);
 100 }
 101
 102 real dhat(real alpha,real k,real h)
 103 {
 104   real dh;
 105
 106   dh = alpha*(sin(k*h)/h) + (1.0-alpha)*(sin(2.0*k*h)/(2.0*h) );
 107
 108   if (fabs(dh) < tol)
 109     return 0.0;
 110   else
 111     return dh;
 112 }
 113
 114 real cufrth(int porder,real k1,real k2,real k3,
 115             real h1,real h2,real h3,int nalias)
 116 {
 117   real kn1,kn2,kn3,tmp;
 118   int  n1,n2,n3;
 119   real ufrth;
 120
 121   real twopi=2*M_PI;
 122
 123   if (nalias == 0) {
 124     tmp   = uhat(porder,k1,k2,k3,h1,h2,h3);
 125     ufrth = tmp*tmp*tmp*tmp;
 126   }
 127   else {
 128     ufrth = 0.0;
 129     for( n1 = -nalias; (n1<= nalias); n1++) {
 130       for( n2 = -nalias; (n2<= nalias); n2++) {
 131         for( n3 = -nalias; (n3<= nalias); n3++) {
 132           kn1 = k1 + n1*twopi/h1;
 133           kn2 = k2 + n2*twopi/h2;
 134           kn3 = k3 + n3*twopi/h3;
 135           tmp = uhat(porder,kn1,kn2,kn3,h1,h2,h3);
 136           ufrth += tmp*tmp*tmp*tmp;
 137         }
 138       }
 139     }
 140   }
 141   if (fabs(ufrth) < tol)
 142     return 0.0;
 143   else
 144     return ufrth;
 145 }
 146
 147 real crsqal(real acut,real r1,real k1,real k2,real k3,
 148             real h1,real h2,real h3,int nalias)
 149 {
 150   real kn1,kn2,kn3;
 151   real h_1,h_2,h_3;
 152   real ksq,kmag,tmp,Rsqal;
 153   int  n1,n2,n3;
 154
 155   real twopi=2*M_PI;
 156
 157   h_1=twopi/h1;
 158   h_2=twopi/h2;
 159   h_3=twopi/h3;
 160   if (nalias==0) {
 161     ksq   = k1*k1 + k2*k2 + k3*k3;
 162     kmag  = sqrt(ksq);
 163     tmp   = shat(acut,kmag,r1);
 164     Rsqal = tmp*tmp/(ksq);
 165   }
 166   else {
 167     Rsqal = 0.0;
 168     for(n1 = -nalias; (n1<= nalias); n1++) {
 169       for( n2 = -nalias; (n2<= nalias); n2++) {
 170         for( n3 = -nalias; (n3<= nalias); n3++) {
 171           kn1   = k1 + n1*h_1;
 172           kn2   = k2 + n2*h_2;
 173           kn3   = k3 + n3*h_3;
 174           ksq   = kn1*kn1 + kn2*kn2 + kn3*kn3;
 175           kmag  = sqrt(ksq);
 176           tmp   = shat(acut,kmag,r1);
 177           Rsqal = Rsqal + tmp*tmp/ksq;
 178         }
 179       }
 180     }
 181   }
 182   return Rsqal;
 183 }
 184
 185
 186 real usqsq(int porder,real k1,real k2,real k3,real h1,real h2,real h3)
 187 {
 188   const real tt=2.0/3.0;
 189   const real tx=2.0/15.0;
 190   real t1,t2,t3,t12,t22,t32,tmp;
 191
 192   t1 = sin(k1*h1*0.5);
 193   t2 = sin(k2*h2*0.5);
 194   t3 = sin(k3*h3*0.5);
 195   t12 = t1*t1;
 196   t22 = t2*t2;
 197   t32 = t3*t3;
 198
 199   if (porder == 1)
 200     tmp = (1.0-tt*t12)*(1-tt*t22)*(1-tt*t32);
 201   else if (porder == 2)
 202     tmp = ( (1.0 - t12 + tx*t12*t12)*
 203             (1.0 - t22 + tx*t22*t22)*
 204             (1.0 - t32 + tx*t32*t32) );
 205   else
 206     fatal_error(0,"porder = %d in usqsq",porder);
 207
 208   return tmp*tmp;
 209 }
 210
 211 real usqsq1D(int porder,real k1,real h1)
 212 {
 213   const real tt=2.0/3.0;
 214   const real tx=2.0/15.0;
 215   real t1,t12,tmp;
 216
 217   t1 = sin(k1*h1*0.5);
 218   t12 = t1*t1;
 219
 220   if (porder == 1)
 221     tmp = (1.0-tt*t12);
 222   else if (porder == 2)
 223     tmp = ( (1.0 - t12 + tx*t12*t12));
 224   else
 225     fatal_error(0,"porder = %d in usqsq",porder);
 226
 227   return tmp*tmp;
 228 }
 229
 230 real  ursum(int term,int porder,real acut,real r1,
 231             real k1,real k2,real k3,real h1,real h2,real h3,int nalias)
 232 {
 233   real kt,ksq,kmag;
 234   /*   real kcutsq; */
 235   real kn1,kn2,kn3,urs,tmp;
 236   real h_1,h_2,h_3;
 237   int  n1,n2,n3;
 238
 239   real twopi=2*M_PI;
 240   h_1=twopi/h1;
 241   h_2=twopi/h2;
 242   h_3=twopi/h3;
 243   /*
 244     c
 245     c     for large enough values of k, the terms become negligable
 246     c     if shat(k) = exp(-k^2/4*acut) < eps
 247     c     kcutsq = 4*alpha* (-ln(eps))
 248     c     eps = 10^-6, -ln(eps) = 14
 249     c     eps = 10^-10, -ln(eps) = 23
 250     c     eps = 10^-20, -ln(eps) = 46
 251     c
 252     c     kcutsq = 4.0*acut*115;
 253   */
 254
 255   if (nalias==0) {
 256     if (term==XX) kt = k1;
 257     if (term==YY) kt = k2;
 258     if (term==ZZ) kt = k3;
 259     ksq = k1*k1 + k2*k2 + k3*k3;
 260     kmag = sqrt(ksq);
 261     tmp = uhat(porder,k1,k2,k3,h1,h2,h3);
 262     urs = tmp*tmp*kt*shat(acut,kmag,r1)/(ksq);
 263   }
 264   else {
 265     urs = 0.0;
 266     for(n1 = -nalias; (n1<= nalias); n1++) {
 267       for( n2 = -nalias; (n2<= nalias); n2++) {
 268         for( n3 = -nalias; (n3<= nalias); n3++) {
 269           kn1 = k1 + n1*h_1;
 270           kn2 = k2 + n2*h_2;
 271           kn3 = k3 + n3*h_3;
 272           ksq = kn1*kn1 + kn2*kn2 + kn3*kn3;
 273
 274           if (term==XX) kt = kn1;
 275           if (term==YY) kt = kn2;
 276           if (term==ZZ) kt = kn3;
 277           if (kt != 0.0) {
 278             kmag = sqrt(ksq);
 279             tmp = uhat(porder,kn1,kn2,kn3,h1,h2,h3);
 280             if (tmp != 0.0)
 281               urs = urs + tmp*tmp*kt*shat(acut,kmag,r1)/ksq;
 282           }
 283         }
 284       }
 285     }
 286   }
 287   return urs;
 288 }
 289
 290 real  ursum1D(int term,int porder,real acut,real r1,real k1,real h1,int nalias)
 291 {
 292   real kt,ksq,kmag;
 293 /*   real kcutsq; */
 294   real kn1,urs,tmp;
 295   real h_1;
 296   int  n1;
 297
 298   real twopi=2*M_PI;
 299   h_1=twopi/h1;
 300   /*
 301     c
 302     c     for large enough values of k, the terms become negligable
 303     c     if shat(k) = exp(-k^2/4*acut) < eps
 304     c     kcutsq = 4*alpha* (-ln(eps))
 305     c     eps = 10^-6, -ln(eps) = 14
 306     c     eps = 10^-10, -ln(eps) = 23
 307     c     eps = 10^-20, -ln(eps) = 46
 308     c
 309     */
 310 /*   kcutsq = 4.0*acut*115; */
 311
 312   if (nalias==0) {
 313     if (term==1) kt = k1;
 314     ksq = k1*k1;
 315     kmag = sqrt(ksq);
 316     tmp = uhat1D(porder,k1,h1);
 317     urs = tmp*tmp*kt*shat(acut,kmag,r1)/(EPSILON0*ksq);
 318   }
 319   else {
 320     urs = 0.0;
 321     for(n1 = -nalias; (n1<= nalias); n1++) {
 322       kn1 = k1 + n1*h_1;
 323       ksq = kn1*kn1;
 324       /*c              if (ksq.lt.kcutsq) then*/
 325       if (term==XX) kt = kn1;
 326       if (kt != 0.0) {
 327         kmag = sqrt(ksq);
 328         tmp = uhat1D(porder,kn1,h1);
 329         if (tmp != 0.0)
 330           urs = urs + tmp*tmp*kt*shat(acut,kmag,r1)/(EPSILON0*ksq);
 331       }
 332       /*c              endif*/
 333     }
 334   }
 335   return urs;
 336 }
 337
 338 real sym(int indx,int maxind)
 339 {
 340   if ( (indx == 0 ) || (indx == maxind/2) )
 341     return 1.0;
 342   else
 343     return 2.0;
 344 }
 345
 346 void calc(bool bSym,bool bVerbose,
 347           const int n1max,const int n2max,const int n3max,
 348           const real h1,const real h2,const real h3,
 349           int nalias,int porder,real acut,real r1,const real alpha,
 350           const bool bSearch,
 351           real ***ghat,real *ppval,real *zzval,real *eeref,real *qqopt)
 352 {
 353   real box1,box2,box3;
 354   real k1,k2,k3,ksq,kmag;
 355   real gnumer,dsq,gdenom,gsq;
 356   real ufrth,rsqal,rsq;
 357   real symfac;
 358   int  l1,l2,l3;
 359   real twopi=2*M_PI;
 360   real d1,d2,d3,u1,u2,u3,ss,gg;
 361   real pval,zval,eref,qopt;
 362   int  N1MAX,N2MAX,N3MAX;
 363
 364   if (bSym) {
 365     N1MAX = n1max/2+1;
 366     N2MAX = n2max/2+1;
 367     N3MAX = n3max/2+1;
 368   }
 369   else {
 370     N1MAX = n1max;
 371     N2MAX = n2max;
 372     N3MAX = n3max;
 373   }
 374
 375   box1 = n1max*h1;
 376   box2 = n2max*h2;
 377   box3 = n3max*h3;
 378
 379   pval = 0.0;
 380   zval = 0.0;
 381   eref = 0.0;
 382   qopt = 0.0;
 383
 384   for(l1=0; (l1<N1MAX); l1++) {
 385     if (bVerbose)
 386       fprintf(stderr,"\rl1=%5d  qopt=%12.6e",l1,qopt);
 387
 388     k1   = twopi*l1/box1;
 389     d1   = dhat(alpha,k1,h1);
 390
 391     for(l2=0; (l2<N2MAX); l2++) {
 392       k2   = twopi*l2/box2;
 393       d2   = dhat(alpha,k2,h2);
 394
 395       for(l3=0; (l3<N3MAX); l3++) {
 396         if (((l1+l2+l3) == 0) /*|| (l1 == n1max/2) || (l2 == n2max/2) ||
 397             (l3 == n3max/2)*/)
 398           ghat[0][0][0] = 0.0;
 399         else {
 400           k3   = twopi*l3/box3;
 401           ksq  = k1*k1 + k2*k2 + k3*k3;
 402           kmag = sqrt(ksq);
 403
 404           d3   = dhat(alpha,k3,h3);
 405           u1   = ursum(XX,porder,acut,r1,k1,k2,k3,h1,h2,h3,nalias);
 406           u2   = ursum(YY,porder,acut,r1,k1,k2,k3,h1,h2,h3,nalias);
 407           u3   = ursum(ZZ,porder,acut,r1,k1,k2,k3,h1,h2,h3,nalias);
 408
 409           gnumer = d1*u1+d2*u2+d3*u3;
 410           dsq    = d1*d1+d2*d2+d3*d3;
 411           gdenom = dsq*usqsq(porder,k1,k2,k3,h1,h2,h3);
 412           if (bSym)
 413             symfac = sym(l1,n1max)*sym(l2,n2max)*sym(l3,n3max);
 414           else
 415             symfac = 1.0;
 416
 417           rsqal  = crsqal(acut,r1,k1,k2,k3,h1,h2,h3,nalias);
 418
 419           if (gdenom != 0)
 420             qopt  += (symfac*(rsqal - sqr(gnumer)/gdenom));
 421           if (debug)
 422             fprintf(debug,"rsqal: %10.3e, gnumer: %10.3e, gdenom: %10.3e, ratio: %10.3e\n",
 423                     rsqal,gnumer,gdenom,gnumer/gdenom);
 424
 425 #ifdef DEBUG
 426           if ((l1 == n1max/2) || (l2 == n2max/2) || (l3 == n3max/2))
 427             printf("L(%2d,%2d,%2d)  D(%10.3e,%10.3e,%10.3e) U(%10.3e,%10.3e,%10.3e) gnumer=%10.3em dsq=%10.3e, gdenom=%10.3e, ghat=%10.3e\n",
 428                    l1,l2,l3,d1,d2,d3,u1,u2,u3,gnumer,dsq,gdenom,
 429                    (gdenom == 0.0) ? 0 : gnumer/gdenom);
 430 #endif
 431           if (!bSearch) {
 432             if (gdenom != 0)
 433               gg = gnumer/gdenom;
 434             else
 435               gg = 0.0;
 436             ghat[l1][l2][l3] = gg/EPSILON0;
 437             gsq = gg*gg;
 438
 439             ufrth = cufrth(porder,k1,k2,k3,h1,h2,h3,nalias);
 440             pval  = pval + symfac*
 441               (dsq*gsq*(usqsq(porder,k1,k2,k3,h1,h2,h3)-ufrth));
 442             if (gdenom != 0)
 443               zval  = zval + symfac*
 444                 (dsq*gsq*ufrth - 2.0*gnumer*gnumer/gdenom + rsqal);
 445             ss    = shat(acut,kmag,r1);
 446             rsq   = ss*ss/ksq;
 447             eref  = eref + symfac* (rsqal - rsq);
 448           }
 449         }
 450       }
 451     }
 452   }
 453   if (bVerbose)
 454     fprintf(stderr,"\n");
 455   *ppval = pval/(box1*box2*box3);
 456   *zzval = zval/(box1*box2*box3);
 457   *eeref = eref/(box1*box2*box3);
 458   *qqopt = qopt/(EPSILON0*box1*box2*box3);
 459 }
 460
 461 void calc1D(bool bSym,bool bVerbose,
 462             const int n1max,const int n2max,const int n3max,
 463             const real h1,const real h2,const real h3,
 464             int nalias,int porder,real acut,real r1,const real alpha,
 465             const bool bSearch,
 466             real ***ghat,real *ppval,real *zzval,real *eeref,real *qqopt)
 467 {
 468   real box1,box2,box3;
 469   real k1,k2,k3;
 470   real gnumer,dsq,gdenom;
 471   real rsqal;
 472   real symfac;
 473   int  l1;
 474   real twopi=2*M_PI;
 475   real d1,u1;
 476   real pval,zval,eref,qopt;
 477   int  N1MAX;
 478 /*   int  N2MAX,N3MAX; */
 479
 480   if (bSym) {
 481     N1MAX = n1max/2+1;
 482 /*     N2MAX = n2max/2+1; */
 483 /*     N3MAX = n3max/2+1; */
 484   }
 485   else {
 486     N1MAX = n1max;
 487 /*     N2MAX = n2max; */
 488 /*     N3MAX = n3max; */
 489   }
 490
 491   box1 = n1max*h1;
 492   box2 = n2max*h2;
 493   box3 = n3max*h3;
 494
 495   pval = 0.0;
 496   zval = 0.0;
 497   eref = 0.0;
 498   qopt = 0.0;
 499
 500   k2 = k3 = 0;
 501
 502   for(l1=0; (l1<N1MAX); l1++) {
 503     if (bVerbose)
 504       fprintf(stderr,"\rl1=%5d  qopt=%12.6e",l1,qopt);
 505
 506     k1   = twopi*l1/box1;
 507     d1   = dhat(alpha,k1,h1);
 508
 509     if (l1 == 0)
 510       ghat[0][0][0] = 0.0;
 511     else {
 512       u1   = ursum1D(XX,porder,acut,r1,k1,h1,nalias);
 513
 514       gnumer = d1*u1;
 515       dsq    = d1*d1;
 516       gdenom = dsq*usqsq(porder,k1,k2,k3,h1,h2,h3);
 517       if (bSym)
 518         symfac = sym(l1,n1max);
 519       else
 520         symfac = 1.0;
 521
 522       rsqal  = crsqal(acut,r1,k1,k2,k3,h1,h2,h3,nalias);
 523
 524       if (gdenom != 0)
 525         qopt  += symfac*(rsqal - (gnumer*gnumer)/gdenom);
 526     }
 527   }
 528   if (bVerbose)
 529     fprintf(stderr,"\n");
 530   *ppval = pval/(box1*box2*box3);
 531   *zzval = zval/(box1*box2*box3);
 532   *eeref = eref/(box1*box2*box3);
 533   *qqopt = qopt/(box1*box2*box3);
 534 }
 535
 536 void read_params(char *fn,t_inputrec *ir,rvec boxs)
 537 {
 538   real   t,lambda;
 539   int    step,natoms,m;
 540   matrix box;
 541
 542   /* Read topology and coordinates */
 543   read_tpx(fn,&step,&t,&lambda,ir,box,&natoms,NULL,NULL,NULL,NULL);
 544   for(m=0; (m<DIM); m++)
 545     boxs[m] = box[m][m];
 546 }
 547
 548 int main(int argc,char *argv[])
 549 {
 550   /* FILE   *fp; */
 551   const  real gold=0.38197;
 552   const  int  mxiter=12;
 553   int    n1max,n2max,n3max;
 554   real   h1,h2,h3;
 555   real   box[3];
 556   int    nalias,porder;
 557   real   acut,alpha,r1;
 558   rvec   beta;
 559   bool   bSearch,bConv;
 560   /* bool   bNL=FALSE; */
 561   real   ***ghat;
 562   real   pval,zval,eref,qopt,norm;
 563   real   alpha0,alpha1,alpha2,alpha3,alptol;
 564   real   q0,q1,q2,q3;
 565   int    niter;
 566   /* int    ii,jj,kk,nn; */
 567   t_inputrec ir;
 568   t_filenm fnm[] = {
 569     { efTPX, NULL, NULL,   ffREAD },
 570     { efHAT, "-o", "ghat", ffWRITE }
 571   };
 572 #define NFILE asize(fnm)
 573   static bool bVerbose=FALSE,bCubic=TRUE,bSym=TRUE;
 574   static t_pargs pa[] = {
 575     { "-v",     FALSE, etBOOL, &bVerbose, "Verbose on"},
 576     { "-cubic", FALSE, etBOOL, &bCubic,   "Force beta to be the same in all directions" },
 577     { "-sym",   FALSE, etBOOL, &bSym,     "HIDDENUse symmetry for the generation of ghat function (turn off for debugging only!)" }
 578   };
 579
 580   CopyRight(stdout,argv[0]);
 581   parse_common_args(&argc,argv,0,TRUE,NFILE,fnm,asize(pa),pa,0,NULL,0,NULL);
 582
 583   read_params(ftp2fn(efTPX,NFILE,fnm),&ir,box);
 584   acut   = ir.rcoulomb;
 585   r1     = ir.rcoulomb_switch;
 586   n1max  = ir.nkx;
 587   n2max  = ir.nky;
 588   n3max  = ir.nkz;
 589   h1     = box[XX]/n1max;
 590   h2     = box[YY]/n2max;
 591   h3     = box[ZZ]/n3max;
 592
 593   /* These are not parameters. They are fixed.
 594    * Luty et al. determined their optimal values to be 2.
 595    */
 596   nalias = 2;
 597   porder = 2;
 598   /*bSym   = TRUE;*/
 599
 600   set_LRconsts(stdout,r1,acut,box,NULL);
 601
 602   printf("Grid cell size is %8.4f x %8.4f x %8.4f nm\n",h1,h2,h3);
 603
 604   ghat=mk_rgrid(n1max,n2max,n3max);
 605
 606   bSearch = TRUE;
 607
 608   niter = 0;
 609   if (!bSearch) {
 610     /* c          alpha = 1.33*/
 611     alpha = 1.0;
 612     calc(bSym,bVerbose,n1max,n2max,n3max,h1,h2,h3,
 613          nalias,porder,acut,r1,alpha,bSearch,
 614          ghat,&pval,&zval,&eref,&qopt);
 615   }
 616   else /* Elsje */ {
 617     alpha0 = 1.0;
 618     alpha1 = 4.0/3.0;
 619     alpha3 = 2.0;
 620     calc(bSym,bVerbose,n1max,n2max,n3max,h1,h2,h3,
 621          nalias,porder,acut,r1,alpha0,bSearch,
 622          ghat,&pval,&zval,&eref,&q0);
 623     calc(bSym,bVerbose,n1max,n2max,n3max,h1,h2,h3,
 624          nalias,porder,acut,r1,alpha1,bSearch,
 625          ghat,&pval,&zval,&eref,&q1);
 626     calc(bSym,bVerbose,n1max,n2max,n3max,h1,h2,h3,
 627          nalias,porder,acut,r1,alpha3,bSearch,
 628          ghat,&pval,&zval,&eref,&q3);
 629
 630     /* if ( (q1 > q0) || (q1 > q3) )
 631       fatal_error(0,"oops q1=%f,q0=%f,q3=%f",q1,q0,q3); */
 632     alpha2 = alpha1 + gold*(alpha3-alpha1);
 633     calc(bSym,bVerbose,n1max,n2max,n3max,h1,h2,h3,
 634          nalias,porder,acut,r1,alpha2,bSearch,
 635          ghat,&pval,&zval,&eref,&q2);
 636     bConv  = FALSE;
 637     alptol = 0.05;
 638
 639     fprintf(stderr,"q0 = %10g, q1= %10g, q2 = %10g, q3 = %10g\n",
 640             q0,q1,q2,q3);
 641
 642     while ((!bConv) && (niter < mxiter)) {
 643       fprintf(stderr,"%2d  %10.4f  %10.4f  %10.4f  %10.4f\n",
 644               niter,alpha0,alpha1,alpha2,alpha3);
 645       niter = niter + 1;
 646       if (q2 < q1) {
 647         alpha0 = alpha1;
 648         q0     = q1;
 649         alpha1 = alpha2;
 650         q1     = q2;
 651         alpha2 = alpha1 + gold*(alpha3-alpha1);
 652         calc(bSym,bVerbose,n1max,n2max,n3max,h1,h2,h3,
 653              nalias,porder,acut,r1,alpha2,bSearch,
 654              ghat,&pval,&zval,&eref,&q2);
 655       }
 656       else {
 657         alpha3 = alpha2;
 658         q3     = q2;
 659         alpha2 = alpha1;
 660         q2     = q1;
 661         alpha1 = alpha2 - gold*(alpha2-alpha0);
 662         calc(bSym,bVerbose,n1max,n2max,n3max,h1,h2,h3,
 663              nalias,porder,acut,r1,alpha1,bSearch,
 664              ghat,&pval,&zval,&eref,&q1);
 665       }
 666       if ((alpha3-alpha0) < alptol)
 667         bConv = TRUE;
 668     }
 669     bSearch = FALSE;
 670     if (q1 < q2) {
 671       alpha = alpha1;
 672       calc(bSym,bVerbose,n1max,n2max,n3max,h1,h2,h3,
 673            nalias,porder,acut,r1,alpha,bSearch,
 674            ghat,&pval,&zval,&eref,&qopt);
 675     }
 676     else {
 677       alpha = alpha2;
 678       calc(bSym,bVerbose,n1max,n2max,n3max,h1,h2,h3,
 679            nalias,porder,acut,r1,alpha,bSearch,
 680            ghat,&pval,&zval,&eref,&qopt);
 681     }
 682
 683     fprintf(stderr,"%10g  %10g  %10g  %10g  %10g  %10g\n",
 684             acut,r1,pval,zval,eref,qopt);
 685     norm = sqr(1.0/(4.0*M_PI*h1*h1))*(4.0/3.0*M_PI*h1*h1*h1);
 686     if (norm != 0) {
 687       pval = sqrt(fabs(pval)/norm)*100.0;
 688       zval = sqrt(fabs(zval)/norm)*100.0;
 689       eref = sqrt(fabs(eref)/norm)*100.0;
 690       qopt = sqrt(fabs(qopt)/norm)*100.0;
 691     }
 692     beta[XX] = beta[YY] = beta[ZZ] = alpha;
 693     wr_ghat(ftp2fn(efHAT,NFILE,fnm),n1max,n2max,n3max,h1,h2,h3,
 694             ghat,nalias,porder,niter,bSym,beta,
 695             r1,acut,pval,zval,eref,qopt);
 696
 697     /*    fp=ftp2FILE(efHAT,NFILE,fnm,"w");
 698           fprintf(fp,"%8d  %8d  %8d  %15.10e  %15.10e %15.10e\n",
 699           n1max,n2max,n3max,h1,h2,h3);
 700           fprintf(fp,"%8d  %8d  %8d  %8d  %15.10e  %15.10e  %15.10e\n",
 701           nalias,porder,niter,bSym,alpha,alpha,alpha);
 702           fprintf(fp,"%10g  %10g  %10g  %10g  %10g  %10g\n",
 703           acut,r1,pval,zval,eref,qopt);
 704           {
 705           int  N1MAX,N2MAX,N3MAX;
 706
 707           if (bSym) {
 708           N1MAX = n1max/2+1;
 709           N2MAX = n2max/2+1;
 710           N3MAX = n3max/2+1;
 711           }
 712           else {
 713           N1MAX = n1max;
 714           N2MAX = n2max;
 715           N3MAX = n3max;
 716           }
 717           for(ii=0; (ii<N1MAX); ii++) {
 718           for(jj=0; (jj<N2MAX); jj++) {
 719           for(kk=0,nn=1; (kk<N3MAX); kk++,nn++) {
 720           bNL=FALSE;
 721           fprintf(fp,"  %12.5e",ghat[ii][jj][kk]);
 722           if ((nn % 6) == 0) {
 723           fprintf(fp,"\n");
 724           bNL=TRUE;
 725           }
 726           }
 727           if (!bNL)
 728           fprintf(fp,"\n");
 729           }
 730           }
 731           fclose(fp);
 732           }
 733     */
 734     pr_scalar_gk("ghat.xvg",n1max,n2max,n3max,box,ghat);
 735   }
 736   return 0;
 737 }
 738