2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2001
5 * BIOSON Research Institute, Dept. of Biophysical Chemistry
6 * University of Groningen, The Netherlands
7 * Copyright (c) 2012,2013, by the GROMACS development team, led by
8 * David van der Spoel, Berk Hess, Erik Lindahl, and including many
9 * others, as listed in the AUTHORS file in the top-level source
10 * directory and at http://www.gromacs.org.
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
44 #include <gmx_fatal.h>
47 #include "powerspect.h"
49 void addtoavgenergy(t_complex *list, real *result, int size, int tsteps){
52 result[i]+=cabs2(list[i])/tsteps;
58 void powerspectavg(real ***intftab, int tsteps, int xbins, int ybins, char **outfiles){
59 /*Fourier plans and output;*/
61 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 */
63 real *pspectavg1; /*power -spectrum 1st interface*/
64 real *pspectavg2; /* -------------- 2nd interface*/
66 FILE *datfile1,*datfile2; /*data-files with interface data*/
68 int fy=ybins/2+1; /* number of (symmetric) fourier y elements; */
69 int rfl=xbins*fy; /*length of real - DFT == Symmetric 2D matrix*/
72 /*Prepare data structures for FFT, with time averaging of power spectrum*/
73 if ( (status=gmx_fft_init_2d_real(&fftp,xbins,ybins,GMX_FFT_FLAG_NONE) )!=0)
76 gmx_fatal(status,__FILE__,__LINE__,"Error allocating FFT");
86 /*Fouriertransform directly (no normalization or anything)*/
87 /*NB! Check carefully indexes here*/
89 for (n=0;n<tsteps;n++){
90 gmx_fft_2d_real(fftp,GMX_FFT_REAL_TO_COMPLEX,intftab[0][n],ftspect1);
91 gmx_fft_2d_real(fftp,GMX_FFT_REAL_TO_COMPLEX,intftab[1][n],ftspect2);
93 /*Add to average for interface 1 here*/
94 addtoavgenergy(ftspect1,pspectavg1,rfl,tsteps);
95 /*Add to average for interface 2 here*/
96 addtoavgenergy(ftspect2,pspectavg2,rfl,tsteps);
98 /*Print out average energy-spectrum to outfiles[0] and outfiles[1];*/
100 datfile1 = ffopen(outfiles[0],"w");
101 datfile2 = ffopen(outfiles[1],"w");
103 /*Filling dat files with spectral data*/
104 fprintf(datfile1,"%s\n","kx\t ky\t\tPower(kx,ky)");
105 fprintf(datfile2,"%s\n","kx\t ky\t\tPower(kx,ky)");
107 fprintf(datfile1,"%d\t%d\t %8.6f\n",(n / fy),(n % fy),pspectavg1[n]);
108 fprintf(datfile2,"%d\t%d\t %8.6f\n",(n /fy),(n % fy),pspectavg2[n]);
118 void powerspectavg_intf(t_interf ***if1, t_interf ***if2, int t, int xb, int yb, char **fnms)
132 surf[0][n][i]=if1[n][i]->Z;
133 surf[1][n][i]=if2[n][i]->Z;
136 powerspectavg(surf, t,xb,yb,fnms);