applications/staticProps/SurfaceDiffusion.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] Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
 * [4] , Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). *
 *  Created by Joseph R. Michalka on Oct 12 2012 
 *  @author  Joseph R. Michalka 
 *  @version $Id: RhoZ.cpp 1665 2011-11-22 20:38:56Z gezelter $
 *
 */

/* Surface Diffusion 
 * Attempting to track/measure the surface diffusion rates of particles on... wait for it..
 * a surface.
 * This program was initially created to track Platinum particles moving around a 557 surface.
 * Hence why we are trying to keep the x and y movement separate. 
 *
 */

#include <algorithm>
#include <fstream>
#include "applications/staticProps/SurfaceDiffusion.hpp"
#include "utils/simError.h"
#include "io/DumpReader.hpp"
#include "primitives/Molecule.hpp"
namespace OpenMD {
  
  SurfaceDiffusion::SurfaceDiffusion(SimInfo* info, const std::string& filename, const std::string& sele, RealType len)
    : StaticAnalyser(info, filename), selectionScript_(sele),  evaluator_(info), seleMan1_(info){

    evaluator_.loadScriptString(sele);
    if (!evaluator_.isDynamic()) {
      seleMan1_.setSelectionSet(evaluator_.evaluate());
    }       
    
    //Depending on the selection 'sele1="select Pt"' need a vector equal to the
    //number of Platinums in the system (for this specific case)
    selectionCount_ = seleMan1_.getSelectionCount();
    cout << "SelectionCount_: " << selectionCount_ << "\n";
    
    moBool_.resize(selectionCount_);
    positions_.resize(selectionCount_);

    filename_ = filename; 
    singleMoveDistance_ = 2.0;
  }

  SurfaceDiffusion::~SurfaceDiffusion(){
  
  }

  void SurfaceDiffusion::process() {
    Molecule* mol;
    RigidBody* rb;
    StuntDouble* sd;
    SimInfo::MoleculeIterator mi;
    Molecule::RigidBodyIterator rbIter;

    DumpReader reader(info_, dumpFilename_);    
    int nFrames = reader.getNFrames();
    frames_ = 0;
    nProcessed_ = nFrames/step_;

    // positions_ and moBool_ are 2D arrays, need the second dimension
    // filled as well
    for(int i = 0; i < selectionCount_; i++){
      moBool_[i].resize(nFrames);
      positions_[i].resize(nFrames);
    }

    int iterator;
    int index = 0; 
    /* Loop over all frames storing the positions in a vec< vec<pos> >
     * At the end, positions.length() should equal seleMan1_.size() or
     * w/e And positions[index].length() should equal nFrames (or
     * nFrames/istep)
     */
    for(int istep = 0; istep < nFrames; istep += step_){
      frames_++;
      reader.readFrame(istep);
      currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();

      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();
        }
      }
      
