OOPSE/libmdtools/NPTi.cpp

#include <iostream>
#include <math.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 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.

NPTi::NPTi ( SimInfo *theInfo, ForceFields* the_ff):
  NPT( theInfo, the_ff )
{
  GenericData* data;
  double *etaArray;
  int test;

  eta = 0.0;
  oldEta = 0.0;

  // retrieve eta array from simInfo if it exists
  data = info->getProperty(ETAVALUE_ID);
  if(data != NULL){
    
    test = data->getDarray(etaArray);
    
    if( test == 9 ){
      
      eta = etaArray[0];
      delete[] etaArray;
    }
    else
      std::cerr << "NPTi error: etaArray is not length 9 (actual = " << test
                << ").\n"
                << "            Simulation wil proceed with eta = 0;\n";
  }
}

NPTi::~NPTi() {
  //nothing for now
}

void NPTi::resetIntegrator() {
  eta = 0.0;
  NPT::resetIntegrator();
}

void NPTi::evolveEtaA() {
  eta += dt2 * ( instaVol * (instaPress - targetPressure) /
                 (p_convert*NkBT*tb2));
  oldEta = eta;
}

void NPTi::evolveEtaB() {

  prevEta = eta;
  eta = oldEta + dt2 * ( instaVol * (instaPress - targetPressure) /
                 (p_convert*NkBT*tb2));
}

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

  for(i=0; i<3; i++) sc[i] = vel[i] * ( chi + eta );
}

void NPTi::getVelScaleB(double sc[3], int index ){
  int i;

  for(i=0; i<3; i++) sc[i] = oldVel[index*3 + i] * ( chi + eta );
}


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

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

  for(j=0; j<3; j++)
    sc[j] *= eta;
}

void NPTi::scaleSimBox( void ){

  double scaleFactor;

  scaleFactor = exp(dt*eta);

  if ((scaleFactor > 1.1) || (scaleFactor < 0.9)) {
    sprintf( painCave.errMsg,
             "NPTi error: Attempting a Box scaling of more than 10 percent"
             " check your tauBarostat, as it is probably too small!\n"
             " eta = %lf, scaleFactor = %lf\n", eta, scaleFactor
             );
    painCave.isFatal = 1;
    simError();
  } else {
    info->scaleBox(scaleFactor);
  }

}

bool NPTi::etaConverged() {

  return ( fabs(prevEta - eta) <= etaTolerance );
}

