Saturday, October 4, 2014

Tight binding DFT

Tight binding DFT (DFTB) is a semi-empirical method with speed and accuracy similar to NDDO-based semiempirical methods such as AM1, PM3, and PM6.  Currently there are three types of DFTB methods called DFT1, DFTB2, and DFTB3. DFTB1 and DFTB2 are sometimes called non-SCC DFTB (non-selfconsistent charge) and SCC-DFTB, respectively.  DFTB3 is generally considered the most accurate for molecules and there are several parameter sets for DFTB2 and DFTB3 for different elements.  Compared to PM6, DFTB has so far been parameterized for relatively few elements.

The closed shell DFTB1 energy is computed from the following equation
$$E^{\text{DFTB1}}=\sum_i^{N/2} \sum_\mu^K \sum_\nu^K 2 C_{\mu i} C_{\nu i} H^0_{\mu\nu}+\sum_A \sum_{B>A} E^{\text{rep}}_{AB}$$
where $C_{\mu i}$ is the molecular orbital coefficient for MO $i$ and basis function $\mu$.
$$H^0_{\mu\nu}= \begin{cases}
\varepsilon^{\text{free atom}}_{\mu\mu} & \text{if } \mu=\nu \\
0 & \text{if } A=B, \mu\ne\nu\\
\langle \chi_\mu | \hat{T}+V_{\text{eff}}[\rho_0^A+\rho_0^B] | \chi_\nu \rangle & \text{if } A\ne B
Here,  $\varepsilon^{\text{free atom}}_{\mu\mu}$ is an orbital energy of a free atom, $\chi$ is a valence Slater-type orbital (STO) or numerical orbital, $\hat{T}$ is the electronic kinetic energy operator, $V_{\text{eff}}$ is the Kohn-Sham potential (electron-nuclear attraction, electron-electron repulsion, and exchange correlation), and $\rho_0^A$ is the electron density of neutral atom $A$.

DFT calculations on free atoms using some functional yield $\left\{\varepsilon^{\text{free atom}}_{\mu\mu} \right\}$, $\left\{\chi\right\}$, and $\rho_0$, which are then used to compute $H^0_{\mu\nu}$ for A-B atom pairs at various separations $R_{AB}$ and stored. When performing DFTB calculations $H^0_{\mu\nu}$ is simply computed for each atom pair A-B by interpolation using this precomputed data set.

Similarly, the overlap matrix $\left\{ \langle \chi_\mu | \chi_\nu \rangle  \right\}$ need to orthonormalize the MOs are computed for various distances and stored for future use.

$E^{\text{rep}}_{AB}$ is an empirical repulsive pairwise atom-atom potential with parameters adjusted to minimize the difference in atomization energies, geometries, and vibrational frequencies computed using DFTB and DFT or electronic structure calculations for set of molecules.

So, a DFTB1 calculation is performed by constructing $\mathbf{H}^0$, diagonalizing it to yield $\mathbf{C}$, and then computing $E^{\text{DFTB1}}$.

$$E^{\text{DFTB2}}=E^{\text{DFTB1}}+\sum_A \sum_{B>A} \gamma_{AB}(R_{AB})\Delta q_A\Delta q_B$$
where $\Delta q_A$ is the Mulliken charge on atom $A$ and $\gamma_{AB}$ is a function of $R_{AB}$ that tends to $1/R_{AB}$ at long distances.

The Mulliken charges depend on $\mathbf{C}$ so a selfconsistent calculation is required:
1. Compute DFTB1 MO coefficients, $\mathbf{C}$
2. Use $\mathbf{C}$ to compute $\left\{ \Delta q \right\}$
3. Construct and diagonalize $H_{\mu \nu}$ to get new MO coefficients, $\mathbf{C}$
$$H_{\mu \nu}=H_{\mu \nu}^0 + \frac{1}{2} S_{\mu\nu} \sum_C (\gamma_{AC}+\gamma_{BC})\Delta q_C, \mu \in A, \nu \in B$$
4. Repeat steps 2 and 3 until selfconsistency.

$$E^{\text{DFTB3}}=E^{\text{DFTB2}}+\sum_A \sum_{B>A} \sum_{C>B>A} \Gamma_{AB}\Delta q_A^2\Delta q_B$$
$\Gamma_{AB}$ is computed using interpolation using precomputed data. A SCF calculation is required.

Parameter sets and availability
DFTB is available in a variety of software packages.  I don't believe DFTB3 is currently in Gaussian and DFTB is also available in CHARMM and CP2K.  DFTB will soon be available in GAMESS.

Note that each user/lab must download the parameter file separately here.  There are several parameter sets.  The most popular sets for molecules are the MIO (materials and biological systems) for DFTB2 and 3OB (DFT3 organic and biological applications).

Dispersion and hydrogen bond corrections
Just like DFT and PM6, the DFTB can be corrected for dispersion and hydrogen bond effect.

This work is licensed under a Creative Commons Attribution 4.0

Saturday, September 6, 2014

Four Practical Comments Regarding RHF Calculations

1. If the program you use to perform a geometry optimization at the HF level complains of convergence problems, the first step is to determine whether it is the geometry optimization that failed to converge (i.e. the program could not find a geometry with a low gradient given a certain number of steps) or whether it is the SCF-procedure at a particular geometry that failed to converge (i.e. the program could not find a set or orbitals that have stopped changing given a certain number of steps). It is my experience that if the SCF does not converge in the default number of steps, then simply increasing the allowed number of SCF steps does not solve the problem. The problem is usually either in the initial guess orbitals or there is an error in the molecular structure.

