Mod:Hunt Research Group/solvent

From wiki
Jump to navigation Jump to search

Introduction to the SMD Model

The SMD model is an implicit solvation model, which treats the solvent as a dielectric continuum. The solute is surrounded by a cavity, the walls of which can be thought of as representing the solvent. The cavity is built by first centering spheres on each solute atom and then combining these spheres. The calculated solvation energy (and the calculated wavefunction) is sensitive to the exact shape of this cavity.

Full details of how the SMD model works (and how it is parameterised) can be found in the original paper (doi: 10.1021/jp810292n). A second important paper (by the same authors) extends the model to ionic liquids, and includes parameters which should be used in the model for some example ionic liquids (doi: 10.1021/jp304365v).

The SMD model attempts to calculate both electrostatic (dependent on the solvent dielectric constant) and non-electrostatic (depends on a range of solvent descriptors, see below) contributions to solvation energy. The electrostatic solvation energy is directly dependent on the electron density of the solute, while the non-electrostatic solvation term depends solely on the solute geometry. The dependence of the electrostatic term on solute density means that using the SMD model will lead to changes in the wavefunction (relative to a gas phase calculation) as the electrostatic term must be taken into account during the scf procedure. If the geometry of a solute is fixed (e.g in a single-point energy calculation) then only the electrostatic terms will affect the calculated wavefunction.

Volumes Two separate cavities are used in the SMD model, one for the electrostatic solvation energy, and one for the non-electrostatic solvation energy.

Carrying Out Calculations With the SMD model

use the keyword scrf (which stands for self-consistent reaction field) is used to carry out SMD calculations.

WARNING Gaussian will not print the parameters (except dielectric constant and refractive index) used to represent pre-defined solvents. There are more internal parameters which you cannot access.

How to simulate a defined solvent environment

The easiest way to do an SMD calculation is by using a pre-defined solvent in gaussian. Gaussian has many previously defined solvent environments. A list is available at the bottom of this page.[1] For example to use the pre-defined water environment simply insert the following keyword into the method line of your input file. The rest of your method line should specify your functional, basis set, optimisation/other type of calculation as usual.

scrf=(smd,solvent=water)

To use a different solvent to water change the solvent=water part to solvent=something else in the list.

For example here is an input file for an SMD calculation using pentane as the solvent.

%chk=Water_Molecule_SMD_Predef_pentane.chk
# hf/3-21g scrf=(smd,solvent=n-pentane)

Title Card Required

0 1
 O                 -0.54187194    1.86371097    0.00000000
 H                  0.41812806    1.86371097    0.00000000
 H                 -0.86232652    2.76864680    0.00000000

How to simulate a generic solvent environment

The SMD model has many parameters. These are already defined inside Gaussian for the list of defined solvents. If you want to use a solvent not on the list e.g. an ionic liquid, you must define these parameters manually. In this case put the following into the method line:

scrf=(smd,solvent=generic)

and append a list of solvent descriptor values after the geometry/connectivity data. For example(for the pentane solvent):

%chk=Water_Molecule_SMD_ReadIn_pentane_all.chk
# hf/3-21g geom=connectivity scrf=(smd,solvent=generic)

Title Card Required

0 1
 O                 -0.54187194    1.86371097    0.00000000
 H                  0.41812806    1.86371097    0.00000000
 H                 -0.86232652    2.76864680    0.00000000

 1 2 1.0 3 1.0
 2
 3

eps=1.837100 epsinf=1.84 HBondAcidity=0.0 HBondBasicity=0.0 SurfaceTensionAtInterface=22.295067 CarbonAromaticity=0.0 ElectronegativeHalogenicity=0.0

Solvent Descriptors

Parameter Symbol Name in Gaussian input file
Dielectric constant ε eps
Index of refraction, squared n2 epsinf
Macroscopic surface tension /cal mol-1 Å-2 γ SurfaceTensionAtInterface
Abraham hydrogen bond acidity parameter Σα2H HBondAcidity
Abraham hydrogen bond basicity parameter Σβ2H HBondBasicity
Fraction of non-hydrogen atoms which are aromatic carbon atoms φ CarbonAromaticity
Fraction of non-hydrogen atoms which are electronegative halogen atoms ψ ElectronegativeHalogenicity

WARNING EPSINF=Refractive index SQUARED

Quality of the approximation

