BAGEL interface for XMS and MS-CASPT2 in FSSH and LZSH (#153)

The BASH interface for BAGEL to perform FSSH and LZSH simulations. It supports XMS-CASPT2, MS-CASPT2 (not advisable) and CASSCF, using analytical NACs. It was tested on cyclopropanone molecule and benchmarked to MOLPRO. Credits to Tomas Jira, Jiri Suchan, Jiri Janos The scripts were also extensively tested on cyclobutanone molecules. The scripts support detection of common issues and automatic restarts with modified thresholds. The scripts were linted by shellcheck by running: shellcheck -a -S info -x bagel.inp r.bagel-sh Co-authored-by: Jiří Janoš <[email protected]> Co-authored-by: Daniel Hollas <[email protected]>
PHOTOX · Jan 15, 2024 · f25a007 · f25a007
1 parent 81fdc94
commit f25a007
Show file tree

Hide file tree

Showing 14 changed files with 1,282 additions and 5 deletions.
diff --git a/dev_scripts/install_mpich.sh b/dev_scripts/install_mpich.sh
@@ -51,6 +51,9 @@ cd $MPICH_DIR/$MPICH_VERSION/src && tar -xzf ../pkg/${TAR_FILE} && cd mpich-${MP
 # Use the two rows below for a debug build/
 # export CFLAGS='-g -O0'
 # --disable-fast --enable-g-option=all \
+# These are needed for newer compilers
+# export FFLAGS=-fallow-argument-mismatch
+# export FCFLAGS=-fallow-argument-mismatch
 ./configure FC=gfortran CC=gcc \
   --enable-fortran=all \
   --with-pm=hydra --with-device=ch3:nemesis \

diff --git a/interfaces/BAGEL-LZ/bagel.inp b/interfaces/BAGEL-LZ/bagel.inp
@@ -0,0 +1,29 @@
+#!/bin/bash
+# NOTE: This script is meant to be sourced, not executed!
+
+### Electronic structure parameters
+# NOTE: Number of states is specified in ABIN input file!
+basis=cc-pvdz
+df_basis=cc-pvdz-jkfit
+nspin=1                 # 0 singlet, 1 dublet, 2 triplet ...
+charge=1
+hf_variant=rohf         # hf/rohf/uhf;  hf stands for restricted hf
+method=xms_caspt2       # xms_caspt2/ms_caspt2/sa_casscf
+nact=4
+nclosed=6
+maxiter_CASSCF=250      # maximum number of CASSCF iterations
+maxiter_CASPT2=250      # maximum number of CASPT2 iterations
+shift=0.25              # imaginary shift for caspt2 (default=0.25)
+
+### Setup BAGEL executable and number of cores 
+# You might need to modify how to launch BAGEL executable,
+# depending on how you compiled it and you cluster setup.
+# NOTE: OpenMP parallelization via BAGEL_NUM_THREADS and MKL_NUM_THREADS
+# does not seem very effective, prefer MPI parallelization via mpirun -n <NPROC>
+# export BAGEL_NUM_THREADS=1
+# export MKL_NUM_THREADS=1
+# export BAGELEXE="mpirun -n 1 BAGEL"
+
+# PHOTOX-specific setting
+# shellcheck source=/dev/null  # See https://www.shellcheck.net/wiki/SC1091
+source SetEnvironment.sh BAGEL
diff --git a/interfaces/BAGEL-LZ/bagel_common.sh b/interfaces/BAGEL-LZ/bagel_common.sh
@@ -0,0 +1,204 @@
+#!/bin/bash
+# Common functions for BAGEL BASH interface
+# NOTE: This file is meant to be sourced, not executed!
+
+function exec_bagel {
+    local input_file=$1
+    local output_file=$2
+    local timestep=$3
+    # BAGEL does not print its input so we will prepend it.
+    cp "$input_file" "$output_file"
+    $BAGELEXE "$input_file" >> "$output_file" 2>&1
+    local exitcode=$?
+
+    # Concatenate all output files
+    all_output="$output_file".all
+    sep="\n############################\n"
+    header="${sep}TIMESTEP = $timestep\n$(date)${sep}"
+    echo -e "${header}" >> "$all_output"
+    cat "$output_file" >> "$all_output"
+    return $exitcode
+}
+
+function bagel_error {
+    local errmsg="$1"
+    local bagel_output="$2"
+    local copy="$3"
+    if [[ -f "$bagel_output" && -n "$copy" ]]; then
+        # NOTE: There might be multiple errors per timestep so concatenate
+        cat "$bagel_output" >> "$copy"
+    fi
+    >&2 echo -e "${errmsg}"
+    >&2 echo "Inspect file $(dirname "$0")/$copy"
+    exit 2
+}
+
+function file_exists {
+    if [[ ! -f $1 ]];then
+        >&2 echo "File \"$1\" does not exist!"
+        exit 2
+    fi
+}
+
+### FUNCTIONS FOR GENERATING BAGEL INPUT
+
+# This function must be called first to start the bagel input
+function specify_molecule {
+    # XYZ coordinates provided by ABIN in geom.dat
+    local xyz=$1
+    local basis=$2
+    local df_basis=$3
+    local inp=$4
+    local natom
+    natom=$(wc -l < "$xyz")
+
+    # NOTE: $basis and $df_basis MUST be defined at the top of the file!
+    cat > "$inp" << EOF
+{"bagel": [
+  {
+    "title": "molecule",
+    "basis": "$basis",
+    "df_basis": "$df_basis",
+    "angstrom": "true",
+    "geometry": [
+EOF
+
+    # Geometry specification in json format
+    awk -v natom="$natom" '{printf "\t{\"atom\": \"%s\", \"xyz\": [ %.10f, %.10f, %.10f ]}", $1, $2, $3, $4} NR != natom {print ","}END{print ""}'  "$xyz" >> "$inp"
+    echo -e "    ]\n  }," >> "$inp"
+}
+
+function generate_hf_orbitals {
+    local hf_variant=$1
+    local charge=$2
+    local nspin=$3
+    local inp=$4
+    cat >> "$inp" << EOF
+  {
+    "title" : "$hf_variant",
+    "charge" : $charge,
+    "nspin" : $nspin
+  },
+EOF
+}
+
+function load_orbitals {
+    local orbfile=$1
+    local inp=$2
+    cat >> "$inp" << EOF
+  {
+    "title": "load_ref",
+    "file": "$orbfile",
+    "continue_geom": false
+  },
+EOF
+}
+
+function save_orbitals {
+    local orbfile=$1
+    local inp=$2
+    cat >> "$inp" << EOF
+  {
+    "title": "save_ref",
+    "file": "$orbfile"
+  },
+EOF
+}
+
+function print_molden {
+    local moldenfile=$1
+    local inp=$2
+    cat >> "$inp" << EOF
+  {
+    "title": "print",
+    "file": "$moldenfile",
+    "orbitals": true
+  }
+EOF
+}
+
+# Specify for which state we calculate the forces
+# (for now we only support one target state)
+# This needs to be called before the CAS section.
+function print_forces {
+    local target_state=$1
+    local inp=$2
+    cat >> "$inp" << EOF
+  {
+    "title" : "forces",
+    "export" : true,
+    "grads" : [
+      {
+        "title": "force",
+        "target": $target_state
+      }
+    ],
+EOF
+}
+
+function print_cas {
+    local method=$1
+    local nspin=$2
+    local charge=$3
+    local nstate=$4
+    local nact=$5
+    local nclosed=$6
+    local thresh_CASSCF=$7
+    local thresh_CASPT2=$8
+    local maxiter_CASSCF=$9
+    local maxiter_CASPT2=${10}
+    local shift=${11}
+    local inp=${12}
+
+    if [[ $method == "xms_caspt2" ]]; then
+
+        print_casscf "caspt2" "$nspin" "$charge" "$nstate" "$nact" "$nclosed" "$thresh_CASSCF" "$maxiter_CASSCF" "$inp"
+        print_caspt2 "true" "$thresh_CASPT2" "$maxiter_CASPT2" "$shift" "$inp"
+
+    elif [[ $method == "ms_caspt2" ]]; then
+
+        print_casscf "caspt2" "$nspin" "$charge" "$nstate" "$nact" "$nclosed" "$thresh_CASSCF" "$maxiter_CASSCF" "$inp"
+        print_caspt2 "false" "$thresh_CASPT2" "$maxiter_CASPT2" "$shift" "$inp"
+
+    elif [[ $method == "sa_casscf" ]]; then
+
+        print_casscf "casscf" "$nspin" "$charge" "$nstate" "$nact" "$nclosed" "$thresh_CASSCF" "$maxiter_CASSCF" "$inp"
+        # Remove extra dangling comma
+        sed -i '$ s/,$//' "$inp"
+
+    else
+
+        >&2 echo "ERROR: Unknown method ($method). Specify one of \"xms_caspt2\", \"ms_caspt2\" or \"sa_casscf\" in bagel.inp"
+        exit 2
+
+    fi
+}
+
+function print_casscf {
+    cat >> "$9" << EOF
+    "method": [{
+      "title": "$1",
+      "nspin": $2,
+      "charge": $3,
+      "nstate": $4,
+      "nact": $5,
+      "nclosed": $6,
+      "thresh": $7,
+      "maxiter": $8,
+EOF
+# NOTE: We have to leave this JSON section unclosed since we might need to append CASPT2 input
+}
+
+function print_caspt2 {
+    # TODO: Make it possible to use a real shift instead of imaginary
+    cat >> "$5" << EOF
+      "smith": {
+        "method": "caspt2",
+        "xms": $1,
+        "thresh": $2,
+        "maxiter": $3,
+        "shift": $4,
+        "shift_imag": true
+      }
+EOF
+}