From a4616d396726876a3701ce059dd69cd2fee8c686 Mon Sep 17 00:00:00 2001
From: "Andrew J. Schultz" <ajschult784@gmail.com>
Date: Mon, 8 Feb 2021 13:10:42 -0500
Subject: [PATCH 01/11] Read initial (lattice) positions into r_p_init

---
 Src/read_write_checkpoint.f90 | 3 +++
 Src/type_definitions.f90      | 2 +-
 2 files changed, 4 insertions(+), 1 deletion(-)

diff --git a/Src/read_write_checkpoint.f90 b/Src/read_write_checkpoint.f90
index 16230d01..48864b60 100755
--- a/Src/read_write_checkpoint.f90
+++ b/Src/read_write_checkpoint.f90
@@ -441,6 +441,9 @@ SUBROUTINE Read_Config(ibox)
                   atom_list(ia,this_im,is)%rxp, &
                   atom_list(ia,this_im,is)%ryp, &
                   atom_list(ia,this_im,is)%rzp
+             atom_list(ia,this_im,is)%rxp_init = atom_list(ia,this_im,is)%rxp
+             atom_list(ia,this_im,is)%ryp_init = atom_list(ia,this_im,is)%ryp
+             atom_list(ia,this_im,is)%rzp_init = atom_list(ia,this_im,is)%rzp
              ! set the frac and exist flags for this atom
              molecule_list(this_im,is)%frac = this_lambda
              atom_list(ia,this_im,is)%exist = .TRUE.
diff --git a/Src/type_definitions.f90 b/Src/type_definitions.f90
index fb140914..50ca2267 100755
--- a/Src/type_definitions.f90
+++ b/Src/type_definitions.f90
@@ -204,7 +204,7 @@ MODULE Type_Definitions
      ! atom_list has dimensions (natoms, max_molecules, nspecies)
 
      REAL(DP) :: rxp, ryp, rzp
-     REAL(DP) :: rxp_nls, ryp_nls, rzp_nls  ! The starting positions for the neighbor list
+     REAL(DP) :: rxp_init, ryp_init, rzp_init  ! The starting positions
      REAL(DP) :: rxp_old, ryp_old, rzp_old
      LOGICAL :: exist
 

From 4b00ccb185cc567a13ecf648cd32f64d16b4d6c8 Mon Sep 17 00:00:00 2001
From: "Andrew J. Schultz" <ajschult784@gmail.com>
Date: Fri, 12 Feb 2021 14:01:01 -0500
Subject: [PATCH 02/11] Compute HMA energy and write out if asked to

---
 Src/accumulate.f90            |  9 ++++++++
 Src/energy_routines.f90       | 36 +++++++++++++++++++++++++++++-
 Src/global_variables.f90      |  7 ++++++
 Src/initialize.f90            |  2 ++
 Src/input_routines.f90        | 21 ++++++++++++++++++
 Src/main.f90                  |  3 +++
 Src/simulation_properties.f90 | 42 +++++++++++++++++++++++++++++++++++
 Src/type_definitions.f90      | 25 +++++++++++++++++++++
 Src/write_properties.f90      | 25 +++++++++++++++++++++
 9 files changed, 169 insertions(+), 1 deletion(-)

diff --git a/Src/accumulate.f90 b/Src/accumulate.f90
index 6b33926c..23eddb88 100755
--- a/Src/accumulate.f90
+++ b/Src/accumulate.f90
@@ -115,6 +115,14 @@ SUBROUTINE Accumulate(ibox)
   mass_density = mass_density / box_list(ibox)%volume
   ac_mass_density(ibox,iblock) = ac_mass_density(ibox,iblock) + mass_density
 
+  IF (need_energy_HMA) THEN
+     IF (energy_HMA(ibox)%last_calc /= i_mcstep) THEN
+        energy_HMA(ibox)%last_calc = i_mcstep
+        CALL Compute_HMA(ibox)
+     END IF
+     ac_energy_HMA(ibox,iblock) = ac_energy_HMA(ibox,iblock) + energy_HMA(ibox)%total
+  END IF
+
   IF (MOD(i_mcstep,block_avg_freq) == 0) THEN
      ac_energy(ibox,iblock)%total      = ac_energy(ibox,iblock)%total      / data_points_per_block
      ac_energy(ibox,iblock)%intra      = ac_energy(ibox,iblock)%intra      / data_points_per_block
@@ -138,5 +146,6 @@ SUBROUTINE Accumulate(ibox)
         ac_density(is,ibox,iblock)     = ac_density(is,ibox,iblock)        / data_points_per_block
         ac_nmols(is,ibox,iblock)       = ac_nmols(is,ibox,iblock)          / data_points_per_block
      END DO
+     ac_energy_HMA(ibox,iblock)        = ac_energy_HMA(ibox,iblock)        / data_points_per_block
   END IF 
 END SUBROUTINE Accumulate
diff --git a/Src/energy_routines.f90 b/Src/energy_routines.f90
index 2f94e7cf..11b2bcee 100755
--- a/Src/energy_routines.f90
+++ b/Src/energy_routines.f90
@@ -2481,6 +2481,7 @@ SUBROUTINE Compute_System_Total_Force(this_box)
    W_tensor_charge(:,:,this_box) = 0.0_DP
    W_tensor_recip(:,:,this_box) = 0.0_DP
    W_tensor_elec(:,:,this_box) =  0.0_DP
+   energy_HMA(this_box)%sum_Fdr = 0.0_DP
 
    DO is = 1, nspecies
       imLOOP1: DO im_1 = 1, nmols(is,this_box)
@@ -2685,7 +2686,7 @@ SUBROUTINE Compute_AtomPair_Force &
     ! Called by: Compute_System_Total_Force
   !-----------------------------------------------------------------------------
     ! Passed to
-    REAL(DP) :: rxij,ryij,rzij,rijsq
+    REAL(DP) :: rijsq
     INTEGER :: is,im,ia,js,jm,ja
     LOGICAL :: get_vdw,get_qq
 
@@ -2706,6 +2707,11 @@ SUBROUTINE Compute_AtomPair_Force &
     ! Coulomb potential
     REAL(DP) :: qi, qj, erfc_val, prefactor
     REAL(DP) :: ewald_constant, exp_const
+    ! HMA
+    REAL(DP) :: dr_dot_product
+    REAL(DP) :: drixp, driyp, drizp, drix, driy, driz
+    REAL(DP) :: drjxp, drjyp, drjzp, drjx, drjy, drjz
+    REAL(DP) :: rxijp, ryijp, rzijp, rxij, ryij, rzij
 
     Wij_vdw = 0.0_DP
     Wij_qq = 0.0_DP
@@ -2817,6 +2823,34 @@ SUBROUTINE Compute_AtomPair_Force &
        ENDIF qq_calc
 
     ENDIF ExistCheck
+
+    HMA_calc: IF (need_energy_HMA) THEN
+       drixp = atom_list(ia,im,is)%rxp - atom_list(ia,im,is)%rxp_init
+       driyp = atom_list(ia,im,is)%ryp - atom_list(ia,im,is)%ryp_init
+       drizp = atom_list(ia,im,is)%rzp - atom_list(ia,im,is)%rzp_init
+       CALL Minimum_Image_Separation(ibox,drixp,driyp,drizp, &
+                  drix,driy,driz)
+
+       drjxp = atom_list(ja,jm,js)%rxp - atom_list(ja,jm,js)%rxp_init
+       drjyp = atom_list(ja,jm,js)%ryp - atom_list(ja,jm,js)%ryp_init
+       drjzp = atom_list(ja,jm,js)%rzp - atom_list(ja,jm,js)%rzp_init
+       CALL Minimum_Image_Separation(ibox,drjxp,drjyp,drjzp, &
+                  drjx,drjy,drjz)
+
+       rxijp = atom_list(ia,im,is)%rxp - atom_list(ja,jm,js)%rxp
+       ryijp = atom_list(ia,im,is)%ryp - atom_list(ja,jm,js)%ryp
+       rzijp = atom_list(ia,im,is)%rzp - atom_list(ja,jm,js)%rzp
+
+       ! Now get the minimum image separation
+       CALL Minimum_Image_Separation(ibox,rxijp,ryijp,rzijp, &
+                  rxij,ryij,rzij)
+
+       dr_dot_product = rxij*(drix-drjx) &
+                      + ryij*(driy-drjy) &
+                      + rzij*(driz-drjz)
+       energy_HMA(ibox)%sum_Fdr = energy_HMA(ibox)%sum_Fdr &
+                                + (Wij_qq+Wij_vdw) * dr_dot_product / rijsq
+    ENDIF HMA_calc
 !------------------------------------------------------------------------------
   CONTAINS
 
diff --git a/Src/global_variables.f90 b/Src/global_variables.f90
index 8d8322bb..c2390abf 100755
--- a/Src/global_variables.f90
+++ b/Src/global_variables.f90
@@ -271,6 +271,12 @@ MODULE Global_Variables
   TYPE(Pressure_Class), DIMENSION(:), ALLOCATABLE, TARGET :: pressure
   LOGICAL :: need_pressure
 
+  !**********************************************************************************
+  ! HMA properties
+
+  TYPE(Energy_HMA_Class), DIMENSION(:), ALLOCATABLE, TARGET :: energy_HMA
+  LOGICAL :: need_energy_HMA
+
   ! **********************************************************************************
   ! system size integers used in memory allocation.
   ! Number of species, molecules, atoms, bonds, angles, dihedrals and impropers should
@@ -423,6 +429,7 @@ MODULE Global_Variables
 
   ! will have dimensions of nbr_boxes
   REAL(DP), DIMENSION(:,:),ALLOCATABLE,TARGET :: ac_volume, ac_enthalpy, ac_pressure, ac_mass_density
+  REAL(DP), DIMENSION(:,:),ALLOCATABLE,TARGET :: ac_energy_HMA
   ! will have dimension of (nspecies,nbr_boxes)
   REAL(DP), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: ac_density, ac_nmols
 
diff --git a/Src/initialize.f90 b/Src/initialize.f90
index cd3712fa..73ad91e1 100755
--- a/Src/initialize.f90
+++ b/Src/initialize.f90
@@ -135,6 +135,7 @@ SUBROUTINE Init_Accumulators
      ALLOCATE(ac_pressure(nbr_boxes,nbr_blocks))
      ALLOCATE(ac_enthalpy(nbr_boxes,nbr_blocks))
      ALLOCATE(ac_mass_density(nbr_boxes,nbr_blocks))
+     ALLOCATE(ac_energy_HMA(nbr_boxes,nbr_blocks))
 
      ac_energy(:,:)%total = 0.0_DP
      ac_energy(:,:)%intra = 0.0_DP
@@ -163,6 +164,7 @@ SUBROUTINE Init_Accumulators
      ac_pressure = 0.0_DP
      ac_enthalpy = 0.0_DP
      ac_mass_density = 0.0_DP
+     ac_energy_HMA = 0.0_DP
   END IF
 
 END SUBROUTINE Init_Accumulators
diff --git a/Src/input_routines.f90 b/Src/input_routines.f90
index a3494c9d..c3886f3e 100755
--- a/Src/input_routines.f90
+++ b/Src/input_routines.f90
@@ -5913,6 +5913,7 @@ SUBROUTINE Get_Property_Info
   max_properties = 0
   cpcollect = .FALSE.
   need_pressure = .FALSE.
+  need_energy_HMA = .FALSE.
 
   DO
      line_nbr = line_nbr + 1
@@ -5983,6 +5984,10 @@ SUBROUTINE Get_Property_Info
            ELSE IF (line_array(1) == 'enthalpy' .OR. line_array(1) == 'Enthalpy') THEN
               nbr_properties = nbr_properties + 1
               IF (int_sim_type /= sim_npt .AND. int_sim_type /= sim_gemc_npt) need_pressure = .TRUE.
+           ELSE IF (line_array(1) == 'energy_HMA' .OR. line_array(1) == 'Energy_HMA') THEN
+              nbr_properties = nbr_properties + 1
+              need_energy_HMA = .TRUE.
+
            ELSE
               ! this is a property for the system
               nbr_properties = nbr_properties + 1
@@ -6161,6 +6166,9 @@ SUBROUTINE Get_Property_Info
               ELSE IF (line_array(1) == 'mass_density' .OR. line_array(1) == 'Mass_Density') THEN
                  nbr_properties = nbr_properties + 1
                  prop_output(nbr_properties,nbr_prop_files(this_box),this_box) = 'Mass_Density'
+              ELSE IF (line_array(1) == 'energy_HMA' .OR. line_array(1) == 'Energy_HMA') THEN
+                 nbr_properties = nbr_properties + 1
+                 prop_output(nbr_properties,nbr_prop_files(this_box),this_box) = 'Energy_HMA'
               ELSE
                 err_msg = ''
                 err_msg(1) = 'Keyword "' // TRIM(line_array(1)) // '" on line ' // &
@@ -6259,6 +6267,19 @@ SUBROUTINE Get_Property_Info
      pressure(:)%last_calc = -1
   END IF
 
+  IF (need_energy_HMA) THEN
+     ALLOCATE(energy_HMA(nbr_boxes))
+     IF (.NOT. ALLOCATED(W_tensor_charge)) ALLOCATE(W_tensor_charge(3,3,nbr_boxes) , W_tensor_recip(3,3,nbr_boxes))
+     IF (.NOT. ALLOCATED(W_tensor_vdw)) ALLOCATE(W_tensor_vdw(3,3,nbr_boxes) , W_tensor_total(3,3,nbr_boxes))
+     IF (.NOT. ALLOCATED(W_tensor_elec)) ALLOCATE(W_tensor_elec(3,3,nbr_boxes), pressure_tensor(3,3,nbr_boxes))
+     IF (.NOT. ALLOCATED(pressure)) ALLOCATE(pressure(nbr_boxes))
+
+     energy_HMA(nbr_boxes)%lattice = 0.0_DP
+     energy_HMA(nbr_boxes)%harmonic = 0.0_DP
+     energy_HMA(nbr_boxes)%anharmonic = 0.0_DP
+     energy_HMA(nbr_boxes)%total = 0.0_DP
+  END IF
+
   WRITE(logunit,'(A80)') '********************************************************************************'
 
 END SUBROUTINE Get_Property_Info
diff --git a/Src/main.f90 b/Src/main.f90
index 1e5f38f4..05afcd5a 100755
--- a/Src/main.f90
+++ b/Src/main.f90
@@ -378,6 +378,9 @@ PROGRAM Main
   WRITE(logunit,'(A80)') '********************************************************************************'
   DO ibox = 1,nbr_boxes
      CALL Check_System_Energy(ibox,.FALSE.)
+     IF (need_energy_HMA) THEN
+        energy_HMA(ibox)%lattice = energy(ibox)%total
+      END IF
   END DO
   WRITE(logunit,'(A80)') '********************************************************************************'
 
diff --git a/Src/simulation_properties.f90 b/Src/simulation_properties.f90
index ca2bbe70..402c6499 100755
--- a/Src/simulation_properties.f90
+++ b/Src/simulation_properties.f90
@@ -204,6 +204,48 @@ SUBROUTINE Compute_Pressure(this_box)
 
   END SUBROUTINE Compute_Pressure
 