Pay attention to the value of the density lying outside the cavity, i.e. inside the dielectric. In G03 this value is labeled as “error on total polarization charges”, and example from an output is given below. As a rule of thumb this value should be less than 0.05 for the calculation to be acceptable.

 Error on total polarization charges =  0.03697

Solvent Descriptors more Detail

Different solvents are specified by using the experimental solvent descriptors below:

Static Dielectric Constant(ε) - Used for calculation of electrostatic terms. Accounts for solute-solvent electrostatic interactions. Higher values lead to more negative solvation energies. A list of static dielectric constants for common solvents can be found on the Gaussian website, under the SCRF keyword. In addition, for many normal materials a list of dielectric constants is here "www.honeywellprocess.com/library/marketing/tech-specs/Dielectric Constant Table.pdf" Wiki on relative permittivity is a good simple description. The dielectric constant is also known as the "static relative permittivity". The relative permittivity varies with frequency! The key point is that this is the value of the relative permittivity at zero frequency.

Refractive Index (n) - Used for calculation of non-electrostatic terms. The refractive index of a solvent is related to the strength of solute-solvent dispersion forces, higher values of n should lead to more negative solvation energies. Values can be found in the CRC handbook for chemistry (section 3 physical constants of organic compounds,https://hbcp.chemnetbase.com). WARNING its actually the refractive index squared which is entered into the gaussian input file (see below)

Relationship between dielectric and refractive index the dielectric ε=n^{2} is the refractive index squared (both are normally expressed as complex numbers). The dielectric depends on the frequency at which it is evaluated. Dielectric is normally reserved for low frequencies (below 100MHz as frequency tends to zero) and refractive index for high frequency (visible light, ≈0.6 THz as frequency tends to infinity) so is also called the optical dielectric. Hence the dielectric of water can be near 78 at room temperature (low frequency) but much lower ≈1.333 (at optical wavelengths).

Abraham Hydrogen Bond Acidity (α) - Used for calculation of non-electrostatic terms. Also affects the shape of the electrostatic cavity if oxygen atoms are present in the solute, and can therefore affect the electrostatic part of the solvation energy. Used to account for solute-solvent hydrogen bonding where the solvent donates the hydrogen atom. Values for many solvents can be found in the supporting information of the original SMD paper (doi: 10.1021/jp810292n)

Abraham Hydrogen Bond Acidity (β) - Used for calculation of non-electrostatic terms. Used to account for solute-solvent hydrogen bonding where the solvent donates the lone pair. Values for many solvents can be found in the supporting information of the original SMD paper (doi: 10.1021/jp810292n)

Abraham H-bonding parameters: Traditional solvents

  • Platt's developed a method for evaluating alpha and beta
  • Abraham alpha =6.562+6.102*EP(Hnuc) in [2]
  • Abraham beta=-0.101+0.011*ΔE(HF) in [3]
  • do check the literature more methods have been evolving since then

Kamlet-Taft vs Abraham H-bonding parameters: Ionic liquids

  • the SMD model requires Abraham H-bondonding parameters (Σα2H, Σβ2H)
  • however Kamlet-Taft (α, β) measurements are more commonly reported for ILs
  • a relationship between the parameters was investigated, giving the following equations:[4]
  • Σα2H = 0.4098α + 0.0064
  • Σβ2H = 0.6138β + 0.0890

Previously the group has developed a simple method for calculating Kamlet-Taft parameters, and the instructions are here.[5]

Surface Tension (γ)- is used for calculation of non-electrostatic terms. Accounts for the energy required to create a cavity in the solvent (for the solute to fit into). Higher values of γ lead to more positive solvation energies (creating an empty cavity in a solvent is unfavorable). Values (in mN/m) can be found in the CRC handbook of chemistry (section 6-182, surface tensions of common liquids) but need to be multipled by 1.43932 to be converted to cal mol-1 Å-2.

  • surface tension is the only parameter with units, those used in SMD are non-standard cal mol-1Å-2
  • the SI units are Jm-2 or Nm-1
  • typical units are dyn cm-1 where 1 dyn = 1 g cm s-2
  • as we tend to work in kJ/mol the energy part of this becomes not J but J/mol
  • 1 dyn cm-1 = 0.001N m-1 = 0.001J m-2
  • 1 m = 1*1010Å and 1J=0.239cal and 1 mol-1=6.022*1023
  • 1 dyn cm-1 = 0.001*0.239cal*6.022*1023mol-1/(1*102*10Å2
  • and if you think about this 1023 on top line cancels with 1020 on bottom line leaving 103 which cancels with the 0.001=10-3 leaving us with 0.239*6.022=1.439
  • 1 dyn cm-1 = 1.439 cal mol-1 Å-2
  • 1 mNm-1 (m for mili 10-3) = 1.439 cal mol-1 Å-2

Fraction of Solvent Atoms which are Aromatic Carbons (Φ) - Used for calculation of non-electrostatic terms. Accounts for systematic differences in solvation energies between aromatic/non-aromatic solvents. Can be worked out from the solvent chemical formula (hydrogen atoms are ignored). Examples of values for ionic liquids can be found in doi:10.1021/jp810292n

Fraction of Solvent Atoms which are Halogens (Ψ) - Used for calculation of non-electrostatic terms. Accounts for systematic differences in solvation energies between halogenated/non-halogenated solvents. Can be worked out from the solvent chemical formula (hydrogen atoms are ignored). Examples of values for ionic liquids can be found in doi:10.1021/jp810292n

Solvent parameters for example solvents

Solvent ε n α β γ φ Ψ
n-pentane 1.8371 1.3575 0.0 0.0 22.30 0.0 0.0
Decanol 7.5305 1.4372 0.37 0.48 41.04 0.0 0.0
[C4C1Im][SCN] 13.7 1.54 0.18 0.52 68.29 0.231 0.0
[C4C1Im][PF6] 11.4 1.4090 0.266 0.216 70.24 0.1765 0.3529

Parameters for n-pentane, Decanol and [C4C1Im][PF6] are taken from the sources mentioned int he introduction. Parameters for [C4C1Im][SCN] are taken from: (doi: 10.1021/je101184s / doi: 10.1016/j.jct.2012.03.024 / doi: 10.1039/B9NJ00751B). More information on the ionic liquid parameters in the next page.

Lactic Acid

also known as 2-Hydroxypropanoic acid

Name in Gaussian input file Value Reference Comments/calculations
eps 5.0 [6] Actual reported values in literature are PLA, Reported dielectric constant of poly-lactic-acid PLA is about 5 at 340K and drops to 3.7 at 365K. range by [6] and 3.7 also by [7] Honeywell Dielectric Constant Table reports 22 (16ºC)[8] Value for propanoic acid is 3.44 (25ºC) and butanoic asid is 2.98 (14 ºC).CRC Handbook.
epsinf 2.071 [9] n=1.4392 squared=2.071 is reported in CRC Handbook. Based on PLA but reported as that for LA, n has range 1.456-1.439, mid-point 1.465 [10] n=1.3980 for butanic acid.
SurfaceTensionAtInterface 54.3 No known surface tension, at 25ºC units in mN m-1: acetic acid (26.60), propanoic acid (26.20), butanoic acid (26.05) and 2-hydroxy butanoic acid (37.7), Fluid Properties, Surface Tension of Common Liquids, CRC Handbook [9] Since there is little variation for these acids, we will take the value for 2-hydroxy butanoic acid (Lactic acid is 2-hydroxy propanoic acid) 37.7 mNm-1 *1.43932 =54.3 cal mol-1 Å-2. surface tension=37.49 butanoic acid.
HBondAcidity (α) 0.65 KT alpha=1.58 Abraham alpha =1.58*0.4098+0.0064=0.6538. alpha=0.60 butanoic acid.
HBondBasicity (β) 0.32 KT beta=0.38 Abraham beta=0.38*0.6138+0.0890=0.32 beta=0.45 butanoic acid.
CarbonAromaticity 0
ElectronegativeHalogenicity 0
eps=5.0 epsinf=2.071 SurfaceTensionAtInterface=54.3 HBondAcidity=1.58 HBondBasicity=0.38 CarbonAromaticity=0 ElectronegativeHalogenicity=0

Propylene Carbonate

Name in Gaussian input file Value Reference Comments/calculations
eps 66.14 [9] static: 64.9 dynamic: 2.02 relative permittivity[11] computational paper, no units provided, no source provided
epsinf 2.013 [9] n=1.4189 squared=2.013
SurfaceTensionAtInterface 54.3 40.9 dynes/cm, at 25ºC, 41.93 dyn/cm at 20ºC [12], 41.1 dynes/cm [13] we have taken the value 41.0 dyn/cm *1.43932 =59.0 cal mol-1 Å-2
HBondAcidity (α) 0.00 [14] from QSPR
HBondBasicity (β) 0.65 [14] from QSPR
CarbonAromaticity 0
ElectronegativeHalogenicity 0
eps=66.14 epsinf=2.013 SurfaceTensionAtInterface=59.0 HBondAcidity=0.0 HBondBasicity=0.64 CarbonAromaticity=0 ElectronegativeHalogenicity=0

Important Information On the Implementation of the SMD model in Gaussian 09

The Treatment of Water

In the SMD model water is treated as a special solvent, in that the non-electrostatic terms are not actually calculated using solvent descriptors (unlike every other solvent). Therefore the correct way to carry out a calculation using the water solvent is to use the gaussian pre-defined solvent (scrf=smd,solvent=water). If water is defined in this way Gaussian will treat it as a special solvent. If ANY parameter is varied for the water solvent then Gaussian will NOT treat the solvent as water. For example, in the input file below (a) uses the pre-defined water solvent, while b) uses the same pre-defined solvent but reads in a dielectric constant (which is actually the one dielectric for the pre-defined water). These two input files lead to different non-electrostatic solvation energies being calculated despite the fact that both use the same solvent descriptors (panel c and d show non-electrostatic solvation energies calculated for input file a and b respectively). The reason is that when Gaussian read in the dielectric constant in b) it changed the solvent type from "water" to "generic" and therefore calculated the solvation energy in the same way it would for any other solvent.

