applications/dynamicProps/ActionCorrFunc.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, 24107 (2008).          
 * [4]  Vardeman & Gezelter, in progress (2009).                        
 */

#include "applications/dynamicProps/ActionCorrFunc.hpp"
#include "utils/PhysicalConstants.hpp"
#include "brains/ForceManager.hpp"
#include "brains/Thermo.hpp"

namespace OpenMD {

  // We need all of the positions, velocities, etc. so that we can
  // recalculate pressures and actions on the fly:
  ActionCorrFunc::ActionCorrFunc(SimInfo* info, const std::string& filename, 
                                 const std::string& sele1, 
                                 const std::string& sele2, 
                                 long long int memSize)
    : FrameTimeCorrFunc(info, filename, sele1, sele2, 
                        DataStorage::dslPosition | 
                        DataStorage::dslVelocity |
                        DataStorage::dslForce, memSize ){

      setCorrFuncType("ActionCorrFunc");
      setOutputName(getPrefix(dumpFilename_) + ".action");
      histogram_.resize(nTimeBins_); 
      count_.resize(nTimeBins_);
    }

  void ActionCorrFunc::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();
    RealType vol1 = snapshot1->getVolume();
    RealType vol2 = snapshot2->getVolume();
       
    int timeBin = int ((time2 - time1) /deltaTime_ + 0.5);

    //std::cerr << "times = " << time1 << " " << time2 << "\n";
    //std::cerr << "vols = " << vol1 << " " << vol2 << "\n";

    int i;
    int j;

    StuntDouble* sd1;

    Mat3x3d actionTensor1(0.0);
    //std::cerr << "at1 = " << actionTensor1 << "\n";
    Mat3x3d actionTensor2(0.0);
    //std::cerr << "at2 = " << actionTensor2 << "\n";

    for (sd1 = seleMan1_.beginSelected(i); sd1 != NULL;
         sd1 = seleMan1_.nextSelected(i)) {
      //std::cerr << "found a SD\n";
      Vector3d r1 = sd1->getPos(frame1);
      //std::cerr << "r1 = " << r1 << "\n";
      Vector3d v1 = sd1->getVel(frame1);
      //std::cerr << "v1 = " << v1 << "\n";
      Vector3d r2 = sd1->getPos(frame2);
      //std::cerr << "r2 = " << r2 << "\n";
      Vector3d v2 = sd1->getVel(frame2);
      //std::cerr << "v2 = " << v2 << "\n";

      RealType m = sd1->getMass();

      //std::cerr << "m = " << m << "\n";

      actionTensor1 += m*outProduct(r1, v1);
      actionTensor2 += m*outProduct(r2, v2);
    }

    actionTensor1 /= vol1;
    //std::cerr << "at1 = " << actionTensor1 << "\n";
    actionTensor2 /= vol2;
    //std::cerr << "at2 = " << actionTensor2 << "\n";

    Mat3x3d corrTensor(0.0);
    //std::cerr << "ct = " << corrTensor << "\n";
    RealType thisTerm;

    for (i = 0; i < 3; i++) {
      for (j = 0; j < 3; j++) {
       
        //std::cerr << "i, j = " << i << " " << j << "\n"; 
        if (i == j) {
          thisTerm = (actionTensor2(i, j) - actionTensor1(i, j) - avePress_ *(time2-time1));
          std::cerr << "at1, at2 = " << actionTensor1(i,j) << " " << actionTensor2(i,j) << " p = " << avePress_ << "\n";
        } else {
          thisTerm = (actionTensor2(i, j) - actionTensor1(i, j));
        }
        
        //std::cerr << "thisTerm = " << thisTerm << "\n"; 
        corrTensor(i, j) += thisTerm * thisTerm;
      }
    }

