2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2010,2011,2013,2014, by the GROMACS development team, led by
5 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
6 * and including many others, as listed in the AUTHORS file in the
7 * top-level source directory and at http://www.gromacs.org.
9 * GROMACS is free software; you can redistribute it and/or
10 * modify it under the terms of the GNU Lesser General Public License
11 * as published by the Free Software Foundation; either version 2.1
12 * of the License, or (at your option) any later version.
14 * GROMACS is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 * Lesser General Public License for more details.
19 * You should have received a copy of the GNU Lesser General Public
20 * License along with GROMACS; if not, see
21 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
22 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
24 * If you want to redistribute modifications to GROMACS, please
25 * consider that scientific software is very special. Version
26 * control is crucial - bugs must be traceable. We will be happy to
27 * consider code for inclusion in the official distribution, but
28 * derived work must not be called official GROMACS. Details are found
29 * in the README & COPYING files - if they are missing, get the
30 * official version at http://www.gromacs.org.
32 * To help us fund GROMACS development, we humbly ask that you cite
33 * the research papers on the package. Check out http://www.gromacs.org.
41 #include "gromacs/fft/fft.h"
42 #include "gmx_fatal.h"
43 #include "gromacs/fileio/futil.h"
45 #include "powerspect.h"
47 void addtoavgenergy(t_complex *list, real *result, int size, int tsteps)
50 for (i = 0; i < size; i++)
52 result[i] += cabs2(list[i])/tsteps;
58 void powerspectavg(real ***intftab, int tsteps, int xbins, int ybins, char **outfiles)
60 /*Fourier plans and output;*/
62 t_complex *ftspect1; /* Spatial FFT of interface for each time frame and interface ftint[time,xycoord][0], ftintf[time,xycoord][1] for interface 1 and 2 respectively */
64 real *pspectavg1; /*power -spectrum 1st interface*/
65 real *pspectavg2; /* -------------- 2nd interface*/
67 FILE *datfile1, *datfile2; /*data-files with interface data*/
69 int fy = ybins/2+1; /* number of (symmetric) fourier y elements; */
70 int rfl = xbins*fy; /*length of real - DFT == Symmetric 2D matrix*/
73 /*Prepare data structures for FFT, with time averaging of power spectrum*/
74 if ( (status = gmx_fft_init_2d_real(&fftp, xbins, ybins, GMX_FFT_FLAG_NONE) ) != 0)
77 gmx_fatal(status, __FILE__, __LINE__, "Error allocating FFT");
84 snew(pspectavg1, rfl);
85 snew(pspectavg2, rfl);
87 /*Fouriertransform directly (no normalization or anything)*/
88 /*NB! Check carefully indexes here*/
90 for (n = 0; n < tsteps; n++)
92 gmx_fft_2d_real(fftp, GMX_FFT_REAL_TO_COMPLEX, intftab[0][n], ftspect1);
93 gmx_fft_2d_real(fftp, GMX_FFT_REAL_TO_COMPLEX, intftab[1][n], ftspect2);
95 /*Add to average for interface 1 here*/
96 addtoavgenergy(ftspect1, pspectavg1, rfl, tsteps);
97 /*Add to average for interface 2 here*/
98 addtoavgenergy(ftspect2, pspectavg2, rfl, tsteps);
100 /*Print out average energy-spectrum to outfiles[0] and outfiles[1];*/
102 datfile1 = gmx_ffopen(outfiles[0], "w");
103 datfile2 = gmx_ffopen(outfiles[1], "w");
105 /*Filling dat files with spectral data*/
106 fprintf(datfile1, "%s\n", "kx\t ky\t\tPower(kx,ky)");
107 fprintf(datfile2, "%s\n", "kx\t ky\t\tPower(kx,ky)");
108 for (n = 0; n < rfl; n++)
110 fprintf(datfile1, "%d\t%d\t %8.6f\n", (n / fy), (n % fy), pspectavg1[n]);
111 fprintf(datfile2, "%d\t%d\t %8.6f\n", (n /fy), (n % fy), pspectavg2[n]);
113 gmx_ffclose(datfile1);
114 gmx_ffclose(datfile2);
121 void powerspectavg_intf(t_interf ***if1, t_interf ***if2, int t, int xb, int yb, char **fnms)
131 for (n = 0; n < t; n++)
133 snew(surf[0][n], xy);
134 snew(surf[1][n], xy);
135 for (i = 0; i < xy; i++)
137 surf[0][n][i] = if1[n][i]->Z;
138 surf[1][n][i] = if2[n][i]->Z;
141 powerspectavg(surf, t, xb, yb, fnms);