Basic Q-Chem Primer

MGCF - College of Chemistry, University of California, Berkeley

This primer is under construction as we change from Gabedit to IQMol.  11/1015.

Q-Chem is a powerful program that can perform many types of molecular wavefunction-based calculations. The recommended graphical user interface is IQMol. In the MGCF, type iqmol to start it. This primer shows the basics of Q-Chem input and analysis, assuming you are already familiar with the fundamental concepts of computational chemistry. Therefore this should NOT be used as a substitute for reading the manual or as a lesson in theory. You can find examples of specific calculations on MGCF workstations in /usr/software/qchem_4.3.1/samples

Contents:

1. Exporting coordinates from Maestro

2. Setting up the calculation in Gabedit

3. Editing/explaining the Q-Chem input file

4. Job submission, output files

5. Using data from a previous calculation

6. Using Gabedit to analyze results

1. Exporting coordinates from Maestro.

Maestro, a user-friendly molecular building program available in this facility (see Maestro tutorials). Any builder that can export pdb, mol2, or xyz format coordinate files can be used at this step. If using Maestro, export the structure by selecting Export Structures from the Project menu. In the Format pulldown, change Maestro to PDB. Name the file in the File field, and click Export.

2. Setting up the calculation in Gabedit.

Start Gabedit by typing gabedit on the command line. When using this program for the first time you may need to go through numerous setup menus, just click continue or OK for all of them.

Note: if you want to be able to submit calculations directly from Gabedit, you need to set the command for the run_qchem script. To do this, select Preferences from the Settings menu, then click the Commands tab. In the Command for execute Q-Chem field, type run_qchem and then type the enter key (important!).

Under the Geometry menu, select Draw. Right click in the black area and select Read -- PDB file, then browse for your file and click Open. Close the Draw Geometry window after making sure the molecule was correctly read in (don't worry about bond orders or lack of drawn bonds, Q-Chem doesn't care about them).

Select File -- New -- Q-Chem Input (or click on the blue icon on the toolbar). Set the charge and multiplicity you want and the job type. You can usually go with the default for SCF Type, but you may want to set Max # of SCF iterations to something like 100. For Exchange Type, select Hartree-Fock or the flavor of DFT you want to use. If you selected an exchange/correlation combination such as B3LYP you should NOT select anything under the Correlation Method pulldown (note that selecting B3LYP actually gives you B3LYP5, which is not the same as the functional used in Gaussian and Jaguar. You need to delete the 5 by hand after making the file if you want the original functional). Finally, select the basis set you want to use (in most cases the Default Mo Guess can be used, see below), and click OK.

The main Gabedit window now displays the input file text. You can edit the input here, or you can save the file and then edit it in a text editor. To save the file, select File -- Save As and enter the name. Gabedit will append the suffix .inp automatically. To submit the file from Gabedit, select Run -- Run a abtinio program (presumably a typo for ab initio). Make sure Program is set to Q-Chem and Server is set to Local (you will actually be submitting the job to the Chemistry Cluster called "Dino" but the script is run locally). Change Folder to the directory where you want your files saved, and then type the jobname in the Save data in file field (note that .inp is automatically added to the filename). Make sure Command to execute is set to run_qchem, then click OK. Your job will be submitted to Dino the same way as if you had submitted it from the command line, see topic 4 below.

3. Editing/explaining the Q-Chem input file.

The series of example calculations presented here are on a FeF2 molecule obtained from replacing the oxygen and hydrogen atoms in water with iron and fluorine atoms. You can either create a basic input file in Gabedit as in Section 2 and modify that, or simply cut/paste the example input into a text editor (save the file as iron1.inp). The first input file is actually a two-job file that will run the second job automatically, with a total computation time of about 3 minutes on one of our machines. The first job in this file is a single point energy calculation, specially set up to build a good set of orbitals to use in the subsequent geometry optimization job using a non-standard, mixed basis set. Line-by-line explanations of the sections and keywords follow.

$comment
    A single point DFT calculation in with basis set projection
$end

$molecule
0   1
Fe1  -0.015333     0.843000    0.193333
F2   -1.718333    -0.442000   -0.141667
F3    1.733667    -0.401000   -0.051667
$end

$rem
    JOBTYPE   sp
    EXCHANGE   b3lyp
    MAX_SCF_CYCLES   200
    XC_GRID   1
    BASIS2   3-21g
    BASIS   gen
    ECP   gen
    PURECART   11
$end

$basis
    F
    aug-cc-pVDZ
    ****
    Fe
    lanl2dz
    ****
$end

$ecp
    Fe
    lanl2dz
    ****
$end

@@@

$comment
    this is the geometry optimization
$end

$molecule
    read
$end

$rem
    JOBTYPE   opt
    SCF_GUESS   read
    EXCHANGE   b3lyp
    XC_GRID   1
    BASIS   gen
    ECP   gen
    PURECART   11
