src/brains/Thermo.cpp

#include <math.h>
#include <iostream>

using namespace std;

#ifdef IS_MPI

#include <mpi.h>
#endif //is_mpi

#include "brains/Thermo.hpp"
#include "primitives/SRI.hpp"
#include "integrators/Integrator.hpp"
#include "utils/simError.h"
#include "math/MatVec3.h"

#ifdef IS_MPI

#define __C
#include "brains/mpiSimulation.hpp"
#endif // is_mpi

namespace oopse {
Thermo::Thermo(SimInfo *the_info) { info = the_info; }

Thermo::~Thermo() { }

double Thermo::getKinetic() {
    const double         e_convert =
                             4.184E - 4; // convert kcal/mol -> (amu A^2)/fs^2
          double         kinetic;
          double         amass;
          double         aVel[3],
                         aJ[3],
                         I[3][3];
          int            i,
                         j,
                         k,
                         kl;

          double         kinetic_global;
    vector<StuntDouble *>integrableObjects = info->integrableObjects;

    kinetic = 0.0;
    kinetic_global = 0.0;

    for( kl = 0; kl < integrableObjects.size(); kl++ ) {
        aVel = integrableObjects[kl]->getVel();
        amass = integrableObjects[kl]->getMass();

        for( j = 0; j < 3; j++ )
        kinetic += amass * aVel[j] * aVel[j];

        if (integrableObjects[kl]->isDirectional()) {
            aJ = integrableObjects[kl]->getJ();
            integrableObjects[kl]->getI(I);

            if (integrableObjects[kl]->isLinear()) {
                i = integrableObjects[kl]->linearAxis();
                j = (i + 1) % 3;
                k = (i + 2) % 3;
                kinetic += aJ[j] * aJ[j] / I[j][j] + aJ[k] * aJ[k] / I[k][k];
            } else {
                for( j = 0; j < 3; j++ )
                kinetic += aJ[j] * aJ[j] / I[j][j];
            }
        }
    }

#ifdef IS_MPI

    MPI_Allreduce(&kinetic, &kinetic_global, 1, MPI_DOUBLE, MPI_SUM,
                  MPI_COMM_WORLD);
    kinetic = kinetic_global;

#endif //is_mpi

    kinetic = kinetic * 0.5 / e_convert;

    return kinetic;
}

double Thermo::getPotential() {
    double    potential_local;
    double    potential;
    int       el,
              nSRI;
    Molecule *molecules;

    molecules = info->molecules;
    nSRI = info->n_SRI;

    potential_local = 0.0;
    potential = 0.0;
    potential_local += info->lrPot;

    for( el = 0; el < info->n_mol; el++ ) {
        potential_local += molecules[el].getPotential();
    }

    // Get total potential for entire system from MPI.

#ifdef IS_MPI

    MPI_Allreduce(&potential_local, &potential, 1, MPI_DOUBLE, MPI_SUM,
                  MPI_COMM_WORLD);

#else

    potential = potential_local;

#endif // is_mpi

    return potential;
}

double Thermo::getTotalE() {
    double total;

    total = this->getKinetic() + this->getPotential();
    return total;
}

double Thermo::getTemperature() {
    const double kb = 1.9872156E - 3; // boltzman's constant in kcal/(mol K)
          double temperature;

    temperature

    = ( 2.0 * this->getKinetic() ) / ((double)info->ndf * kb );

    return temperature;
}

double Thermo::getVolume() { return info->boxVol; }

double Thermo::getPressure() {

    // Relies on the calculation of the full molecular pressure tensor

    const double p_convert = 1.63882576e8;
          double press[3][3];
          double pressure;

    this->getPressureTensor(press);

    pressure = p_convert * (press[0][0] + press[1][1] + press[2][2]) / 3.0;

    return pressure;
}

double Thermo::getPressureX() {

    // Relies on the calculation of the full molecular pressure tensor

    const double p_convert = 1.63882576e8;
          double press[3][3];
          double pressureX;

    this->getPressureTensor(press);

    pressureX = p_convert * press[0][0];

    return pressureX;
}

double Thermo::getPressureY() {

    // Relies on the calculation of the full molecular pressure tensor

    const double p_convert = 1.63882576e8;
          double press[3][3];
          double pressureY;

    this->getPressureTensor(press);

    pressureY = p_convert * press[1][1];

    return pressureY;
}

double Thermo::getPressureZ() {

    // Relies on the calculation of the full molecular pressure tensor

    const double p_convert = 1.63882576e8;
          double press[3][3];
          double pressureZ;

    this->getPressureTensor(press);

    pressureZ = p_convert * press[2][2];

    return pressureZ;
}

void Thermo::getPressureTensor(double press[3][3]) {
    // returns pressure tensor in units amu*fs^-2*Ang^-1
    // routine derived via viral theorem description in:
    // Paci, E. and Marchi, M. J.Phys.Chem. 1996, 100, 4314-4322

    const double e_convert = 4.184e - 4;

          double molmass,
                 volume;
          double vcom[3];
          double p_local[9],
                 p_global[9];
          int    i,
                 j,
                 k;

    for( i = 0; i < 9; i++ ) {
        p_local[i] = 0.0;
        p_global[i] = 0.0;
    }

    // use velocities of integrableObjects and their masses:  

    for( i = 0; i < info->integrableObjects.size(); i++ ) {
        molmass = info->integrableObjects[i]->getMass();

        vcom = info->integrableObjects[i]->getVel();

        p_local[0] += molmass * (vcom[0] * vcom[0]);
        p_local[1] += molmass * (vcom[0] * vcom[1]);
        p_local[2] += molmass * (vcom[0] * vcom[2]);
        p_local[3] += molmass * (vcom[1] * vcom[0]);
        p_local[4] += molmass * (vcom[1] * vcom[1]);
        p_local[5] += molmass * (vcom[1] * vcom[2]);
        p_local[6] += molmass * (vcom[2] * vcom[0]);
        p_local[7] += molmass * (vcom[2] * vcom[1]);
        p_local[8] += molmass * (vcom[2] * vcom[2]);
    }

    // Get total for entire system from MPI.

#ifdef IS_MPI

    MPI_Allreduce(p_local, p_global, 9, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);

#else

    for( i = 0; i < 9; i++ ) {
        p_global[i] = p_local[i];
    }

#endif // is_mpi

    volume = this->getVolume();

    for( i = 0; i < 3; i++ ) {
        for( j = 0; j < 3; j++ ) {
            k = 3 * i + j;

            press

            [i][j] = (p_global[k] + info->tau[k]*e_convert) / volume;
        }
    }
}
} //end namespace oopse
Revision:	1725
Committed:	Wed Nov 10 22:01:06 2004 UTC (21 years ago) by tim
File size:	6186 byte(s)
Log Message:	another painful day (1) SimCreator, SimInfo, mpiSimulation (2) DumpReader, DumpWriter (InitializeFrom File will be removed) (3) ForceField (at least LJ) and BondType, BendType, TorsionType (4)Integrator (5)oopse.cpp (6)visitors & Dump2XYZ (7)SimpleBuilder (8)Constraint & ZConstraint
#	User	Rev	Content
1	tim	1725	#include <math.h>
2			#include <iostream>
3
4			using namespace std;
5
6			#ifdef IS_MPI
7
8			#include <mpi.h>
9			#endif //is_mpi
10
11			#include "brains/Thermo.hpp"
12			#include "primitives/SRI.hpp"
13			#include "integrators/Integrator.hpp"
14			#include "utils/simError.h"
15			#include "math/MatVec3.h"
16
17			#ifdef IS_MPI
18
19			#define __C
20			#include "brains/mpiSimulation.hpp"
21			#endif // is_mpi
22
23			namespace oopse {
24			Thermo::Thermo(SimInfo *the_info) { info = the_info; }
25
26			Thermo::~Thermo() { }
27
28			double Thermo::getKinetic() {
29			const double e_convert =
30			4.184E - 4; // convert kcal/mol -> (amu A^2)/fs^2
31			double kinetic;
32			double amass;
33			double aVel[3],
34			aJ[3],
35			I[3][3];
36			int i,
37			j,
38			k,
39			kl;
40
41			double kinetic_global;
42			vector<StuntDouble *>integrableObjects = info->integrableObjects;
43
44			kinetic = 0.0;
45			kinetic_global = 0.0;
46
47			for( kl = 0; kl < integrableObjects.size(); kl++ ) {
48			aVel = integrableObjects[kl]->getVel();
49			amass = integrableObjects[kl]->getMass();
50
51			for( j = 0; j < 3; j++ )
52			kinetic += amass * aVel[j] * aVel[j];
53
54			if (integrableObjects[kl]->isDirectional()) {
55			aJ = integrableObjects[kl]->getJ();
56			integrableObjects[kl]->getI(I);
57
58			if (integrableObjects[kl]->isLinear()) {
59			i = integrableObjects[kl]->linearAxis();
60			j = (i + 1) % 3;
61			k = (i + 2) % 3;
62			kinetic += aJ[j] * aJ[j] / I[j][j] + aJ[k] * aJ[k] / I[k][k];
63			} else {
64			for( j = 0; j < 3; j++ )
65			kinetic += aJ[j] * aJ[j] / I[j][j];
66			}
67			}
68			}
69
70			#ifdef IS_MPI
71
72			MPI_Allreduce(&kinetic, &kinetic_global, 1, MPI_DOUBLE, MPI_SUM,
73			MPI_COMM_WORLD);
74			kinetic = kinetic_global;
75
76			#endif //is_mpi
77
78			kinetic = kinetic * 0.5 / e_convert;
79
80			return kinetic;
81			}
82
83			double Thermo::getPotential() {
84			double potential_local;
85			double potential;
86			int el,
87			nSRI;
88			Molecule *molecules;
89
90			molecules = info->molecules;
91			nSRI = info->n_SRI;
92
93			potential_local = 0.0;
94			potential = 0.0;
95			potential_local += info->lrPot;
96
97			for( el = 0; el < info->n_mol; el++ ) {
98			potential_local += molecules[el].getPotential();
99			}
100
101			// Get total potential for entire system from MPI.
102
103			#ifdef IS_MPI
104
105			MPI_Allreduce(&potential_local, &potential, 1, MPI_DOUBLE, MPI_SUM,
106			MPI_COMM_WORLD);
107
108			#else
109
110			potential = potential_local;
111
112			#endif // is_mpi
113
114			return potential;
115			}
116
117			double Thermo::getTotalE() {
118			double total;
119
120			total = this->getKinetic() + this->getPotential();
121			return total;
122			}
123
124			double Thermo::getTemperature() {
125			const double kb = 1.9872156E - 3; // boltzman's constant in kcal/(mol K)
126			double temperature;
127
128			temperature
129
130			= ( 2.0 * this->getKinetic() ) / ((double)info->ndf * kb );
131
132			return temperature;
133			}
134
135			double Thermo::getVolume() { return info->boxVol; }
136
137			double Thermo::getPressure() {
138
139			// Relies on the calculation of the full molecular pressure tensor
140
141			const double p_convert = 1.63882576e8;
142			double press[3][3];
143			double pressure;
144
145			this->getPressureTensor(press);
146
147			pressure = p_convert * (press[0][0] + press[1][1] + press[2][2]) / 3.0;
148
149			return pressure;
150			}
151
152			double Thermo::getPressureX() {
153
154			// Relies on the calculation of the full molecular pressure tensor
155
156			const double p_convert = 1.63882576e8;
157			double press[3][3];
158			double pressureX;
159
160			this->getPressureTensor(press);
161
162			pressureX = p_convert * press[0][0];
163
164			return pressureX;
165			}
166
167			double Thermo::getPressureY() {
168
169			// Relies on the calculation of the full molecular pressure tensor
170
171			const double p_convert = 1.63882576e8;
172			double press[3][3];
173			double pressureY;
174
175			this->getPressureTensor(press);
176
177			pressureY = p_convert * press[1][1];
178
179			return pressureY;
180			}
181
182			double Thermo::getPressureZ() {
183
184			// Relies on the calculation of the full molecular pressure tensor
185
186			const double p_convert = 1.63882576e8;
187			double press[3][3];
188			double pressureZ;
189
190			this->getPressureTensor(press);
191
192			pressureZ = p_convert * press[2][2];
193
194			return pressureZ;
195			}
196
197			void Thermo::getPressureTensor(double press[3][3]) {
198			// returns pressure tensor in units amufs^-2Ang^-1
199			// routine derived via viral theorem description in:
200			// Paci, E. and Marchi, M. J.Phys.Chem. 1996, 100, 4314-4322
201
202			const double e_convert = 4.184e - 4;
203
204			double molmass,
205			volume;
206			double vcom[3];
207			double p_local[9],
208			p_global[9];
209			int i,
210			j,
211			k;
212
213			for( i = 0; i < 9; i++ ) {
214			p_local[i] = 0.0;
215			p_global[i] = 0.0;
216			}
217
218			// use velocities of integrableObjects and their masses:
219
220			for( i = 0; i < info->integrableObjects.size(); i++ ) {
221			molmass = info->integrableObjects[i]->getMass();
222
223			vcom = info->integrableObjects[i]->getVel();
224
225			p_local[0] += molmass * (vcom[0] * vcom[0]);
226			p_local[1] += molmass * (vcom[0] * vcom[1]);
227			p_local[2] += molmass * (vcom[0] * vcom[2]);
228			p_local[3] += molmass * (vcom[1] * vcom[0]);
229			p_local[4] += molmass * (vcom[1] * vcom[1]);
230			p_local[5] += molmass * (vcom[1] * vcom[2]);
231			p_local[6] += molmass * (vcom[2] * vcom[0]);
232			p_local[7] += molmass * (vcom[2] * vcom[1]);
233			p_local[8] += molmass * (vcom[2] * vcom[2]);
234			}
235
236			// Get total for entire system from MPI.
237
238			#ifdef IS_MPI
239
240			MPI_Allreduce(p_local, p_global, 9, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
241
242			#else
243
244			for( i = 0; i < 9; i++ ) {
245			p_global[i] = p_local[i];
246			}
247
248			#endif // is_mpi
249
250			volume = this->getVolume();
251
252			for( i = 0; i < 3; i++ ) {
253			for( j = 0; j < 3; j++ ) {
254			k = 3 * i + j;
255
256			press
257
258			[i][j] = (p_global[k] + info->tau[k]*e_convert) / volume;
259			}
260			}
261			}
262			} //end namespace oopse