applications/staticProps/Hxy.cpp

/*
 * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
 *
 * The University of Notre Dame grants you ("Licensee") a
 * non-exclusive, royalty free, license to use, modify and
 * redistribute this software in source and binary code form, provided
 * that the following conditions are met:
 *
 * 1. Acknowledgement of the program authors must be made in any
 *    publication of scientific results based in part on use of the
 *    program.  An acceptable form of acknowledgement is citation of
 *    the article in which the program was described (Matthew
 *    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
 *    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
 *    Parallel Simulation Engine for Molecular Dynamics,"
 *    J. Comput. Chem. 26, pp. 252-271 (2005))
 *
 * 2. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 *
 * 3. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the
 *    distribution.
 *
 * This software is provided "AS IS," without a warranty of any
 * kind. All express or implied conditions, representations and
 * warranties, including any implied warranty of merchantability,
 * fitness for a particular purpose or non-infringement, are hereby
 * excluded.  The University of Notre Dame and its licensors shall not
 * be liable for any damages suffered by licensee as a result of
 * using, modifying or distributing the software or its
 * derivatives. In no event will the University of Notre Dame or its
 * licensors be liable for any lost revenue, profit or data, or for
 * direct, indirect, special, consequential, incidental or punitive
 * damages, however caused and regardless of the theory of liability,
 * arising out of the use of or inability to use software, even if the
 * University of Notre Dame has been advised of the possibility of
 * such damages.
 *
 *
 *  Hxy.cpp
 *  OOPSE-2.0
 *
 *  Created by Xiuquan Sun on 05/09/06.
 *  @author  Xiuquan Sun 
 *  @version $Id: Hxy.cpp,v 1.1 2006-05-12 21:34:43 xsun Exp $
 *
 */

/* Calculates the undulation spectrum of the lipid membrance. */

#include <algorithm>
#include <fstream>
#include "applications/staticProps/Hxy.hpp"
#include "utils/simError.h"
#include "io/DumpReader.hpp"
#include "primitives/Molecule.hpp"
#include<stdio.h>
#include<string.h>
#include<stdlib.h>
#include<math.h>
#include<fftw3.h>
#include<mkl_lapack64.h>

namespace oopse {
  
  Hxy::Hxy(SimInfo* info, const std::string& filename, const std::string& sele, int nbins_x, int nbins_y, int nrbins)
    : StaticAnalyser(info, filename), selectionScript_(sele),  evaluator_(info), seleMan_(info), nBinsX_(nbins_x), nBinsY_(nbins_y), nbins_(nrbins){

    evaluator_.loadScriptString(sele);
    if (!evaluator_.isDynamic()) {
      seleMan_.setSelectionSet(evaluator_.evaluate());
    }

    gridsample_.resize(nBinsX_*nBinsY_);
    gridZ_.resize(nBinsX_*nBinsY_);

    sum_bin.resize(nbins);
    avg_bin.resize(nbins);
    errbin_sum.resize(nbins);
    errbin.resize(nbins);
    sum_bin_sq.resize(nbins);
    avg_bin_sq.resize(nbins);
    errbin_sum_sq.resize(nbins);
    errbin_sq.resize(nbins);

    setOutputName(getPrefix(filename) + ".Hxy");
  }