$end

$basis
    F
    aug-cc-pVDZ
    ****
    Fe
    lanl2dz
    ****
$end

$ecp
    Fe
    lanl2dz
    ****
$end

The Q-Chem input file normally consists of one text file containing all of the information needed for the calculation. Each section starts with $section_name and ends with $end. Input is in free format and is not case sensitive (except for atom symbols); the indentation pattern shown above is just for ease of viewing. Blank lines between sections as well as any text on a line after the necessary input will be ignored as a comment. The $comment section contains a text string that is echoed in Q-Chem output files but no data is read from it. The $molecule section is where the coordinates are specified, in angstroms/degrees (in Cartesian or Z-matrix format, see the manual), on which the first line states the charge and multiplicity of the system (0 1 means neutral singlet). The $rem section lists all of the non-default settings that the program will use in that job. The order in which keywords are specified in this section does not matter. The $molecule and $rem sections are mandatory for all calculations.

For our first calculation JOBTYPE is set to sp for single point and EXCHANGE is set to B3LYP (the popular hybrid DFT functional). Notice that no CORRELATION keyword is specified because we already specified the whole XC functional in the EXCHANGE keyword. To make sure we have enough SCF iterations to achieve convergence, MAX_SCF_CYCLES is set to 200. Additionally, setting the XC_GRID to 1 (0 is default) calls for a finer DFT grid which in this case is necessary for convergence.

The rest of the keywords in this $rem section relate to basis sets. This example shows the use a custom basis set, which introduces a complication: the high-quality SAD initial orbital guess does not work with a custom basis set. To work around this, the program first obtains a wavefunction in a standard basis set via the BASIS2 keyword and then a basis set projection is automatically performed to convert the orbitals from the first calculation to an initial guess for the custom basis set calculation. The 3-21G basis set is usually the best choice because it is available for all elements up through Xe and usually gives better convergence than STO-3G. The reason why we broke this input file up into two jobs is that if we used basis set projection in the geometry optimization it would do it at every geometry step, yielding undesirable results. The BASIS and ECP (effective core potential or pseudopotential) sections both have gen (or general) for an argument in place of a basis set name, meaning input sections for both follow.

