From a70a2cdd5d086e7c23f301ca43a82b94a7dfe02a Mon Sep 17 00:00:00 2001
From: Daniel Hollas <daniel.hollas@bristol.ac.uk>
Date: Thu, 25 Jan 2024 14:23:34 +0000
Subject: [PATCH] Bring back TeraChem file interface (#165)

The Amber MPI interface in TC was removed in 2021.
We've used this interface for ground state MD.
Instead, we now have the TCPB interface, but still
in an experimental stage. Therefore, for now we bring back
the simple file interface as a quick hassle free way
to run MD with TeraChem.
WARNING: For small / medium sized molecules, you will
lose a lot of time during repeated GPU initialization!
---
 interfaces/TERACHEM/r.terachem | 86 ++++++++++++++++++++++++++++++++++
 src/force_abin.F90             |  2 +-
 src/init.F90                   |  2 +-
 utils/abin-randomint.f90       |  1 +
 4 files changed, 89 insertions(+), 2 deletions(-)
 create mode 100755 interfaces/TERACHEM/r.terachem

diff --git a/interfaces/TERACHEM/r.terachem b/interfaces/TERACHEM/r.terachem
new file mode 100755
index 00000000..f74577c5
--- /dev/null
+++ b/interfaces/TERACHEM/r.terachem
@@ -0,0 +1,86 @@
+#!/bin/bash
+# This script is called by ABIN as:
+# TERACHEM/r.terachem <timestep> <ibead>
+
+# SETUP the TeraChem environment
+# $TERAEXE (used below) should hold the path to the terachem binary
+# The SetEnvironment.sh script that handles this is specific to PHOTOX clusters.
+source ./SetEnvironment.sh TERACHEM
+
+set -u
+
+cd "$(dirname "$0")" || exit 1
+ibead=$2
+input="input$ibead.com"
+geom="../geom.dat.$ibead"
+natom=$(wc -l < "$geom")
+
+### USER INPUT FOR  TERACHEM
+cat > "$input" << EOF
+basis 		6-31g*
+charge   	0
+spinmult	1
+method 		pbe
+convthre        1e-6
+threall         1e-13
+EOF
+
+parallel=0 # =1 if we execute ABIN in parallel for PIMD
+numgpus=1  # number of gpus per single point calculation
+
+### END OF USER INPUT
+
+scrdir="./scratch$ibead"
+
+if [[ $parallel -eq "1" ]];then
+    gpuid=""
+    (( gpu0=ibead-1 ))
+    (( gpu0=gpu0*numgpus ))
+else
+    gpu0=0
+fi
+
+for ((i=0;i<numgpus;i++ )) {
+    (( gpuindex=gpu0+i ))
+    gpuid="${gpuid-} $gpuindex"
+}
+
+### Pass the wave function from previous step if it exists
+guess_string=""
+if [[ -f $scrdir/c0 ]];then
+    guess_string="guess $scrdir/c0"
+fi
+if [[ -f $scrdir/ca0 && -f $scrdir/cb0 ]];then
+    guess_string="guess $scrdir/ca0 $scrdir/cb0"
+fi
+
+cat >> "$input" << EOF
+scrdir        $scrdir
+keep_scr      yes
+coordinates   input$ibead.xyz
+units         angstrom
+gpus          $numgpus $gpuid
+run           gradient
+$guess_string
+end
+EOF
+
+### Create XYZ geometry file
+echo -ne "$natom\n\n" > "input$ibead.xyz"
+head -n "$natom" "$geom" >> "input$ibead.xyz"
+
+### Launch TeraChem
+export OMP_NUM_THREADS=$numgpus
+
+if $TERAEXE "$input" > "$input.out" 2>&1; then
+    cp "$input.out" "$input.out.old"
+else
+    echo "ERROR: TeraChem calculation probably failed."
+    echo "See TERACHEM/$input.out.error"
+    cp "$input.out" "$input.out.error"
+    exit 2
+fi
+
+### Extract energy and forces
+grep 'FINAL ENERGY' "$input.out" | awk '{print $3}' > "../engrad.dat.$ibead"
+grep -A"$natom" 'dE/dX' "$input.out" | tail -"$natom" >> "../engrad.dat.$ibead" 
diff --git a/src/force_abin.F90 b/src/force_abin.F90
index 26d03f0c..e6d49bca 100644
--- a/src/force_abin.F90
+++ b/src/force_abin.F90
@@ -122,7 +122,7 @@ subroutine call_shell(shellscript, it, iw, ipimd, abort)
       ! Passing arguments to bash script
       ! First argument is time step
       ! Second argument is the bead index, neccessary for parallel calculations
-      write (call_cmd, '(A,I0,I4.3)') './'//trim(shellscript)//' ', it, iw
+      write (call_cmd, '(A,I0,I4.3)') trim(shellscript)//' ', it, iw
       call_cmd = append_rank(call_cmd)
 
       ! For SH, pass the 4th parameter: precision of forces as 10^(-force_accu1)
diff --git a/src/init.F90 b/src/init.F90
index 14647dd9..f1543bb3 100644
--- a/src/init.F90
+++ b/src/init.F90
@@ -116,7 +116,7 @@ subroutine init(dt)
       ! - general: Most basic MD settings + misc
       ! - nhcopt:  parameters for thermostats
       ! - remd:    parameters for Replica Exchange MD
-      ! - sh:      parameters for Surface Hopping
+      ! - sh:      parameters for Surface Hopping, moved to mod_sh module
       ! - system:  system-specific parameters for model potentials, masses, SHAKE constraints...
       ! - lz:      parameters for Landau-Zener excited state dynamics.
       ! - qmmm:    parameters for internal QMMM (not really tested).
diff --git a/utils/abin-randomint.f90 b/utils/abin-randomint.f90
index 0ad7a083..8b599960 100644
--- a/utils/abin-randomint.f90
+++ b/utils/abin-randomint.f90
@@ -39,6 +39,7 @@ program abin_randomint
    do i = 1, nran
       write (*, '(I0)') irans(i)
    end do
+   deallocate (irans)
 end program
 
 subroutine get_cmdline(iseed, nran)