--- trunk/OOPSE/libmdtools/NPTf.cpp 2003/10/28 16:03:37 829 +++ trunk/OOPSE/libmdtools/NPTf.cpp 2003/10/29 00:19:10 837 @@ -7,7 +7,7 @@ #include "Thermo.hpp" #include "ReadWrite.hpp" #include "Integrator.hpp" -#include "simError.h" +#include "simError.h" #ifdef IS_MPI #include "mpiSimulation.hpp" @@ -17,25 +17,46 @@ // 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 NPTf::NPTf ( SimInfo *theInfo, ForceFields* the_ff): T( theInfo, the_ff ) { - + GenericData* data; + DoubleArrayData * etaValue; + vector etaArray; int i,j; - + for(i = 0; i < 3; i++){ for (j = 0; j < 3; j++){ - + eta[i][j] = 0.0; oldEta[i][j] = 0.0; } } + + // retrieve eta array from simInfo if it exists + data = info->getProperty(ETAVALUE_ID); + if(data){ + etaValue = dynamic_cast(data); + + if(etaValue){ + etaArray = etaValue->getData(); + + for(i = 0; i < 3; i++){ + for (j = 0; j < 3; j++){ + eta[i][j] = etaArray[3*i+j]; + oldEta[i][j] = eta[i][j]; + } + } + + } + } + } template NPTf::~NPTf() { @@ -44,37 +65,37 @@ template void NPTf::resetIntegrator() { } template void NPTf::resetIntegrator() { - + int i, j; - + for(i = 0; i < 3; i++) for (j = 0; j < 3; j++) eta[i][j] = 0.0; - + T::resetIntegrator(); } template void NPTf::evolveEtaA() { - + int i, j; - + for(i = 0; i < 3; i ++){ for(j = 0; j < 3; j++){ if( i == j) - eta[i][j] += dt2 * instaVol * + eta[i][j] += dt2 * instaVol * (press[i][j] - targetPressure/p_convert) / (NkBT*tb2); else eta[i][j] += dt2 * instaVol * press[i][j] / (NkBT*tb2); } } - + for(i = 0; i < 3; i++) for (j = 0; j < 3; j++) oldEta[i][j] = eta[i][j]; } template void NPTf::evolveEtaB() { - + int i,j; for(i = 0; i < 3; i++) @@ -84,7 +105,7 @@ template void NPTf::evolveEtaB() { for(i = 0; i < 3; i ++){ for(j = 0; j < 3; j++){ if( i == j) { - eta[i][j] = oldEta[i][j] + dt2 * instaVol * + eta[i][j] = oldEta[i][j] + dt2 * instaVol * (press[i][j] - targetPressure/p_convert) / (NkBT*tb2); } else { eta[i][j] = oldEta[i][j] + dt2 * instaVol * press[i][j] / (NkBT*tb2); @@ -100,13 +121,13 @@ template void NPTf::getVelScaleA(double for (i = 0; i < 3; i++ ) { for (j = 0; j < 3; j++ ) { vScale[i][j] = eta[i][j]; - + if (i == j) { - vScale[i][j] += chi; - } + vScale[i][j] += chi; + } } } - + info->matVecMul3( vScale, vel, sc ); } @@ -118,20 +139,20 @@ template void NPTf::getVelScaleB(double for (i = 0; i < 3; i++ ) { for (j = 0; j < 3; j++ ) { vScale[i][j] = eta[i][j]; - + if (i == j) { - vScale[i][j] += chi; - } + vScale[i][j] += chi; + } } } - + for (j = 0; j < 3; j++) myVel[j] = oldVel[3*index + j]; info->matVecMul3( vScale, myVel, sc ); } -template void NPTf::getPosScale(double pos[3], double COM[3], +template void NPTf::getPosScale(double pos[3], double COM[3], int index, double sc[3]){ int j; double rj[3]; @@ -149,29 +170,29 @@ template void NPTf::scaleSimBox( void ) double eta2ij; double bigScale, smallScale, offDiagMax; double hm[3][3], hmnew[3][3]; - + // Scale the box after all the positions have been moved: - + // Use a taylor expansion for eta products: Hmat = Hmat . exp(dt * etaMat) // Hmat = Hmat . ( Ident + dt * etaMat + dt^2 * etaMat*etaMat / 2) - + bigScale = 1.0; smallScale = 1.0; offDiagMax = 0.0; - + for(i=0; i<3; i++){ for(j=0; j<3; j++){ - + // Calculate the matrix Product of the eta array (we only need // the ij element right now): - + eta2ij = 0.0; for(k=0; k<3; k++){ eta2ij += eta[i][k] * eta[k][j]; } - + scaleMat[i][j] = 0.0; // identity matrix (see above): if (i == j) scaleMat[i][j] = 1.0; @@ -179,14 +200,14 @@ template void NPTf::scaleSimBox( void ) scaleMat[i][j] += dt*eta[i][j] + 0.5*dt*dt*eta2ij; if (i != j) - if (fabs(scaleMat[i][j]) > offDiagMax) + if (fabs(scaleMat[i][j]) > offDiagMax) offDiagMax = fabs(scaleMat[i][j]); } if (scaleMat[i][i] > bigScale) bigScale = scaleMat[i][i]; if (scaleMat[i][i] < smallScale) smallScale = scaleMat[i][i]; } - + if ((bigScale > 1.1) || (smallScale < 0.9)) { sprintf( painCave.errMsg, "NPTf error: Attempting a Box scaling of more than 10 percent.\n" @@ -224,15 +245,15 @@ template bool NPTf::etaConverged() { sumEta = 0; for(i = 0; i < 3; i++) - sumEta += pow(prevEta[i][i] - eta[i][i], 2); - + sumEta += pow(prevEta[i][i] - eta[i][i], 2); + diffEta = sqrt( sumEta / 3.0 ); - + return ( diffEta <= etaTolerance ); } template double NPTf::getConservedQuantity(void){ - + double conservedQuantity; double totalEnergy; double thermostat_kinetic; @@ -244,7 +265,7 @@ template double NPTf::getConservedQuant totalEnergy = tStats->getTotalE(); - thermostat_kinetic = fkBT * tt2 * chi * chi / + thermostat_kinetic = fkBT * tt2 * chi * chi / (2.0 * eConvert); thermostat_potential = fkBT* integralOfChidt / eConvert; @@ -253,22 +274,39 @@ template double NPTf::getConservedQuant info->matMul3(a, eta, b); trEta = info->matTrace3(b); - barostat_kinetic = NkBT * tb2 * trEta / + barostat_kinetic = NkBT * tb2 * trEta / (2.0 * eConvert); - - barostat_potential = (targetPressure * tStats->getVolume() / p_convert) / + + barostat_potential = (targetPressure * tStats->getVolume() / p_convert) / eConvert; conservedQuantity = totalEnergy + thermostat_kinetic + thermostat_potential + barostat_kinetic + barostat_potential; - + // cout.width(8); // cout.precision(8); -// cerr << info->getTime() << "\t" << Energy << "\t" << thermostat_kinetic << -// "\t" << thermostat_potential << "\t" << barostat_kinetic << +// cerr << info->getTime() << "\t" << Energy << "\t" << thermostat_kinetic << +// "\t" << thermostat_potential << "\t" << barostat_kinetic << // "\t" << barostat_potential << "\t" << conservedQuantity << endl; - return conservedQuantity; - + return conservedQuantity; + } + +template string NPTf::getAdditionalParameters(void){ + string parameters; + const int BUFFERSIZE = 2000; // size of the read buffer + char buffer[BUFFERSIZE]; + + sprintf(buffer,"\t%f\t%f;", chi, integralOfChidt); + parameters += buffer; + + for(int i = 0; i < 3; i++){ + sprintf(buffer,"\t%g\t%g\t%g;", eta[3*i], eta[3*i+1], eta[3*i+2]); + parameters += buffer; + } + + return parameters; + +}