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.
#	Content
1	/*
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
100	// filled as well
101	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	* 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	*/
113	for(int istep = 0; istep < nFrames; istep += step_){
114	frames_++;
115	reader.readFrame(istep);
116	currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
117
118	for(mol = info_->beginMolecule(mi); mol != NULL;
119	mol = info_->nextMolecule(mi)){
120	//change the positions of atoms which belong to the rigidbodies
121	for(rb = mol->beginRigidBody(rbIter); rb != NULL;
122	rb = mol->nextRigidBody(rbIter)){
123	rb->updateAtoms();
124	}
125	}
126
127	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	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	for(std::size_t i = 0; i < positions_.size(); i++){
141	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	for(std::size_t i = 0; i < moBool_.size(); i++){
160	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	moBool_[i][0] = mobileAtom; // is true if any value later in the
166	// array is true, false otherwise
167	if(mobileAtom){
168	mobileAtomCount++;
169	}
170	}
171
172	cout << "Mobile atom count: " << mobileAtomCount << "\n";
173
174	// 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	positions2_.resize(mobileAtomCount);
179	moBool2_.resize(mobileAtomCount);
180	int pos2index = 0;
181	for(std::size_t i = 0; i < positions_.size(); i++){
182	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	diffStream << "#time, <x(t)-x(0)>, <y(t)-y(0)>, <r(t)-r(0)>\n";
211
212	for(std::size_t i = 0; i < xHist_.size(); i++){
213	diffStream << i << ", " << xHist_[i] << ", " << yHist_[i] << ", "
214	<< rHist_[i] << "\n";
215	}
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	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	for(int j = 0; j < frames; j++){
246	// if the particle is mobile between j and j + 1, then count
247	// it for all timeShifts
248	if(moBool2_[i][j+1]){
249	for(std::size_t k = j; k < positions2_[0].size(); k++){
250	//<x(t)-x(0)> <y(t)-y(0)> <r(t)-r(0)>
251	//The positions stored are not wrapped, thus I don't need
252	//to worry about pbc
253	//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	for(std::size_t i = 0; i < xHist_.size(); i++){
279	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	}