src/tools/powerspect.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2001
   5  * BIOSON Research Institute, Dept. of Biophysical Chemistry
   6  * University of Groningen, The Netherlands
   7  * Copyright (c) 2012, 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.
  11  *
  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.
  16  *
  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.
  21  *
  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.
  26  *
  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.
  34  *
  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.
  37  */
  38
  39 #ifdef HAVE_CONFIG_H
  40 #include <config.h>
  41 #endif
  42 #include <smalloc.h>
  43 #include <gmx_fft.h>
  44 #include <gmx_fatal.h>
  45 #include <futil.h>
  46 #include "interf.h"
  47 #include "powerspect.h"
  48
  49 void addtoavgenergy(t_complex *list, real *result, int size, int tsteps){
  50   int i;
  51         for (i=0;i<size;i++){
  52                 result[i]+=cabs2(list[i])/tsteps;
  53         }
  54
  55 }
  56
  57
  58 void powerspectavg(real ***intftab, int tsteps, int xbins, int ybins, char **outfiles){
  59   /*Fourier plans and output;*/
  60   gmx_fft_t  fftp;
  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 */
  62   t_complex *ftspect2;
  63   real *pspectavg1; /*power -spectrum 1st interface*/
  64   real *pspectavg2; /* -------------- 2nd interface*/
  65   real *temp;
  66   FILE *datfile1,*datfile2; /*data-files with interface data*/
  67   int n; /*time index*/
  68   int fy=ybins/2+1; /* number of (symmetric) fourier y elements; */
  69   int rfl=xbins*fy; /*length of real - DFT == Symmetric 2D matrix*/
  70   int status;
  71
  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)
  74   {
  75         free(fftp);
  76         gmx_fatal(status,__FILE__,__LINE__,"Error allocating FFT");
  77     }
  78
  79 /*Initialize arrays*/
  80 snew(ftspect1,rfl);
  81 snew(ftspect2,rfl);
  82 snew(temp,rfl);
  83 snew(pspectavg1,rfl);
  84 snew(pspectavg2,rfl);
  85
  86 /*Fouriertransform directly (no normalization or anything)*/
  87 /*NB! Check carefully indexes here*/
  88
  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);
  92
  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);
  97         }
  98         /*Print out average energy-spectrum to outfiles[0] and outfiles[1];*/
  99
 100 datfile1 = ffopen(outfiles[0],"w");
 101 datfile2 = ffopen(outfiles[1],"w");
 102
 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)");
 106         for(n=0;n<rfl;n++){
 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]);
 109         }
 110 ffclose(datfile1);
 111 ffclose(datfile2);
 112
 113 free(ftspect1);
 114 free(ftspect2);
 115
 116 }
 117
 118 void powerspectavg_intf(t_interf ***if1, t_interf ***if2, int t, int xb, int yb, char **fnms)
 119 {
 120         real ***surf;
 121
 122         int xy=xb*yb;
 123         int i,n;
 124
 125         snew(surf,2);
 126         snew(surf[0],t);
 127         snew(surf[1],t);
 128         for (n=0;n<t;n++){
 129                 snew(surf[0][n],xy);
 130                 snew(surf[1][n],xy);
 131                 for(i=0;i<xy;i++){
 132                         surf[0][n][i]=if1[n][i]->Z;
 133                         surf[1][n][i]=if2[n][i]->Z;
 134                 }
 135         }
 136         powerspectavg(surf, t,xb,yb,fnms);
 137 }
 138