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(int 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(int 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(int 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(int 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(int 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(int 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(int 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:	2032
Committed:	Fri Oct 31 18:57:19 2014 UTC (10 years, 6 months ago) by jmichalk
File size:	9649 byte(s)
Log Message:	Forgot to svn add the files. Additionally, examined the files where changes were committed that I didn't realize were going to be changed. Primary one to examine SimplePreprocessor, bufferSize=8192 instead of 1024... everything else was commented cerr/cout or line spacing issues
#	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			//positions_ and moBool_ are 2D arrays, need the second dimension filled as well
100			for(int i = 0; i < selectionCount_; i++){
101			moBool_[i].resize(nFrames);
102			positions_[i].resize(nFrames);
103			}
104
105			int iterator;
106			int index = 0;
107			/* Loop over all frames storing the positions in a vec< vec<pos> >
108			* At the end, positions.length() should equal seleMan1_.size() or w/e
109			* And positions[index].length() should equal nFrames (or nFrames/istep)
110			*/
111			for(int istep = 0; istep < nFrames; istep += step_){
112			frames_++;
113			reader.readFrame(istep);
114			currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
115
116			for(mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)){
117			//change the positions of atoms which belong to the rigidbodies
118			for(rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)){
119			rb->updateAtoms();
120			}
121			}
122
123			index = 0; //count over atoms since iterators aren't the most friendly for such plebian things
124			for(sd = seleMan1_.beginSelected(iterator); sd != NULL; sd = seleMan1_.nextSelected(iterator)){
125			Vector3d pos = sd->getPos();
126			positions_[index][istep] = pos;
127			index++;
128			}
129			}
130
131
132			cout << "Position Array size: " << positions_.size() << "\n";
133			cout << "Frames analyzed: " << positions_[0].size() << "\n";
134
135
136			for(int i = 0; i < positions_.size(); i++){
137			int frameIndex = positions_[i].size();
138			for(int j = 1; j < frameIndex; j++){
139			Vector3d posF1 = positions_[i][j-1];
140			Vector3d posF2 = positions_[i][j];
141			Vector3d diff = posF2 - posF1;
142			if(usePeriodicBoundaryConditions_){
143			currentSnapshot_->wrapVector(diff);
144			}
145			double dist = diff.length();
146			if(dist > singleMoveDistance_){
147			moBool_[i][j] = true;
148			}else{
149			moBool_[i][j] = false;
150			}
151			}
152			}
153
154			int mobileAtomCount = 0;
155			for(int i = 0; i < moBool_.size(); i++){
156			int frameIndex = moBool_[i].size();
157			bool mobileAtom = false;
158			for(int j = 0; j < frameIndex; j++){
159			mobileAtom = mobileAtom \|\| moBool_[i][j];
160			}
161			moBool_[i][0] = mobileAtom; //is true if any value later in the array is true, false otherwise
162			if(mobileAtom){
163			mobileAtomCount++;
164			}
165			}
166
167			cout << "Mobile atom count: " << mobileAtomCount << "\n";
168
169			//Here I shrink the size of the arrays, why look through 3888, when you only need ~800.
170			//Additionally, all of these are mobile at some point in time, the others aren't, dead weight and
171			//memory
172			positions2_.resize(mobileAtomCount);
173			moBool2_.resize(mobileAtomCount);
174			int pos2index = 0;
175			for(int i = 0; i < positions_.size(); i++){
176			int frameCount = positions_[i].size();
177			if(moBool_[i][0]){
178			for(int j = 0; j < frameCount; j++){
179			positions2_[pos2index].push_back(positions_[i][j]);
180			moBool2_[pos2index].push_back(moBool_[i][j]);
181			}
182			pos2index++;
183			}
184			}
185
186			positions_.clear();
187			moBool_.clear();
188
189			cout << "positions_ has been cleared: " << positions_.size() << "\n";
190			cout << "positions2_ has been filled: " << positions2_.size() << "\n";
191			cout << "positions2_ has " << positions2_[0].size() << " frames\n";
192
193			//The important one!
194			positionCorrelation();
195
196
197			//Write out my data
198			std::ofstream diffStream;
199			setOutputName(getPrefix(filename_) + ".Mdiffusion");
200			diffStream.open(outputFilename_.c_str());
201			diffStream << "#X&Y diffusion amounts\n";
202			diffStream << "#singleMoveDistance_: " << singleMoveDistance_ << "\n";
203
204			diffStream << "#Number of mobile atoms: " << positions2_.size() << "\n";
205
206			diffStream << "#time, <x(t)-x(0)>, <y(t)-y(0)>, <r(t)-r(0)>\n";
207			for(int i = 0; i < xHist_.size(); i++){
208			diffStream << i << ", " << xHist_[i] << ", " << yHist_[i] << ", " << rHist_[i] << "\n";
209			}
210			diffStream.close();
211
212			}
213
214			void SurfaceDiffusion::positionCorrelation(){
215			RealType xDist = 0.0;
216			RealType yDist = 0.0;
217			RealType rDist = 0.0;
218			int timeShift = 0;
219			Vector3d kPos;
220			Vector3d jPos;
221			//biggest timeShift is positions2_[0].size() - 1?
222			xHist_.clear();
223			yHist_.clear();
224			rHist_.clear();
225			count_.clear();
226			int frameResize = positions2_[0].size();
227			xHist_.resize(frameResize);
228			yHist_.resize(frameResize);
229			rHist_.resize(frameResize);
230			count_.resize(frameResize);
231			//loop over particles
232			// loop over frames starting at j
233			// loop over frames starting at k = j (time shift of 0)
234			for(int i = 0; i < positions2_.size(); i++){
235			int frames = positions2_[i].size() - 1; //for counting properly, otherwise moBool2_[i][j+1] will
236			// go over
237			for(int j = 0; j < frames; j++){
238			//if the particle is mobile between j and j + 1, then count it for all timeShifts
239			if(moBool2_[i][j+1]){
240			for(int k = j; k < positions2_[0].size(); k++){
241			//<x(t)-x(0)> <y(t)-y(0)> <r(t)-r(0)>
242			//The positions stored are not wrapped, thus I don't need to worry about pbc
243			//Mean square displacement
244			//So I do want the squared distances
245
246			kPos = positions2_[i][k];
247			jPos = positions2_[i][j];
248			xDist = kPos.x() - jPos.x();
249			xDist = xDist*xDist;
250
251			yDist = kPos.y() - jPos.y();
252			yDist = yDist*yDist;
253
254			rDist = (kPos - jPos).lengthSquare();
255
256
257			timeShift = k - j;
258			xHist_[timeShift] += xDist;
259			yHist_[timeShift] += yDist;
260			rHist_[timeShift] += rDist;
261			count_[timeShift]++;
262			}
263			}
264			}
265			}
266			cout << "X, Y, R calculated\n";
267
268			for(int i = 0; i < xHist_.size(); i++){
269			xHist_[i] = xHist_[i]/(count_[i]);
270			yHist_[i] = yHist_[i]/(count_[i]);
271			rHist_[i] = rHist_[i]/(count_[i]);
272			}
273			cout << "X, Y, R normalized\n";
274			}
275
276			}