3
$\begingroup$

I am running a single point calculation on a fragment of a molecule using restricted open-shell theory. The fragment is in the triplet spin-state and of pentagonal planar structure (D5h). I ran the calculation with rohf/LANL2DZ and was able to achieve convergence quite easily, however with rob3ylp/LANL2DZ I keep running into the following error:

Convergence failure -- run terminated.

I have added the following keywords to the route section in an attempt to help convergence to no avail:

#P pop=full rob3lyp/LANL2DZ scf=(maxcycles=2048,tight,fermi) use=L506 symmetry=(loose) integral=(superfine) IOp(3/33=1) 

Are there any other keywords that might be useful? I need to use the current geometry of the fragment (it is already optimized for the whole molecule), and if possible not change the basis set.

Any help is greatly appreciated! Thank you!

$\endgroup$
1
  • 1
    $\begingroup$ Does the ROB3LYP wavefunction break the spatial symmetry of the molecule? If yes, then it may be simply impossible to converge the SCF under D5h symmetry, and you may want to try its largest Abelian subgroup (C2v) instead. $\endgroup$
    – wzkchem5
    Aug 2, 2023 at 7:39

1 Answer 1

6
$\begingroup$

Before introducing solutions to solve this problem, let us analyze unnecessary keywords:

pop=full IOp(3/33=1) use=L506 does not affect the SCF convergence at all, but increases the output. Often there is no need to write these two keywords (unless you have some special calculation purposes, e.g. you need the AO overlap integrals).

Increasing the SCF maxcycles will work only when the SCF is close to the convergence threshold and the electronic energy becomes lower gradually in the last few cycles. In such case 50 extra cycles would usually be sufficient. maxcycles=2048 seems crazy and it is possible that the user did not monitor the SCF procedure. One can monitor the SCF convergence by running

grep 'Delta-E=' xxx.log

in Linux command line.

integral=(superfine) usually does some help. I don't think symmetry=(loose) or scf=fermi would be useful.

Using the basis set and ECP LANL2DZ for all atoms is not recommended. For non-metal atoms (H, C, etc), better basis sets like 6-311G(d,p) or def2-TZVP are recommended. If your molecule is large, you can use 6-31G(d) or def2-SVP. While for transition metal atoms, def2-TZVP or cc-pVTZ(-PP) is recommended. If you insist on using the LANL2, you can use LANL2TZ(f), which is much more accurate than LANL2DZ.

Here comes the solutions:

(1) Reading molecular orbitals (MOs) from a previous job. For example, reading MOs of ROHF/LANL2DZ. You can also read MOs from another density functional (e.g. BLYP) if SCF can be converged using that functional. The Gaussian input is like

%oldchk=previous_job.chk
%chk=current_job.chk
#p ROB3LYP/LANL2DZ guess=read geom=allcheck (other keywords you want)

Note that %oldchk can only be used in Gaussian16. For Gaussian09 users, you need to copy/backup a previous chk file and use %chk= to specify the filename. By specifying geom=allcheck, Gaussian program will read Title Card, charge, spin multiplicity, and Cartesian coordinates from the previous chk file, so there is no need to write them in the gjf file.

(2) Try other keywords in Gaussian, e.g. scf(vshift=300).

(3) Read converged ROB3LYP MOs of a smaller basis set. This trick is often used for cc-pVDZ -> cc-pVTZ, def2-SVP -> def2-TZVP, etc. Since you already use a small basis set LANL2DZ, this trick can hardly be used.

(4) Reading converged ROB3LYP MOs from other quantum chemistry packages, e.g. PySCF. Solutions (1), (2) and (3) require some trial and error. Using this solution, however, you will probably succeed by trying only once. Firstly, you need to install the PySCF and MOKIT packages:

https://github.com/pyscf/pyscf
https://pyscf.org
https://gitlab.com/jxzou/mokit

Both of these two packages can be installed via conda install if connected to network. The we write a Python script (since you do not show coordinates of your molecule, let us take the triplet O2 as an example: O2.py):

from pyscf import dft, lib
from mokit.lib.gaussian import load_mol_from_fch
from mokit.lib.py2fch_direct import fchk

lib.num_threads(2) # 2 CPU cores
fchname = 'O2.fch'
mol = load_mol_from_fch(fchname)

mf = dft.ROKS(mol)
mf.verbose = 4
mf.xc = 'b3lypg'
mf.grids.atom_grid = (99,590) # ultrafine
mf.kernel()
# SCF may not converged here. We add an orbital-optimization step using the second-order method

mf = mf.newton()
mf.kernel()
fchk(mf, 'O2_ROB3LYP.fch')

Running this script is simple:

python O2.py >O2.out 2>&1

Then you get the file O2_ROB3LYP.fch. Run

unfchk O2_ROB3LYP.fch O2_ROB3LYP.chk

to obtain the chk file. Now you obtain converged MOs from another package, this will make Gaussian SCF converge immediately (using %chk to read MOs mentioned in Solution (1)).

Moreover, you can do something that Gaussian does not support currently. For example, check the wave function stability of RODFT in PySCF, before transferring MOs to Gaussian:

mf.kernel()
mf.stability()

which makes your wave function more reliable.

$\endgroup$

You must log in to answer this question.