From 984c71026fe08df083c3ab3098d36f227a3eacbd Mon Sep 17 00:00:00 2001
From: Daniel Hollas <daniel.hollas@bristol.ac.uk>
Date: Fri, 17 May 2024 17:26:34 +0100
Subject: [PATCH 01/10] WIP: DFTB+ interface upgrade

---
 interfaces/DFTB/dftb_in.hsd |  7 +++--
 interfaces/DFTB/r.dftb      | 55 +++++++++++++++++++++++++++----------
 2 files changed, 44 insertions(+), 18 deletions(-)

diff --git a/interfaces/DFTB/dftb_in.hsd b/interfaces/DFTB/dftb_in.hsd
index 482fa5b0..a2e7e59c 100644
--- a/interfaces/DFTB/dftb_in.hsd
+++ b/interfaces/DFTB/dftb_in.hsd
@@ -4,10 +4,11 @@ Geometry = GenFormat {
 
 Hamiltonian = DFTB {
   SCC = Yes
+  MaxSCCIterations = 100
   charge = 0.0
   #ReadInitialCharges = yes
   SlaterKosterFiles = Type2FileNames {
-   Prefix = "/home/hollas/3ob-2-1/"
+   Prefix = "/home/hollas/3ob-3-1/"
    Separator = "-"
    Suffix = ".skf"
   }
@@ -17,8 +18,8 @@ Hamiltonian = DFTB {
   }
 }
 
-Options {
-CalculateForces = Yes
+Analysis {
+  CalculateForces = Yes
 }
 
 ParserOptions {
diff --git a/interfaces/DFTB/r.dftb b/interfaces/DFTB/r.dftb
index 9aa97cfc..1e5a6fa0 100755
--- a/interfaces/DFTB/r.dftb
+++ b/interfaces/DFTB/r.dftb
@@ -1,6 +1,15 @@
 #!/bin/bash
+# File interface to DFTB+ program, tested for version 24.1
 cd $(dirname $0)
-source ../SetEnvironment.sh DFTB
+set -u
+
+if [[ -f ../SetEnvironment.sh ]];then
+  # This is specific to Prague clusters
+  source ../SetEnvironment.sh DFTB
+else
+  # We assume dftb+ is in PATH already
+  DFTBEXE=dftb+
+fi
 
 timestep=$1
 ibead=$2
@@ -8,13 +17,15 @@ input=input$ibead
 geom=../geom.dat.$ibead
 
 natom=$(wc -l < $geom)
-WRKDIR=OUT$ibead.$natom
-rm -rf ${WRKDIR}.old
-if [[ -d $WRKDIR ]];then
-  mv $WRKDIR ${WRKDIR}.old
+CHARGES_FILE="charges.bin"
+
+# Working directory for the DFTB calculation
+WORKDIR=OUT$ibead.$natom
+rm -rf ${WORKDIR}.previous
+if [[ -d $WORKDIR ]];then
+    mv $WORKDIR ${WORKDIR}.previous
 fi
-mkdir -p $WRKDIR
-cd $WRKDIR
+mkdir -p $WORKDIR
 
 # You have to specify, which elements are present in your system
 # i.e. define array id[x]="element"
@@ -39,16 +50,17 @@ awk -v natom="$natom" 'BEGIN{
          print NR, i, $2, $3, $4
       }
    }
-}'  ../$geom > geom_in.gen
+}' $geom > $WORKDIR/geom_in.gen
 
 # Read initial charge distribution
-if [ -e charges.bin ];then
-   sed 's/#ReadInitialCharges/ReadInitialCharges/' ../dftb_in.hsd > dftb_in.hsd
+if [[ -f $CHARGES_FILE ]]; then
+   cp $CHARGES_FILE $WORKDIR/
+   sed 's/#ReadInitialCharges/ReadInitialCharges/' dftb_in.hsd > $WORKDIR/dftb_in.hsd
 else
-   cp ../dftb_in.hsd .
+   cp dftb_in.hsd $WORKDIR/
 fi
 
-rm -f detailed.out
+cd $WORKDIR || exit 2
 
 $DFTBEXE &> $input.out
 if [[ $? -eq 0 ]];then
@@ -60,7 +72,20 @@ else
    exit 2
 fi
 
