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(gmx_bool bSym,gmx_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 gmx_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(gmx_bool bSym,gmx_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 gmx_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   gmx_bool   bSearch,bConv;
 567   /* gmx_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 gmx_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