Note: if you are using a standard basis set for the whole calculation you do not need the BASIS2 keyword, but you must specify the same basis set name for BASIS and ECP. Finally, the PURECART keyword specifies which orbital types use pure (PURECART=1 for 5 d, 7 f, etc.) vs. Cartesian functions (PURECART=2 for 6 d, 10 f, etc.), with the highest angular momentum listed first (here we've requested pure f and d functions, respectively). This keyword should always be used when using custom basis sets.

The $basis section specifies the user-input basis set. Within the section, the first line contains only the element using the first basis set. It is possible to list basis sets for individual atoms but the format is different, see the manual. The second line contains the name of a standard basis set, or exponents and coefficients for any general basis. Each basis set specification ends with a line containing only the character string ****. The $ecp section is constructed in an identical manner but it only contains the atoms using the ECP.

The character string @@@ signifies the end of the first job specification. The temporary files created for that job containing the molecular orbital coefficients, density matrix, etc., will be saved until they are read in during the next part of the job (after that they are overwritten). Note that the next job will run regardless of whether the first job completed successfully.

The second job input also contains a $molecule section, but read is entered instead of the molecule specification, meaning all of the information is to be read from temporary (restart/checkpoint) files created by the previous job. In the $rem section, the JOBTYPE keyword now has the opt argument, for a geometry optimization. Because we want to use the orbitals from the previous calculation, the SCF_GUESS keyword is used with read as its argument. The remaining input is the same as in the previous job but BASIS2 is NOT specified.

4. Job submission, output files.

To submit your input file to the Chemistry Cluster called "Dino" from the command line, type

run_qchem iron1 save

with or without the .inp suffix. The save setting is optional, but it is needed if you want to use data such as molecular stuctures and orbital coefficients for later calculations (see below). The run command creates the file iron1.job and the job is automatically submitted to cluster's queuing system. While the calculation is running, the only output file generated is iron1.out (standard output) which contains information about the progress of the job and, in the case of geometry optimizations, the atomic coordinates for each structure. The output starts with the Q-Chem reference, then echoes the input. What follows depends on the job, but normally includes atomic coordinates in the standard orientation, information about memory settings, SCF-related output, then MO energies and dipole moments, etc. Output for the job ends with "Thank you very much for using Q-Chem;" however, just because the job completed does not mean it was successful!. You can monitor the optimization while it is running or analyze the results in Gabedit, see section 6, part A.

If you used the save option, a folder called iron1-saved is created in your working directory (the directory from where the job was submitted) when the job finishes. This contains files needed for a restart or plots or a subsequent job using data from the current one. For the example input presented here, this folder will appear after the first part (SP calculation) finishes, and then be overwritten by the second part. See the following section for its usage.

 5. Using data from a previous calculation

Note: To use data from a jobname-saved folder for a subsequent calculation you MUST use the save option when submitting the new job.

(note from Kathy - this changed in the current version so we are in the process of modifying these instructions to match the qchem and run_qchem behavior. This page will be updated shortly).

Part A: orbital plots.

The folder iron1-saved created at the end our optimization contains everything needed to restart a calculation, included the current geometry, molecular orbital coefficients, and the force constant matrix or Hessian. Note: some of these are binary files that cannot be opened with a text editor. To use these data, first make a copy of the folder using the name that you will call the next job:

cp -r iron1-saved iron_orb-saved

where the -r is required because these are folders, not files.

Next, make a copy of iron1.inp called iron_orb.inp, open iron_orb.inp in a text editor and delete the first job specification (also the @@@ characters). In the $rem section, change JOBTYPE to sp and add two lines: PRINT_GENERAL_BASIS with the true argument, and PRINT_ORBITALS also with the true argument. We didn't specify these keywords in the optimization because the extra information would be printed for every optimization iteration.

Change the read line from

$molecule
    read
$end

to point to the location of the molecule file (for example):

$molecule
    read /home/kathy/tests/qchem/iron_orb/molecule
$end

This input file will read the geometry and orbitals from the iron_orb folder and generate everything Gabedit needs to make orbital plots (see section 6, part B). Submit iron_orb.inp to the cluster using: run_qchem iron_orb.inp save

Part B: vibrational frequencies, thermodynamics.

Make a copy of the iron1-saved folder called iron_freq, as in part A. Again, make a copy of iron1.inp called iron_freq.inp and delete the first job specification. Now, change JOBTYPE to freq. Q-Chem will calculate the frequencies and infrared intensities. If you want Raman intensities also use the DORAMAN keyword with the true argument, but keep in mind that this will take much longer and ECPs are not supported. Submit iron_freq.inp to the cluster using run_qchem.

6. Analyze results.

Part A: Geometry optimization.

In the main Gabedit window, select File -- Open and browse for iron1.inp (click on the Home button to get to your top-level directory). This opens the input file and also the output file. Click on the output tab (right above the whitespace in the main Gabedit window) to display the output panel. Note: if you are monitoring a job that is still running, you should click the Update button to the right of the output panel so you have the most up-to-date output. Next, click on the Geom. Conv. button. This displays a graph of energy vs. iteration number, which is useful in analyzing the quality or progress of an optimization. Click on the Draw button in the Geometry Convergence window and the geometry at the selected iteration will be displayed (the graph minimizes when you do this, retrieve it from the task bar). Click on each point on the graph and see how the structure changes throughout the optimization. Before moving on, look at the text output for iron1 and make sure the job completed successfully.

Part B: Visualizing molecular orbitals.

This part assumes you have already done part 5, section A of this primer. From the output panel, click on the Dens. Orb. button (you do not need to have any files open yet). Right click in the black space and select Orbitals -- Read geometry and orbitals from a Q-Chem output file and browse for iron_orb.out. Select an orbital from the Orbitals window that pops up, and click OK. The next dialog to come up can generally use the default values, but you probably want to adjust the Value in the following Iso-Value dialog. Smaller numbers generate larger orbital lobes; 0.05 works well for this example. To change the isovalue after the orbital is rendered, right click and select Surfaces -- reset isovalue.

When you want to display another orbital, you must first delete the current one as there is no toggle feature. To do this, right click and select Surfaces -- delete all. If you have a large molecule and you think you want to view that orbital again, save a Gabedit cube file (right click, Cube -- Save) which can be quickly read in later through the Cube menu. Next, right click and select Orbitals -- Selection and then repeat as above. The orbital surface can be rendered transparent by selecting Render -- Surface -- Transparency.

Part C: Visualizing vibrational modes.

This section assumes you have already done part 5, section B. From the output panel, click on the Dens. Orb. button (you do not need to have any files open yet). Right click in the black space and select Animation -- Vibration. This opens the Vibration window, in which you can read in iron_freq.out through the File menu. Select the mode you wish to visualize and click the Play button. The Scale factor field controls the magnitude of the displacements, which also adjusts the length of the displayed displacement vectors.

You can create an animation of a mode by selecting Create a film before pressing Play (you should probably change BMP to something smaller like JPG). For this purpose a good number of steps by cycle is 16. This will generate 32 frames called gab##.jpg. Next, change the names of the first 9 frames from gab#.jpg to gab0#.jpg. To convert the frames to a gif movie, type the following on the command line:

convert -delay 10 -loop 1000 gab*.jpg iron_animate.gif

This gives you a gif-format movie (called iron_animate.gif in this case) that shows each frame for 10 milliseconds and loops 1000 times before stopping. You can adjust these to suit your purposes.