double NPTi::getConservedQuantity(void){

  double conservedQuantity;
  double Energy;
  double thermostat_kinetic;
  double thermostat_potential;
  double barostat_kinetic;
  double barostat_potential;

  Energy = tStats->getTotalE();

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

  thermostat_potential = fkBT* integralOfChidt / eConvert;


  barostat_kinetic = 3.0 * NkBT * tb2 * eta * eta /
    (2.0 * eConvert);

  barostat_potential = (targetPressure * tStats->getVolume() / p_convert) /
    eConvert;

  conservedQuantity = Energy + 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* NPTi::getAdditionalParameters(void){

  sprintf(addParamBuffer,
          "\t%G\t%G;"
          "\t%G\t%0.0\t%0.0;"
          "\t%0.0\t%G\t%0.0;"
          "\t%0.0\t%0.0\t%G;",
          chi, integralOfChidt,
          eta, eta, eta
          );

  return addParamBuffer;
}
Revision:	851
Committed:	Wed Nov 5 19:18:17 2003 UTC (21 years, 8 months ago) by mmeineke
File size:	3772 byte(s)
Log Message:	some work on trying to find the compression bug
#	User	Rev	Content
1	mmeineke	851	#include <iostream>
2	gezelter	829	#include <math.h>
3	mmeineke	851
4
5	gezelter	574	#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	574
15	tim	763	#ifdef IS_MPI
16			#include "mpiSimulation.hpp"
17			#endif
18	gezelter	574
19			// Basic isotropic thermostating and barostating via the Melchionna
20			// 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	574	//
25			// and
26	tim	837	//
27	gezelter	574	// Hoover, W. G., 1986, Phys. Rev. A, 34, 2499.
28
29	mmeineke	849	NPTi::NPTi ( SimInfo theInfo, ForceFields the_ff):
30	mmeineke	850	NPT( theInfo, the_ff )
31	gezelter	574	{
32	tim	837	GenericData* data;
33	mmeineke	850	double *etaArray;
34	mmeineke	851	int test;
35	tim	837
36	gezelter	574	eta = 0.0;
37	mmeineke	778	oldEta = 0.0;
38	tim	837
39	mmeineke	850	// retrieve eta array from simInfo if it exists
40	tim	837	data = info->getProperty(ETAVALUE_ID);
41	mmeineke	850	if(data != NULL){
42
43	mmeineke	851	test = data->getDarray(etaArray);
44	mmeineke	850
45			if( test == 9 ){
46
47	tim	837	eta = etaArray[0];
48	mmeineke	850	delete[] etaArray;
49	tim	837	}
50	mmeineke	850	else
51			std::cerr << "NPTi error: etaArray is not length 9 (actual = " << test
52			<< ").\n"
53			<< " Simulation wil proceed with eta = 0;\n";
54	tim	837	}
55	gezelter	574	}
56
57	mmeineke	849	NPTi::~NPTi() {
58	mmeineke	778	//nothing for now
59	tim	763	}
60
61	mmeineke	849	void NPTi::resetIntegrator() {
62	mmeineke	778	eta = 0.0;
63	mmeineke	850	NPT::resetIntegrator();
64	mmeineke	778	}
65	gezelter	600
66	mmeineke	849	void NPTi::evolveEtaA() {
67	tim	837	eta += dt2 * ( instaVol * (instaPress - targetPressure) /
68	mmeineke	778	(p_convertNkBTtb2));
69			oldEta = eta;
70			}
71	gezelter	574
72	mmeineke	849	void NPTi::evolveEtaB() {
73	tim	837
74	mmeineke	778	prevEta = eta;
75	tim	837	eta = oldEta + dt2 * ( instaVol * (instaPress - targetPressure) /
76	mmeineke	778	(p_convertNkBTtb2));
77			}
78	gezelter	574
79	mmeineke	849	void NPTi::getVelScaleA(double sc[3], double vel[3]) {
80	mmeineke	778	int i;
81	gezelter	574
82	mmeineke	778	for(i=0; i<3; i++) sc[i] = vel[i] * ( chi + eta );
83			}
84	gezelter	574
85	mmeineke	849	void NPTi::getVelScaleB(double sc[3], int index ){
86	mmeineke	778	int i;
87	gezelter	600
88	mmeineke	778	for(i=0; i<3; i++) sc[i] = oldVel[index3 + i] ( chi + eta );
89			}
90	gezelter	574
91	tim	763
92	mmeineke	849	void NPTi::getPosScale(double pos[3], double COM[3],
93	mmeineke	778	int index, double sc[3]){
94			int j;
95	gezelter	574
96	mmeineke	778	for(j=0; j<3; j++)
97			sc[j] = ( oldPos[index*3+j] + pos[j]) / 2.0 - COM[j];
98	gezelter	600
99	mmeineke	778	for(j=0; j<3; j++)
100			sc[j] *= eta;
101			}
102	gezelter	600
103	mmeineke	849	void NPTi::scaleSimBox( void ){
104	gezelter	600
105	mmeineke	778	double scaleFactor;
106	gezelter	600
107	gezelter	611	scaleFactor = exp(dt*eta);
108
109	mmeineke	614	if ((scaleFactor > 1.1) \|\| (scaleFactor < 0.9)) {
110	gezelter	611	sprintf( painCave.errMsg,
111			"NPTi error: Attempting a Box scaling of more than 10 percent"
112			" check your tauBarostat, as it is probably too small!\n"
113			" eta = %lf, scaleFactor = %lf\n", eta, scaleFactor
114			);
115			painCave.isFatal = 1;
116			simError();
117	tim	837	} else {
118			info->scaleBox(scaleFactor);
119			}
120	mmeineke	614
121	gezelter	574	}
122
123	mmeineke	849	bool NPTi::etaConverged() {
124	tim	763
125	mmeineke	778	return ( fabs(prevEta - eta) <= etaTolerance );
126	gezelter	574	}
127
128	mmeineke	849	double NPTi::getConservedQuantity(void){
129	tim	763
130			double conservedQuantity;
131	tim	769	double Energy;
132			double thermostat_kinetic;
133			double thermostat_potential;
134			double barostat_kinetic;
135			double barostat_potential;
136	tim	837
137	tim	769	Energy = tStats->getTotalE();
138	tim	763
139	tim	837	thermostat_kinetic = fkBT* tt2 * chi * chi /
140	tim	769	(2.0 * eConvert);
141	tim	763
142	tim	769	thermostat_potential = fkBT* integralOfChidt / eConvert;
143	tim	763
144
145	tim	837	barostat_kinetic = 3.0 * NkBT * tb2 * eta * eta /
146	tim	769	(2.0 * eConvert);
147	tim	837
148			barostat_potential = (targetPressure * tStats->getVolume() / p_convert) /
149	tim	769	eConvert;
150	tim	767
151	tim	769	conservedQuantity = Energy + thermostat_kinetic + thermostat_potential +
152			barostat_kinetic + barostat_potential;
153	tim	837
154	mmeineke	778	// cout.width(8);
155			// cout.precision(8);
156	tim	763
157	tim	837	// cerr << info->getTime() << "\t" << Energy << "\t" << thermostat_kinetic <<
158			// "\t" << thermostat_potential << "\t" << barostat_kinetic <<
159	mmeineke	778	// "\t" << barostat_potential << "\t" << conservedQuantity << endl;
160	tim	837	return conservedQuantity;
161	tim	763	}
162	tim	837
163	mmeineke	851	char* NPTi::getAdditionalParameters(void){
164	tim	837
165	mmeineke	851	sprintf(addParamBuffer,
166			"\t%G\t%G;"
167			"\t%G\t%0.0\t%0.0;"
168			"\t%0.0\t%G\t%0.0;"
169			"\t%0.0\t%0.0\t%G;",
170			chi, integralOfChidt,
171			eta, eta, eta
172			);
173	tim	837
174	mmeineke	851	return addParamBuffer;
175	tim	837	}