OOPSE/libmdtools/NPTf.cpp

#include <stdlib.h>
#include <math.h>
#include <string.h>

#include "Atom.hpp"
#include "SRI.hpp"
#include "AbstractClasses.hpp"
#include "SimInfo.hpp"
#include "ForceFields.hpp"
#include "Thermo.hpp"
#include "ReadWrite.hpp"
#include "Integrator.hpp"
#include "simError.h"

#ifdef IS_MPI
#include "mpiSimulation.hpp"
#endif

// Basic non-isotropic thermostating and barostating via the Melchionna
// modification of the Hoover algorithm:
//
//    Melchionna, S., Ciccotti, G., and Holian, B. L., 1993,
//       Molec. Phys., 78, 533.
//
//           and
//
//    Hoover, W. G., 1986, Phys. Rev. A, 34, 2499.

NPTf::NPTf ( SimInfo *theInfo, ForceFields* the_ff):
  NPT( theInfo, the_ff )
{
  GenericData* data;
  double *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 != NULL){
    
    int test = data->getDarray(etaArray);
    
    if( test == 9 ){
      
      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];
        }
      }    
      delete[] etaArray;
    }
    else
      std::cerr << "NPTf error: etaArray is not length 9 (actual = " << test
                << ").\n"
                << "            Simulation wil proceed with eta = 0;\n";
  }
}

NPTf::~NPTf() {

  // empty for now
}

void NPTf::resetIntegrator() {

  int i, j;

  for(i = 0; i < 3; i++)
    for (j = 0; j < 3; j++)
      eta[i][j] = 0.0;

  NPT::resetIntegrator();
}

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 *
          (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];
}

void NPTf::evolveEtaB() {

  int i,j;

  for(i = 0; i < 3; i++)
    for (j = 0; j < 3; j++)
      prevEta[i][j] = eta[i][j];

  for(i = 0; i < 3; i ++){
    for(j = 0; j < 3; j++){
      if( i == j) {
        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);
      }
    }
  }
}

void NPTf::getVelScaleA(double sc[3], double vel[3]) {
  int i,j;
  double vScale[3][3];

  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;
      }
    }
  }

  info->matVecMul3( vScale, vel, sc );
}

void NPTf::getVelScaleB(double sc[3], int index ){
  int i,j;
  double myVel[3];
  double vScale[3][3];

//   std::cerr << "velScaleB chi = " << chi << "\n";

  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;
      }
    }
  }

  for (j = 0; j < 3; j++)
    myVel[j] = oldVel[3*index + j];

//   std::cerr << "velScaleB = \n" 
//          << "[ " << vScale[0][0] << " , " << vScale[0][1] << " , " << vScale[0][2] << "]\n"
//          << "[ " << vScale[1][0] << " , " << vScale[1][1] << " , " << vScale[1][2] << "]\n"
//          << "[ " << vScale[2][0] << " , " << vScale[2][1] << " , " << vScale[2][2] << "]\n\n";


 //  std::cerr << "myVel " << index << " in => " 
//          << myVel[0] << ", " << myVel[1] << ", " << myVel[2] << "\n";

  info->matVecMul3( vScale, myVel, sc );

 //  std::cerr << "sc " << index << " out => " 
//          << sc[0] << ", " << sc[1] << ", " << sc[2] << "\n";
}

void NPTf::getPosScale(double pos[3], double COM[3],
                                               int index, double sc[3]){
  int j;
  double rj[3];

  for(j=0; j<3; j++)
    rj[j] = ( oldPos[index*3+j] + pos[j]) / 2.0 - COM[j];

  info->matVecMul3( eta, rj, sc );
}

