--- trunk/src/applications/hydrodynamics/HydrodynamicsModel.cpp 2006/02/24 21:17:05 893 +++ trunk/src/applications/hydrodynamics/HydrodynamicsModel.cpp 2006/03/15 20:58:47 900 @@ -77,7 +77,16 @@ bool HydrodynamicsModel::calcHydrodyanmicsProps() { std::cout << "can not create beads" << std::endl; return false; } - + + //calcResistanceTensor(); + calcDiffusionTensor(); + return true; +} + +void HydrodynamicsModel::calcResistanceTensor() { +} + +void HydrodynamicsModel::calcDiffusionTensor() { int nbeads = beads_.size(); DynamicRectMatrix B(3*nbeads, 3*nbeads); DynamicRectMatrix C(3*nbeads, 3*nbeads); @@ -105,18 +114,12 @@ bool HydrodynamicsModel::calcHydrodyanmicsProps() { Tij(2, 2) = constant; } B.setSubMatrix(i*3, j*3, Tij); - std::cout << Tij << std::endl; } } - std::cout << "B=\n" - << B << std::endl; //invert B Matrix invertMatrix(B, C); - std::cout << "C=\n" - << C << std::endl; - //prepare U Matrix relative to arbitrary origin O(0.0, 0.0, 0.0) std::vector U; for (int i = 0; i < nbeads; ++i) { @@ -134,7 +137,7 @@ bool HydrodynamicsModel::calcHydrodyanmicsProps() { double volume = 0.0; for (std::vector::iterator iter = beads_.begin(); iter != beads_.end(); ++iter) { - volume = 4.0/3.0 * NumericConstant::PI * pow((*iter).radius,3); + volume += 4.0/3.0 * NumericConstant::PI * pow((*iter).radius,3); } for (std::size_t i = 0; i < nbeads; ++i) { @@ -144,22 +147,16 @@ bool HydrodynamicsModel::calcHydrodyanmicsProps() { Xitt += Cij; Xitr += U[i] * Cij; - Xirr += -U[i] * Cij * U[j]; - //Xirr += -U[i] * Cij * U[j] + (0.166*6 * viscosity_ * volume) * I; + Xirr += -U[i] * Cij * U[j] + (6 * viscosity_ * volume) * I; } } - //invert Xi to get Diffusion Tensor at arbitrary origin O - RectMatrix Xi; - RectMatrix Do; - Xi.setSubMatrix(0, 0, Xitt); - Xi.setSubMatrix(0, 3, Xitr.transpose()); - Xi.setSubMatrix(3, 0, Xitr); - Xi.setSubMatrix(3, 3, Xirr); - //invertMatrix(Xi, Do); - double kt = OOPSEConstant::kB * temperature_ * 1.66E-2; - //Do *= kt; + const double convertConstant = 6.023; //convert poise.angstrom to amu/fs + Xitt *= convertConstant; + Xitr *= convertConstant; + Xirr *= convertConstant; + double kt = OOPSEConstant::kB * temperature_; Mat3x3d Dott; //translational diffusion tensor at arbitrary origin O Mat3x3d Dorr; //rotational diffusion tensor at arbitrary origin O @@ -170,10 +167,6 @@ bool HydrodynamicsModel::calcHydrodyanmicsProps() { Mat3x3d XittInv(0.0); XittInv = Xitt.inverse(); - //Xirr may not be inverted,if it one of the diagonal element is zero, for example - //( a11 a12 0) - //( a21 a22 0) - //( 0 0 0) Mat3x3d XirrInv; XirrInv = Xirr.inverse(); @@ -182,18 +175,14 @@ bool HydrodynamicsModel::calcHydrodyanmicsProps() { tmp = Xitt - Xitr.transpose() * XirrInv * Xitr; tmpInv = tmp.inverse(); - Dott = kt * tmpInv; - Dotr = -kt*XirrInv * Xitr * tmpInv* 1.0E8; - + Dott = kt*tmpInv; + Dotr = -kt*XirrInv * Xitr * tmpInv; + tmp = Xirr - Xitr * XittInv * Xitr.transpose(); tmpInv = tmp.inverse(); - Dorr = kt * tmpInv*1.0E16; + Dorr = kt * tmpInv; - //Do.getSubMatrix(0, 0 , Dott); - //Do.getSubMatrix(3, 0, Dotr); - //Do.getSubMatrix(3, 3, Dorr); - //calculate center of diffusion tmp(0, 0) = Dorr(1, 1) + Dorr(2, 2); tmp(0, 1) = - Dorr(0, 1); @@ -225,13 +214,50 @@ bool HydrodynamicsModel::calcHydrodyanmicsProps() { Ddtt = Dott - Uod * Dorr * Uod + Dotr.transpose() * Uod - Uod * Dotr; Ddrr = Dorr; Ddtr = Dotr + Dorr * Uod; - + props_.diffCenter = rod; - props_.transDiff = Ddtt; - props_.transRotDiff = Ddtr; - props_.rotDiff = Ddrr; + props_.Ddtt = Ddtt; + props_.Ddtr = Ddtr; + props_.Ddrr = Ddrr; - return true; + SquareMatrix Dd; + Dd.setSubMatrix(0, 0, Ddtt); + Dd.setSubMatrix(0, 3, Ddtr.transpose()); + Dd.setSubMatrix(3, 0, Ddtr); + Dd.setSubMatrix(3, 3, Ddrr); + SquareMatrix Xid; + invertMatrix(Dd, Xid); + + //Xidtt in units of kcal*fs*mol^-1*Ang^-2 + Xid *= OOPSEConstant::kb*temperature_; + + Xid.getSubMatrix(0, 0, props_.Xidtt); + Xid.getSubMatrix(0, 3, props_.Xidrt); + Xid.getSubMatrix(3, 0, props_.Xidtr); + Xid.getSubMatrix(3, 3, props_.Xidrr); + + + std::cout << "center of diffusion :" << std::endl; + std::cout << rod << std::endl; + std::cout << "diffusion tensor at center of diffusion " << std::endl; + std::cout << "translation(A^2/fs) :" << std::endl; + std::cout << Ddtt << std::endl; + std::cout << "translation-rotation(A^3/fs):" << std::endl; + std::cout << Ddtr << std::endl; + std::cout << "rotation(A^4/fs):" << std::endl; + std::cout << Ddrr << std::endl; + + std::cout << "resistance tensor at center of diffusion " << std::endl; + std::cout << "translation(kcal*fs*mol^-1*Ang^-2) :" << std::endl; + std::cout << props_.Xidtt << std::endl; + std::cout << "rotation-translation (kcal*fs*mol^-1*Ang^-3):" << std::endl; + std::cout << props_.Xidrt << std::endl; + std::cout << "translation-rotation(kcal*fs*mol^-1*Ang^-3):" << std::endl; + std::cout << props_.Xidtr << std::endl; + std::cout << "rotation(kcal*fs*mol^-1*Ang^-4):" << std::endl; + std::cout << props_.Xidrr << std::endl; + + } void HydrodynamicsModel::writeBeads(std::ostream& os) { @@ -245,39 +271,38 @@ void HydrodynamicsModel::writeDiffCenterAndDiffTensor( } void HydrodynamicsModel::writeDiffCenterAndDiffTensor(std::ostream& os) { - os << "//viscosity = " << viscosity_ << std::endl; - os << "//temperature = " << temperature_<< std::endl; - std::vector::iterator iter; - os << sd_->getType() << "\n"; - os << "//diffusion center" << std::endl; - os << props_.diffCenter << std::endl; + os << sd_->getType() << "\t"; + os << props_.diffCenter[0] << "\t" << props_.diffCenter[1] << "\t" << props_.diffCenter[2] << "\t"; - os << "//translational diffusion tensor" << std::endl; - os << props_.transDiff << std::endl; + os << props_.Ddtt(0, 0) << "\t" << props_.Ddtt(0, 1) << "\t" << props_.Ddtt(0, 2) << "\t" + << props_.Ddtt(1, 0) << "\t" << props_.Ddtt(1, 1) << "\t" << props_.Ddtt(1, 2) << "\t" + << props_.Ddtt(2, 0) << "\t" << props_.Ddtt(2, 1) << "\t" << props_.Ddtt(2, 2) << "\t"; + + os << props_.Ddtr(0, 0) << "\t" << props_.Ddtr(0, 1) << "\t" << props_.Ddtr(0, 2) << "\t" + << props_.Ddtr(1, 0) << "\t" << props_.Ddtr(1, 1) << "\t" << props_.Ddtr(1, 2) << "\t" + << props_.Ddtr(2, 0) << "\t" << props_.Ddtr(2, 1) << "\t" << props_.Ddtr(2, 2) << "\t"; - os << "//translation-rotation coupling diffusion tensor" << std::endl; - os << props_.transRotDiff << std::endl; + os << props_.Ddrr(0, 0) << "\t" << props_.Ddrr(0, 1) << "\t" << props_.Ddrr(0, 2) << "\t" + << props_.Ddrr(1, 0) << "\t" << props_.Ddrr(1, 1) << "\t" << props_.Ddrr(1, 2) << "\t" + << props_.Ddrr(2, 0) << "\t" << props_.Ddrr(2, 1) << "\t" << props_.Ddrr(2, 2) <<"\t"; - os << "//rotational diffusion tensor" << std::endl; - os << props_.rotDiff << std::endl; - - /* - os << props_.diffCenter[0] << "\t" << props_.diffCenter[1] << "\t" << props_.diffCenter[2] << "\n" + os << props_.Xidtt(0, 0) << "\t" << props_.Xidtt(0, 1) << "\t" << props_.Xidtt(0, 2) << "\t" + << props_.Xidtt(1, 0) << "\t" << props_.Xidtt(1, 1) << "\t" << props_.Xidtt(1, 2) << "\t" + << props_.Xidtt(2, 0) << "\t" << props_.Xidtt(2, 1) << "\t" << props_.Xidtt(2, 2) << "\t"; - os << props_.transDiff(0, 0) << "\t" << props_.transDiff(0, 1) << "\t" << props_.transDiff(0, 2) << "\t" - << props_.transDiff(1, 0) << "\t" << props_.transDiff(1, 1) << "\t" << props_.transDiff(1, 2) << "\t" - << props_.transDiff(2, 0) << "\t" << props_.transDiff(2, 1) << "\t" << props_.transDiff(2, 2) << "\n"; + os << props_.Xidrt(0, 0) << "\t" << props_.Xidrt(0, 1) << "\t" << props_.Xidrt(0, 2) << "\t" + << props_.Xidrt(1, 0) << "\t" << props_.Xidrt(1, 1) << "\t" << props_.Xidrt(1, 2) << "\t" + << props_.Xidrt(2, 0) << "\t" << props_.Xidrt(2, 1) << "\t" << props_.Xidrt(2, 2) << "\t"; - os << props_.transRotDiff(0, 0) << "\t" << props_.transRotDiff(0, 1) << "\t" << props_.transRotDiff(0, 2) << "\t" - << props_.transRotDiff(1, 0) << "\t" << props_.transRotDiff(1, 1) << "\t" << props_.transRotDiff(1, 2) << "\t" - << props_.transRotDiff(2, 0) << "\t" << props_.transRotDiff(2, 1) << "\t" << props_.transRotDiff(2, 2) << "\t" + os << props_.Xidtr(0, 0) << "\t" << props_.Xidtr(0, 1) << "\t" << props_.Xidtr(0, 2) << "\t" + << props_.Xidtr(1, 0) << "\t" << props_.Xidtr(1, 1) << "\t" << props_.Xidtr(1, 2) << "\t" + << props_.Xidtr(2, 0) << "\t" << props_.Xidtr(2, 1) << "\t" << props_.Xidtr(2, 2) << "\t"; - os << props_.rotDiff(0, 0) << "\t" << props_.rotDiff(0, 1) << "\t" << props_.rotDiff(0, 2) << "\t" - << props_.rotDiff(1, 0) << "\t" << props_.rotDiff(1, 1) << "\t" << props_.rotDiff(1, 2) << "\t" - << props_.rotDiff(2, 0) << "\t" << props_.rotDiff(2, 1) << "\t" << props_.rotDiff(2, 2) << ";" - << std::endl; - */ + os << props_.Xidrr(0, 0) << "\t" << props_.Xidrr(0, 1) << "\t" << props_.Xidrr(0, 2) << "\t" + << props_.Xidrr(1, 0) << "\t" << props_.Xidrr(1, 1) << "\t" << props_.Xidrr(1, 2) << "\t" + << props_.Xidrr(2, 0) << "\t" << props_.Xidrr(2, 1) << "\t" << props_.Xidrr(2, 2) << std::endl; + } }