  void Hxy::process() {
    DumpReader reader(info_, dumpFilename_);    
    int nFrames = reader.getNFrames();
    nProcessed_ = nFrames/step_;

    std::vector<double> mag, newmag;
    double lenX_, lenY_;
    double gridX_, gridY_;
    double halfBoxX_, halfBoxY_;
    int binNoX, binNoY;
    double interpsum, value;
    int ninterp, px, py, newp;
    int newx, newy, newindex, index;
    int new_i, new_j, new_index;
    double freq_x, freq_y, zero_freq_x, zero_freq_y, freq;
    double maxfreqx, maxfreqy, maxfreq, dfreq;
    int whichbin;
    
    for (int istep = 0; istep < nFrames; istep += step_) {
      
      int nMolecules = reader.getNMolecules();

      fftw_complex *in, *out;
      fftw_plan p;

      in = fftw_malloc(sizeof(fftw_complex) * (nBinsX_*nBinsY_));
      out = fftw_malloc(sizeof(fftw_complex) *(nBinsX_*nBinsY_));
      p =  fftw_plan_dft_2d(nBinsX_, 
                            nBinsY_, 
                            in, out,
                            FFTW_FORWARD, 
                            FFTW_ESTIMATE);
      
      int i, j;   

      for(i=0; i < nBinsX_*nBinsY_; i++){
        gridsample_[i].clear();
        gridZ_[i].clear();
      }

      for(i=0; i < nbins; i++){
        sum_bin[i].clear();
        avg_bin[i].clear();
        errbin_sum[i].clear();
        errbin[i].clear();
        sum_bin_sq[i].clear();
        avg_bin_sq[i].clear();
        errbin_sum_sq[i].clear();
        errbin_sq[i].clear();
      }

      mag.resize(nBinsX_*nBinsY_);
      newmag.resize(nBinsX_*nBinsY_);

      for(i=0; i < nBinsX_*nBinsY_; i++){
        mag[i].clear();
        newmag[i].clear();
      }
      
      StuntDouble* sd;
      reader.readFrame(istep);
      currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
        
      Mat3x3d hmat = currentSnapshot_->getHmat();

      lenX_ = hmat(0,0);
      lenY_ = hmat(1,1);

      gridX_ = lenX_ /(nBinsX_);
      gridY_ = lenY_ /(nBinsY_);

      double halfBoxX_ = lenX_ / 2.0;      
      double halfBoxY_ = lenY_ / 2.0;      
        
      if (evaluator_.isDynamic()) {
        seleMan_.setSelectionSet(evaluator_.evaluate());
      }
      
      //wrap the stuntdoubles into a cell     
      for (sd = seleMan_.beginSelected(i); sd != NULL; sd = seleMan_.nextSelected(i)) {
        Vector3d pos = sd->getPos();
        currentSnapshot_->wrapVector(pos);
        sd->setPos(pos);
      }
      
      //determine which atom belongs to which grid
      for (sd = seleMan_.beginSelected(i); sd != NULL; sd = seleMan_.nextSelected(i)) {
        Vector3d pos = sd->getPos();
        //int binNo = (pos.z() /deltaR_) - 1;
        int binNoX = (pos.x() + halfBoxX_) /gridX_; 
        int binNoY = (pos.y() + halfBoxY_) /gridY_;
        //std::cout << "pos.z = " << pos.z() << " halfBoxZ_ = " << halfBoxZ_ << " deltaR_ = "  << deltaR_ << " binNo = " << binNo << "\n";
        gridZ_[binNoX*nBinsY_+binNoY] += pos.z();
        gridsample_[binNoX*nBinsY_+binNoY]++;
      }
      
      // FFT stuff depends on nx and ny, so delay allocation until we have
      // that information
      
      for (i=0; i< nBinsX_; i++) {
        for(j=0; j< nBinsY_; j++) {
          newindex = i*nBinsY_ + j;
          mag[newindex] = 0.0;
          newmag[newindex] = 0.0;
        }
      }
      
      for(i = 0; i < nBinsX_; i++){
        for(j = 0; j < nBinsY_; j++){
          newindex = i * nBinsY_ + j;
          if(gridsample_[newindex] > 0){
            gridZ_[newindex] = gridZ_[newindex] / (double)gridsample_[newindex];
          }
        }
      }
        
      for (i=0; i< nBinsX_; i++) {
        for(j=0; j< nBinsY_; j++) {
          newindex = i*nBinsY_ + j;
          if (gridsample_[newindex] == 0) {
            // interpolate from surrounding points:
            
            interpsum = 0.0;
            ninterp = 0;
            
            //point1 = bottom;
            
            px = i;
            py = j - 1;
            newp = px*nBinsY_ + py;
            if ((py >= 0) && (gridsample_[newp] > 0)) {
              interpsum += gridZ_[newp];
              ninterp++;
            } 
            
            //point2 = top;
            
            px = i;
            py = j + 1;
            newp = px*nBinsY_ + py;
            if ((py < nBinsY_) && (gridsample_[newp] > 0)) {
              interpsum += gridZ_[newp];
              ninterp++;
            } 
            
            //point3 = left;
            
            px = i - 1;
            py = j;
            newp = px*nBinsY_ + py;
            if ((px >= 0) && (gridsample_[newp] > 0)) {
              interpsum += gridZ_[newp];
              ninterp++;
            }
            
            //point4 = right;
            
            px = i + 1;
            py = j;
            newp = px*nBinsY_ + py;
            if ((px < nBinsX_) && (gridsample_[newp] > 0)) {
              interpsum += gridZ_[newp];
              ninterp++;
            } 
            
            value = interpsum / (double)ninterp;
            
            gridZ_[newindex] = value;
          }
        }
      }
      
      for (i=0; i < nBinsX_; i++) {
        for (j=0; j < nBinsY_; j++) {
          newindex = i*nBinsY_ + j;
          
          c_re(in[newindex]) = gridZ_[newindex];
          c_im(in[newindex]) = 0.0;
        } 
      }
      
      fftw_execute(p); 
      
      for (i=0; i< nBinsX_; i++) {
        for(j=0; j< nBinsY_; j++) {
          newindex = i*nBinsY_ + j;
          mag[newindex] = pow(c_re(out[newindex]),2) + pow(c_im(out[newindex]),2);
        }
      }
      
      fftw_destroy_plan(p);
      fftw_free(out);
      fftw_free(in);

      for (i=0; i< (nBinsX_/2); i++) {
        for(j=0; j< (nBinsY_/2); j++) {
          index = i*nBinsY_ + j;
          new_i = i + (nBinsX_/2);
          new_j = j + (nBinsY_/2);
          new_index = new_i*nBinsY_ + new_j;
          newmag[new_index] = mag[index];
        }
      }
      
      for (i=(nBinsX_/2); i< nBinsX_; i++) {
        for(j=0; j< (nBinsY_/2); j++) {
          index = i*nBinsY_ + j;
          new_i = i - (nBinsX_/2);
          new_j = j + (nBinsY_/2);
          new_index = new_i*nBinsY_ + new_j;
          newmag[new_index] = mag[index];
        }
      }
      
      for (i=0; i< (nBinsX_/2); i++) {
        for(j=(nBinsY_/2); j< nBinsY_; j++) {
          index = i*nBinsY_ + j;
          new_i = i + (nBinsX_/2);
          new_j = j - (nBinsY_/2);
          new_index = new_i*nBinsY_ + new_j;
          newmag[new_index] = mag[index];
        }
      }
      
      for (i=(nBinsX_/2); i< nBinsX_; i++) {
        for(j=(nBinsY_/2); j< nBinsY_; j++) {
          index = i*nBinsY_ + j;
          new_i = i - (nBinsX_/2);
          new_j = j - (nBinsY_/2);
          new_index = new_i*nBinsY_ + new_j;
          newmag[new_index] = mag[index];
        }
      }
    
      maxfreqx = 1.0 / gridX_;
      maxfreqy = 1.0 / gridY_;
      
      //  printf("%lf\t%lf\t%lf\t%lf\n", dx, dy, maxfreqx, maxfreqy);
      
      maxfreq = sqrt(maxfreqx*maxfreqx + maxfreqy*maxfreqy);
      dfreq = maxfreq/(double)(nbins-1);
    
      //printf("%lf\n", dfreq);
      
      zero_freq_x = nBinsX_/2; 
      zero_freq_y = nBinsY_/2; 
      
      for (i=0; i< nBinsX_; i++) {
        for(j=0; j< nBinsY_; j++) {
          
          freq_x = (double)(i - zero_freq_x)*maxfreqx*2 / nBinsX_;
          freq_y = (double)(j - zero_freq_y)*maxfreqy*2 / nBinsY_;
          
          freq = sqrt(freq_x*freq_x + freq_y*freq_y);
          
          whichbin = (int) (freq / dfreq);
          newindex = i*nBinsY_ + j;
          //    printf("%d %d %lf %lf\n", whichbin, newindex, freq, dfreq);
          bin[whichbin][istep] += newmag[newindex];
          samples[whichbin][istep]++;
        }
      }
      
      for ( i = 0; i < nbins; i++) {
        if ( samples[i][istep] > 0) {
          bin[i][istep] = 4.0 * sqrt(bin[i][istep] / (double)samples[i][istep]) / (double)nMolecules;
        }
      }    
      
    }
  
    for (i = 0; i < nbins; i++) {
      for (j = 0; j < nFrames; j++) {
        sum_bin[i] += bin[i][j];
        sum_bin_sq[i] += bin[i][j] * bin[i][j];
      }
      avg_bin[i] = sum_bin[i] / (double)nFrames;
      avg_bin_sq[i] = sum_bin_sq[i] / (double)nFrames;
      for (j = 0; j < nFrames; j++) {
        errbin_sum[i] += pow((bin[i][j] - avg_bin[i]), 2);
        errbin_sum_sq[i] += pow((bin[i][j] * bin[i][j] - avg_bin_sq[i]), 2);
      }
      errbin[i] = sqrt( errbin_sum[i] / (double)nFrames );
      errbin_sq[i] = sqrt( errbin_sum_sq[i] / (double)nFrames );
    }

    printSpectrum();
  }
    