Water SMD Picture.jpg

NOTE: To generate the exact non-electrostatic solvation energy of input file b) by defining all terms yourself the following input file is used.

%chk=Water_Molecule_SMD_ReadIn_Water_GenericSolvent.chk
# hf/3-21g geom=connectivity scrf=(smd,solvent=generic)

Title Card Required

0 1
 O                 -0.54187194    1.86371097    0.00000000
 H                  0.41812806    1.86371097    0.00000000
 H                 -0.86232652    2.76864680    0.00000000

 1 2 1.0 3 1.0
 2
 3

eps=78.3553 epsinf=1.777849 HBondAcidity=0.82 HBondBasicity=0.35 SurfaceTensionAtInterface=0.0 CarbonAromaticity=0.0 ElectronegativeHalogenicity=0.0

Meaning Gaussian 09 does not have a surface tension value for the pre-defined water solvent (103.63 cal mol-1 Å-2 is the surface tension of water).

The importance of using the "special" method of calculating non-electrostatic solvation energy for the water solvent, compared to reading in the correct parameters and treating it as a generic solvent, has been tested for multiple molecules (table below). Single point energy calculations were carried out for each structure (optimised at the HF/3-21 level in the gas phase) using both a user-defined water solvent (all relevant parameters read in, and the solvent treated as "generic") and the pre-defined Gaussian solvent (treated as a "special" solvent).

