13
$\begingroup$

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.

$\endgroup$
2
  • 3
    $\begingroup$ The fastest way is probably to use Bash scripts. It takes you at most a few hours to learn Bash to the extent that you can write that script by yourself. $\endgroup$
    – wzkchem5
    May 23 at 20:07
  • 1
    $\begingroup$ Which software do you want to run the calculations in? Then I might be able to write a script for that. $\endgroup$
    – S R Maiti
    May 24 at 9:03
13
$\begingroup$

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.

$\endgroup$
1
  • 1
    $\begingroup$ Thank you, in the past I wrote some bash script but I didn't know sed program to modify text files. $\endgroup$
    – NickZ
    May 24 at 21:57
12
$\begingroup$

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

$\endgroup$
4
  • $\begingroup$ Thank you, I knew that a python script would be one solution but unfortunately I don't know anything about python programming. Does the last part of the code act also as a queuing system to run code sequentially when the previous finished? $\endgroup$
    – NickZ
    May 24 at 22:05
  • $\begingroup$ @NickZ The part where the QM program is executed is inside the for loop and the QM program is called with os.system() which invokes the shell/command prompt. The os.system() will pause the program until the terminal returns. So, whether the running is sequential or not depends on the program. For most programs I have seen like GAMESS/Orca, the terminal doesn't return until the calculation is complete, in those cases, the Python script will run sequentially. However, if you are using a job scheduler and submit scripts etc. then the loop won't stop. I have edited the answer to include this. $\endgroup$
    – S R Maiti
    May 25 at 9:04
  • 1
    $\begingroup$ Do you want to try showing what the input and output should look like on GitHub like I did here: github.com/HPQC-LABS/Modeling_Matters/tree/master/5040/sed ? I made folders for my sed solution and my MATLAB solution, so you could make a folder called Python and show us what the files look like before and after. Maybe in the next 70 years of our lives, many people will benefit from being able to download those when they want to run a batch of calculations! $\endgroup$ May 25 at 15:37
  • $\begingroup$ @NikeDattani I have added a github link for the input and output files, although I don't know if Python will even exist after 70 years!! $\endgroup$
    – S R Maiti
    May 26 at 9:38
8
$\begingroup$

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$Pictorial representation of snakemake workflow

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)
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.