-# Extract energy and gradients
+# Extract energy
 grep 'Total energy:' detailed.out | awk '{print $3}' > ../../engrad.dat.$ibead
-awk -v natom=$natom '{if ($2=="Forces"){for (i=1;i<=natom;i++){getline;printf"%3.15e %3.15e %3.15e \n",-$1,-$2,-$3}}}' \
- detailed.out >> ../../engrad.dat.$ibead
+
+# Extract gradients (note the conversion from forces to gradients)
+awk -v natom=$natom '{
+  if ($2 == "Forces") {
+    for (i=0; i < natom; i++) {
+      getline;
+      printf("%3.15e %3.15e %3.15e\n", -$2, -$3, -$4)
+    }
+  }
+}' detailed.out >> ../../engrad.dat.$ibead
+
+# Copy the charges for next step
+if [[ -f $CHARGES_FILE ]];then
+  cp $CHARGES_FILE ../
+fi

From 3037dffc1cf3db3280b272899aeaa8819a67deec Mon Sep 17 00:00:00 2001
From: Daniel Hollas <daniel.hollas@bristol.ac.uk>
Date: Sat, 18 May 2024 22:31:06 +0100
Subject: [PATCH 02/10] More cleanup

---
 interfaces/DFTB/r.dftb | 89 +++++++++++++++++++++---------------------
 1 file changed, 44 insertions(+), 45 deletions(-)

diff --git a/interfaces/DFTB/r.dftb b/interfaces/DFTB/r.dftb
index 1e5a6fa0..23c4455f 100755
--- a/interfaces/DFTB/r.dftb
+++ b/interfaces/DFTB/r.dftb
@@ -1,9 +1,11 @@
 #!/bin/bash
-# File interface to DFTB+ program, tested for version 24.1
-cd $(dirname $0)
+# File interface to DFTB+ program, tested with version 24.1
+# NOTE: You need to tweak the awk routine for converting XYZ geometries
+# from ABIN to dftb+ compatible input file.
+cd "$(dirname "$0")" || exit 2
 set -u
 
-if [[ -f ../SetEnvironment.sh ]];then
+if [[ -f ../SetEnvironment.sh ]]; then
   # This is specific to Prague clusters
   source ../SetEnvironment.sh DFTB
 else
@@ -11,81 +13,78 @@ else
   DFTBEXE=dftb+
 fi
 
-timestep=$1
-ibead=$2
-input=input$ibead
-geom=../geom.dat.$ibead
+timestep="$1"
+ibead="$2"
+geom=../geom.dat."$ibead"
 
-natom=$(wc -l < $geom)
+natom=$(wc -l < "$geom")
 CHARGES_FILE="charges.bin"
 
 # Working directory for the DFTB calculation
-WORKDIR=OUT$ibead.$natom
-rm -rf ${WORKDIR}.previous
-if [[ -d $WORKDIR ]];then
-    mv $WORKDIR ${WORKDIR}.previous
+WORKDIR=CALC."$ibead"
+rm -rf "${WORKDIR}".previous
+if [[ -d "$WORKDIR" ]];then
+    mv "$WORKDIR" "${WORKDIR}".previous
 fi
-mkdir -p $WORKDIR
+mkdir -p "$WORKDIR"
 
-# You have to specify, which elements are present in your system
-# i.e. define array id[x]="element"
+# You have to specify, which (unique) elements are present in your system
+# i.e. define array el[i]="element"
 # No extra elements are allowed.
 awk -v natom="$natom" 'BEGIN{
-   id[1]="O"
-   id[2]="H"
-   nid=2    # number of different elements
+   el[1]="O"
+   el[2]="H"
+   num_el=2    # number of unique elements
 