Solute Non-Electrostatic Solvation energy

(User Defined water, kJ mol-1)

Non-Electrostatic Solvation energy

(Pre-defined water, kJ mol-1)

User Defined - Pre-Defined

Non-Electrostatic Solvation energy (kJ mol-1)

[C4C1Im][SCN] ion pair 40.1 29.3 10.8
[C4C1Im]+ Cation 31.5 24.3 7.2
[SCN]- Anion 20.7 14.3 6.4
Methane 13.2 11.4 1.8
Water 19.7 6.1 13.1

Kamlett Taft Vs Abraham Parameters

Values for Abraham hydrogen bond acidity and abraham hydrogen bond basicity are sometimes not available, while Kamlett-Taft values are (e.g the case for ionic liquids). As Abraham and Kamelett-Taft parameters both describe the strength of solute-solvent hydrogen bonds they can be related to each other. An empirical relation ship between Kamlett-Taft and Abraham hydrogen bonding parameters (see doi: 10.1021/jp304365v) is:

αAbraham = 0.4098 αKT + 0.0064

βAbraham = 0.6138 βKT + 0.0890

Where αAbraham and βAbraham are the Abraham hydrogen bond acidity and basicity respectively. αKT and βKT are the Kamlett-Taft hydrogen bond acidity and basicity respectively. Note that the relationship between Kamlet-Taft and Abraham parameters is not perfect, and it is therefore best to Abraham parameters which have been directly measured by experiment rather than converting Kamlet-Taft parameters to Abraham parameters.

Modifying the non-electrostatic Cavity

