Difference between revisions of "Tutorial:Polarizable Continuum Model (PCM)"

From OctopusWiki
Jump to navigation Jump to search
(15 intermediate revisions by the same user not shown)
Line 1: Line 1:
In this tutorial we will show how to set up ground state and time-dependent calculations for a molecule in solution, using the PCM/TD-PCM Octopus features. For this we will use two examples Hydrogen Fluoride and Methane, the Hydrogen Fluoride for the ground-state calculations and Methane for time-dependent and absorption spectrum calculations.
 
  
By now, PCM is only implemented for TheoryLevel = DFT, and it has been tested properly for CalculationMode = gs and td.
+
<span style="color: red>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 part2: "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</span>
  
== HF diatomic molecule in water: ground state ==
+
= Solvation energy of a molecule in solution within the Polarizable Continuum Model =
  
The important input keywords to activate a gs PCM run are {{variable|PCMCalculation|PCM}} and {{variable|PCMStaticEpsilon|PCM}}.
+
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).
  
PCMCalculation just activates the PCM mode of calculation. PCMStaticEpsilon should be put to the value of the static dielectric constant of the solvent medium (in this particular example, we are trying to mimic water).
+
At the time being PCM is only implemented for TheoryLevel = DFT.
  
By adding the latter keywords to an input for a gs calculation of a molecule in vacuo, one sets the calculation for the molecule in solution. Now, we will check a PCM calculation in practice. To that aim, we use the following input:
+
== Input ==
 +
 
 +
The required input keywords to activate a gs PCM run are {{variable|PCMCalculation|PCM}} and {{variable|PCMStaticEpsilon|PCM}}.
 +
 
 +
{{variable|PCMCalculation|PCM}} just enables the PCM part of the calculation. {{variable|PCMStaticEpsilon|PCM}} contains the value of the relative static dielectric constant of the solvent medium (in this particular example water with <math>\epsilon = 78.39</math>).
 +
 
 +
The corresponding minimal input for the ground state calculation in solution is
  
 
  {{variable|CalculationMode|Calculation_Modes}} = gs
 
  {{variable|CalculationMode|Calculation_Modes}} = gs
Line 28: Line 33:
 
  {{variable|PCMStaticEpsilon|PCM}} = 78.39
 
  {{variable|PCMStaticEpsilon|PCM}} = 78.39
  
The latter input example is intended to make you familiar with the running of Octopus for PCM calculations, it does not provide any benchmark for HF molecule. Actually, the geometry might not be optimal and the gs energy might not be converged with respect to the Radius and Spacing variables. N.B. A new convergence analysis w.r.t. Radius and Spacing (different from the one in vacuo) is in principle due when PCM is activated. N.B. Geometry optimization for solvated molecules within Octopus is still under development; it might work properly only when there are minor adjustments in the geometry when passing from vacuuum to the solvent. On the other hand, the selection of the PseudopotentialSet and the implicit choice of LDA xc functional is arbitrary.
+
=== Warnings ===
 +
 
 +
* Please note that the latter input example is intended to make you familiar with the running of Octopus for PCM calculations, it does not provide any benchmark for HF molecule.  
 +
* Actually, the geometry might not be optimal and the gs energy might not be converged with respect to the Radius and Spacing variables. A new convergence analysis w.r.t. Radius and Spacing (independent of the one in vacuo) is in principle due when PCM is activated.  
 +
* Geometry optimization for solvated molecules within Octopus is still under development; it might work properly only when there are minor adjustments in the geometry when passing from vacuuum to the solvent.  
 +
* The selection of the PseudopotentialSet and the implicit choice of LDA xc functional for this particular example is arbitrary.
  
When we perform the calculation of HF in water using the previous input, a new folder called pcm will be created inside the work directory where input file lies. This pcm folder contains the following files:
+
== Output ==
  
* cavity_mol.xyz - useful way to visualize the molecule inside the cavity tessellation. The cavity is plot as a collection of fictitious Hydrogen atoms. Checking up this file is good to be sure that the molecule fits well inside the cavity and that the tessellation doesn't contain artifacts.
+
<span style="color: red>Carlo: add a few lines explaining the main output of this calculation (i.e. the solvation energy). Point to the info file.</span>
  
