From edf014d74ace834ce90f4a1aed26f0181f98b668 Mon Sep 17 00:00:00 2001
From: mkirsz <s1351949@sms.ed.ac.uk>
Date: Thu, 26 Dec 2024 23:26:45 +0000
Subject: [PATCH] Added density calculation
---
include/tadah/mlip/structure.h | 3 +++
src/structure.cpp | 8 ++++++++
2 files changed, 11 insertions(+)
diff --git a/include/tadah/mlip/structure.h b/include/tadah/mlip/structure.h
index 306ca2c..c36baf2 100644
--- a/include/tadah/mlip/structure.h
+++ b/include/tadah/mlip/structure.h
@@ -137,6 +137,9 @@ struct Structure {
/** @return volume of this structure. */
double get_volume() const;
+ /** @return density of this structure in g/cm^3 */
+ double get_density() const;
+
/** @return virial pressure calculated from the stress tensor.
*
* Units: energy/distance^3
diff --git a/src/structure.cpp b/src/structure.cpp
index 2d870bb..9db209c 100644
--- a/src/structure.cpp
+++ b/src/structure.cpp
@@ -173,6 +173,14 @@ size_t Structure::get_nn_iindex(const size_t i, const size_t j, const size_t jj)
double Structure::get_volume() const {
return cell.row(0)*(cell.row(1).cross(cell.row(2)));
}
+double Structure::get_density() const {
+ double V = cell.row(0)*(cell.row(1).cross(cell.row(2)));
+ V*=1e-24; // convert to cm^3
+ double amu = 1.66053906660e-24; // g
+ double mass = 0;
+ for (const auto& a:atoms) mass += PeriodicTable::get_mass(a.Z);
+ return amu*mass/V;
+}
double Structure::get_virial_pressure() const {
return stress.trace()/get_volume()/3;
--
GitLab