DFTB+ interface upgrade

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://git.luolix.topmunity/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.

interfaces/DFTB/dftb_in.hsd

@@ -1,13 +1,14 @@
-Geometry = GenFormat {
- <<< "geom_in.gen"
+Geometry = xyzFormat {
+ <<< "geometry.xyz"
 }
 
 Hamiltonian = DFTB {
   SCC = Yes
+  MaxSCCIterations = 100
   charge = 0.0
   #ReadInitialCharges = yes
   SlaterKosterFiles = Type2FileNames {
-   Prefix = "/home/hollas/3ob-2-1/"
+   Prefix = "/path/to/parameters/3ob-3-1/"
    Separator = "-"
    Suffix = ".skf"
   }
@@ -17,8 +18,8 @@ Hamiltonian = DFTB {
   }
 }
 
-Options {
-CalculateForces = Yes
+Analysis {
+  CalculateForces = Yes
 }
 
 ParserOptions {

interfaces/DFTB/r.dftb

@@ -1,66 +1,79 @@
 #!/bin/bash
-cd $(dirname $0)
-source ../SetEnvironment.sh DFTB
+# 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.
+#
+# 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.
+cd "$(dirname "$0")" || exit 2
+set -u
 
-timestep=$1
-ibead=$2
-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
+if [[ -f ../SetEnvironment.sh ]]; then
+  # This is specific to Prague clusters
+  source ../SetEnvironment.sh DFTB
+else
+  # We assume dftb+ is in PATH already. If not, add it here.
+  DFTBEXE=dftb+
 fi
-mkdir -p $WRKDIR
-cd $WRKDIR
 
-# You have to specify, which elements are present in your system
-# i.e. define array id[x]="element"
-# No extra elements are allowed.
-awk -v natom="$natom" 'BEGIN{
-   id[1]="O"
-   id[2]="H"
-   nid=2    # number of different elements
+timestep="$1"
+ibead="$2"
+geom=../geom.dat."$ibead"
+natom=$(wc -l < "$geom")
+
+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"
 
-# END OF USER INPUT
-   print natom,"C"
-   for (i=1;i<=nid;i++) {
-      printf"%s ",id[i]
-   }
-   print ""
-   print ""
+    echo -e "$atom\n" > "$WORKDIR/geometry.xyz"
+    cat $geom >> "$WORKDIR/geometry.xyz"
+
+    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=1;i<=nid;i++) {
-      if ( $1 == id[i] ) {
-         print NR, i, $2, $3, $4
+
+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)
+        }
       }
-   }
-}'  ../$geom > geom_in.gen
+    }' $dftb_out >> $engrad_file
+}
 
-# Read initial charge distribution
-if [ -e charges.bin ];then
-   sed 's/#ReadInitialCharges/ReadInitialCharges/' ../dftb_in.hsd > dftb_in.hsd
-else
-   cp ../dftb_in.hsd .
-fi
+##### LET'S GO! #####
+cd "$WORKDIR" || exit 2
 
-rm -f detailed.out
+prepare_dftb_inputs
 
-$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
-   exit 2
+$DFTBEXE &> dftb.out
+if [[ $? -ne 0 ]]; then
+    echo "ERROR: DFTB calculation failed."
+    echo "See file '$PWD/dftb.out' to find out why."
+    exit 2
 fi
 
-# Extract energy and gradients
-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_energy_and_gradients ../../engrad.dat."$ibead"