applications/dynamicProps/LegendreCorrFunc.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.
 */

#include "applications/dynamicProps/LegendreCorrFunc.hpp"
#include "math/LegendrePolynomial.hpp"
#include "utils/simError.h"

namespace oopse {
  LegendreCorrFunc::LegendreCorrFunc(SimInfo* info, const std::string& filename, const std::string& sele1, const std::string& sele2, int order)
    : ParticleTimeCorrFunc(info, filename, sele1, sele2, DataStorage::dslAmat | DataStorage::dslElectroFrame){

      setCorrFuncType("Legendre Correlation Function");
      setOutputName(getPrefix(dumpFilename_) + ".lcorr");
      histogram_.resize(nTimeBins_); 
      count_.resize(nTimeBins_);
      nSelected_ =   seleMan1_.getSelectionCount();  
      std::cout << "sele1 = " << sele1 << "\n";
      std::cout << "sele2 = " << sele2 << "\n";
      std::cout << "nSelected = " << nSelected_ << "\n";
      assert(  nSelected_ == seleMan2_.getSelectionCount());

      LegendrePolynomial polynomial(order);
      legendre_ = polynomial.getLegendrePolynomial(order);
    }

  void LegendreCorrFunc::correlateFrames(int frame1, int frame2) {
    Snapshot* snapshot1 = bsMan_->getSnapshot(frame1);
    Snapshot* snapshot2 = bsMan_->getSnapshot(frame2);
    assert(snapshot1 && snapshot2);

    RealType time1 = snapshot1->getTime();
    RealType time2 = snapshot2->getTime();

    int timeBin = int ((time2 - time1) /deltaTime_ + 0.5);
    count_[timeBin] += nSelected_;    

    int i;
    int j;
    StuntDouble* sd1;
    StuntDouble* sd2;

    for (sd1 = seleMan1_.beginSelected(i), sd2 = seleMan2_.beginSelected(j);
         sd1 != NULL && sd2 != NULL;
         sd1 = seleMan1_.nextSelected(i), sd2 = seleMan2_.nextSelected(j)) {

      Vector3d corrVals = calcCorrVals(frame1, frame2, sd1, sd2);
      histogram_[timeBin] += corrVals; 
    }
    
  }

  void LegendreCorrFunc::postCorrelate() {
    for (int i =0 ; i < nTimeBins_; ++i) {
      if (count_[i] > 0) {
        histogram_[i] /= count_[i];
      }
    }
  }

  void LegendreCorrFunc::preCorrelate() {
    // Fill the histogram with empty Vector3d:
    std::fill(histogram_.begin(), histogram_.end(), Vector3d(0.0));
    // count array set to zero
    std::fill(count_.begin(), count_.end(), 0);
  }


  Vector3d LegendreCorrFunc::calcCorrVals(int frame1, int frame2, StuntDouble* sd1,  StuntDouble* sd2) {
    Vector3d v1x = sd1->getA(frame1).getColumn(0);
    Vector3d v2x = sd2->getA(frame2).getColumn(0);
    Vector3d v1y = sd1->getA(frame1).getColumn(1);
    Vector3d v2y = sd2->getA(frame2).getColumn(1);
    Vector3d v1z = sd1->getA(frame1).getColumn(2);
    Vector3d v2z = sd2->getA(frame2).getColumn(2);

    RealType uxprod = legendre_.evaluate(dot(v1x, v2x)/(v1x.length()*v2x.length()));
    RealType uyprod = legendre_.evaluate(dot(v1y, v2y)/(v1y.length()*v2y.length()));
    RealType uzprod = legendre_.evaluate(dot(v1z, v2z)/(v1z.length()*v2z.length()));

    return Vector3d(uxprod, uyprod, uzprod);

  }


  void LegendreCorrFunc::validateSelection(const SelectionManager& seleMan) {
    StuntDouble* sd;
    int i;    
    for (sd = seleMan1_.beginSelected(i); sd != NULL; sd = seleMan1_.nextSelected(i)) {
      if (!sd->isDirectionalAtom()) {
        sprintf(painCave.errMsg,
                "LegendreCorrFunc::validateSelection Error: selected atoms are not Directional\n");
        painCave.isFatal = 1;
        simError();        
      }
    }
    
  }

