From 925196a24fec54926cf6fe1c29f1797ef4afb463 Mon Sep 17 00:00:00 2001
From: Steve Plimpton <sjplimp@sandia.gov>
Date: Wed, 27 Jul 2022 09:38:23 -0600
Subject: [PATCH] virial adjustments for MDI stress

---
 examples/mdi/LATTE/README | 49 +++++++++++----------------------------
 src/MDI/fix_mdi_qm.cpp    |  9 +++----
 src/MDI/mdi_engine.cpp    |  5 ++++
 3 files changed, 23 insertions(+), 40 deletions(-)

diff --git a/examples/mdi/LATTE/README b/examples/mdi/LATTE/README
index 8955ec1a75..34c9e20330 100644
--- a/examples/mdi/LATTE/README
+++ b/examples/mdi/LATTE/README
@@ -1,4 +1,4 @@
-Examples in this directoy use LAMMPS as an MDI driver code and the
+Examples in this directory use LAMMPS as an MDI driver code and the
 LATTE tight-binding code as an MDI engine to perform ab initio MD and
 evaluate the energy/forces/virial for a series of MD snapshots or
 conformations with LATTE.
@@ -7,41 +7,8 @@ Talk about LATTE package and fix latte command
 
 ---------
 
-Build LAMMPS with MDI support via the MDI package:
-
-% cd lammps/lib/mdi
-% python Install.py -m mpi    # creates libmdi.a (not shared)
-% cd lammps/src
-% make yes-mdi
-% make mpi                    # builds LAMMPS with MDI
-
----------
-
-Build LATTE in its own Git repo:
-
-
-cd latte/git
-  cp makefiles/makefile.CHOICES.mdi makefile.CHOICES
-  toggle settings:
-    MAKELIB = OFF
-    SHARED = ON
-    MDI = ON
-  for MDI, also need to set MDI_PATH to mdi/git/build/MDI_Library
-  also commented out some LIB lines for Cuda
-  then just type "make"
-  produces exe = LATTE_DOUBLE and libfile = liblatte_mdi.so
-  the exe is not the non-MDI LATTE main, but a new one that goes into engine mode
-
-for building with fix latte:
-
-cmake -D PKG_MOLECULE=yes -D PKG_LATTE=yes -D USE_INTERNAL_LINALG=on ../cmake
-
----------
-
 copy lmp_mpi and LATTE_DOUBLE into this dir
-set LD_LIB_PATH to ~/latte/git for liblatte_mdi.so file in plugin mode
-LATTE always uses in.latte as input file, contains no element info
-LATTE always uses latte.in
+LATTE always uses latte.in as input file, contains no element info
 
 this version of LATTE does not use MPI, only OpenMP or GPU
 run LATTE with OpenMP: setenv OMP_NUM_THREADS max-core-count
@@ -59,7 +26,9 @@ mpirun -np 1 lmp_mpi -mdi "-name LMP -role DRIVER -method MPI" \
   -log log.aimd.driver.mpi.1 -in in.mdi.aimd : \
   -np 1 LATTE_DOUBLE -mdi "-name LATTE -role ENGINE -method MPI"
 
-mpirun -np 1 valgrind lmp_mpi -mdi "-name LMP -role DRIVER -method MPI" \
+++ Run with MPI: 2 procs for LAMMPS, 1 for LATTE
+
+mpirun -np 2 lmp_mpi -mdi "-name LMP -role DRIVER -method MPI" \
   -log log.aimd.driver.mpi.1 -in in.mdi.aimd : \
   -np 1 LATTE_DOUBLE -mdi "-name LATTE -role ENGINE -method MPI"
 
@@ -69,6 +38,14 @@ lmp_mpi -mdi \
   "-name LMP -role DRIVER -method LINK -plugin_path /home/sjplimp/latte/git" \
   -log log.aimd.driver.plugin.1 -in in.mdi.aimd.plugin
 
+mpirun -np 2 lmp_mpi -mdi \
+  "-name LMP -role DRIVER -method LINK -plugin_path /home/sjplimp/latte/git" \
+  -log log.aimd.driver.plugin.1 -in in.mdi.aimd.plugin
+
+++ Run with fix latte: 1 proc
+
+lmp_mpi < in.aimd
+
 ---------
 
 Run sequence of configs:
diff --git a/src/MDI/fix_mdi_qm.cpp b/src/MDI/fix_mdi_qm.cpp
index 25e6870bd8..d5e60c11e0 100644
--- a/src/MDI/fix_mdi_qm.cpp
+++ b/src/MDI/fix_mdi_qm.cpp
@@ -371,7 +371,8 @@ void FixMDIQM::post_force(int vflag)
     }
   }
 
-  // optionally request stress tensor from MDI engine, convert to virial
+  // optionally request stress tensor from MDI engine, convert to 6-value virial
+  // MDI defines virial tensor as intensive (divided by volume), LAMMPS does not
   // qm_virial = fix output for global QM virial
 
   if (virialflag) {
@@ -607,18 +608,18 @@ void FixMDIQM::unit_conversions()
     mdi2lmp_force = angstrom_to_bohr / ev_to_hartree;
   }
 
-  // pressure or stress units = force/area = energy/volume
+  // stress units = force/area = energy/volume
 
   mdi2lmp_pressure = 1.0;
   lmp2mdi_pressure = 1.0;
 
   if (lmpunits == REAL) {
     lmp2mdi_pressure = (kelvin_to_hartree / force->boltz) /
-        (angstrom_to_bohr * angstrom_to_bohr * angstrom_to_bohr) / force->nktv2p;
+      (angstrom_to_bohr * angstrom_to_bohr * angstrom_to_bohr);
     mdi2lmp_pressure = 1.0 / lmp2mdi_pressure;
   } else if (lmpunits == METAL) {
     lmp2mdi_pressure =
-        ev_to_hartree / (angstrom_to_bohr * angstrom_to_bohr * angstrom_to_bohr) / force->nktv2p;
+      ev_to_hartree / (angstrom_to_bohr * angstrom_to_bohr * angstrom_to_bohr);
     mdi2lmp_pressure = 1.0 / lmp2mdi_pressure;
   }
 }
diff --git a/src/MDI/mdi_engine.cpp b/src/MDI/mdi_engine.cpp
index 91774659ff..aeb18e856c 100644
--- a/src/MDI/mdi_engine.cpp
+++ b/src/MDI/mdi_engine.cpp
@@ -1493,6 +1493,9 @@ void MDIEngine::send_pe()
 /* ----------------------------------------------------------------------
    <STRESS command
    send 9-component stress tensor (no kinetic energy term)
+   should be intensive quantity (divided by volume in pressure compute)
+   MDI stress tensor units are energy/volume,
+     so conversion factor includes nktv2p to convert pressure back to virial
 ---------------------------------------------------------------------- */
 
 void MDIEngine::send_stress()
@@ -1837,6 +1840,8 @@ void MDIEngine::unit_conversions()
   }
 
   // pressure or stress units = force/area = energy/volume
+  // MDI energy/volume = Hartree/Bohr^3,
+  //   so need to remove LAMMPS nktv2p from pressure
 
   mdi2lmp_pressure = 1.0;
   lmp2mdi_pressure = 1.0;