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 amu*fs^-2*Ang^-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 |