    //std::cerr << "ct = " << corrTensor << "\n";
    //std::cerr << "hist = " << histogram_[timeBin] << "\n";
    histogram_[timeBin] += corrTensor;    
    count_[timeBin]++;    
    
  }

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


  void ActionCorrFunc::preCorrelate() {
    // Fill the histogram with empty 3x3 matrices:
    std::fill(histogram_.begin(), histogram_.end(), Mat3x3d(0.0));
    // count array set to zero
    std::fill(count_.begin(), count_.end(), 0);

    SimInfo::MoleculeIterator mi;
    Molecule* mol;
    Molecule::AtomIterator ai;
    Atom* atom;

    // We'll need the force manager to compute forces for the average pressure
    ForceManager* forceMan = new ForceManager(info_);
    forceMan->init();

    // We'll need thermo to compute the pressures from the virial
    Thermo* thermo =  new Thermo(info_);

    // prepare the averages
    RealType pSum = 0.0;
    RealType vSum = 0.0;
    int nsamp = 0;

    // dump files can be enormous, so read them in block-by-block:
    int nblocks = bsMan_->getNBlocks();
    for (int i = 0; i < nblocks; ++i) {
      std::cerr << "block = " << i << "\n";
      bsMan_->loadBlock(i);
      assert(bsMan_->isBlockActive(i));      
      SnapshotBlock block1 = bsMan_->getSnapshotBlock(i);
      for (int j = block1.first; j < block1.second; ++j) {

        // go snapshot-by-snapshot through this block:
        Snapshot* snap = bsMan_->getSnapshot(j);
       
        // update the positions and velocities of the atoms belonging
        // to rigid bodies:

        updateFrame(j);

        // do the forces:
        forceMan->calcForces(true, true);
        // call thermo to get the pressure and volume.
        pSum += thermo->getPressure();
        vSum += thermo->getVolume();
        nsamp++;
        
      }
      bsMan_->unloadBlock(i);
    }

    avePress_ = pSum / ( PhysicalConstants::pressureConvert * (RealType)nsamp);
    aveVol_ = vSum / (RealType)nsamp;
    std::cout << "pAve = " << avePress_ << " vAve = " << aveVol_ << "\n";
  }   


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

    if (ofs.is_open()) {

      ofs << "#" << getCorrFuncType() << "\n";
      ofs << "#time\tcorrTensor\txx\txy\txz\tyx\tyy\tyz\tzx\tzy\tzz\n";

      for (int i = 0; i < nTimeBins_; ++i) {
        ofs << time_[i] << "\t" << 
          histogram_[i](0,0) << "\t" <<
          histogram_[i](0,1) << "\t" <<
          histogram_[i](0,2) << "\t" <<
          histogram_[i](1,0) << "\t" <<
          histogram_[i](1,1) << "\t" <<
          histogram_[i](1,2) << "\t" <<
          histogram_[i](2,0) << "\t" <<
          histogram_[i](2,1) << "\t" <<
          histogram_[i](2,2) << "\t" << "\n";
      }
            
    } else {
      sprintf(painCave.errMsg,
              "ActionCorrFunc::writeCorrelate Error: fail to open %s\n", getOutputFileName().c_str());
      painCave.isFatal = 1;
      simError();        
    }

    ofs.close();    
  }

}
Revision:	1596
Committed:	Mon Jul 25 17:30:53 2011 UTC (14 years, 1 month ago) by gezelter
File size:	8257 byte(s)
Log Message:	Updated the BlockSnapshotManager to use a specified memory footprint in constructor and not to rely on physmem and residentMem to figure out free memory. DynamicProps is the only program that uses the BlockSnapshotManager, so substantial changes were needed there as well.
#	User	Rev	Content
1	xsun	1183	/*
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	gezelter	1390	* 1. Redistributions of source code must retain the above copyright
10	xsun	1183	* notice, this list of conditions and the following disclaimer.
11			*
12	gezelter	1390	* 2. Redistributions in binary form must reproduce the above copyright
13	xsun	1183	* 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	gezelter	1390	*
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, 24107 (2008).
39			* [4] Vardeman & Gezelter, in progress (2009).
40	xsun	1183	*/
41
42			#include "applications/dynamicProps/ActionCorrFunc.hpp"
43	gezelter	1390	#include "utils/PhysicalConstants.hpp"
44	xsun	1183	#include "brains/ForceManager.hpp"
45			#include "brains/Thermo.hpp"
46
47	gezelter	1390	namespace OpenMD {
48	xsun	1183
49			// We need all of the positions, velocities, etc. so that we can
50			// recalculate pressures and actions on the fly:
51			ActionCorrFunc::ActionCorrFunc(SimInfo* info, const std::string& filename,
52			const std::string& sele1,
53	gezelter	1596	const std::string& sele2,
54			long long int memSize)
55	xsun	1183	: FrameTimeCorrFunc(info, filename, sele1, sele2,
56			DataStorage::dslPosition \|
57			DataStorage::dslVelocity \|
58	gezelter	1596	DataStorage::dslForce, memSize ){
59	xsun	1183
60			setCorrFuncType("ActionCorrFunc");
61			setOutputName(getPrefix(dumpFilename_) + ".action");
62			histogram_.resize(nTimeBins_);
63			count_.resize(nTimeBins_);
64			}
65
66			void ActionCorrFunc::correlateFrames(int frame1, int frame2) {
67			Snapshot* snapshot1 = bsMan_->getSnapshot(frame1);
68			Snapshot* snapshot2 = bsMan_->getSnapshot(frame2);
69			assert(snapshot1 && snapshot2);
70
71			RealType time1 = snapshot1->getTime();
72			RealType time2 = snapshot2->getTime();
73			RealType vol1 = snapshot1->getVolume();
74			RealType vol2 = snapshot2->getVolume();
75
76			int timeBin = int ((time2 - time1) /deltaTime_ + 0.5);
77
78			//std::cerr << "times = " << time1 << " " << time2 << "\n";
79			//std::cerr << "vols = " << vol1 << " " << vol2 << "\n";
80
81			int i;
82			int j;
83
84			StuntDouble* sd1;
85
86			Mat3x3d actionTensor1(0.0);
87			//std::cerr << "at1 = " << actionTensor1 << "\n";
88			Mat3x3d actionTensor2(0.0);
89			//std::cerr << "at2 = " << actionTensor2 << "\n";
90
91			for (sd1 = seleMan1_.beginSelected(i); sd1 != NULL;
92			sd1 = seleMan1_.nextSelected(i)) {
93			//std::cerr << "found a SD\n";
94			Vector3d r1 = sd1->getPos(frame1);
95			//std::cerr << "r1 = " << r1 << "\n";
96			Vector3d v1 = sd1->getVel(frame1);
97			//std::cerr << "v1 = " << v1 << "\n";
98			Vector3d r2 = sd1->getPos(frame2);
99			//std::cerr << "r2 = " << r2 << "\n";
100			Vector3d v2 = sd1->getVel(frame2);
101			//std::cerr << "v2 = " << v2 << "\n";
102
103			RealType m = sd1->getMass();
104
105			//std::cerr << "m = " << m << "\n";
106
107			actionTensor1 += m*outProduct(r1, v1);
108			actionTensor2 += m*outProduct(r2, v2);
109			}
110
111			actionTensor1 /= vol1;
112			//std::cerr << "at1 = " << actionTensor1 << "\n";
113			actionTensor2 /= vol2;
114			//std::cerr << "at2 = " << actionTensor2 << "\n";
115
116			Mat3x3d corrTensor(0.0);
117			//std::cerr << "ct = " << corrTensor << "\n";
118			RealType thisTerm;
119
120			for (i = 0; i < 3; i++) {
121			for (j = 0; j < 3; j++) {
122
123			//std::cerr << "i, j = " << i << " " << j << "\n";
124			if (i == j) {
125			thisTerm = (actionTensor2(i, j) - actionTensor1(i, j) - avePress_ *(time2-time1));
126			std::cerr << "at1, at2 = " << actionTensor1(i,j) << " " << actionTensor2(i,j) << " p = " << avePress_ << "\n";
127			} else {
128			thisTerm = (actionTensor2(i, j) - actionTensor1(i, j));
129			}
130
131			//std::cerr << "thisTerm = " << thisTerm << "\n";
132			corrTensor(i, j) += thisTerm * thisTerm;
133			}
134			}
135
136			//std::cerr << "ct = " << corrTensor << "\n";
137			//std::cerr << "hist = " << histogram_[timeBin] << "\n";
138			histogram_[timeBin] += corrTensor;
139			count_[timeBin]++;
140
141			}
142
143			void ActionCorrFunc::postCorrelate() {
144			for (int i =0 ; i < nTimeBins_; ++i) {
145			if (count_[i] > 0) {
146			histogram_[i] /= count_[i];
147			}
148			}
149			}
150
151
152			void ActionCorrFunc::preCorrelate() {
153			// Fill the histogram with empty 3x3 matrices:
154			std::fill(histogram_.begin(), histogram_.end(), Mat3x3d(0.0));
155			// count array set to zero
156			std::fill(count_.begin(), count_.end(), 0);
157
158			SimInfo::MoleculeIterator mi;
159			Molecule* mol;
160			Molecule::AtomIterator ai;
161			Atom* atom;
162
163			// We'll need the force manager to compute forces for the average pressure
164			ForceManager* forceMan = new ForceManager(info_);
165			forceMan->init();
166
167			// We'll need thermo to compute the pressures from the virial
168			Thermo* thermo = new Thermo(info_);
169
170			// prepare the averages
171			RealType pSum = 0.0;
172			RealType vSum = 0.0;
173			int nsamp = 0;
174
175			// dump files can be enormous, so read them in block-by-block:
176			int nblocks = bsMan_->getNBlocks();
177			for (int i = 0; i < nblocks; ++i) {
178			std::cerr << "block = " << i << "\n";
179			bsMan_->loadBlock(i);
180			assert(bsMan_->isBlockActive(i));
181			SnapshotBlock block1 = bsMan_->getSnapshotBlock(i);
182			for (int j = block1.first; j < block1.second; ++j) {
183
184			// go snapshot-by-snapshot through this block:
185			Snapshot* snap = bsMan_->getSnapshot(j);
186
187			// update the positions and velocities of the atoms belonging
188			// to rigid bodies:
189
190			updateFrame(j);
191
192			// do the forces:
193			forceMan->calcForces(true, true);
194			// call thermo to get the pressure and volume.
195			pSum += thermo->getPressure();
196			vSum += thermo->getVolume();
197			nsamp++;
198
199			}
200			bsMan_->unloadBlock(i);
201			}
202
203	gezelter	1390	avePress_ = pSum / ( PhysicalConstants::pressureConvert * (RealType)nsamp);
204	xsun	1183	aveVol_ = vSum / (RealType)nsamp;
205			std::cout << "pAve = " << avePress_ << " vAve = " << aveVol_ << "\n";
206			}
207
208
209			void ActionCorrFunc::writeCorrelate() {
210			std::ofstream ofs(getOutputFileName().c_str());
211
212			if (ofs.is_open()) {
213
214			ofs << "#" << getCorrFuncType() << "\n";
215			ofs << "#time\tcorrTensor\txx\txy\txz\tyx\tyy\tyz\tzx\tzy\tzz\n";
216
217			for (int i = 0; i < nTimeBins_; ++i) {
218			ofs << time_[i] << "\t" <<
219			histogram_[i](0,0) << "\t" <<
220			histogram_[i](0,1) << "\t" <<
221			histogram_[i](0,2) << "\t" <<
222			histogram_[i](1,0) << "\t" <<
223			histogram_[i](1,1) << "\t" <<
224			histogram_[i](1,2) << "\t" <<
225			histogram_[i](2,0) << "\t" <<
226			histogram_[i](2,1) << "\t" <<
227			histogram_[i](2,2) << "\t" << "\n";
228			}
229
230			} else {
231			sprintf(painCave.errMsg,
232			"ActionCorrFunc::writeCorrelate Error: fail to open %s\n", getOutputFileName().c_str());
233			painCave.isFatal = 1;
234			simError();
235			}
236
237			ofs.close();
238			}
239
240			}