3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
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.
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.
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.
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.
30 * For more info, check our website at http://www.gromacs.org
33 * GROwing Monsters And Cloning Shrimps
35 /* This file is completely threadsafe - keep it that way! */
52 static void calc_k(rvec lll,int ix,int iy,int iz,int nx,int ny,int nz,rvec k)
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]);
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.
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.
78 for(i=0; (i<=nx/2); i++) {
80 for(j=0; (j<=ny/2); j++) {
82 for(k=0; (k<=nz/2); 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;
97 void mk_ghat(FILE *fp,int nx,int ny,int nz,real ***ghat,
98 rvec box,real r1,real rc,gmx_bool bSym,gmx_bool bOld)
101 int ixmax,iymax,izmax;
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);
123 if ((ix == 0) && (iy == 0) && (iz == 0))
127 ggg = gk(sqrt(k2),rc,r1)/(k2*EPSILON0);
129 ggg = gknew(sqrt(k2),rc,r1)/(k2*EPSILON0);
131 ghat[ix][iy][iz]=ggg;
136 symmetrize_ghat(nx,ny,nz,ghat);
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,gmx_bool bSym,rvec beta,
142 real r1,real rc,real pval,real zval,real eref,real qopt)
145 int N1MAX,N2MAX,N3MAX;
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);
168 for(ii=0; (ii<N1MAX); ii++) {
169 for(jj=0; (jj<N2MAX); jj++) {
170 for(kk=0,nn=1; (kk<N3MAX); kk++,nn++) {
172 fprintf(fp," %12.5e",ghat[ii][jj][kk]);
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]);
198 void pr_scalar_gk(const char *fn,const output_env_t oenv,int nx,int ny,int nz,
199 rvec box,real ***ghat)
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);
214 fprintf(fp,"%10g %10g\n",k1,ghat[ii][jj][kk]);
221 static real ***mk_rgrid(int nx,int ny,int nz)
233 for(i=0; (i<nx); i++) {
235 for(j=0; (j<ny); j++,n2++) {
236 ptr2[n2] = &(ptr1[n3]);
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)
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;
253 in=gmx_fio_fopen(fn,"r");
254 if(6 != fscanf(in,"%d%d%d%lf%lf%lf",&ix,&iy,&iz,&gx,&gy,&gz))
256 gmx_fatal(FARGS,"Error reading from file %s",fn);
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))
264 gmx_fatal(FARGS,"Error reading from file %s",fn);
267 if(6 != fscanf(in,"%lf%lf%lf%lf%lf%lf",&acut,&r11,&pval,&zval,&eref,&qopt))
269 gmx_fatal(FARGS,"Error reading from file %s",fn);
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);
289 gh = mk_rgrid(ix,iy,iz);
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))
306 gmx_fatal(FARGS,"Error reading from file %s",fn);
309 gh[ix][iy][iz] = ddd;
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);
318 symmetrize_ghat(igrid[XX],igrid[YY],igrid[ZZ],gh);
320 fprintf(log,"\nSuccessfully read ghat function!\n");