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, 9 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)
#	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	#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	: ParticleTimeCorrFunc(info, filename, sele1, sele2, DataStorage::dslAmat \| DataStorage::dslElectroFrame){
49
50	setCorrFuncType("Legendre Correlation Function");
51	setOutputName(getPrefix(dumpFilename_) + ".lcorr");
52	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	LegendrePolynomial polynomial(order);
61	legendre_ = polynomial.getLegendrePolynomial(order);
62	}
63
64	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	}
87
88	}
89
90	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
98	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	}
104
105
106
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	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	"LegendreCorrFunc::validateSelection Error: selected atoms are not Directional\n");
131	painCave.isFatal = 1;
132	simError();
133	}
134	}
135
136	}
137
138	void LegendreCorrFunc::writeCorrelate() {
139	std::ofstream ofs(getOutputFileName().c_str());
140
141	if (ofs.is_open()) {
142
143	ofs << "#" << getCorrFuncType() << "\n";
144	ofs << "#time\tPn(costheta_x)\tPn(costheta_y)\tPn(costheta_z)\n";
145
146	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	}