The cavity can be improved/refined in several direction if needed by the specific calculation. You can increase the number of points in the cavity by changing the default value of the variable PCMTessSubdivider from 1 to 4. In practice, the tessellation of each sphere centered at each atomic position (but Hydrogen's) from 60 to 240 points. You can also relax the constrain on the minimum distance PCMTessMinDistance (=0.1 A by default) between tesserae so as to obtain a more smooth cavity. You can construct the molecular cavity by putting spheres also in Hydrogen, i.e., by setting PCMSpheresOnH = yes. You can change the radius of the spheres used to build up the cavity by manipulating PCMRadiusScaling. Moreover, there is as well the possibility to read the cavity from file (instead of being generated inside Octopus, default) using PCMCavity keyword, whose value should be the complete path to the cavity file. But keep in mind that usually the default values for these quantities are well supported by extensive numerical checks and published results in the literature.
+
* '''info'''
  
* pcm_info.out - a summary of gs pcm calculation containing information on the tessellation and a table of PCM energetic contributions and total polarization charges due electrons and nuclei for all iterations of the SCF cycle. The PCM energetic contributions are splitted in the interaction energy of: 1) electrons and polarization charges due to electrons E_ee, 2) electrons and polarization charges due to nuclei E_en, 3) nuclei and polarization charges due to nuclei E_nn, and finally, 4) nuclei and polarization charges due to electrons.
+
<span style="color: red>Carlo: add the first few lines of each output file</span>
  
In this regard, one important check is to ascertain wether Gauss theorem obtains or not for the total polarization charges, i.e., if <math>Q_{pol}=-(\epsilon-1)/\epsilon \times Q_M</math> is valid or not, where <math>Q_M</math> is the nominal charge of the molecule. If the latter fails, by setting PCMRenormCharges = yes and manipulating PCMQtotTol Gauss theorem is recover at each SCF iteration (or time-step).
+
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:
  
* pcm_matrix.out - the PCM response matrix. It is printed for advanced debug and specially to check against GAMESS PCM matrix in the same conditions (same solvent dielectric constant, same geometry of the cavity).
+
* '''cavity_mol.xyz''' - useful way to visualize the molecule inside the cavity tessellation. The cavity is plotted as a collection of fictitious Hydrogen atoms. Checking up this file is good to be sure that the molecule fits well inside the cavity and that the tessellation doesn't contain artifacts.
  
* ASC_e.dat - polarization charges due electrons.
+
* '''pcm_info.out''' - a summary of gs PCM calculation containing information on the tessellation and a table of PCM energetic contributions and total polarization charges due electrons and nuclei for all iterations of the SCF cycle. The PCM energetic contributions are split into the interaction energy of: 1) electrons and polarization charges due to electrons E<sub>ee</sub>, 2) electrons and polarization charges due to nuclei E<sub>en</sub>, 3) nuclei and polarization charges due to nuclei E<sub>nn</sub>, and finally, 4) nuclei and polarization charges due to electrons.
* ASC_n.dat - polarization charges due nuclei.
+
 
* ASC_sol.dat - polarization charges due to the full molecule.
+
* '''pcm_matrix.out''' - the PCM response matrix. It is printed for advanced debug and specially to check against GAMESS PCM matrix in the same conditions (same solvent dielectric constant, same geometry of the cavity).
 +
 
 +
* '''ASC_e.dat''' - polarization charges due to the electrons.
 +
 
 +
* '''ASC_n.dat''' - polarization charges due to the nuclei.
 +
 
 +
* '''ASC_sol.dat''' - polarization charges due to the full molecule.
 +
 
 +
<span style="color: red>Carlo: a word of warning about the signs of these charges here</span>
  
 
The format of the polarization charge file is:
 
The format of the polarization charge file is:
  
cavity position x | cavity position y | cavity position z | polarization charge | label  
+
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).
 
where | indicates column separation (actually, replace by spaces in the ASC_*.dat files).
  
Now, we look at the convergence of the total energy in solution in Figure 1. What it is most important to notice in the case of Figure 1 is the fact that there is a stabilization of the system in solvent with respect to the in vacuo case. A clear sign of this is the finite value and sign of the energy difference between the both cases obtained upon convergence (see inset plotted in logscale and different interval).
+
== Exercises ==
 +
 
 +