-# END OF USER INPUT
-   print natom,"C"
-   for (i=1;i<=nid;i++) {
-      printf"%s ",id[i]
+   # END OF USER INPUT
+   print natom, "C"
+   for (i = 0; i < num_el; i++) {
+      printf"%s ", el[i]
    }
    print ""
    print ""
 }
 #conversion of xyz input to dftb geom
 {
-   for (i=1;i<=nid;i++) {
-      if ( $1 == id[i] ) {
+   for (i = 0; i < num_el; i++) {
+      if ( $1 == el[i] ) {
          print NR, i, $2, $3, $4
       }
    }
-}' $geom > $WORKDIR/geom_in.gen
+}' "$geom" > "$WORKDIR"/geom_in.gen
 
 # Read initial charge distribution
-if [[ -f $CHARGES_FILE ]]; then
-   cp $CHARGES_FILE $WORKDIR/
-   sed 's/#ReadInitialCharges/ReadInitialCharges/' dftb_in.hsd > $WORKDIR/dftb_in.hsd
+if [[ -f "$CHARGES_FILE" ]]; then
+   cp "$CHARGES_FILE" "$WORKDIR"/
+   sed 's/#ReadInitialCharges/ReadInitialCharges/' dftb_in.hsd > "$WORKDIR"/dftb_in.hsd
 else
-   cp dftb_in.hsd $WORKDIR/
+   cp dftb_in.hsd "$WORKDIR"/
 fi
 
-cd $WORKDIR || exit 2
+cd "$WORKDIR" || exit 2
 
-$DFTBEXE &> $input.out
-if [[ $? -eq 0 ]];then
-   cp $input.out $input.out.old
-else
-   echo "ERROR: DFTB calculation probably failed."
-   echo "See $input.out.error" 
-   cp $input.out $input.out.error
+# Everything is ready, run dftb+!
+$DFTBEXE &> dftb.out
+if [[ $? -ne 0 ]]; then
+   echo "ERROR: DFTB calculation failed."
+   echo "See $PWD/dftb.out" 
    exit 2
 fi
 
-# Extract energy
-grep 'Total energy:' detailed.out | awk '{print $3}' > ../../engrad.dat.$ibead
+# Extract energy from dftb+ output
+grep 'Total energy:' detailed.out | awk '{print $3}' > ../../engrad.dat."$ibead"
 
 # Extract gradients (note the conversion from forces to gradients)
-awk -v natom=$natom '{
+awk -v natom="$natom" '{
   if ($2 == "Forces") {
-    for (i=0; i < natom; i++) {
+    for (i = 0; i < natom; i++) {
       getline;
       printf("%3.15e %3.15e %3.15e\n", -$2, -$3, -$4)
     }
   }
-}' detailed.out >> ../../engrad.dat.$ibead
+}' detailed.out >> ../../engrad.dat."$ibead"
 
-# Copy the charges for next step
-if [[ -f $CHARGES_FILE ]];then
-  cp $CHARGES_FILE ../
+# Copy the charges for the next step
+if [[ -f "$CHARGES_FILE" ]]; then
+  cp "$CHARGES_FILE" ../
 fi

From d1408f7e50026759ce3b927129c493ce07224dfe Mon Sep 17 00:00:00 2001
From: Daniel Hollas <daniel.hollas@bristol.ac.uk>
Date: Sat, 18 May 2024 22:31:34 +0100
Subject: [PATCH 03/10] Tweak CONTRIBUTING.md

---
 CONTRIBUTING.md | 25 ++++++++++---------------
 1 file changed, 10 insertions(+), 15 deletions(-)

diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md
index 53eb2b5e..af9f2018 100644
--- a/CONTRIBUTING.md
+++ b/CONTRIBUTING.md
@@ -1,6 +1,6 @@
 # How to contribute to ABIN
 
-Thanks for contributing to our code! 💜
+Thanks for contributing to our code!💜
 Here's a couple of guidelines that you should keep in mind.
 
 ## Setup your environment
@@ -32,15 +32,15 @@ Here's a quick summary of our code style that we try to adhere to:
  - variables and subroutines in modules should be private by default, use `private` attribute,e.g.
 ```fortran
 module mod_my_module
-  private
-  integer :: public_var
-  integer :: private_var
-  public :: public_var
-  public :: public_subroutine
+   private
+   integer :: public_var
+   integer :: private_var
+   public :: public_var
+   public :: public_subroutine
 contains
- subroutine public_subroutine
-    ...
- end subroutine public_subroutine
+   subroutine public_subroutine
+      ...
+   end subroutine public_subroutine
 end module mod_my_module
 ```
 
@@ -68,7 +68,7 @@ Here's a quick summary of our formatting style, as it is defined in `.fprettify.
 
 ### Inspecting Git history
 
