ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/applications/dynamicProps/LegendreCorrFuncZ.cpp
Revision: 1915
Committed: Mon Jul 29 17:55:17 2013 UTC (12 years, 1 month ago) by gezelter
File size: 6938 byte(s)
Log Message:
Added Legendre Correlation function (as a function of Z), working on architecture for Ewald
Fixed some bugs in FlucQ


File Contents

# User Rev Content
1 gezelter 1915 /*
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. Redistributions of source code must retain the above copyright
10     * notice, this list of conditions and the following disclaimer.
11     *
12     * 2. Redistributions in binary form must reproduce the above copyright
13     * 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     * 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, 234107 (2008).
39     * [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010).
40     * [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41     */
42    
43     #include "applications/dynamicProps/LegendreCorrFuncZ.hpp"
44     #include "math/LegendrePolynomial.hpp"
45     #include "utils/simError.h"
46    
47     namespace OpenMD {
48     LegendreCorrFuncZ::LegendreCorrFuncZ(SimInfo* info,
49     const std::string& filename,
50     const std::string& sele1,
51     const std::string& sele2,
52     int order, int nZbins,
53     long long int memSize)
54     : ParticleTimeCorrFunc(info, filename, sele1, sele2,
55     DataStorage::dslAmat, memSize), nZBins_(nZbins) {
56    
57     setCorrFuncType("Legendre Correlation Function of Z");
58     setOutputName(getPrefix(dumpFilename_) + ".lcorrZ");
59     histogram_.resize(nTimeBins_);
60     counts_.resize(nTimeBins_);
61     for (int i = 0; i < nTimeBins_; i++) {
62     histogram_[i].resize(nZBins_);
63     counts_[i].resize(nZBins_);
64     }
65     LegendrePolynomial polynomial(order);
66     legendre_ = polynomial.getLegendrePolynomial(order);
67     }
68    
69     void LegendreCorrFuncZ::correlateFrames(int frame1, int frame2) {
70     Snapshot* snapshot1 = bsMan_->getSnapshot(frame1);
71     Snapshot* snapshot2 = bsMan_->getSnapshot(frame2);
72     assert(snapshot1 && snapshot2);
73    
74     Mat3x3d hmat = snapshot1->getHmat();
75     RealType halfBoxZ_ = hmat(2,2) / 2.0;
76    
77     RealType time1 = snapshot1->getTime();
78     RealType time2 = snapshot2->getTime();
79    
80     int timeBin = int ((time2 - time1) /deltaTime_ + 0.5);
81    
82     int i;
83     int j;
84     StuntDouble* sd1;
85     StuntDouble* sd2;
86    
87     for (sd1 = seleMan1_.beginSelected(i), sd2 = seleMan2_.beginSelected(j);
88     sd1 != NULL && sd2 != NULL;
89     sd1 = seleMan1_.nextSelected(i), sd2 = seleMan2_.nextSelected(j)) {
90    
91     Vector3d pos = sd1->getPos();
92     if (info_->getSimParams()->getUsePeriodicBoundaryConditions())
93     snapshot1->wrapVector(pos);
94     int zBin = int(nZBins_ * (halfBoxZ_ + pos.z()) / hmat(2,2));
95    
96     Vector3d corrVals = calcCorrVals(frame1, frame2, sd1, sd2);
97     histogram_[timeBin][zBin] += corrVals;
98     counts_[timeBin][zBin]++;
99     }
100    
101     }
102    
103     void LegendreCorrFuncZ::postCorrelate() {
104     for (int i =0 ; i < nTimeBins_; ++i) {
105     for (int j = 0; j < nZBins_; ++j) {
106     if (counts_[i][j] > 0) {
107     histogram_[i][j] /= counts_[i][j];
108     }
109     }
110     }
111     }
112    
113     void LegendreCorrFuncZ::preCorrelate() {
114     for (int i = 0; i < nTimeBins_; i++) {
115     std::fill(histogram_[i].begin(), histogram_[i].end(), Vector3d(0.0));
116     std::fill(counts_[i].begin(), counts_[i].end(), 0);
117     }
118     }
119    
120    
121    
122     Vector3d LegendreCorrFuncZ::calcCorrVals(int frame1, int frame2, StuntDouble* sd1, StuntDouble* sd2) {
123    
124     // The lab frame vector corresponding to the body-fixed
125     // z-axis is simply the second column of A.transpose()
126     // or, identically, the second row of A itself.
127     // Similar identites give the 0th and 1st rows of A for
128     // the lab vectors associated with body-fixed x- and y- axes.
129    
130     Vector3d v1x = sd1->getA(frame1).getRow(0);
131     Vector3d v2x = sd2->getA(frame2).getRow(0);
132     Vector3d v1y = sd1->getA(frame1).getRow(1);
133     Vector3d v2y = sd2->getA(frame2).getRow(1);
134     Vector3d v1z = sd1->getA(frame1).getRow(2);
135     Vector3d v2z = sd2->getA(frame2).getRow(2);
136    
137     RealType uxprod = legendre_.evaluate(dot(v1x, v2x)/(v1x.length()*v2x.length()));
138     RealType uyprod = legendre_.evaluate(dot(v1y, v2y)/(v1y.length()*v2y.length()));
139     RealType uzprod = legendre_.evaluate(dot(v1z, v2z)/(v1z.length()*v2z.length()));
140    
141     return Vector3d(uxprod, uyprod, uzprod);
142    
143     }
144    
145    
146     void LegendreCorrFuncZ::validateSelection(const SelectionManager& seleMan) {
147     StuntDouble* sd;
148     int i;
149     for (sd = seleMan1_.beginSelected(i); sd != NULL; sd = seleMan1_.nextSelected(i)) {
150     if (!sd->isDirectionalAtom()) {
151     sprintf(painCave.errMsg,
152     "LegendreCorrFunc::validateSelection Error: selected atoms are not Directional\n");
153     painCave.isFatal = 1;
154     simError();
155     }
156     }
157    
158     }
159    
160     void LegendreCorrFuncZ::writeCorrelate() {
161     std::ofstream ofs(getOutputFileName().c_str());
162    
163     if (ofs.is_open()) {
164    
165     ofs << "#" << getCorrFuncType() << "\n";
166     ofs << "#time\tPn(costheta_z)\n";
167    
168     for (int i = 0; i < nTimeBins_; ++i) {
169    
170     ofs << time_[i];
171    
172     for (int j = 0; j < nZBins_; ++j) {
173     ofs << "\t" << histogram_[i][j](2);
174     }
175     ofs << "\n";
176     }
177    
178     } else {
179     sprintf(painCave.errMsg,
180     "LegendreCorrFuncZ::writeCorrelate Error: fail to open %s\n", getOutputFileName().c_str());
181     painCave.isFatal = 1;
182     simError();
183     }
184     ofs.close();
185     }
186     }

Properties

Name Value
svn:eol-style native
svn:executable *