[[Image:PCM_gs_convergence.png|thumb|420px|Ground state convergence of the total energy of HF in solution|right]]
 +
 
 +
* Look first at the convergence of the total energy of the system. <span style="color: red>Carlo: clarify if you mean the energy of the molecule or the solvation energy, and make the figure accordingly</span> It should look like in the Figure. What it is most important to notice is the fact that there is a stabilization of the system in solvent with respect to the ''in vacuo'' case. Find the final value and sign of the energy difference between the both cases obtained upon convergence.
 +
 
 +
<span style="color: red>Carlo: Please rename to something meaningful such as PCM_gs_convergence.png and make the figure larger. Mediawiki will thumbnail it automatically. Also I am not sure if the inset is useful. Maybe suggest plotting in log scale as an exercise</span>
 +
 
 +
 
 +
* Another important check is to ascertain whether Gauss' theorem is fulfilled for the total polarization charges, i.e., if <math>Q_{pol}=-(\epsilon-1)/\epsilon \times Q_M</math> is valid or not, where <math>Q_M</math> is the nominal charge of the molecule. If the latter fails, by setting {{variable|PCMRenormCharges|PCM}} = yes and manipulating {{variable|PCMQtotTol|PCM}} Gauss' theorem is recovered at each SCF iteration (or time-step).
 +
 
 +
== Advanced settings ==
 +
 
 +
* The cavity can be improved/refined in several ways.
 +
** You can increase the number of points in the cavity by changing the default value of the variable {{variable|PCMTessSubdivider|PCM}} from 1 to 4. In practice, the tessellation of each sphere centered at each atomic position (but Hydrogen's) from 60 to 240 points. Check the convergence and the final value of the total energy with respect to different tessellations.
 +
** You can also relax the constrain on the minimum distance {{variable|PCMTessMinDistance|PCM}} (=0.1 A by default) between tesserae so as to obtain a more smooth cavity. This might be useful when the geometry of the molecule is complicated. Check if and how the total energy changes.
 +
** You can construct the molecular cavity by putting spheres also in Hydrogen, i.e., by setting {{variable|PCMSpheresOnH|PCM}} = yes.
 +
** You can change the radius of the spheres used to build up the cavity by manipulating {{variable|PCMRadiusScaling|PCM}}.
 +
** Moreover, by using {{variable|PCMCavity|PCM}} = 'full path to cavity file' keyword the cavity geometry can be read from a file (instead of being generated inside Octopus, which is the default).
 +
 
 +
= 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.
 +
 
 +
<span style="color: red>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. </span>
 +
 
 +
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 <math>\epsilon(\omega)=\epsilon_0</math>. 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</span>
 +
 
 +
== Input ==
 +
 
 +
{{variable|CalculationMode|Calculation_Modes}} = td
 +
{{variable|UnitsOutput|Execution}} = eV_Angstrom
 +
 +
{{variable|Radius|Mesh}} = 3.5*angstrom
 +
{{variable|Spacing|Mesh}} = 0.18*angstrom
 +
 
 +
{{variable|PseudopotentialSet|System}} = hgh_lda
 +
 +
FH = 0.917*angstrom
 +
%{{variable|Coordinates|System}}
 +
  "H" |  0 | 0 | 0
 +
  "F" |  FH | 0 | 0
 +
%
 +
 
 +
{{variable|PCMCalculation|PCM}} = yes
 +
{{variable|PCMStaticEpsilon|PCM}} = 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., {{variable|PCMTDLevel|PCM}} = neq and {{variable|PCMDynamicEpsilon|PCM}}, the latter to the dynamic dielectric constant value.
  
[[Image:fig1.png|thumb|420px|<sub>4</sub>]]
+
== Input ==
  
== CH4 case: time-dependent and absorption spectrum ==
+
== Output ==
  
Now we move to a PCM td calculation. As in the case of gas-phase calculations, you should perform always the gs calculation before the td one.
+
= Equation of motion PCM =
  
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. The first and most elementary approximation would be to consider an evolution where the solvent is able to equilibrate with the solute density fluctuations in real time. Formally within PCM, this means that the dielectric function is constant and equal to its static (zero-frequency) value for all frequencies. In practice, this approximation requires minimal changes to the gs input file: just replacing gs by td in the CalculationMode.  
+
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.
  
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 ==
  
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.
+
== Output ==
  
  

Revision as of 11:44, 28 November 2019

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 part2: "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 ).

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

  • Please note that the latter input example is intended to make you familiar with the running of Octopus for PCM calculations, it does not provide any benchmark for HF molecule.
  • Actually, the geometry might not be optimal and the gs energy might not be converged with respect to the Radius and Spacing variables. A new convergence analysis w.r.t. Radius and Spacing (independent of the one in vacuo) is in principle due when PCM is activated.
  • Geometry optimization for solvated molecules within Octopus is still under development; it might work properly only when there are minor adjustments in the geometry when passing from vacuuum to the solvent.
  • The selection of the PseudopotentialSet and the implicit choice of LDA xc functional for this particular example is arbitrary.

Output

Carlo: add a few lines explaining the main output of this calculation (i.e. the solvation energy). Point to the info file.

  • info

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:

  • cavity_mol.xyz - useful way to visualize the molecule inside the cavity tessellation. The cavity is plotted as a collection of fictitious Hydrogen atoms. Checking up this file is good to be sure that the molecule fits well inside the cavity and that the tessellation doesn't contain artifacts.
  • pcm_info.out - a summary of gs PCM calculation containing information on the tessellation and a table of PCM energetic contributions and total polarization charges due electrons and nuclei for all iterations of the SCF cycle. The PCM energetic contributions are split into the interaction energy of: 1) electrons and polarization charges due to electrons Eee, 2) electrons and polarization charges due to nuclei Een, 3) nuclei and polarization charges due to nuclei Enn, and finally, 4) nuclei and polarization charges due to electrons.
  • pcm_matrix.out - the PCM response matrix. It is printed for advanced debug and specially to check against GAMESS PCM matrix in the same conditions (same solvent dielectric constant, same geometry of the cavity).
  • ASC_e.dat - polarization charges due to the electrons.
  • ASC_n.dat - polarization charges due to the nuclei.
  • ASC_sol.dat - polarization charges due to the full molecule.

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

