applications/sequentialProps/ContactAngle2.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, 234107 (2008).          
 * [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
 * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
 */

#include <algorithm>
#include <functional>
#include "applications/sequentialProps/ContactAngle2.hpp"
#include "utils/simError.h"
#include "io/DumpReader.hpp"
#include "primitives/Molecule.hpp"
#include "utils/NumericConstant.hpp"
#include "utils/PhysicalConstants.hpp"
#include "math/Polynomial.hpp"

namespace OpenMD {

  ContactAngle2::ContactAngle2(SimInfo* info, const std::string& filename, 
                               const std::string& sele, RealType solidZ,
                               RealType threshDens, int nrbins, int nzbins)
    : SequentialAnalyzer(info, filename), selectionScript_(sele), 
      evaluator_(info), seleMan_(info), solidZ_(solidZ),
      threshDens_(threshDens), nRBins_(nrbins), nZBins_(nzbins) {
    
    setOutputName(getPrefix(filename) + ".ca2");
    
    evaluator_.loadScriptString(sele);
    
    if (!evaluator_.isDynamic()) {
      seleMan_.setSelectionSet(evaluator_.evaluate());
    }            
  }

  void ContactAngle2::doFrame() {
    StuntDouble* sd;
    int i;

    // set up the bins for density analysis

    Mat3x3d hmat = info_->getSnapshotManager()->getCurrentSnapshot()->getHmat();
    RealType len = std::min(hmat(0, 0), hmat(1, 1));
    RealType zLen = hmat(2,2);
    RealType dr = len / (RealType) nRBins_;
    RealType dz = zLen / (RealType) nZBins_;

    std::vector<std::vector<RealType> > histo;
    histo.resize(nRBins_);
    for (int i = 0 ; i < nRBins_; ++i) {
      histo[i].resize(nZBins_);
    }
    for (unsigned int i = 0; i < histo.size(); ++i){
      std::fill(histo[i].begin(), histo[i].end(), 0.0);
    }      
        
    if (evaluator_.isDynamic()) {
      seleMan_.setSelectionSet(evaluator_.evaluate());
    }
    

    RealType mtot = 0.0;
    Vector3d com(V3Zero);
    RealType mass;
    
    for (sd = seleMan_.beginSelected(i); sd != NULL;
         sd = seleMan_.nextSelected(i)) {      
      mass = sd->getMass();
      mtot += mass;
      com += sd->getPos() * mass;
    }

    com /= mtot;

    // now that we have the centroid, we can make cylindrical density maps
    Vector3d pos;
    RealType r;
    RealType z;
    
    for (sd = seleMan_.beginSelected(i); sd != NULL;
         sd = seleMan_.nextSelected(i)) {      
      pos = sd->getPos() - com;
      
      r = sqrt(pow(pos.x(), 2) + pow(pos.y(), 2));
      z = pos.z() - solidZ_;
      
      int whichRBin = int(r / dr);
      int whichZBin = int(z/ dz);
      
      if ((r <= len) && (z <= zLen)) 
        histo[whichRBin][whichZBin] += sd->getMass();
      
    }
    
    for(unsigned int i = 0 ; i < histo.size(); ++i){

      RealType rL = i * dr;
      RealType rU = rL + dr;
      RealType volSlice = NumericConstant::PI * dz * (( rU*rU ) - ( rL*rL ));

      for (unsigned int j = 0; j < histo[i].size(); ++j){
        histo[i][j] *= PhysicalConstants::densityConvert / volSlice;
      }
    }

    for (unsigned int j = 0; j < nZBins_;  ++j) {
      RealType thez = dz * (j + 0.5);
      bool aboveThresh = false;
      for (unsigned int i = 0; i < nRBins_;  ++i) {
        RealType ther = dr * (i + 0.5);
        if (histo[i][j] >= threshDens_) aboveThresh = true;

        if (aboveThresh && (histo[i][j] <= threshDens_)) {
          cerr << thez << "\t" << ther << "\n";
          break;         
        }
      }
    }
    
    // values_.push_back( acos(maxct)*(180.0/M_PI) );
    
  }   
}


Revision:	2036
Committed:	Tue Nov 4 16:20:31 2014 UTC (10 years, 5 months ago) by gezelter
File size:	5549 byte(s)
Log Message:	A scan fix for CA2
#	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, 234107 (2008).
39	* [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010).
40	* [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41	*/
42
43	#include <algorithm>
44	#include <functional>
45	#include "applications/sequentialProps/ContactAngle2.hpp"
46	#include "utils/simError.h"
47	#include "io/DumpReader.hpp"
48	#include "primitives/Molecule.hpp"
49	#include "utils/NumericConstant.hpp"
50	#include "utils/PhysicalConstants.hpp"
51	#include "math/Polynomial.hpp"
52
53	namespace OpenMD {
54
55	ContactAngle2::ContactAngle2(SimInfo* info, const std::string& filename,
56	const std::string& sele, RealType solidZ,
57	RealType threshDens, int nrbins, int nzbins)
58	: SequentialAnalyzer(info, filename), selectionScript_(sele),
59	evaluator_(info), seleMan_(info), solidZ_(solidZ),
60	threshDens_(threshDens), nRBins_(nrbins), nZBins_(nzbins) {
61
62	setOutputName(getPrefix(filename) + ".ca2");
63
64	evaluator_.loadScriptString(sele);
65
66	if (!evaluator_.isDynamic()) {
67	seleMan_.setSelectionSet(evaluator_.evaluate());
68	}
69	}
70
71	void ContactAngle2::doFrame() {
72	StuntDouble* sd;
73	int i;
74
75	// set up the bins for density analysis
76
77	Mat3x3d hmat = info_->getSnapshotManager()->getCurrentSnapshot()->getHmat();
78	RealType len = std::min(hmat(0, 0), hmat(1, 1));
79	RealType zLen = hmat(2,2);
80	RealType dr = len / (RealType) nRBins_;
81	RealType dz = zLen / (RealType) nZBins_;
82
83	std::vector<std::vector<RealType> > histo;
84	histo.resize(nRBins_);
85	for (int i = 0 ; i < nRBins_; ++i) {
86	histo[i].resize(nZBins_);
87	}
88	for (unsigned int i = 0; i < histo.size(); ++i){
89	std::fill(histo[i].begin(), histo[i].end(), 0.0);
90	}
91
92	if (evaluator_.isDynamic()) {
93	seleMan_.setSelectionSet(evaluator_.evaluate());
94	}
95
96
97	RealType mtot = 0.0;
98	Vector3d com(V3Zero);
99	RealType mass;
100
101	for (sd = seleMan_.beginSelected(i); sd != NULL;
102	sd = seleMan_.nextSelected(i)) {
103	mass = sd->getMass();
104	mtot += mass;
105	com += sd->getPos() * mass;
106	}
107
108	com /= mtot;
109
110	// now that we have the centroid, we can make cylindrical density maps
111	Vector3d pos;
112	RealType r;
113	RealType z;
114
115	for (sd = seleMan_.beginSelected(i); sd != NULL;
116	sd = seleMan_.nextSelected(i)) {
117	pos = sd->getPos() - com;
118
119	r = sqrt(pow(pos.x(), 2) + pow(pos.y(), 2));
120	z = pos.z() - solidZ_;
121
122	int whichRBin = int(r / dr);
123	int whichZBin = int(z/ dz);
124
125	if ((r <= len) && (z <= zLen))
126	histo[whichRBin][whichZBin] += sd->getMass();
127
128	}
129
130	for(unsigned int i = 0 ; i < histo.size(); ++i){
131
132	RealType rL = i * dr;
133	RealType rU = rL + dr;
134	RealType volSlice = NumericConstant::PI * dz * (( rUrU ) - ( rLrL ));
135
136	for (unsigned int j = 0; j < histo[i].size(); ++j){
137	histo[i][j] *= PhysicalConstants::densityConvert / volSlice;
138	}
139	}
140
141	for (unsigned int j = 0; j < nZBins_; ++j) {
142	RealType thez = dz * (j + 0.5);
143	bool aboveThresh = false;
144	for (unsigned int i = 0; i < nRBins_; ++i) {
145	RealType ther = dr * (i + 0.5);
146	if (histo[i][j] >= threshDens_) aboveThresh = true;
147
148	if (aboveThresh && (histo[i][j] <= threshDens_)) {
149	cerr << thez << "\t" << ther << "\n";
150	break;
151	}
152	}
153	}
154
155	// values_.push_back( acos(maxct)*(180.0/M_PI) );
156
157	}
158	}
159
160