--- trunk/OOPSE/libmdtools/NPTxyz.cpp 2003/10/28 16:03:37 829 +++ trunk/OOPSE/libmdtools/NPTxyz.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,45 @@ // 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 NPTxyz::NPTxyz ( 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 NPTxyz::~NPTxyz() { @@ -44,37 +64,37 @@ template void NPTxyz::resetIntegrator() } template void NPTxyz::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 NPTxyz::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] = 0.0; } } - + for(i = 0; i < 3; i++) for (j = 0; j < 3; j++) oldEta[i][j] = eta[i][j]; } template void NPTxyz::evolveEtaB() { - + int i,j; for(i = 0; i < 3; i++) @@ -84,7 +104,7 @@ template void NPTxyz::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] = 0.0; @@ -100,13 +120,13 @@ template void NPTxyz::getVelScaleA(doub 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 +138,20 @@ template void NPTxyz::getVelScaleB(doub 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 NPTxyz::getPosScale(double pos[3], double COM[3], +template void NPTxyz::getPosScale(double pos[3], double COM[3], int index, double sc[3]){ int j; double rj[3]; @@ -149,31 +169,31 @@ template void NPTxyz::scaleSimBox( void double eta2ij, scaleFactor; 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++){ scaleMat[i][j] = 0.0; if(i==j) scaleMat[i][j] = 1.0; } } - + for(i=0;i<3;i++){ - + // calculate the scaleFactors - + scaleFactor = exp(dt*eta[i][i]); - + scaleMat[i][i] = scaleFactor; if (scaleMat[i][i] > bigScale) bigScale = scaleMat[i][i]; @@ -182,15 +202,15 @@ template void NPTxyz::scaleSimBox( void // 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; @@ -198,14 +218,14 @@ template void NPTxyz::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, "NPTxyz error: Attempting a Box scaling of more than 10 percent.\n" @@ -231,15 +251,15 @@ template bool NPTxyz::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 NPTxyz::getConservedQuantity(void){ - + double conservedQuantity; double totalEnergy; double thermostat_kinetic; @@ -251,7 +271,7 @@ template double NPTxyz::getConservedQua totalEnergy = tStats->getTotalE(); - thermostat_kinetic = fkBT * tt2 * chi * chi / + thermostat_kinetic = fkBT * tt2 * chi * chi / (2.0 * eConvert); thermostat_potential = fkBT* integralOfChidt / eConvert; @@ -260,22 +280,39 @@ template double NPTxyz::getConservedQua 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 NPTxyz::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; + +}