+!  !****************************************************************************
+
+  SUBROUTINE Compute_HMA(this_box)
+    !***************************************************************************
+    !
+    ! This subroutine calculates the pressure of the the box
+    !
+    ! CALLED BY:
+    !       Write_Properties
+    !
+    ! CALLS :
+    !       Compute_System_Total_Force
+    !
+    ! Written by Andrew Schultz on 02/10/21
+    !
+    !***************************************************************************
+
+    IMPLICIT NONE
+
+    INTEGER :: this_box
+
+    REAL(DP) :: dxp, dyp, dzp
+    REAL(DP) :: dx, dy, dz
+    REAL(DP) :: Fdr
+
+    INTEGER :: is, mol_is, im, ia
+
+    energy_HMA(this_box)%harmonic = 1.5_DP / beta(this_box) &
+                                    * (SUM(nmols(:,this_box))-1)
+
+    ! compute Fdr
+    CALL Compute_System_Total_Force(this_box)
+
+    energy_HMA(this_box)%anharmonic = 0.5_DP * energy_HMA(this_box)%sum_Fdr &
+                                    + energy(this_box)%total &
+                                    - energy_HMA(this_box)%lattice
+    energy_HMA(this_box)%total = 0.5_DP * energy_HMA(this_box)%sum_Fdr &
+                               + energy(this_box)%total &
+                               + energy_HMA(this_box)%harmonic
+
+  END SUBROUTINE Compute_HMA
+
 !  !****************************************************************************
 
 END MODULE Simulation_Properties
diff --git a/Src/type_definitions.f90 b/Src/type_definitions.f90
index 50ca2267..de5ce258 100755
--- a/Src/type_definitions.f90
+++ b/Src/type_definitions.f90
@@ -593,6 +593,31 @@ MODULE Type_Definitions
 
  END TYPE Pressure_Class
 
+!-------------------------------------------------------------------------------------------------
+
+ TYPE Energy_HMA_Class
+    ! This class holds the HMA energy values for each box
+
+    ! Lattice energy -- energy at the beginning of the simulation
+    REAL(DP) :: lattice
+
+    ! Harmonic energy: (3/2) (N-1) kT
+    REAL(DP) :: harmonic
+
+    ! Anharmonic energy (from HMA)
+    REAL(DP) :: anharmonic
+
+    ! Total energy (from HMA): U_lat + H_harm + U_anh
+    REAL(DP) :: total
+
+    ! Sum of Fi dot dri over all atoms
+    REAL(DP) :: sum_Fdr
+
+    ! last calculation
+    INTEGER :: last_calc
+
+ END TYPE Energy_HMA_Class
+
 !-------------------------------------------------------------------------------------------------
 
   TYPE Rotation_Class
diff --git a/Src/write_properties.f90 b/Src/write_properties.f90
index a8c51104..59b9d669 100755
--- a/Src/write_properties.f90
+++ b/Src/write_properties.f90
@@ -148,6 +148,10 @@ SUBROUTINE Write_Header(file_number)
 
           prop_unit(ii) = '(kJ/mol)'
 
+       ELSE IF (prop_to_write(1:10) == 'Energy_HMA') THEN
+
+          prop_unit(ii) = '(kJ/mol)-Ext'
+
        END IF
 
        WRITE(this_unit,'(A16,2X)',ADVANCE='NO') (TRIM(prop_unit(ii)))
@@ -418,6 +422,18 @@ SUBROUTINE Write_Properties_Buffer(file_number)
             write_buff(ii+1) = mass_density / box_list(this_box)%volume
          END IF
          write_buff(ii+1) = write_buff(ii+1) * atomic_to_kgm3
+
+      ELSE IF (prop_written == 'Energy_HMA') THEN
+         IF (block_avg) THEN
+            write_buff(ii+1) = ac_energy_HMA(this_box,iblock)
+         ELSE
+            IF (energy_HMA(this_box)%last_calc /= i_mcstep) THEN
+               energy_HMA(this_box)%last_calc = i_mcstep
+               CALL Compute_HMA(this_box)
+            END IF
+            write_buff(ii+1) = energy_HMA(this_box)%total
+         END IF
+         write_buff(ii+1) = write_buff(ii+1) * atomic_to_kJmol
       END IF
       
       ! At the end increment property counter by 1
@@ -785,6 +801,15 @@ SUBROUTINE Write_Mean_Error(ibox)
             write_mean(ii+1) = write_mean(ii+1) * atomic_to_kgm3
             write_err(ii+1)  = sqrt(write_err(ii+1) / nbr_blocks) * atomic_to_kgm3
 
+         ELSE IF (prop_written == 'Energy_HMA') THEN
+
+            write_mean(ii+1) = SUM(ac_energy_HMA(ibox,:)) / nbr_blocks
+            DO iblock = 1, nbr_blocks
+               write_err(ii+1) = write_err(ii+1) + (ac_energy_HMA(ibox,iblock) - write_mean(ii+1))**2
+            END DO
+            write_mean(ii+1) = write_mean(ii+1) * atomic_to_kJmol
+            write_err(ii+1)  = sqrt(write_err(ii+1) / nbr_blocks) * atomic_to_kJmol
+
          END IF
          
          ! At the end increment property counter by 1

From 41d38cfe6566ce056984f67555d7a3cd44abc41b Mon Sep 17 00:00:00 2001
From: "Andrew J. Schultz" <ajschult784@gmail.com>
Date: Fri, 19 Feb 2021 09:52:55 -0500
Subject: [PATCH 03/11] Add HMA pressure, HMA section for input file

---
 Src/accumulate.f90            |   4 +-
 Src/energy_routines.f90       |   2 +-
 Src/global_variables.f90      |   5 +-
 Src/input_routines.f90        | 132 ++++++++++++++++++++++++++++------
 Src/main.f90                  |   4 +-
 Src/nvtmc_control.f90         |   2 +
 Src/simulation_properties.f90 |  15 +++-
 Src/type_definitions.f90      |  22 ++++++
 Src/write_properties.f90      |  24 +++++++
 9 files changed, 182 insertions(+), 28 deletions(-)

diff --git a/Src/accumulate.f90 b/Src/accumulate.f90
index 23eddb88..8413ed32 100755
--- a/Src/accumulate.f90
+++ b/Src/accumulate.f90
@@ -115,12 +115,14 @@ SUBROUTINE Accumulate(ibox)
   mass_density = mass_density / box_list(ibox)%volume
   ac_mass_density(ibox,iblock) = ac_mass_density(ibox,iblock) + mass_density
 
-  IF (need_energy_HMA) THEN
+  IF (need_HMA) THEN
      IF (energy_HMA(ibox)%last_calc /= i_mcstep) THEN
         energy_HMA(ibox)%last_calc = i_mcstep
+        pressure_HMA(ibox)%last_calc = i_mcstep
         CALL Compute_HMA(ibox)
      END IF
      ac_energy_HMA(ibox,iblock) = ac_energy_HMA(ibox,iblock) + energy_HMA(ibox)%total
+     ac_pressure_HMA(ibox,iblock) = ac_pressure_HMA(ibox,iblock) + pressure_HMA(ibox)%total
   END IF
 
   IF (MOD(i_mcstep,block_avg_freq) == 0) THEN
diff --git a/Src/energy_routines.f90 b/Src/energy_routines.f90
index 11b2bcee..c8d96383 100755
--- a/Src/energy_routines.f90
+++ b/Src/energy_routines.f90
@@ -2824,7 +2824,7 @@ SUBROUTINE Compute_AtomPair_Force &
 
     ENDIF ExistCheck
 
-    HMA_calc: IF (need_energy_HMA) THEN
+    HMA_calc: IF (need_HMA) THEN
        drixp = atom_list(ia,im,is)%rxp - atom_list(ia,im,is)%rxp_init
        driyp = atom_list(ia,im,is)%ryp - atom_list(ia,im,is)%ryp_init
        drizp = atom_list(ia,im,is)%rzp - atom_list(ia,im,is)%rzp_init
diff --git a/Src/global_variables.f90 b/Src/global_variables.f90
index c2390abf..52e89c88 100755
--- a/Src/global_variables.f90
+++ b/Src/global_variables.f90
@@ -275,7 +275,8 @@ MODULE Global_Variables
   ! HMA properties
 
   TYPE(Energy_HMA_Class), DIMENSION(:), ALLOCATABLE, TARGET :: energy_HMA
-  LOGICAL :: need_energy_HMA
+  TYPE(Pressure_HMA_Class), DIMENSION(:), ALLOCATABLE, TARGET :: pressure_HMA
+  LOGICAL :: need_HMA
 
   ! **********************************************************************************
   ! system size integers used in memory allocation.
@@ -429,7 +430,7 @@ MODULE Global_Variables
 
   ! will have dimensions of nbr_boxes
   REAL(DP), DIMENSION(:,:),ALLOCATABLE,TARGET :: ac_volume, ac_enthalpy, ac_pressure, ac_mass_density
-  REAL(DP), DIMENSION(:,:),ALLOCATABLE,TARGET :: ac_energy_HMA
+  REAL(DP), DIMENSION(:,:),ALLOCATABLE,TARGET :: ac_energy_HMA, ac_pressure_HMA
   ! will have dimension of (nspecies,nbr_boxes)
   REAL(DP), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: ac_density, ac_nmols
 
diff --git a/Src/input_routines.f90 b/Src/input_routines.f90
index c3886f3e..d5adf107 100755
--- a/Src/input_routines.f90
+++ b/Src/input_routines.f90
@@ -5913,7 +5913,6 @@ SUBROUTINE Get_Property_Info
   max_properties = 0
   cpcollect = .FALSE.
   need_pressure = .FALSE.
-  need_energy_HMA = .FALSE.
 
   DO
      line_nbr = line_nbr + 1
@@ -5984,9 +5983,6 @@ SUBROUTINE Get_Property_Info
            ELSE IF (line_array(1) == 'enthalpy' .OR. line_array(1) == 'Enthalpy') THEN
               nbr_properties = nbr_properties + 1
               IF (int_sim_type /= sim_npt .AND. int_sim_type /= sim_gemc_npt) need_pressure = .TRUE.
-           ELSE IF (line_array(1) == 'energy_HMA' .OR. line_array(1) == 'Energy_HMA') THEN
-              nbr_properties = nbr_properties + 1
-              need_energy_HMA = .TRUE.
 
            ELSE
               ! this is a property for the system
@@ -6167,8 +6163,15 @@ SUBROUTINE Get_Property_Info
                  nbr_properties = nbr_properties + 1
                  prop_output(nbr_properties,nbr_prop_files(this_box),this_box) = 'Mass_Density'
               ELSE IF (line_array(1) == 'energy_HMA' .OR. line_array(1) == 'Energy_HMA') THEN
-                 nbr_properties = nbr_properties + 1
-                 prop_output(nbr_properties,nbr_prop_files(this_box),this_box) = 'Energy_HMA'
+                 IF (need_HMA) THEN
+                    nbr_properties = nbr_properties + 1
+                    prop_output(nbr_properties,nbr_prop_files(this_box),this_box) = 'Energy_HMA'
+                  END IF
+              ELSE IF (line_array(1) == 'pressure_HMA' .OR. line_array(1) == 'Pressure_HMA') THEN
+                 IF (need_HMA) THEN
+                    nbr_properties = nbr_properties + 1
+                    prop_output(nbr_properties,nbr_prop_files(this_box),this_box) = 'Pressure_HMA'
+                 END IF
               ELSE
                 err_msg = ''
                 err_msg(1) = 'Keyword "' // TRIM(line_array(1)) // '" on line ' // &
@@ -6178,7 +6181,7 @@ SUBROUTINE Get_Property_Info
                 err_msg(4) = '  energy_angle, energy_dihedral, energy_improper, energy_intravdw, energy_intraq'
                 err_msg(5) = '  energy_intervdw, energy_interq, energy_lrc, energy_recip, energy_self,'
                 err_msg(6) = '  enthalpy, pressure, pressure_xx, pressure_yy, pressure_zz, volume, density,'
-                err_msg(7) = '  nmols, mass_density'
+                err_msg(7) = '  nmols, mass_density, energy_HMA, pressure_HMA'
                 CALL Clean_Abort(err_msg,'Get_Property_Info')
               END IF
 
@@ -6251,7 +6254,7 @@ SUBROUTINE Get_Property_Info
 
   END IF
 
-  IF (need_pressure) THEN
+  IF (need_pressure .OR. need_HMA) THEN
      ALLOCATE(W_tensor_charge(3,3,nbr_boxes) , W_tensor_recip(3,3,nbr_boxes))
      ALLOCATE(W_tensor_vdw(3,3,nbr_boxes) , W_tensor_total(3,3,nbr_boxes))
      ALLOCATE(W_tensor_elec(3,3,nbr_boxes), pressure_tensor(3,3,nbr_boxes))
@@ -6267,19 +6270,6 @@ SUBROUTINE Get_Property_Info
      pressure(:)%last_calc = -1
   END IF
 
-  IF (need_energy_HMA) THEN
-     ALLOCATE(energy_HMA(nbr_boxes))
-     IF (.NOT. ALLOCATED(W_tensor_charge)) ALLOCATE(W_tensor_charge(3,3,nbr_boxes) , W_tensor_recip(3,3,nbr_boxes))
-     IF (.NOT. ALLOCATED(W_tensor_vdw)) ALLOCATE(W_tensor_vdw(3,3,nbr_boxes) , W_tensor_total(3,3,nbr_boxes))
-     IF (.NOT. ALLOCATED(W_tensor_elec)) ALLOCATE(W_tensor_elec(3,3,nbr_boxes), pressure_tensor(3,3,nbr_boxes))
-     IF (.NOT. ALLOCATED(pressure)) ALLOCATE(pressure(nbr_boxes))
-
-     energy_HMA(nbr_boxes)%lattice = 0.0_DP
-     energy_HMA(nbr_boxes)%harmonic = 0.0_DP
-     energy_HMA(nbr_boxes)%anharmonic = 0.0_DP
-     energy_HMA(nbr_boxes)%total = 0.0_DP
-  END IF
-
   WRITE(logunit,'(A80)') '********************************************************************************'
 
 END SUBROUTINE Get_Property_Info
@@ -6724,4 +6714,104 @@ SUBROUTINE Get_Verbosity_Info
 
 END SUBROUTINE Get_Verbosity_Info
 