-To ignore bulk whitespace changes in blame history, use:
+To ignore bulk whitespace changes in git blame history, use:
 ```sh
 git blame --ignore-revs-file .git-blame-ignore-revs
 ```
@@ -78,11 +78,6 @@ or to do it automatically:
 git config blame.ignoreRevsFile .git-blame-ignore-revs
 ```
 
-Unfortunately, this is not yet supported in
-[the GitHub UI](https://github.community/t/support-ignore-revs-file-in-githubs-blame-view/3256),
-but Github UI already allows to browse git blame a bit.
-
-
 ## Submitting code changes
 
 Last but not least, to get your code merged to the main repository, please open a Pull Request (PR) on Github.

From eb1bf4e14dad7f2a2729864e3550c41cbe2d3680 Mon Sep 17 00:00:00 2001
From: Daniel Hollas <daniel.hollas@bristol.ac.uk>
Date: Tue, 15 Oct 2024 13:17:25 +0100
Subject: [PATCH 04/10] Use XYZ geometry as input

---
 interfaces/DFTB/dftb_in.hsd |   6 +-
 interfaces/DFTB/r.dftb      | 112 ++++++++++++++++--------------------
 2 files changed, 52 insertions(+), 66 deletions(-)

diff --git a/interfaces/DFTB/dftb_in.hsd b/interfaces/DFTB/dftb_in.hsd
index a2e7e59c..6beda72b 100644
--- a/interfaces/DFTB/dftb_in.hsd
+++ b/interfaces/DFTB/dftb_in.hsd
@@ -1,5 +1,5 @@
-Geometry = GenFormat {
- <<< "geom_in.gen"
+Geometry = xyzFormat {
+ <<< "geometry.xyz"
 }
 
 Hamiltonian = DFTB {
@@ -8,7 +8,7 @@ Hamiltonian = DFTB {
   charge = 0.0
   #ReadInitialCharges = yes
   SlaterKosterFiles = Type2FileNames {
-   Prefix = "/home/hollas/3ob-3-1/"
+   Prefix = "/path/to/parameters/3ob-3-1/"
    Separator = "-"
    Suffix = ".skf"
   }
diff --git a/interfaces/DFTB/r.dftb b/interfaces/DFTB/r.dftb
index 23c4455f..9fee2fda 100755
--- a/interfaces/DFTB/r.dftb
+++ b/interfaces/DFTB/r.dftb
@@ -1,7 +1,12 @@
 #!/bin/bash
-# File interface to DFTB+ program, tested with version 24.1
-# NOTE: You need to tweak the awk routine for converting XYZ geometries
-# from ABIN to dftb+ compatible input file.
+# File interface to DFTB+ program, https://dftbplus.org/
+# This file typically does not need to be modified.
+# The input DFTB parameters are given in a file `dftb_in.hsd` in this directory,
+# which you need to modify for your needs.
+#
+# Tested with version 24.1. Other versions might work, but always test
+# by verifying energy conservation in a short NVE simulation.
+# Versions older than 20.1 will not work since they do not support XYZ geometry input.
 cd "$(dirname "$0")" || exit 2
 set -u
 
@@ -9,82 +14,63 @@ if [[ -f ../SetEnvironment.sh ]]; then
   # This is specific to Prague clusters
   source ../SetEnvironment.sh DFTB
 else
-  # We assume dftb+ is in PATH already
+  # We assume dftb+ is in PATH already. If not, add it here.
   DFTBEXE=dftb+
 fi
 
 timestep="$1"
 ibead="$2"
 geom=../geom.dat."$ibead"
-
 natom=$(wc -l < "$geom")
-CHARGES_FILE="charges.bin"
 
-# Working directory for the DFTB calculation
-WORKDIR=CALC."$ibead"
-rm -rf "${WORKDIR}".previous
-if [[ -d "$WORKDIR" ]];then
-    mv "$WORKDIR" "${WORKDIR}".previous
-fi
-mkdir -p "$WORKDIR"
+function prepare_dftb_inputs() {
+    # Working directory for the DFTB calculation
+    WORKDIR="CALC.$ibead"
+    rm -rf "${WORKDIR}".previous
+    if [[ -d "$WORKDIR" ]];then
+        mv "$WORKDIR" "${WORKDIR}".previous
+    fi
+    mkdir -p "$WORKDIR"
 
-# You have to specify, which (unique) elements are present in your system
-# i.e. define array el[i]="element"
-# No extra elements are allowed.
-awk -v natom="$natom" 'BEGIN{
-   el[1]="O"
-   el[2]="H"
-   num_el=2    # number of unique elements
+    echo -e "$atom\n" > "$WORKDIR/geometry.xyz"
+    cat $geom >> "$WORKDIR/geometry.xyz"
 
-   # END OF USER INPUT
-   print natom, "C"
-   for (i = 0; i < num_el; i++) {
-      printf"%s ", el[i]
-   }
-   print ""
-   print ""
+    cp dftb_in.hsd "$WORKDIR"
+    # Read charge distribution from previous step if available
+    charges_file="$WORKDIR.previous/charges.bin"
+    if [[ -f "$charges_file" ]]; then
+       cp "$charges_file" "$WORKDIR"
+       sed -i 's/#ReadInitialCharges/ReadInitialCharges/' "$WORKDIR/dftb_in.hsd"
+    fi
 }
-#conversion of xyz input to dftb geom
-{
-   for (i = 0; i < num_el; i++) {
-      if ( $1 == el[i] ) {
-         print NR, i, $2, $3, $4
-      }
-   }
-}' "$geom" > "$WORKDIR"/geom_in.gen
 
-# Read initial charge distribution
-if [[ -f "$CHARGES_FILE" ]]; then
-   cp "$CHARGES_FILE" "$WORKDIR"/
-   sed 's/#ReadInitialCharges/ReadInitialCharges/' dftb_in.hsd > "$WORKDIR"/dftb_in.hsd
-else
-   cp dftb_in.hsd "$WORKDIR"/
-fi
+function extract_energy_and_gradients() {
+    # Extract energy dftb+ output file
+    dftb_out="detailed.out"
+    engrad_file=$1
+    grep 'Total energy:' $dftb_out | awk '{print $3}' > $engrad_file
+
+    # Extract gradients (note the conversion from forces to gradients)
+    awk -v natom="$natom" '{
+      if ($2 == "Forces") {
+        for (i = 0; i < natom; i++) {
+          getline;
+          printf("%3.15e %3.15e %3.15e\n", -$2, -$3, -$4)
+        }
+      }
+    }' $dftb_out >> $engrad_file
+}
 
+##### LET'S GO! #####
 cd "$WORKDIR" || exit 2
 
-# Everything is ready, run dftb+!
+prepare_dftb_inputs
+
 $DFTBEXE &> dftb.out
 if [[ $? -ne 0 ]]; then
-   echo "ERROR: DFTB calculation failed."
-   echo "See $PWD/dftb.out" 
-   exit 2
+    echo "ERROR: DFTB calculation failed."
+    echo "See file '$PWD/dftb.out' to find out why."
+    exit 2
 fi
 
-# Extract energy from dftb+ output
-grep 'Total energy:' detailed.out | awk '{print $3}' > ../../engrad.dat."$ibead"
-
-# Extract gradients (note the conversion from forces to gradients)
-awk -v natom="$natom" '{
-  if ($2 == "Forces") {
-    for (i = 0; i < natom; i++) {
-      getline;
-      printf("%3.15e %3.15e %3.15e\n", -$2, -$3, -$4)
-    }
-  }
-}' detailed.out >> ../../engrad.dat."$ibead"
-
-# Copy the charges for the next step
-if [[ -f "$CHARGES_FILE" ]]; then
-  cp "$CHARGES_FILE" ../
-fi
+extract_energy_and_gradients ../../engrad.dat."$ibead"

From 5742634042709b88fb8020aa26ed98ca34aca955 Mon Sep 17 00:00:00 2001
From: Daniel Hollas <daniel.hollas@bristol.ac.uk>
Date: Tue, 15 Oct 2024 13:32:40 +0100
Subject: [PATCH 05/10] Add note about the dftb+ executable

---
 interfaces/DFTB/r.dftb | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/interfaces/DFTB/r.dftb b/interfaces/DFTB/r.dftb
index 9fee2fda..08415694 100755
--- a/interfaces/DFTB/r.dftb
+++ b/interfaces/DFTB/r.dftb
@@ -4,6 +4,9 @@
 # The input DFTB parameters are given in a file `dftb_in.hsd` in this directory,
 # which you need to modify for your needs.
 #
+# NOTE: This script assumes that the 'dftb+' executable is already in your PATH.
+# If that is not the case, modify the variable 'DFTBEXE' below accordingly.
+#
 # Tested with version 24.1. Other versions might work, but always test
 # by verifying energy conservation in a short NVE simulation.
 # Versions older than 20.1 will not work since they do not support XYZ geometry input.

From 7c285b84195f3dc4ca3cf5852821704dfde83784 Mon Sep 17 00:00:00 2001
From: Daniel Hollas <daniel.hollas@bristol.ac.uk>
Date: Tue, 15 Oct 2024 14:53:18 +0100
Subject: [PATCH 06/10] Fixes

---
 interfaces/DFTB/dftb_in.hsd | 26 +++++++++++++++++++++++---
 interfaces/DFTB/r.dftb      | 23 +++++++++++++----------
 2 files changed, 36 insertions(+), 13 deletions(-)

diff --git a/interfaces/DFTB/dftb_in.hsd b/interfaces/DFTB/dftb_in.hsd
index 6beda72b..3a033e5d 100644
--- a/interfaces/DFTB/dftb_in.hsd
+++ b/interfaces/DFTB/dftb_in.hsd
@@ -3,19 +3,39 @@ Geometry = xyzFormat {
 }
 
 Hamiltonian = DFTB {
+  charge = 0.0
   SCC = Yes
   MaxSCCIterations = 100
-  charge = 0.0
-  #ReadInitialCharges = yes
+
+  # Slater-Koster files need to be downloaded separately from dftb.org
+  # https://dftb.org/parameters/download/
   SlaterKosterFiles = Type2FileNames {
-   Prefix = "/path/to/parameters/3ob-3-1/"
+   Prefix = "/path/to/slater_koster_params/3ob-3-1/"
    Separator = "-"
    Suffix = ".skf"
   }
+
+  # NOTE: These need to be defined for your system and depend on the SK parameters used,
+  # see documentation for the MaxAngularMomentum for more info.
   MaxAngularMomentum {
     O = "p"
     H = "s"
   }
+
+  # NOTE: The parameters below are specific for the 3ob-3-1 SK parameters!
+  # https://dftb.org/parameters/download/3ob/3ob-3-1-cc
+  ThirdOrderFull = Yes
+  HubbardDerivs = {
+    O = -0.1575
+    H = -0.1857
+  }
+  HCorrection = Damping {
+    Exponent = 4.00
+  }
+
+  # Keep the line below! It is needed to automatically
+  # load initial guess charges from previous time step, see r.dftb
+  #ReadInitialCharges = yes
 }
 
 Analysis {
diff --git a/interfaces/DFTB/r.dftb b/interfaces/DFTB/r.dftb
index 08415694..826b6824 100755
--- a/interfaces/DFTB/r.dftb
+++ b/interfaces/DFTB/r.dftb
@@ -21,22 +21,25 @@ else
   DFTBEXE=dftb+
 fi
 
-timestep="$1"
+# timestep var not used in this script
+# timestep="$1"
+# Bead index in PIMD, "001" for classical MD
 ibead="$2"
-geom=../geom.dat."$ibead"
+
+geom="../geom.dat.$ibead"
 natom=$(wc -l < "$geom")
+WORKDIR="CALC.$ibead"
 
 function prepare_dftb_inputs() {
     # Working directory for the DFTB calculation
-    WORKDIR="CALC.$ibead"
     rm -rf "${WORKDIR}".previous
     if [[ -d "$WORKDIR" ]];then
         mv "$WORKDIR" "${WORKDIR}".previous
     fi
     mkdir -p "$WORKDIR"
 
-    echo -e "$atom\n" > "$WORKDIR/geometry.xyz"
-    cat $geom >> "$WORKDIR/geometry.xyz"
+    echo -e "$natom\n" > "$WORKDIR/geometry.xyz"
+    cat "$geom" >> "$WORKDIR/geometry.xyz"
 
     cp dftb_in.hsd "$WORKDIR"
     # Read charge distribution from previous step if available
@@ -51,7 +54,7 @@ function extract_energy_and_gradients() {
     # Extract energy dftb+ output file
     dftb_out="detailed.out"
     engrad_file=$1
-    grep 'Total energy:' $dftb_out | awk '{print $3}' > $engrad_file
+    grep 'Total energy:' $dftb_out | awk '{print $3}' > "$engrad_file"
 
     # Extract gradients (note the conversion from forces to gradients)
     awk -v natom="$natom" '{
@@ -61,14 +64,14 @@ function extract_energy_and_gradients() {
           printf("%3.15e %3.15e %3.15e\n", -$2, -$3, -$4)
         }
       }
-    }' $dftb_out >> $engrad_file
+    }' $dftb_out >> "$engrad_file"
 }
 
 ##### LET'S GO! #####
-cd "$WORKDIR" || exit 2
-
 prepare_dftb_inputs
 
+cd "$WORKDIR" || exit 2
+
 $DFTBEXE &> dftb.out
 if [[ $? -ne 0 ]]; then
     echo "ERROR: DFTB calculation failed."
@@ -76,4 +79,4 @@ if [[ $? -ne 0 ]]; then
     exit 2
 fi
 
-extract_energy_and_gradients ../../engrad.dat."$ibead"
+extract_energy_and_gradients "../../engrad.dat.$ibead"

From 2768a50fdc139f3d0fbc2c1ae1fc354e225f6450 Mon Sep 17 00:00:00 2001
From: Daniel Hollas <daniel.hollas@bristol.ac.uk>
Date: Wed, 16 Oct 2024 16:25:29 +0100
Subject: [PATCH 07/10] Add check and output example

---
 interfaces/DFTB/r.dftb | 31 +++++++++++++++++++++++++------
 1 file changed, 25 insertions(+), 6 deletions(-)

diff --git a/interfaces/DFTB/r.dftb b/interfaces/DFTB/r.dftb
index 826b6824..6d4cf917 100755
--- a/interfaces/DFTB/r.dftb
+++ b/interfaces/DFTB/r.dftb
@@ -11,7 +11,7 @@
 # by verifying energy conservation in a short NVE simulation.
 # Versions older than 20.1 will not work since they do not support XYZ geometry input.
 cd "$(dirname "$0")" || exit 2
-set -u
+set -uo pipefail
 
 if [[ -f ../SetEnvironment.sh ]]; then
   # This is specific to Prague clusters
@@ -51,20 +51,39 @@ function prepare_dftb_inputs() {
 }
 
 function extract_energy_and_gradients() {
-    # Extract energy dftb+ output file
     dftb_out="detailed.out"
     engrad_file=$1
-    grep 'Total energy:' $dftb_out | awk '{print $3}' > "$engrad_file"
+
+    # Extract energy dftb+ output file
+    # We're matching the third column from this line:
+    # ```
+    # Total energy:  -4.0698318096 H     -110.7458 eV
+    # ```
+    match_regex="^Total energy"
+    # The matched line should appear only once in the DFTB+ output
+    match_count=$(grep -E -c "$match_regex" "$dftb_out")
+    if [[ $match_count -ne 1 ]]; then
+        echo "ERROR: Unexpected DFTB output in file '$dftb_out'"
+        echo "Regular expression '$match_regex' was matched $match_count times"
+        echo "This likely indicates a problem with the calculation or incompatible DFTB+ version"
+        exit 2
+    fi
+    grep -E "$match_regex" "$dftb_out" | awk '{print $3}' > "$engrad_file"
 
     # Extract gradients (note the conversion from forces to gradients)
