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