+
+!******************************************************************************
+SUBROUTINE Get_HMA_Info
+!******************************************************************************
+! This routine opens the input file and determines if HMA properties will be
+! computed and (optionally) the harmonic pressure.
+!
+! # HMA
+! TRUE
+! harmonic_pressure 9
+!
+! The routine searches for the keyword "# HMA" and then reads the necessary
+! information underneath the key word.
+!******************************************************************************
+
+  INTEGER :: ierr,line_nbr,nbr_entries
+  CHARACTER(STRING_LEN) :: line_string, line_array(60)
+
+!******************************************************************************
+  WRITE(logunit,*)
+  WRITE(logunit,'(A)') 'HMA'
+  WRITE(logunit,'(A80)') '********************************************************************************'
+
+  need_HMA = .FALSE.
+
+  REWIND(inputunit)
+
+  ierr = 0
+  line_nbr = 0
+
+  DO
+     line_nbr = line_nbr + 1
+     CALL Read_String(inputunit,line_string,ierr)
+
+     IF (ierr .NE. 0) THEN
+        EXIT
+     END IF
+
+     IF (line_string(1:13) == '# HMA') THEN
+        line_nbr = line_nbr + 1
+        CALL Parse_String(inputunit,line_nbr,1,nbr_entries,line_array,ierr)
+        IF (line_array(1) == 'TRUE' .OR. line_array(1) == '.TRUE.' .OR. &
+            line_array(1) == 'true' .OR. line_array(1) == '.true.') THEN
+           need_HMA = .TRUE.
+
+           ALLOCATE(energy_HMA(nbr_boxes))
+           ALLOCATE(pressure_HMA(nbr_boxes))
+
+           energy_HMA(:)%lattice = 0.0_DP
+           energy_HMA(:)%harmonic = 0.0_DP
+           energy_HMA(:)%anharmonic = 0.0_DP
+           energy_HMA(:)%total = 0.0_DP
+
+           pressure_HMA(:)%lattice = 0.0_DP
+           pressure_HMA(:)%harmonic = 0.0_DP
+           pressure_HMA(:)%anharmonic = 0.0_DP
+           pressure_HMA(:)%total = 0.0_DP
+
+        ELSE IF (line_array(1) == 'FALSE' .OR. line_array(1) == '.FALSE.' .OR. &
+                 line_array(1) == 'false' .OR. line_array(1) == '.false.') THEN
+           need_HMA = .FALSE.
+           EXIT
+        ELSE
+           err_msg = ''
+           err_msg(1) = 'Keyword ' // TRIM(line_array(1)) // ' on line number ' // &
+                        TRIM(Int_To_String(line_nbr)) // ' of the input file is not supported'
+           err_msg(2) = 'Supported keywords are: true, false'
+           CALL Clean_Abort(err_msg,'Get_Verbosity_Info')
+        END IF
+
+        line_nbr = line_nbr + 1
+        CALL Parse_String(inputunit,line_nbr,0,nbr_entries,line_array,ierr)
+
+        IF (nbr_entries == 0 .OR. line_array(1)(1:1) == '!') THEN
+           EXIT
+        END IF
+
+        IF (line_array(1) == 'harmonic_pressure' .OR. line_array(1) == 'Harmonic_Pressure') THEN
+           ! assign the pressures
+           pressure_HMA(:)%harmonic = String_To_Double(line_array(2))
+           WRITE(logunit,'(A,X,F9.3,X,A)',ADVANCE='NO') 'HMA harmonic pressure is', pressure_HMA(1)%harmonic,  'bar'
+
+           ! convert pressure into atomic units
+           pressure_HMA(:)%harmonic = pressure_HMA(:)%harmonic / atomic_to_bar
+           WRITE(logunit,'(X,A,X,E13.6,X,A)') '=', pressure_HMA(1)%harmonic,  'amu / (A ps^2)'
+        END IF
+
+        EXIT
+     ENDIF
+
+  ENDDO
+
+  IF (need_HMA) THEN
+     WRITE(logunit,'(A)') 'HMA enabled'
+  ELSE
+     WRITE(logunit,'(A)') 'HMA not enabled'
+  END IF
+  WRITE(logunit,'(A80)') '********************************************************************************'
+
+END SUBROUTINE Get_HMA_Info
 END MODULE Input_Routines
diff --git a/Src/main.f90 b/Src/main.f90
index 05afcd5a..d40fada5 100755
--- a/Src/main.f90
+++ b/Src/main.f90
@@ -378,8 +378,10 @@ PROGRAM Main
   WRITE(logunit,'(A80)') '********************************************************************************'
   DO ibox = 1,nbr_boxes
      CALL Check_System_Energy(ibox,.FALSE.)
-     IF (need_energy_HMA) THEN
+     IF (need_HMA) THEN
+        CALL Compute_Pressure(ibox)
         energy_HMA(ibox)%lattice = energy(ibox)%total
+        pressure_HMA(ibox)%lattice = pressure(ibox)%computed - pressure(ibox)%ideal
       END IF
   END DO
   WRITE(logunit,'(A80)') '********************************************************************************'
diff --git a/Src/nvtmc_control.f90 b/Src/nvtmc_control.f90
index 769b5167..a72a8f5c 100755
--- a/Src/nvtmc_control.f90
+++ b/Src/nvtmc_control.f90
@@ -117,6 +117,8 @@ SUBROUTINE NVTMC_Control
   ! Write Widom insertion info to log file, which requires units from Get_Simulation_Length_Info
   CALL Log_Widom_Info
 
+  CALL Get_HMA_Info
+
   CALL Get_Property_Info
 
   CALL Get_Rcutoff_Low
diff --git a/Src/simulation_properties.f90 b/Src/simulation_properties.f90
index 402c6499..01f4356e 100755
--- a/Src/simulation_properties.f90
+++ b/Src/simulation_properties.f90
@@ -227,12 +227,15 @@ SUBROUTINE Compute_HMA(this_box)
 
     REAL(DP) :: dxp, dyp, dzp
     REAL(DP) :: dx, dy, dz
-    REAL(DP) :: Fdr
+    Real(DP) :: fV
 
     INTEGER :: is, mol_is, im, ia
+    INTEGER :: this_nmols
+
+    this_nmols = SUM(nmols(:,this_box))
 
     energy_HMA(this_box)%harmonic = 1.5_DP / beta(this_box) &
-                                    * (SUM(nmols(:,this_box))-1)
+                                    * (this_nmols-1)
 
     ! compute Fdr
     CALL Compute_System_Total_Force(this_box)
@@ -244,6 +247,14 @@ SUBROUTINE Compute_HMA(this_box)
                                + energy(this_box)%total &
                                + energy_HMA(this_box)%harmonic
 
+    fV = (beta(this_box)*pressure_HMA(this_box)%harmonic - this_nmols / box_list(this_box)%volume) &
+       / (3*(this_nmols)-1)
+
+    pressure_HMA(this_box)%anharmonic = pressure(this_box)%computed - pressure(this_box)%ideal &
+                                      - pressure_HMA(this_box)%lattice + fV * energy_HMA(this_box)%sum_Fdr
+    pressure_HMA(this_box)%total = pressure_HMA(this_box)%harmonic + pressure(this_box)%computed &
+                                 - pressure(this_box)%ideal + fV * energy_HMA(this_box)%sum_Fdr
+
   END SUBROUTINE Compute_HMA
 
 !  !****************************************************************************
diff --git a/Src/type_definitions.f90 b/Src/type_definitions.f90
index de5ce258..1f74401b 100755
--- a/Src/type_definitions.f90
+++ b/Src/type_definitions.f90
@@ -618,6 +618,28 @@ MODULE Type_Definitions
 
  END TYPE Energy_HMA_Class
 
+!-------------------------------------------------------------------------------------------------
+
+ TYPE Pressure_HMA_Class
+    ! This class holds the HMA energy values for each box
+
+    ! Lattice pressure -- non-ideal pressure at the beginning of the simulation
+    REAL(DP) :: lattice
+
+    ! Harmonic pressure: read from input file
+    REAL(DP) :: harmonic
+
+    ! Anharmonic pressure (from HMA)
+    REAL(DP) :: anharmonic
+
+    ! Total pressure (from HMA): P_lat + P_harm + P_anh
+    REAL(DP) :: total
+
+    ! last calculation
+    INTEGER :: last_calc
+
+ END TYPE Pressure_HMA_Class
+
 !-------------------------------------------------------------------------------------------------
 
   TYPE Rotation_Class
diff --git a/Src/write_properties.f90 b/Src/write_properties.f90
index 59b9d669..d79fad8f 100755
--- a/Src/write_properties.f90
+++ b/Src/write_properties.f90
@@ -152,6 +152,10 @@ SUBROUTINE Write_Header(file_number)
 
           prop_unit(ii) = '(kJ/mol)-Ext'
 
+       ELSE IF (prop_to_write(1:12) == 'Pressure_HMA') THEN
+
+          prop_unit(ii) = '(bar)'
+
        END IF
 
        WRITE(this_unit,'(A16,2X)',ADVANCE='NO') (TRIM(prop_unit(ii)))
@@ -434,6 +438,17 @@ SUBROUTINE Write_Properties_Buffer(file_number)
             write_buff(ii+1) = energy_HMA(this_box)%total
          END IF
          write_buff(ii+1) = write_buff(ii+1) * atomic_to_kJmol
+      ELSE IF (prop_written == 'Pressure_HMA') THEN
+         IF (block_avg) THEN
+            write_buff(ii+1) = ac_pressure_HMA(this_box,iblock)
+         ELSE
+            IF (pressure_HMA(this_box)%last_calc /= i_mcstep) THEN
+               pressure_HMA(this_box)%last_calc = i_mcstep
+               CALL Compute_HMA(this_box)
+            END IF
+            write_buff(ii+1) = pressure_HMA(this_box)%total
+         END IF
+         write_buff(ii+1) = write_buff(ii+1) * atomic_to_bar
       END IF
       
       ! At the end increment property counter by 1
@@ -810,6 +825,15 @@ SUBROUTINE Write_Mean_Error(ibox)
             write_mean(ii+1) = write_mean(ii+1) * atomic_to_kJmol
             write_err(ii+1)  = sqrt(write_err(ii+1) / nbr_blocks) * atomic_to_kJmol
 
+         ELSE IF (prop_written == 'Pressure_HMA') THEN
+
+            write_mean(ii+1) = SUM(ac_pressure_HMA(ibox,:)) / nbr_blocks
+            DO iblock = 1, nbr_blocks
+               write_err(ii+1) = write_err(ii+1) + (ac_pressure_HMA(ibox,iblock) - write_mean(ii+1))**2
+            END DO
+            write_mean(ii+1) = write_mean(ii+1) * atomic_to_bar
+            write_err(ii+1)  = sqrt(write_err(ii+1) / nbr_blocks) * atomic_to_bar
+
          END IF
          
          ! At the end increment property counter by 1

From e0c393f0453cb6d2a104e530718d8105fce1378c Mon Sep 17 00:00:00 2001
From: "Andrew J. Schultz" <ajschult784@gmail.com>
Date: Fri, 26 Feb 2021 09:53:03 -0500
Subject: [PATCH 04/11] Use lattice-site cutoff when doing HMA

---
 Src/energy_routines.f90 | 72 +++++++++++++++++++++++++++++++++++------
 1 file changed, 62 insertions(+), 10 deletions(-)

