applications/dynamicProps/TimeCorrFunc.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/TimeCorrFunc.hpp"
#include "utils/simError.h"
#include "primitives/Molecule.hpp"
namespace OpenMD {

  TimeCorrFunc::TimeCorrFunc(SimInfo* info, const std::string& filename, 
                             const std::string& sele1, const std::string& sele2,
                             int storageLayout, long long int memSize)
    : info_(info), storageLayout_(storageLayout), memSize_(memSize),
      dumpFilename_(filename), selectionScript1_(sele1), 
      selectionScript2_(sele2), evaluator1_(info), evaluator2_(info), 
      seleMan1_(info), seleMan2_(info){

      int nAtoms = info->getNGlobalAtoms();
      int nRigidBodies = info->getNGlobalRigidBodies();
      int nStuntDoubles =   nAtoms + nRigidBodies;

      std::set<AtomType*> atomTypes = info->getUniqueAtomTypes();
      std::set<AtomType*>::iterator i;
      bool hasDirectionalAtom = false;
      bool hasMultipole = false;    
      for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
        if ((*i)->isDirectional()){
          hasDirectionalAtom = true;
        }
        if ((*i)->isMultipole()){
          hasMultipole = true;
        }
      }
    
      if (nRigidBodies > 0 || hasDirectionalAtom) {
        storageLayout_ |= DataStorage::dslAmat;
        if(storageLayout_ & DataStorage::dslVelocity) {
          storageLayout_ |= DataStorage::dslAngularMomentum;
        }
        if (storageLayout_ & DataStorage::dslForce) {
          storageLayout_ |= DataStorage::dslTorque;
        }
      }
      if (hasMultipole) {
        storageLayout_ |= DataStorage::dslElectroFrame;
      }
        
      bsMan_ = new BlockSnapshotManager(info, dumpFilename_, storageLayout_, memSize_);
      info_->setSnapshotManager(bsMan_);

      evaluator1_.loadScriptString(selectionScript1_);
      evaluator2_.loadScriptString(selectionScript2_);
    
      //if selection is static, we only need to evaluate it once
      if (!evaluator1_.isDynamic()) {
        seleMan1_.setSelectionSet(evaluator1_.evaluate());
        validateSelection(seleMan1_);
      } else {
        sprintf(painCave.errMsg,
                "TimeCorrFunc Error: dynamic selection is not supported\n");
        painCave.isFatal = 1;
        simError();  
      }
    
      if (!evaluator2_.isDynamic()) {
        seleMan2_.setSelectionSet(evaluator2_.evaluate());
        validateSelection(seleMan2_);
      } else {
        sprintf(painCave.errMsg,
                "TimeCorrFunc Error: dynamic selection is not supported\n");
        painCave.isFatal = 1;
        simError();  
      }
    

      /**@todo Fix Me */
      Globals* simParams = info_->getSimParams();
      if (simParams->haveSampleTime()){
        deltaTime_ = simParams->getSampleTime();
      } else {
        sprintf(painCave.errMsg,
                "TimeCorrFunc::writeCorrelate Error: can not figure out deltaTime\n");
        painCave.isFatal = 1;
        simError();  
      }

      int nframes =  bsMan_->getNFrames();
      nTimeBins_ = nframes;
      histogram_.resize(nTimeBins_);
      count_.resize(nTimeBins_);
      time_.resize(nTimeBins_);

