applications/staticProps/RadialDistrFunc.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).
 * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
 */

#include <algorithm>

#include "RadialDistrFunc.hpp"
#include "io/DumpReader.hpp"
#include "primitives/Molecule.hpp"

namespace OpenMD {

  RadialDistrFunc::RadialDistrFunc(SimInfo* info, 
                                   const std::string& filename, 
                                   const std::string& sele1, 
                                   const std::string& sele2)
    : StaticAnalyser(info, filename), selectionScript1_(sele1), 
      selectionScript2_(sele2), evaluator1_(info), evaluator2_(info), 
      seleMan1_(info), seleMan2_(info), common_(info), 
      sele1_minus_common_(info), sele2_minus_common_(info) {
          
      evaluator1_.loadScriptString(sele1);
      evaluator2_.loadScriptString(sele2);

      if (!evaluator1_.isDynamic()) {
        seleMan1_.setSelectionSet(evaluator1_.evaluate());
        validateSelection1(seleMan1_);
      }
      if (!evaluator2_.isDynamic()) {
        seleMan2_.setSelectionSet(evaluator2_.evaluate());
        validateSelection2(seleMan2_);
      }

      if (!evaluator1_.isDynamic() && !evaluator2_.isDynamic()) {
        // If all selections are static, we can precompute the number
        // of real pairs.
        common_ = seleMan1_ & seleMan2_;
        sele1_minus_common_ = seleMan1_ - common_;
        sele2_minus_common_ = seleMan2_ - common_;      

        int nSelected1 = seleMan1_.getSelectionCount();
        int nSelected2 = seleMan2_.getSelectionCount();
        int nIntersect = common_.getSelectionCount();
        
        nPairs_ = nSelected1 * nSelected2 - (nIntersect +1) * nIntersect/2;  
      }
    
    }

  void RadialDistrFunc::process() {
    Molecule* mol;
    RigidBody* rb;
    SimInfo::MoleculeIterator mi;
    Molecule::RigidBodyIterator rbIter;
    
    preProcess();
    
    DumpReader reader(info_, dumpFilename_);    
    int nFrames = reader.getNFrames();
    nProcessed_ = nFrames / step_;

    for (int i = 0; i < nFrames; i += step_) {
      reader.readFrame(i);
      currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();

      if (evaluator1_.isDynamic()) {
        seleMan1_.setSelectionSet(evaluator1_.evaluate());
        validateSelection1(seleMan1_);
      }
      if (evaluator2_.isDynamic()) {
        seleMan2_.setSelectionSet(evaluator2_.evaluate());
        validateSelection2(seleMan2_);
      }

      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();
        }
      }
        
      initializeHistogram();
        
      // Selections may overlap, and we need a bit of logic to deal
      // with this.
      //
      // |     s1    |
      // | s1 -c | c |
      //         | c | s2 - c |
      //         |    s2      |
      //
      // s1 : Set of StuntDoubles in selection1
      // s2 : Set of StuntDoubles in selection2
      // c  : Intersection of selection1 and selection2
      // 
      // When we loop over the pairs, we can divide the looping into 3
      // stages:
      //
      // Stage 1 :     [s1-c]      [s2]
      // Stage 2 :     [c]         [s2 - c]
      // Stage 3 :     [c]         [c]
      // Stages 1 and 2 are completely non-overlapping.
      // Stage 3 is completely overlapping.

      if (evaluator1_.isDynamic() || evaluator2_.isDynamic()) {
        common_ = seleMan1_ & seleMan2_;
        sele1_minus_common_ = seleMan1_ - common_;
        sele2_minus_common_ = seleMan2_ - common_;            
        int nSelected1 = seleMan1_.getSelectionCount();
        int nSelected2 = seleMan2_.getSelectionCount();
        int nIntersect = common_.getSelectionCount();
            
        nPairs_ = nSelected1 * nSelected2 - (nIntersect +1) * nIntersect/2;
      }
      
      processNonOverlapping(sele1_minus_common_, seleMan2_);
      processNonOverlapping(common_,             sele2_minus_common_);
      processOverlapping(common_);
      