void NPTf::scaleSimBox( void ){

  int i,j,k;
  double scaleMat[3][3];
  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;
      // Taylor expansion for the exponential truncated at second order:
      scaleMat[i][j] += dt*eta[i][j]  + 0.5*dt*dt*eta2ij;

      if (i != j)
        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"
             " Check your tauBarostat, as it is probably too small!\n\n"
             " scaleMat = [%lf\t%lf\t%lf]\n"
             "            [%lf\t%lf\t%lf]\n"
             "            [%lf\t%lf\t%lf]\n",
             scaleMat[0][0],scaleMat[0][1],scaleMat[0][2],
             scaleMat[1][0],scaleMat[1][1],scaleMat[1][2],
             scaleMat[2][0],scaleMat[2][1],scaleMat[2][2]);
    painCave.isFatal = 1;
    simError();
  } else if (offDiagMax > 0.1) {
    sprintf( painCave.errMsg,
             "NPTf error: Attempting an off-diagonal Box scaling of more than 10 percent.\n"
             " Check your tauBarostat, as it is probably too small!\n\n"
             " scaleMat = [%lf\t%lf\t%lf]\n"
             "            [%lf\t%lf\t%lf]\n"
             "            [%lf\t%lf\t%lf]\n",
             scaleMat[0][0],scaleMat[0][1],scaleMat[0][2],
             scaleMat[1][0],scaleMat[1][1],scaleMat[1][2],
             scaleMat[2][0],scaleMat[2][1],scaleMat[2][2]);
    painCave.isFatal = 1;
    simError();
  } else {
    info->getBoxM(hm);
    info->matMul3(hm, scaleMat, hmnew);
    info->setBoxM(hmnew);
  }
}

bool NPTf::etaConverged() {
  int i;
  double diffEta, sumEta;

  sumEta = 0;
  for(i = 0; i < 3; i++)
    sumEta += pow(prevEta[i][i] - eta[i][i], 2);

  diffEta = sqrt( sumEta / 3.0 );

  return ( diffEta <= etaTolerance );
}

double NPTf::getConservedQuantity(void){

  double conservedQuantity;
  double totalEnergy;
  double thermostat_kinetic;
  double thermostat_potential;
  double barostat_kinetic;
  double barostat_potential;
  double trEta;
  double a[3][3], b[3][3];

  totalEnergy = tStats->getTotalE();

  thermostat_kinetic = fkBT * tt2 * chi * chi /
    (2.0 * eConvert);

  thermostat_potential = fkBT* integralOfChidt / eConvert;

  info->transposeMat3(eta, a);
  info->matMul3(a, eta, b);
  trEta = info->matTrace3(b);

  barostat_kinetic = NkBT * tb2 * trEta /
    (2.0 * eConvert);

  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 <<
//       "\t" << barostat_potential << "\t" << conservedQuantity << endl;

  return conservedQuantity;

}

