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

 /* Uses the Helfand-moment method for calculating thermal
  * conductivity using the relation kappa = (N,V)lim(t)->inf 1/(2*k_B*T^2*V*t) <[G_K(t)-G_K(0)]^2>
  * where G_K is the Helfand moment for thermal conductivity definded as
  * G_K(t) = sum_{a=1}{^N} x_a(E_a-<E_a>) and E_a is defined to be
  *     E_a = p_2^2/(2*m)+1/2 sum_{b.ne.a} u(r_ab) where p is momentum and u is pot energy for the 
  * particle pair a-b. This routine calculates E_a, <E_a> and does the correlation
  * <[G_K(t)-G_K(0)]^2>.
  * See Viscardy et al. JCP 126, 184513 (2007)
  */


#include "applications/dynamicProps/EnergyCorrFunc.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:
  EnergyCorrFunc::EnergyCorrFunc(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 |
                        DataStorage::dslTorque |                        
                        DataStorage::dslParticlePot,
                        memSize){

      setCorrFuncType("EnergyCorrFunc");
      setOutputName(getPrefix(dumpFilename_) + ".moment");
      histogram_.resize(nTimeBins_); 
      count_.resize(nTimeBins_);
    }

  void EnergyCorrFunc::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);

    Vector3d G_t_frame1 = G_t_[frame1];
    Vector3d G_t_frame2 = G_t_[frame2];
    
    
    RealType diff_x = G_t_frame1.x()-G_t_frame2.x();
    RealType diff_x_sq = diff_x * diff_x;
    
    RealType diff_y = G_t_frame1.y()-G_t_frame2.y();
    RealType diff_y_sq = diff_y * diff_y;
    
    RealType diff_z = G_t_frame1.z()-G_t_frame2.z();
    RealType diff_z_sq = diff_z*diff_z;    
    
    histogram_[timeBin][0] += diff_x_sq;
    histogram_[timeBin][1] += diff_y_sq;
    histogram_[timeBin][2] += diff_z_sq;
    
    count_[timeBin]++;
    
  }

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

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

    SimInfo::MoleculeIterator mi;
    Molecule::IntegrableObjectIterator mj;
    Molecule* mol;
    Molecule::AtomIterator ai;
    Atom* atom;
    std::vector<RealType > particleEnergies;

    // 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();
    bool firsttime = true;
    int junkframe = 0;
    for (int i = 0; i < nblocks; ++i) {
      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(); 
        
        int index = 0;
        
        for (mol = info_->beginMolecule(mi); mol != NULL; 
             mol = info_->nextMolecule(mi)) {
          for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
            RealType mass = atom->getMass();
            Vector3d vel = atom->getVel();
            RealType kinetic = mass * (vel[0]*vel[0] + vel[1]*vel[1] + 
                                       vel[2]*vel[2]) / PhysicalConstants::energyConvert;
            RealType potential =  atom->getParticlePot();
            RealType eatom = (kinetic + potential)/2.0;
            particleEnergies.push_back(eatom);
            if(firsttime)
              {
                AvgE_a_.push_back(eatom);
              } else {
              /* We assume the the number of atoms does not change.*/
              AvgE_a_[index] += eatom;
            }
            index++;
          }      
        }
        firsttime = false;
        E_a_.push_back(particleEnergies);
      }
      
      bsMan_->unloadBlock(i);
    }
    
    int nframes =  bsMan_->getNFrames();
    for (int i = 0; i < AvgE_a_.size(); i++){
      AvgE_a_[i] /= nframes;
    }
    
    int frame = 0;
    
    // Do it again to compute G^(kappa)(t) for x,y,z
    for (int i = 0; i < nblocks; ++i) {
      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);        

        // this needs to be updated to the frame value:
        particleEnergies = E_a_[j];
        
        int thisAtom = 0;
        Vector3d G_t;
        
        for (mol = info_->beginMolecule(mi); mol != NULL; 
             mol = info_->nextMolecule(mi)) {
          for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
                        
            Vector3d pos = atom->getPos();

            G_t[0] += pos.x()*(particleEnergies[thisAtom]-AvgE_a_[thisAtom]);
            G_t[1] += pos.y()*(particleEnergies[thisAtom]-AvgE_a_[thisAtom]);
            G_t[2] += pos.z()*(particleEnergies[thisAtom]-AvgE_a_[thisAtom]);
            
            thisAtom++;                    
          }
        }

        G_t_.push_back(G_t);
        //std::cerr <<"Frame: " << j <<"\t" << G_t << std::endl;
      } 
      bsMan_->unloadBlock(i);
    }
  }   


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

    if (ofs.is_open()) {

      ofs << "#" << getCorrFuncType() << "\n";
      ofs << "#time\tK_x\tK_y\tK_z\n";

      for (int i = 0; i < nTimeBins_; ++i) {
        ofs << time_[i] << "\t" << 
              histogram_[i].x() << "\t" <<
              histogram_[i].y() << "\t" <<
              histogram_[i].z() << "\t" << "\n";
      }
            
    } else {
      sprintf(painCave.errMsg,
              "EnergyCorrFunc::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:	9400 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	chuckv	1246	/*
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	chuckv	1246	* notice, this list of conditions and the following disclaimer.
11			*
12	gezelter	1390	* 2. Redistributions in binary form must reproduce the above copyright
13	chuckv	1246	* 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	chuckv	1246	*/
41
42			/* Uses the Helfand-moment method for calculating thermal
43			* conductivity using the relation kappa = (N,V)lim(t)->inf 1/(2k_BT^2Vt) <[G_K(t)-G_K(0)]^2>
44			* where G_K is the Helfand moment for thermal conductivity definded as
45			* G_K(t) = sum_{a=1}{^N} x_a(E_a-<E_a>) and E_a is defined to be
46			* E_a = p_2^2/(2*m)+1/2 sum_{b.ne.a} u(r_ab) where p is momentum and u is pot energy for the
47			* particle pair a-b. This routine calculates E_a, <E_a> and does the correlation
48			* <[G_K(t)-G_K(0)]^2>.
49			* See Viscardy et al. JCP 126, 184513 (2007)
50			*/
51
52
53
54			#include "applications/dynamicProps/EnergyCorrFunc.hpp"
55	gezelter	1390	#include "utils/PhysicalConstants.hpp"
56	chuckv	1246	#include "brains/ForceManager.hpp"
57			#include "brains/Thermo.hpp"
58
59	gezelter	1390	namespace OpenMD {
60	chuckv	1246
61			// We need all of the positions, velocities, etc. so that we can
62			// recalculate pressures and actions on the fly:
63			EnergyCorrFunc::EnergyCorrFunc(SimInfo* info, const std::string& filename,
64			const std::string& sele1,
65	gezelter	1596	const std::string& sele2,
66			long long int memSize)
67	chuckv	1246	: FrameTimeCorrFunc(info, filename, sele1, sele2,
68			DataStorage::dslPosition \|
69			DataStorage::dslVelocity \|
70			DataStorage::dslForce \|
71			DataStorage::dslTorque \|
72	gezelter	1596	DataStorage::dslParticlePot,
73			memSize){
74	chuckv	1246
75			setCorrFuncType("EnergyCorrFunc");
76			setOutputName(getPrefix(dumpFilename_) + ".moment");
77			histogram_.resize(nTimeBins_);
78			count_.resize(nTimeBins_);
79			}
80
81			void EnergyCorrFunc::correlateFrames(int frame1, int frame2) {
82			Snapshot* snapshot1 = bsMan_->getSnapshot(frame1);
83			Snapshot* snapshot2 = bsMan_->getSnapshot(frame2);
84			assert(snapshot1 && snapshot2);
85
86			RealType time1 = snapshot1->getTime();
87			RealType time2 = snapshot2->getTime();
88
89			int timeBin = int ((time2 - time1) /deltaTime_ + 0.5);
90
91	gezelter	1313	Vector3d G_t_frame1 = G_t_[frame1];
92			Vector3d G_t_frame2 = G_t_[frame2];
93
94
95			RealType diff_x = G_t_frame1.x()-G_t_frame2.x();
96			RealType diff_x_sq = diff_x * diff_x;
97
98			RealType diff_y = G_t_frame1.y()-G_t_frame2.y();
99			RealType diff_y_sq = diff_y * diff_y;
100
101			RealType diff_z = G_t_frame1.z()-G_t_frame2.z();
102			RealType diff_z_sq = diff_z*diff_z;
103
104			histogram_[timeBin][0] += diff_x_sq;
105			histogram_[timeBin][1] += diff_y_sq;
106			histogram_[timeBin][2] += diff_z_sq;
107
108			count_[timeBin]++;
109
110	chuckv	1246	}
111
112			void EnergyCorrFunc::postCorrelate() {
113			for (int i =0 ; i < nTimeBins_; ++i) {
114			if (count_[i] > 0) {
115			histogram_[i] /= count_[i];
116			}
117			}
118			}
119
120			void EnergyCorrFunc::preCorrelate() {
121			// Fill the histogram with empty 3x3 matrices:
122			std::fill(histogram_.begin(), histogram_.end(), 0.0);
123			// count array set to zero
124			std::fill(count_.begin(), count_.end(), 0);
125
126			SimInfo::MoleculeIterator mi;
127			Molecule::IntegrableObjectIterator mj;
128			Molecule* mol;
129			Molecule::AtomIterator ai;
130			Atom* atom;
131			std::vector<RealType > particleEnergies;
132
133			// We'll need the force manager to compute forces for the average pressure
134			ForceManager* forceMan = new ForceManager(info_);
135
136			forceMan->init();
137	gezelter	1313
138	chuckv	1246	// We'll need thermo to compute the pressures from the virial
139			Thermo* thermo = new Thermo(info_);
140
141			// prepare the averages
142			RealType pSum = 0.0;
143			RealType vSum = 0.0;
144			int nsamp = 0;
145
146			// dump files can be enormous, so read them in block-by-block:
147			int nblocks = bsMan_->getNBlocks();
148			bool firsttime = true;
149	chuckv	1247	int junkframe = 0;
150	chuckv	1246	for (int i = 0; i < nblocks; ++i) {
151			bsMan_->loadBlock(i);
152			assert(bsMan_->isBlockActive(i));
153			SnapshotBlock block1 = bsMan_->getSnapshotBlock(i);
154			for (int j = block1.first; j < block1.second; ++j) {
155
156			// go snapshot-by-snapshot through this block:
157			Snapshot* snap = bsMan_->getSnapshot(j);
158	gezelter	1249
159	chuckv	1246	// update the positions and velocities of the atoms belonging
160			// to rigid bodies:
161	gezelter	1249
162	chuckv	1246	updateFrame(j);
163	gezelter	1249
164	chuckv	1246	// do the forces:
165	gezelter	1249
166	gezelter	1464	forceMan->calcForces();
167	gezelter	1249
168	chuckv	1246	int index = 0;
169	gezelter	1249
170	chuckv	1246	for (mol = info_->beginMolecule(mi); mol != NULL;
171	gezelter	1249	mol = info_->nextMolecule(mi)) {
172	chuckv	1246	for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
173			RealType mass = atom->getMass();
174			Vector3d vel = atom->getVel();
175	gezelter	1249	RealType kinetic = mass * (vel[0]vel[0] + vel[1]vel[1] +
176	gezelter	1390	vel[2]*vel[2]) / PhysicalConstants::energyConvert;
177	gezelter	1313	RealType potential = atom->getParticlePot();
178			RealType eatom = (kinetic + potential)/2.0;
179	chuckv	1246	particleEnergies.push_back(eatom);
180			if(firsttime)
181	gezelter	1249	{
182			AvgE_a_.push_back(eatom);
183			} else {
184	chuckv	1246	/* We assume the the number of atoms does not change.*/
185			AvgE_a_[index] += eatom;
186			}
187			index++;
188			}
189			}
190			firsttime = false;
191			E_a_.push_back(particleEnergies);
192			}
193	chuckv	1247
194	chuckv	1246	bsMan_->unloadBlock(i);
195			}
196
197	gezelter	1249	int nframes = bsMan_->getNFrames();
198			for (int i = 0; i < AvgE_a_.size(); i++){
199			AvgE_a_[i] /= nframes;
200			}
201	chuckv	1246
202	gezelter	1249	int frame = 0;
203	chuckv	1246
204	gezelter	1249	// Do it again to compute G^(kappa)(t) for x,y,z
205			for (int i = 0; i < nblocks; ++i) {
206			bsMan_->loadBlock(i);
207			assert(bsMan_->isBlockActive(i));
208			SnapshotBlock block1 = bsMan_->getSnapshotBlock(i);
209			for (int j = block1.first; j < block1.second; ++j) {
210	chuckv	1246
211	gezelter	1249	// go snapshot-by-snapshot through this block:
212			Snapshot* snap = bsMan_->getSnapshot(j);
213
214			// update the positions and velocities of the atoms belonging
215			// to rigid bodies:
216
217			updateFrame(j);
218
219			// this needs to be updated to the frame value:
220			particleEnergies = E_a_[j];
221
222			int thisAtom = 0;
223			Vector3d G_t;
224
225			for (mol = info_->beginMolecule(mi); mol != NULL;
226			mol = info_->nextMolecule(mi)) {
227			for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
228
229			Vector3d pos = atom->getPos();
230
231			G_t[0] += pos.x()*(particleEnergies[thisAtom]-AvgE_a_[thisAtom]);
232			G_t[1] += pos.y()*(particleEnergies[thisAtom]-AvgE_a_[thisAtom]);
233			G_t[2] += pos.z()*(particleEnergies[thisAtom]-AvgE_a_[thisAtom]);
234
235			thisAtom++;
236			}
237			}
238
239			G_t_.push_back(G_t);
240	gezelter	1313	//std::cerr <<"Frame: " << j <<"\t" << G_t << std::endl;
241	gezelter	1249	}
242			bsMan_->unloadBlock(i);
243			}
244	chuckv	1246	}
245
246
247			void EnergyCorrFunc::writeCorrelate() {
248			std::ofstream ofs(getOutputFileName().c_str());
249
250			if (ofs.is_open()) {
251
252			ofs << "#" << getCorrFuncType() << "\n";
253			ofs << "#time\tK_x\tK_y\tK_z\n";
254
255			for (int i = 0; i < nTimeBins_; ++i) {
256			ofs << time_[i] << "\t" <<
257			histogram_[i].x() << "\t" <<
258			histogram_[i].y() << "\t" <<
259			histogram_[i].z() << "\t" << "\n";
260			}
261
262			} else {
263			sprintf(painCave.errMsg,
264			"EnergyCorrFunc::writeCorrelate Error: fail to open %s\n", getOutputFileName().c_str());
265			painCave.isFatal = 1;
266			simError();
267			}
268
269			ofs.close();
270
271			}
272
273			}