--- trunk/src/applications/hydrodynamics/HydrodynamicsModel.cpp 2006/02/23 23:16:43 892 +++ 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); @@ -95,7 +104,7 @@ bool HydrodynamicsModel::calcHydrodyanmicsProps() { double rij2 = rij * rij; double sumSigma2OverRij2 = ((beads_[i].radius*beads_[i].radius) + (beads_[i].radius*beads_[i].radius)) / rij2; Mat3x3d tmpMat; - tmpMat = outProduct(beads_[i].pos, beads_[j].pos) / rij2; + tmpMat = outProduct(Rij, Rij) / rij2; double constant = 8.0 * NumericConstant::PI * viscosity_ * rij; Tij = ((1.0 + sumSigma2OverRij2/3.0) * I + (1.0 - sumSigma2OverRij2) * tmpMat ) / constant; }else { @@ -110,6 +119,7 @@ bool HydrodynamicsModel::calcHydrodyanmicsProps() { //invert B Matrix invertMatrix(B, C); + //prepare U Matrix relative to arbitrary origin O(0.0, 0.0, 0.0) std::vector U; for (int i = 0; i < nbeads; ++i) { @@ -122,6 +132,13 @@ bool HydrodynamicsModel::calcHydrodyanmicsProps() { Mat3x3d Xitt; Mat3x3d Xirr; Mat3x3d Xitr; + + //calculate the total volume + + 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); + } for (std::size_t i = 0; i < nbeads; ++i) { for (std::size_t j = 0; j < nbeads; ++j) { @@ -129,59 +146,43 @@ bool HydrodynamicsModel::calcHydrodyanmicsProps() { C.getSubMatrix(i*3, j*3, Cij); Xitt += Cij; - Xirr += U[i] * Cij; - Xitr += U[i] * Cij * U[j]; + Xitr += U[i] * Cij; + 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 Mat3x3d Dotr; //translation-rotation couplingl diffusion tensor at arbitrary origin O - Mat3x3d XirrInv(0.0); - Mat3x3d XirrCopy; - XirrCopy = Xirr; + const static Mat3x3d zeroMat(0.0); Mat3x3d XittInv(0.0); - Mat3x3d XittCopy; - XittCopy = Xitt; - invertMatrix(XittCopy, XittInv); + XittInv = Xitt.inverse(); + + Mat3x3d XirrInv; + XirrInv = Xirr.inverse(); Mat3x3d tmp; Mat3x3d tmpInv; tmp = Xitt - Xitr.transpose() * XirrInv * Xitr; + tmpInv = tmp.inverse(); - const static Mat3x3d zeroMat(0.0); - if (!invertMatrix(tmp, tmpInv)) { - tmpInv = zeroMat; - } - - Dott = kt * tmpInv; + Dott = kt*tmpInv; Dotr = -kt*XirrInv * Xitr * tmpInv; - - tmp = Xirr - Xitr * XittInv * Xitr.transpose(); - if(!invertMatrix(tmp, tmpInv)) { - tmpInv = zeroMat; - } + tmp = Xirr - Xitr * XittInv * Xitr.transpose(); + tmpInv = tmp.inverse(); + 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); @@ -198,9 +199,7 @@ bool HydrodynamicsModel::calcHydrodyanmicsProps() { tmpVec[1] = Dotr(2, 0) - Dotr(0, 2); tmpVec[2] = Dotr(0, 1) - Dotr(1, 0); - if(!invertMatrix(tmp, tmpInv)) { - tmpInv = zeroMat; - } + tmpInv = tmp.inverse(); Vector3d rod = tmpInv * tmpVec; @@ -215,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) { @@ -235,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 << "//translational 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; + } }