      index = 0; // count over atoms since iterators aren't the most
                 // friendly for such plebian things
      for(sd = seleMan1_.beginSelected(iterator); sd != NULL; 
          sd = seleMan1_.nextSelected(iterator)){
        Vector3d pos = sd->getPos();
        positions_[index][istep] = pos;
        index++;
      }      
    }

    cout << "Position Array size: " << positions_.size() << "\n";
    cout << "Frames analyzed: " << positions_[0].size() << "\n";
      
    for(std::size_t i = 0; i < positions_.size(); i++){
      int frameIndex = positions_[i].size();
      for(int j = 1; j < frameIndex; j++){
        Vector3d posF1 = positions_[i][j-1];
        Vector3d posF2 = positions_[i][j];
        Vector3d diff = posF2 - posF1;
        if(usePeriodicBoundaryConditions_){
          currentSnapshot_->wrapVector(diff);
        }
        double dist = diff.length();
        if(dist > singleMoveDistance_){
          moBool_[i][j] = true;
        }else{
          moBool_[i][j] = false;
        }
      }
    } 

    int mobileAtomCount = 0;
    for(std::size_t i = 0; i < moBool_.size(); i++){
      int frameIndex = moBool_[i].size();
      bool mobileAtom = false;
      for(int j = 0; j < frameIndex; j++){
        mobileAtom = mobileAtom || moBool_[i][j];
      }
      moBool_[i][0] = mobileAtom; // is true if any value later in the
                                  // array is true, false otherwise
      if(mobileAtom){
        mobileAtomCount++;
      }
    } 

    cout << "Mobile atom count: " << mobileAtomCount << "\n";

    // Here I shrink the size of the arrays, why look through 3888,
    // when you only need ~800.  Additionally, all of these are mobile
    // at some point in time, the others aren't, dead weight and
    // memory
    positions2_.resize(mobileAtomCount);
    moBool2_.resize(mobileAtomCount);
    int pos2index = 0;
    for(std::size_t i = 0; i < positions_.size(); i++){
      int frameCount = positions_[i].size();
      if(moBool_[i][0]){
        for(int j = 0; j < frameCount; j++){
          positions2_[pos2index].push_back(positions_[i][j]);
          moBool2_[pos2index].push_back(moBool_[i][j]);
        }
        pos2index++;
      }
    }   
    
    positions_.clear();
    moBool_.clear();

    cout << "positions_ has been cleared: " << positions_.size() << "\n";
    cout << "positions2_ has been filled: " << positions2_.size() << "\n";
    cout << "positions2_ has " << positions2_[0].size() << " frames\n";
  
    //The important one! 
    positionCorrelation();
 

    //Write out my data
    std::ofstream diffStream;
    setOutputName(getPrefix(filename_) + ".Mdiffusion");
    diffStream.open(outputFilename_.c_str());
    diffStream << "#X&Y diffusion amounts\n";
    diffStream << "#singleMoveDistance_: " << singleMoveDistance_ << "\n";
    diffStream << "#Number of mobile atoms: " << positions2_.size() << "\n";
    diffStream << "#time, <x(t)-x(0)>, <y(t)-y(0)>, <r(t)-r(0)>\n";

    for(std::size_t i = 0; i < xHist_.size(); i++){
      diffStream << i << ", " << xHist_[i] << ", " << yHist_[i] << ", " 
                 << rHist_[i] << "\n";
    }
    diffStream.close();

  }
 
  void SurfaceDiffusion::positionCorrelation(){
    RealType xDist = 0.0;
    RealType yDist = 0.0;
    RealType rDist = 0.0;
    int timeShift = 0;
    Vector3d kPos;
    Vector3d jPos;
    //biggest timeShift is positions2_[0].size() - 1?
    xHist_.clear();
    yHist_.clear();
    rHist_.clear();
    count_.clear();
    int frameResize = positions2_[0].size();
    xHist_.resize(frameResize);
    yHist_.resize(frameResize);
    rHist_.resize(frameResize);
    count_.resize(frameResize);
    //loop over particles
    // loop over frames starting at j
    //  loop over frames starting at k = j (time shift of 0)
    for(std::size_t i = 0; i < positions2_.size(); i++){
      int frames = positions2_[i].size() - 1; // for counting
                                              // properly, otherwise
                                              // moBool2_[i][j+1] will
                                              // go over
      for(int j = 0; j < frames; j++){
        // if the particle is mobile between j and j + 1, then count
        // it for all timeShifts
        if(moBool2_[i][j+1]){
          for(std::size_t k = j; k < positions2_[0].size(); k++){
            //<x(t)-x(0)>  <y(t)-y(0)>  <r(t)-r(0)>
            //The positions stored are not wrapped, thus I don't need
            //to worry about pbc          
            //Mean square displacement
            //So I do want the squared distances
          
            kPos = positions2_[i][k];
            jPos = positions2_[i][j];
            xDist = kPos.x() - jPos.x();
            xDist = xDist*xDist;
          
            yDist = kPos.y() - jPos.y();
            yDist = yDist*yDist;

            rDist = (kPos - jPos).lengthSquare();
  

            timeShift = k - j;
            xHist_[timeShift] += xDist;
            yHist_[timeShift] += yDist;  
            rHist_[timeShift] += rDist;
            count_[timeShift]++;
          }
        }
      }
    }
    cout << "X, Y, R calculated\n";

    for(std::size_t i = 0; i < xHist_.size(); i++){
      xHist_[i] = xHist_[i]/(count_[i]);
      yHist_[i] = yHist_[i]/(count_[i]);
      rHist_[i] = rHist_[i]/(count_[i]);
    }
    cout << "X, Y, R normalized\n";
  }

}
Revision:	2071
Committed:	Sat Mar 7 21:41:51 2015 UTC (10 years, 7 months ago) by gezelter
File size:	9806 byte(s)
Log Message:	Reducing the number of warnings when using g++ to compile.
#	User	Rev	Content
1	jmichalk	2032	/*
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, 24107 (2008).
39			* [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010).
40			* [4] , Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). *
41			* Created by Joseph R. Michalka on Oct 12 2012
42			* @author Joseph R. Michalka
43			* @version $Id: RhoZ.cpp 1665 2011-11-22 20:38:56Z gezelter $
44			*
45			*/
46
47			/* Surface Diffusion
48			* Attempting to track/measure the surface diffusion rates of particles on... wait for it..
49			* a surface.
50			* This program was initially created to track Platinum particles moving around a 557 surface.
51			* Hence why we are trying to keep the x and y movement separate.
52			*
53			*/
54
55			#include <algorithm>
56			#include <fstream>
57			#include "applications/staticProps/SurfaceDiffusion.hpp"
58			#include "utils/simError.h"
59			#include "io/DumpReader.hpp"
60			#include "primitives/Molecule.hpp"
61			namespace OpenMD {
62
63			SurfaceDiffusion::SurfaceDiffusion(SimInfo* info, const std::string& filename, const std::string& sele, RealType len)
64			: StaticAnalyser(info, filename), selectionScript_(sele), evaluator_(info), seleMan1_(info){
65
66			evaluator_.loadScriptString(sele);
67			if (!evaluator_.isDynamic()) {
68			seleMan1_.setSelectionSet(evaluator_.evaluate());
69			}
70
71			//Depending on the selection 'sele1="select Pt"' need a vector equal to the
72			//number of Platinums in the system (for this specific case)
73			selectionCount_ = seleMan1_.getSelectionCount();
74			cout << "SelectionCount_: " << selectionCount_ << "\n";
75
76			moBool_.resize(selectionCount_);
77			positions_.resize(selectionCount_);
78
79			filename_ = filename;
80			singleMoveDistance_ = 2.0;
81			}
82
83			SurfaceDiffusion::~SurfaceDiffusion(){
84
85			}
86
87			void SurfaceDiffusion::process() {
88			Molecule* mol;
89			RigidBody* rb;
90			StuntDouble* sd;
91			SimInfo::MoleculeIterator mi;
92			Molecule::RigidBodyIterator rbIter;
93
94			DumpReader reader(info_, dumpFilename_);
95			int nFrames = reader.getNFrames();
96			frames_ = 0;
97			nProcessed_ = nFrames/step_;
98
99	gezelter	2071	// positions_ and moBool_ are 2D arrays, need the second dimension
100			// filled as well
101	jmichalk	2032	for(int i = 0; i < selectionCount_; i++){
102			moBool_[i].resize(nFrames);
103			positions_[i].resize(nFrames);
104			}
105
106			int iterator;
107			int index = 0;
108			/* Loop over all frames storing the positions in a vec< vec<pos> >
109	gezelter	2071	* At the end, positions.length() should equal seleMan1_.size() or
110			* w/e And positions[index].length() should equal nFrames (or
111			* nFrames/istep)
112	jmichalk	2032	*/
113			for(int istep = 0; istep < nFrames; istep += step_){
114			frames_++;
115			reader.readFrame(istep);
116			currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
117
118	gezelter	2071	for(mol = info_->beginMolecule(mi); mol != NULL;
119			mol = info_->nextMolecule(mi)){
120	jmichalk	2032	//change the positions of atoms which belong to the rigidbodies
121	gezelter	2071	for(rb = mol->beginRigidBody(rbIter); rb != NULL;
122			rb = mol->nextRigidBody(rbIter)){
123	jmichalk	2032	rb->updateAtoms();
124			}
125			}
126
127	gezelter	2071	index = 0; // count over atoms since iterators aren't the most
128			// friendly for such plebian things
129			for(sd = seleMan1_.beginSelected(iterator); sd != NULL;
130			sd = seleMan1_.nextSelected(iterator)){
131	jmichalk	2032	Vector3d pos = sd->getPos();
132			positions_[index][istep] = pos;
133			index++;
134			}
135			}
136
137			cout << "Position Array size: " << positions_.size() << "\n";
138			cout << "Frames analyzed: " << positions_[0].size() << "\n";
139
140	gezelter	2071	for(std::size_t i = 0; i < positions_.size(); i++){
141	jmichalk	2032	int frameIndex = positions_[i].size();
142			for(int j = 1; j < frameIndex; j++){
143			Vector3d posF1 = positions_[i][j-1];
144			Vector3d posF2 = positions_[i][j];
145			Vector3d diff = posF2 - posF1;
146			if(usePeriodicBoundaryConditions_){
147			currentSnapshot_->wrapVector(diff);
148			}
149			double dist = diff.length();
150			if(dist > singleMoveDistance_){
151			moBool_[i][j] = true;
152			}else{
153			moBool_[i][j] = false;
154			}
155			}
156			}
157
158			int mobileAtomCount = 0;
159	gezelter	2071	for(std::size_t i = 0; i < moBool_.size(); i++){
160	jmichalk	2032	int frameIndex = moBool_[i].size();
161			bool mobileAtom = false;
162			for(int j = 0; j < frameIndex; j++){
163			mobileAtom = mobileAtom \|\| moBool_[i][j];
164			}
165	gezelter	2071	moBool_[i][0] = mobileAtom; // is true if any value later in the
166			// array is true, false otherwise
167	jmichalk	2032	if(mobileAtom){
168			mobileAtomCount++;
169			}
170			}
171
172			cout << "Mobile atom count: " << mobileAtomCount << "\n";
173
174	gezelter	2071	// Here I shrink the size of the arrays, why look through 3888,
175			// when you only need ~800. Additionally, all of these are mobile
176			// at some point in time, the others aren't, dead weight and
177			// memory
178	jmichalk	2032	positions2_.resize(mobileAtomCount);
179			moBool2_.resize(mobileAtomCount);
180			int pos2index = 0;
181	gezelter	2071	for(std::size_t i = 0; i < positions_.size(); i++){
182	jmichalk	2032	int frameCount = positions_[i].size();
183			if(moBool_[i][0]){
184			for(int j = 0; j < frameCount; j++){
185			positions2_[pos2index].push_back(positions_[i][j]);
186			moBool2_[pos2index].push_back(moBool_[i][j]);
187			}
188			pos2index++;
189			}
190			}
191
192			positions_.clear();
193			moBool_.clear();
194
195			cout << "positions_ has been cleared: " << positions_.size() << "\n";
196			cout << "positions2_ has been filled: " << positions2_.size() << "\n";
197			cout << "positions2_ has " << positions2_[0].size() << " frames\n";
198
199			//The important one!
200			positionCorrelation();
201
202
203			//Write out my data
204			std::ofstream diffStream;
205			setOutputName(getPrefix(filename_) + ".Mdiffusion");
206			diffStream.open(outputFilename_.c_str());
207			diffStream << "#X&Y diffusion amounts\n";
208			diffStream << "#singleMoveDistance_: " << singleMoveDistance_ << "\n";
209			diffStream << "#Number of mobile atoms: " << positions2_.size() << "\n";
210	gezelter	2071	diffStream << "#time, <x(t)-x(0)>, <y(t)-y(0)>, <r(t)-r(0)>\n";
211	jmichalk	2032
212	gezelter	2071	for(std::size_t i = 0; i < xHist_.size(); i++){
213			diffStream << i << ", " << xHist_[i] << ", " << yHist_[i] << ", "
214			<< rHist_[i] << "\n";
215	jmichalk	2032	}
216			diffStream.close();
217
218			}
219
220			void SurfaceDiffusion::positionCorrelation(){
221			RealType xDist = 0.0;
222			RealType yDist = 0.0;
223			RealType rDist = 0.0;
224			int timeShift = 0;
225			Vector3d kPos;
226			Vector3d jPos;
227			//biggest timeShift is positions2_[0].size() - 1?
228			xHist_.clear();
229			yHist_.clear();
230			rHist_.clear();
231			count_.clear();
232			int frameResize = positions2_[0].size();
233			xHist_.resize(frameResize);
234			yHist_.resize(frameResize);
235			rHist_.resize(frameResize);
236			count_.resize(frameResize);
237			//loop over particles
238			// loop over frames starting at j
239			// loop over frames starting at k = j (time shift of 0)
240	gezelter	2071	for(std::size_t i = 0; i < positions2_.size(); i++){
241			int frames = positions2_[i].size() - 1; // for counting
242			// properly, otherwise
243			// moBool2_[i][j+1] will
244			// go over
245	jmichalk	2032	for(int j = 0; j < frames; j++){
246	gezelter	2071	// if the particle is mobile between j and j + 1, then count
247			// it for all timeShifts
248	jmichalk	2032	if(moBool2_[i][j+1]){
249	gezelter	2071	for(std::size_t k = j; k < positions2_[0].size(); k++){
250	jmichalk	2032	//<x(t)-x(0)> <y(t)-y(0)> <r(t)-r(0)>
251	gezelter	2071	//The positions stored are not wrapped, thus I don't need
252			//to worry about pbc
253	jmichalk	2032	//Mean square displacement
254			//So I do want the squared distances
255
256			kPos = positions2_[i][k];
257			jPos = positions2_[i][j];
258			xDist = kPos.x() - jPos.x();
259			xDist = xDist*xDist;
260
261			yDist = kPos.y() - jPos.y();
262			yDist = yDist*yDist;
263
264			rDist = (kPos - jPos).lengthSquare();
265
266
267			timeShift = k - j;
268			xHist_[timeShift] += xDist;
269			yHist_[timeShift] += yDist;
270			rHist_[timeShift] += rDist;
271			count_[timeShift]++;
272			}
273			}
274			}
275			}
276			cout << "X, Y, R calculated\n";
277
278	gezelter	2071	for(std::size_t i = 0; i < xHist_.size(); i++){
279	jmichalk	2032	xHist_[i] = xHist_[i]/(count_[i]);
280			yHist_[i] = yHist_[i]/(count_[i]);
281			rHist_[i] = rHist_[i]/(count_[i]);
282			}
283			cout << "X, Y, R normalized\n";
284			}
285
286			}