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.5 2006-05-22 15:30:42 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>

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_);
    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(int i=0; i < bin.size(); i++)
        bin[i].clear();
      for(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(int k=0; k < bin.size(); k++)
      bin[k].resize(nFrames);
    for(int k=0; k < samples.size(); k++)
      samples[k].resize(nFrames);

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

    int binNoX, binNoY;
    RealType interpsum, value;
    int ninterp, px, py, newp;
    int newx, newy, 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;
    int nMolecules;

    for (int istep = 0; istep < nFrames; istep += step_) {
      
      reader.readFrame(istep);
      currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
      nMolecules = info_->getNGlobalMolecules();
      
      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
      
      int i, j;   
      
      std::fill(gridsample_.begin(), gridsample_.end(), 0);
      std::fill(gridZ_.begin(), gridZ_.end(), 0.0);
      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);
      std::fill(mag.begin(), mag.end(), 0.0);
      std::fill(newmag.begin(), newmag.end(), 0.0);

      for(i=0; i < bin.size(); i++)
        std::fill(bin[i].begin(), bin[i].end(), 0.0);

      for(i=0; i < samples.size(); i++)
        std::fill(samples[i].begin(), samples[i].end(), 0);
      
      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();
        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)nMolecules;
        }
      }    
    }

    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:	2763
