Script to generate input files for benchmark purpose

Question

I would like to run some benchmark of different DFT functionals. The most time consuming part is obviously the generation of the input files. I'm not a programmer, but I think this could be done by writing a generic input file with geometry and parameters except for the dft functional, then another text file with a list of the functionals one wants to test. Then I think a script could generate every single input file, maybe appending to it the name of the functional. Hope you have some scripts to share or other ideas to solve my problem.

Answer 1

Consider the MOLPRO input file below template.inp, obtained from here, which was used for a previous MMSE answer about the energy of a zinc atom, except now we're using DFT and I have just entered the word functional instead of the actual functional being used:

 ***,Zn
memory 5880, M
geometry={Zn}
basis = cc-pVDZ
wf,charge=0
{df-rks,functional}

Now consider the file named batch_script.x:

sed s/functional/B3LYP/ B3LYP.inp;
sed s/functional/PBE/   PBE.inp;
sed s/functional/HSE06/ HSE06.inp;

If you run the following commands in the folder containing template.inp and batch_script.x (the first is to make the script file executable, and the second is to run the script):

chmod +x batch_script.x
./batch_script.x

your folder will contain the following three input files:

B3LYP.inp
PBE.inp
HSE06.inp

which look exactly like the template input file shown in the beginning of this answer, but with functional replaced by B3LYP, PBE and HSE06 respectively.

Explanation of the commands in the script file:

• sed s/functional/B3LYP/ B3LYP.inp runs the program sed on B3LYP.inp:

• s/functional/B3LYP/ means substitute functional with B3LYP

I've made the above files available on GitHub: They are here in the folder called sed. You'll find template.inp, batch_script.x and the three output files that get produced when you run the .x script. There's also a folder called MATLAB containing the skeleton of an alternative solution I wrote.

Answer 2

This is an alternative approach with Python scripting. For bash scripting, @NikeDattani has already posted an excellent answer.

Let's say you are trying to run the input files with GAMESS. Start with a template file that may look like this:

!  Methane with B3LYP/6-31G(d)
 $BASIS GBASIS=N31 NGAUSS=6 NDFUNC=1 $END
 $CONTRL SCFTYP=RHF RUNTYP=ENERGY DFTTYP=B3LYP $END

 $DATA 
Title
C1
C     6.0    -5.29801     1.38743     0.00000
H     1.0    -4.22801     1.38743     0.00000
H     1.0    -5.65468     2.34391     0.32068
H     1.0    -5.65468     0.63147     0.66800
H     1.0    -5.65468     1.18690    -0.98867
 $END

Now save the template file as template.inp.

This Python script will read the template file and replace "B3LYP" with names of other functionals and write new input files with the name of the functional appended.

import os
import sys

# list of density functionals
dft_list = ["wB97X-D", "M06-2X", "MN15", "REVM11","TPSS"]

# optionally, read this in from a text file where the functionals
# are separated by space

#with open("dft_list.txt", "r") as dft_file_hnd:
#    dft_list=dft_file_hnd.read().replace("\n","").replace("\r","").split(" ")

# Iterate through dft_list
if os.path.isfile("template.inp")==False:
    print("Unable to find template file template.inp")
    print("Template file template.inp is required!")
    print("Exiting...")
    sys.exit()
for functional in dft_list:
    # write the new input files as input_M06-2X.inp etc.
    # Warning! the name of the functional must not contain unusual
    # characters which cannot be used as file names
    new_file_name = "input_" + functional + ".inp"
    if os.path.isfile(new_file_name):
        print("File:",new_file_name," already exists!")
        print("Please delete it before proceeding")
        print("Exiting...")
        break
    new_file_hnd = open(new_file_name,"a")
    with open("template.inp", "r") as temp_file_hnd:
        for line in temp_file_hnd:
            # search and replace B3LYP with new functional
            line = line.replace("B3LYP",functional)
            # write line into new input file
            new_file_hnd.write(line)
    new_file_hnd.close()
    # run the input files with a QM program
    #output files are named output_M06-2X.log etc.
    print("")
    print("Written ",new_file_name)
    out_file_name = "output_" + functional + ".log"
    print("Running input file:",new_file_name)
    command_string="rungms " + new_file_name + " 00 1 >" + out_file_name
    os.system(command_string)
    # This runs input files with gamess.00.exe with 1 processes (uses rungms)

Note that you can also choose to run the input files inside the same script. There are some lines for that at the bottom. The command_string has to be changed depending on your QM software and OS. The QM software is invoked with os.system() of python which calls the command line terminal of that OS (shell for Linux, cmd for Windows). The os.system() pauses the Python script execution until the terminal returns.

So, if the QM program or command you are running doesn't return until the calculation finishes then the script will run the files sequentially, as the os.system() is inside the for loop. However, if you are running a job submit script in an HPC or something similar to that, then the program will not pause to allow the calculation to finish, it will submit all the jobs in quick succession. Another issue with os.system() is that you would have to escape backslashes for path names if you are using it on Windows (because it is based on C system() IIRC). A better option might then be to use os.subprocess.run().

Also note that the names of the functionals might be slightly different in different programs. For example in Gaussian $\text{M06-2X}$ is represented as m062x or M062X, whereas in GAMESS it is represented as m06-2x or M06-2X.

The advantages of using Python is that it is OS-agnostic (for the most part), so it can be run on all systems. (Bash is not natively available on Windows). Another benefit is that Python is more flexible than bash, so if you need some further modification to the input files, you can do that easily. The disadvantage is obviously that you have to have Python installed.

The python script, template file and sample input and output files can be found here: https://github.com/ShoubhikRaj/molecular-modelling/tree/main/matter-modelling/dft-benchmark

Answer 3

To take @Shoubhik's suggestion of Python one step further, you could also consider setting up your benchmark using Snakemake, a Python program/extension designed for building rule-based data analysis workflows. An example of this approach is the ACCDB (A Collection of Chemistry DataBases) project$^{1,2}$, which in addition to collecting geometries from a number of databases provides a simple interface to rerun these calculations with different basis sets or functionals.

The disadvantage of this approach is the need to learn how to implement a snakemake procedure, which relies on knowing some Python plus additional Snakemake specific syntax. The advantage is that it provides a convenient way to quickly extend your benchmark set to include new molecules/functionals/basis sets/etc. It also can be used to control the job submission and data collection processes, entirely automating the creation of the benchmark data.

$\qquad\qquad\qquad\qquad\qquad$

From personal experience, with a decent Python background, it wasn't terribly difficult to make my own snakemake routine.

1. https://github.com/peverati/ACCDB
2. Morgante, P., Peverati, R.. J. Comput. Chem. 2019, 40, 839– 848. DOI: 10.1002/jcc.25761 (or https://arxiv.org/ftp/arxiv/papers/1809/1809.01707.pdf)