BAGEL interface for XMS and MS-CASPT2 in FSSH and LZSH

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 \

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

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
+}