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
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
jsub tm big_anion_numfq_cosmo jNumForce -ri -central -cosmo
-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
x2t si9_sncy3_2_2m.xyz > coord
define and skip the first two questions with
In the geometry menu, read in
coord file, determine the point group symmetry (C2v) and generate redundant internal coordinates:
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:
In the general menu, first turn on DFT, choose PBE0 hybrid functional, and m4 integration grid:
Next, turn on resolution of identity (RI) and set RI memory to 0 MB:
Finally, turn on multipole-accelerated RI-J:
You can now exit
define with command *.
Next, add keyword
control file. The fastest way is with
vi, using shortcut
vic (vi control)
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
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:
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:
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.