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 i = 0; i < histo.size(); ++i) {
      RealType ther = dr * (i + 0.5);
      for(unsigned int j = 0; j < histo[i].size(); ++j) {
        if (histo[i][j] <= threshDens_) {
          RealType thez = dz * (j + 0.5);
          cerr << ther << "\t" << thez << "\n";
          break;
        }
      }
    }

    // values_.push_back( acos(maxct)*(180.0/M_PI) );
    
  }   
}


Revision:	2035
Committed:	Tue Nov 4 15:31:51 2014 UTC (10 years, 5 months ago) by gezelter
File size:	5438 byte(s)
Log Message:	First checkin of the second contact angle method for sequentialProps.
#	User	Rev	Content
1	gezelter	2035	/*
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 i = 0; i < histo.size(); ++i) {
142			RealType ther = dr * (i + 0.5);
143			for(unsigned int j = 0; j < histo[i].size(); ++j) {
144			if (histo[i][j] <= threshDens_) {
145			RealType thez = dz * (j + 0.5);
146			cerr << ther << "\t" << thez << "\n";
147			break;
148			}
149			}
150			}
151
152			// values_.push_back( acos(maxct)*(180.0/M_PI) );
153
154			}
155			}
156
157