--- trunk/src/applications/hydrodynamics/AnalyticalModel.cpp 2006/03/17 23:20:35 906 +++ 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,109 +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; + } } - -/** - * calculate the ratio of friction coeffiction constant between ellipsoid and spheric - * with same volume. - * @param m - * @param n - * @note - * Reference: - * - * (1) Victor A. Bloomfield, On-Line Biophysics Textbook, Volume: Separations and Hydrodynamics - * Chapter 1,Survey of Biomolecular Hydrodynamics - * http://www.biophysics.org/education/vbloomfield.pdf - * (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 ft; - double fra; - double frb; - double a = ellipsoid->getA(); - double b = ellipsoid->getB(); - double q = a/b; //? - if (q > 1.0) {//prolate - ft = sqrt(1-q*q)/(pow(q, 2.0/3.0)*log((1 + sqrt(1-q*q))/q)); - fra = 4*(1-q*q)/(3*(2 - 2*pow(q, 4.0/3.0)/ft)); //not sure - frb = 4*(1-q*q*q*q) /(3*q*q*(2*pow(q, -2.0/3.0)*(2-q*q)/ft-2)); - } else {//oblate - ft = sqrt(1-q*q)/(pow(q, 2.0/3.0)*atan(sqrt(q*q-1))); - fra = 4*(1-q*q)/(3*(2 - 2*pow(q, 4.0/3.0)/ft)); //not sure - frb = 4*(1-q*q*q*q) /(3*q*q*(2*pow(q, -2.0/3.0)*(2-q*q)/ft-2)); - } - - double radius = pow(a*a*b, 1.0/3.0); - HydroProps props; - 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; - - const double convertConstant = 6.023; //convert poise.angstrom to amu/fs - props.Xi *= convertConstant; - props.Xi(0,0) *= ft; - props.Xi(1,1) *= ft; - props.Xi(2,2) *= ft; - props.Xi(3,3) *= fra; - props.Xi(4,4) *= fra; - props.Xi(5,5) *= frb; - - Mat6x6d XiCopy = props.Xi; - XiCopy /= OOPSEConstant::kb * temperature; - invertMatrix(XiCopy, props.D); - double kt = OOPSEConstant::kB * temperature; - props.D *= kt; - - setCR(props); - setCD(props); - - return true; -} - -bool AnalyticalModel::calcHydroProps(CompositeShape* compositexShape, double viscosity, double temperature) { - return false; -} - - - -}