Ground state convergence of the total energy of HF in solution
  • Look first at the convergence of the total energy of the system. Carlo: clarify if you mean the energy of the molecule or the solvation energy, and make the figure accordingly It should look like in the Figure. What it is most important to notice is the fact that there is a stabilization of the system in solvent with respect to the in vacuo case. Find the final value and sign of the energy difference between the both cases obtained upon convergence.

Carlo: Please rename to something meaningful such as PCM_gs_convergence.png and make the figure larger. Mediawiki will thumbnail it automatically. Also I am not sure if the inset is useful. Maybe suggest plotting in log scale as an exercise


  • Another important check is to ascertain whether Gauss' theorem is fulfilled for the total polarization charges, i.e., if is valid or not, where is the nominal charge of the molecule. If the latter fails, by setting PCMRenormCharges = yes and manipulating PCMQtotTol Gauss' theorem is recovered at each SCF iteration (or time-step).

Advanced settings

  • The cavity can be improved/refined in several ways.
    • You can increase the number of points in the cavity by changing the default value of the variable PCMTessSubdivider from 1 to 4. In practice, the tessellation of each sphere centered at each atomic position (but Hydrogen's) from 60 to 240 points. Check the convergence and the final value of the total energy with respect to different tessellations.
    • You can also relax the constrain on the minimum distance PCMTessMinDistance (=0.1 A by default) between tesserae so as to obtain a more smooth cavity. This might be useful when the geometry of the molecule is complicated. Check if and how the total energy changes.
    • You can construct the molecular cavity by putting spheres also in Hydrogen, i.e., by setting PCMSpheresOnH = yes.
    • You can change the radius of the spheres used to build up the cavity by manipulating PCMRadiusScaling.
    • Moreover, by using PCMCavity = 'full path to cavity file' keyword the cavity geometry can be read from a file (instead of being generated inside Octopus, which is the default).

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 . 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. 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



Back to Tutorials