/**************************************/ /* Spectral method to calculate the */ /* scalar thermal conductivity and */ /* shear viscosity. */ /* */ /* Version 1.1 Jan.22, 1999 */ /* Ju Li (MIT) */ /**************************************/ #include "smile.h" static void realft(double data[], unsigned long n, int isign); static void four1(double data[], unsigned long nn, int isign); void realft (double data[], unsigned long n, int isign) { void four1(double data[], unsigned long nn, int isign); unsigned long i,i1,i2,i3,i4,np3; double c1=0.5,c2,h1r,h1i,h2r,h2i; double wr,wi,wpr,wpi,wtemp,theta; theta=3.141592653589793/(n>>1); if (isign == 1) { c2 = -0.5; four1(data,n>>1,1); } else { c2=0.5; theta = -theta; } wtemp=sin(0.5*theta); wpr = -2.0*wtemp*wtemp; wpi=sin(theta); wr=1.0+wpr; wi=wpi; np3=n+3; for (i=2;i<=(n>>2);i++) { i4=1+(i3=np3-(i2=1+(i1=i+i-1))); h1r=c1*(data[i1]+data[i3]); h1i=c1*(data[i2]-data[i4]); h2r = -c2*(data[i2]+data[i4]); h2i=c2*(data[i1]-data[i3]); data[i1]=h1r+wr*h2r-wi*h2i; data[i2]=h1i+wr*h2i+wi*h2r; data[i3]=h1r-wr*h2r+wi*h2i; data[i4] = -h1i+wr*h2i+wi*h2r; wr=(wtemp=wr)*wpr-wi*wpi+wr; wi=wi*wpr+wtemp*wpi+wi; } if (isign == 1) { data[1] = (h1r=data[1])+data[2]; data[2] = h1r-data[2]; } else { data[1]=c1*((h1r=data[1])+data[2]); data[2]=c1*(h1r-data[2]); four1(data,n>>1,-1); } /* Added by Ju Li: factor of 2/n is */ /* multiplied for inverse realft(); */ if (isign == -1) for (i=1; i<=n; i++) data[i] *= 2./n; return; } /* end realft() */ #define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr void four1(double data[], unsigned long nn, int isign) { unsigned long n,mmax,m,j,istep,i; double wtemp,wr,wpr,wpi,wi,theta; double tempr,tempi; n=nn << 1; j=1; for (i=1;i<n;i+=2) { if (j > i) { SWAP(data[j],data[i]); SWAP(data[j+1],data[i+1]); } m=n >> 1; while (m >= 2 && j > m) { j -= m; m >>= 1; } j += m; } mmax=2; while (n > mmax) { istep=mmax << 1; theta=isign*(6.28318530717959/mmax); wtemp=sin(0.5*theta); wpr = -2.0*wtemp*wtemp; wpi=sin(theta); wr=1.0; wi=0.0; for (m=1;m<mmax;m+=2) { for (i=m;i<=n;i+=istep) { j=i+mmax; tempr=wr*data[j]-wi*data[j+1]; tempi=wr*data[j+1]+wi*data[j]; data[j]=data[i]-tempr; data[j+1]=data[i+1]-tempi; data[i] += tempr; data[i+1] += tempi; } wr=(wtemp=wr)*wpr-wi*wpi+wr; wi=wi*wpr+wtemp*wpi+wi; } mmax=istep; } } /* end four1() */ #undef SWAP /*****************************************************/ /* Calculate the power spectrum and auto-correlation */ /* function of a "current" averaged over three */ /* orthogonal components and write result to OUTPUT. */ /*****************************************************/ double smile (int w, FILE *output) { int i, j, k, l; unsigned long size, n, ncorr, nfreq; double delta, total, first_dip, *data, *corr, freq[SMILE_FREQ_CHANNELS] = {0}; FILE *current[3], *freq_output, *corr_output; if ( (w!=SMILE_HEAT) && (w!=SMILE_STRESS) ) { fprintf (stderr, "smile: illegal option w = %d.\n", w); exit(1); } for (i=0; i<3; i++) if ((current[i]=fopen(fn_flux[w][i],"r")) == NULL) { fprintf (stderr, "smile: cannot open file \"%s\".\n", fn_flux[w][i]); exit(1); } /* the first double precision number is delta in ps */ fread (&delta, sizeof(double), 1, current[0]); printf ("Found delta = %.4f ps\n", delta); fseek (current[0], 0L, SEEK_END); size = ftell(current[0]); if (size/sizeof(double)*sizeof(double)!=size) { fprintf (stderr, "smile: size of \"%s\" = %ld.\n", fn_flux[w][0], size); exit(1); } size = size/sizeof(double)-1; /* FFT only uses the chunk that is integer power of 2 */ for (n=2; n<=size; n*=2); n/=2; printf ("\nUsable %s current length = %ld,\n", fn_flux[w][0], n); if ( (data=(double*)malloc(n*sizeof(double))) == NULL ) { fprintf (stderr, "smile: not enough memory to load in %s\n.", fn_flux[w][0]); exit(1); } /* conform to Fortran array convention */ data--; nfreq = (SMILE_FREQ_CHANNELS>n/2)?n/2:SMILE_FREQ_CHANNELS; for (j=0; j<3; j++) { printf ("Loading in \"%s\"...\n", fn_flux[w][j]); fseek (current[j], -n*sizeof(double), SEEK_END); fread (&data[1], sizeof(double), n, current[j]); fclose (current[j]); printf ("raw data loaded, FFT this current...\n"); if (n>=2) realft (data, n, 1); printf ("FFT complete, getting the power spectrum...\n"); /* get rid of net drift */ data[1] = 0.; if (n>=2) data[2] = data[2]*data[2]/n; for (i=3; i<=n; i+=2) { data[i] = (data[i]*data[i]+data[i+1]*data[i+1])/n; data[i+1] = 0.; } /* divide by 2 to be the single sided integral */ for (i=0; i<nfreq; i++) freq[i] += data[2*i+1]/3./2.; printf ("Inverse FFT the power spectrum...\n"); if (n>=2) realft (data, n, -1); printf ("complete, accumulate to counters.\n\n"); if (j==0) { /* search approx. location corr. turns negative */ for (ncorr=0; (ncorr<n)&&(data[1+ncorr/SMILE_CORR_EXTENSION]>0.); ncorr++); corr = (double *) malloc((ncorr+1)*sizeof(double)); corr--; for (i=1; i<=ncorr; i++) corr[i] = data[i]/3.; } else for (i=1; i<=ncorr; i++) corr[i] += data[i]/3.; } /* loop over principle directions */ /* look for the first point where the */ /* correlation function turns negative */ for (i=1; (i<=ncorr)&&(corr[i]>0.); i++); for (total=0.,l=1; l<i; l++) total += corr[l]; fprintf (output, "\n Based on a piece of %.1f ps current data (x 3),\n", n * delta); fprintf (output, (w==SMILE_HEAT)? " <J(0)^2>/kT^2V = %f W/M/K/ps,\n": " <s_xy(0)^2>V/kT = %e Pa,\n", corr[1]/delta); fprintf (output, " the first dip in correlation function occurs at %.3f ps\n", i*delta); if (w == SMILE_HEAT) fprintf (output, " and by then the thermal conductivity is %.3f mW/m/K.\n\n", 1000*total); else fprintf (output, " and by then the shear viscosity is %.3f uPa.second.\n\n", 1000000*total); /* save the correlation function */ corr_output = fopen (fn_corr[w], "w+"); for (i=1; i<=ncorr; i++) fprintf (corr_output, "%.4f %e\n", (i-1)*delta, corr[i]); fclose (corr_output); /* save the power spectrum */ freq_output = fopen (fn_freq[w], "w+"); for (i=0; i<nfreq; i++) fprintf (freq_output, "%d %f\n", i, freq[i]); fclose (freq_output); /* free the allocations */ free(data+1); free(corr+1); return (total); } /* end smile() */ #ifdef SMILE_TEST int main (int argc, char *argv[]) { if (argc!=2) { printf ("Usage: smile 0 (conductivity), or 1 (viscosity)\n"); return(1); } smile (atoi(argv[1]), stdout); return(0); } /* end main() */ #endif