OpenMD 3.1
Molecular Dynamics in the Open
Loading...
Searching...
No Matches
Thermo.cpp
1/*
2 * Copyright (c) 2004-present, The University of Notre Dame. All rights
3 * reserved.
4 *
5 * Redistribution and use in source and binary forms, with or without
6 * modification, are permitted provided that the following conditions are met:
7 *
8 * 1. Redistributions of source code must retain the above copyright notice,
9 * this list of conditions and the following disclaimer.
10 *
11 * 2. Redistributions in binary form must reproduce the above copyright notice,
12 * this list of conditions and the following disclaimer in the documentation
13 * and/or other materials provided with the distribution.
14 *
15 * 3. Neither the name of the copyright holder nor the names of its
16 * contributors may be used to endorse or promote products derived from
17 * this software without specific prior written permission.
18 *
19 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
20 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
21 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
22 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
23 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
24 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
25 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
26 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
27 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
28 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
29 * POSSIBILITY OF SUCH DAMAGE.
30 *
31 * SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your
32 * research, please cite the appropriate papers when you publish your
33 * work. Good starting points are:
34 *
35 * [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).
36 * [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).
37 * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).
38 * [4] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
39 * [5] Kuang & Gezelter, Mol. Phys., 110, 691-701 (2012).
40 * [6] Lamichhane, Gezelter & Newman, J. Chem. Phys. 141, 134109 (2014).
41 * [7] Lamichhane, Newman & Gezelter, J. Chem. Phys. 141, 134110 (2014).
42 * [8] Bhattarai, Newman & Gezelter, Phys. Rev. B 99, 094106 (2019).
43 */
44
45#include "brains/Thermo.hpp"
46
47#include <cmath>
48#include <iostream>
49
50#ifdef IS_MPI
51#include <mpi.h>
52#endif
53
55#include "types/FixedChargeAdapter.hpp"
56#include "types/FluctuatingChargeAdapter.hpp"
57#include "types/MultipoleAdapter.hpp"
58#include "utils/Constants.hpp"
59#include "utils/simError.h"
60#ifdef HAVE_QHULL
61#include "math/AlphaHull.hpp"
62#include "math/ConvexHull.hpp"
63#endif
64
65using namespace std;
66namespace OpenMD {
67
68 RealType Thermo::getTranslationalKinetic() {
70
71 if (!snap->hasTranslationalKineticEnergy) {
72 SimInfo::MoleculeIterator miter;
73 vector<StuntDouble*>::iterator iiter;
74 Molecule* mol;
75 StuntDouble* sd;
76 Vector3d vel;
77 RealType mass;
78 RealType kinetic(0.0);
79
80 for (mol = info_->beginMolecule(miter); mol != NULL;
81 mol = info_->nextMolecule(miter)) {
82 for (sd = mol->beginIntegrableObject(iiter); sd != NULL;
83 sd = mol->nextIntegrableObject(iiter)) {
84 mass = sd->getMass();
85 vel = sd->getVel();
86
87 kinetic +=
88 mass * (vel[0] * vel[0] + vel[1] * vel[1] + vel[2] * vel[2]);
89 }
90 }
91
92#ifdef IS_MPI
93 MPI_Allreduce(MPI_IN_PLACE, &kinetic, 1, MPI_REALTYPE, MPI_SUM,
94 MPI_COMM_WORLD);
95#endif
96
97 kinetic = kinetic * 0.5 / Constants::energyConvert;
98
99 snap->setTranslationalKineticEnergy(kinetic);
100 }
101 return snap->getTranslationalKineticEnergy();
102 }
103
104 RealType Thermo::getRotationalKinetic() {
106
107 if (!snap->hasRotationalKineticEnergy) {
108 SimInfo::MoleculeIterator miter;
109 vector<StuntDouble*>::iterator iiter;
110 Molecule* mol;
111 StuntDouble* sd;
112 Vector3d angMom;
113 Mat3x3d I;
114 int i, j, k;
115 RealType kinetic(0.0);
116
117 for (mol = info_->beginMolecule(miter); mol != NULL;
118 mol = info_->nextMolecule(miter)) {
119 for (sd = mol->beginIntegrableObject(iiter); sd != NULL;
120 sd = mol->nextIntegrableObject(iiter)) {
121 if (sd->isDirectional()) {
122 angMom = sd->getJ();
123 I = sd->getI();
124
125 if (sd->isLinear()) {
126 i = sd->linearAxis();
127 j = (i + 1) % 3;
128 k = (i + 2) % 3;
129 kinetic += angMom[j] * angMom[j] / I(j, j) +
130 angMom[k] * angMom[k] / I(k, k);
131 } else {
132 kinetic += angMom[0] * angMom[0] / I(0, 0) +
133 angMom[1] * angMom[1] / I(1, 1) +
134 angMom[2] * angMom[2] / I(2, 2);
135 }
136 }
137 }
138 }
139
140#ifdef IS_MPI
141 MPI_Allreduce(MPI_IN_PLACE, &kinetic, 1, MPI_REALTYPE, MPI_SUM,
142 MPI_COMM_WORLD);
143#endif
144
145 kinetic = kinetic * 0.5 / Constants::energyConvert;
146
147 snap->setRotationalKineticEnergy(kinetic);
148 }
149 return snap->getRotationalKineticEnergy();
150 }
151
152 RealType Thermo::getKinetic() {
154
155 if (!snap->hasKineticEnergy) {
156 RealType ke = getTranslationalKinetic() + getRotationalKinetic() +
157 getElectronicKinetic();
158
159 snap->setKineticEnergy(ke);
160 }
161 return snap->getKineticEnergy();
162 }
163
164 RealType Thermo::getPotential() {
165 // ForceManager computes the potential and stores it in the
166 // Snapshot. All we have to do is report it.
167
169 return snap->getPotentialEnergy();
170 }
171
172 potVec Thermo::getSelectionPotentials() {
173 // ForceManager computes the selection potentials and stores them
174 // in the Snapshot. All we have to do is report them.
175
177 return snap->getSelectionPotentials();
178 }
179
180 RealType Thermo::getTotalEnergy() {
182
183 if (!snap->hasTotalEnergy) {
184 snap->setTotalEnergy(this->getKinetic() + this->getPotential());
185 }
186
187 return snap->getTotalEnergy();
188 }
189
190 /*
191 * Returns only the nuclear portion of the temperature - see
192 * getElectronicTemperature for the electronic portion */
193 RealType Thermo::getTemperature() {
195
196 if (!snap->hasTemperature) {
197 RealType nuclearKE =
198 this->getTranslationalKinetic() + this->getRotationalKinetic();
199
200 RealType temperature {};
201
202 // With no degrees of freedom, T is not well defined, but we'll
203 // set to zero.
204
205 if (info_->getNdf() > 0)
206 temperature = (2.0 * nuclearKE) / (info_->getNdf() * Constants::kb);
207 else
208 temperature = 0.0;
209
210 snap->setTemperature(temperature);
211 }
212
213 return snap->getTemperature();
214 }
215
216 RealType Thermo::getElectronicKinetic() {
218
219 if (!snap->hasElectronicKineticEnergy) {
220 SimInfo::MoleculeIterator miter;
221 vector<Atom*>::iterator iiter;
222 Molecule* mol;
223 Atom* atom;
224 RealType cvel;
225 RealType cmass;
226 RealType kinetic(0.0);
227
228 for (mol = info_->beginMolecule(miter); mol != NULL;
229 mol = info_->nextMolecule(miter)) {
230 for (atom = mol->beginFluctuatingCharge(iiter); atom != NULL;
231 atom = mol->nextFluctuatingCharge(iiter)) {
232 cmass = atom->getChargeMass();
233 cvel = atom->getFlucQVel();
234 kinetic += cmass * cvel * cvel;
235 }
236 }
237
238#ifdef IS_MPI
239 MPI_Allreduce(MPI_IN_PLACE, &kinetic, 1, MPI_REALTYPE, MPI_SUM,
240 MPI_COMM_WORLD);
241#endif
242
243 kinetic *= 0.5;
244 snap->setElectronicKineticEnergy(kinetic);
245 }
246
247 return snap->getElectronicKineticEnergy();
248 }
249
250 RealType Thermo::getElectronicTemperature() {
252
253 if (!snap->hasElectronicTemperature) {
254 RealType eTemp = (2.0 * this->getElectronicKinetic()) /
255 (info_->getNFluctuatingCharges() * Constants::kb);
256
257 snap->setElectronicTemperature(eTemp);
258 }
259
260 return snap->getElectronicTemperature();
261 }
262
263 RealType Thermo::getNetCharge() {
265
266 if (!snap->hasNetCharge) {
267 SimInfo::MoleculeIterator miter;
268 vector<Atom*>::iterator aiter;
269 Molecule* mol;
270 Atom* atom;
271 RealType charge;
272 RealType netCharge(0.0);
273
274 for (mol = info_->beginMolecule(miter); mol != NULL;
275 mol = info_->nextMolecule(miter)) {
276 for (atom = mol->beginAtom(aiter); atom != NULL;
277 atom = mol->nextAtom(aiter)) {
278 charge = 0.0;
279
281 if (fca.isFixedCharge()) { charge = fca.getCharge(); }
282
285 if (fqa.isFluctuatingCharge()) { charge += atom->getFlucQPos(); }
286
287 netCharge += charge;
288 }
289 }
290
291#ifdef IS_MPI
292 MPI_Allreduce(MPI_IN_PLACE, &netCharge, 1, MPI_REALTYPE, MPI_SUM,
293 MPI_COMM_WORLD);
294#endif
295
296 snap->setNetCharge(netCharge);
297 }
298
299 return snap->getNetCharge();
300 }
301
302 RealType Thermo::getChargeMomentum() {
304
305 if (!snap->hasChargeMomentum) {
306 SimInfo::MoleculeIterator miter;
307 vector<Atom*>::iterator iiter;
308 Molecule* mol;
309 Atom* atom;
310 RealType cvel;
311 RealType cmass;
312 RealType momentum(0.0);
313
314 for (mol = info_->beginMolecule(miter); mol != NULL;
315 mol = info_->nextMolecule(miter)) {
316 for (atom = mol->beginFluctuatingCharge(iiter); atom != NULL;
317 atom = mol->nextFluctuatingCharge(iiter)) {
318 cmass = atom->getChargeMass();
319 cvel = atom->getFlucQVel();
320
321 momentum += cmass * cvel;
322 }
323 }
324
325#ifdef IS_MPI
326 MPI_Allreduce(MPI_IN_PLACE, &momentum, 1, MPI_REALTYPE, MPI_SUM,
327 MPI_COMM_WORLD);
328#endif
329
330 snap->setChargeMomentum(momentum);
331 }
332
333 return snap->getChargeMomentum();
334 }
335
336 std::vector<Vector3d> Thermo::getCurrentDensity() {
338 AtomTypeSet simTypes = info_->getSimulatedAtomTypes();
339
340 SimInfo::MoleculeIterator miter;
341 std::vector<Atom*>::iterator iiter;
342 std::vector<RigidBody*>::iterator ri;
343 Molecule* mol;
344 RigidBody* rb;
345 Atom* atom;
346 AtomType* atype;
347 AtomTypeSet::iterator at;
348 Vector3d Jc(0.0);
349 std::vector<Vector3d> typeJc(simTypes.size(), V3Zero);
350
351 for (mol = info_->beginMolecule(miter); mol != NULL;
352 mol = info_->nextMolecule(miter)) {
353 // change the velocities of atoms which belong to the rigidbodies
354 for (rb = mol->beginRigidBody(ri); rb != NULL;
355 rb = mol->nextRigidBody(ri)) {
356 rb->updateAtomVel();
357 }
358
359 for (atom = mol->beginAtom(iiter); atom != NULL;
360 atom = mol->nextAtom(iiter)) {
361 Vector3d v = atom->getVel();
362 RealType q = 0.0;
363 int typeIndex(-1);
364
365 atype = atom->getAtomType();
367 if (fca.isFixedCharge()) q = fca.getCharge();
369 if (fqa.isFluctuatingCharge()) q += atom->getFlucQPos();
370
371 typeIndex = -1;
372 at = std::find(simTypes.begin(), simTypes.end(), atype);
373 if (at != simTypes.end()) {
374 typeIndex = std::distance(simTypes.begin(), at);
375 }
376
377 if (typeIndex != -1) { typeJc[typeIndex] += q * v; }
378 Jc += q * v;
379 }
380 }
381
382#ifdef IS_MPI
383 MPI_Allreduce(MPI_IN_PLACE, &Jc[0], 3, MPI_REALTYPE, MPI_SUM,
384 MPI_COMM_WORLD);
385 for (unsigned int j = 0; j < simTypes.size(); j++) {
386 MPI_Allreduce(MPI_IN_PLACE, &typeJc[j][0], 3, MPI_REALTYPE, MPI_SUM,
387 MPI_COMM_WORLD);
388 }
389#endif
390
391 RealType vol = snap->getVolume();
392
393 Jc /= (vol * Constants::currentDensityConvert);
394 for (unsigned int j = 0; j < simTypes.size(); j++) {
395 typeJc[j] /= (vol * Constants::currentDensityConvert);
396 }
397
398 std::vector<Vector3d> result;
399 result.clear();
400
401 result.push_back(Jc);
402 for (unsigned int j = 0; j < simTypes.size(); j++)
403 result.push_back(typeJc[j]);
404
405 return result;
406 }
407
408 RealType Thermo::getVolume() {
410 return snap->getVolume();
411 }
412
413 RealType Thermo::getPressure() {
415
416 if (!snap->hasPressure) {
417 // Relies on the calculation of the full molecular pressure tensor
418
419 Mat3x3d tensor;
420 RealType pressure;
421
422 tensor = getPressureTensor();
423
424 pressure = Constants::pressureConvert *
425 (tensor(0, 0) + tensor(1, 1) + tensor(2, 2)) / 3.0;
426
427 snap->setPressure(pressure);
428 }
429
430 return snap->getPressure();
431 }
432
433 RealType Thermo::getPressure(Snapshot* snap) {
434 if (!snap->hasPressure) {
435 // Relies on the calculation of the full molecular pressure tensor
436 Mat3x3d tensor;
437 RealType pressure;
438
439 tensor = getPressureTensor(snap);
440
441 pressure = Constants::pressureConvert *
442 (tensor(0, 0) + tensor(1, 1) + tensor(2, 2)) / 3.0;
443
444 snap->setPressure(pressure);
445 }
446
447 return snap->getPressure();
448 }
449
451 // returns pressure tensor in units amu*fs^-2*Ang^-1
452 // routine derived via viral theorem description in:
453 // Paci, E. and Marchi, M. J.Phys.Chem. 1996, 100, 4314-4322
455
456 if (!snap->hasPressureTensor) {
457 Mat3x3d pressureTensor;
458 Mat3x3d p_tens(0.0);
459 RealType mass;
460 Vector3d vcom;
461
462 SimInfo::MoleculeIterator i;
463 vector<StuntDouble*>::iterator j;
464 Molecule* mol;
465 StuntDouble* sd;
466 for (mol = info_->beginMolecule(i); mol != NULL;
467 mol = info_->nextMolecule(i)) {
468 for (sd = mol->beginIntegrableObject(j); sd != NULL;
469 sd = mol->nextIntegrableObject(j)) {
470 mass = sd->getMass();
471 vcom = sd->getVel();
472 p_tens += mass * outProduct(vcom, vcom);
473 }
474 }
475
476#ifdef IS_MPI
477 MPI_Allreduce(MPI_IN_PLACE, p_tens.getArrayPointer(), 9, MPI_REALTYPE,
478 MPI_SUM, MPI_COMM_WORLD);
479#endif
480
481 RealType volume = this->getVolume();
482 Mat3x3d virialTensor = snap->getVirialTensor();
483
484 pressureTensor =
485 (p_tens + Constants::energyConvert * virialTensor) / volume;
486
487 snap->setPressureTensor(pressureTensor);
488 }
489 return snap->getPressureTensor();
490 }
491
493 // returns pressure tensor in units amu*fs^-2*Ang^-1
494 // routine derived via viral theorem description in:
495 // Paci, E. and Marchi, M. J.Phys.Chem. 1996, 100, 4314-4322
496 if (!snap->hasPressureTensor) {
497 Mat3x3d pressureTensor;
498 Mat3x3d p_tens(0.0);
499 RealType mass;
500 Vector3d vcom;
501
502 SimInfo::MoleculeIterator i;
503 vector<StuntDouble*>::iterator j;
504 Molecule* mol;
505 StuntDouble* sd;
506 for (mol = info_->beginMolecule(i); mol != NULL;
507 mol = info_->nextMolecule(i)) {
508 for (sd = mol->beginIntegrableObject(j); sd != NULL;
509 sd = mol->nextIntegrableObject(j)) {
510 mass = sd->getMass();
511 vcom = sd->getVel(snap);
512 p_tens += mass * outProduct(vcom, vcom);
513 }
514 }
515
516#ifdef IS_MPI
517 MPI_Allreduce(MPI_IN_PLACE, p_tens.getArrayPointer(), 9, MPI_REALTYPE,
518 MPI_SUM, MPI_COMM_WORLD);
519#endif
520
521 RealType volume = this->getVolume();
522 Mat3x3d virialTensor = snap->getVirialTensor();
523
524 pressureTensor =
525 (p_tens + Constants::energyConvert * virialTensor) / volume;
526
527 snap->setPressureTensor(pressureTensor);
528 }
529 return snap->getPressureTensor();
530 }
531
534
535 if (!snap->hasSystemDipole) {
536 SimInfo::MoleculeIterator miter;
537 vector<Atom*>::iterator aiter;
538 Molecule* mol;
539 Atom* atom;
540 RealType charge;
541 Vector3d ri(0.0);
542 Vector3d dipoleVector(0.0);
543 Vector3d nPos(0.0);
544 Vector3d pPos(0.0);
545 RealType nChg(0.0);
546 RealType pChg(0.0);
547 int nCount = 0;
548 int pCount = 0;
549
550 RealType chargeToC = 1.60217733e-19;
551 RealType angstromToM = 1.0e-10;
552 RealType debyeToCm = 3.33564095198e-30;
553
554 for (mol = info_->beginMolecule(miter); mol != NULL;
555 mol = info_->nextMolecule(miter)) {
556 for (atom = mol->beginAtom(aiter); atom != NULL;
557 atom = mol->nextAtom(aiter)) {
558 charge = 0.0;
559
561 if (fca.isFixedCharge()) { charge = fca.getCharge(); }
562
565 if (fqa.isFluctuatingCharge()) { charge += atom->getFlucQPos(); }
566
567 charge *= chargeToC;
568
569 ri = atom->getPos();
570 snap->wrapVector(ri);
571 ri *= angstromToM;
572
573 if (charge < 0.0) {
574 nPos += ri;
575 nChg -= charge;
576 nCount++;
577 } else if (charge > 0.0) {
578 pPos += ri;
579 pChg += charge;
580 pCount++;
581 }
582
583 if (atom->isDipole()) {
584 dipoleVector += atom->getDipole() * debyeToCm;
585 }
586 }
587 }
588
589#ifdef IS_MPI
590 MPI_Allreduce(MPI_IN_PLACE, &pChg, 1, MPI_REALTYPE, MPI_SUM,
591 MPI_COMM_WORLD);
592 MPI_Allreduce(MPI_IN_PLACE, &nChg, 1, MPI_REALTYPE, MPI_SUM,
593 MPI_COMM_WORLD);
594
595 MPI_Allreduce(MPI_IN_PLACE, &pCount, 1, MPI_INTEGER, MPI_SUM,
596 MPI_COMM_WORLD);
597 MPI_Allreduce(MPI_IN_PLACE, &nCount, 1, MPI_INTEGER, MPI_SUM,
598 MPI_COMM_WORLD);
599
600 MPI_Allreduce(MPI_IN_PLACE, pPos.getArrayPointer(), 3, MPI_REALTYPE,
601 MPI_SUM, MPI_COMM_WORLD);
602 MPI_Allreduce(MPI_IN_PLACE, nPos.getArrayPointer(), 3, MPI_REALTYPE,
603 MPI_SUM, MPI_COMM_WORLD);
604
605 MPI_Allreduce(MPI_IN_PLACE, dipoleVector.getArrayPointer(), 3,
606 MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
607#endif
608
609 // first load the accumulated dipole moment (if dipoles were present)
610 Vector3d boxDipole = dipoleVector;
611 // now include the dipole moment due to charges
612 // use the lesser of the positive and negative charge totals
613 RealType chg_value = nChg <= pChg ? nChg : pChg;
614
615 // find the average positions
616 if (pCount > 0 && nCount > 0) {
617 pPos /= pCount;
618 nPos /= nCount;
619 }
620
621 // dipole is from the negative to the positive (physics notation)
622 boxDipole += (pPos - nPos) * chg_value;
623 snap->setSystemDipole(boxDipole);
624 }
625
626 return snap->getSystemDipole();
627 }
628
631
632 if (!snap->hasSystemQuadrupole) {
633 SimInfo::MoleculeIterator miter;
634 vector<Atom*>::iterator aiter;
635 Molecule* mol;
636 Atom* atom;
637 RealType charge;
638 Vector3d ri(0.0);
639 Vector3d dipole(0.0);
640 Mat3x3d qpole(0.0);
641
642 RealType chargeToC = 1.60217733e-19;
643 RealType angstromToM = 1.0e-10;
644 RealType debyeToCm = 3.33564095198e-30;
645
646 for (mol = info_->beginMolecule(miter); mol != NULL;
647 mol = info_->nextMolecule(miter)) {
648 for (atom = mol->beginAtom(aiter); atom != NULL;
649 atom = mol->nextAtom(aiter)) {
650 ri = atom->getPos();
651 snap->wrapVector(ri);
652 ri *= angstromToM;
653
654 charge = 0.0;
655
657 if (fca.isFixedCharge()) { charge = fca.getCharge(); }
658
661 if (fqa.isFluctuatingCharge()) { charge += atom->getFlucQPos(); }
662
663 charge *= chargeToC;
664
665 qpole += 0.5 * charge * outProduct(ri, ri);
666
668
669 if (ma.isDipole()) {
670 dipole = atom->getDipole() * debyeToCm;
671 qpole += 0.5 * outProduct(dipole, ri);
672 qpole += 0.5 * outProduct(ri, dipole);
673 }
674
675 if (ma.isQuadrupole()) {
676 qpole += atom->getQuadrupole() * debyeToCm * angstromToM;
677 }
678 }
679 }
680
681#ifdef IS_MPI
682 MPI_Allreduce(MPI_IN_PLACE, qpole.getArrayPointer(), 9, MPI_REALTYPE,
683 MPI_SUM, MPI_COMM_WORLD);
684#endif
685
686 snap->setSystemQuadrupole(qpole);
687 }
688
689 return snap->getSystemQuadrupole();
690 }
691
692 // Returns the Heat Flux Vector for the system
693 Vector3d Thermo::getHeatFlux() {
694 Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
695 SimInfo::MoleculeIterator miter;
696 vector<StuntDouble*>::iterator iiter;
697 Molecule* mol;
698 StuntDouble* sd;
699 RigidBody::AtomIterator ai;
700 Atom* atom;
701 Vector3d vel;
702 Vector3d angMom;
703 Mat3x3d I;
704 int i;
705 int j;
706 int k;
707 RealType mass;
708
709 Vector3d x_a;
710 RealType kinetic;
711 RealType potential;
712 RealType eatom;
713 // Convective portion of the heat flux
714 Vector3d heatFluxJc = V3Zero;
715
716 /* Calculate convective portion of the heat flux */
717 for (mol = info_->beginMolecule(miter); mol != NULL;
718 mol = info_->nextMolecule(miter)) {
719 for (sd = mol->beginIntegrableObject(iiter); sd != NULL;
720 sd = mol->nextIntegrableObject(iiter)) {
721 mass = sd->getMass();
722 vel = sd->getVel();
723
724 kinetic = mass * (vel[0] * vel[0] + vel[1] * vel[1] + vel[2] * vel[2]);
725
726 if (sd->isDirectional()) {
727 angMom = sd->getJ();
728 I = sd->getI();
729
730 if (sd->isLinear()) {
731 i = sd->linearAxis();
732 j = (i + 1) % 3;
733 k = (i + 2) % 3;
734 kinetic += angMom[j] * angMom[j] / I(j, j) +
735 angMom[k] * angMom[k] / I(k, k);
736 } else {
737 kinetic += angMom[0] * angMom[0] / I(0, 0) +
738 angMom[1] * angMom[1] / I(1, 1) +
739 angMom[2] * angMom[2] / I(2, 2);
740 }
741 }
742
743 potential = 0.0;
744
745 if (sd->isRigidBody()) {
746 RigidBody* rb = dynamic_cast<RigidBody*>(sd);
747 for (atom = rb->beginAtom(ai); atom != NULL;
748 atom = rb->nextAtom(ai)) {
749 potential += atom->getParticlePot();
750 }
751 } else {
752 potential = sd->getParticlePot();
753 }
754
755 potential *= Constants::energyConvert; // amu A^2/fs^2
756 // The potential may not be a 1/2 factor
757 eatom = (kinetic + potential) / 2.0; // amu A^2/fs^2
758 heatFluxJc[0] += eatom * vel[0]; // amu A^3/fs^3
759 heatFluxJc[1] += eatom * vel[1]; // amu A^3/fs^3
760 heatFluxJc[2] += eatom * vel[2]; // amu A^3/fs^3
761 }
762 }
763
764 /* The J_v vector is reduced in the forceManager so everyone has
765 * the global Jv. Jc is computed over the local atoms and must be
766 * reduced among all processors.
767 */
768#ifdef IS_MPI
769 MPI_Allreduce(MPI_IN_PLACE, &heatFluxJc[0], 3, MPI_REALTYPE, MPI_SUM,
770 MPI_COMM_WORLD);
771#endif
772
773 // (kcal/mol * A/fs) * conversion => (amu A^3)/fs^3
774
775 Vector3d heatFluxJv =
776 currSnapshot->getConductiveHeatFlux() * Constants::energyConvert;
777
778 // Correct for the fact the flux is 1/V (Jc + Jv)
779 return (heatFluxJv + heatFluxJc) / this->getVolume(); // amu / fs^3
780 }
781
784
785 if (!snap->hasCOMvel) {
786 SimInfo::MoleculeIterator i;
787 Molecule* mol;
788
789 Vector3d comVel(0.0);
790 RealType totalMass(0.0);
791
792 for (mol = info_->beginMolecule(i); mol != NULL;
793 mol = info_->nextMolecule(i)) {
794 RealType mass = mol->getMass();
795 totalMass += mass;
796 comVel += mass * mol->getComVel();
797 }
798
799#ifdef IS_MPI
800 MPI_Allreduce(MPI_IN_PLACE, &totalMass, 1, MPI_REALTYPE, MPI_SUM,
801 MPI_COMM_WORLD);
802 MPI_Allreduce(MPI_IN_PLACE, comVel.getArrayPointer(), 3, MPI_REALTYPE,
803 MPI_SUM, MPI_COMM_WORLD);
804#endif
805
806 comVel /= totalMass;
807 snap->setCOMvel(comVel);
808 }
809 return snap->getCOMvel();
810 }
811
814
815 if (!snap->hasCOM) {
816 SimInfo::MoleculeIterator i;
817 Molecule* mol;
818
819 Vector3d com(0.0);
820 RealType totalMass(0.0);
821
822 for (mol = info_->beginMolecule(i); mol != NULL;
823 mol = info_->nextMolecule(i)) {
824 RealType mass = mol->getMass();
825 totalMass += mass;
826 com += mass * mol->getCom();
827 }
828
829#ifdef IS_MPI
830 MPI_Allreduce(MPI_IN_PLACE, &totalMass, 1, MPI_REALTYPE, MPI_SUM,
831 MPI_COMM_WORLD);
832 MPI_Allreduce(MPI_IN_PLACE, com.getArrayPointer(), 3, MPI_REALTYPE,
833 MPI_SUM, MPI_COMM_WORLD);
834#endif
835
836 com /= totalMass;
837 snap->setCOM(com);
838 }
839 return snap->getCOM();
840 }
841
842 /**
843 * Returns center of mass and center of mass velocity in one
844 * function call.
845 */
846 void Thermo::getComAll(Vector3d& com, Vector3d& comVel) {
848
849 if (!(snap->hasCOM && snap->hasCOMvel)) {
850 SimInfo::MoleculeIterator i;
851 Molecule* mol;
852
853 RealType totalMass(0.0);
854
855 com = 0.0;
856 comVel = 0.0;
857
858 for (mol = info_->beginMolecule(i); mol != NULL;
859 mol = info_->nextMolecule(i)) {
860 RealType mass = mol->getMass();
861 totalMass += mass;
862 com += mass * mol->getCom();
863 comVel += mass * mol->getComVel();
864 }
865
866#ifdef IS_MPI
867 MPI_Allreduce(MPI_IN_PLACE, &totalMass, 1, MPI_REALTYPE, MPI_SUM,
868 MPI_COMM_WORLD);
869 MPI_Allreduce(MPI_IN_PLACE, com.getArrayPointer(), 3, MPI_REALTYPE,
870 MPI_SUM, MPI_COMM_WORLD);
871 MPI_Allreduce(MPI_IN_PLACE, comVel.getArrayPointer(), 3, MPI_REALTYPE,
872 MPI_SUM, MPI_COMM_WORLD);
873#endif
874
875 com /= totalMass;
876 comVel /= totalMass;
877 snap->setCOM(com);
878 snap->setCOMvel(comVel);
879 }
880 com = snap->getCOM();
881 comVel = snap->getCOMvel();
882 return;
883 }
884
885 /**
886 * \brief Return inertia tensor for entire system and angular momentum
887 * Vector.
888 *
889 *
890 *
891 * [ Ixx -Ixy -Ixz ]
892 * I =| -Iyx Iyy -Iyz |
893 * [ -Izx -Iyz Izz ]
894 */
896 Vector3d& angularMomentum) {
898 Molecule::IntegrableObjectIterator j;
899 StuntDouble* sd;
900
901 if (!(snap->hasInertiaTensor && snap->hasCOMw)) {
902 RealType xx = 0.0;
903 RealType yy = 0.0;
904 RealType zz = 0.0;
905 RealType xy = 0.0;
906 RealType xz = 0.0;
907 RealType yz = 0.0;
908 Vector3d com(0.0);
909 Vector3d comVel(0.0);
910
911 getComAll(com, comVel);
912
913 SimInfo::MoleculeIterator i;
914 Molecule* mol;
915
916 Vector3d r(0.0);
917 Vector3d v(0.0);
918 Mat3x3d Itmp{};
919
920 RealType m = 0.0;
921
922 for (mol = info_->beginMolecule(i); mol != NULL;
923 mol = info_->nextMolecule(i)) {
924 for (sd = mol->beginIntegrableObject(j); sd != NULL;
925 sd = mol->nextIntegrableObject(j)) {
926 r = sd->getPos() - com;
927 v = sd->getVel() - comVel;
928
929 m = sd->getMass();
930
931 // Compute moment of intertia coefficients.
932 xx += r[0] * r[0] * m;
933 yy += r[1] * r[1] * m;
934 zz += r[2] * r[2] * m;
935
936 // compute products of intertia
937 xy += r[0] * r[1] * m;
938 xz += r[0] * r[2] * m;
939 yz += r[1] * r[2] * m;
940
941 angularMomentum += cross(r, v) * m;
942
943 if (sd->isDirectional()) {
944 RotMat3x3d A = sd->getA();
945 Itmp += A.transpose() * sd->getI() * A;
946 angularMomentum += sd->getJ();
947 }
948 }
949 }
950
951 inertiaTensor(0, 0) = yy + zz;
952 inertiaTensor(0, 1) = -xy;
953 inertiaTensor(0, 2) = -xz;
954 inertiaTensor(1, 0) = -xy;
955 inertiaTensor(1, 1) = xx + zz;
956 inertiaTensor(1, 2) = -yz;
957 inertiaTensor(2, 0) = -xz;
958 inertiaTensor(2, 1) = -yz;
959 inertiaTensor(2, 2) = xx + yy;
960
961 inertiaTensor += Itmp;
962
963#ifdef IS_MPI
964 MPI_Allreduce(MPI_IN_PLACE, inertiaTensor.getArrayPointer(), 9,
965 MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
966 MPI_Allreduce(MPI_IN_PLACE, angularMomentum.getArrayPointer(), 3,
967 MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
968#endif
969
970 snap->setCOMw(angularMomentum);
971 snap->setInertiaTensor(inertiaTensor);
972 }
973
974 angularMomentum = snap->getCOMw();
975 inertiaTensor = snap->getInertiaTensor();
976
977 return;
978 }
979
982
983 if (!(snap->hasBoundingBox)) {
984 SimInfo::MoleculeIterator i;
985 Molecule::RigidBodyIterator ri;
986 Molecule::AtomIterator ai;
987 Molecule* mol;
988 RigidBody* rb;
989 Atom* atom;
990 Vector3d pos, bMax, bMin;
991 int index = 0;
992
993 for (mol = info_->beginMolecule(i); mol != NULL;
994 mol = info_->nextMolecule(i)) {
995 // change the positions of atoms which belong to the rigidbodies
996 for (rb = mol->beginRigidBody(ri); rb != NULL;
997 rb = mol->nextRigidBody(ri)) {
998 rb->updateAtoms();
999 }
1000
1001 for (atom = mol->beginAtom(ai); atom != NULL;
1002 atom = mol->nextAtom(ai)) {
1003 pos = atom->getPos();
1004
1005 if (index == 0) {
1006 bMax = pos;
1007 bMin = pos;
1008 } else {
1009 for (int i = 0; i < 3; i++) {
1010 bMax[i] = max(bMax[i], pos[i]);
1011 bMin[i] = min(bMin[i], pos[i]);
1012 }
1013 }
1014 index++;
1015 }
1016 }
1017
1018#ifdef IS_MPI
1019 MPI_Allreduce(MPI_IN_PLACE, &bMax[0], 3, MPI_REALTYPE, MPI_MAX,
1020 MPI_COMM_WORLD);
1021
1022 MPI_Allreduce(MPI_IN_PLACE, &bMin[0], 3, MPI_REALTYPE, MPI_MIN,
1023 MPI_COMM_WORLD);
1024#endif
1025 Mat3x3d bBox = Mat3x3d(0.0);
1026 for (int i = 0; i < 3; i++) {
1027 bBox(i, i) = bMax[i] - bMin[i];
1028 }
1029 snap->setBoundingBox(bBox);
1030 }
1031
1032 return snap->getBoundingBox();
1033 }
1034
1035 // Returns the angular momentum of the system
1038
1039 if (!snap->hasCOMw) {
1040 Vector3d com(0.0);
1041 Vector3d comVel(0.0);
1042 Vector3d angularMomentum(0.0);
1043
1044 getComAll(com, comVel);
1045
1046 SimInfo::MoleculeIterator i;
1047 Molecule* mol;
1048
1049 Vector3d thisr(0.0);
1050 Vector3d thisp(0.0);
1051
1052 RealType thisMass;
1053
1054 for (mol = info_->beginMolecule(i); mol != NULL;
1055 mol = info_->nextMolecule(i)) {
1056 thisMass = mol->getMass();
1057 thisr = mol->getCom() - com;
1058 thisp = (mol->getComVel() - comVel) * thisMass;
1059
1060 angularMomentum += cross(thisr, thisp);
1061 }
1062
1063#ifdef IS_MPI
1064 MPI_Allreduce(MPI_IN_PLACE, angularMomentum.getArrayPointer(), 3,
1065 MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1066#endif
1067
1068 snap->setCOMw(angularMomentum);
1069 }
1070
1071 return snap->getCOMw();
1072 }
1073
1074 /**
1075 * Returns the Volume of the system based on a ellipsoid with
1076 * semi-axes based on the radius of gyration V=4/3*Pi*R_1*R_2*R_3
1077 * where R_i are related to the principle inertia moments
1078 * R_i = sqrt(C*I_i/N), this reduces to
1079 * V = 4/3*Pi*(C/N)^3/2*sqrt(det(I)).
1080 * See S.E. Baltazar et. al. Comp. Mat. Sci. 37 (2006) 526-536.
1081 */
1084
1085 if (!snap->hasGyrationalVolume) {
1086 Mat3x3d intTensor;
1087 RealType det;
1088 Vector3d dummyAngMom;
1089 RealType sysconstants;
1090 RealType geomCnst;
1091 RealType volume;
1092
1093 geomCnst = 3.0 / 2.0;
1094 /* Get the inertial tensor and angular momentum for free*/
1095 getInertiaTensor(intTensor, dummyAngMom);
1096
1097 det = intTensor.determinant();
1098 sysconstants =
1099 geomCnst / (RealType)(info_->getNGlobalIntegrableObjects());
1100 volume =
1101 4.0 / 3.0 * Constants::PI * pow(sysconstants, geomCnst) * sqrt(det);
1102
1103 snap->setGyrationalVolume(volume);
1104 }
1105 return snap->getGyrationalVolume();
1106 }
1107
1108 void Thermo::getGyrationalVolume(RealType& volume, RealType& detI) {
1110
1111 if (!(snap->hasInertiaTensor && snap->hasGyrationalVolume)) {
1112 Mat3x3d intTensor;
1113 Vector3d dummyAngMom;
1114 RealType sysconstants;
1115 RealType geomCnst;
1116
1117 geomCnst = 3.0 / 2.0;
1118 /* Get the inertia tensor and angular momentum for free*/
1119 this->getInertiaTensor(intTensor, dummyAngMom);
1120
1121 detI = intTensor.determinant();
1122 sysconstants =
1123 geomCnst / (RealType)(info_->getNGlobalIntegrableObjects());
1124 volume =
1125 4.0 / 3.0 * Constants::PI * pow(sysconstants, geomCnst) * sqrt(detI);
1126 snap->setGyrationalVolume(volume);
1127 } else {
1128 volume = snap->getGyrationalVolume();
1129 detI = snap->getInertiaTensor().determinant();
1130 }
1131 return;
1132 }
1133
1134 RealType Thermo::getTaggedAtomPairDistance() {
1135 Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
1136 Globals* simParams = info_->getSimParams();
1137
1138 if (simParams->haveTaggedAtomPair() &&
1139 simParams->havePrintTaggedPairDistance()) {
1140 if (simParams->getPrintTaggedPairDistance()) {
1141 pair<int, int> tap = simParams->getTaggedAtomPair();
1142 Vector3d pos1, pos2, rab;
1143
1144#ifdef IS_MPI
1145 int mol1 = info_->getGlobalMolMembership(tap.first);
1146 int mol2 = info_->getGlobalMolMembership(tap.second);
1147
1148 int proc1 = info_->getMolToProc(mol1);
1149 int proc2 = info_->getMolToProc(mol2);
1150
1151 RealType data[3];
1152 if (proc1 == worldRank) {
1153 StuntDouble* sd1 = info_->getIOIndexToIntegrableObject(tap.first);
1154 pos1 = sd1->getPos();
1155 data[0] = pos1.x();
1156 data[1] = pos1.y();
1157 data[2] = pos1.z();
1158 MPI_Bcast(data, 3, MPI_REALTYPE, proc1, MPI_COMM_WORLD);
1159 } else {
1160 MPI_Bcast(data, 3, MPI_REALTYPE, proc1, MPI_COMM_WORLD);
1161 pos1 = Vector3d(data);
1162 }
1163
1164 if (proc2 == worldRank) {
1165 StuntDouble* sd2 = info_->getIOIndexToIntegrableObject(tap.second);
1166 pos2 = sd2->getPos();
1167 data[0] = pos2.x();
1168 data[1] = pos2.y();
1169 data[2] = pos2.z();
1170 MPI_Bcast(data, 3, MPI_REALTYPE, proc2, MPI_COMM_WORLD);
1171 } else {
1172 MPI_Bcast(data, 3, MPI_REALTYPE, proc2, MPI_COMM_WORLD);
1173 pos2 = Vector3d(data);
1174 }
1175#else
1176 StuntDouble* at1 = info_->getIOIndexToIntegrableObject(tap.first);
1177 StuntDouble* at2 = info_->getIOIndexToIntegrableObject(tap.second);
1178 pos1 = at1->getPos();
1179 pos2 = at2->getPos();
1180#endif
1181 rab = pos2 - pos1;
1182 currSnapshot->wrapVector(rab);
1183 return rab.length();
1184 }
1185 return 0.0;
1186 }
1187 return 0.0;
1188 }
1189
1190 RealType Thermo::getHullVolume() {
1191#ifdef HAVE_QHULL
1193 if (!snap->hasHullVolume) {
1194 Hull* surfaceMesh_;
1195
1196 Globals* simParams = info_->getSimParams();
1197 const std::string ht = simParams->getHULL_Method();
1198
1199 if (ht == "Convex") {
1200 surfaceMesh_ = new ConvexHull();
1201 } else if (ht == "AlphaShape") {
1202 surfaceMesh_ = new AlphaHull(simParams->getAlpha());
1203 } else {
1204 return 0.0;
1205 }
1206
1207 // Build a vector of stunt doubles to determine if they are
1208 // surface atoms
1209 std::vector<StuntDouble*> localSites_;
1210 Molecule* mol;
1211 StuntDouble* sd;
1212 SimInfo::MoleculeIterator i;
1213 Molecule::IntegrableObjectIterator j;
1214
1215 for (mol = info_->beginMolecule(i); mol != NULL;
1216 mol = info_->nextMolecule(i)) {
1217 for (sd = mol->beginIntegrableObject(j); sd != NULL;
1218 sd = mol->nextIntegrableObject(j)) {
1219 localSites_.push_back(sd);
1220 }
1221 }
1222
1223 // Compute surface Mesh
1224 surfaceMesh_->computeHull(localSites_);
1225 snap->setHullVolume(surfaceMesh_->getVolume());
1226
1227 delete surfaceMesh_;
1228 }
1229
1230 return snap->getHullVolume();
1231#else
1232 return 0.0;
1233#endif
1234 }
1235
1236} // namespace OpenMD
Compute alpha complex or alpha shape.
Definition AlphaHull.hpp:79
AtomType * getAtomType()
Returns the AtomType of this Atom.
Definition A.hpp:93
AtomType is what OpenMD looks to for unchanging data about an atom.
Definition AtomType.hpp:66
RealType getMass()
get total mass of this molecule
Definition Molecule.cpp:266
Vector3d getComVel()
Returns the velocity of center of mass of this molecule.
Definition Molecule.cpp:351
Vector3d getCom()
Returns the current center of mass position of this molecule.
Definition Molecule.cpp:298
Real * getArrayPointer()
Returns the pointer of internal array.
void updateAtoms()
update the positions of atoms belong to this rigidbody
int getNdf()
Returns the number of degrees of freedom.
Definition SimInfo.hpp:220
int getNGlobalIntegrableObjects()
Returns the total number of integrable objects (total number of rigid bodies plus the total number of...
Definition SimInfo.hpp:139
int getMolToProc(int globalIndex)
Finds the processor where a molecule resides.
Definition SimInfo.hpp:667
Molecule * beginMolecule(MoleculeIterator &i)
Returns the first molecule in this SimInfo and intialize the iterator.
Definition SimInfo.cpp:240
int getNFluctuatingCharges()
Returns the total number of fluctuating charges that are present.
Definition SimInfo.hpp:217
Molecule * nextMolecule(MoleculeIterator &i)
Returns the next avaliable Molecule based on the iterator.
Definition SimInfo.cpp:245
SnapshotManager * getSnapshotManager()
Returns the snapshot manager.
Definition SimInfo.hpp:248
StuntDouble * getIOIndexToIntegrableObject(int index)
return an integral objects by its global index.
Definition SimInfo.cpp:1027
AtomTypeSet getSimulatedAtomTypes()
Returns the set of atom types present in this simulation.
Definition SimInfo.cpp:712
The Snapshot class is a repository storing dynamic data during a Simulation.
Definition Snapshot.hpp:147
Mat3x3d getBoundingBox()
Returns the Bounding Box.
Definition Snapshot.cpp:281
void setBoundingBox(const Mat3x3d &m)
Sets the Bounding Box.
Definition Snapshot.cpp:284
void wrapVector(Vector3d &v)
Wrapping the vector according to periodic boundary condition.
Definition Snapshot.cpp:337
Snapshot * getCurrentSnapshot()
Returns the pointer of current snapshot.
Real determinant() const
Returns the determinant of this matrix.
"Don't move, or you're dead! Stand up! Captain, we've got them!"
RotMat3x3d getA()
Returns the current rotation matrix of this stuntDouble.
Vector3d getVel()
Returns the current velocity of this stuntDouble.
int linearAxis()
Returns the linear axis of the rigidbody, atom and directional atom will always return -1.
RealType getMass()
Returns the mass of this stuntDouble.
Mat3x3d getQuadrupole()
Returns the current quadrupole tensor of this stuntDouble.
virtual Mat3x3d getI()=0
Returns the inertia tensor of this stuntDouble.
bool isLinear()
Tests the if this stuntDouble is a linear rigidbody.
Vector3d getPos()
Returns the current position of this stuntDouble.
RealType getFlucQPos()
Returns the current fluctuating charge of this stuntDouble.
Vector3d getDipole()
Returns the current dipole vector of this stuntDouble.
RealType getParticlePot()
Returns the current particlePot of this stuntDouble.
bool isRigidBody()
Tests if this stuntDouble is a rigid body.
Vector3d getJ()
Returns the current angular momentum of this stuntDouble (body -fixed).
bool isDirectional()
Tests if this stuntDouble is a directional one.
RealType getFlucQVel()
Returns the current charge velocity of this stuntDouble.
Mat3x3d getPressureTensor()
gives the pressure tensor in amu*fs^-2*Ang^-1
Definition Thermo.cpp:450
Vector3d getSystemDipole()
accumulate and return the simulation box dipole moment in C*m
Definition Thermo.cpp:532
void getComAll(Vector3d &com, Vector3d &comVel)
Returns the center of the mass and Center of Mass velocity of the whole system.
Definition Thermo.cpp:846
void getInertiaTensor(Mat3x3d &inertiaTensor, Vector3d &angularMomentum)
Returns the inertia tensor and the total angular momentum for for the entire system.
Definition Thermo.cpp:895
RealType getGyrationalVolume()
Returns volume of system as estimated by an ellipsoid defined by the radii of gyration.
Definition Thermo.cpp:1082
Vector3d getCom()
Returns the center of the mass of the whole system.
Definition Thermo.cpp:812
Mat3x3d getSystemQuadrupole()
accumulate and return the simulation box dipole moment in debye Angstroms
Definition Thermo.cpp:629
Vector3d getAngularMomentum()
Returns system angular momentum.
Definition Thermo.cpp:1036
Vector3d getComVel()
Returns the velocity of center of mass of the whole system.
Definition Thermo.cpp:782
Mat3x3d getBoundingBox()
Returns the Axis-aligned bounding box for the current system.
Definition Thermo.cpp:980
Real & z()
Returns reference of the third element of Vector3.
Definition Vector3.hpp:120
Real & x()
Returns reference of the first element of Vector3.
Definition Vector3.hpp:96
Real & y()
Returns reference of the second element of Vector3.
Definition Vector3.hpp:108
Real length()
Returns the length of this vector.
Definition Vector.hpp:393
Real * getArrayPointer()
Returns the pointer of internal array.
Definition Vector.hpp:171
This basic Periodic Table class was originally taken from the data.cpp file in OpenBabel.
Vector3< Real > cross(const Vector3< Real > &v1, const Vector3< Real > &v2)
Returns the cross product of two Vectors.
Definition Vector3.hpp:136