--- trunk/OOPSE/libmdtools/NPT.cpp 2003/09/19 20:00:27 778 +++ trunk/OOPSE/libmdtools/NPT.cpp 2004/04/12 20:32:20 1097 @@ -1,4 +1,5 @@ -#include +#include + #include "Atom.hpp" #include "SRI.hpp" #include "AbstractClasses.hpp" @@ -7,7 +8,7 @@ #include "Thermo.hpp" #include "ReadWrite.hpp" #include "Integrator.hpp" -#include "simError.h" +#include "simError.h" #ifdef IS_MPI #include "mpiSimulation.hpp" @@ -18,15 +19,22 @@ // modification of the Hoover algorithm: // // Melchionna, S., Ciccotti, G., and Holian, B. L., 1993, -// Molec. Phys., 78, 533. +// Molec. Phys., 78, 533. // // and -// +// // Hoover, W. G., 1986, Phys. Rev. A, 34, 2499. template NPT::NPT ( SimInfo *theInfo, ForceFields* the_ff): T( theInfo, the_ff ) { + GenericData* data; + DoubleData * chiValue; + DoubleData * integralOfChidtValue; + + chiValue = NULL; + integralOfChidtValue = NULL; + chi = 0.0; integralOfChidt = 0.0; have_tau_thermostat = 0; @@ -37,9 +45,26 @@ template NPT::NPT ( SimInfo *theInfo, F have_eta_tolerance = 0; have_pos_iter_tolerance = 0; - oldPos = new double[3*nAtoms]; - oldVel = new double[3*nAtoms]; - oldJi = new double[3*nAtoms]; + // retrieve chi and integralOfChidt from simInfo + data = info->getProperty(CHIVALUE_ID); + if(data){ + chiValue = dynamic_cast(data); + } + + data = info->getProperty(INTEGRALOFCHIDT_ID); + if(data){ + integralOfChidtValue = dynamic_cast(data); + } + + // chi and integralOfChidt should appear by pair + if(chiValue && integralOfChidtValue){ + chi = chiValue->getData(); + integralOfChidt = integralOfChidtValue->getData(); + } + + oldPos = new double[3*integrableObjects.size()]; + oldVel = new double[3*integrableObjects.size()]; + oldJi = new double[3*integrableObjects.size()]; #ifdef IS_MPI Nparticles = mpiSim->getTotAtoms(); #else @@ -55,10 +80,9 @@ template void NPT::moveA() { } template void NPT::moveA() { - + //new version of NPT int i, j, k; - DirectionalAtom* dAtom; double Tb[3], ji[3]; double mass; double vel[3], pos[3], frc[3]; @@ -66,54 +90,56 @@ template void NPT::moveA() { double COM[3]; instaTemp = tStats->getTemperature(); - instaPress = tStats->getPressure(); + tStats->getPressureTensor( press ); + instaPress = p_convert * (press[0][0] + press[1][1] + press[2][2]) / 3.0; instaVol = tStats->getVolume(); - + tStats->getCOM(COM); - + //evolve velocity half step - for( i=0; igetVel( vel ); - atoms[i]->getFrc( frc ); + calcVelScale(); + + for( i=0; igetMass(); + integrableObjects[i]->getVel( vel ); + integrableObjects[i]->getFrc( frc ); + mass = integrableObjects[i]->getMass(); + getVelScaleA( sc, vel ); for (j=0; j < 3; j++) { - + // velocity half step (use chi from previous step here): vel[j] += dt2 * ((frc[j] / mass ) * eConvert - sc[j]); - + } - atoms[i]->setVel( vel ); - - if( atoms[i]->isDirectional() ){ + integrableObjects[i]->setVel( vel ); - dAtom = (DirectionalAtom *)atoms[i]; + if( integrableObjects[i]->isDirectional() ){ // get and convert the torque to body frame - - dAtom->getTrq( Tb ); - dAtom->lab2Body( Tb ); - + + integrableObjects[i]->getTrq( Tb ); + integrableObjects[i]->lab2Body( Tb ); + // get the angular momentum, and propagate a half step - dAtom->getJ( ji ); + integrableObjects[i]->getJ( ji ); - for (j=0; j < 3; j++) + for (j=0; j < 3; j++) ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi); - - this->rotationPropagation( dAtom, ji ); - - dAtom->setJ( ji ); - } + + this->rotationPropagation( integrableObjects[i], ji ); + + integrableObjects[i]->setJ( ji ); + } } // evolve chi and eta half step - + evolveChiA(); evolveEtaA(); @@ -121,64 +147,61 @@ template void NPT::moveA() { integralOfChidt += dt2*chi; //save the old positions - for(i = 0; i < nAtoms; i++){ - atoms[i]->getPos(pos); + for(i = 0; i < integrableObjects.size(); i++){ + integrableObjects[i]->getPos(pos); for(j = 0; j < 3; j++) oldPos[i*3 + j] = pos[j]; } - - //the first estimation of r(t+dt) is equal to r(t) - + + //the first estimation of r(t+dt) is equal to r(t) + for(k = 0; k < 5; k ++){ - for(i =0 ; i < nAtoms; i++){ + for(i =0 ; i < integrableObjects.size(); i++){ - atoms[i]->getVel(vel); - atoms[i]->getPos(pos); + integrableObjects[i]->getVel(vel); + integrableObjects[i]->getPos(pos); this->getPosScale( pos, COM, i, sc ); - + for(j = 0; j < 3; j++) pos[j] = oldPos[i*3 + j] + dt*(vel[j] + sc[j]); - atoms[i]->setPos( pos ); + integrableObjects[i]->setPos( pos ); } - + if (nConstrained){ constrainA(); } } - + // Scale the box after all the positions have been moved: - + this->scaleSimBox(); } template void NPT::moveB( void ){ - + //new version of NPT int i, j, k; - DirectionalAtom* dAtom; double Tb[3], ji[3], sc[3]; double vel[3], frc[3]; double mass; - + // Set things up for the iteration: - for( i=0; igetVel( vel ); + integrableObjects[i]->getVel( vel ); for (j=0; j < 3; j++) oldVel[3*i + j] = vel[j]; - if( atoms[i]->isDirectional() ){ + if( integrableObjects[i]->isDirectional() ){ - dAtom = (DirectionalAtom *)atoms[i]; + integrableObjects[i]->getJ( ji ); - dAtom->getJ( ji ); - for (j=0; j < 3; j++) oldJi[3*i + j] = ji[j]; @@ -188,9 +211,9 @@ template void NPT::moveB( void ){ // do the iteration: instaVol = tStats->getVolume(); - + for (k=0; k < 4; k++) { - + instaTemp = tStats->getTemperature(); instaPress = tStats->getPressure(); @@ -198,42 +221,41 @@ template void NPT::moveB( void ){ this->evolveChiB(); this->evolveEtaB(); - - for( i=0; icalcVelScale(); - atoms[i]->getFrc( frc ); - atoms[i]->getVel(vel); - - mass = atoms[i]->getMass(); - + for( i=0; igetFrc( frc ); + integrableObjects[i]->getVel(vel); + + mass = integrableObjects[i]->getMass(); + getVelScaleB( sc, i ); // velocity half step - for (j=0; j < 3; j++) + for (j=0; j < 3; j++) vel[j] = oldVel[3*i+j] + dt2 * ((frc[j] / mass ) * eConvert - sc[j]); - - atoms[i]->setVel( vel ); - - if( atoms[i]->isDirectional() ){ - dAtom = (DirectionalAtom *)atoms[i]; - - // get and convert the torque to body frame - - dAtom->getTrq( Tb ); - dAtom->lab2Body( Tb ); - - for (j=0; j < 3; j++) + integrableObjects[i]->setVel( vel ); + + if( integrableObjects[i]->isDirectional() ){ + + // get and convert the torque to body frame + + integrableObjects[i]->getTrq( Tb ); + integrableObjects[i]->lab2Body( Tb ); + + for (j=0; j < 3; j++) ji[j] = oldJi[3*i + j] + dt2 * (Tb[j] * eConvert - oldJi[3*i+j]*chi); - - dAtom->setJ( ji ); + + integrableObjects[i]->setJ( ji ); } } - + if (nConstrained){ constrainB(); - } - + } + if ( this->chiConverged() && this->etaConverged() ) break; } @@ -254,14 +276,14 @@ template void NPT::evolveChiB() { } template void NPT::evolveChiB() { - + prevChi = chi; chi = oldChi + dt2 * ( instaTemp / targetTemp - 1.0) / tt2; } template bool NPT::chiConverged() { - - return ( fabs( prevChi - chi ) <= chiTolerance ); + + return ( fabs( prevChi - chi ) <= chiTolerance ); } template int NPT::readyCheck() { @@ -269,10 +291,10 @@ template int NPT::readyCheck() { //check parent's readyCheck() first if (T::readyCheck() == -1) return -1; - - // First check to see if we have a target temperature. - // Not having one is fatal. - + + // First check to see if we have a target temperature. + // Not having one is fatal. + if (!have_target_temp) { sprintf( painCave.errMsg, "NPT error: You can't use the NPT integrator\n" @@ -292,9 +314,9 @@ template int NPT::readyCheck() { simError(); return -1; } - + // We must set tauThermostat. - + if (!have_tau_thermostat) { sprintf( painCave.errMsg, "NPT error: If you use the NPT\n" @@ -302,10 +324,10 @@ template int NPT::readyCheck() { painCave.isFatal = 1; simError(); return -1; - } + } // We must set tauBarostat. - + if (!have_tau_barostat) { sprintf( painCave.errMsg, "NPT error: If you use the NPT\n" @@ -313,7 +335,7 @@ template int NPT::readyCheck() { painCave.isFatal = 1; simError(); return -1; - } + } if (!have_chi_tolerance) { sprintf( painCave.errMsg, @@ -322,7 +344,7 @@ template int NPT::readyCheck() { have_chi_tolerance = 1; painCave.isFatal = 0; simError(); - } + } if (!have_eta_tolerance) { sprintf( painCave.errMsg, @@ -331,18 +353,18 @@ template int NPT::readyCheck() { have_eta_tolerance = 1; painCave.isFatal = 0; simError(); - } - + } + // We need NkBT a lot, so just set it here: This is the RAW number // of particles, so no subtraction or addition of constraints or // orientational degrees of freedom: - + NkBT = (double)Nparticles * kB * targetTemp; - + // fkBT is used because the thermostat operates on more degrees of freedom // than the barostat (when there are particles with orientational degrees // of freedom). ndf = 3 * (n_atoms + n_oriented -1) - n_constraint - nZcons - + fkBT = (double)info->ndf * kB * targetTemp; tt2 = tauThermostat * tauThermostat;