src/mdlib/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  * Getting the Right Output Means no Artefacts in Calculating Stuff
  31  */
  32
  33 /* This file is completely threadsafe - keep it that way! */
  34
  35 #include <stdio.h>
  36 #include "typedefs.h"
  37 #include "futil.h"
  38 #include "vec.h"
  39 #include "physics.h"
  40 #include "shift_util.h"
  41 #include "fftgrid.h"
  42 #include "xvgr.h"
  43
  44 void symmetrize_ghat(int nx,int ny,int nz,real ***ghat)
  45 {
  46   int  i,j,k;
  47   int  iip,jjp,kkp;
  48   real ggg;
  49
  50   /* fprintf(stderr,"Symmetrizing Ghat function\n");   */
  51   /* Only the lower octant of the rectangle has been saved,
  52    * so we must construct the other 7 octants by symmetry operations.
  53    */
  54   for(i=0; (i<=nx/2); i++) {
  55     iip = (nx-i) % nx;
  56     for(j=0; (j<=ny/2); j++) {
  57       jjp = (ny-j) % ny;
  58       for(k=0; (k<=nz/2); k++) {
  59         kkp = (nz-k) % nz;
  60         ggg                 = ghat[i][j][k];
  61         ghat[i]  [jjp][k]   = ggg;
  62         ghat[i]  [j]  [kkp] = ggg;
  63         ghat[i]  [jjp][kkp] = ggg;
  64         ghat[iip][j]  [k]   = ggg;
  65         ghat[iip][jjp][k]   = ggg;
  66         ghat[iip][j]  [kkp] = ggg;
  67         ghat[iip][jjp][kkp] = ggg;
  68       }
  69     }
  70   }
  71 }
  72
  73 void mk_ghat(FILE *fp,int nx,int ny,int nz,real ***ghat,
  74              rvec box,real r1,real rc,bool bSym,bool bOld)
  75 {
  76   int  ix,iy,iz;
  77   int  ixmax,iymax,izmax;
  78   real k2,ggg;
  79   rvec k,lll;
  80
  81   calc_lll(box,lll);
  82
  83   if (bSym) {
  84     ixmax=nx/2+1;
  85     iymax=ny/2+1;
  86     izmax=nz/2+1;
  87   }
  88   else {
  89     ixmax=nx;
  90     iymax=ny;
  91     izmax=nz;
  92   }
  93   /* Loop over lattice vectors in fourier space */
  94   for(ix=0; (ix < ixmax); ix++) {
  95     for(iy=0; (iy < iymax); iy++) {
  96       for(iz=0; (iz < izmax); iz++) {
  97         calc_k(lll,ix,iy,iz,nx,ny,nz,k);
  98         k2 = iprod(k,k);
  99         if ((ix == 0) && (iy == 0) && (iz == 0))
 100           ggg = 0.0;
 101         else {
 102           if (bOld)
 103             ggg = gk(sqrt(k2),rc,r1)/(k2*EPSILON0);
 104           else
 105             ggg = gknew(sqrt(k2),rc,r1)/(k2*EPSILON0);
 106         }
 107         ghat[ix][iy][iz]=ggg;
 108       }
 109     }
 110   }
 111   if (bSym)
 112     symmetrize_ghat(nx,ny,nz,ghat);
 113 }
 114
 115 void wr_ghat(char *fn,int n1max,int n2max,int n3max,real h1,real h2,real h3,
 116              real ***ghat,int nalias,int porder,int niter,bool bSym,rvec beta,
 117              real r1,real rc,real pval,real zval,real eref,real qopt)
 118 {
 119   FILE *fp;
 120   int  N1MAX,N2MAX,N3MAX;
 121   bool bNL=FALSE;
 122   real rx,ry,rz;
 123   int  ii,jj,kk,nn;
 124
 125   fp = ffopen(fn,"w");
 126   fprintf(fp,"%8d  %8d  %8d  %15.10e  %15.10e %15.10e\n",
 127           n1max,n2max,n3max,h1,h2,h3);
 128   fprintf(fp,"%8d  %8d  %8d  %8d  %15.10e  %15.10e  %15.10e\n",
 129           nalias,porder,niter,bSym,beta[XX],beta[YY],beta[ZZ]);
 130   fprintf(fp,"%10g  %10g  %10g  %10g  %10g  %10g\n",
 131           rc,r1,pval,zval,eref,qopt);
 132
 133   if (bSym) {
 134     N1MAX = n1max/2+1;
 135     N2MAX = n2max/2+1;
 136     N3MAX = n3max/2+1;
 137   }
 138   else {
 139     N1MAX = n1max;
 140     N2MAX = n2max;
 141     N3MAX = n3max;
 142   }
 143   for(ii=0; (ii<N1MAX); ii++) {
 144     for(jj=0; (jj<N2MAX); jj++) {
 145       for(kk=0,nn=1; (kk<N3MAX); kk++,nn++) {
 146         bNL=FALSE;
 147         fprintf(fp,"  %12.5e",ghat[ii][jj][kk]);
 148         if ((nn % 6) == 0) {
 149           fprintf(fp,"\n");
 150           bNL=TRUE;
 151         }
 152       }
 153       if (!bNL)
 154         fprintf(fp,"\n");
 155     }
 156   }
 157   fclose(fp);
 158
 159   fp=xvgropen("ghat.xvg","G-Hat","k","gk");
 160   for(ii=0; (ii<N1MAX); ii++) {
 161     rx=sqr((real)(ii*h1));
 162     for(jj=0; (jj<N2MAX); jj++) {
 163       ry=rx+sqr((real)(jj*h2));
 164       for(kk=0; (kk<N3MAX); kk++) {
 165         rz=ry+sqr((real)(kk*h3));
 166         fprintf(fp,"%10g  %10g\n",sqrt(rz),ghat[ii][jj][kk]);
 167       }
 168     }
 169   }
 170   fclose(fp);
 171 }
 172
 173 void pr_scalar_gk(char *fn,int nx,int ny,int nz,rvec box,real ***ghat)
 174 {
 175   FILE *fp;
 176   int  ii,jj,kk;
 177   real k1;
 178   rvec k,lll;
 179
 180   calc_lll(box,lll);
 181
 182   fp=xvgropen(fn,"G-Hat","k","gk");
 183   for(ii=0; (ii<nx); ii++) {
 184     for(jj=0; (jj<ny); jj++) {
 185       for(kk=0; (kk<nz); kk++) {
 186         calc_k(lll,ii,jj,kk,nx,ny,nz,k);
 187         k1 = norm(k);
 188         fprintf(fp,"%10g  %10g\n",k1,ghat[ii][jj][kk]);
 189       }
 190     }
 191   }
 192   fclose(fp);
 193 }
 194
 195 real ***rd_ghat(FILE *log,char *fn,ivec igrid,rvec gridspace,
 196                 rvec beta,int *porder,real *r1,real *rc)
 197 {
 198   FILE   *in;
 199   real   ***gh;
 200   double gx,gy,gz,alX,alY,alZ,ddd;
 201   double acut,pval,zval,eref,qopt,r11;
 202   int    nalias,niter,bSym;
 203   int    ix,iy,iz,ixmax,iymax,izmax;
 204
 205   in=ffopen(fn,"r");
 206   fscanf(in,"%d%d%d%lf%lf%lf",&ix,&iy,&iz,&gx,&gy,&gz);
 207   igrid[XX]=ix, igrid[YY]=iy, igrid[ZZ]=iz;
 208   gridspace[XX]=gx,  gridspace[YY]=gy,  gridspace[ZZ]=gz;
 209   fscanf(in,"%d%d%d%d%lf%lf%lf",&nalias,porder,&niter,&bSym,&alX,&alY,&alZ);
 210   fscanf(in,"%lf%lf%lf%lf%lf%lf",&acut,&r11,&pval,&zval,&eref,&qopt);
 211   *r1 = r11;
 212   *rc = acut;
 213
 214   fprintf(log,"\nOpened %s for reading ghat function\n",fn);
 215   fprintf(log,"gridsize: %10d %10d %10d\n",ix,iy,iz);
 216   fprintf(log,"spacing:  %10g %10g %10g\n",gx,gy,gz);
 217   fprintf(log,"    nalias    porder     niter      bSym      beta[X-Z]\n"
 218           "%10d%10d%10d%10d%10g%10g%10g\n",
 219           nalias,*porder,niter,bSym,alX,alY,alZ);
 220   fprintf(log,"      acut        r1      pval      zval      eref      qopt\n"
 221           "%10g%10g%10g%10g%10g%10g\n",acut,*r1,pval,zval,eref,qopt);
 222   fflush(log);
 223
 224   beta[XX] = alX;
 225   beta[YY] = alY;
 226   beta[ZZ] = alZ;
 227
 228   gh      = mk_rgrid(ix,iy,iz);
 229   if (bSym) {
 230     ixmax=igrid[XX]/2+1;
 231     iymax=igrid[YY]/2+1;
 232     izmax=igrid[ZZ]/2+1;
 233   }
 234   else {
 235     ixmax=igrid[XX];
 236     iymax=igrid[YY];
 237     izmax=igrid[ZZ];
 238   }
 239   fprintf(log,"Reading ghat of %d %d %d\n",ixmax,iymax,izmax);
 240   for(ix=0; (ix<ixmax); ix++)
 241     for(iy=0; (iy<iymax); iy++)
 242       for(iz=0; (iz<izmax); iz++) {
 243         fscanf(in,"%lf",&ddd);
 244         gh[ix][iy][iz] = ddd;
 245       }
 246   ffclose(in);
 247
 248   wr_ghat("output.hat",igrid[XX],igrid[YY],igrid[ZZ],gx,gy,gz,gh,
 249           nalias,*porder,niter,bSym,beta,
 250           *r1,*rc,pval,zval,eref,qopt);
 251
 252   if (bSym)
 253     symmetrize_ghat(igrid[XX],igrid[YY],igrid[ZZ],gh);
 254
 255   fprintf(log,"\nSuccessfully read ghat function!\n");
 256
 257
 258   return gh;
 259 }
 260