-    awk -v natom="$natom" '{
-      if ($2 == "Forces") {
+    # ```
+    # Total Forces
+    #    1      0.009551273894      0.004605933524      0.000709843407
+    #    2      0.010527153681      0.006652360906      0.002907870190
+    # ```
+    awk -v natom="$natom" '
+      $1 == "Total" && $2 == "Forces" {
         for (i = 0; i < natom; i++) {
           getline;
           printf("%3.15e %3.15e %3.15e\n", -$2, -$3, -$4)
         }
       }
-    }' $dftb_out >> "$engrad_file"
+    ' "$dftb_out" >> "$engrad_file"
 }
 
 ##### LET'S GO! #####

From 12596fa0114f1872610f346678f57df1d479396f Mon Sep 17 00:00:00 2001
From: Daniel Hollas <daniel.hollas@bristol.ac.uk>
Date: Wed, 16 Oct 2024 16:58:33 +0100
Subject: [PATCH 08/10] README fixes

---
 README.md | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/README.md b/README.md
index 9d23c20a..1eec94dd 100644
--- a/README.md
+++ b/README.md
@@ -5,9 +5,9 @@
 ## What is ABIN?
 
 ABIN is a program for performing ab initio molecular dynamics.
-It is a general purpose program that was initially designed to deal with nuclear quantum effects (NQE).
+It is a general purpose program that was initially designed to model nuclear quantum effects (NQE).
 NQE can be most rigirously captured with path integral MD (PIMD), but also within the Quantum Thermostat based on General Langevin Equation framework developed by Michele Cerriotti.
