src/mdlib/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  * GROwing Monsters And Cloning Shrimps
  34  */
  35 /* This file is completely threadsafe - keep it that way! */
  36 #ifdef HAVE_CONFIG_H
  37 #include <config.h>
  38 #endif
  39
  40 #include <stdio.h>
  41 #include "typedefs.h"
  42 #include "futil.h"
  43 #include "vec.h"
  44 #include "physics.h"
  45 #include "coulomb.h"
  46 #include "pppm.h"
  47 #include "xvgr.h"
  48 #include "gmxfio.h"
  49 #include "pppm.h"
  50 #include "smalloc.h"
  51
  52 static void calc_k(rvec lll,int ix,int iy,int iz,int nx,int ny,int nz,rvec k)
  53 {
  54 #define IDX(i,n,x)  (i<=n/2) ? (i*x) : ((i-n)*x)
  55   k[XX] = IDX(ix,nx,lll[XX]);
  56   k[YY] = IDX(iy,ny,lll[YY]);
  57   k[ZZ] = IDX(iz,nz,lll[ZZ]);
  58 #undef IDX
  59 }
  60
  61 void symmetrize_ghat(int nx,int ny,int nz,real ***ghat)
  62 /* Symmetrize the Ghat function. It is assumed that the
  63  * first octant of the Ghat function is either read or generated
  64  * (all k-vectors from 0..nx/2 0..ny/2 0..nz/2).
  65  * Since Gk depends on the absolute value of k only,
  66  * symmetry operations may shorten the time to generate it.
  67  */
  68
  69 {
  70   int  i,j,k;
  71   int  iip,jjp,kkp;
  72   real ggg;
  73
  74   /* fprintf(stderr,"Symmetrizing Ghat function\n");   */
  75   /* Only the lower octant of the rectangle has been saved,
  76    * so we must construct the other 7 octants by symmetry operations.
  77    */
  78   for(i=0; (i<=nx/2); i++) {
  79     iip = (nx-i) % nx;
  80     for(j=0; (j<=ny/2); j++) {
  81       jjp = (ny-j) % ny;
  82       for(k=0; (k<=nz/2); k++) {
  83         kkp = (nz-k) % nz;
  84         ggg                 = ghat[i][j][k];
  85         ghat[i]  [jjp][k]   = ggg;
  86         ghat[i]  [j]  [kkp] = ggg;
  87         ghat[i]  [jjp][kkp] = ggg;
  88         ghat[iip][j]  [k]   = ggg;
  89         ghat[iip][jjp][k]   = ggg;
  90         ghat[iip][j]  [kkp] = ggg;
  91         ghat[iip][jjp][kkp] = ggg;
  92       }
  93     }
  94   }
  95 }
  96
  97 void mk_ghat(FILE *fp,int nx,int ny,int nz,real ***ghat,
  98              rvec box,real r1,real rc,bool bSym,bool bOld)
  99 {
 100   int  ix,iy,iz;
 101   int  ixmax,iymax,izmax;
 102   real k2,ggg;
 103   rvec k,lll;
 104
 105   calc_lll(box,lll);
 106
 107   if (bSym) {
 108     ixmax=nx/2+1;
 109     iymax=ny/2+1;
 110     izmax=nz/2+1;
 111   }
 112   else {
 113     ixmax=nx;
 114     iymax=ny;
 115     izmax=nz;
 116   }
 117   /* Loop over lattice vectors in fourier space */
 118   for(ix=0; (ix < ixmax); ix++) {
 119     for(iy=0; (iy < iymax); iy++) {
 120       for(iz=0; (iz < izmax); iz++) {
 121         calc_k(lll,ix,iy,iz,nx,ny,nz,k);
 122         k2 = iprod(k,k);
 123         if ((ix == 0) && (iy == 0) && (iz == 0))
 124           ggg = 0.0;
 125         else {
 126           if (bOld)
 127             ggg = gk(sqrt(k2),rc,r1)/(k2*EPSILON0);
 128           else
 129             ggg = gknew(sqrt(k2),rc,r1)/(k2*EPSILON0);
 130         }
 131         ghat[ix][iy][iz]=ggg;
 132       }
 133     }
 134   }
 135   if (bSym)
 136     symmetrize_ghat(nx,ny,nz,ghat);
 137 }
 138
 139 void wr_ghat(const char *fn,const output_env_t oenv,
 140              int n1max,int n2max,int n3max,real h1,real h2,real h3,
 141              real ***ghat,int nalias,int porder,int niter,bool bSym,rvec beta,
 142              real r1,real rc,real pval,real zval,real eref,real qopt)
 143 {
 144   FILE *fp;
 145   int  N1MAX,N2MAX,N3MAX;
 146   bool bNL=FALSE;
 147   real rx,ry,rz;
 148   int  ii,jj,kk,nn;
 149
 150   fp = gmx_fio_fopen(fn,"w");
 151   fprintf(fp,"%8d  %8d  %8d  %15.10e  %15.10e %15.10e\n",
 152           n1max,n2max,n3max,h1,h2,h3);
 153   fprintf(fp,"%8d  %8d  %8d  %8d  %15.10e  %15.10e  %15.10e\n",
 154           nalias,porder,niter,bSym,beta[XX],beta[YY],beta[ZZ]);
 155   fprintf(fp,"%10g  %10g  %10g  %10g  %10g  %10g\n",
 156           rc,r1,pval,zval,eref,qopt);
 157
 158   if (bSym) {
 159     N1MAX = n1max/2+1;
 160     N2MAX = n2max/2+1;
 161     N3MAX = n3max/2+1;
 162   }
 163   else {
 164     N1MAX = n1max;
 165     N2MAX = n2max;
 166     N3MAX = n3max;
 167   }
 168   for(ii=0; (ii<N1MAX); ii++) {
 169     for(jj=0; (jj<N2MAX); jj++) {
 170       for(kk=0,nn=1; (kk<N3MAX); kk++,nn++) {
 171         bNL=FALSE;
 172         fprintf(fp,"  %12.5e",ghat[ii][jj][kk]);
 173         if ((nn % 6) == 0) {
 174           fprintf(fp,"\n");
 175           bNL=TRUE;
 176         }
 177       }
 178       if (!bNL)
 179         fprintf(fp,"\n");
 180     }
 181   }
 182   gmx_fio_fclose(fp);
 183
 184   fp=xvgropen("ghat.xvg","G-Hat","k","gk",oenv);
 185   for(ii=0; (ii<N1MAX); ii++) {
 186     rx=sqr((real)(ii*h1));
 187     for(jj=0; (jj<N2MAX); jj++) {
 188       ry=rx+sqr((real)(jj*h2));
 189       for(kk=0; (kk<N3MAX); kk++) {
 190         rz=ry+sqr((real)(kk*h3));
 191         fprintf(fp,"%10g  %10g\n",sqrt(rz),ghat[ii][jj][kk]);
 192       }
 193     }
 194   }
 195   gmx_fio_fclose(fp);
 196 }
 197
 198 void pr_scalar_gk(const char *fn,const output_env_t oenv,int nx,int ny,int nz,
 199                   rvec box,real ***ghat)
 200 {
 201   FILE *fp;
 202   int  ii,jj,kk;
 203   real k1;
 204   rvec k,lll;
 205
 206   calc_lll(box,lll);
 207
 208   fp=xvgropen(fn,"G-Hat","k","gk",oenv);
 209   for(ii=0; (ii<nx); ii++) {
 210     for(jj=0; (jj<ny); jj++) {
 211       for(kk=0; (kk<nz); kk++) {
 212         calc_k(lll,ii,jj,kk,nx,ny,nz,k);
 213         k1 = norm(k);
 214         fprintf(fp,"%10g  %10g\n",k1,ghat[ii][jj][kk]);
 215       }
 216     }
 217   }
 218   gmx_fio_fclose(fp);
 219 }
 220
 221 static real ***mk_rgrid(int nx,int ny,int nz)
 222 {
 223   real *ptr1;
 224   real **ptr2;
 225   real ***ptr3;
 226   int  i,j,n2,n3;
 227
 228   snew(ptr1,nx*ny*nz);
 229   snew(ptr2,nx*ny);
 230   snew(ptr3,nx);
 231
 232   n2=n3=0;
 233   for(i=0; (i<nx); i++) {
 234     ptr3[i]=&(ptr2[n2]);
 235     for(j=0; (j<ny); j++,n2++) {
 236       ptr2[n2] = &(ptr1[n3]);
 237       n3 += nz;
 238     }
 239   }
 240   return ptr3;
 241 }
 242
 243 real ***rd_ghat(FILE *log,const output_env_t oenv,char *fn,ivec igrid,
 244                 rvec gridspace, rvec beta,int *porder,real *r1,real *rc)
 245 {
 246   FILE   *in;
 247   real   ***gh;
 248   double gx,gy,gz,alX,alY,alZ,ddd;
 249   double acut,pval,zval,eref,qopt,r11;
 250   int    nalias,niter,bSym;
 251   int    ix,iy,iz,ixmax,iymax,izmax;
 252
 253   in=gmx_fio_fopen(fn,"r");
 254   if(6 != fscanf(in,"%d%d%d%lf%lf%lf",&ix,&iy,&iz,&gx,&gy,&gz))
 255   {
 256       gmx_fatal(FARGS,"Error reading from file %s",fn);
 257   }
 258
 259
 260   igrid[XX]=ix, igrid[YY]=iy, igrid[ZZ]=iz;
 261   gridspace[XX]=gx,  gridspace[YY]=gy,  gridspace[ZZ]=gz;
 262   if(7 != fscanf(in,"%d%d%d%d%lf%lf%lf",&nalias,porder,&niter,&bSym,&alX,&alY,&alZ))
 263   {
 264       gmx_fatal(FARGS,"Error reading from file %s",fn);
 265   }
 266
 267   if(6 != fscanf(in,"%lf%lf%lf%lf%lf%lf",&acut,&r11,&pval,&zval,&eref,&qopt))
 268   {
 269     gmx_fatal(FARGS,"Error reading from file %s",fn);
 270   }
 271
 272   *r1 = r11;
 273   *rc = acut;
 274
 275   fprintf(log,"\nOpened %s for reading ghat function\n",fn);
 276   fprintf(log,"gridsize: %10d %10d %10d\n",ix,iy,iz);
 277   fprintf(log,"spacing:  %10g %10g %10g\n",gx,gy,gz);
 278   fprintf(log,"    nalias    porder     niter      bSym      beta[X-Z]\n"
 279           "%10d%10d%10d%10d%10g%10g%10g\n",
 280           nalias,*porder,niter,bSym,alX,alY,alZ);
 281   fprintf(log,"      acut        r1      pval      zval      eref      qopt\n"
 282           "%10g%10g%10g%10g%10g%10g\n",acut,*r1,pval,zval,eref,qopt);
 283   fflush(log);
 284
 285   beta[XX] = alX;
 286   beta[YY] = alY;
 287   beta[ZZ] = alZ;
 288
 289   gh      = mk_rgrid(ix,iy,iz);
 290   if (bSym) {
 291     ixmax=igrid[XX]/2+1;
 292     iymax=igrid[YY]/2+1;
 293     izmax=igrid[ZZ]/2+1;
 294   }
 295   else {
 296     ixmax=igrid[XX];
 297     iymax=igrid[YY];
 298     izmax=igrid[ZZ];
 299   }
 300   fprintf(log,"Reading ghat of %d %d %d\n",ixmax,iymax,izmax);
 301   for(ix=0; (ix<ixmax); ix++)
 302     for(iy=0; (iy<iymax); iy++)
 303       for(iz=0; (iz<izmax); iz++) {
 304         if( 1 != fscanf(in,"%lf",&ddd))
 305         {
 306             gmx_fatal(FARGS,"Error reading from file %s",fn);
 307         }
 308
 309         gh[ix][iy][iz] = ddd;
 310       }
 311   gmx_fio_fclose(in);
 312
 313   wr_ghat("output.hat",oenv,igrid[XX],igrid[YY],igrid[ZZ],gx,gy,gz,gh,
 314           nalias,*porder,niter,bSym,beta,
 315           *r1,*rc,pval,zval,eref,qopt);
 316
 317   if (bSym)
 318     symmetrize_ghat(igrid[XX],igrid[YY],igrid[ZZ],gh);
 319
 320   fprintf(log,"\nSuccessfully read ghat function!\n");
 321
 322
 323   return gh;
 324 }
 325