  void LegendreCorrFunc::writeCorrelate() {
    std::ofstream ofs(getOutputFileName().c_str());

    if (ofs.is_open()) {

      ofs << "#" << getCorrFuncType() << "\n";
      ofs << "#time\tPn(costheta_x)\tPn(costheta_y)\tPn(costheta_z)\n";

      for (int i = 0; i < nTimeBins_; ++i) {
        ofs << time_[i] << "\t" << 
          histogram_[i](0) << "\t" <<
          histogram_[i](1) << "\t" <<
          histogram_[i](2) << "\t" << "\n";
      }
            
    } else {
      sprintf(painCave.errMsg,
              "LegendreCorrFunc::writeCorrelate Error: fail to open %s\n", getOutputFileName().c_str());
      painCave.isFatal = 1;
      simError();        
    }
    ofs.close();    
  }
}
Revision:	3235
Committed:	Tue Oct 2 12:05:25 2007 UTC (16 years, 11 months ago) by xsun
File size:	6086 byte(s)
Log Message:	Added: Action Correlation function (for viscosity), Directional RCorr (for non-spherical bodies) Modified: LCorr (for multiple axes)
#	User	Rev	Content
1	tim	2576	/*
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			#include "applications/dynamicProps/LegendreCorrFunc.hpp"
43			#include "math/LegendrePolynomial.hpp"
44			#include "utils/simError.h"
45
46			namespace oopse {
47			LegendreCorrFunc::LegendreCorrFunc(SimInfo* info, const std::string& filename, const std::string& sele1, const std::string& sele2, int order)
48	xsun	3235	: ParticleTimeCorrFunc(info, filename, sele1, sele2, DataStorage::dslAmat \| DataStorage::dslElectroFrame){
49	tim	2576
50			setCorrFuncType("Legendre Correlation Function");
51			setOutputName(getPrefix(dumpFilename_) + ".lcorr");
52	xsun	3235	histogram_.resize(nTimeBins_);
53			count_.resize(nTimeBins_);
54			nSelected_ = seleMan1_.getSelectionCount();
55			std::cout << "sele1 = " << sele1 << "\n";
56			std::cout << "sele2 = " << sele2 << "\n";
57			std::cout << "nSelected = " << nSelected_ << "\n";
58			assert( nSelected_ == seleMan2_.getSelectionCount());
59
60	tim	2576	LegendrePolynomial polynomial(order);
61			legendre_ = polynomial.getLegendrePolynomial(order);
62	xsun	3235	}
63	tim	2576
64	xsun	3235	void LegendreCorrFunc::correlateFrames(int frame1, int frame2) {
65			Snapshot* snapshot1 = bsMan_->getSnapshot(frame1);
66			Snapshot* snapshot2 = bsMan_->getSnapshot(frame2);
67			assert(snapshot1 && snapshot2);
68
69			RealType time1 = snapshot1->getTime();
70			RealType time2 = snapshot2->getTime();
71
72			int timeBin = int ((time2 - time1) /deltaTime_ + 0.5);
73			count_[timeBin] += nSelected_;
74
75			int i;
76			int j;
77			StuntDouble* sd1;
78			StuntDouble* sd2;
79
80			for (sd1 = seleMan1_.beginSelected(i), sd2 = seleMan2_.beginSelected(j);
81			sd1 != NULL && sd2 != NULL;
82			sd1 = seleMan1_.nextSelected(i), sd2 = seleMan2_.nextSelected(j)) {
83
84			Vector3d corrVals = calcCorrVals(frame1, frame2, sd1, sd2);
85			histogram_[timeBin] += corrVals;
86	tim	2576	}
87	xsun	3235
88			}
89	tim	2576
90	xsun	3235	void LegendreCorrFunc::postCorrelate() {
91			for (int i =0 ; i < nTimeBins_; ++i) {
92			if (count_[i] > 0) {
93			histogram_[i] /= count_[i];
94			}
95			}
96			}
97	tim	2576
98	xsun	3235	void LegendreCorrFunc::preCorrelate() {
99			// Fill the histogram with empty Vector3d:
100			std::fill(histogram_.begin(), histogram_.end(), Vector3d(0.0));
101			// count array set to zero
102			std::fill(count_.begin(), count_.end(), 0);
103	tim	2576	}
104
105
106	xsun	3235
107			Vector3d LegendreCorrFunc::calcCorrVals(int frame1, int frame2, StuntDouble* sd1, StuntDouble* sd2) {
108			Vector3d v1x = sd1->getA(frame1).getColumn(0);
109			Vector3d v2x = sd2->getA(frame2).getColumn(0);
110			Vector3d v1y = sd1->getA(frame1).getColumn(1);
111			Vector3d v2y = sd2->getA(frame2).getColumn(1);
112			Vector3d v1z = sd1->getA(frame1).getColumn(2);
113			Vector3d v2z = sd2->getA(frame2).getColumn(2);
114
115			RealType uxprod = legendre_.evaluate(dot(v1x, v2x)/(v1x.length()*v2x.length()));
116			RealType uyprod = legendre_.evaluate(dot(v1y, v2y)/(v1y.length()*v2y.length()));
117			RealType uzprod = legendre_.evaluate(dot(v1z, v2z)/(v1z.length()*v2z.length()));
118
119			return Vector3d(uxprod, uyprod, uzprod);
120
121			}
122
123
124	tim	2576	void LegendreCorrFunc::validateSelection(const SelectionManager& seleMan) {
125			StuntDouble* sd;
126			int i;
127			for (sd = seleMan1_.beginSelected(i); sd != NULL; sd = seleMan1_.nextSelected(i)) {
128			if (!sd->isDirectionalAtom()) {
129			sprintf(painCave.errMsg,
130	xsun	3235	"LegendreCorrFunc::validateSelection Error: selected atoms are not Directional\n");
131	tim	2576	painCave.isFatal = 1;
132			simError();
133			}
134			}
135
136			}
137
138	xsun	3235	void LegendreCorrFunc::writeCorrelate() {
139			std::ofstream ofs(getOutputFileName().c_str());
140	tim	2576
141	xsun	3235	if (ofs.is_open()) {
142	tim	2576
143	xsun	3235	ofs << "#" << getCorrFuncType() << "\n";
144			ofs << "#time\tPn(costheta_x)\tPn(costheta_y)\tPn(costheta_z)\n";
145	tim	2576
146	xsun	3235	for (int i = 0; i < nTimeBins_; ++i) {
147			ofs << time_[i] << "\t" <<
148			histogram_[i](0) << "\t" <<
149			histogram_[i](1) << "\t" <<
150			histogram_[i](2) << "\t" << "\n";
151			}
152
153			} else {
154			sprintf(painCave.errMsg,
155			"LegendreCorrFunc::writeCorrelate Error: fail to open %s\n", getOutputFileName().c_str());
156			painCave.isFatal = 1;
157			simError();
158			}
159			ofs.close();
160			}
161			}