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. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 *
 * 2. 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.
 *
 * SUPPORT OPEN SCIENCE!  If you use OpenMD or its source code in your
 * research, please cite the appropriate papers when you publish your
 * work.  Good starting points are:
 *                                                                      
 * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).             
 * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).          
 * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).          
 * [4] Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
 * [4] , Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). *
 *
 *  Created by Xiuquan Sun on 05/09/06.
 *  @author  Xiuquan Sun 
 *  @version $Id$
 *
 */

/* 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>

namespace OpenMD {
  
  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_);
    mag.resize(nBinsX_*nBinsY_);
    newmag.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_);

    bin.resize(nbins_);
    samples.resize(nbins_);

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

  Hxy::~Hxy(){
      gridsample_.clear();
      gridZ_.clear();
      sum_bin.clear();
      avg_bin.clear();
      errbin_sum.clear();
      errbin.clear();
      sum_bin_sq.clear();
      avg_bin_sq.clear();
      errbin_sum_sq.clear();
      errbin_sq.clear();
      
      for(unsigned int i=0; i < bin.size(); i++)
        bin[i].clear();
      for(unsigned int i=0; i < samples.size(); i++)
        samples[i].clear();

      mag.clear();
      newmag.clear();
  }

  void Hxy::process() {
#if defined(HAVE_FFTW_H) || defined(HAVE_DFFTW_H) || defined(HAVE_FFTW3_H)
    DumpReader reader(info_, dumpFilename_);    
    int nFrames = reader.getNFrames();
    nProcessed_ = nFrames/step_;
    
    for(unsigned int k=0; k < bin.size(); k++)
      bin[k].resize(nFrames);
    for(unsigned int k=0; k < samples.size(); k++)
      samples[k].resize(nFrames);

    RealType lenX_, lenY_;
    RealType gridX_, gridY_;
    RealType halfBoxX_, halfBoxY_;

    RealType interpsum, value;
    int ninterp, px, py, newp;
    int newindex, index;
    int new_i, new_j, new_index;

    RealType freq_x, freq_y, zero_freq_x, zero_freq_y, freq;
    RealType maxfreqx, maxfreqy, maxfreq;

    int whichbin;

    std::fill(sum_bin.begin(), sum_bin.end(), 0.0);
    std::fill(avg_bin.begin(), avg_bin.end(), 0.0);
    std::fill(errbin_sum.begin(), errbin_sum.end(), 0.0);
    std::fill(errbin.begin(), errbin.end(), 0.0);
    std::fill(sum_bin_sq.begin(), sum_bin_sq.end(), 0.0);
    std::fill(avg_bin_sq.begin(), avg_bin_sq.end(), 0.0);
    std::fill(errbin_sum_sq.begin(), errbin_sum_sq.end(), 0.0);
    std::fill(errbin_sq.begin(), errbin_sq.end(), 0.0);
    
    for(unsigned int i=0; i < bin.size(); i++)
      std::fill(bin[i].begin(), bin[i].end(), 0.0);
    
    for(unsigned int i=0; i < samples.size(); i++)
      std::fill(samples[i].begin(), samples[i].end(), 0);
    
    for (int istep = 0; istep < nFrames; istep += step_) {
      
      reader.readFrame(istep);
      currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
      
      Mat3x3d hmat = currentSnapshot_->getHmat();
      
#ifdef HAVE_FFTW3_H
      fftw_plan p;
#else
      fftwnd_plan p;
#endif
      fftw_complex *in, *out;
      
      in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * (nBinsX_*nBinsY_));
      out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) *(nBinsX_*nBinsY_));

#ifdef HAVE_FFTW3_H
      p = fftw_plan_dft_2d(nBinsX_, nBinsY_, in, out, 
                           FFTW_FORWARD, FFTW_ESTIMATE); 
#else
      p = fftw2d_create_plan(nBinsX_, nBinsY_, FFTW_FORWARD, FFTW_ESTIMATE);
#endif

      std::fill(gridsample_.begin(), gridsample_.end(), 0);
      std::fill(gridZ_.begin(), gridZ_.end(), 0.0);
      std::fill(mag.begin(), mag.end(), 0.0);
      std::fill(newmag.begin(), newmag.end(), 0.0);

      int i, j;   
      
      StuntDouble* sd;
      
      lenX_ = hmat(0,0);
      lenY_ = hmat(1,1);
      
      gridX_ = lenX_ /(nBinsX_);
      gridY_ = lenY_ /(nBinsY_);
      
      halfBoxX_ = lenX_ / 2.0;      
      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();
        if (usePeriodicBoundaryConditions_)
          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 = (int) ((pos.x() + halfBoxX_) / gridX_);
        int binNoY = (int) ((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;
          if(gridsample_[newindex] > 0){
            gridZ_[newindex] = gridZ_[newindex] / (RealType)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 / (RealType)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;
        } 
      }

#ifdef HAVE_FFTW3_H
      fftw_execute(p);
#else
      fftwnd_one(p, in, out);
#endif
      
      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);
        }
      }

#ifdef HAVE_FFTW3_H
      fftw_destroy_plan(p);
#else
      fftwnd_destroy_plan(p);
#endif      
      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/(RealType)(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 = (RealType)(i - zero_freq_x)*maxfreqx*2 / nBinsX_;
          freq_y = (RealType)(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] / (RealType)samples[i][istep]) / (RealType)nBinsX_ / (RealType)nBinsY_;
        }
      }    
    }

    for (int i = 0; i < nbins_; i++) {
      for (int 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] / (RealType)nFrames;
      avg_bin_sq[i] = sum_bin_sq[i] / (RealType)nFrames;
      for (int 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] / (RealType)nFrames );
      errbin_sq[i] = sqrt( errbin_sum_sq[i] / (RealType)nFrames );
    }

    printSpectrum();

#else
    sprintf(painCave.errMsg, "Hxy: FFTW support was not compiled in!\n");
    painCave.isFatal = 1;
    simError();  

#endif
  }
    
  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 << (RealType)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:	1796
Committed:	Mon Sep 10 18:38:44 2012 UTC (12 years, 7 months ago) by gezelter
File size:	12972 byte(s)
Log Message:	Updating MPI calls to C++, fixing a DumpWriter bug, cleaning cruft
#	User	Rev	Content
1	xsun	955	/*
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	gezelter	1390	* 1. Redistributions of source code must retain the above copyright
10	xsun	955	* notice, this list of conditions and the following disclaimer.
11			*
12	gezelter	1390	* 2. Redistributions in binary form must reproduce the above copyright
13	xsun	955	* notice, this list of conditions and the following disclaimer in the
14			* documentation and/or other materials provided with the
15			* distribution.
16			*
17			* This software is provided "AS IS," without a warranty of any
18			* kind. All express or implied conditions, representations and
19			* warranties, including any implied warranty of merchantability,
20			* fitness for a particular purpose or non-infringement, are hereby
21			* excluded. The University of Notre Dame and its licensors shall not
22			* be liable for any damages suffered by licensee as a result of
23			* using, modifying or distributing the software or its
24			* derivatives. In no event will the University of Notre Dame or its
25			* licensors be liable for any lost revenue, profit or data, or for
26			* direct, indirect, special, consequential, incidental or punitive
27			* damages, however caused and regardless of the theory of liability,
28			* arising out of the use of or inability to use software, even if the
29			* University of Notre Dame has been advised of the possibility of
30			* such damages.
31			*
32	gezelter	1390	* SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your
33			* research, please cite the appropriate papers when you publish your
34			* work. Good starting points are:
35			*
36			* [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).
37			* [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).
38			* [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).
39	gezelter	1782	* [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010).
40			* [4] , Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). *
41	xsun	955	*
42			* Created by Xiuquan Sun on 05/09/06.
43			* @author Xiuquan Sun
44	gezelter	1442	* @version $Id$
45	xsun	955	*
46			*/
47
48			/* Calculates the undulation spectrum of the lipid membrance. */
49
50			#include <algorithm>
51			#include <fstream>
52			#include "applications/staticProps/Hxy.hpp"
53			#include "utils/simError.h"
54			#include "io/DumpReader.hpp"
55			#include "primitives/Molecule.hpp"
56			#include<stdio.h>
57			#include<string.h>
58			#include<stdlib.h>
59			#include<math.h>
60
61	gezelter	1390	namespace OpenMD {
62	xsun	955
63			Hxy::Hxy(SimInfo* info, const std::string& filename, const std::string& sele, int nbins_x, int nbins_y, int nrbins)
64			: StaticAnalyser(info, filename), selectionScript_(sele), evaluator_(info), seleMan_(info), nBinsX_(nbins_x), nBinsY_(nbins_y), nbins_(nrbins){
65
66			evaluator_.loadScriptString(sele);
67			if (!evaluator_.isDynamic()) {
68			seleMan_.setSelectionSet(evaluator_.evaluate());
69			}
70
71			gridsample_.resize(nBinsX_*nBinsY_);
72			gridZ_.resize(nBinsX_*nBinsY_);
73	xsun	967	mag.resize(nBinsX_*nBinsY_);
74			newmag.resize(nBinsX_*nBinsY_);
75	xsun	955
76	gezelter	956	sum_bin.resize(nbins_);
77			avg_bin.resize(nbins_);
78			errbin_sum.resize(nbins_);
79			errbin.resize(nbins_);
80			sum_bin_sq.resize(nbins_);
81			avg_bin_sq.resize(nbins_);
82			errbin_sum_sq.resize(nbins_);
83			errbin_sq.resize(nbins_);
84	xsun	955
85	xsun	967	bin.resize(nbins_);
86			samples.resize(nbins_);
87
88	xsun	955	setOutputName(getPrefix(filename) + ".Hxy");
89			}
90
91	xsun	967	Hxy::~Hxy(){
92			gridsample_.clear();
93			gridZ_.clear();
94			sum_bin.clear();
95			avg_bin.clear();
96			errbin_sum.clear();
97			errbin.clear();
98			sum_bin_sq.clear();
99			avg_bin_sq.clear();
100			errbin_sum_sq.clear();
101			errbin_sq.clear();
102
103	gezelter	1782	for(unsigned int i=0; i < bin.size(); i++)
104	xsun	967	bin[i].clear();
105	gezelter	1782	for(unsigned int i=0; i < samples.size(); i++)
106	xsun	967	samples[i].clear();
107
108			mag.clear();
109			newmag.clear();
110			}
111
112	xsun	955	void Hxy::process() {
113	gezelter	957	#if defined(HAVE_FFTW_H) \|\| defined(HAVE_DFFTW_H) \|\| defined(HAVE_FFTW3_H)
114	xsun	955	DumpReader reader(info_, dumpFilename_);
115			int nFrames = reader.getNFrames();
116			nProcessed_ = nFrames/step_;
117	gezelter	956
118	gezelter	1782	for(unsigned int k=0; k < bin.size(); k++)
119	xsun	967	bin[k].resize(nFrames);
120	gezelter	1782	for(unsigned int k=0; k < samples.size(); k++)
121	xsun	967	samples[k].resize(nFrames);
122
123	tim	963	RealType lenX_, lenY_;
124			RealType gridX_, gridY_;
125			RealType halfBoxX_, halfBoxY_;
126	xsun	967
127	tim	963	RealType interpsum, value;
128	xsun	955	int ninterp, px, py, newp;
129	gezelter	1796	int newindex, index;
130	xsun	955	int new_i, new_j, new_index;
131	xsun	967
132	tim	963	RealType freq_x, freq_y, zero_freq_x, zero_freq_y, freq;
133	xsun	967	RealType maxfreqx, maxfreqy, maxfreq;
134
135	xsun	955	int whichbin;
136	xsun	967
137	xsun	969	std::fill(sum_bin.begin(), sum_bin.end(), 0.0);
138			std::fill(avg_bin.begin(), avg_bin.end(), 0.0);
139			std::fill(errbin_sum.begin(), errbin_sum.end(), 0.0);
140			std::fill(errbin.begin(), errbin.end(), 0.0);
141			std::fill(sum_bin_sq.begin(), sum_bin_sq.end(), 0.0);
142			std::fill(avg_bin_sq.begin(), avg_bin_sq.end(), 0.0);
143			std::fill(errbin_sum_sq.begin(), errbin_sum_sq.end(), 0.0);
144			std::fill(errbin_sq.begin(), errbin_sq.end(), 0.0);
145
146	gezelter	1782	for(unsigned int i=0; i < bin.size(); i++)
147	xsun	969	std::fill(bin[i].begin(), bin[i].end(), 0.0);
148
149	gezelter	1782	for(unsigned int i=0; i < samples.size(); i++)
150	xsun	969	std::fill(samples[i].begin(), samples[i].end(), 0);
151
152	xsun	955	for (int istep = 0; istep < nFrames; istep += step_) {
153
154	gezelter	956	reader.readFrame(istep);
155			currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
156	gezelter	957
157	gezelter	956	Mat3x3d hmat = currentSnapshot_->getHmat();
158
159	gezelter	957	#ifdef HAVE_FFTW3_H
160			fftw_plan p;
161			#else
162			fftwnd_plan p;
163			#endif
164	xsun	955	fftw_complex in, out;
165	gezelter	956
166			in = (fftw_complex) fftw_malloc(sizeof(fftw_complex) (nBinsX_*nBinsY_));
167			out = (fftw_complex) fftw_malloc(sizeof(fftw_complex) (nBinsX_*nBinsY_));
168	gezelter	957
169			#ifdef HAVE_FFTW3_H
170			p = fftw_plan_dft_2d(nBinsX_, nBinsY_, in, out,
171			FFTW_FORWARD, FFTW_ESTIMATE);
172			#else
173			p = fftw2d_create_plan(nBinsX_, nBinsY_, FFTW_FORWARD, FFTW_ESTIMATE);
174			#endif
175	xsun	969
176	xsun	967	std::fill(gridsample_.begin(), gridsample_.end(), 0);
177			std::fill(gridZ_.begin(), gridZ_.end(), 0.0);
178			std::fill(mag.begin(), mag.end(), 0.0);
179			std::fill(newmag.begin(), newmag.end(), 0.0);
180
181	xsun	969	int i, j;
182	gezelter	956
183	xsun	955	StuntDouble* sd;
184	gezelter	956
185	xsun	955	lenX_ = hmat(0,0);
186			lenY_ = hmat(1,1);
187	gezelter	956
188	xsun	955	gridX_ = lenX_ /(nBinsX_);
189			gridY_ = lenY_ /(nBinsY_);
190	gezelter	956
191	xsun	967	halfBoxX_ = lenX_ / 2.0;
192			halfBoxY_ = lenY_ / 2.0;
193	gezelter	956
194	xsun	955	if (evaluator_.isDynamic()) {
195			seleMan_.setSelectionSet(evaluator_.evaluate());
196			}
197
198			//wrap the stuntdoubles into a cell
199			for (sd = seleMan_.beginSelected(i); sd != NULL; sd = seleMan_.nextSelected(i)) {
200			Vector3d pos = sd->getPos();
201	gezelter	1078	if (usePeriodicBoundaryConditions_)
202			currentSnapshot_->wrapVector(pos);
203	xsun	955	sd->setPos(pos);
204	xsun	967	}
205	xsun	955
206			//determine which atom belongs to which grid
207			for (sd = seleMan_.beginSelected(i); sd != NULL; sd = seleMan_.nextSelected(i)) {
208			Vector3d pos = sd->getPos();
209			//int binNo = (pos.z() /deltaR_) - 1;
210	xsun	967	int binNoX = (int) ((pos.x() + halfBoxX_) / gridX_);
211			int binNoY = (int) ((pos.y() + halfBoxY_) / gridY_);
212	xsun	955	//std::cout << "pos.z = " << pos.z() << " halfBoxZ_ = " << halfBoxZ_ << " deltaR_ = " << deltaR_ << " binNo = " << binNo << "\n";
213			gridZ_[binNoX*nBinsY_+binNoY] += pos.z();
214			gridsample_[binNoX*nBinsY_+binNoY]++;
215			}
216
217			// FFT stuff depends on nx and ny, so delay allocation until we have
218			// that information
219
220			for(i = 0; i < nBinsX_; i++){
221			for(j = 0; j < nBinsY_; j++){
222			newindex = i * nBinsY_ + j;
223			if(gridsample_[newindex] > 0){
224	tim	963	gridZ_[newindex] = gridZ_[newindex] / (RealType)gridsample_[newindex];
225	xsun	955	}
226			}
227			}
228	gezelter	956
229	xsun	955	for (i=0; i< nBinsX_; i++) {
230			for(j=0; j< nBinsY_; j++) {
231			newindex = i*nBinsY_ + j;
232			if (gridsample_[newindex] == 0) {
233			// interpolate from surrounding points:
234
235			interpsum = 0.0;
236			ninterp = 0;
237
238			//point1 = bottom;
239
240			px = i;
241			py = j - 1;
242			newp = px*nBinsY_ + py;
243			if ((py >= 0) && (gridsample_[newp] > 0)) {
244			interpsum += gridZ_[newp];
245			ninterp++;
246			}
247
248			//point2 = top;
249
250			px = i;
251			py = j + 1;
252			newp = px*nBinsY_ + py;
253			if ((py < nBinsY_) && (gridsample_[newp] > 0)) {
254			interpsum += gridZ_[newp];
255			ninterp++;
256			}
257
258			//point3 = left;
259
260			px = i - 1;
261			py = j;
262			newp = px*nBinsY_ + py;
263			if ((px >= 0) && (gridsample_[newp] > 0)) {
264			interpsum += gridZ_[newp];
265			ninterp++;
266			}
267
268			//point4 = right;
269
270			px = i + 1;
271			py = j;
272			newp = px*nBinsY_ + py;
273	gezelter	956	if ( (px < nBinsX_ ) && ( gridsample_[newp] > 0 )) {
274			interpsum += gridZ_[newp];
275			ninterp++;
276			}
277
278	tim	963	value = interpsum / (RealType)ninterp;
279	xsun	955
280			gridZ_[newindex] = value;
281			}
282			}
283			}
284
285			for (i=0; i < nBinsX_; i++) {
286			for (j=0; j < nBinsY_; j++) {
287			newindex = i*nBinsY_ + j;
288
289			c_re(in[newindex]) = gridZ_[newindex];
290			c_im(in[newindex]) = 0.0;
291			}
292			}
293	xsun	967
294	gezelter	957	#ifdef HAVE_FFTW3_H
295			fftw_execute(p);
296			#else
297			fftwnd_one(p, in, out);
298			#endif
299	xsun	955
300			for (i=0; i< nBinsX_; i++) {
301			for(j=0; j< nBinsY_; j++) {
302			newindex = i*nBinsY_ + j;
303			mag[newindex] = pow(c_re(out[newindex]),2) + pow(c_im(out[newindex]),2);
304			}
305			}
306	xsun	967
307	gezelter	957	#ifdef HAVE_FFTW3_H
308	xsun	955	fftw_destroy_plan(p);
309	gezelter	957	#else
310			fftwnd_destroy_plan(p);
311			#endif
312	xsun	955	fftw_free(out);
313			fftw_free(in);
314
315			for (i=0; i< (nBinsX_/2); i++) {
316			for(j=0; j< (nBinsY_/2); 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			for (i=(nBinsX_/2); i< nBinsX_; i++) {
326			for(j=0; j< (nBinsY_/2); j++) {
327			index = i*nBinsY_ + j;
328			new_i = i - (nBinsX_/2);
329			new_j = j + (nBinsY_/2);
330			new_index = new_i*nBinsY_ + new_j;
331			newmag[new_index] = mag[index];
332			}
333			}
334
335			for (i=0; i< (nBinsX_/2); i++) {
336			for(j=(nBinsY_/2); j< nBinsY_; j++) {
337			index = i*nBinsY_ + j;
338			new_i = i + (nBinsX_/2);
339			new_j = j - (nBinsY_/2);
340			new_index = new_i*nBinsY_ + new_j;
341			newmag[new_index] = mag[index];
342			}
343			}
344
345			for (i=(nBinsX_/2); i< nBinsX_; i++) {
346			for(j=(nBinsY_/2); j< nBinsY_; j++) {
347			index = i*nBinsY_ + j;
348			new_i = i - (nBinsX_/2);
349			new_j = j - (nBinsY_/2);
350			new_index = new_i*nBinsY_ + new_j;
351			newmag[new_index] = mag[index];
352			}
353			}
354
355			maxfreqx = 1.0 / gridX_;
356			maxfreqy = 1.0 / gridY_;
357
358			// printf("%lf\t%lf\t%lf\t%lf\n", dx, dy, maxfreqx, maxfreqy);
359
360			maxfreq = sqrt(maxfreqxmaxfreqx + maxfreqymaxfreqy);
361	tim	963	dfreq = maxfreq/(RealType)(nbins_-1);
362	xsun	955
363			//printf("%lf\n", dfreq);
364
365			zero_freq_x = nBinsX_/2;
366			zero_freq_y = nBinsY_/2;
367
368			for (i=0; i< nBinsX_; i++) {
369			for(j=0; j< nBinsY_; j++) {
370
371	tim	963	freq_x = (RealType)(i - zero_freq_x)maxfreqx2 / nBinsX_;
372			freq_y = (RealType)(j - zero_freq_y)maxfreqy2 / nBinsY_;
373	xsun	955
374			freq = sqrt(freq_xfreq_x + freq_yfreq_y);
375
376			whichbin = (int) (freq / dfreq);
377			newindex = i*nBinsY_ + j;
378			// printf("%d %d %lf %lf\n", whichbin, newindex, freq, dfreq);
379			bin[whichbin][istep] += newmag[newindex];
380			samples[whichbin][istep]++;
381			}
382			}
383
384	gezelter	956	for ( i = 0; i < nbins_; i++) {
385	xsun	955	if ( samples[i][istep] > 0) {
386	xsun	969	bin[i][istep] = 4.0 * sqrt(bin[i][istep] / (RealType)samples[i][istep]) / (RealType)nBinsX_ / (RealType)nBinsY_;
387	xsun	955	}
388			}
389			}
390	xsun	967
391	gezelter	956	for (int i = 0; i < nbins_; i++) {
392			for (int j = 0; j < nFrames; j++) {
393	xsun	955	sum_bin[i] += bin[i][j];
394			sum_bin_sq[i] += bin[i][j] * bin[i][j];
395			}
396	tim	963	avg_bin[i] = sum_bin[i] / (RealType)nFrames;
397			avg_bin_sq[i] = sum_bin_sq[i] / (RealType)nFrames;
398	gezelter	956	for (int j = 0; j < nFrames; j++) {
399	xsun	955	errbin_sum[i] += pow((bin[i][j] - avg_bin[i]), 2);
400			errbin_sum_sq[i] += pow((bin[i][j] * bin[i][j] - avg_bin_sq[i]), 2);
401			}
402	tim	963	errbin[i] = sqrt( errbin_sum[i] / (RealType)nFrames );
403			errbin_sq[i] = sqrt( errbin_sum_sq[i] / (RealType)nFrames );
404	xsun	955	}
405	xsun	967
406	gezelter	956	printSpectrum();
407	xsun	967
408	gezelter	956	#else
409			sprintf(painCave.errMsg, "Hxy: FFTW support was not compiled in!\n");
410			painCave.isFatal = 1;
411			simError();
412	xsun	955
413	gezelter	956	#endif
414	xsun	967	}
415	xsun	955
416			void Hxy::printSpectrum() {
417			std::ofstream rdfStream(outputFilename_.c_str());
418			if (rdfStream.is_open()) {
419	xsun	967
420			for (int i = 0; i < nbins_; ++i) {
421	xsun	955	if ( avg_bin[i] > 0 ){
422	xsun	967	rdfStream << (RealType)i * dfreq << "\t"
423	xsun	955	<<pow(avg_bin[i], 2)<<"\t"
424			<<errbin_sq[i]<<"\t"
425			<<avg_bin[i]<<"\t"
426			<<errbin[i]<<"\n";
427			}
428			}
429			} else {
430
431			sprintf(painCave.errMsg, "Hxy: unable to open %s\n", outputFilename_.c_str());
432			painCave.isFatal = 1;
433			simError();
434			}
435	xsun	967
436	xsun	955	rdfStream.close();
437	xsun	967
438	xsun	955	}
439
440			}