--- trunk/src/applications/hydrodynamics/AnalyticalModel.cpp 2006/03/21 00:26:55 909 +++ trunk/src/applications/hydrodynamics/AnalyticalModel.cpp 2009/11/25 20:02:06 1390 @@ -6,19 +6,10 @@ * redistribute this software in source and binary code form, provided * that the following conditions are met: * - * 1. Acknowledgement of the program authors must be made in any - * publication of scientific results based in part on use of the - * program. An acceptable form of acknowledgement is citation of - * the article in which the program was described (Matthew - * A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher - * J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented - * Parallel Simulation Engine for Molecular Dynamics," - * J. Comput. Chem. 26, pp. 252-271 (2005)) - * - * 2. Redistributions of source code must retain the above copyright + * 1. Redistributions of source code must retain the above copyright * notice, this list of conditions and the following disclaimer. * - * 3. Redistributions in binary form must reproduce the above copyright + * 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. @@ -37,102 +28,63 @@ * 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] Vardeman & Gezelter, in progress (2009). */ #include "applications/hydrodynamics/AnalyticalModel.hpp" -#include "applications/hydrodynamics/Spheric.hpp" -#include "applications/hydrodynamics/Ellipsoid.hpp" +#include "hydrodynamics/Sphere.hpp" +#include "hydrodynamics/Ellipsoid.hpp" #include "applications/hydrodynamics/CompositeShape.hpp" #include "math/LU.hpp" -namespace oopse { -bool AnalyticalModel::calcHydroProps(Spheric* spheric, double viscosity, double temperature) { - - double radius = spheric->getRadius(); - HydroProps props; - props.center =V3Zero; - double Xitt = 6.0 * NumericConstant::PI * viscosity * radius; - double Xirr = 8.0 * NumericConstant::PI * viscosity * radius * radius * radius; - props.Xi(0, 0) = Xitt; - props.Xi(1, 1) = Xitt; - props.Xi(2, 2) = Xitt; - props.Xi(3, 3) = Xirr; - props.Xi(4, 4) = Xirr; - props.Xi(5, 5) = Xirr; +namespace OpenMD { + + bool AnalyticalModel::calcHydroProps(Shape* shape, RealType viscosity, RealType temperature) { - const double convertConstant = 6.023; //convert poise.angstrom to amu/fs - props.Xi *= convertConstant; - Mat6x6d XiCopy = props.Xi; - invertMatrix(XiCopy, props.D); - double kt = OOPSEConstant::kB * temperature; - props.D *= kt; - props.Xi *= OOPSEConstant::kb * temperature; - - setCR(props); - setCD(props); - - return true; + HydroProp* props; + Sphere* sphere = dynamic_cast(shape); + if (sphere != NULL) { + props = sphere->getHydroProp(viscosity, temperature); + setCR(props); + setCD(props); + return true; + } else { + Ellipsoid* ellipsoid = dynamic_cast(shape); + if (ellipsoid != NULL) { + props = ellipsoid->getHydroProp(viscosity, temperature); + setCR(props); + setCD(props); + return true; + } else { + CompositeShape* composite = dynamic_cast(shape); + if (composite != NULL) { +// props = composite->getHydroProp(viscosity, temperature); +// setCR(props); +// setCD(props); +// return true; + return false; + } else { + sprintf( painCave.errMsg, + "Could not figure out what kind of shape this is!\n"); + painCave.severity = OPENMD_ERROR; + painCave.isFatal = 1; + simError(); + return false; + } + } + } + } + void AnalyticalModel::writeBeads(std::ostream& os) { + os << "1\n"; + os << "Generated by Hydro\n"; + Vector3d pos = sd_->getPos(); + os << sd_->getType() << "\t" << pos[0] << "\t" << pos[1] << "\t" << pos[2] << std::endl; + } } - -/** - * Reference: - * (2) F. Perrin , J. Phys. Radium, [7] 5, 497-511, 1934 - * (3) F. Perrin, J. Phys. Radium, [7] 7, 1-11, 1936 - */ -bool AnalyticalModel::calcHydroProps(Ellipsoid* ellipsoid, double viscosity, double temperature) { - - double rMajor = ellipsoid->getRMajor(); - double rMinor = ellipsoid->getRMinor(); - - double a = rMinor; - double b = rMajor; - double a2 = a * a; - double b2 = b* b; - - double p = a /b; - double S; - if (p > 1.0) { //prolate - S = 2.0/sqrt(a2 - b2) * log((a + sqrt(a2-b2))/b); - } else { //oblate - S = 2.0/sqrt(b2 - a2) * atan(sqrt(b2-a2)/a); - } - - double P = 1.0/(a2 - b2) * (S - 2.0/a); - double Q = 0.5/(a2-b2) * (2.0*a/b2 - S); - - double transMinor = 16.0 * NumericConstant::PI * viscosity * (a2 - b2) /((2.0*a2-b2)*S -2.0*a); - double transMajor = 32.0 * NumericConstant::PI * viscosity * (a2 - b2) /((2.0*a2-3.0*b2)*S +2.0*a); - double rotMinor = 32.0/3.0 * NumericConstant::PI * viscosity *(a2 - b2) * b2 /(2.0*a -b2*S); - double rotMajor = 32.0/3.0 * NumericConstant::PI * viscosity *(a2*a2 - b2*b2)/((2.0*a2-b2)*S-2.0*a); - - - HydroProps props; - - props.Xi(0,0) = transMajor; - props.Xi(1,1) = transMajor; - props.Xi(2,2) = transMinor; - props.Xi(3,3) = rotMajor; - props.Xi(4,4) = rotMajor; - props.Xi(5,5) = rotMinor; - - const double convertConstant = 6.023; //convert poise.angstrom to amu/fs - props.Xi *= convertConstant; - - Mat6x6d XiCopy = props.Xi; - invertMatrix(XiCopy, props.D); - double kt = OOPSEConstant::kB * temperature; - props.D *= kt; - props.Xi *= OOPSEConstant::kb * temperature; - - setCR(props); - setCD(props); - - return true; -} - -bool AnalyticalModel::calcHydroProps(CompositeShape* compositexShape, double viscosity, double temperature) { - return false; -} - - - -}