The cavity used to calculated electrostatic solvation energy can be modified using the normal cavity control keywords on the gaussian website (scrf keyword). The main reason i've done this is just to figure out what Gaussian 09 uses as defaults. Generally you want to use the default values, but there may be a situtation (although i cant think of one) where you'd want to modify the non-electrostatic cavity. To change the atomic radii used to construct the non-electrostatic cavity the keyword "NElRadii=" is used (and appended to the bottom of the input file. For example to use atomic radii from the UFF model the input file is:

%chk=Water_Molecule_SMD_ReadIn_NECavity_Testing3.chk
# hf/3-21g geom=connectivity scrf=(smd,solvent=n-pentane)

Title Card Required

0 1
 O                 -0.54187194    1.86371097    0.00000000
 H                  0.41812806    1.86371097    0.00000000
 H                 -0.86232652    2.76864680    0.00000000

 1 2 1.0 3 1.0
 2
 3

NElRadii=UFF

Note: Gaussian seems to use Bondi Radii as default (as is suggested in the 2009 SMD paper) for the non-electrostatic cavity.

The non-electrostatic cavity is (by default, and as recommended in the 2009 SMD paper) a solvent accesible surface, meaning the the solvent radius is added to each of the atomic centered sphere when generating the cavity. For practical purposes this means that the shape of the non-electrostatic cavity depends on the solvent radius. The solvent radius used for the non-electrostatic cavity is modified with the keyword "SMDRSolvCDS=" using scrf=(smd,solvent=xxxxx,read) and inputting the value at the bottom of an input file. For example:

%chk=Water_Molecule_SMD_pentane_CDSRSolv_0pt3.chk
# hf/3-21g geom=connectivity scrf=(smd,solvent=n-pentane,read)

Title Card Required

0 1
 O                 -0.54187194    1.86371097    0.00000000
 H                  0.41812806    1.86371097    0.00000000
 H                 -0.86232652    2.76864680    0.00000000

 1 2 1.0 3 1.0
 2
 3

SMDRSolvCDS=0.3

Note: Gaussian seems to use a non-electrostatic solvent radius of 0.4Å for all solvents (as is recommended in the 2009 SMD paper).

  1. http://www.gaussian.com/g_tech/g_ur/k_scrf.htm
  2. Phys. Chem. Chem. Phys., 2000, 2, 973-980
  3. Phys. Chem. Chem. Phys., 2000, 2, 3115-3120
  4. Bernales 2012 http://pubs.acs.org/doi/abs/10.1021/jp304365v
  5. http://www.huntresearchgroup.org.uk/research/research_il_alpha_beta_intro.html
  6. 6.0 6.1 Claudius Dichtl, Pit Sippel, and Stephan Krohns, Advances in Materials Science and Engineering Volume, 2017, Article ID 6913835, https://doi.org/10.1155/2017/6913835
  7. Giovanni Spinelli, Rumiana Kotsilkova, Evgeni Ivanov, Radost Ivanova, Carlo Naddeo and Vittorio Romano, Polymers 2020, 12(10), 2414; https://doi.org/10.3390/polym12102414
  8. Honeywell https://www.google.com/url?sa=t&rct=j&q=&esrc=s&source=web&cd=&ved=2ahUKEwjKntHato70AhXEb30KHWG4Ct0QFnoECA8QAQ&url=https%3A%2F%2Fwww.honeywellprocess.com%2Flibrary%2Fmarketing%2Ftech-specs%2FDielectric%2520Constant%2520Table.pdf&usg=AOvVaw1xdJAaUad4jOLmfh2bmZI6
  9. 9.0 9.1 9.2 9.3 CRC_Handbook https://hbcp.chemnetbase.com/faces/documents/06_32/06_32_0001.html
  10. Garlotta Garlotta D; Journal of Polymers and the Environment 9(2): 63-84 (2001). https://link.springer.com/article/10.1023/A:1020200822435
  11. K. Hernández-Burgos , S. E. Burkhardt , G. G. Rodríguez-Calero , R. G. Hennig and H. D. Abruña, J. Phys. Chem. C 2014, 118, 12, 6046–6051
  12. https://macro.lsu.edu/howto/solvents/Propylene%20Carbonate.htm
  13. A. W. Adamson and A. P. Gast, in Physical chemistry of surfaces (Wiley) 6th ed.
  14. 14.0 14.1 J. Chem. Inf. Comput. Sci. 2002, 42, 6, 1320–1331