      processHistogram();
        
    }

    postProcess();

    writeRdf();
  }

  void RadialDistrFunc::processNonOverlapping( SelectionManager& sman1, 
                                               SelectionManager& sman2) {
    StuntDouble* sd1;
    StuntDouble* sd2;
    int i;    
    int j;
    
    // This is the same as a non-overlapping pairwise loop structure:
    // for (int i = 0;  i < ni ; ++i ) {
    //   for (int j = 0; j < nj; ++j) {} 
    // }

    for (sd1 = sman1.beginSelected(i); sd1 != NULL; 
         sd1 = sman1.nextSelected(i)) {
      for (sd2 = sman2.beginSelected(j); sd2 != NULL; 
           sd2 = sman2.nextSelected(j)) {
        collectHistogram(sd1, sd2);
      }
    }
  }

  void RadialDistrFunc::processOverlapping( SelectionManager& sman) {
    StuntDouble* sd1;
    StuntDouble* sd2;
    int i;    
    int j;

    // This is the same as a pairwise loop structure:
    // for (int i = 0;  i < n-1 ; ++i ) {
    //   for (int j = i + 1; j < n; ++j) {} 
    // }
    
    for (sd1 = sman.beginSelected(i); sd1 != NULL; 
         sd1 = sman.nextSelected(i)) {                    
      for (j  = i, sd2 = sman.nextSelected(j); sd2 != NULL; 
           sd2 = sman.nextSelected(j)) {
        collectHistogram(sd1, sd2);
      }            
    }
  }
}
Revision:	1785
Committed:	Wed Aug 22 18:43:27 2012 UTC (13 years, 2 months ago) by jmichalk
File size:	7432 byte(s)
Log Message:	Trunk: The changes in this commit are confined to applications/staticProps and for the most part deal with a misspelling of initialize. The one other change took place in StaticProps.cpp and deals with the default treatment of sele2. It had previously been set to 'select all' which seems to go against what would be desired by not specifying it with regard to proper operations of many of the analysis programs ( g of r's especially)
#	User	Rev	Content
1	tim	306	/*
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	306	* 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	306	* 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	gezelter	1782	* [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010).
40			* [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41	tim	306	*/
42
43			#include <algorithm>
44
45			#include "RadialDistrFunc.hpp"
46	tim	311	#include "io/DumpReader.hpp"
47			#include "primitives/Molecule.hpp"
48	gezelter	1782
49	gezelter	1390	namespace OpenMD {
50	tim	306
51	gezelter	1782	RadialDistrFunc::RadialDistrFunc(SimInfo* info,
52			const std::string& filename,
53			const std::string& sele1,
54			const std::string& sele2)
55			: StaticAnalyser(info, filename), selectionScript1_(sele1),
56			selectionScript2_(sele2), evaluator1_(info), evaluator2_(info),
57			seleMan1_(info), seleMan2_(info), common_(info),
58			sele1_minus_common_(info), sele2_minus_common_(info) {
59	tim	306
60	gezelter	507	evaluator1_.loadScriptString(sele1);
61			evaluator2_.loadScriptString(sele2);
62	tim	306
63	gezelter	507	if (!evaluator1_.isDynamic()) {
64	tim	361	seleMan1_.setSelectionSet(evaluator1_.evaluate());
65			validateSelection1(seleMan1_);
66	gezelter	507	}
67			if (!evaluator2_.isDynamic()) {
68	tim	361	seleMan2_.setSelectionSet(evaluator2_.evaluate());
69			validateSelection2(seleMan2_);
70	gezelter	507	}
71	tim	306
72	gezelter	507	if (!evaluator1_.isDynamic() && !evaluator2_.isDynamic()) {
73	gezelter	1782	// If all selections are static, we can precompute the number
74			// of real pairs.
75	tim	353	common_ = seleMan1_ & seleMan2_;
76			sele1_minus_common_ = seleMan1_ - common_;
77			sele2_minus_common_ = seleMan2_ - common_;
78	tim	347
79	tim	348	int nSelected1 = seleMan1_.getSelectionCount();
80			int nSelected2 = seleMan2_.getSelectionCount();
81	tim	353	int nIntersect = common_.getSelectionCount();
82
83	gezelter	1782	nPairs_ = nSelected1 * nSelected2 - (nIntersect +1) * nIntersect/2;
84	gezelter	507	}
85
86	tim	347	}
87	tim	306
88	gezelter	507	void RadialDistrFunc::process() {
89	tim	311	Molecule* mol;
90			RigidBody* rb;
91			SimInfo::MoleculeIterator mi;
92			Molecule::RigidBodyIterator rbIter;
93	tim	318
94	tim	306	preProcess();
95
96			DumpReader reader(info_, dumpFilename_);
97	tim	311	int nFrames = reader.getNFrames();
98	tim	351	nProcessed_ = nFrames / step_;
99
100	tim	306	for (int i = 0; i < nFrames; i += step_) {
101	gezelter	507	reader.readFrame(i);
102			currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
103	tim	306
104	gezelter	507	if (evaluator1_.isDynamic()) {
105			seleMan1_.setSelectionSet(evaluator1_.evaluate());
106			validateSelection1(seleMan1_);
107			}
108			if (evaluator2_.isDynamic()) {
109			seleMan2_.setSelectionSet(evaluator2_.evaluate());
110			validateSelection2(seleMan2_);
111			}
112	tim	306
113	gezelter	1782	for (mol = info_->beginMolecule(mi); mol != NULL;
114			mol = info_->nextMolecule(mi)) {
115	tim	311
116	gezelter	1782	// Change the positions of atoms which belong to the RigidBodies
117			for (rb = mol->beginRigidBody(rbIter); rb != NULL;
118			rb = mol->nextRigidBody(rbIter)) {
119	gezelter	507	rb->updateAtoms();
120			}
121			}
122	tim	311
123	jmichalk	1785	initializeHistogram();
124	tim	353
125	gezelter	1782	// Selections may overlap, and we need a bit of logic to deal
126			// with this.
127	gezelter	507	//
128	gezelter	1782	// \| s1 \|
129			// \| s1 -c \| c \|
130			// \| c \| s2 - c \|
131			// \| s2 \|
132	gezelter	507	//
133	gezelter	1782	// s1 : Set of StuntDoubles in selection1
134			// s2 : Set of StuntDoubles in selection2
135			// c : Intersection of selection1 and selection2
136			//
137			// When we loop over the pairs, we can divide the looping into 3
138			// stages:
139			//
140			// Stage 1 : [s1-c] [s2]
141			// Stage 2 : [c] [s2 - c]
142			// Stage 3 : [c] [c]
143			// Stages 1 and 2 are completely non-overlapping.
144			// Stage 3 is completely overlapping.
145	tim	353
146	gezelter	507	if (evaluator1_.isDynamic() \|\| evaluator2_.isDynamic()) {
147			common_ = seleMan1_ & seleMan2_;
148			sele1_minus_common_ = seleMan1_ - common_;
149			sele2_minus_common_ = seleMan2_ - common_;
150			int nSelected1 = seleMan1_.getSelectionCount();
151			int nSelected2 = seleMan2_.getSelectionCount();
152			int nIntersect = common_.getSelectionCount();
153	tim	353
154	xsun	1213	nPairs_ = nSelected1 * nSelected2 - (nIntersect +1) * nIntersect/2;
155	gezelter	507	}
156	gezelter	1782
157	gezelter	507	processNonOverlapping(sele1_minus_common_, seleMan2_);
158	gezelter	1782	processNonOverlapping(common_, sele2_minus_common_);
159	gezelter	507	processOverlapping(common_);
160	xsun	1213
161	gezelter	507	processHistogram();
162	tim	306
163			}
164
165			postProcess();
166
167	tim	307	writeRdf();
168	gezelter	507	}
169	tim	306
170	gezelter	1782	void RadialDistrFunc::processNonOverlapping( SelectionManager& sman1,
171			SelectionManager& sman2) {
172	tim	353	StuntDouble* sd1;
173			StuntDouble* sd2;
174			int i;
175			int j;
176
177	gezelter	1782	// This is the same as a non-overlapping pairwise loop structure:
178			// for (int i = 0; i < ni ; ++i ) {
179			// for (int j = 0; j < nj; ++j) {}
180			// }
181	tim	347
182	gezelter	1782	for (sd1 = sman1.beginSelected(i); sd1 != NULL;
183			sd1 = sman1.nextSelected(i)) {
184			for (sd2 = sman2.beginSelected(j); sd2 != NULL;
185			sd2 = sman2.nextSelected(j)) {
186	gezelter	507	collectHistogram(sd1, sd2);
187	gezelter	1782	}
188	tim	353	}
189	gezelter	507	}
190	tim	347
191	gezelter	507	void RadialDistrFunc::processOverlapping( SelectionManager& sman) {
192	tim	353	StuntDouble* sd1;
193			StuntDouble* sd2;
194			int i;
195			int j;
196
197	gezelter	1782	// This is the same as a pairwise loop structure:
198			// for (int i = 0; i < n-1 ; ++i ) {
199			// for (int j = i + 1; j < n; ++j) {}
200			// }
201	tim	353
202	gezelter	1782	for (sd1 = sman.beginSelected(i); sd1 != NULL;
203			sd1 = sman.nextSelected(i)) {
204			for (j = i, sd2 = sman.nextSelected(j); sd2 != NULL;
205			sd2 = sman.nextSelected(j)) {
206	gezelter	507	collectHistogram(sd1, sd2);
207			}
208	tim	347	}
209	gezelter	507	}
210	tim	347	}