applications/dynamicProps/cOHz.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, 234107 (2008).          
 * [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
 * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
 */

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

namespace OpenMD {
  COHZ::COHZ(SimInfo* info, 
             const std::string& filename, 
             const std::string& sele1, 
             const std::string& sele2, 
             int order, int nZbins,
             long long int memSize)
    : ParticleTimeCorrFunc(info, filename, sele1, sele2, 
                           DataStorage::dslAmat, memSize), nZBins_(nZbins) {

      setCorrFuncType("Legendre Correlation Function for OH bond vector of Z");
      setOutputName1(getPrefix(dumpFilename_) + ".cohZ");
      setOutputName2(getPrefix(dumpFilename_) + ".lcorrZ");
      histogram_.resize(nTimeBins_);
      counts_.resize(nTimeBins_);
      for (int i = 0; i < nTimeBins_; i++) {
        histogram_[i].resize(nZBins_);
        counts_[i].resize(nZBins_);
      }
      LegendrePolynomial polynomial(order);
      legendre_ = polynomial.getLegendrePolynomial(order);
    }

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

    Mat3x3d hmat = snapshot1->getHmat();
    RealType halfBoxZ_ = hmat(2,2) / 2.0;      

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

    int timeBin = int ((time2 - time1) /deltaTime_ + 0.5);

    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 pos = sd1->getPos();
      if (info_->getSimParams()->getUsePeriodicBoundaryConditions())
        snapshot1->wrapVector(pos);
      int zBin = int(nZBins_ * (halfBoxZ_ + pos.z()) / hmat(2,2));

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

  void COHZ::postCorrelate() {
    for (int i =0 ; i < nTimeBins_; ++i) {
      for (int j = 0; j < nZBins_; ++j) {
        if (counts_[i][j] > 0) {
          histogram_[i][j] /= counts_[i][j];
        }
      }
    }
  }

  void COHZ::preCorrelate() {
    for (int i = 0; i < nTimeBins_; i++) {
      std::fill(histogram_[i].begin(), histogram_[i].end(), Vector3d(0.0));
      std::fill(counts_[i].begin(), counts_[i].end(), 0);
    }
  }

  Vector3d COHZ::calcCorrVals(int frame1, int frame2, StuntDouble* sd1,  StuntDouble* sd2) {
    
    // Vectors v1x, v1y, and v1z are the body-fixed axes on the
    // molecule in frame 1 in the laboratory frame.

    // Vectors v2x, v2y, and v2z are the body-fixed axes on the
    // molecule in frame 2 in the laboratory frame.

    // Vectors u1 & u2 are the first OH bond vector in frames 1 & 2
    // respectively.  Here we assume SPC/E geometry.

    // Vectors w1 & w2 are the second OH bond vector in frames 1 & 2
    // respectively.  Here we assume SPC/E geometry.

    Vector3d v1x = sd1->getA(frame1).getRow(0);
    Vector3d v2x = sd2->getA(frame2).getRow(0);

    Vector3d v1y = sd1->getA(frame1).getRow(1);
    Vector3d v2y = sd2->getA(frame2).getRow(1);

    Vector3d v1z = sd1->getA(frame1).getRow(2);
    Vector3d v2z = sd2->getA(frame2).getRow(2);

    Vector3d u1 = 0.81649 * v1y + 0.57736 * v1z;
    Vector3d u2 = 0.81649 * v2y + 0.57736 * v2z;

    Vector3d w1 = -0.81649 * v1y + 0.57736 * v1z;
    Vector3d w2 = -0.81649 * v2y + 0.57736 * v2z;

    
    RealType vprod = legendre_.evaluate(dot(v1z, v2z)/(v1z.length()*v2z.length()));
    RealType uprod = legendre_.evaluate(dot(u1, u2)/(u1.length()*u2.length()));
    RealType wprod = legendre_.evaluate(dot(w1, w2)/(w1.length()*w2.length()));

    return Vector3d(vprod, uprod, wprod);

  }


  void COHZ::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 COHZ::writeCorrelate() {
    std::ofstream ofs1(getOutputFileName1().c_str());
    std::ofstream ofs2(getOutputFileName2().c_str());

    if (ofs1.is_open()) {

      ofs1 << "#" << getCorrFuncType() << "\n";
      ofs1 << "#time\tPn(costheta_z)\n";

      for (int i = 0; i < nTimeBins_; ++i) {

        ofs1 << time_[i];

        for (int j = 0; j < nZBins_; ++j) {          
          ofs1 << "\t" << 0.5*(histogram_[i][j](1) +  histogram_[i][j](2));
        }
        ofs1 << "\n";
      }
            
    } else {
      sprintf(painCave.errMsg,
              "cOHz::writeCorrelate Error: fail to open %s\n", getOutputFileName1().c_str());
      painCave.isFatal = 1;
      simError();        
    }
    ofs1.close();    

    if (ofs2.is_open()) {

      ofs2 << "#" << getCorrFuncType() << "\n";
      ofs2 << "#time\tPn(costheta_z)\n";

      for (int i = 0; i < nTimeBins_; ++i) {

        ofs2 << time_[i];

        for (int j = 0; j < nZBins_; ++j) {          
          ofs2 << "\t" << histogram_[i][j](0);
        }
        ofs2 << "\n";
      }
            
    } else {
      sprintf(painCave.errMsg,
              "cOHz::writeCorrelate Error: fail to open %s\n", getOutputFileName2().c_str());
      painCave.isFatal = 1;
      simError();        
    }
    ofs2.close();    
  }
}
Revision:	1934
Committed:	Mon Aug 26 14:15:09 2013 UTC (12 years, 2 months ago) by gezelter
File size:	7785 byte(s)
Log Message:	Added OH bond time correlation function
#	User	Rev	Content
1	gezelter	1934	/*
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/cOHz.hpp"
44			#include "math/LegendrePolynomial.hpp"
45			#include "utils/simError.h"
46
47			namespace OpenMD {
48			COHZ::COHZ(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 for OH bond vector of Z");
58			setOutputName1(getPrefix(dumpFilename_) + ".cohZ");
59			setOutputName2(getPrefix(dumpFilename_) + ".lcorrZ");
60			histogram_.resize(nTimeBins_);
61			counts_.resize(nTimeBins_);
62			for (int i = 0; i < nTimeBins_; i++) {
63			histogram_[i].resize(nZBins_);
64			counts_[i].resize(nZBins_);
65			}
66			LegendrePolynomial polynomial(order);
67			legendre_ = polynomial.getLegendrePolynomial(order);
68			}
69
70			void COHZ::correlateFrames(int frame1, int frame2) {
71			Snapshot* snapshot1 = bsMan_->getSnapshot(frame1);
72			Snapshot* snapshot2 = bsMan_->getSnapshot(frame2);
73			assert(snapshot1 && snapshot2);
74
75			Mat3x3d hmat = snapshot1->getHmat();
76			RealType halfBoxZ_ = hmat(2,2) / 2.0;
77
78			RealType time1 = snapshot1->getTime();
79			RealType time2 = snapshot2->getTime();
80
81			int timeBin = int ((time2 - time1) /deltaTime_ + 0.5);
82
83			int i;
84			int j;
85			StuntDouble* sd1;
86			StuntDouble* sd2;
87
88			for (sd1 = seleMan1_.beginSelected(i), sd2 = seleMan2_.beginSelected(j);
89			sd1 != NULL && sd2 != NULL;
90			sd1 = seleMan1_.nextSelected(i), sd2 = seleMan2_.nextSelected(j)) {
91
92			Vector3d pos = sd1->getPos();
93			if (info_->getSimParams()->getUsePeriodicBoundaryConditions())
94			snapshot1->wrapVector(pos);
95			int zBin = int(nZBins_ * (halfBoxZ_ + pos.z()) / hmat(2,2));
96
97			Vector3d corrVals = calcCorrVals(frame1, frame2, sd1, sd2);
98			histogram_[timeBin][zBin] += corrVals;
99			counts_[timeBin][zBin]++;
100			}
101
102			}
103
104			void COHZ::postCorrelate() {
105			for (int i =0 ; i < nTimeBins_; ++i) {
106			for (int j = 0; j < nZBins_; ++j) {
107			if (counts_[i][j] > 0) {
108			histogram_[i][j] /= counts_[i][j];
109			}
110			}
111			}
112			}
113
114			void COHZ::preCorrelate() {
115			for (int i = 0; i < nTimeBins_; i++) {
116			std::fill(histogram_[i].begin(), histogram_[i].end(), Vector3d(0.0));
117			std::fill(counts_[i].begin(), counts_[i].end(), 0);
118			}
119			}
120
121			Vector3d COHZ::calcCorrVals(int frame1, int frame2, StuntDouble* sd1, StuntDouble* sd2) {
122
123			// Vectors v1x, v1y, and v1z are the body-fixed axes on the
124			// molecule in frame 1 in the laboratory frame.
125
126			// Vectors v2x, v2y, and v2z are the body-fixed axes on the
127			// molecule in frame 2 in the laboratory frame.
128
129			// Vectors u1 & u2 are the first OH bond vector in frames 1 & 2
130			// respectively. Here we assume SPC/E geometry.
131
132			// Vectors w1 & w2 are the second OH bond vector in frames 1 & 2
133			// respectively. Here we assume SPC/E geometry.
134
135			Vector3d v1x = sd1->getA(frame1).getRow(0);
136			Vector3d v2x = sd2->getA(frame2).getRow(0);
137
138			Vector3d v1y = sd1->getA(frame1).getRow(1);
139			Vector3d v2y = sd2->getA(frame2).getRow(1);
140
141			Vector3d v1z = sd1->getA(frame1).getRow(2);
142			Vector3d v2z = sd2->getA(frame2).getRow(2);
143
144			Vector3d u1 = 0.81649 * v1y + 0.57736 * v1z;
145			Vector3d u2 = 0.81649 * v2y + 0.57736 * v2z;
146
147			Vector3d w1 = -0.81649 * v1y + 0.57736 * v1z;
148			Vector3d w2 = -0.81649 * v2y + 0.57736 * v2z;
149
150
151			RealType vprod = legendre_.evaluate(dot(v1z, v2z)/(v1z.length()*v2z.length()));
152			RealType uprod = legendre_.evaluate(dot(u1, u2)/(u1.length()*u2.length()));
153			RealType wprod = legendre_.evaluate(dot(w1, w2)/(w1.length()*w2.length()));
154
155			return Vector3d(vprod, uprod, wprod);
156
157			}
158
159
160			void COHZ::validateSelection(const SelectionManager& seleMan) {
161			StuntDouble* sd;
162			int i;
163			for (sd = seleMan1_.beginSelected(i); sd != NULL; sd = seleMan1_.nextSelected(i)) {
164			if (!sd->isDirectionalAtom()) {
165			sprintf(painCave.errMsg,
166			"LegendreCorrFunc::validateSelection Error: selected atoms are not Directional\n");
167			painCave.isFatal = 1;
168			simError();
169			}
170			}
171
172			}
173
174			void COHZ::writeCorrelate() {
175			std::ofstream ofs1(getOutputFileName1().c_str());
176			std::ofstream ofs2(getOutputFileName2().c_str());
177
178			if (ofs1.is_open()) {
179
180			ofs1 << "#" << getCorrFuncType() << "\n";
181			ofs1 << "#time\tPn(costheta_z)\n";
182
183			for (int i = 0; i < nTimeBins_; ++i) {
184
185			ofs1 << time_[i];
186
187			for (int j = 0; j < nZBins_; ++j) {
188			ofs1 << "\t" << 0.5*(histogram_[i][j](1) + histogram_[i][j](2));
189			}
190			ofs1 << "\n";
191			}
192
193			} else {
194			sprintf(painCave.errMsg,
195			"cOHz::writeCorrelate Error: fail to open %s\n", getOutputFileName1().c_str());
196			painCave.isFatal = 1;
197			simError();
198			}
199			ofs1.close();
200
201			if (ofs2.is_open()) {
202
203			ofs2 << "#" << getCorrFuncType() << "\n";
204			ofs2 << "#time\tPn(costheta_z)\n";
205
206			for (int i = 0; i < nTimeBins_; ++i) {
207
208			ofs2 << time_[i];
209
210			for (int j = 0; j < nZBins_; ++j) {
211			ofs2 << "\t" << histogram_[i][j](0);
212			}
213			ofs2 << "\n";
214			}
215
216			} else {
217			sprintf(painCave.errMsg,
218			"cOHz::writeCorrelate Error: fail to open %s\n", getOutputFileName2().c_str());
219			painCave.isFatal = 1;
220			simError();
221			}
222			ofs2.close();
223			}
224			}