ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/applications/hydrodynamics/AnalyticalModel.cpp
(Generate patch)

Comparing trunk/src/applications/hydrodynamics/AnalyticalModel.cpp (file contents):
Revision 908 by tim, Mon Mar 20 19:12:14 2006 UTC vs.
Revision 1879 by gezelter, Sun Jun 16 15:15:42 2013 UTC

# Line 6 | Line 6
6   * redistribute this software in source and binary code form, provided
7   * that the following conditions are met:
8   *
9 < * 1. Acknowledgement of the program authors must be made in any
10 < *    publication of scientific results based in part on use of the
11 < *    program.  An acceptable form of acknowledgement is citation of
12 < *    the article in which the program was described (Matthew
13 < *    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14 < *    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15 < *    Parallel Simulation Engine for Molecular Dynamics,"
16 < *    J. Comput. Chem. 26, pp. 252-271 (2005))
17 < *
18 < * 2. Redistributions of source code must retain the above copyright
9 > * 1. Redistributions of source code must retain the above copyright
10   *    notice, this list of conditions and the following disclaimer.
11   *
12 < * 3. Redistributions in binary form must reproduce the above copyright
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.
# Line 37 | Line 28
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   #include "applications/hydrodynamics/AnalyticalModel.hpp"
43 < #include "applications/hydrodynamics/Spheric.hpp"
44 < #include "applications/hydrodynamics/Ellipsoid.hpp"
43 > #include "hydrodynamics/Sphere.hpp"
44 > #include "hydrodynamics/Ellipsoid.hpp"
45   #include "applications/hydrodynamics/CompositeShape.hpp"
46   #include "math/LU.hpp"
47 < namespace oopse {
48 < bool AnalyticalModel::calcHydroProps(Spheric* spheric, double viscosity, double temperature) {
49 <
49 <    double radius = spheric->getRadius();
50 <    HydroProps props;
51 <    props.center =V3Zero;
52 <    double Xitt  = 6.0 * NumericConstant::PI * viscosity * radius;
53 <    double Xirr = 8.0 * NumericConstant::PI * viscosity * radius * radius * radius;
54 <    props.Xi(0, 0) = Xitt;
55 <    props.Xi(1, 1) = Xitt;
56 <    props.Xi(2, 2) = Xitt;
57 <    props.Xi(3, 3) = Xirr;
58 <    props.Xi(4, 4) = Xirr;
59 <    props.Xi(5, 5) = Xirr;
47 > namespace OpenMD {
48 >  
49 >  bool AnalyticalModel::calcHydroProps(Shape* shape, RealType viscosity, RealType temperature) {
50      
51 <    const double convertConstant = 6.023; //convert poise.angstrom to amu/fs
52 <    props.Xi *= convertConstant;
53 <    Mat6x6d XiCopy = props.Xi;
54 <    invertMatrix(XiCopy, props.D);
55 <    double kt = OOPSEConstant::kB * temperature;
56 <    props.D *= kt;
57 <    props.Xi *= OOPSEConstant::kb * temperature;
58 <
59 <    setCR(props);
60 <    setCD(props);
61 <
62 <    return true;
51 >    HydroProp* props;
52 >    Sphere* sphere = dynamic_cast<Sphere*>(shape);        
53 >    if (sphere != NULL) {
54 >      props = sphere->getHydroProp(viscosity, temperature);
55 >      setCR(props);
56 >      setCD(props);
57 >      return true;
58 >    } else {
59 >      Ellipsoid* ellipsoid = dynamic_cast<Ellipsoid*>(shape);        
60 >      if (ellipsoid != NULL) {
61 >        props = ellipsoid->getHydroProp(viscosity, temperature);
62 >        setCR(props);
63 >        setCD(props);
64 >        return true;
65 >      } else {
66 >        CompositeShape* composite = dynamic_cast<CompositeShape*>(shape);
67 >        if (composite != NULL) {
68 > //           props = composite->getHydroProp(viscosity, temperature);
69 > //           setCR(props);
70 > //           setCD(props);
71 > //           return true;
72 >          return false;
73 >        } else {
74 >          sprintf( painCave.errMsg,
75 >                   "Could not figure out what kind of shape this is!\n");
76 >          painCave.severity = OPENMD_ERROR;
77 >          painCave.isFatal = 1;
78 >          simError();    
79 >          return false;
80 >        }
81 >      }
82 >    }          
83 >  }
84      
85 +  void AnalyticalModel::writeBeads(std::ostream& os) {
86 +    os << "1\n";
87 +    os << "Generated by Hydro\n";
88 +    Vector3d pos = sd_->getPos();
89 +    os << sd_->getType() << "\t" << pos[0] << "\t" << pos[1] << "\t" << pos[2] << std::endl;
90 +  }  
91   }
75
76 /**
77 * Reference:
78 * (2) F. Perrin , J. Phys. Radium, [7] 5, 497-511, 1934
79 * (3) F. Perrin, J. Phys. Radium, [7] 7, 1-11, 1936
80 */        
81 bool AnalyticalModel::calcHydroProps(Ellipsoid* ellipsoid, double viscosity, double temperature) {
82
83    double rMajor = ellipsoid->getRMajor();
84    double rMinor = ellipsoid->getRMinor();
85
86    double a = rMinor;
87    double b = rMajor;
88    double a2 = a * a;
89    double b2 = b* b;
90    
91    double p = a /b;
92    double S;
93    if (p > 1.0) { //prolate
94        S = 2.0/sqrt(a2 - b2) * log((a + sqrt(a2-b2))/b);
95    } { //oblate
96        S = 2.0/sqrt(b2 - a2) * atan(sqrt(b2-a2)/a);
97    }
98
99    double P = 1.0/(a2 - b2) * (S - 2.0/a);
100    double Q = 0.5/(a2-b2) * (2.0*a/b2 - S);
101
102    double transMinor = 16.0 * NumericConstant::PI * viscosity * (a2 - b2) /((2.0*a2-b2)*S -2.0*a);
103    double transMajor = 32.0 * NumericConstant::PI * viscosity * (a2 - b2) /((2.0*a2-3.0*b2)*S +2.0*a);
104    double rotMinor = 32.0/3.0 * NumericConstant::PI * viscosity *(a2 - b2) * b2 /(2.0*a -b2*S);
105    double rotMajor = 32.0/3.0 * NumericConstant::PI * viscosity *(a2*a2 - b2*b2)/((2.0*a2-b2)*S-2.0*a);
106    
107        
108    HydroProps props;
109
110    props.Xi(0,0) = transMajor;
111    props.Xi(1,1) = transMajor;
112    props.Xi(2,2) = transMinor;
113    props.Xi(3,3) = rotMajor;
114    props.Xi(4,4) = rotMajor;
115    props.Xi(5,5) = rotMinor;
116    
117    const double convertConstant = 6.023; //convert poise.angstrom to amu/fs
118    props.Xi *= convertConstant;    
119    
120    Mat6x6d XiCopy = props.Xi;
121    invertMatrix(XiCopy, props.D);
122    double kt = OOPSEConstant::kB * temperature;
123    props.D *= kt;
124    props.Xi *= OOPSEConstant::kb * temperature;
125    
126    setCR(props);
127    setCD(props);
128
129    return true;
130 }
131
132 bool AnalyticalModel::calcHydroProps(CompositeShape* compositexShape, double viscosity, double temperature) {
133    return false;
134 }
135        
136
137
138 }

Comparing trunk/src/applications/hydrodynamics/AnalyticalModel.cpp (property svn:keywords):
Revision 908 by tim, Mon Mar 20 19:12:14 2006 UTC vs.
Revision 1879 by gezelter, Sun Jun 16 15:15:42 2013 UTC

# Line 0 | Line 1
1 + Author Id Revision Date

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines