Page tree
Skip to end of metadata
Go to start of metadata

Table of contents

What is COSMO?

TURBOMOLE includes the Conductor-like Screening Model (COSMO) which is a continuum solvation model. In this approach, the solvent is represented by a dielectric continuum with permittivity ε and the solute molecule forms a cavity within this continuum (see TURBOMOLE manual chapter 18 for more details and original literature references). 

Using COSMO for anionic systems

One very important use case for the COSMO solvent field is anionic systems with charge higher than –1. Kohn-Sham DFT has major issues in describing such anions and typically yields positive energies for the highest occupied orbitals (that is, the electrons are not actually bound). This issue can be solved by using COSMO solvent field to counter the anionic charge of the system. Adding the dielectric continuum around the anion normally leads in a better description of the electronic structure.

The permittivity ε could in principle set to the literature value of any solvent. However, in particular when the solvent field is only used to counter the anionic charge of the system, it is recommended to keep it in the default value (ε = ∞).

How to use COSMO in a geometry optimization?

To turn on COSMO, first run define normally. After that, add the following keyword to the beginning of the control file:

$cosmo

This will turn on COSMO with default settings. You can submit the geometry optimization normally.

Frequency calculations with COSMO

If you are using COSMO and want to run frequency calculation, it is recommended to run a numerical frequency calculation with NumForce instead of aoforce. You need to use -cosmo switch and you cannot use -c switch:

jsub tm big_anion_numfq_cosmo jNumForce -ri -central -cosmo

The -c switch does not work with -cosmo because COSMO needs data from the gradient check run.

Complete example: [Si9(SnCy3)2]2–

This example is based on [Si9(SnCy3)2]2– published in Chem. Sci., 2019, 10, 9130 (DOI). XYZ coordinates for the example system are available here: si9_sncy3_2_2m.xyz. These coordinates relatively close to the final optimized structure.

Creating input files with define

Start by creating a new empty directory and copy the XYZ file there. Then, convert the XYZ file to a coord file:

x2t si9_sncy3_2_2m.xyz > coord

Run define and skip the first two questions with <Enter>.

In the geometry menu, read in coord file, determine the point group symmetry (C2v) and generate redundant internal coordinates:

a coord
desy 1d-3
ired
*

In the basis set menu, select def-TZVP for C and H atoms and then select def2-TZVP for Si and Sn:

b "c","h" def-TZVP
b "si","sn" def2-TZVP
*

We use larger def2-TZVP basis set for the Si9 cluster and the connected Sn atoms, while for the cyclohexanyl ligands, it is enough to use the slightly smaller def-TZVP basis set. def-TZVP is still clearly better than using a SVP-quality basis set for C and H.

In the occupation number and molecular orbital definition menu, choose EHT guess, accept the default parameters, use -2 as the total charge, and accept the EHT guess:

eht
<Enter>
-2
<Enter>

In the general menu, first turn on DFT, choose PBE0 hybrid functional, and m4 integration grid:

dft
on
func pbe0
grid m4
<Enter>

Next, turn on resolution of identity (RI) and set RI memory to 0 MB:

ri
on
m 0
<Enter>

Finally, turn on multipole-accelerated RI-J:

marij
<Enter>

You can now exit define with command *.

Next, add keyword $cosmo to control file. The fastest way is with vi, using shortcut vic (vi control)

vic
o
$cosmo
<Esc>
:wq
<Enter>

But you can edit the file in any way you like!

Submitting the geometry optimization job

The job can now be submitted:

jsub -np 12 -smp -mem 2G tm si9_sncy3_2_2m_pbe0_tzvp_opt jobex -ri

Warning! The job took three hours to run with 12 CPUs (Intel(R) Xeon(R) CPU E5-2630 v2 @ 2.60GHz).

(If you do not want to actually run the job, you can download the results calculated with TURBOMOLE 7.5: si9_sncy3_2_2m_pbe0_tzvp.tar.gz. Unpack the results with tar xfz si9_sncy3_2_2m_pbe0_tzvp.tar.gz)

While the job is running, it is possible to check how it proceeds by executing

grep cycle gradient

This will print out the total energy and gradient for each geometry optimization step (cycle)

cycle =  1 SCF energy = -4443.7873926390 |dE/dxyz| = 0.001296
...
cycle = 48 SCF energy = -4443.7865256710 |dE/dxyz| = 0.000060

In this case, the optimization started from a rather good geometry, but still took many steps because the potential energy surface is rather flat. COSMO can also make the optimization bit more time-consuming. Another thing to note is that the energy actually slightly increased during the optimization (but gradient decreased). But the increase of 2 kJ/mol is practically nothing for a system of this size and not a reason to worry. This can happen in optimizations with COSMO. 

Analysis of the geometry optimization job

Confirm that you have file GEO_OPT_CONVERGED in the directory. It is also a good practice to check how the optimization has behaved with

grep cycle gradient

Next, run eiger:

eiger > holu.out

and check from holu.out that the HOMO-LUMO gap and orbital energies are reasonable (HOMO-LUMO gap must be positive, energies of occupied orbitals must be negative). For example, in this case everything looks fine:

holu.out
 HOMO-LUMO Separation
   HOMO:   227.    61 b1    -0.17941415 H =     -4.88211 eV
   LUMO:   228.    62 b1    -0.02024713 H =     -0.55095 eV
   Gap :                    +0.15916702 H =     +4.33116 eV

Number of MOs=   1513, Electrons=    454.00, Symmetry: c2v

   Nr.   Orbital    Occupation       Energy
  244.    66 b1                    +0.052851 H =        +1.438 eV
  (some virtual orbitals cut out from here)
  228.    62 b1                    -0.020247 H =        -0.551 eV
  227.    61 b1       2.000        -0.179414 H =        -4.882 eV
  226.    52 b2       2.000        -0.183585 H =        -4.996 eV

Submitting a frequency calculation job

The next task is to run a numerical frequency calculation to confirm that the system is a true local minimum. Submit a numerical frequency calculation with COSMO:

jsub -np 12 -smp -mem 2G tm si9_sncy3_2_2m_pbe0_tzvp_fq jNumForce -ri -central -cosmo

Analysis of the frequency calculation job

After the job has finished, you should check that vibspectrum file does not show any imaginary frequencies:

vibspectrum
$vibrational spectrum
#  mode     symmetry     wave number   IR intensity    selection rules
#                         cm**(-1)        km/mol         IR     RAMAN
     1                       -0.00         0.00000        -       -
     2                       -0.00         0.00000        -       -
     3                       -0.00         0.00000        -       -
     4                       -0.00         0.00000        -       -
     5                        0.00         0.00000        -       -
     6                        0.00         0.00000        -       -
     7        b2             15.34         0.11795       YES     YES

In this case, the system is a true local minimum. Sometimes you can see small imaginary modes with frequencies between, say, –1 and –50 cm−1. You should visualize these modes with for example Jmol or TMoleX. If the modes correspond for example to rotation of some methyl groups or ammonia ligands, there is no reason to worry about them. It might be possible to get rid of them by reoptimizing with even tighter convergence criteria  (jobex -gcart 4), but often they persist even after that, especially if you are also using COSMO. It is OK to transparently report such rotational modes of ligands (and the imaginary frequencies) for example in computational details in a publication.

  • No labels