char* NPTf::getAdditionalParameters(void){

  sprintf(addParamBuffer,
          "\t%G\t%G;"
          "\t%G\t%G\t%G;"
          "\t%G\t%G\t%G;"
          "\t%G\t%G\t%G;",
          chi, integralOfChidt,
          eta[0][0], eta[0][1], eta[0][2],
          eta[1][0], eta[1][1], eta[1][2],
          eta[2][0], eta[2][1], eta[2][2]
          );

  return addParamBuffer;
}
Revision:	851
Committed:	Wed Nov 5 19:18:17 2003 UTC (22 years ago) by mmeineke
File size:	7942 byte(s)
Log Message:	some work on trying to find the compression bug
#	User	Rev	Content
1	mmeineke	850	#include <stdlib.h>
2	gezelter	829	#include <math.h>
3	mmeineke	851	#include <string.h>
4	mmeineke	850
5	gezelter	576	#include "Atom.hpp"
6			#include "SRI.hpp"
7			#include "AbstractClasses.hpp"
8			#include "SimInfo.hpp"
9			#include "ForceFields.hpp"
10			#include "Thermo.hpp"
11			#include "ReadWrite.hpp"
12			#include "Integrator.hpp"
13	tim	837	#include "simError.h"
14	gezelter	576
15	gezelter	772	#ifdef IS_MPI
16			#include "mpiSimulation.hpp"
17			#endif
18	gezelter	576
19	gezelter	578	// Basic non-isotropic thermostating and barostating via the Melchionna
20	gezelter	576	// modification of the Hoover algorithm:
21			//
22			// Melchionna, S., Ciccotti, G., and Holian, B. L., 1993,
23	tim	837	// Molec. Phys., 78, 533.
24	gezelter	576	//
25			// and
26	tim	837	//
27	gezelter	576	// Hoover, W. G., 1986, Phys. Rev. A, 34, 2499.
28
29	mmeineke	849	NPTf::NPTf ( SimInfo theInfo, ForceFields the_ff):
30	mmeineke	850	NPT( theInfo, the_ff )
31	gezelter	576	{
32	tim	837	GenericData* data;
33	mmeineke	850	double *etaArray;
34	mmeineke	780	int i,j;
35	tim	837
36	mmeineke	780	for(i = 0; i < 3; i++){
37			for (j = 0; j < 3; j++){
38	tim	837
39	gezelter	588	eta[i][j] = 0.0;
40	mmeineke	780	oldEta[i][j] = 0.0;
41			}
42			}
43	tim	837
44	mmeineke	850	// retrieve eta array from simInfo if it exists
45			data = info->getProperty(ETAVALUE_ID);
46			if(data != NULL){
47
48			int test = data->getDarray(etaArray);
49
50			if( test == 9 ){
51
52			for(i = 0; i < 3; i++){
53			for (j = 0; j < 3; j++){
54			eta[i][j] = etaArray[3*i+j];
55			oldEta[i][j] = eta[i][j];
56			}
57			}
58			delete[] etaArray;
59	tim	837	}
60	mmeineke	850	else
61			std::cerr << "NPTf error: etaArray is not length 9 (actual = " << test
62			<< ").\n"
63			<< " Simulation wil proceed with eta = 0;\n";
64			}
65	mmeineke	780	}
66	gezelter	588
67	mmeineke	849	NPTf::~NPTf() {
68	tim	767
69	mmeineke	780	// empty for now
70			}
71	tim	767
72	mmeineke	849	void NPTf::resetIntegrator() {
73	tim	837
74	mmeineke	780	int i, j;
75	tim	837
76	mmeineke	780	for(i = 0; i < 3; i++)
77			for (j = 0; j < 3; j++)
78			eta[i][j] = 0.0;
79	tim	837
80	mmeineke	850	NPT::resetIntegrator();
81	gezelter	576	}
82
83	mmeineke	849	void NPTf::evolveEtaA() {
84	tim	837
85	mmeineke	780	int i, j;
86	tim	837
87	mmeineke	780	for(i = 0; i < 3; i ++){
88			for(j = 0; j < 3; j++){
89			if( i == j)
90	tim	837	eta[i][j] += dt2 * instaVol *
91	mmeineke	780	(press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
92			else
93			eta[i][j] += dt2 * instaVol * press[i][j] / (NkBT*tb2);
94			}
95			}
96	tim	837
97	mmeineke	780	for(i = 0; i < 3; i++)
98			for (j = 0; j < 3; j++)
99			oldEta[i][j] = eta[i][j];
100	tim	767	}
101
102	mmeineke	849	void NPTf::evolveEtaB() {
103	tim	837
104	mmeineke	780	int i,j;
105	gezelter	772
106	mmeineke	780	for(i = 0; i < 3; i++)
107			for (j = 0; j < 3; j++)
108			prevEta[i][j] = eta[i][j];
109	mmeineke	778
110	mmeineke	780	for(i = 0; i < 3; i ++){
111			for(j = 0; j < 3; j++){
112			if( i == j) {
113	tim	837	eta[i][j] = oldEta[i][j] + dt2 * instaVol *
114	mmeineke	780	(press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
115			} else {
116			eta[i][j] = oldEta[i][j] + dt2 * instaVol * press[i][j] / (NkBT*tb2);
117			}
118			}
119			}
120			}
121	gezelter	600
122	mmeineke	849	void NPTf::getVelScaleA(double sc[3], double vel[3]) {
123	mmeineke	780	int i,j;
124			double vScale[3][3];
125	gezelter	576
126	gezelter	588	for (i = 0; i < 3; i++ ) {
127			for (j = 0; j < 3; j++ ) {
128	tim	767	vScale[i][j] = eta[i][j];
129	tim	837
130	gezelter	588	if (i == j) {
131	tim	837	vScale[i][j] += chi;
132			}
133	gezelter	588	}
134			}
135	tim	837
136	mmeineke	780	info->matVecMul3( vScale, vel, sc );
137			}
138	gezelter	600
139	mmeineke	849	void NPTf::getVelScaleB(double sc[3], int index ){
140	mmeineke	780	int i,j;
141			double myVel[3];
142			double vScale[3][3];
143	gezelter	600
144	mmeineke	851	// std::cerr << "velScaleB chi = " << chi << "\n";
145
146	mmeineke	780	for (i = 0; i < 3; i++ ) {
147			for (j = 0; j < 3; j++ ) {
148			vScale[i][j] = eta[i][j];
149	tim	837
150	mmeineke	780	if (i == j) {
151	tim	837	vScale[i][j] += chi;
152			}
153	tim	767	}
154			}
155	tim	837
156	mmeineke	780	for (j = 0; j < 3; j++)
157			myVel[j] = oldVel[3*index + j];
158	tim	767
159	mmeineke	851	// std::cerr << "velScaleB = \n"
160			// << "[ " << vScale[0][0] << " , " << vScale[0][1] << " , " << vScale[0][2] << "]\n"
161			// << "[ " << vScale[1][0] << " , " << vScale[1][1] << " , " << vScale[1][2] << "]\n"
162			// << "[ " << vScale[2][0] << " , " << vScale[2][1] << " , " << vScale[2][2] << "]\n\n";
163
164
165			// std::cerr << "myVel " << index << " in => "
166			// << myVel[0] << ", " << myVel[1] << ", " << myVel[2] << "\n";
167
168	mmeineke	780	info->matVecMul3( vScale, myVel, sc );
169	mmeineke	851
170			// std::cerr << "sc " << index << " out => "
171			// << sc[0] << ", " << sc[1] << ", " << sc[2] << "\n";
172	mmeineke	780	}
173	tim	767
174	mmeineke	849	void NPTf::getPosScale(double pos[3], double COM[3],
175	mmeineke	780	int index, double sc[3]){
176			int j;
177			double rj[3];
178	tim	767
179	mmeineke	780	for(j=0; j<3; j++)
180			rj[j] = ( oldPos[index*3+j] + pos[j]) / 2.0 - COM[j];
181	tim	767
182	mmeineke	780	info->matVecMul3( eta, rj, sc );
183			}
184	tim	767
185	mmeineke	849	void NPTf::scaleSimBox( void ){
186	tim	767
187	mmeineke	780	int i,j,k;
188			double scaleMat[3][3];
189			double eta2ij;
190			double bigScale, smallScale, offDiagMax;
191			double hm[3][3], hmnew[3][3];
192	tim	767
193	mmeineke	780
194	tim	837
195	gezelter	577	// Scale the box after all the positions have been moved:
196	tim	837
197	gezelter	578	// Use a taylor expansion for eta products: Hmat = Hmat . exp(dt * etaMat)
198			// Hmat = Hmat . ( Ident + dt * etaMat + dt^2 * etaMat*etaMat / 2)
199	tim	837
200	gezelter	617	bigScale = 1.0;
201			smallScale = 1.0;
202			offDiagMax = 0.0;
203	tim	837
204	gezelter	578	for(i=0; i<3; i++){
205			for(j=0; j<3; j++){
206	tim	837
207	gezelter	588	// Calculate the matrix Product of the eta array (we only need
208			// the ij element right now):
209	tim	837
210	gezelter	588	eta2ij = 0.0;
211	gezelter	578	for(k=0; k<3; k++){
212	gezelter	588	eta2ij += eta[i][k] * eta[k][j];
213	gezelter	578	}
214	tim	837
215	gezelter	588	scaleMat[i][j] = 0.0;
216			// identity matrix (see above):
217			if (i == j) scaleMat[i][j] = 1.0;
218			// Taylor expansion for the exponential truncated at second order:
219			scaleMat[i][j] += dteta[i][j] + 0.5dtdteta2ij;
220	gezelter	617
221			if (i != j)
222	tim	837	if (fabs(scaleMat[i][j]) > offDiagMax)
223	gezelter	617	offDiagMax = fabs(scaleMat[i][j]);
224	gezelter	578	}
225	gezelter	617
226			if (scaleMat[i][i] > bigScale) bigScale = scaleMat[i][i];
227			if (scaleMat[i][i] < smallScale) smallScale = scaleMat[i][i];
228	gezelter	578	}
229	tim	837
230	gezelter	617	if ((bigScale > 1.1) \|\| (smallScale < 0.9)) {
231			sprintf( painCave.errMsg,
232			"NPTf error: Attempting a Box scaling of more than 10 percent.\n"
233			" Check your tauBarostat, as it is probably too small!\n\n"
234			" scaleMat = [%lf\t%lf\t%lf]\n"
235			" [%lf\t%lf\t%lf]\n"
236			" [%lf\t%lf\t%lf]\n",
237			scaleMat[0][0],scaleMat[0][1],scaleMat[0][2],
238			scaleMat[1][0],scaleMat[1][1],scaleMat[1][2],
239			scaleMat[2][0],scaleMat[2][1],scaleMat[2][2]);
240			painCave.isFatal = 1;
241			simError();
242			} else if (offDiagMax > 0.1) {
243			sprintf( painCave.errMsg,
244			"NPTf error: Attempting an off-diagonal Box scaling of more than 10 percent.\n"
245			" Check your tauBarostat, as it is probably too small!\n\n"
246			" scaleMat = [%lf\t%lf\t%lf]\n"
247			" [%lf\t%lf\t%lf]\n"
248			" [%lf\t%lf\t%lf]\n",
249			scaleMat[0][0],scaleMat[0][1],scaleMat[0][2],
250			scaleMat[1][0],scaleMat[1][1],scaleMat[1][2],
251			scaleMat[2][0],scaleMat[2][1],scaleMat[2][2]);
252			painCave.isFatal = 1;
253			simError();
254			} else {
255			info->getBoxM(hm);
256			info->matMul3(hm, scaleMat, hmnew);
257			info->setBoxM(hmnew);
258			}
259	gezelter	576	}
260
261	mmeineke	849	bool NPTf::etaConverged() {
262	mmeineke	780	int i;
263			double diffEta, sumEta;
264	gezelter	600
265	mmeineke	780	sumEta = 0;
266	tim	767	for(i = 0; i < 3; i++)
267	tim	837	sumEta += pow(prevEta[i][i] - eta[i][i], 2);
268
269	mmeineke	780	diffEta = sqrt( sumEta / 3.0 );
270	tim	837
271	mmeineke	780	return ( diffEta <= etaTolerance );
272	gezelter	576	}
273
274	mmeineke	849	double NPTf::getConservedQuantity(void){
275	tim	837
276	tim	763	double conservedQuantity;
277	mmeineke	782	double totalEnergy;
278	gezelter	772	double thermostat_kinetic;
279			double thermostat_potential;
280			double barostat_kinetic;
281			double barostat_potential;
282			double trEta;
283			double a[3][3], b[3][3];
284	tim	763
285	mmeineke	782	totalEnergy = tStats->getTotalE();
286	tim	763
287	tim	837	thermostat_kinetic = fkBT * tt2 * chi * chi /
288	gezelter	772	(2.0 * eConvert);
289	tim	763
290	gezelter	772	thermostat_potential = fkBT* integralOfChidt / eConvert;
291	tim	763
292	gezelter	772	info->transposeMat3(eta, a);
293			info->matMul3(a, eta, b);
294			trEta = info->matTrace3(b);
295	tim	767
296	tim	837	barostat_kinetic = NkBT * tb2 * trEta /
297	gezelter	772	(2.0 * eConvert);
298	tim	837
299			barostat_potential = (targetPressure * tStats->getVolume() / p_convert) /
300	gezelter	772	eConvert;
301	tim	767
302	mmeineke	782	conservedQuantity = totalEnergy + thermostat_kinetic + thermostat_potential +
303	gezelter	772	barostat_kinetic + barostat_potential;
304	tim	837
305	mmeineke	780	// cout.width(8);
306			// cout.precision(8);
307	tim	767
308	tim	837	// cerr << info->getTime() << "\t" << Energy << "\t" << thermostat_kinetic <<
309			// "\t" << thermostat_potential << "\t" << barostat_kinetic <<
310	mmeineke	780	// "\t" << barostat_potential << "\t" << conservedQuantity << endl;
311	gezelter	772
312	tim	837	return conservedQuantity;
313
314	tim	763	}
315	tim	837
316	mmeineke	851	char* NPTf::getAdditionalParameters(void){
317	tim	837
318	mmeineke	851	sprintf(addParamBuffer,
319			"\t%G\t%G;"
320			"\t%G\t%G\t%G;"
321			"\t%G\t%G\t%G;"
322			"\t%G\t%G\t%G;",
323			chi, integralOfChidt,
324			eta[0][0], eta[0][1], eta[0][2],
325			eta[1][0], eta[1][1], eta[1][2],
326			eta[2][0], eta[2][1], eta[2][2]
327			);
328	tim	837
329	mmeineke	851	return addParamBuffer;
330	tim	837	}