-ABIN can also simulate non-adiabatic events using Surface-hoping algorithm, using either the classical fewest-switches algorithm (FSSH) or simpler Landay-Zener approach which does not require non-adiabatic couplings. The LZ approach can also capture singlet-triplet transitions.
+ABIN can also simulate non-adiabatic events using Surface-hoping algorithm, using either the classical fewest-switches algorithm (FSSH) or simpler Landau-Zener approach which does not require non-adiabatic couplings. The LZ approach can also capture singlet-triplet transitions.
 
 The basic philosophy of ABIN program is simple — 
 while the program itself handles the propagation of the system according to the equations of motion,

From 3a1be80e1c5bc1104f7a97af7006f79b3bf2b32d Mon Sep 17 00:00:00 2001
From: Daniel Hollas <daniel.hollas@bristol.ac.uk>
Date: Wed, 16 Oct 2024 19:23:45 +0100
Subject: [PATCH 09/10] Check forces output format

---
 interfaces/DFTB/r.dftb | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/interfaces/DFTB/r.dftb b/interfaces/DFTB/r.dftb
index 6d4cf917..493929d5 100755
--- a/interfaces/DFTB/r.dftb
+++ b/interfaces/DFTB/r.dftb
@@ -80,6 +80,11 @@ function extract_energy_and_gradients() {
       $1 == "Total" && $2 == "Forces" {
         for (i = 0; i < natom; i++) {
           getline;
+          if ($1 != i || NF != 4) {
+            print "ERROR: Unexpected line in the DFTB+ output file"
+            print $0
+            exit 2
+          }
           printf("%3.15e %3.15e %3.15e\n", -$2, -$3, -$4)
         }
       }

From ecc7288063af868795bf2e2d0ee4c6ff06863925 Mon Sep 17 00:00:00 2001
From: Daniel Hollas <daniel.hollas@bristol.ac.uk>
Date: Wed, 16 Oct 2024 19:31:02 +0100
Subject: [PATCH 10/10] Fix!

---
 interfaces/DFTB/r.dftb | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/interfaces/DFTB/r.dftb b/interfaces/DFTB/r.dftb
index 493929d5..d45869ce 100755
--- a/interfaces/DFTB/r.dftb
+++ b/interfaces/DFTB/r.dftb
@@ -76,19 +76,19 @@ function extract_energy_and_gradients() {
     #    1      0.009551273894      0.004605933524      0.000709843407
     #    2      0.010527153681      0.006652360906      0.002907870190
     # ```
-    awk -v natom="$natom" '
+    awk -v natom="$natom" -v out="$engrad_file" '
       $1 == "Total" && $2 == "Forces" {
-        for (i = 0; i < natom; i++) {
-          getline;
+        for (i = 1; i <= natom; i++) {
+          getline
           if ($1 != i || NF != 4) {
             print "ERROR: Unexpected line in the DFTB+ output file"
             print $0
             exit 2
           }
-          printf("%3.15e %3.15e %3.15e\n", -$2, -$3, -$4)
+          printf("%3.15e %3.15e %3.15e\n", -$2, -$3, -$4) >> out
         }
       }
-    ' "$dftb_out" >> "$engrad_file"
+    ' "$dftb_out"
 }
 
 ##### LET'S GO! #####