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. 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, diff --git a/interfaces/DFTB/dftb_in.hsd b/interfaces/DFTB/dftb_in.hsd index 482fa5b0..3a033e5d 100644 --- a/interfaces/DFTB/dftb_in.hsd +++ b/interfaces/DFTB/dftb_in.hsd @@ -1,24 +1,45 @@ -Geometry = GenFormat { - <<< "geom_in.gen" +Geometry = xyzFormat { + <<< "geometry.xyz" } Hamiltonian = DFTB { - SCC = Yes charge = 0.0 - #ReadInitialCharges = yes + SCC = Yes + MaxSCCIterations = 100 + + # Slater-Koster files need to be downloaded separately from dftb.org + # https://dftb.org/parameters/download/ SlaterKosterFiles = Type2FileNames { - Prefix = "/home/hollas/3ob-2-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 } -Options { -CalculateForces = Yes +Analysis { + CalculateForces = Yes } ParserOptions { diff --git a/interfaces/DFTB/r.dftb b/interfaces/DFTB/r.dftb index 9aa97cfc..d45869ce 100755 --- a/interfaces/DFTB/r.dftb +++ b/interfaces/DFTB/r.dftb @@ -1,66 +1,106 @@ #!/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 -uo pipefail -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 var not used in this script +# timestep="$1" +# Bead index in PIMD, "001" for classical MD +ibead="$2" + +geom="../geom.dat.$ibead" +natom=$(wc -l < "$geom") +WORKDIR="CALC.$ibead" -# END OF USER INPUT - print natom,"C" - for (i=1;i<=nid;i++) { - printf"%s ",id[i] - } - print "" - print "" +function prepare_dftb_inputs() { + # Working directory for the DFTB calculation + rm -rf "${WORKDIR}".previous + if [[ -d "$WORKDIR" ]];then + mv "$WORKDIR" "${WORKDIR}".previous + fi + mkdir -p "$WORKDIR" + + 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 + 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() { + dftb_out="detailed.out" + engrad_file=$1 + + # 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) + # ``` + # Total Forces + # 1 0.009551273894 0.004605933524 0.000709843407 + # 2 0.010527153681 0.006652360906 0.002907870190 + # ``` + awk -v natom="$natom" -v out="$engrad_file" ' + $1 == "Total" && $2 == "Forces" { + 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) >> out + } } - } -}' ../$geom > geom_in.gen + ' "$dftb_out" +} -# 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! ##### +prepare_dftb_inputs -rm -f detailed.out +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 - 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"