  void Hxy::printSpectrum() {
    std::ofstream rdfStream(outputFilename_.c_str());
    if (rdfStream.is_open()) {
      
      for (int i = 0; i < nbins; i++) {
        if ( avg_bin[i] > 0 ){
          rdfStream << i*dfreq << "\t"
                    <<pow(avg_bin[i], 2)<<"\t"
                    <<errbin_sq[i]<<"\t"
                    <<avg_bin[i]<<"\t"
                    <<errbin[i]<<"\n";
        }
      }
    } else {
      
      sprintf(painCave.errMsg, "Hxy: unable to open %s\n", outputFilename_.c_str());
      painCave.isFatal = 1;
      simError();  
    }
    
    rdfStream.close();
  }
  
}
Revision:	2751
Committed:	Fri May 12 21:34:43 2006 UTC (18 years, 3 months ago) by xsun
File size:	11430 byte(s)
Log Message:	Changes to calculate undulation spectrum
#	User	Rev	Content
1	xsun	2751	/*
2			* Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3			*
4			* The University of Notre Dame grants you ("Licensee") a
5			* non-exclusive, royalty free, license to use, modify and
6			* redistribute this software in source and binary code form, provided
7			* that the following conditions are met:
8			*
9			* 1. Acknowledgement of the program authors must be made in any
10			* publication of scientific results based in part on use of the
11			* program. An acceptable form of acknowledgement is citation of
12			* the article in which the program was described (Matthew
13			* A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14			* J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15			* Parallel Simulation Engine for Molecular Dynamics,"
16			* J. Comput. Chem. 26, pp. 252-271 (2005))
17			*
18			* 2. Redistributions of source code must retain the above copyright
19			* notice, this list of conditions and the following disclaimer.
20			*
21			* 3. Redistributions in binary form must reproduce the above copyright
22			* notice, this list of conditions and the following disclaimer in the
23			* documentation and/or other materials provided with the
24			* distribution.
25			*
26			* This software is provided "AS IS," without a warranty of any
27			* kind. All express or implied conditions, representations and
28			* warranties, including any implied warranty of merchantability,
29			* fitness for a particular purpose or non-infringement, are hereby
30			* excluded. The University of Notre Dame and its licensors shall not
31			* be liable for any damages suffered by licensee as a result of
32			* using, modifying or distributing the software or its
33			* derivatives. In no event will the University of Notre Dame or its
34			* licensors be liable for any lost revenue, profit or data, or for
35			* direct, indirect, special, consequential, incidental or punitive
36			* damages, however caused and regardless of the theory of liability,
37			* arising out of the use of or inability to use software, even if the
38			* University of Notre Dame has been advised of the possibility of
39			* such damages.
40			*
41			*
42			* Hxy.cpp
43			* OOPSE-2.0
44			*
45			* Created by Xiuquan Sun on 05/09/06.
46			* @author Xiuquan Sun
47			* @version $Id: Hxy.cpp,v 1.1 2006-05-12 21:34:43 xsun Exp $
48			*
49			*/
50
51			/* Calculates the undulation spectrum of the lipid membrance. */
52
53			#include <algorithm>
54			#include <fstream>
55			#include "applications/staticProps/Hxy.hpp"
56			#include "utils/simError.h"
57			#include "io/DumpReader.hpp"
58			#include "primitives/Molecule.hpp"
59			#include<stdio.h>
60			#include<string.h>
61			#include<stdlib.h>
62			#include<math.h>
63			#include<fftw3.h>
64			#include<mkl_lapack64.h>
65
66			namespace oopse {
67
68			Hxy::Hxy(SimInfo* info, const std::string& filename, const std::string& sele, int nbins_x, int nbins_y, int nrbins)
69			: StaticAnalyser(info, filename), selectionScript_(sele), evaluator_(info), seleMan_(info), nBinsX_(nbins_x), nBinsY_(nbins_y), nbins_(nrbins){
70
71			evaluator_.loadScriptString(sele);
72			if (!evaluator_.isDynamic()) {
73			seleMan_.setSelectionSet(evaluator_.evaluate());
74			}
75
76			gridsample_.resize(nBinsX_*nBinsY_);
77			gridZ_.resize(nBinsX_*nBinsY_);
78
79			sum_bin.resize(nbins);
80			avg_bin.resize(nbins);
81			errbin_sum.resize(nbins);
82			errbin.resize(nbins);
83			sum_bin_sq.resize(nbins);
84			avg_bin_sq.resize(nbins);
85			errbin_sum_sq.resize(nbins);
86			errbin_sq.resize(nbins);
87
88			setOutputName(getPrefix(filename) + ".Hxy");
89			}
90
91			void Hxy::process() {
92			DumpReader reader(info_, dumpFilename_);
93			int nFrames = reader.getNFrames();
94			nProcessed_ = nFrames/step_;
95
96			std::vector<double> mag, newmag;
97			double lenX_, lenY_;
98			double gridX_, gridY_;
99			double halfBoxX_, halfBoxY_;
100			int binNoX, binNoY;
101			double interpsum, value;
102			int ninterp, px, py, newp;
103			int newx, newy, newindex, index;
104			int new_i, new_j, new_index;
105			double freq_x, freq_y, zero_freq_x, zero_freq_y, freq;
106			double maxfreqx, maxfreqy, maxfreq, dfreq;
107			int whichbin;
108
109			for (int istep = 0; istep < nFrames; istep += step_) {
110
111			int nMolecules = reader.getNMolecules();
112
113			fftw_complex in, out;
114			fftw_plan p;
115
116			in = fftw_malloc(sizeof(fftw_complex) * (nBinsX_*nBinsY_));
117			out = fftw_malloc(sizeof(fftw_complex) (nBinsX_nBinsY_));
118			p = fftw_plan_dft_2d(nBinsX_,
119			nBinsY_,
120			in, out,
121			FFTW_FORWARD,
122			FFTW_ESTIMATE);
123
124			int i, j;
125
126			for(i=0; i < nBinsX_*nBinsY_; i++){
127			gridsample_[i].clear();
128			gridZ_[i].clear();
129			}
130
131			for(i=0; i < nbins; i++){
132			sum_bin[i].clear();
133			avg_bin[i].clear();
134			errbin_sum[i].clear();
135			errbin[i].clear();
136			sum_bin_sq[i].clear();
137			avg_bin_sq[i].clear();
138			errbin_sum_sq[i].clear();
139			errbin_sq[i].clear();
140			}
141
142			mag.resize(nBinsX_*nBinsY_);
143			newmag.resize(nBinsX_*nBinsY_);
144
145			for(i=0; i < nBinsX_*nBinsY_; i++){
146			mag[i].clear();
147			newmag[i].clear();
148			}
149
150			StuntDouble* sd;
151			reader.readFrame(istep);
152			currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
153
154			Mat3x3d hmat = currentSnapshot_->getHmat();
155
156			lenX_ = hmat(0,0);
157			lenY_ = hmat(1,1);
158
159			gridX_ = lenX_ /(nBinsX_);
160			gridY_ = lenY_ /(nBinsY_);
161
162			double halfBoxX_ = lenX_ / 2.0;
163			double halfBoxY_ = lenY_ / 2.0;
164
165			if (evaluator_.isDynamic()) {
166			seleMan_.setSelectionSet(evaluator_.evaluate());
167			}
168
169			//wrap the stuntdoubles into a cell
170			for (sd = seleMan_.beginSelected(i); sd != NULL; sd = seleMan_.nextSelected(i)) {
171			Vector3d pos = sd->getPos();
172			currentSnapshot_->wrapVector(pos);
173			sd->setPos(pos);
174			}
175
176			//determine which atom belongs to which grid
177			for (sd = seleMan_.beginSelected(i); sd != NULL; sd = seleMan_.nextSelected(i)) {
178			Vector3d pos = sd->getPos();
179			//int binNo = (pos.z() /deltaR_) - 1;
180			int binNoX = (pos.x() + halfBoxX_) /gridX_;
181			int binNoY = (pos.y() + halfBoxY_) /gridY_;
182			//std::cout << "pos.z = " << pos.z() << " halfBoxZ_ = " << halfBoxZ_ << " deltaR_ = " << deltaR_ << " binNo = " << binNo << "\n";
183			gridZ_[binNoX*nBinsY_+binNoY] += pos.z();
184			gridsample_[binNoX*nBinsY_+binNoY]++;
185			}
186
187			// FFT stuff depends on nx and ny, so delay allocation until we have
188			// that information
189
190			for (i=0; i< nBinsX_; i++) {
191			for(j=0; j< nBinsY_; j++) {
192			newindex = i*nBinsY_ + j;
193			mag[newindex] = 0.0;
194			newmag[newindex] = 0.0;
195			}
196			}
197
198			for(i = 0; i < nBinsX_; i++){
199			for(j = 0; j < nBinsY_; j++){
200			newindex = i * nBinsY_ + j;
201			if(gridsample_[newindex] > 0){
202			gridZ_[newindex] = gridZ_[newindex] / (double)gridsample_[newindex];
203			}
204			}
205			}
206
207			for (i=0; i< nBinsX_; i++) {
208			for(j=0; j< nBinsY_; j++) {
209			newindex = i*nBinsY_ + j;
210			if (gridsample_[newindex] == 0) {
211			// interpolate from surrounding points:
212
213			interpsum = 0.0;
214			ninterp = 0;
215
216			//point1 = bottom;
217
218			px = i;
219			py = j - 1;
220			newp = px*nBinsY_ + py;
221			if ((py >= 0) && (gridsample_[newp] > 0)) {
222			interpsum += gridZ_[newp];
223			ninterp++;
224			}
225
226			//point2 = top;
227
228			px = i;
229			py = j + 1;
230			newp = px*nBinsY_ + py;
231			if ((py < nBinsY_) && (gridsample_[newp] > 0)) {
232			interpsum += gridZ_[newp];
233			ninterp++;
234			}
235
236			//point3 = left;
237
238			px = i - 1;
239			py = j;
240			newp = px*nBinsY_ + py;
241			if ((px >= 0) && (gridsample_[newp] > 0)) {
242			interpsum += gridZ_[newp];
243			ninterp++;
244			}
245
246			//point4 = right;
247
248			px = i + 1;
249			py = j;
250			newp = px*nBinsY_ + py;
251			if ((px < nBinsX_) && (gridsample_[newp] > 0)) {
252			interpsum += gridZ_[newp];
253			ninterp++;
254			}
255
256			value = interpsum / (double)ninterp;
257
258			gridZ_[newindex] = value;
259			}
260			}
261			}
262
263			for (i=0; i < nBinsX_; i++) {
264			for (j=0; j < nBinsY_; j++) {
265			newindex = i*nBinsY_ + j;
266
267			c_re(in[newindex]) = gridZ_[newindex];
268			c_im(in[newindex]) = 0.0;
269			}
270			}
271
272			fftw_execute(p);
273
274			for (i=0; i< nBinsX_; i++) {
275			for(j=0; j< nBinsY_; j++) {
276			newindex = i*nBinsY_ + j;
277			mag[newindex] = pow(c_re(out[newindex]),2) + pow(c_im(out[newindex]),2);
278			}
279			}
280
281			fftw_destroy_plan(p);
282			fftw_free(out);
283			fftw_free(in);
284
285			for (i=0; i< (nBinsX_/2); i++) {
286			for(j=0; j< (nBinsY_/2); j++) {
287			index = i*nBinsY_ + j;
288			new_i = i + (nBinsX_/2);
289			new_j = j + (nBinsY_/2);
290			new_index = new_i*nBinsY_ + new_j;
291			newmag[new_index] = mag[index];
292			}
293			}
294
295			for (i=(nBinsX_/2); i< nBinsX_; i++) {
296			for(j=0; j< (nBinsY_/2); j++) {
297			index = i*nBinsY_ + j;
298			new_i = i - (nBinsX_/2);
299			new_j = j + (nBinsY_/2);
300			new_index = new_i*nBinsY_ + new_j;
301			newmag[new_index] = mag[index];
302			}
303			}
304
305			for (i=0; i< (nBinsX_/2); i++) {
306			for(j=(nBinsY_/2); j< nBinsY_; j++) {
307			index = i*nBinsY_ + j;
308			new_i = i + (nBinsX_/2);
309			new_j = j - (nBinsY_/2);
310			new_index = new_i*nBinsY_ + new_j;
311			newmag[new_index] = mag[index];
312			}
313			}
314
315			for (i=(nBinsX_/2); i< nBinsX_; i++) {
316			for(j=(nBinsY_/2); j< nBinsY_; j++) {
317			index = i*nBinsY_ + j;
318			new_i = i - (nBinsX_/2);
319			new_j = j - (nBinsY_/2);
320			new_index = new_i*nBinsY_ + new_j;
321			newmag[new_index] = mag[index];
322			}
323			}
324
325			maxfreqx = 1.0 / gridX_;
326			maxfreqy = 1.0 / gridY_;
327
328			// printf("%lf\t%lf\t%lf\t%lf\n", dx, dy, maxfreqx, maxfreqy);
329
330			maxfreq = sqrt(maxfreqxmaxfreqx + maxfreqymaxfreqy);
331			dfreq = maxfreq/(double)(nbins-1);
332
333			//printf("%lf\n", dfreq);
334
335			zero_freq_x = nBinsX_/2;
336			zero_freq_y = nBinsY_/2;
337
338			for (i=0; i< nBinsX_; i++) {
339			for(j=0; j< nBinsY_; j++) {
340
341			freq_x = (double)(i - zero_freq_x)maxfreqx2 / nBinsX_;
342			freq_y = (double)(j - zero_freq_y)maxfreqy2 / nBinsY_;
343
344			freq = sqrt(freq_xfreq_x + freq_yfreq_y);
345
346			whichbin = (int) (freq / dfreq);
347			newindex = i*nBinsY_ + j;
348			// printf("%d %d %lf %lf\n", whichbin, newindex, freq, dfreq);
349			bin[whichbin][istep] += newmag[newindex];
350			samples[whichbin][istep]++;
351			}
352			}
353
354			for ( i = 0; i < nbins; i++) {
355			if ( samples[i][istep] > 0) {
356			bin[i][istep] = 4.0 * sqrt(bin[i][istep] / (double)samples[i][istep]) / (double)nMolecules;
357			}
358			}
359
360			}
361
362			for (i = 0; i < nbins; i++) {
363			for (j = 0; j < nFrames; j++) {
364			sum_bin[i] += bin[i][j];
365			sum_bin_sq[i] += bin[i][j] * bin[i][j];
366			}
367			avg_bin[i] = sum_bin[i] / (double)nFrames;
368			avg_bin_sq[i] = sum_bin_sq[i] / (double)nFrames;
369			for (j = 0; j < nFrames; j++) {
370			errbin_sum[i] += pow((bin[i][j] - avg_bin[i]), 2);
371			errbin_sum_sq[i] += pow((bin[i][j] * bin[i][j] - avg_bin_sq[i]), 2);
372			}
373			errbin[i] = sqrt( errbin_sum[i] / (double)nFrames );
374			errbin_sq[i] = sqrt( errbin_sum_sq[i] / (double)nFrames );
375			}
376
377			printSpectrum();
378			}
379
380			void Hxy::printSpectrum() {
381			std::ofstream rdfStream(outputFilename_.c_str());
382			if (rdfStream.is_open()) {
383
384			for (int i = 0; i < nbins; i++) {
385			if ( avg_bin[i] > 0 ){
386			rdfStream << i*dfreq << "\t"
387			<<pow(avg_bin[i], 2)<<"\t"
388			<<errbin_sq[i]<<"\t"
389			<<avg_bin[i]<<"\t"
390			<<errbin[i]<<"\n";
391			}
392			}
393			} else {
394
395			sprintf(painCave.errMsg, "Hxy: unable to open %s\n", outputFilename_.c_str());
396			painCave.isFatal = 1;
397			simError();
398			}
399
400			rdfStream.close();
401			}
402
403			}