Polarizable Continuum Model (PCM)

Carlo: How about splitting this tutorial in two parts, so it better fits the categorization or the other tutorials? I propose part 1: "Solvation energy of a molecule in solution within the Polarizable Continuum Model", and part 2: "Time propagation in a solvent within the Polarizable Continuum Model". I have already modified the intros, but I haven't split the page yet, so that Gabriel can finish it first

Solvation energy of a molecule in solution within the Polarizable Continuum Model

In this tutorial we show how calculate the solvation energy of the Hydrogen Fluoride molecule in water solution by setting up a ground state calculation using the Polarizable Continnum Model (PCM).

At the time being PCM is only implemented for TheoryLevel = DFT.

Input

The required input keywords to activate a gs PCM run are PCMCalculation and PCMStaticEpsilon.

PCMCalculation just enables the PCM part of the calculation. PCMStaticEpsilon contains the value of the relative static dielectric constant of the solvent medium (in this particular example water with $\epsilon = 78.39$).

The corresponding minimal input for the ground state calculation in solution is


 CalculationMode = gs
 UnitsOutput = eV_Angstrom

 Radius = 3.5*angstrom
 Spacing = 0.18*angstrom

 PseudopotentialSet = hgh_lda

 FH = 0.917*angstrom
 %Coordinates
  "H" |   0 | 0 | 0
  "F" |  FH | 0 | 0
 %

 PCMCalculation = yes
 PCMStaticEpsilon = 78.39

Warnings

Output

The main output from such a gs PCM calculation is the solvation energy of the system, which can be extracted from the file

The actual info file produced from a calculation with the previous output contains the following energetic contributions:


 Energy [eV]:
       Total       =      -680.55034366
       Free        =      -680.55034366
       -----------
       Ion-ion     =       109.92094773
       Eigenvalues =      -120.75757382
       Hartree     =       708.30289945
       Int[n*v_xc] =      -176.84277019
       Exchange    =      -121.31897823
       Correlation =       -13.43829837
       vanderWaals =         0.00000000
       Delta XC    =         0.00000000
       Entropy     =         0.00000000
       -TS         =        -0.00000000
       E_e-solvent =         1.44141920
       E_n-solvent =        -2.34199070
       E_M-solvent =        -0.90057151
       Kinetic     =       556.30947697
       External    =     -1919.42310315
       Non-local   =        42.69395221

E_M-solvent is the solvation energy, which is in turn split in a term due to electrons and another one due to the nuclei, E_e-solvent and E_n-solvent, respectively.

Carlo: add the first few lines of each output file

When we perform the calculation of HF in water using the previous input, a new folder called pcm is created inside the working directory where input file lies. This pcm folder contains the following files:

The first lines of the files are:


    62

   H        1.46477445     1.44399256    -0.00000000
   H        1.08627161     1.44399256     0.52096446
   H        0.47384116     1.44399256     0.32197374

Below, an example of the file:


# Configuration: Molecule + Solvent
###################################
# Epsilon(Solvent) =       78.390
#
# Number of interlocking spheres =    1
#
# SPHERE   ELEMENT               CENTER  (X,Y,Z) (A)                    RADIUS (A)
# ------   -------    -------------------------------------------       ----------
#    1       F        0.91700000      0.00000000      0.00000000        1.54440000
#
# Number of tesserae / sphere = 60
#
#
# Total number of tesserae =   60
# Cavity surface area (A^2) =    29.972947
# Cavity volume (A^3) =    15.430073
#
#
#    iter              E_ee                     E_en                     E_nn                     E_ne                     E_M-solv                 Q_pol                    deltaQ^e                 Q_pol                    deltaQ^n
       1            -293.25112588             294.71995750            -296.89466500             294.71995750              -0.70587588              -7.87251916              -0.02542700               7.89664876              -0.00129741
       2            -293.13633571             294.54433718            -296.89466500             294.54433718              -0.94232634              -7.87358521              -0.02436096               7.89664876              -0.00129741
       3            -293.08877935             294.51502795            -296.89466500             294.51502795              -0.95338845              -7.87305159              -0.02489458               7.89664876              -0.00129741

Carlo: a word of warning about the signs of these charges here

The format of the polarization charge file is:


 cavity position x | cavity position y | cavity position z | polarization charge | label

where | indicates column separation (actually, replace by spaces in the ASC_*.dat files).

Exercises

Convergence of the solvation energy of HF in water solution as a function of the number of SCF iterations.
Convergence of the solvation energy of HF in water solution as a function of the number of SCF iterations.

Advanced settings

Time propagation in a solvent within the Polarizable Continuum Model

In this tutorial we show how perform a time propagation of the Methene molecule in a solvent by setting up a time-dependent calculation and using the Time-dependent Polarizable Continnum Model (TD-PCM)

As in the case of gas-phase calculations, you should perform always the gs calculation before the td one. At the time being PCM is only implemented for TheoryLevel = DFT.

Carlo: I put the different methods into different sections. Used the same naming as in the octopus paper. Pleas add and comment inputs and outputs.

There are several ways to set up a PCM td calculation in Octopus, corresponding to different kinds of approximation for the coupled evolution of the solute-solvent system.

Equilibrium TD-PCM

The first and most elementary approximation would be to consider an evolution where the solvent is able to instantaneously equilibrate with the solute density fluctuations. Formally within PCM, this means that the dielectric function is constant and equal to its static (zero-frequency) value for all frequencies $\epsilon(\omega)=\epsilon_0$. In practice, this approximation requires minimal changes to the gs input file: just replacing gs by td in the CalculationMode leaving the other PCM related variables unchanged. <span style=“color: red>Carlo: Gabriel, please check that this is corresponds with what you mean

Input


 CalculationMode = td
 UnitsOutput = eV_Angstrom

 Radius = 3.5*angstrom
 Spacing = 0.18*angstrom

 PseudopotentialSet = hgh_lda

 FH = 0.917*angstrom
 %Coordinates
  "H" |   0 | 0 | 0
  "F" |  FH | 0 | 0
 %

 PCMCalculation = yes
 PCMStaticEpsilon = 78.39

Output

Inertial TD-PCM

The first nonequilibrium approximation is to consider that there are fast and slow degrees of freedom of the solvent, of which only the former are able to equilibrate in real time with the solute. The slow degrees of freedom of the solvent that are frozen in time, in equilibrium with the initial (gs) state of the propagation. Formally, it means to consider two dielectric functions for the zero-frequency and high-frequency regime (namely, static and dynamic dielectric constants). In practice, this approximation is finally activated by setting extra input keywords, i.e., PCMTDLevel = neq and PCMDynamicEpsilon, the latter to the dynamic dielectric constant value.

Input

Output

Equation of motion PCM

Finally, in the last –and most accurate– nonequilibrium approximation there is an equation of motion for the evolution of the solvent polarization charges, rendering it history-dependent. Formally, this requires to use a specific model for the frequency-dependence for the dielectric function.

Input

Output