diff --git a/Src/energy_routines.f90 b/Src/energy_routines.f90
index c8d96383..f6e83d2d 100755
--- a/Src/energy_routines.f90
+++ b/Src/energy_routines.f90
@@ -634,12 +634,22 @@ SUBROUTINE Compute_Atom_Nonbond_Energy(this_atom,this_molecule,this_species, &
                 ! Intermolecular energy so apply pbc.
 
                 ! First compute the parent separation
-                rxijp = atom_list(ia,im,is)%rxp &
-                      - atom_list(this_atom,this_molecule,this_species)%rxp
-                ryijp = atom_list(ia,im,is)%ryp &
-                      - atom_list(this_atom,this_molecule,this_species)%ryp
-                rzijp = atom_list(ia,im,is)%rzp &
-                      - atom_list(this_atom,this_molecule,this_species)%rzp
+                IF (need_HMA) THEN
+                   ! HMA uses lattice sites for cutoff
+                   rxijp = atom_list(ia,im,is)%rxp_init &
+                         - atom_list(this_atom,this_molecule,this_species)%rxp_init
+                   ryijp = atom_list(ia,im,is)%ryp_init &
+                         - atom_list(this_atom,this_molecule,this_species)%ryp_init
+                   rzijp = atom_list(ia,im,is)%rzp_init &
+                         - atom_list(this_atom,this_molecule,this_species)%rzp_init
+                ELSE
+                   rxijp = atom_list(ia,im,is)%rxp &
+                         - atom_list(this_atom,this_molecule,this_species)%rxp
+                   ryijp = atom_list(ia,im,is)%ryp &
+                         - atom_list(this_atom,this_molecule,this_species)%ryp
+                   rzijp = atom_list(ia,im,is)%rzp &
+                         - atom_list(this_atom,this_molecule,this_species)%rzp
+                END IF
 
                 ! Now get the minimum image separation
                 CALL Minimum_Image_Separation(this_box,rxijp,ryijp,rzijp, &
@@ -656,6 +666,21 @@ SUBROUTINE Compute_Atom_Nonbond_Energy(this_atom,this_molecule,this_species, &
 
              CALL Check_AtomPair_Cutoff(rijsq,get_vdw,get_qq,this_box)
 
+             IF (need_HMA) THEN
+                rxijp = atom_list(ia,im,is)%rxp &
+                      - atom_list(this_atom,this_molecule,this_species)%rxp
+                ryijp = atom_list(ia,im,is)%ryp &
+                      - atom_list(this_atom,this_molecule,this_species)%ryp
+                rzijp = atom_list(ia,im,is)%rzp &
+                      - atom_list(this_atom,this_molecule,this_species)%rzp
+
+                ! Now get the minimum image separation
+                CALL Minimum_Image_Separation(this_box,rxijp,ryijp,rzijp, &
+                     rxij,ryij,rzij)
+
+                rijsq = rxij*rxij + ryij*ryij + rzij*rzij
+             END IF
+
              ! Compute vdw and q-q energy using if required
              IF (get_vdw .OR. get_qq) THEN
 
@@ -2301,6 +2326,11 @@ SUBROUTINE Check_MoleculePair_Cutoff(im_1,is_1,im_2,is_2,get_interaction, &
 
     get_interaction = .TRUE.
 
+    IF (need_HMA) THEN
+       ! HMA uses an atomic cutoff on lattice sites
+       RETURN
+    END IF
+
     ! Figure out the box to be used later.
 
     this_box = molecule_list(im_1,is_1)%which_box
@@ -2481,7 +2511,9 @@ SUBROUTINE Compute_System_Total_Force(this_box)
    W_tensor_charge(:,:,this_box) = 0.0_DP
    W_tensor_recip(:,:,this_box) = 0.0_DP
    W_tensor_elec(:,:,this_box) =  0.0_DP
-   energy_HMA(this_box)%sum_Fdr = 0.0_DP
+   IF (need_HMA) THEN
+      energy_HMA(this_box)%sum_Fdr = 0.0_DP
+   END IF
 
    DO is = 1, nspecies
       imLOOP1: DO im_1 = 1, nmols(is,this_box)
@@ -2613,9 +2645,16 @@ SUBROUTINE Compute_MoleculePair_Force(im,is,jm,js,this_box,tens_vdw,tens_charge,
        DO ja = 1, natoms(js)
 
           ! Obtain the minimum image separation
-          rxijp = atom_list(ia,im,is)%rxp - atom_list(ja,jm,js)%rxp
-          ryijp = atom_list(ia,im,is)%ryp - atom_list(ja,jm,js)%ryp
-          rzijp = atom_list(ia,im,is)%rzp - atom_list(ja,jm,js)%rzp
+          IF (need_HMA) THEN
+             ! HMA uses lattice sites for cutoff
+             rxijp = atom_list(ia,im,is)%rxp_init - atom_list(ja,jm,js)%rxp_init
+             ryijp = atom_list(ia,im,is)%ryp_init - atom_list(ja,jm,js)%ryp_init
+             rzijp = atom_list(ia,im,is)%rzp_init - atom_list(ja,jm,js)%rzp_init
+          ELSE
+             rxijp = atom_list(ia,im,is)%rxp - atom_list(ja,jm,js)%rxp
+             ryijp = atom_list(ia,im,is)%ryp - atom_list(ja,jm,js)%ryp
+             rzijp = atom_list(ia,im,is)%rzp - atom_list(ja,jm,js)%rzp
+          END IF
 
           ! Now get the minimum image separation
           CALL Minimum_Image_Separation(this_box,rxijp,ryijp,rzijp, &
@@ -2626,6 +2665,19 @@ SUBROUTINE Compute_MoleculePair_Force(im,is,jm,js,this_box,tens_vdw,tens_charge,
           ! Now figure out what needs to be computed, then call pair_energy
           CALL Check_AtomPair_Cutoff(rijsq,get_vdw,get_qq,this_box)
 
+          IF (need_HMA) THEN
+             ! now compute actual atom-atom separation
+             rxijp = atom_list(ia,im,is)%rxp - atom_list(ja,jm,js)%rxp
+             ryijp = atom_list(ia,im,is)%ryp - atom_list(ja,jm,js)%ryp
+             rzijp = atom_list(ia,im,is)%rzp - atom_list(ja,jm,js)%rzp
+
+             ! Now get the minimum image separation
+             CALL Minimum_Image_Separation(this_box,rxijp,ryijp,rzijp, &
+                  rxij,ryij,rzij)
+
+             rijsq = rxij*rxij + ryij*ryij + rzij*rzij
+          END IF
+
           ! Compute vdw and q-q energy using if required
           IF (get_vdw .OR. get_qq) THEN
 

From 7c7f409ccc83765d0fb9dbb74687c57124153275 Mon Sep 17 00:00:00 2001
From: "Andrew J. Schultz" <ajschult784@gmail.com>
Date: Fri, 26 Feb 2021 10:42:02 -0500
Subject: [PATCH 05/11] Drop HMA section, just list harmonic pressure next to
 prop

---
 Src/input_routines.f90 | 144 ++++++++++-------------------------------
 Src/nvtmc_control.f90  |   2 -
 2 files changed, 35 insertions(+), 111 deletions(-)

diff --git a/Src/input_routines.f90 b/Src/input_routines.f90
index d5adf107..2dcab975 100755
--- a/Src/input_routines.f90
+++ b/Src/input_routines.f90
@@ -5984,6 +5984,37 @@ SUBROUTINE Get_Property_Info
               nbr_properties = nbr_properties + 1
               IF (int_sim_type /= sim_npt .AND. int_sim_type /= sim_gemc_npt) need_pressure = .TRUE.
 
+           ELSE IF (line_array(1) == 'energy_HMA' .OR. line_array(1) == 'Energy_HMA' .OR. &
+                    line_array(1) == 'pressure_HMA' .OR. line_array(1) == 'Pressure_HMA') THEN
+
+              IF (.NOT. need_HMA) THEN
+                 ALLOCATE(energy_HMA(nbr_boxes))
+                 ALLOCATE(pressure_HMA(nbr_boxes))
+
+                 energy_HMA(:)%lattice = 0.0_DP
+                 energy_HMA(:)%harmonic = 0.0_DP
+                 energy_HMA(:)%anharmonic = 0.0_DP
+                 energy_HMA(:)%total = 0.0_DP
+
+                 pressure_HMA(:)%lattice = 0.0_DP
+                 pressure_HMA(:)%harmonic = 0.0_DP
+                 pressure_HMA(:)%anharmonic = 0.0_DP
+                 pressure_HMA(:)%total = 0.0_DP
+              END IF
+
+              need_HMA = .TRUE.
+
+              IF (line_array(1) == 'pressure_HMA' .OR. line_array(1) == 'pressure_HMA') THEN
+                 ! harmonic pressure is the second value on the line
+                 pressure_HMA(this_box)%harmonic = String_To_Double(line_array(2))
+                 WRITE(logunit,'(A,X,F9.3,X,A)',ADVANCE='NO') 'HMA harmonic pressure of box', this_box, ' is', pressure_HMA(this_box)%harmonic,  'bar'
+
+                 ! convert pressure into atomic units
+                 pressure_HMA(this_box)%harmonic = pressure_HMA(this_box)%harmonic / atomic_to_bar
+                 WRITE(logunit,'(X,A,X,E13.6,X,A)') '=', pressure_HMA(this_box)%harmonic,  'amu / (A ps^2)'
+              END IF
+
+              nbr_properties = nbr_properties + 1
            ELSE
               ! this is a property for the system
               nbr_properties = nbr_properties + 1
@@ -6163,15 +6194,11 @@ SUBROUTINE Get_Property_Info
                  nbr_properties = nbr_properties + 1
                  prop_output(nbr_properties,nbr_prop_files(this_box),this_box) = 'Mass_Density'
               ELSE IF (line_array(1) == 'energy_HMA' .OR. line_array(1) == 'Energy_HMA') THEN
-                 IF (need_HMA) THEN
-                    nbr_properties = nbr_properties + 1
-                    prop_output(nbr_properties,nbr_prop_files(this_box),this_box) = 'Energy_HMA'
-                  END IF
+                 nbr_properties = nbr_properties + 1
+                 prop_output(nbr_properties,nbr_prop_files(this_box),this_box) = 'Energy_HMA'
               ELSE IF (line_array(1) == 'pressure_HMA' .OR. line_array(1) == 'Pressure_HMA') THEN
-                 IF (need_HMA) THEN
-                    nbr_properties = nbr_properties + 1
-                    prop_output(nbr_properties,nbr_prop_files(this_box),this_box) = 'Pressure_HMA'
-                 END IF
+                 nbr_properties = nbr_properties + 1
+                 prop_output(nbr_properties,nbr_prop_files(this_box),this_box) = 'Pressure_HMA'
               ELSE
                 err_msg = ''
                 err_msg(1) = 'Keyword "' // TRIM(line_array(1)) // '" on line ' // &
@@ -6713,105 +6740,4 @@ SUBROUTINE Get_Verbosity_Info
   WRITE(logunit,'(A80)') '********************************************************************************'
 
 END SUBROUTINE Get_Verbosity_Info
-
-
-!******************************************************************************
-SUBROUTINE Get_HMA_Info
-!******************************************************************************
-! This routine opens the input file and determines if HMA properties will be
-! computed and (optionally) the harmonic pressure.
-!
-! # HMA
-! TRUE
-! harmonic_pressure 9
-!
-! The routine searches for the keyword "# HMA" and then reads the necessary
-! information underneath the key word.
-!******************************************************************************
-
-  INTEGER :: ierr,line_nbr,nbr_entries
-  CHARACTER(STRING_LEN) :: line_string, line_array(60)
-
-!******************************************************************************
-  WRITE(logunit,*)
-  WRITE(logunit,'(A)') 'HMA'
-  WRITE(logunit,'(A80)') '********************************************************************************'
-
-  need_HMA = .FALSE.
-
-  REWIND(inputunit)
-
-  ierr = 0
-  line_nbr = 0
-
-  DO
-     line_nbr = line_nbr + 1
-     CALL Read_String(inputunit,line_string,ierr)
-
-     IF (ierr .NE. 0) THEN
-        EXIT
-     END IF
-
-     IF (line_string(1:13) == '# HMA') THEN
-        line_nbr = line_nbr + 1
-        CALL Parse_String(inputunit,line_nbr,1,nbr_entries,line_array,ierr)
-        IF (line_array(1) == 'TRUE' .OR. line_array(1) == '.TRUE.' .OR. &
-            line_array(1) == 'true' .OR. line_array(1) == '.true.') THEN
-           need_HMA = .TRUE.
-
-           ALLOCATE(energy_HMA(nbr_boxes))
-           ALLOCATE(pressure_HMA(nbr_boxes))
-
-           energy_HMA(:)%lattice = 0.0_DP
-           energy_HMA(:)%harmonic = 0.0_DP
-           energy_HMA(:)%anharmonic = 0.0_DP
-           energy_HMA(:)%total = 0.0_DP
-
-           pressure_HMA(:)%lattice = 0.0_DP
-           pressure_HMA(:)%harmonic = 0.0_DP
-           pressure_HMA(:)%anharmonic = 0.0_DP
-           pressure_HMA(:)%total = 0.0_DP
-
-        ELSE IF (line_array(1) == 'FALSE' .OR. line_array(1) == '.FALSE.' .OR. &
-                 line_array(1) == 'false' .OR. line_array(1) == '.false.') THEN
-           need_HMA = .FALSE.
-           EXIT
-        ELSE
-           err_msg = ''
-           err_msg(1) = 'Keyword ' // TRIM(line_array(1)) // ' on line number ' // &
-                        TRIM(Int_To_String(line_nbr)) // ' of the input file is not supported'
-           err_msg(2) = 'Supported keywords are: true, false'
-           CALL Clean_Abort(err_msg,'Get_Verbosity_Info')
-        END IF
-
-        line_nbr = line_nbr + 1
-        CALL Parse_String(inputunit,line_nbr,0,nbr_entries,line_array,ierr)
-
-        IF (nbr_entries == 0 .OR. line_array(1)(1:1) == '!') THEN
-           EXIT
-        END IF
-
-        IF (line_array(1) == 'harmonic_pressure' .OR. line_array(1) == 'Harmonic_Pressure') THEN
-           ! assign the pressures
-           pressure_HMA(:)%harmonic = String_To_Double(line_array(2))
-           WRITE(logunit,'(A,X,F9.3,X,A)',ADVANCE='NO') 'HMA harmonic pressure is', pressure_HMA(1)%harmonic,  'bar'
-
-           ! convert pressure into atomic units
-           pressure_HMA(:)%harmonic = pressure_HMA(:)%harmonic / atomic_to_bar
-           WRITE(logunit,'(X,A,X,E13.6,X,A)') '=', pressure_HMA(1)%harmonic,  'amu / (A ps^2)'
-        END IF
-
-        EXIT
-     ENDIF
-
-  ENDDO
-
-  IF (need_HMA) THEN
-     WRITE(logunit,'(A)') 'HMA enabled'
-  ELSE
-     WRITE(logunit,'(A)') 'HMA not enabled'
-  END IF
-  WRITE(logunit,'(A80)') '********************************************************************************'
-
-END SUBROUTINE Get_HMA_Info
 END MODULE Input_Routines
diff --git a/Src/nvtmc_control.f90 b/Src/nvtmc_control.f90
index a72a8f5c..769b5167 100755
--- a/Src/nvtmc_control.f90
+++ b/Src/nvtmc_control.f90
@@ -117,8 +117,6 @@ SUBROUTINE NVTMC_Control
   ! Write Widom insertion info to log file, which requires units from Get_Simulation_Length_Info
   CALL Log_Widom_Info
 
-  CALL Get_HMA_Info
-
   CALL Get_Property_Info
 
   CALL Get_Rcutoff_Low

From 0a9ea11d8909859ede74ee76c7d6d4e5cfab0505 Mon Sep 17 00:00:00 2001
From: "Andrew J. Schultz" <ajschult784@gmail.com>
Date: Fri, 26 Feb 2021 13:51:22 -0500
Subject: [PATCH 06/11] If we compute HMA pressure, we need standard pressure
 too

print box for each harmonic pressure
---
 Src/input_routines.f90 | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/Src/input_routines.f90 b/Src/input_routines.f90
index 2dcab975..c200290a 100755
--- a/Src/input_routines.f90
+++ b/Src/input_routines.f90
@@ -6005,9 +6005,10 @@ SUBROUTINE Get_Property_Info
               need_HMA = .TRUE.
 
               IF (line_array(1) == 'pressure_HMA' .OR. line_array(1) == 'pressure_HMA') THEN
+                 need_pressure = .TRUE.
                  ! harmonic pressure is the second value on the line
                  pressure_HMA(this_box)%harmonic = String_To_Double(line_array(2))
-                 WRITE(logunit,'(A,X,F9.3,X,A)',ADVANCE='NO') 'HMA harmonic pressure of box', this_box, ' is', pressure_HMA(this_box)%harmonic,  'bar'
+                 WRITE(logunit,'(A,X,I1,X,A,X,F9.3,X,A)',ADVANCE='NO') 'HMA harmonic pressure of box', this_box, 'is', pressure_HMA(this_box)%harmonic,  'bar'
 
                  ! convert pressure into atomic units
                  pressure_HMA(this_box)%harmonic = pressure_HMA(this_box)%harmonic / atomic_to_bar

From 486a974f1577bc350ee0b9486bacc33d0767bee8 Mon Sep 17 00:00:00 2001
From: "Andrew J. Schultz" <ajschult784@gmail.com>
Date: Fri, 26 Feb 2021 13:52:18 -0500
Subject: [PATCH 07/11] Fix lattice truncation

---
 Src/energy_routines.f90 | 79 ++++++++++++++++++++++++++---------------
 1 file changed, 50 insertions(+), 29 deletions(-)

diff --git a/Src/energy_routines.f90 b/Src/energy_routines.f90
index f6e83d2d..0a144e20 100755
--- a/Src/energy_routines.f90
+++ b/Src/energy_routines.f90
@@ -666,23 +666,23 @@ SUBROUTINE Compute_Atom_Nonbond_Energy(this_atom,this_molecule,this_species, &
 
              CALL Check_AtomPair_Cutoff(rijsq,get_vdw,get_qq,this_box)
 
-             IF (need_HMA) THEN
-                rxijp = atom_list(ia,im,is)%rxp &
-                      - atom_list(this_atom,this_molecule,this_species)%rxp
-                ryijp = atom_list(ia,im,is)%ryp &
-                      - atom_list(this_atom,this_molecule,this_species)%ryp
-                rzijp = atom_list(ia,im,is)%rzp &
-                      - atom_list(this_atom,this_molecule,this_species)%rzp
+             ! Compute vdw and q-q energy using if required
+             IF (get_vdw .OR. get_qq) THEN
+                IF (need_HMA) THEN
+                   rxijp = atom_list(ia,im,is)%rxp &
+                         - atom_list(this_atom,this_molecule,this_species)%rxp
+                   ryijp = atom_list(ia,im,is)%ryp &
+                         - atom_list(this_atom,this_molecule,this_species)%ryp
+                   rzijp = atom_list(ia,im,is)%rzp &
+                         - atom_list(this_atom,this_molecule,this_species)%rzp
 
-                ! Now get the minimum image separation
-                CALL Minimum_Image_Separation(this_box,rxijp,ryijp,rzijp, &
-                     rxij,ryij,rzij)
+                   ! Now get the minimum image separation
+                   CALL Minimum_Image_Separation(this_box,rxijp,ryijp,rzijp, &
+                        rxij,ryij,rzij)
 
-                rijsq = rxij*rxij + ryij*ryij + rzij*rzij
-             END IF
+                   rijsq = rxij*rxij + ryij*ryij + rzij*rzij
+                END IF
 
-             ! Compute vdw and q-q energy using if required
-             IF (get_vdw .OR. get_qq) THEN
 
                 CALL Compute_AtomPair_Energy(rxij,ryij,rzij,rijsq, &
                      is,im,ia,this_species,this_molecule,this_atom,&
@@ -1131,9 +1131,16 @@ SUBROUTINE Compute_MoleculePair_Energy(im,is,jm,js,this_box, &
         IF ( .NOT. atom_list(ja,jm,js)%exist) CYCLE
 
         ! Obtain the minimum image separation
-        rxijp = atom_list(ia,im,is)%rxp - atom_list(ja,jm,js)%rxp
-        ryijp = atom_list(ia,im,is)%ryp - atom_list(ja,jm,js)%ryp
-        rzijp = atom_list(ia,im,is)%rzp - atom_list(ja,jm,js)%rzp
+        IF (need_HMA) THEN
+          ! HMA uses lattice sites for cutoff
+          rxijp = atom_list(ia,im,is)%rxp_init - atom_list(ja,jm,js)%rxp_init
+          ryijp = atom_list(ia,im,is)%ryp_init - atom_list(ja,jm,js)%ryp_init
+          rzijp = atom_list(ia,im,is)%rzp_init - atom_list(ja,jm,js)%rzp_init
+        ELSE
+          rxijp = atom_list(ia,im,is)%rxp - atom_list(ja,jm,js)%rxp
+          ryijp = atom_list(ia,im,is)%ryp - atom_list(ja,jm,js)%ryp
+          rzijp = atom_list(ia,im,is)%rzp - atom_list(ja,jm,js)%rzp
+        END IF
 
         ! Now get the minimum image separation
         CALL Minimum_Image_Separation(this_box,rxijp,ryijp,rzijp,rxij,ryij,rzij)
@@ -1154,6 +1161,16 @@ SUBROUTINE Compute_MoleculePair_Energy(im,is,jm,js,this_box, &
 
         ! Compute vdw and q-q energy, if required
         IF (get_vdw .OR. get_qq) THEN
+          IF (need_HMA) THEN
+            ! HMA uses lattice sites for cutoff
+            rxijp = atom_list(ia,im,is)%rxp - atom_list(ja,jm,js)%rxp
+            ryijp = atom_list(ia,im,is)%ryp - atom_list(ja,jm,js)%ryp
+            rzijp = atom_list(ia,im,is)%rzp - atom_list(ja,jm,js)%rzp
+
+            CALL Minimum_Image_Separation(this_box,rxijp,ryijp,rzijp,rxij,ryij,rzij)
+
+            rijsq = rxij*rxij + ryij*ryij + rzij*rzij
+          END IF
 
           CALL Compute_AtomPair_Energy(rxij,ryij,rzij,rijsq, &
                is,im,ia,js,jm,ja,get_vdw,get_qq, &
@@ -2665,21 +2682,25 @@ SUBROUTINE Compute_MoleculePair_Force(im,is,jm,js,this_box,tens_vdw,tens_charge,
           ! Now figure out what needs to be computed, then call pair_energy
           CALL Check_AtomPair_Cutoff(rijsq,get_vdw,get_qq,this_box)
 
-          IF (need_HMA) THEN
-             ! now compute actual atom-atom separation
-             rxijp = atom_list(ia,im,is)%rxp - atom_list(ja,jm,js)%rxp
-             ryijp = atom_list(ia,im,is)%ryp - atom_list(ja,jm,js)%ryp
-             rzijp = atom_list(ia,im,is)%rzp - atom_list(ja,jm,js)%rzp
+          ! Compute vdw and q-q energy using if required
+          IF (get_vdw .OR. get_qq) THEN
 
-             ! Now get the minimum image separation
-             CALL Minimum_Image_Separation(this_box,rxijp,ryijp,rzijp, &
-                  rxij,ryij,rzij)
+             IF (need_HMA) THEN
+                ! now compute actual atom-atom separation
+                rxijp = atom_list(ia,im,is)%rxp - atom_list(ja,jm,js)%rxp
+                ryijp = atom_list(ia,im,is)%ryp - atom_list(ja,jm,js)%ryp
+                rzijp = atom_list(ia,im,is)%rzp - atom_list(ja,jm,js)%rzp
 
-             rijsq = rxij*rxij + ryij*ryij + rzij*rzij
-          END IF
+                ! Now get the minimum image separation
+                CALL Minimum_Image_Separation(this_box,rxijp,ryijp,rzijp, &
+                     rxij,ryij,rzij)
+
+                rabx = rxij
+                raby = ryij
+                rabz = rzij
+                rijsq = rxij*rxij + ryij*ryij + rzij*rzij
+             END IF
 
-          ! Compute vdw and q-q energy using if required
-          IF (get_vdw .OR. get_qq) THEN
 
              CALL Compute_AtomPair_Force(rijsq,is,im,ia,js,jm,ja,&
                   get_vdw,get_qq,Wij_vdw,Wij_qq)

From 53a3c9912b2442f371f477dfc3cf9cd2c4f50b0b Mon Sep 17 00:00:00 2001
From: "Andrew J. Schultz" <ajschult784@gmail.com>
Date: Wed, 3 Mar 2021 14:16:29 -0500
Subject: [PATCH 08/11] Move lattice cutoff into its own method

---
 Src/energy_routines.f90 | 143 +++++++++++++++++++---------------------
 1 file changed, 68 insertions(+), 75 deletions(-)

diff --git a/Src/energy_routines.f90 b/Src/energy_routines.f90
index 0a144e20..8c0b378a 100755
--- a/Src/energy_routines.f90
+++ b/Src/energy_routines.f90
@@ -134,6 +134,11 @@ MODULE Energy_Routines
   !
   ! Check_AtomPair_Cutoff:
   !
+  ! Check_AtomPair_Lattice_Cutoff: Uses lattice sites to determine truncation
+  !                       instead of molecules (or atoms) so that atoms do not
+  !                       fluctuate in and out of the truncation during the
+  !                       simulation.
+  !
   ! Get_Molecule_Energy: Computes the intra- and inter-molecular energy of
   !                       a given molecule interacting with all other molecules.
   !
@@ -634,22 +639,12 @@ SUBROUTINE Compute_Atom_Nonbond_Energy(this_atom,this_molecule,this_species, &
                 ! Intermolecular energy so apply pbc.
 
                 ! First compute the parent separation
-                IF (need_HMA) THEN
-                   ! HMA uses lattice sites for cutoff
-                   rxijp = atom_list(ia,im,is)%rxp_init &
-                         - atom_list(this_atom,this_molecule,this_species)%rxp_init
-                   ryijp = atom_list(ia,im,is)%ryp_init &
-                         - atom_list(this_atom,this_molecule,this_species)%ryp_init
-                   rzijp = atom_list(ia,im,is)%rzp_init &
-                         - atom_list(this_atom,this_molecule,this_species)%rzp_init
-                ELSE
-                   rxijp = atom_list(ia,im,is)%rxp &
-                         - atom_list(this_atom,this_molecule,this_species)%rxp
-                   ryijp = atom_list(ia,im,is)%ryp &
-                         - atom_list(this_atom,this_molecule,this_species)%ryp
-                   rzijp = atom_list(ia,im,is)%rzp &
-                         - atom_list(this_atom,this_molecule,this_species)%rzp
-                END IF
+                rxijp = atom_list(ia,im,is)%rxp &
+                      - atom_list(this_atom,this_molecule,this_species)%rxp
+                ryijp = atom_list(ia,im,is)%ryp &
+                      - atom_list(this_atom,this_molecule,this_species)%ryp
+                rzijp = atom_list(ia,im,is)%rzp &
+                      - atom_list(this_atom,this_molecule,this_species)%rzp
 
                 ! Now get the minimum image separation
                 CALL Minimum_Image_Separation(this_box,rxijp,ryijp,rzijp, &
@@ -664,25 +659,16 @@ SUBROUTINE Compute_Atom_Nonbond_Energy(this_atom,this_molecule,this_species, &
 
              ENDIF
 
-             CALL Check_AtomPair_Cutoff(rijsq,get_vdw,get_qq,this_box)
+             IF (need_HMA) THEN
+                CALL Check_AtomPair_Lattice_Cutoff(atom_list(ia,im,is), &
+                        atom_list(this_atom,this_molecule,this_species), &
+                        get_vdw,get_qq,this_box)
+             ELSE
+                CALL Check_AtomPair_Cutoff(rijsq,get_vdw,get_qq,this_box)
+             END IF
 
              ! Compute vdw and q-q energy using if required
              IF (get_vdw .OR. get_qq) THEN
-                IF (need_HMA) THEN
-                   rxijp = atom_list(ia,im,is)%rxp &
-                         - atom_list(this_atom,this_molecule,this_species)%rxp
-                   ryijp = atom_list(ia,im,is)%ryp &
-                         - atom_list(this_atom,this_molecule,this_species)%ryp
-                   rzijp = atom_list(ia,im,is)%rzp &
-                         - atom_list(this_atom,this_molecule,this_species)%rzp
-
-                   ! Now get the minimum image separation
-                   CALL Minimum_Image_Separation(this_box,rxijp,ryijp,rzijp, &
-                        rxij,ryij,rzij)
-
-                   rijsq = rxij*rxij + ryij*ryij + rzij*rzij
-                END IF
-
 
                 CALL Compute_AtomPair_Energy(rxij,ryij,rzij,rijsq, &
                      is,im,ia,this_species,this_molecule,this_atom,&
@@ -1131,16 +1117,9 @@ SUBROUTINE Compute_MoleculePair_Energy(im,is,jm,js,this_box, &
         IF ( .NOT. atom_list(ja,jm,js)%exist) CYCLE
 
         ! Obtain the minimum image separation
-        IF (need_HMA) THEN
-          ! HMA uses lattice sites for cutoff
-          rxijp = atom_list(ia,im,is)%rxp_init - atom_list(ja,jm,js)%rxp_init
-          ryijp = atom_list(ia,im,is)%ryp_init - atom_list(ja,jm,js)%ryp_init
-          rzijp = atom_list(ia,im,is)%rzp_init - atom_list(ja,jm,js)%rzp_init
-        ELSE
-          rxijp = atom_list(ia,im,is)%rxp - atom_list(ja,jm,js)%rxp
-          ryijp = atom_list(ia,im,is)%ryp - atom_list(ja,jm,js)%ryp
-          rzijp = atom_list(ia,im,is)%rzp - atom_list(ja,jm,js)%rzp
-        END IF
+        rxijp = atom_list(ia,im,is)%rxp - atom_list(ja,jm,js)%rxp
+        ryijp = atom_list(ia,im,is)%ryp - atom_list(ja,jm,js)%ryp
+        rzijp = atom_list(ia,im,is)%rzp - atom_list(ja,jm,js)%rzp
 
         ! Now get the minimum image separation
         CALL Minimum_Image_Separation(this_box,rxijp,ryijp,rzijp,rxij,ryij,rzij)
@@ -1153,7 +1132,12 @@ SUBROUTINE Compute_MoleculePair_Energy(im,is,jm,js,this_box, &
         END IF
 
         ! Now figure out what needs to be computed, then call pair_energy
-        CALL Check_AtomPair_Cutoff(rijsq,get_vdw,get_qq,this_box)
+        IF (need_HMA) THEN
+           CALL Check_AtomPair_Lattice_Cutoff(atom_list(ia,im,is), &
+                   atom_list(ja,jm,js),get_vdw,get_qq,this_box)
+        ELSE
+           CALL Check_AtomPair_Cutoff(rijsq,get_vdw,get_qq,this_box)
+        END IF
 
         IF(cbmc_flag .AND. (.NOT. species_list(is)%L_Coul_CBMC)) THEN
           get_qq=.FALSE.
@@ -1161,16 +1145,6 @@ SUBROUTINE Compute_MoleculePair_Energy(im,is,jm,js,this_box, &
 
         ! Compute vdw and q-q energy, if required
         IF (get_vdw .OR. get_qq) THEN
-          IF (need_HMA) THEN
-            ! HMA uses lattice sites for cutoff
-            rxijp = atom_list(ia,im,is)%rxp - atom_list(ja,jm,js)%rxp
-            ryijp = atom_list(ia,im,is)%ryp - atom_list(ja,jm,js)%ryp
-            rzijp = atom_list(ia,im,is)%rzp - atom_list(ja,jm,js)%rzp
-
-            CALL Minimum_Image_Separation(this_box,rxijp,ryijp,rzijp,rxij,ryij,rzij)
-
-            rijsq = rxij*rxij + ryij*ryij + rzij*rzij
-          END IF
 
           CALL Compute_AtomPair_Energy(rxij,ryij,rzij,rijsq, &
                is,im,ia,js,jm,ja,get_vdw,get_qq, &
@@ -2493,6 +2467,36 @@ SUBROUTINE Check_AtomPair_Cutoff(rijsq,get_vdw,get_qq,this_box)
 
  END SUBROUTINE Check_AtomPair_Cutoff
 
+!*******************************************************************************
+ SUBROUTINE Check_AtomPair_Lattice_Cutoff(atom1,atom2,get_vdw,get_qq,this_box)
+   !***************************************************************************
+   ! The subroutine determines if the interaction between atom1 and atom2
+   ! should be truncated or not based on the distance between their lattice
+   ! sites.  This routine computes the lattice distance and uses
+   ! Check_AtomPair_Cutoff to decide on the actual truncation.
+   !***************************************************************************
+   TYPE(Atom_Class) :: atom1, atom2
+   INTEGER          :: this_box
+   LOGICAL          :: get_vdw, get_qq
+
+   REAL(DP) :: rxijp, ryijp, rzijp, rxij, ryij, rzij, rijsq
+
+   rxijp = atom1%rxp_init - atom2%rxp_init
+   ryijp = atom1%ryp_init - atom2%ryp_init
+   rzijp = atom1%rzp_init - atom2%rzp_init
+
+   ! Now get the minimum image separation
+   CALL Minimum_Image_Separation(this_box,rxijp,ryijp,rzijp,rxij,ryij,rzij)
+
+   rijsq = rxij*rxij + ryij*ryij + rzij*rzij
+
+   ! Now figure out what needs to be computed, then call pair_energy
+   CALL Check_AtomPair_Cutoff(rijsq,get_vdw,get_qq,this_box)
+
+   RETURN
+
+ END SUBROUTINE Check_AtomPair_Lattice_Cutoff
+
  SUBROUTINE Compute_System_Total_Force(this_box)
 
    !****************************************************************************
@@ -2662,16 +2666,9 @@ SUBROUTINE Compute_MoleculePair_Force(im,is,jm,js,this_box,tens_vdw,tens_charge,
        DO ja = 1, natoms(js)
 
           ! Obtain the minimum image separation
-          IF (need_HMA) THEN
-             ! HMA uses lattice sites for cutoff
-             rxijp = atom_list(ia,im,is)%rxp_init - atom_list(ja,jm,js)%rxp_init
-             ryijp = atom_list(ia,im,is)%ryp_init - atom_list(ja,jm,js)%ryp_init
-             rzijp = atom_list(ia,im,is)%rzp_init - atom_list(ja,jm,js)%rzp_init
-          ELSE
-             rxijp = atom_list(ia,im,is)%rxp - atom_list(ja,jm,js)%rxp
-             ryijp = atom_list(ia,im,is)%ryp - atom_list(ja,jm,js)%ryp
-             rzijp = atom_list(ia,im,is)%rzp - atom_list(ja,jm,js)%rzp
-          END IF
+          rxijp = atom_list(ia,im,is)%rxp - atom_list(ja,jm,js)%rxp
+          ryijp = atom_list(ia,im,is)%ryp - atom_list(ja,jm,js)%ryp
+          rzijp = atom_list(ia,im,is)%rzp - atom_list(ja,jm,js)%rzp
 
           ! Now get the minimum image separation
           CALL Minimum_Image_Separation(this_box,rxijp,ryijp,rzijp, &
@@ -2680,28 +2677,24 @@ SUBROUTINE Compute_MoleculePair_Force(im,is,jm,js,this_box,tens_vdw,tens_charge,
           rijsq = rxij*rxij + ryij*ryij + rzij*rzij
 
           ! Now figure out what needs to be computed, then call pair_energy
-          CALL Check_AtomPair_Cutoff(rijsq,get_vdw,get_qq,this_box)
+          IF (need_HMA) THEN
+             CALL Check_AtomPair_Lattice_Cutoff(atom_list(ia,im,is), &
+                     atom_list(ja,jm,js),get_vdw,get_qq,this_box)
+          ELSE
+             CALL Check_AtomPair_Cutoff(rijsq,get_vdw,get_qq,this_box)
+          END IF
 
           ! Compute vdw and q-q energy using if required
           IF (get_vdw .OR. get_qq) THEN
 
              IF (need_HMA) THEN
-                ! now compute actual atom-atom separation
-                rxijp = atom_list(ia,im,is)%rxp - atom_list(ja,jm,js)%rxp
-                ryijp = atom_list(ia,im,is)%ryp - atom_list(ja,jm,js)%ryp
-                rzijp = atom_list(ia,im,is)%rzp - atom_list(ja,jm,js)%rzp
-
-                ! Now get the minimum image separation
-                CALL Minimum_Image_Separation(this_box,rxijp,ryijp,rzijp, &
-                     rxij,ryij,rzij)
-
+                ! molecule-molecule vector was not computed before, but we
+                ! know the molecules have only one atom each
                 rabx = rxij
                 raby = ryij
                 rabz = rzij
-                rijsq = rxij*rxij + ryij*ryij + rzij*rzij
              END IF
 
-
              CALL Compute_AtomPair_Force(rijsq,is,im,ia,js,jm,ja,&
                   get_vdw,get_qq,Wij_vdw,Wij_qq)
 

From 6691feae4073da7ed4a67be0ce6ad74f4381baa8 Mon Sep 17 00:00:00 2001
From: "Andrew J. Schultz" <ajschult784@gmail.com>
Date: Tue, 18 May 2021 12:52:09 -0400
Subject: [PATCH 09/11] Add HMA example

---
 Examples/NVT/HMA/README     |  19 ++
 Examples/NVT/HMA/fcc500.xyz | 502 ++++++++++++++++++++++++++++++++++++
 Examples/NVT/HMA/inp500_1.0 |  79 ++++++
 Examples/NVT/HMA/lj.mcf     |  34 +++
 4 files changed, 634 insertions(+)
 create mode 100644 Examples/NVT/HMA/README
 create mode 100644 Examples/NVT/HMA/fcc500.xyz
 create mode 100644 Examples/NVT/HMA/inp500_1.0
 create mode 100644 Examples/NVT/HMA/lj.mcf

diff --git a/Examples/NVT/HMA/README b/Examples/NVT/HMA/README
new file mode 100644
index 00000000..8a31e420
--- /dev/null
+++ b/Examples/NVT/HMA/README
@@ -0,0 +1,19 @@
+This directory contains example HMA input files to run a 500-atom FCC
+simulation using harmonically mapped averaging (HMA) with a Lennard-Jones
+potential.
+
+More information about HMA is available in these publications:
+
+A. J. Schultz, D. A. Kofke, “Comprehensive high-precision high-accuracy
+equation of state and coexistence properties for classical Lennard-Jones
+crystals and low-temperature fluid phases”, J. Chem. Phys. 149, 204508 (2018)
+https://doi.org/10.1063/1.5053714
+
+S. G. Moustafa, A. J. Schultz, D. A. Kofke, “Harmonically Assisted Methods for
+Computing the Free Energy of Classical Crystals by Molecular Simulation: A
+Comparative Study”, J. Chem. Theory Comput. 13, 825-834 (2017)
+https://doi.org/10.1021/acs.jctc.6b01082
+
+S. G. Moustafa, A. J. Schultz, D. A. Kofke, “Very fast averaging of thermal
+properties of crystals by molecular simulation”, Phys. Rev. E 92, 043303 (2015)
+https://doi.org/10.1103/PhysRevE.92.043303
diff --git a/Examples/NVT/HMA/fcc500.xyz b/Examples/NVT/HMA/fcc500.xyz
new file mode 100644
index 00000000..bd43c226
--- /dev/null
+++ b/Examples/NVT/HMA/fcc500.xyz
@@ -0,0 +1,502 @@
+500
+# lattice
+A   -3.571652366928449   -3.571652366928449   -3.571652366928449
+A   -2.777951840944349   -2.777951840944349   -3.571652366928449
+A   -2.777951840944349   -3.571652366928449   -2.777951840944349
+A   -3.571652366928449   -2.777951840944349   -2.777951840944349
+A   -3.571652366928449   -3.571652366928449   -1.984251314960249
+A   -2.777951840944349   -2.777951840944349   -1.984251314960249
+A   -2.777951840944349   -3.571652366928449   -1.190550788976150
+A   -3.571652366928449   -2.777951840944349   -1.190550788976150
+A   -3.571652366928449   -3.571652366928449   -0.396850262992050
+A   -2.777951840944349   -2.777951840944349   -0.396850262992050
+A   -2.777951840944349   -3.571652366928449    0.396850262992050
+A   -3.571652366928449   -2.777951840944349    0.396850262992050
+A   -3.571652366928449   -3.571652366928449    1.190550788976150
+A   -2.777951840944349   -2.777951840944349    1.190550788976150
+A   -2.777951840944349   -3.571652366928449    1.984251314960249
+A   -3.571652366928449   -2.777951840944349    1.984251314960249
+A   -3.571652366928449   -3.571652366928449    2.777951840944349
+A   -2.777951840944349   -2.777951840944349    2.777951840944349
+A   -2.777951840944349   -3.571652366928449    3.571652366928448
+A   -3.571652366928449   -2.777951840944349    3.571652366928448
+A   -3.571652366928449   -1.984251314960249   -3.571652366928449
+A   -2.777951840944349   -1.190550788976150   -3.571652366928449
+A   -2.777951840944349   -1.984251314960249   -2.777951840944349
+A   -3.571652366928449   -1.190550788976150   -2.777951840944349
+A   -3.571652366928449   -1.984251314960249   -1.984251314960249
+A   -2.777951840944349   -1.190550788976150   -1.984251314960249
+A   -2.777951840944349   -1.984251314960249   -1.190550788976150
+A   -3.571652366928449   -1.190550788976150   -1.190550788976150
+A   -3.571652366928449   -1.984251314960249   -0.396850262992050
+A   -2.777951840944349   -1.190550788976150   -0.396850262992050
+A   -2.777951840944349   -1.984251314960249    0.396850262992050
+A   -3.571652366928449   -1.190550788976150    0.396850262992050
+A   -3.571652366928449   -1.984251314960249    1.190550788976150
+A   -2.777951840944349   -1.190550788976150    1.190550788976150
+A   -2.777951840944349   -1.984251314960249    1.984251314960249
+A   -3.571652366928449   -1.190550788976150    1.984251314960249
+A   -3.571652366928449   -1.984251314960249    2.777951840944349
+A   -2.777951840944349   -1.190550788976150    2.777951840944349
+A   -2.777951840944349   -1.984251314960249    3.571652366928448
+A   -3.571652366928449   -1.190550788976150    3.571652366928448
+A   -3.571652366928449   -0.396850262992050   -3.571652366928449
+A   -2.777951840944349    0.396850262992050   -3.571652366928449
+A   -2.777951840944349   -0.396850262992050   -2.777951840944349
+A   -3.571652366928449    0.396850262992050   -2.777951840944349
+A   -3.571652366928449   -0.396850262992050   -1.984251314960249
+A   -2.777951840944349    0.396850262992050   -1.984251314960249
+A   -2.777951840944349   -0.396850262992050   -1.190550788976150
+A   -3.571652366928449    0.396850262992050   -1.190550788976150
+A   -3.571652366928449   -0.396850262992050   -0.396850262992050
+A   -2.777951840944349    0.396850262992050   -0.396850262992050
+A   -2.777951840944349   -0.396850262992050    0.396850262992050
+A   -3.571652366928449    0.396850262992050    0.396850262992050
+A   -3.571652366928449   -0.396850262992050    1.190550788976150
+A   -2.777951840944349    0.396850262992050    1.190550788976150
+A   -2.777951840944349   -0.396850262992050    1.984251314960249
+A   -3.571652366928449    0.396850262992050    1.984251314960249
+A   -3.571652366928449   -0.396850262992050    2.777951840944349
+A   -2.777951840944349    0.396850262992050    2.777951840944349
+A   -2.777951840944349   -0.396850262992050    3.571652366928448
+A   -3.571652366928449    0.396850262992050    3.571652366928448
+A   -3.571652366928449    1.190550788976150   -3.571652366928449
+A   -2.777951840944349    1.984251314960249   -3.571652366928449
+A   -2.777951840944349    1.190550788976150   -2.777951840944349
+A   -3.571652366928449    1.984251314960249   -2.777951840944349
+A   -3.571652366928449    1.190550788976150   -1.984251314960249
+A   -2.777951840944349    1.984251314960249   -1.984251314960249
+A   -2.777951840944349    1.190550788976150   -1.190550788976150
+A   -3.571652366928449    1.984251314960249   -1.190550788976150
+A   -3.571652366928449    1.190550788976150   -0.396850262992050
+A   -2.777951840944349    1.984251314960249   -0.396850262992050
+A   -2.777951840944349    1.190550788976150    0.396850262992050
+A   -3.571652366928449    1.984251314960249    0.396850262992050
+A   -3.571652366928449    1.190550788976150    1.190550788976150
+A   -2.777951840944349    1.984251314960249    1.190550788976150
+A   -2.777951840944349    1.190550788976150    1.984251314960249
+A   -3.571652366928449    1.984251314960249    1.984251314960249
+A   -3.571652366928449    1.190550788976150    2.777951840944349
+A   -2.777951840944349    1.984251314960249    2.777951840944349
+A   -2.777951840944349    1.190550788976150    3.571652366928448
+A   -3.571652366928449    1.984251314960249    3.571652366928448
+A   -3.571652366928449    2.777951840944349   -3.571652366928449
+A   -2.777951840944349    3.571652366928448   -3.571652366928449
+A   -2.777951840944349    2.777951840944349   -2.777951840944349
+A   -3.571652366928449    3.571652366928448   -2.777951840944349
+A   -3.571652366928449    2.777951840944349   -1.984251314960249
+A   -2.777951840944349    3.571652366928448   -1.984251314960249
+A   -2.777951840944349    2.777951840944349   -1.190550788976150
+A   -3.571652366928449    3.571652366928448   -1.190550788976150
+A   -3.571652366928449    2.777951840944349   -0.396850262992050
+A   -2.777951840944349    3.571652366928448   -0.396850262992050
+A   -2.777951840944349    2.777951840944349    0.396850262992050
+A   -3.571652366928449    3.571652366928448    0.396850262992050
+A   -3.571652366928449    2.777951840944349    1.190550788976150
+A   -2.777951840944349    3.571652366928448    1.190550788976150
+A   -2.777951840944349    2.777951840944349    1.984251314960249
+A   -3.571652366928449    3.571652366928448    1.984251314960249
+A   -3.571652366928449    2.777951840944349    2.777951840944349
+A   -2.777951840944349    3.571652366928448    2.777951840944349
+A   -2.777951840944349    2.777951840944349    3.571652366928448
+A   -3.571652366928449    3.571652366928448    3.571652366928448
+A   -1.984251314960249   -3.571652366928449   -3.571652366928449
+A   -1.190550788976150   -2.777951840944349   -3.571652366928449
+A   -1.190550788976150   -3.571652366928449   -2.777951840944349
+A   -1.984251314960249   -2.777951840944349   -2.777951840944349
+A   -1.984251314960249   -3.571652366928449   -1.984251314960249
+A   -1.190550788976150   -2.777951840944349   -1.984251314960249
+A   -1.190550788976150   -3.571652366928449   -1.190550788976150
+A   -1.984251314960249   -2.777951840944349   -1.190550788976150
+A   -1.984251314960249   -3.571652366928449   -0.396850262992050
+A   -1.190550788976150   -2.777951840944349   -0.396850262992050
+A   -1.190550788976150   -3.571652366928449    0.396850262992050
+A   -1.984251314960249   -2.777951840944349    0.396850262992050
+A   -1.984251314960249   -3.571652366928449    1.190550788976150
+A   -1.190550788976150   -2.777951840944349    1.190550788976150
+A   -1.190550788976150   -3.571652366928449    1.984251314960249
+A   -1.984251314960249   -2.777951840944349    1.984251314960249
+A   -1.984251314960249   -3.571652366928449    2.777951840944349
+A   -1.190550788976150   -2.777951840944349    2.777951840944349
+A   -1.190550788976150   -3.571652366928449    3.571652366928448
+A   -1.984251314960249   -2.777951840944349    3.571652366928448
+A   -1.984251314960249   -1.984251314960249   -3.571652366928449
+A   -1.190550788976150   -1.190550788976150   -3.571652366928449
+A   -1.190550788976150   -1.984251314960249   -2.777951840944349
+A   -1.984251314960249   -1.190550788976150   -2.777951840944349
+A   -1.984251314960249   -1.984251314960249   -1.984251314960249
+A   -1.190550788976150   -1.190550788976150   -1.984251314960249
+A   -1.190550788976150   -1.984251314960249   -1.190550788976150
+A   -1.984251314960249   -1.190550788976150   -1.190550788976150
+A   -1.984251314960249   -1.984251314960249   -0.396850262992050
+A   -1.190550788976150   -1.190550788976150   -0.396850262992050
+A   -1.190550788976150   -1.984251314960249    0.396850262992050
+A   -1.984251314960249   -1.190550788976150    0.396850262992050
+A   -1.984251314960249   -1.984251314960249    1.190550788976150
+A   -1.190550788976150   -1.190550788976150    1.190550788976150
+A   -1.190550788976150   -1.984251314960249    1.984251314960249
+A   -1.984251314960249   -1.190550788976150    1.984251314960249
+A   -1.984251314960249   -1.984251314960249    2.777951840944349
+A   -1.190550788976150   -1.190550788976150    2.777951840944349
+A   -1.190550788976150   -1.984251314960249    3.571652366928448
+A   -1.984251314960249   -1.190550788976150    3.571652366928448
+A   -1.984251314960249   -0.396850262992050   -3.571652366928449
+A   -1.190550788976150    0.396850262992050   -3.571652366928449
+A   -1.190550788976150   -0.396850262992050   -2.777951840944349
+A   -1.984251314960249    0.396850262992050   -2.777951840944349
+A   -1.984251314960249   -0.396850262992050   -1.984251314960249
+A   -1.190550788976150    0.396850262992050   -1.984251314960249
+A   -1.190550788976150   -0.396850262992050   -1.190550788976150
+A   -1.984251314960249    0.396850262992050   -1.190550788976150
+A   -1.984251314960249   -0.396850262992050   -0.396850262992050
+A   -1.190550788976150    0.396850262992050   -0.396850262992050
+A   -1.190550788976150   -0.396850262992050    0.396850262992050
+A   -1.984251314960249    0.396850262992050    0.396850262992050
+A   -1.984251314960249   -0.396850262992050    1.190550788976150
+A   -1.190550788976150    0.396850262992050    1.190550788976150
+A   -1.190550788976150   -0.396850262992050    1.984251314960249
+A   -1.984251314960249    0.396850262992050    1.984251314960249
+A   -1.984251314960249   -0.396850262992050    2.777951840944349
+A   -1.190550788976150    0.396850262992050    2.777951840944349
+A   -1.190550788976150   -0.396850262992050    3.571652366928448
+A   -1.984251314960249    0.396850262992050    3.571652366928448
+A   -1.984251314960249    1.190550788976150   -3.571652366928449
+A   -1.190550788976150    1.984251314960249   -3.571652366928449
+A   -1.190550788976150    1.190550788976150   -2.777951840944349
+A   -1.984251314960249    1.984251314960249   -2.777951840944349
+A   -1.984251314960249    1.190550788976150   -1.984251314960249
+A   -1.190550788976150    1.984251314960249   -1.984251314960249
+A   -1.190550788976150    1.190550788976150   -1.190550788976150
+A   -1.984251314960249    1.984251314960249   -1.190550788976150
+A   -1.984251314960249    1.190550788976150   -0.396850262992050
+A   -1.190550788976150    1.984251314960249   -0.396850262992050
+A   -1.190550788976150    1.190550788976150    0.396850262992050
+A   -1.984251314960249    1.984251314960249    0.396850262992050
+A   -1.984251314960249    1.190550788976150    1.190550788976150
+A   -1.190550788976150    1.984251314960249    1.190550788976150
+A   -1.190550788976150    1.190550788976150    1.984251314960249
+A   -1.984251314960249    1.984251314960249    1.984251314960249
+A   -1.984251314960249    1.190550788976150    2.777951840944349
+A   -1.190550788976150    1.984251314960249    2.777951840944349
+A   -1.190550788976150    1.190550788976150    3.571652366928448
+A   -1.984251314960249    1.984251314960249    3.571652366928448
+A   -1.984251314960249    2.777951840944349   -3.571652366928449
+A   -1.190550788976150    3.571652366928448   -3.571652366928449
+A   -1.190550788976150    2.777951840944349   -2.777951840944349
+A   -1.984251314960249    3.571652366928448   -2.777951840944349
+A   -1.984251314960249    2.777951840944349   -1.984251314960249
+A   -1.190550788976150    3.571652366928448   -1.984251314960249
+A   -1.190550788976150    2.777951840944349   -1.190550788976150
+A   -1.984251314960249    3.571652366928448   -1.190550788976150
+A   -1.984251314960249    2.777951840944349   -0.396850262992050
+A   -1.190550788976150    3.571652366928448   -0.396850262992050
+A   -1.190550788976150    2.777951840944349    0.396850262992050
+A   -1.984251314960249    3.571652366928448    0.396850262992050
+A   -1.984251314960249    2.777951840944349    1.190550788976150
+A   -1.190550788976150    3.571652366928448    1.190550788976150
+A   -1.190550788976150    2.777951840944349    1.984251314960249
+A   -1.984251314960249    3.571652366928448    1.984251314960249
+A   -1.984251314960249    2.777951840944349    2.777951840944349
+A   -1.190550788976150    3.571652366928448    2.777951840944349
+A   -1.190550788976150    2.777951840944349    3.571652366928448
+A   -1.984251314960249    3.571652366928448    3.571652366928448
+A   -0.396850262992050   -3.571652366928449   -3.571652366928449
+A    0.396850262992050   -2.777951840944349   -3.571652366928449
+A    0.396850262992050   -3.571652366928449   -2.777951840944349
+A   -0.396850262992050   -2.777951840944349   -2.777951840944349
+A   -0.396850262992050   -3.571652366928449   -1.984251314960249
+A    0.396850262992050   -2.777951840944349   -1.984251314960249
+A    0.396850262992050   -3.571652366928449   -1.190550788976150
+A   -0.396850262992050   -2.777951840944349   -1.190550788976150
+A   -0.396850262992050   -3.571652366928449   -0.396850262992050
+A    0.396850262992050   -2.777951840944349   -0.396850262992050
+A    0.396850262992050   -3.571652366928449    0.396850262992050
+A   -0.396850262992050   -2.777951840944349    0.396850262992050
+A   -0.396850262992050   -3.571652366928449    1.190550788976150
+A    0.396850262992050   -2.777951840944349    1.190550788976150
+A    0.396850262992050   -3.571652366928449    1.984251314960249
+A   -0.396850262992050   -2.777951840944349    1.984251314960249
+A   -0.396850262992050   -3.571652366928449    2.777951840944349
+A    0.396850262992050   -2.777951840944349    2.777951840944349
+A    0.396850262992050   -3.571652366928449    3.571652366928448
+A   -0.396850262992050   -2.777951840944349    3.571652366928448
+A   -0.396850262992050   -1.984251314960249   -3.571652366928449
+A    0.396850262992050   -1.190550788976150   -3.571652366928449
+A    0.396850262992050   -1.984251314960249   -2.777951840944349
+A   -0.396850262992050   -1.190550788976150   -2.777951840944349
+A   -0.396850262992050   -1.984251314960249   -1.984251314960249
+A    0.396850262992050   -1.190550788976150   -1.984251314960249
+A    0.396850262992050   -1.984251314960249   -1.190550788976150
+A   -0.396850262992050   -1.190550788976150   -1.190550788976150
+A   -0.396850262992050   -1.984251314960249   -0.396850262992050
+A    0.396850262992050   -1.190550788976150   -0.396850262992050
+A    0.396850262992050   -1.984251314960249    0.396850262992050
+A   -0.396850262992050   -1.190550788976150    0.396850262992050
+A   -0.396850262992050   -1.984251314960249    1.190550788976150
+A    0.396850262992050   -1.190550788976150    1.190550788976150
+A    0.396850262992050   -1.984251314960249    1.984251314960249
+A   -0.396850262992050   -1.190550788976150    1.984251314960249
+A   -0.396850262992050   -1.984251314960249    2.777951840944349
+A    0.396850262992050   -1.190550788976150    2.777951840944349
+A    0.396850262992050   -1.984251314960249    3.571652366928448
+A   -0.396850262992050   -1.190550788976150    3.571652366928448
+A   -0.396850262992050   -0.396850262992050   -3.571652366928449
+A    0.396850262992050    0.396850262992050   -3.571652366928449
+A    0.396850262992050   -0.396850262992050   -2.777951840944349
+A   -0.396850262992050    0.396850262992050   -2.777951840944349
+A   -0.396850262992050   -0.396850262992050   -1.984251314960249
+A    0.396850262992050    0.396850262992050   -1.984251314960249
+A    0.396850262992050   -0.396850262992050   -1.190550788976150
+A   -0.396850262992050    0.396850262992050   -1.190550788976150
+A   -0.396850262992050   -0.396850262992050   -0.396850262992050
+A    0.396850262992050    0.396850262992050   -0.396850262992050
+A    0.396850262992050   -0.396850262992050    0.396850262992050
+A   -0.396850262992050    0.396850262992050    0.396850262992050
+A   -0.396850262992050   -0.396850262992050    1.190550788976150
+A    0.396850262992050    0.396850262992050    1.190550788976150
+A    0.396850262992050   -0.396850262992050    1.984251314960249
+A   -0.396850262992050    0.396850262992050    1.984251314960249
+A   -0.396850262992050   -0.396850262992050    2.777951840944349
+A    0.396850262992050    0.396850262992050    2.777951840944349
+A    0.396850262992050   -0.396850262992050    3.571652366928448
+A   -0.396850262992050    0.396850262992050    3.571652366928448
+A   -0.396850262992050    1.190550788976150   -3.571652366928449
+A    0.396850262992050    1.984251314960249   -3.571652366928449
+A    0.396850262992050    1.190550788976150   -2.777951840944349
+A   -0.396850262992050    1.984251314960249   -2.777951840944349
+A   -0.396850262992050    1.190550788976150   -1.984251314960249
+A    0.396850262992050    1.984251314960249   -1.984251314960249
+A    0.396850262992050    1.190550788976150   -1.190550788976150
+A   -0.396850262992050    1.984251314960249   -1.190550788976150
+A   -0.396850262992050    1.190550788976150   -0.396850262992050
+A    0.396850262992050    1.984251314960249   -0.396850262992050
+A    0.396850262992050    1.190550788976150    0.396850262992050
+A   -0.396850262992050    1.984251314960249    0.396850262992050
+A   -0.396850262992050    1.190550788976150    1.190550788976150
+A    0.396850262992050    1.984251314960249    1.190550788976150
+A    0.396850262992050    1.190550788976150    1.984251314960249
+A   -0.396850262992050    1.984251314960249    1.984251314960249
+A   -0.396850262992050    1.190550788976150    2.777951840944349
+A    0.396850262992050    1.984251314960249    2.777951840944349
+A    0.396850262992050    1.190550788976150    3.571652366928448
+A   -0.396850262992050    1.984251314960249    3.571652366928448
+A   -0.396850262992050    2.777951840944349   -3.571652366928449
+A    0.396850262992050    3.571652366928448   -3.571652366928449
+A    0.396850262992050    2.777951840944349   -2.777951840944349
+A   -0.396850262992050    3.571652366928448   -2.777951840944349
+A   -0.396850262992050    2.777951840944349   -1.984251314960249
+A    0.396850262992050    3.571652366928448   -1.984251314960249
+A    0.396850262992050    2.777951840944349   -1.190550788976150
+A   -0.396850262992050    3.571652366928448   -1.190550788976150
+A   -0.396850262992050    2.777951840944349   -0.396850262992050
+A    0.396850262992050    3.571652366928448   -0.396850262992050
+A    0.396850262992050    2.777951840944349    0.396850262992050
+A   -0.396850262992050    3.571652366928448    0.396850262992050
+A   -0.396850262992050    2.777951840944349    1.190550788976150
+A    0.396850262992050    3.571652366928448    1.190550788976150
+A    0.396850262992050    2.777951840944349    1.984251314960249
+A   -0.396850262992050    3.571652366928448    1.984251314960249
+A   -0.396850262992050    2.777951840944349    2.777951840944349
+A    0.396850262992050    3.571652366928448    2.777951840944349
+A    0.396850262992050    2.777951840944349    3.571652366928448
+A   -0.396850262992050    3.571652366928448    3.571652366928448
+A    1.190550788976150   -3.571652366928449   -3.571652366928449
+A    1.984251314960249   -2.777951840944349   -3.571652366928449
+A    1.984251314960249   -3.571652366928449   -2.777951840944349
+A    1.190550788976150   -2.777951840944349   -2.777951840944349
+A    1.190550788976150   -3.571652366928449   -1.984251314960249
+A    1.984251314960249   -2.777951840944349   -1.984251314960249
+A    1.984251314960249   -3.571652366928449   -1.190550788976150
+A    1.190550788976150   -2.777951840944349   -1.190550788976150
+A    1.190550788976150   -3.571652366928449   -0.396850262992050
+A    1.984251314960249   -2.777951840944349   -0.396850262992050
+A    1.984251314960249   -3.571652366928449    0.396850262992050
+A    1.190550788976150   -2.777951840944349    0.396850262992050
+A    1.190550788976150   -3.571652366928449    1.190550788976150
+A    1.984251314960249   -2.777951840944349    1.190550788976150
+A    1.984251314960249   -3.571652366928449    1.984251314960249
+A    1.190550788976150   -2.777951840944349    1.984251314960249
+A    1.190550788976150   -3.571652366928449    2.777951840944349
+A    1.984251314960249   -2.777951840944349    2.777951840944349
+A    1.984251314960249   -3.571652366928449    3.571652366928448
+A    1.190550788976150   -2.777951840944349    3.571652366928448
+A    1.190550788976150   -1.984251314960249   -3.571652366928449
+A    1.984251314960249   -1.190550788976150   -3.571652366928449
+A    1.984251314960249   -1.984251314960249   -2.777951840944349
+A    1.190550788976150   -1.190550788976150   -2.777951840944349
+A    1.190550788976150   -1.984251314960249   -1.984251314960249
+A    1.984251314960249   -1.190550788976150   -1.984251314960249
+A    1.984251314960249   -1.984251314960249   -1.190550788976150
+A    1.190550788976150   -1.190550788976150   -1.190550788976150
+A    1.190550788976150   -1.984251314960249   -0.396850262992050
+A    1.984251314960249   -1.190550788976150   -0.396850262992050
+A    1.984251314960249   -1.984251314960249    0.396850262992050
+A    1.190550788976150   -1.190550788976150    0.396850262992050
+A    1.190550788976150   -1.984251314960249    1.190550788976150
+A    1.984251314960249   -1.190550788976150    1.190550788976150
+A    1.984251314960249   -1.984251314960249    1.984251314960249
+A    1.190550788976150   -1.190550788976150    1.984251314960249
+A    1.190550788976150   -1.984251314960249    2.777951840944349
+A    1.984251314960249   -1.190550788976150    2.777951840944349
+A    1.984251314960249   -1.984251314960249    3.571652366928448
+A    1.190550788976150   -1.190550788976150    3.571652366928448
+A    1.190550788976150   -0.396850262992050   -3.571652366928449
+A    1.984251314960249    0.396850262992050   -3.571652366928449
+A    1.984251314960249   -0.396850262992050   -2.777951840944349
+A    1.190550788976150    0.396850262992050   -2.777951840944349
+A    1.190550788976150   -0.396850262992050   -1.984251314960249
+A    1.984251314960249    0.396850262992050   -1.984251314960249
+A    1.984251314960249   -0.396850262992050   -1.190550788976150
+A    1.190550788976150    0.396850262992050   -1.190550788976150
+A    1.190550788976150   -0.396850262992050   -0.396850262992050
+A    1.984251314960249    0.396850262992050   -0.396850262992050
+A    1.984251314960249   -0.396850262992050    0.396850262992050
+A    1.190550788976150    0.396850262992050    0.396850262992050
+A    1.190550788976150   -0.396850262992050    1.190550788976150
+A    1.984251314960249    0.396850262992050    1.190550788976150
+A    1.984251314960249   -0.396850262992050    1.984251314960249
+A    1.190550788976150    0.396850262992050    1.984251314960249
+A    1.190550788976150   -0.396850262992050    2.777951840944349
+A    1.984251314960249    0.396850262992050    2.777951840944349
+A    1.984251314960249   -0.396850262992050    3.571652366928448
+A    1.190550788976150    0.396850262992050    3.571652366928448
+A    1.190550788976150    1.190550788976150   -3.571652366928449
+A    1.984251314960249    1.984251314960249   -3.571652366928449
+A    1.984251314960249    1.190550788976150   -2.777951840944349
+A    1.190550788976150    1.984251314960249   -2.777951840944349
+A    1.190550788976150    1.190550788976150   -1.984251314960249
+A    1.984251314960249    1.984251314960249   -1.984251314960249
+A    1.984251314960249    1.190550788976150   -1.190550788976150
+A    1.190550788976150    1.984251314960249   -1.190550788976150
+A    1.190550788976150    1.190550788976150   -0.396850262992050
+A    1.984251314960249    1.984251314960249   -0.396850262992050
+A    1.984251314960249    1.190550788976150    0.396850262992050
+A    1.190550788976150    1.984251314960249    0.396850262992050
+A    1.190550788976150    1.190550788976150    1.190550788976150
+A    1.984251314960249    1.984251314960249    1.190550788976150
+A    1.984251314960249    1.190550788976150    1.984251314960249
+A    1.190550788976150    1.984251314960249    1.984251314960249
+A    1.190550788976150    1.190550788976150    2.777951840944349
+A    1.984251314960249    1.984251314960249    2.777951840944349
+A    1.984251314960249    1.190550788976150    3.571652366928448
+A    1.190550788976150    1.984251314960249    3.571652366928448
+A    1.190550788976150    2.777951840944349   -3.571652366928449
+A    1.984251314960249    3.571652366928448   -3.571652366928449
+A    1.984251314960249    2.777951840944349   -2.777951840944349
+A    1.190550788976150    3.571652366928448   -2.777951840944349
+A    1.190550788976150    2.777951840944349   -1.984251314960249
+A    1.984251314960249    3.571652366928448   -1.984251314960249
+A    1.984251314960249    2.777951840944349   -1.190550788976150
+A    1.190550788976150    3.571652366928448   -1.190550788976150
+A    1.190550788976150    2.777951840944349   -0.396850262992050
+A    1.984251314960249    3.571652366928448   -0.396850262992050
+A    1.984251314960249    2.777951840944349    0.396850262992050
+A    1.190550788976150    3.571652366928448    0.396850262992050
+A    1.190550788976150    2.777951840944349    1.190550788976150
+A    1.984251314960249    3.571652366928448    1.190550788976150
+A    1.984251314960249    2.777951840944349    1.984251314960249
+A    1.190550788976150    3.571652366928448    1.984251314960249
+A    1.190550788976150    2.777951840944349    2.777951840944349
+A    1.984251314960249    3.571652366928448    2.777951840944349
+A    1.984251314960249    2.777951840944349    3.571652366928448
+A    1.190550788976150    3.571652366928448    3.571652366928448
+A    2.777951840944349   -3.571652366928449   -3.571652366928449
+A    3.571652366928448   -2.777951840944349   -3.571652366928449
+A    3.571652366928448   -3.571652366928449   -2.777951840944349
+A    2.777951840944349   -2.777951840944349   -2.777951840944349
+A    2.777951840944349   -3.571652366928449   -1.984251314960249
+A    3.571652366928448   -2.777951840944349   -1.984251314960249
+A    3.571652366928448   -3.571652366928449   -1.190550788976150
+A    2.777951840944349   -2.777951840944349   -1.190550788976150
+A    2.777951840944349   -3.571652366928449   -0.396850262992050
+A    3.571652366928448   -2.777951840944349   -0.396850262992050
+A    3.571652366928448   -3.571652366928449    0.396850262992050
+A    2.777951840944349   -2.777951840944349    0.396850262992050
+A    2.777951840944349   -3.571652366928449    1.190550788976150
+A    3.571652366928448   -2.777951840944349    1.190550788976150
+A    3.571652366928448   -3.571652366928449    1.984251314960249
+A    2.777951840944349   -2.777951840944349    1.984251314960249
+A    2.777951840944349   -3.571652366928449    2.777951840944349
+A    3.571652366928448   -2.777951840944349    2.777951840944349
+A    3.571652366928448   -3.571652366928449    3.571652366928448
+A    2.777951840944349   -2.777951840944349    3.571652366928448
+A    2.777951840944349   -1.984251314960249   -3.571652366928449
+A    3.571652366928448   -1.190550788976150   -3.571652366928449
+A    3.571652366928448   -1.984251314960249   -2.777951840944349
+A    2.777951840944349   -1.190550788976150   -2.777951840944349
+A    2.777951840944349   -1.984251314960249   -1.984251314960249
+A    3.571652366928448   -1.190550788976150   -1.984251314960249
+A    3.571652366928448   -1.984251314960249   -1.190550788976150
+A    2.777951840944349   -1.190550788976150   -1.190550788976150
+A    2.777951840944349   -1.984251314960249   -0.396850262992050
+A    3.571652366928448   -1.190550788976150   -0.396850262992050
+A    3.571652366928448   -1.984251314960249    0.396850262992050
+A    2.777951840944349   -1.190550788976150    0.396850262992050
+A    2.777951840944349   -1.984251314960249    1.190550788976150
+A    3.571652366928448   -1.190550788976150    1.190550788976150
+A    3.571652366928448   -1.984251314960249    1.984251314960249
+A    2.777951840944349   -1.190550788976150    1.984251314960249
+A    2.777951840944349   -1.984251314960249    2.777951840944349
+A    3.571652366928448   -1.190550788976150    2.777951840944349
+A    3.571652366928448   -1.984251314960249    3.571652366928448
+A    2.777951840944349   -1.190550788976150    3.571652366928448
+A    2.777951840944349   -0.396850262992050   -3.571652366928449
+A    3.571652366928448    0.396850262992050   -3.571652366928449
+A    3.571652366928448   -0.396850262992050   -2.777951840944349
+A    2.777951840944349    0.396850262992050   -2.777951840944349
+A    2.777951840944349   -0.396850262992050   -1.984251314960249
+A    3.571652366928448    0.396850262992050   -1.984251314960249
+A    3.571652366928448   -0.396850262992050   -1.190550788976150
+A    2.777951840944349    0.396850262992050   -1.190550788976150
+A    2.777951840944349   -0.396850262992050   -0.396850262992050
+A    3.571652366928448    0.396850262992050   -0.396850262992050
+A    3.571652366928448   -0.396850262992050    0.396850262992050
+A    2.777951840944349    0.396850262992050    0.396850262992050
+A    2.777951840944349   -0.396850262992050    1.190550788976150
+A    3.571652366928448    0.396850262992050    1.190550788976150
+A    3.571652366928448   -0.396850262992050    1.984251314960249
+A    2.777951840944349    0.396850262992050    1.984251314960249
+A    2.777951840944349   -0.396850262992050    2.777951840944349
+A    3.571652366928448    0.396850262992050    2.777951840944349
+A    3.571652366928448   -0.396850262992050    3.571652366928448
+A    2.777951840944349    0.396850262992050    3.571652366928448
+A    2.777951840944349    1.190550788976150   -3.571652366928449
+A    3.571652366928448    1.984251314960249   -3.571652366928449
+A    3.571652366928448    1.190550788976150   -2.777951840944349
+A    2.777951840944349    1.984251314960249   -2.777951840944349
+A    2.777951840944349    1.190550788976150   -1.984251314960249
+A    3.571652366928448    1.984251314960249   -1.984251314960249
+A    3.571652366928448    1.190550788976150   -1.190550788976150
+A    2.777951840944349    1.984251314960249   -1.190550788976150
+A    2.777951840944349    1.190550788976150   -0.396850262992050
+A    3.571652366928448    1.984251314960249   -0.396850262992050
+A    3.571652366928448    1.190550788976150    0.396850262992050
+A    2.777951840944349    1.984251314960249    0.396850262992050
+A    2.777951840944349    1.190550788976150    1.190550788976150
+A    3.571652366928448    1.984251314960249    1.190550788976150
+A    3.571652366928448    1.190550788976150    1.984251314960249
+A    2.777951840944349    1.984251314960249    1.984251314960249
+A    2.777951840944349    1.190550788976150    2.777951840944349
+A    3.571652366928448    1.984251314960249    2.777951840944349
+A    3.571652366928448    1.190550788976150    3.571652366928448
+A    2.777951840944349    1.984251314960249    3.571652366928448
+A    2.777951840944349    2.777951840944349   -3.571652366928449
+A    3.571652366928448    3.571652366928448   -3.571652366928449
+A    3.571652366928448    2.777951840944349   -2.777951840944349
+A    2.777951840944349    3.571652366928448   -2.777951840944349
+A    2.777951840944349    2.777951840944349   -1.984251314960249
+A    3.571652366928448    3.571652366928448   -1.984251314960249
+A    3.571652366928448    2.777951840944349   -1.190550788976150
+A    2.777951840944349    3.571652366928448   -1.190550788976150
+A    2.777951840944349    2.777951840944349   -0.396850262992050
+A    3.571652366928448    3.571652366928448   -0.396850262992050
+A    3.571652366928448    2.777951840944349    0.396850262992050
+A    2.777951840944349    3.571652366928448    0.396850262992050
+A    2.777951840944349    2.777951840944349    1.190550788976150
+A    3.571652366928448    3.571652366928448    1.190550788976150
+A    3.571652366928448    2.777951840944349    1.984251314960249
+A    2.777951840944349    3.571652366928448    1.984251314960249
+A    2.777951840944349    2.777951840944349    2.777951840944349
+A    3.571652366928448    3.571652366928448    2.777951840944349
+A    3.571652366928448    2.777951840944349    3.571652366928448
+A    2.777951840944349    3.571652366928448    3.571652366928448
diff --git a/Examples/NVT/HMA/inp500_1.0 b/Examples/NVT/HMA/inp500_1.0
new file mode 100644
index 00000000..60396c08
--- /dev/null
+++ b/Examples/NVT/HMA/inp500_1.0
@@ -0,0 +1,79 @@
+# Run_Name
+out500
+!------------------------------------------------------------------------------
+
+# Sim_Type
+nvt
+!------------------------------------------------------------------------------
+
+# Nbr_Species
+1
+!------------------------------------------------------------------------------
+
+# Molecule_Files
+lj.mcf 500
+!------------------------------------------------------------------------------
+
+# Box_Info
+1
+cubic
+7.937005259840997
+!------------------------------------------------------------------------------
+
+# VDW_Style
+lj cut 3.
+!---------------------------------------------------------one line for each box
+
+# Temperature_Info
+1.
+!---------------------------------------------------one entry/line for each box
+
+# Move_Probability_Info
+
+# Prob_Translation
+100.
+0.1
+
+# Done_Probability_Info
+!---------------------one line for each box and one entry/line for each species
+
+# Start_Type
+read_config 500 fcc500.xyz
+!---------------------one line for each species and one entry/line for each box
+
+# Run_Type
+production   50000 100
+!------------------------------------------------------------------------------
+
+# Simulation_Length_Info
+units        sweeps
+prop_freq    10
+coord_freq   200
+run          2000
+!------------------------------------------------------------------------------
+
+# Property_Info 1 
+energy_total
+energy_HMA
+pressure
+pressure_HMA 1311.62
+!------------------------------------------------------------------------------
+
+# CBMC_Info
+kappa_ins 12
+rcut_cbmc 6.5
+!------------------------------------------------------------------------------
+
+# Rcutoff_Low
+0.8
+!------------------------------------------------------------------------------
+
+# Pair_Energy
+true
+!------------------------------------------------------------------------------
+
+# Seed_Info
+9799807 23445206
+!------------------------------------------------------------------------------
+
+END
diff --git a/Examples/NVT/HMA/lj.mcf b/Examples/NVT/HMA/lj.mcf
new file mode 100644
index 00000000..b004374f
--- /dev/null
+++ b/Examples/NVT/HMA/lj.mcf
@@ -0,0 +1,34 @@
+!*******************************************************************************
+!Molecular connectivity file for n2.pdb
+!*******************************************************************************
+!Atom Format
+!index type element mass charge vdw_type parameters
+!vdw_type="LJ", parms=epsilon sigma
+!vdw_type="Mie", parms=epsilon sigma repulsion_exponent dispersion_exponent
+
+# Atom_Info
+1
+1     A       A     1.000   0.000  LJ     1.000     1.000
+
+!Bond Format
+!index i j type parameters
+!type="fixed", parms=bondLength
+
+# Bond_Info
+0
+
+!Angle Format
+!index i j k type parameters
+!type="fixed", parms=equilibrium_angle
+!type="harmonic", parms=force_constant equilibrium_angle
+
+# Angle_Info
+0
+
+# Dihedral_Info
+0
+
+# Improper_Info
+0
+
+END

From f68ef7039f6d915656fa683a97d0602f4bad69f7 Mon Sep 17 00:00:00 2001
From: "Andrew J. Schultz" <ajschult784@gmail.com>
Date: Fri, 21 May 2021 14:05:03 -0400
Subject: [PATCH 10/11] Add bits to describe additional properties in docs

---
 Documentation/source/guides/input_files.rst | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/Documentation/source/guides/input_files.rst b/Documentation/source/guides/input_files.rst
index 15aa30cf..f9c123b4 100644
--- a/Documentation/source/guides/input_files.rst
+++ b/Documentation/source/guides/input_files.rst
@@ -1353,6 +1353,8 @@ Additional supported keywords are:
 * ``nmols``: Number of molecules of each species
 * ``density``: Density of each species in Å\ :sup:`3`
 * ``mass_density``: Density of the system in kg/m\ :sup:`3`
+* ``energy_HMA``: Energy of the system measured with HMA
+* ``pressure_HMA``: Pressure of the system measured with HMA, in bar.  The harmonic pressure must be provided as a second argument on this line.
 
 For example, if you would like total energy, volume and pressure of a one box
 system to be written, you may specify the following:

From cdc2fc125b2eea5b67e8b083e3a225e527bb62cd Mon Sep 17 00:00:00 2001
From: "Andrew J. Schultz" <ajschult784@gmail.com>
Date: Fri, 21 May 2021 15:34:57 -0400
Subject: [PATCH 11/11] Add new page in guides for HMA

---
 Documentation/source/guides/hma.rst         | 42 +++++++++++++++++++++
 Documentation/source/guides/input_files.rst |  4 +-
 Documentation/source/index.rst              |  1 +
 3 files changed, 45 insertions(+), 2 deletions(-)
 create mode 100644 Documentation/source/guides/hma.rst

diff --git a/Documentation/source/guides/hma.rst b/Documentation/source/guides/hma.rst
new file mode 100644
index 00000000..7c20e40d
--- /dev/null
+++ b/Documentation/source/guides/hma.rst
@@ -0,0 +1,42 @@
+.. highlight:: none
+
+.. _ch:hma:
+
+Harmonically Mapped Averaging (HMA)
+===================================
+
+.. _sec:hma:
+
+Background
+----------
+
+Harmonically mapped averaging (HMA) is a technique that improves the precision
+of measured themodynamic properites of crystalline systems.  Beyond precision,
+other benefits include smaller potential-truncation effects, finite-size
+effects, faster equilibration and shorter decorrelation time.
+
+Limitations
+-----------
+
+HMA can only be used in the NVT ensemble and when the atoms only vibrate around
+their lattice sites (without diffusing).  HMA can be used in Cassandra with any
+pair potential, but not with Ewald summation.
+
+Truncation
+----------
+
+In order to supress property fluctuations due to atoms moving into and out of
+the truncation radius, the potential in Cassandra is truncated based on the
+lattice site distance instead of the distance between the atoms.  This will lead
+to a small difference in the measured properties, but the difference will
+disappear in the limit of long truncation.
+
+Example
+-------
+
+An simuation that uses HMA is included in `Examples/NVT/HMA/ <https://github.com/MaginnGroup/Cassandra/tree/master/Examples/NVT/HMA>`_ in the source tree.  The example simulation runs a Lennard-Jones FCC crystal near its melting point.  Energy and pressure are computed using both conventional and HMA and reported in the property output.
+
+References
+----------
+
+Sabry G. Moustafa, Andrew J. Schultz, and David A. Kofke, Very fast averaging of thermal properties of crystals by molecular simulation, `Phys. Rev. E [92], 043303 (2015) <https://dx.doi.org/10.1103/PhysRevE.92.043303>`_
diff --git a/Documentation/source/guides/input_files.rst b/Documentation/source/guides/input_files.rst
index f9c123b4..49b40db0 100644
--- a/Documentation/source/guides/input_files.rst
+++ b/Documentation/source/guides/input_files.rst
@@ -1353,8 +1353,8 @@ Additional supported keywords are:
 * ``nmols``: Number of molecules of each species
 * ``density``: Density of each species in Å\ :sup:`3`
 * ``mass_density``: Density of the system in kg/m\ :sup:`3`
-* ``energy_HMA``: Energy of the system measured with HMA
-* ``pressure_HMA``: Pressure of the system measured with HMA, in bar.  The harmonic pressure must be provided as a second argument on this line.
+* ``energy_HMA``: Energy of the system measured with :ref:`HMA <ch:hma>`
+* ``pressure_HMA``: Pressure of the system measured with :ref:`HMA <ch:hma>`, in bar.  The harmonic pressure must be provided as a second argument on this line.
 
 For example, if you would like total energy, volume and pressure of a one box
 system to be written, you may specify the following:
diff --git a/Documentation/source/index.rst b/Documentation/source/index.rst
index 7e2cdade..7fd6f447 100644
--- a/Documentation/source/index.rst
+++ b/Documentation/source/index.rst
@@ -24,6 +24,7 @@ Cassandra Monte Carlo Software
    guides/input_files
    guides/utilities
    guides/solids
+   guides/hma
 
 .. toctree::
    :hidden: