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