      for (int i = 0; i < nTimeBins_; ++i) {
        time_[i] = i * deltaTime_;
      }
    
    }


  void TimeCorrFunc::doCorrelate() {
    preCorrelate();

    int nblocks = bsMan_->getNBlocks();

    for (int i = 0; i < nblocks; ++i) {
      bsMan_->loadBlock(i);

      for (int j = i; j < nblocks; ++j) {
        bsMan_->loadBlock(j);
        correlateBlocks(i, j);
        bsMan_->unloadBlock(j);
      }
        
      bsMan_->unloadBlock(i);
    }
    
    postCorrelate();

    writeCorrelate();
  }

  void TimeCorrFunc::correlateBlocks(int block1, int block2) {

    int jstart, jend;

    assert(bsMan_->isBlockActive(block1) && bsMan_->isBlockActive(block2));

    SnapshotBlock snapshotBlock1 = bsMan_->getSnapshotBlock(block1);
    SnapshotBlock snapshotBlock2 = bsMan_->getSnapshotBlock(block2);

    jend = snapshotBlock2.second;

    for (int i = snapshotBlock1.first; i < snapshotBlock1.second; ++i) {

      //update the position or velocity of the atoms belong to rigid bodies
      updateFrame(i);

      // if the two blocks are the same, we don't want to correlate
      // backwards in time, so start j at the same frame as i
      if (block1 == block2) {
        jstart = i;
      } else {
        jstart = snapshotBlock2.first;
      }
        
      for(int j  = jstart; j < jend; ++j) {
        //update the position or velocity of the atoms belong to rigid bodies
        updateFrame(j);

        correlateFrames(i, j);
      }
    }
  }

  void TimeCorrFunc::updateFrame(int frame){
    Molecule* mol;
    RigidBody* rb;
    SimInfo::MoleculeIterator mi;
    Molecule::RigidBodyIterator rbIter;

    /** @todo need improvement */    
    if (storageLayout_ & DataStorage::dslPosition) {
      for (mol = info_->beginMolecule(mi); mol != NULL; 
           mol = info_->nextMolecule(mi)) {

        //change the positions of atoms which belong to the rigidbodies
        for (rb = mol->beginRigidBody(rbIter); rb != NULL; 
             rb = mol->nextRigidBody(rbIter)) {
          rb->updateAtoms(frame);
        }
      }        
    }

    if (storageLayout_ & DataStorage::dslVelocity) {
      for (mol = info_->beginMolecule(mi); mol != NULL; 
           mol = info_->nextMolecule(mi)) {
        
        //change the positions of atoms which belong to the rigidbodies
        for (rb = mol->beginRigidBody(rbIter); rb != NULL; 
             rb = mol->nextRigidBody(rbIter)) {
          rb->updateAtomVel(frame);
        }
      }
      
    }
    
  }


  void TimeCorrFunc::preCorrelate() {
    std::fill(histogram_.begin(), histogram_.end(), 0.0);
    std::fill(count_.begin(), count_.end(), 0);
  }

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


  void TimeCorrFunc::writeCorrelate() {
    std::ofstream ofs(outputFilename_.c_str());

    if (ofs.is_open()) {

      ofs << "#" << getCorrFuncType() << "\n";
      ofs << "#selection script1: \"" << selectionScript1_ <<"\"\tselection script2: \"" << selectionScript2_ << "\"\n";
      ofs << "#extra information: " << extra_ << "\n";
      ofs << "#time\tcorrVal\n";

      for (int i = 0; i < nTimeBins_; ++i) {
        ofs << time_[i] << "\t" << histogram_[i] << "\n";
      }
            
    } else {
      sprintf(painCave.errMsg,
              "TimeCorrFunc::writeCorrelate Error: fail to open %s\n", outputFilename_.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:	8463 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	tim	333	/*
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	tim	333	* notice, this list of conditions and the following disclaimer.
11			*
12	gezelter	1390	* 2. Redistributions in binary form must reproduce the above copyright
13	tim	333	* 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	tim	333	*/
41
42			#include "applications/dynamicProps/TimeCorrFunc.hpp"
43			#include "utils/simError.h"
44			#include "primitives/Molecule.hpp"
45	gezelter	1390	namespace OpenMD {
46	tim	333
47	gezelter	507	TimeCorrFunc::TimeCorrFunc(SimInfo* info, const std::string& filename,
48	gezelter	1596	const std::string& sele1, const std::string& sele2,
49			int storageLayout, long long int memSize)
50			: info_(info), storageLayout_(storageLayout), memSize_(memSize),
51			dumpFilename_(filename), selectionScript1_(sele1),
52			selectionScript2_(sele2), evaluator1_(info), evaluator2_(info),
53			seleMan1_(info), seleMan2_(info){
54	tim	333
55	gezelter	507	int nAtoms = info->getNGlobalAtoms();
56			int nRigidBodies = info->getNGlobalRigidBodies();
57			int nStuntDoubles = nAtoms + nRigidBodies;
58	tim	333
59	gezelter	507	std::set<AtomType*> atomTypes = info->getUniqueAtomTypes();
60			std::set<AtomType*>::iterator i;
61			bool hasDirectionalAtom = false;
62			bool hasMultipole = false;
63			for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
64	tim	333	if ((*i)->isDirectional()){
65	gezelter	507	hasDirectionalAtom = true;
66	tim	333	}
67			if ((*i)->isMultipole()){
68	gezelter	507	hasMultipole = true;
69	tim	333	}
70	gezelter	507	}
71	tim	333
72	gezelter	507	if (nRigidBodies > 0 \|\| hasDirectionalAtom) {
73	tim	333	storageLayout_ \|= DataStorage::dslAmat;
74	xsun	1183	if(storageLayout_ & DataStorage::dslVelocity) {
75			storageLayout_ \|= DataStorage::dslAngularMomentum;
76			}
77			if (storageLayout_ & DataStorage::dslForce) {
78			storageLayout_ \|= DataStorage::dslTorque;
79			}
80	gezelter	507	}
81			if (hasMultipole) {
82	tim	333	storageLayout_ \|= DataStorage::dslElectroFrame;
83	gezelter	507	}
84	tim	333
85	gezelter	1596	bsMan_ = new BlockSnapshotManager(info, dumpFilename_, storageLayout_, memSize_);
86	gezelter	507	info_->setSnapshotManager(bsMan_);
87	tim	333
88	gezelter	507	evaluator1_.loadScriptString(selectionScript1_);
89			evaluator2_.loadScriptString(selectionScript2_);
90	tim	334
91	gezelter	507	//if selection is static, we only need to evaluate it once
92			if (!evaluator1_.isDynamic()) {
93	tim	333	seleMan1_.setSelectionSet(evaluator1_.evaluate());
94			validateSelection(seleMan1_);
95	gezelter	507	} else {
96			sprintf(painCave.errMsg,
97	tim	333	"TimeCorrFunc Error: dynamic selection is not supported\n");
98	gezelter	507	painCave.isFatal = 1;
99			simError();
100			}
101	tim	333
102	gezelter	507	if (!evaluator2_.isDynamic()) {
103	tim	333	seleMan2_.setSelectionSet(evaluator2_.evaluate());
104			validateSelection(seleMan2_);
105	gezelter	507	} else {
106			sprintf(painCave.errMsg,
107	tim	333	"TimeCorrFunc Error: dynamic selection is not supported\n");
108	gezelter	507	painCave.isFatal = 1;
109			simError();
110			}
111	tim	333
112
113
114	gezelter	507	/*@todo Fix Me /
115			Globals* simParams = info_->getSimParams();
116			if (simParams->haveSampleTime()){
117	tim	333	deltaTime_ = simParams->getSampleTime();
118	gezelter	507	} else {
119			sprintf(painCave.errMsg,
120	tim	333	"TimeCorrFunc::writeCorrelate Error: can not figure out deltaTime\n");
121	gezelter	507	painCave.isFatal = 1;
122			simError();
123			}
124	tim	333
125	gezelter	507	int nframes = bsMan_->getNFrames();
126			nTimeBins_ = nframes;
127			histogram_.resize(nTimeBins_);
128			count_.resize(nTimeBins_);
129			time_.resize(nTimeBins_);
130	tim	333
131	gezelter	507	for (int i = 0; i < nTimeBins_; ++i) {
132	tim	333	time_[i] = i * deltaTime_;
133	gezelter	507	}
134
135	tim	333	}
136
137
138	gezelter	507	void TimeCorrFunc::doCorrelate() {
139	tim	333	preCorrelate();
140
141			int nblocks = bsMan_->getNBlocks();
142
143			for (int i = 0; i < nblocks; ++i) {
144	gezelter	507	bsMan_->loadBlock(i);
145	tim	351
146	gezelter	507	for (int j = i; j < nblocks; ++j) {
147			bsMan_->loadBlock(j);
148			correlateBlocks(i, j);
149			bsMan_->unloadBlock(j);
150			}
151	tim	351
152	gezelter	507	bsMan_->unloadBlock(i);
153	tim	333	}
154
155			postCorrelate();
156
157			writeCorrelate();
158	gezelter	507	}
159	tim	333
160	gezelter	507	void TimeCorrFunc::correlateBlocks(int block1, int block2) {
161	tim	333
162	tim	334	int jstart, jend;
163	tim	333
164	tim	334	assert(bsMan_->isBlockActive(block1) && bsMan_->isBlockActive(block2));
165	tim	333
166	tim	334	SnapshotBlock snapshotBlock1 = bsMan_->getSnapshotBlock(block1);
167			SnapshotBlock snapshotBlock2 = bsMan_->getSnapshotBlock(block2);
168	tim	333
169	tim	334	jend = snapshotBlock2.second;
170
171			for (int i = snapshotBlock1.first; i < snapshotBlock1.second; ++i) {
172
173	gezelter	507	//update the position or velocity of the atoms belong to rigid bodies
174			updateFrame(i);
175	tim	334
176	gezelter	507	// if the two blocks are the same, we don't want to correlate
177			// backwards in time, so start j at the same frame as i
178			if (block1 == block2) {
179	tim	334	jstart = i;
180	gezelter	507	} else {
181			jstart = snapshotBlock2.first;
182			}
183	tim	334
184	gezelter	507	for(int j = jstart; j < jend; ++j) {
185			//update the position or velocity of the atoms belong to rigid bodies
186			updateFrame(j);
187	tim	334
188	gezelter	507	correlateFrames(i, j);
189			}
190	tim	333	}
191	gezelter	507	}
192	tim	333
193	gezelter	507	void TimeCorrFunc::updateFrame(int frame){
194	tim	333	Molecule* mol;
195			RigidBody* rb;
196			SimInfo::MoleculeIterator mi;
197			Molecule::RigidBodyIterator rbIter;
198
199			/** @todo need improvement */
200			if (storageLayout_ & DataStorage::dslPosition) {
201	xsun	1183	for (mol = info_->beginMolecule(mi); mol != NULL;
202			mol = info_->nextMolecule(mi)) {
203	tim	333
204	gezelter	507	//change the positions of atoms which belong to the rigidbodies
205	xsun	1183	for (rb = mol->beginRigidBody(rbIter); rb != NULL;
206			rb = mol->nextRigidBody(rbIter)) {
207	gezelter	507	rb->updateAtoms(frame);
208			}
209			}
210	tim	333	}
211
212			if (storageLayout_ & DataStorage::dslVelocity) {
213	xsun	1183	for (mol = info_->beginMolecule(mi); mol != NULL;
214			mol = info_->nextMolecule(mi)) {
215
216	gezelter	507	//change the positions of atoms which belong to the rigidbodies
217	xsun	1183	for (rb = mol->beginRigidBody(rbIter); rb != NULL;
218			rb = mol->nextRigidBody(rbIter)) {
219	gezelter	507	rb->updateAtomVel(frame);
220			}
221			}
222	xsun	1183
223	tim	333	}
224
225	gezelter	507	}
226	tim	333
227
228	gezelter	507	void TimeCorrFunc::preCorrelate() {
229	tim	333	std::fill(histogram_.begin(), histogram_.end(), 0.0);
230			std::fill(count_.begin(), count_.end(), 0);
231	gezelter	507	}
232	tim	333
233	gezelter	507	void TimeCorrFunc::postCorrelate() {
234	tim	333	for (int i =0 ; i < nTimeBins_; ++i) {
235	gezelter	507	if (count_[i] > 0) {
236			histogram_[i] /= count_[i];
237			}
238	tim	333	}
239	gezelter	507	}
240	tim	333
241
242	gezelter	507	void TimeCorrFunc::writeCorrelate() {
243	tim	333	std::ofstream ofs(outputFilename_.c_str());
244
245			if (ofs.is_open()) {
246
247	gezelter	507	ofs << "#" << getCorrFuncType() << "\n";
248			ofs << "#selection script1: \"" << selectionScript1_ <<"\"\tselection script2: \"" << selectionScript2_ << "\"\n";
249			ofs << "#extra information: " << extra_ << "\n";
250			ofs << "#time\tcorrVal\n";
251	tim	333
252	gezelter	507	for (int i = 0; i < nTimeBins_; ++i) {
253			ofs << time_[i] << "\t" << histogram_[i] << "\n";
254			}
255	tim	333
256			} else {
257	gezelter	507	sprintf(painCave.errMsg,
258			"TimeCorrFunc::writeCorrelate Error: fail to open %s\n", outputFilename_.c_str());
259			painCave.isFatal = 1;
260			simError();
261	tim	333	}
262
263			ofs.close();
264	gezelter	507	}
265	tim	333
266			}