Committed:	Mon May 22 15:30:42 2006 UTC (18 years, 1 month ago) by xsun
File size:	12931 byte(s)
Log Message:	fixed the bugs of Hxy.cpp.
#	Content
1	/*
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.5 2006-05-22 15:30:42 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
64	namespace oopse {
65
66	Hxy::Hxy(SimInfo* info, const std::string& filename, const std::string& sele, int nbins_x, int nbins_y, int nrbins)
67	: StaticAnalyser(info, filename), selectionScript_(sele), evaluator_(info), seleMan_(info), nBinsX_(nbins_x), nBinsY_(nbins_y), nbins_(nrbins){
68
69	evaluator_.loadScriptString(sele);
70	if (!evaluator_.isDynamic()) {
71	seleMan_.setSelectionSet(evaluator_.evaluate());
72	}
73
74	gridsample_.resize(nBinsX_*nBinsY_);
75	gridZ_.resize(nBinsX_*nBinsY_);
76	mag.resize(nBinsX_*nBinsY_);
77	newmag.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	bin.resize(nbins_);
89	samples.resize(nbins_);
90
91	setOutputName(getPrefix(filename) + ".Hxy");
92	}
93
94	Hxy::~Hxy(){
95	gridsample_.clear();
96	gridZ_.clear();
97	sum_bin.clear();
98	avg_bin.clear();
99	errbin_sum.clear();
100	errbin.clear();
101	sum_bin_sq.clear();
102	avg_bin_sq.clear();
103	errbin_sum_sq.clear();
104	errbin_sq.clear();
105
106	for(int i=0; i < bin.size(); i++)
107	bin[i].clear();
108	for(int i=0; i < samples.size(); i++)
109	samples[i].clear();
110
111	mag.clear();
112	newmag.clear();
113	}
114
115	void Hxy::process() {
116	#if defined(HAVE_FFTW_H) \|\| defined(HAVE_DFFTW_H) \|\| defined(HAVE_FFTW3_H)
117	DumpReader reader(info_, dumpFilename_);
118	int nFrames = reader.getNFrames();
119	nProcessed_ = nFrames/step_;
120
121	for(int k=0; k < bin.size(); k++)
122	bin[k].resize(nFrames);
123	for(int k=0; k < samples.size(); k++)
124	samples[k].resize(nFrames);
125
126	RealType lenX_, lenY_;
127	RealType gridX_, gridY_;
128	RealType halfBoxX_, halfBoxY_;
129
130	int binNoX, binNoY;
131	RealType interpsum, value;
132	int ninterp, px, py, newp;
133	int newx, newy, newindex, index;
134	int new_i, new_j, new_index;
135
136	RealType freq_x, freq_y, zero_freq_x, zero_freq_y, freq;
137	RealType maxfreqx, maxfreqy, maxfreq;
138
139	int whichbin;
140	int nMolecules;
141
142	for (int istep = 0; istep < nFrames; istep += step_) {
143
144	reader.readFrame(istep);
145	currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
146	nMolecules = info_->getNGlobalMolecules();
147
148	Mat3x3d hmat = currentSnapshot_->getHmat();
149
150	#ifdef HAVE_FFTW3_H
151	fftw_plan p;
152	#else
153	fftwnd_plan p;
154	#endif
155	fftw_complex in, out;
156
157	in = (fftw_complex) fftw_malloc(sizeof(fftw_complex) (nBinsX_*nBinsY_));
158	out = (fftw_complex) fftw_malloc(sizeof(fftw_complex) (nBinsX_*nBinsY_));
159
160	#ifdef HAVE_FFTW3_H
161	p = fftw_plan_dft_2d(nBinsX_, nBinsY_, in, out,
162	FFTW_FORWARD, FFTW_ESTIMATE);
163	#else
164	p = fftw2d_create_plan(nBinsX_, nBinsY_, FFTW_FORWARD, FFTW_ESTIMATE);
165	#endif
166
167	int i, j;
168
169	std::fill(gridsample_.begin(), gridsample_.end(), 0);
170	std::fill(gridZ_.begin(), gridZ_.end(), 0.0);
171	std::fill(sum_bin.begin(), sum_bin.end(), 0.0);
172	std::fill(avg_bin.begin(), avg_bin.end(), 0.0);
173	std::fill(errbin_sum.begin(), errbin_sum.end(), 0.0);
174	std::fill(errbin.begin(), errbin.end(), 0.0);
175	std::fill(sum_bin_sq.begin(), sum_bin_sq.end(), 0.0);
176	std::fill(avg_bin_sq.begin(), avg_bin_sq.end(), 0.0);
177	std::fill(errbin_sum_sq.begin(), errbin_sum_sq.end(), 0.0);
178	std::fill(errbin_sq.begin(), errbin_sq.end(), 0.0);
179	std::fill(mag.begin(), mag.end(), 0.0);
180	std::fill(newmag.begin(), newmag.end(), 0.0);
181
182	for(i=0; i < bin.size(); i++)
183	std::fill(bin[i].begin(), bin[i].end(), 0.0);
184
185	for(i=0; i < samples.size(); i++)
186	std::fill(samples[i].begin(), samples[i].end(), 0);
187
188	StuntDouble* sd;
189
190	lenX_ = hmat(0,0);
191	lenY_ = hmat(1,1);
192
193	gridX_ = lenX_ /(nBinsX_);
194	gridY_ = lenY_ /(nBinsY_);
195
196	halfBoxX_ = lenX_ / 2.0;
197	halfBoxY_ = lenY_ / 2.0;
198
199	if (evaluator_.isDynamic()) {
200	seleMan_.setSelectionSet(evaluator_.evaluate());
201	}
202
203	//wrap the stuntdoubles into a cell
204	for (sd = seleMan_.beginSelected(i); sd != NULL; sd = seleMan_.nextSelected(i)) {
205	Vector3d pos = sd->getPos();
206	currentSnapshot_->wrapVector(pos);
207	sd->setPos(pos);
208	}
209
210	//determine which atom belongs to which grid
211	for (sd = seleMan_.beginSelected(i); sd != NULL; sd = seleMan_.nextSelected(i)) {
212	Vector3d pos = sd->getPos();
213	//int binNo = (pos.z() /deltaR_) - 1;
214	int binNoX = (int) ((pos.x() + halfBoxX_) / gridX_);
215	int binNoY = (int) ((pos.y() + halfBoxY_) / gridY_);
216	//std::cout << "pos.z = " << pos.z() << " halfBoxZ_ = " << halfBoxZ_ << " deltaR_ = " << deltaR_ << " binNo = " << binNo << "\n";
217	gridZ_[binNoX*nBinsY_+binNoY] += pos.z();
218	gridsample_[binNoX*nBinsY_+binNoY]++;
219	}
220
221	// FFT stuff depends on nx and ny, so delay allocation until we have
222	// that information
223
224	for(i = 0; i < nBinsX_; i++){
225	for(j = 0; j < nBinsY_; j++){
226	newindex = i * nBinsY_ + j;
227	if(gridsample_[newindex] > 0){
228	gridZ_[newindex] = gridZ_[newindex] / (RealType)gridsample_[newindex];
229	}
230	}
231	}
232
233	for (i=0; i< nBinsX_; i++) {
234	for(j=0; j< nBinsY_; j++) {
235	newindex = i*nBinsY_ + j;
236	if (gridsample_[newindex] == 0) {
237	// interpolate from surrounding points:
238
239	interpsum = 0.0;
240	ninterp = 0;
241
242	//point1 = bottom;
243
244	px = i;
245	py = j - 1;
246	newp = px*nBinsY_ + py;
247	if ((py >= 0) && (gridsample_[newp] > 0)) {
248	interpsum += gridZ_[newp];
249	ninterp++;
250	}
251
252	//point2 = top;
253
254	px = i;
255	py = j + 1;
256	newp = px*nBinsY_ + py;
257	if ((py < nBinsY_) && (gridsample_[newp] > 0)) {
258	interpsum += gridZ_[newp];
259	ninterp++;
260	}
261
262	//point3 = left;
263
264	px = i - 1;
265	py = j;
266	newp = px*nBinsY_ + py;
267	if ((px >= 0) && (gridsample_[newp] > 0)) {
268	interpsum += gridZ_[newp];
269	ninterp++;
270	}
271
272	//point4 = right;
273
274	px = i + 1;
275	py = j;
276	newp = px*nBinsY_ + py;
277	if ( (px < nBinsX_ ) && ( gridsample_[newp] > 0 )) {
278	interpsum += gridZ_[newp];
279	ninterp++;
280	}
281
282	value = interpsum / (RealType)ninterp;
283
284	gridZ_[newindex] = value;
285	}
286	}
287	}
288
289	for (i=0; i < nBinsX_; i++) {
290	for (j=0; j < nBinsY_; j++) {
291	newindex = i*nBinsY_ + j;
292
293	c_re(in[newindex]) = gridZ_[newindex];
294	c_im(in[newindex]) = 0.0;
295	}
296	}
297
298	#ifdef HAVE_FFTW3_H
299	fftw_execute(p);
300	#else
301	fftwnd_one(p, in, out);
302	#endif
303
304	for (i=0; i< nBinsX_; i++) {
305	for(j=0; j< nBinsY_; j++) {
306	newindex = i*nBinsY_ + j;
307	mag[newindex] = pow(c_re(out[newindex]),2) + pow(c_im(out[newindex]),2);
308	}
309	}
310
311	#ifdef HAVE_FFTW3_H
312	fftw_destroy_plan(p);
313	#else
314	fftwnd_destroy_plan(p);
315	#endif
316	fftw_free(out);
317	fftw_free(in);
318
319	for (i=0; i< (nBinsX_/2); i++) {
320	for(j=0; j< (nBinsY_/2); j++) {
321	index = i*nBinsY_ + j;
322	new_i = i + (nBinsX_/2);
323	new_j = j + (nBinsY_/2);
324	new_index = new_i*nBinsY_ + new_j;
325	newmag[new_index] = mag[index];
326	}
327	}
328
329	for (i=(nBinsX_/2); i< nBinsX_; i++) {
330	for(j=0; j< (nBinsY_/2); j++) {
331	index = i*nBinsY_ + j;
332	new_i = i - (nBinsX_/2);
333	new_j = j + (nBinsY_/2);
334	new_index = new_i*nBinsY_ + new_j;
335	newmag[new_index] = mag[index];
336	}
337	}
338
339	for (i=0; i< (nBinsX_/2); i++) {
340	for(j=(nBinsY_/2); j< nBinsY_; j++) {
341	index = i*nBinsY_ + j;
342	new_i = i + (nBinsX_/2);
343	new_j = j - (nBinsY_/2);
344	new_index = new_i*nBinsY_ + new_j;
345	newmag[new_index] = mag[index];
346	}
347	}
348
349	for (i=(nBinsX_/2); i< nBinsX_; i++) {
350	for(j=(nBinsY_/2); j< nBinsY_; j++) {
351	index = i*nBinsY_ + j;
352	new_i = i - (nBinsX_/2);
353	new_j = j - (nBinsY_/2);
354	new_index = new_i*nBinsY_ + new_j;
355	newmag[new_index] = mag[index];
356	}
357	}
358
359	maxfreqx = 1.0 / gridX_;
360	maxfreqy = 1.0 / gridY_;
361
362	// printf("%lf\t%lf\t%lf\t%lf\n", dx, dy, maxfreqx, maxfreqy);
363
364	maxfreq = sqrt(maxfreqxmaxfreqx + maxfreqymaxfreqy);
365	dfreq = maxfreq/(RealType)(nbins_-1);
366
367	//printf("%lf\n", dfreq);
368
369	zero_freq_x = nBinsX_/2;
370	zero_freq_y = nBinsY_/2;
371
372	for (i=0; i< nBinsX_; i++) {
373	for(j=0; j< nBinsY_; j++) {
374
375	freq_x = (RealType)(i - zero_freq_x)maxfreqx2 / nBinsX_;
376	freq_y = (RealType)(j - zero_freq_y)maxfreqy2 / nBinsY_;
377
378	freq = sqrt(freq_xfreq_x + freq_yfreq_y);
379
380	whichbin = (int) (freq / dfreq);
381	newindex = i*nBinsY_ + j;
382	// printf("%d %d %lf %lf\n", whichbin, newindex, freq, dfreq);
383	bin[whichbin][istep] += newmag[newindex];
384	samples[whichbin][istep]++;
385	}
386	}
387
388	for ( i = 0; i < nbins_; i++) {
389	if ( samples[i][istep] > 0) {
390	bin[i][istep] = 4.0 * sqrt(bin[i][istep] / (RealType)samples[i][istep]) / (RealType)nMolecules;
391	}
392	}
393	}
394
395	for (int i = 0; i < nbins_; i++) {
396	for (int j = 0; j < nFrames; j++) {
397	sum_bin[i] += bin[i][j];
398	sum_bin_sq[i] += bin[i][j] * bin[i][j];
399	}
400	avg_bin[i] = sum_bin[i] / (RealType)nFrames;
401	avg_bin_sq[i] = sum_bin_sq[i] / (RealType)nFrames;
402	for (int j = 0; j < nFrames; j++) {
403	errbin_sum[i] += pow((bin[i][j] - avg_bin[i]), 2);
404	errbin_sum_sq[i] += pow((bin[i][j] * bin[i][j] - avg_bin_sq[i]), 2);
405	}
406	errbin[i] = sqrt( errbin_sum[i] / (RealType)nFrames );
407	errbin_sq[i] = sqrt( errbin_sum_sq[i] / (RealType)nFrames );
408	}
409
410	printSpectrum();
411
412	#else
413	sprintf(painCave.errMsg, "Hxy: FFTW support was not compiled in!\n");
414	painCave.isFatal = 1;
415	simError();
416
417	#endif
418	}
419
420	void Hxy::printSpectrum() {
421	std::ofstream rdfStream(outputFilename_.c_str());
422	if (rdfStream.is_open()) {
423
424	for (int i = 0; i < nbins_; ++i) {
425	if ( avg_bin[i] > 0 ){
426	rdfStream << (RealType)i * dfreq << "\t"
427	<<pow(avg_bin[i], 2)<<"\t"
428	<<errbin_sq[i]<<"\t"
429	<<avg_bin[i]<<"\t"
430	<<errbin[i]<<"\n";
431	}
432	}
433	} else {
434
435	sprintf(painCave.errMsg, "Hxy: unable to open %s\n", outputFilename_.c_str());
436	painCave.isFatal = 1;
437	simError();
438	}
439
440	rdfStream.close();
441
442	}
443
444	}