2. Most quantum chemistry programs will have RHF as the default method, which implies a singlet wave function. So if you provide the program with a molecule that has an odd number of electrons (e.g. NH4) the program will complain that it cannot do a singlet calculation on a molecule with an odd number of electrons, and asks you to change either the charge or the multiplicity, i.e. to specify whether you want to do an UHF or ROHF calculation on the doublet state of the NH4 radical, or whether you meant to do a singlet RHF calculation on NH4+. You will get very different results depending on your choice, so the object here is not to use trial and error to find a combination of charge and multiplicity that satisfies the program. Rather you have to decide what makes sense chemically, before you do the calculation.

3. Once you have defined the position of the nuclei and the number of electrons (by specifying the overall charge of the molecule) the SCF procedure does the rest. Therefore, you don’t have to specify that the positive charge “is on the NH3 group” in CH3NH3+. Each orbital has access to all the basis functions in the molecule and the electrons generally “go where they want”. Also, changing the CN bond to a double bond in a graphics program does not affect the result (as long as the CN bond length is the same).

4. This should be obvious, but just in case. It makes no sense to compare the total energy (E) of two molecules (or collection of molecules) that have a different number and kinds of atoms. For example, substituting a hydrogen atom in a molecule with a methyl group will lead to a considerable lowering of the RHF energy. Most of this energy lowering will some from the electron-nuclear attraction energy of the electron in the 1s orbital of the C atom, but says nothing about the relative stability of the two molecules. In the Introduction to Molecular Modeling course I used to teach, I used to subtract 20 points (out of a hundred) for such a comparison in a laboratory report

Tuesday, September 2, 2014

Computational Chemistry Highlights: August issue

The August issue of Computational Chemistry Highlights is out.

CCH is an overlay journal that identifies the most important papers in computational and theoretical chemistry published in the last 1-2 years. CCH is not affiliated with any publisher: it is a free resource run by scientists for scientists. You can read more about it here.

Table of content for this issue features contributions from CCH editors Steven Bachrach, David Bowler, Mario Barbatti, and Jan Jensen:

Interactive Chemical Reactivity Exploration

Saturday, August 9, 2014

Eight practical comments regarding geometry optimizations.

 1. An energy minimization can stop either because a stationary point is found or the program ran out of steps. It’s important that you know why it stopped. If it ran out of steps you have to do another energy minimization starting with the geometry that has the lowest energy. This is not always the last geometry!

2. Larger molecules will require more steps to reach a minimum than smaller molecules. The relationship is roughly linear, so if you double the size of the molecules you will probably take at least twice as many steps, and maybe more. The reason is that each step explores one degree of freedom the molecules, all degrees of freedom must be explored to find the minimum (unless you get lucky), and the number of degrees of freedom increases linearly with the number of atoms (3N-6). This is not quite true for Newton-Raphson energy minimizations if the Hessian is computed, but this is rarely done in practice (see point 5).

3. Floppy molecules with relative flat PESs require more steps than rigid molecules, e.g. a linear hydrocarbon takes longer to optimize than the corresponding cyclic structure. The reason is that flat surfaces mean small gradients and therefore small steps. Also, flat surfaces tend to be non-quadratic so the Newton-Raphson steps have to be scaled back.

4. It is sometimes more efficient to use internal coordinates (bond lengths and angles) for the optimization than Cartesian coordinates, since internal coordinates are more quadratic than Cartesian coordinates. However, if the molecule has many fused rings it might be better to use Cartesians coordinates because of linear-dependence issues.

5. For Newton-Raphson minimizations the Hessian is usually estimated rather than computed exactly because it is expensive to calculate. Attempts are usually made to improve or update the Hessian after each step. So, if you run 50 Newton-Raphson steps from an initial geometry, take the last geometry as the new guess, and run another 50 Newton-Raphson steps you will take a different path towards the minimum than if you run 100 Newton-Raphson steps from the initial geometry. This is because you start with a new guess at the Hessian in step “51” in the first approach. This usually works better in my experience.

6. If you want to find a transition state you must use the Newton-Raphson method, since only the Hessian can tell you in which direction you have to go uphill (the direction is defined by the normal mode corresponding to the imaginary frequency). For the same reason you must compute the Hessian for the initial guess geometry of your transition state.

7. For molecules with more than ~1000 atoms the algorithm needed to invert the Hessian (i.e. produce the inverse of H) becomes very expensive, so steepest descent or conjugate gradient methods are usually used to energy minimize such big molecules.

8. If your starting geometry for ammonium is planar (e.g. in the xy plane), it will stay planar. The reason for this is that the gradient perpendicular to the plane is zero because of symmetry. In essence there is two opposing ways of making ammonia non-planar and the net effect is a cancellation of the force out-of-the plane. Notice that the gradient components in the x and y directions are small (below the optimization threshold) but non-zero, while the gradient in the z-direction is exactly zero. Of course, a zero gradient means no displacement in the out-of-plane direction and ammonia stays planar during the energy minimization.

More generally, the molecule will retain the symmetry of the starting guess structure during the minimization, whether that symmetry corresponds to an energy minimum or not.