A PCM calculation is requested by setting SOLVENT_METHOD = PCM in the $rem section. As mentioned above, there are a variety of different theoretical models that fall within the PCM family, so additional fine-tuning may be required, as described below.
Most PCM job control is accomplished via options specified in the $pcm input section, which allows the user to specify which flavor of PCM will be used, which algorithm will be used to solve the PCM equations, and other options. The format of the $pcm section is analogous to that of the $rem section:
$pcm <Keyword> <parameter/option> $end
Note: The following job control variables belong only in the $pcm section. Do not place them in the $rem section.
Construction of the solute cavity is an important part of the model and users should consult the literature in this capacity, especially with regard to the radii used for the atomic spheres. The default values provided in Q-Chem correspond to the consensus choice that has emerged over several decades, namely, to use vdW radii scaled by a factor of 1.2. The most widely-used set of vdW radii are those determined from crystallographic data by Bondi,106 although the radius for hydrogen was later adjusted to 1.1 Å,795 and radii for those main-group elements not addressed by Bondi were provided later.600 This extended set of vdW is used by default in Q-Chem, and for simplicity we call these “Bondi radii” regardless of whether they come from Bondi’s original paper or the later work. Alternatively, atomic radii from the Universal Force Field (UFF) are available.766 The main appeal of UFF radii is that they are defined for all atoms of the periodic table, though the quality of these radii for PCM applications is unclear. Finally, the user may specify his or her own radii for cavity construction using a $van_der_waals input section, the format for which is described in Section 12.2.9. No scaling factor is applied to user-defined radii. Note that is allowed for a particular atomic radius, in which case the atom in question is not used to construct the cavity surface. This feature facilitates the construction of “united atom” cavities,55 in which the hydrogen atoms do not get their own spheres and the heavy-atom radii are increased to compensate Finally, since the solvent molecules should not be able to penetrate all the way to the atomic vdW radii of the solute, it is traditional either to scale the atomic radii (vdW surface construction) or else to augment them with an assumed radius of a spherical solvent molecule (SAS construction), but not both.
Historically, discretization of the cavity surface has involved “tessellation” methods that divide the cavity surface area into finite polygonal “tesserae”. (The GEPOL algorithm698 is perhaps the most widely-used tessellation scheme.) Tessellation methods, however, suffer not only from discontinuities in the cavity surface area and solvation energy as a function of the nuclear coordinates, but in addition they lead to analytic energy gradients that are complicated to derive and implement. To avoid these problems, Q-Chem’s SWIG PCM implementation519, 520, 521 uses Lebedev grids to discretize the atomic spheres. These are atom-centered grids with icosahedral symmetry, and may consist of anywhere from 26 to 5294 grid points per atomic sphere. The default values used by Q-Chem were selected based on extensive numerical tests,520, 521 and they are looser for MM atoms (in MM/PCM or QM/MM/PCM jobs) than they are for QM atoms, reflecting the more complicated electrostatic potential that is generated by a QM density as compared to an MM point charge. For QM atoms, the default is to use points for hydrogen atoms and points for all other atoms, whereas for MM atoms the default is for hydrogen and for non-hydrogen. These default values exhibit good rotational invariance ( kcal/mol differences in when the molecule is rotated520) and absolute solvation energies that typically lie within kcal/mol of the limit,520, 521 at least for charge-neutral solutes.
Note that earlier versions of Q-Chem used denser grids by default. (In versions up to and including Q-Chem v. 4.2, the default was for all QM atoms, but was switched to beginning with v. 4.2.1. The defaults mentioned above are the current ones starting with v. 5.2.) However, grid errors of 0.5 kcal/mol are well within the intrinsic accuracy of for these models. Grid errors in solvatochromatic shifts for excitation energies tend to be 0.01–0.02 eV, which is well within the intrinsic accuracy of nearly any excited-state methodology. If questions about grid accuracy arise, we suggest using as a high-quality option and as an essentially converged option.520, 521
Note: The acceptable values for the number of Lebedev points per sphere are , 50, 110, 194, 302, 434, 590, 770, 974, 1202, 1454, 1730, 2030, 2354, 2702, 3074, 3470, 3890, 4334, 4802, or 5294.
Especially for complicated molecules, the user may want to visualize the cavity surface. This can be accomplished by setting PrintLevel , which will trigger the generation of several .PQR files that describe the cavity surface. (These are written to the Q-Chem output file.) The .PQR format is similar to the common .PDB (Protein Data Bank) format, but also contains charge and radius information for each atom. One of the output .PQR files contains the charges computed in the PCM calculation and radii (in Å) that are half of the square root of the surface area represented by each surface grid point. Thus, in examining this representation of the surface, larger discretization points are associated with larger surface areas. A second .PQR file contains the solute’s electrostatic potential (in atomic units), in place of the charge information, and uses uniform radii for the grid points. These .PQR files can be visualized using various third-party software, including the freely-available Visual Molecular Dynamics (VMD) program,W. Humphrey, A. Dalke, and K. Schulten (1996), 17 which is particularly useful for coloring the .PQR surface grid points according to their charge, and sizing them according to their contribution to the molecular surface area. (Examples of such visualizations can be found in Ref. 519.)
Finally, note that setting Method to Spherical in the $pcm input selection requests the construction of a solute cavity consisting of a single, fixed sphere. This is generally not recommended but is occasionally useful for making contact with the results of Born models in the literature, or the Kirkwood-Onsager model discussed in Section 12.2.1. In this case, the cavity radius and its center must also be specified in the $pcm section. The keyword HeavyPoints controls the number of Lebedev grid points used to discretize the surface.
The following example shows a very basic PCM job. The solvent dielectric is specified in the $solvent section, which is described below.
$molecule 0 1 C -0.245826 -0.351674 -0.019873 C 0.244003 0.376569 1.241371 O 0.862012 -0.527016 2.143243 F 0.776783 -0.909300 -0.666009 F -0.858739 0.511576 -0.827287 F -1.108290 -1.303001 0.339419 H -0.587975 0.878499 1.736246 H 0.963047 1.147195 0.961639 H 0.191283 -1.098089 2.489052 $end $rem JOBTYPE OPT BASIS 6-31G* METHOD B3LYP SOLVENT_METHOD PCM $end $pcm Theory CPCM Method SWIG Solver Inversion HeavyPoints 194 HPoints 194 Radii Bondi vdwScale 1.2 $end $solvent Dielectric 78.39 $end
$molecule 0 1 O 0.00000000 0.00000000 0.11722303 H -0.75908339 0.00000000 -0.46889211 H 0.75908339 0.00000000 -0.46889211 $end $rem METHOD HF BASIS 6-31g** SOLVENT_METHOD pcm $end $pcm method spherical ! single spherical cavity with 590 discretization points HeavyPoints 590 CavityRadius 1.8 ! Solute Radius, in Angstrom CavityCenter 0.0 0.0 0.0 ! Will be at center of Standard Nuclear Orientation Theory SSVPE $end $solvent Dielectric 35.9 ! Acetonitrile $end
Finally, we consider an example of a united-atom cavity. Note that a user-defined van der Waals radius is supplied only for carbon, so the hydrogen radius is taken to be zero and thus the hydrogen atoms are not used to construct the cavity surface. (As mentioned above, the format for the $van_der_waals input section is discussion in Section 12.2.9).
$comment Benzene (in benzene), with a united-atom cavity construction R = 2.28 A for carbon, R = 0 for hydrogen $end $molecule 0 1 C 1.38620 0.000000 0.000000 C 0.69310 1.200484 0.000000 C -0.69310 1.200484 0.000000 C -1.38620 0.000000 0.000000 C -0.69310 -1.200484 0.000000 C 0.69310 -1.200484 0.000000 H 2.46180 0.000000 0.000000 H 1.23090 2.131981 0.000000 H -1.23090 2.131981 0.000000 H -2.46180 0.000000 0.000000 H -1.23090 -2.131981 0.000000 H 1.23090 -2.131981 0.000000 $end $rem EXCHANGE hf BASIS 6-31G* SOLVENT_METHOD pcm $end $pcm theory iefpcm ! this is a synonym for ssvpe method swig printlevel 1 radii read $end $solvent dielectric 2.27 $end $van_der_waals 1 6 2.28 1 0.00 $end
The solvent for PCM calculations is specified using the $solvent section, as documented below. In addition, the $solvent section can be used to incorporate non-electrostatic interaction terms into the solvation energy. (The Theory keyword in the $pcm section specifies only how the electrostatic interactions are handled.) The general form of the $solvent input section is shown below. The $solvent section was used above to specify parameters for the Kirkwood-Onsager SCRF model, and will be used again below to specify the solvent for SM calculations (Section 12.2.8); in each case, the particular options that can be listed in the $solvent section depend upon the value of the $rem variable SOLVENT_METHOD.
$solvent NonEls <Option> NSolventAtoms <Number unique of solvent atoms> SolventAtom <Number1> <Number2> <Number3> <SASradius> SolventAtom <Number1> <Number2> <Number3> <SASradius> . . . <Keyword> <parameter/option> . . . $end
The keyword SolventAtom requires multiple parameters, whereas all other keywords require only a single parameter. In addition to any (optional) non-electrostatic parameters, the $solvent section is also used to specify the solvent’s dielectric constant. If non-electrostatic interactions are ignored, then this is the only keyword that is necessary in the $solvent section. For nonequilibrium TDDFT/C-PCM calculations (Section 220.127.116.11), the optical dielectric constant should be specified in the $solvent section as well.
The non-electrostatic interactions currently available in Q-Chem are based on the work of Cossi et al.,193 and are computed outside of the SCF procedure used to determine the electrostatic interactions. The non-electrostatic energy is highly dependent on the input parameters and can be extremely sensitive to the radii chosen to define the solute cavity. Accordingly, the inclusion of non-electrostatic interactions is highly empirical and probably needs to be considered on a case-by-case basis. Following Ref. 193, the cavitation energy is computed using the same solute cavity that is used to compute the electrostatic energy, whereas the dispersion/repulsion energy is computed using a solvent-accessible surface.
The following keywords (in the $solvent section) are used to define non-electrostatic parameters for PCM calculations.
The following example illustrates the use of the non-electrostatic interactions.
$molecule 0 1 C -0.245826 -0.351674 -0.019873 23 C 0.244003 0.376569 1.241371 22 O 0.862012 -0.527016 2.143243 24 F 0.776783 -0.909300 -0.666009 26 F -0.858739 0.511576 -0.827287 26 F -1.108290 -1.303001 0.339419 26 H -0.587975 0.878499 1.736246 27 H 0.963047 1.147195 0.961639 27 H 0.191283 -1.098089 2.489052 25 $end $rem JOBTYPE OPT BASIS 6-31G* METHOD B3LYP SOLVENT_METHOD PCM FORCE_FIELD OPLSAA $end $pcm Theory CPCM Method SWIG Solver Inversion HeavyPoints 194 HPoints 194 Radii Bondi vdwScale 1.2 $end $solvent NonEls LJCav NSolventAtoms 2 SolventAtom 8 1 186 1.30 SolventAtom 1 2 187 0.01 SolventRadius 1.35 Temperature 298.15 Pressure 1.0 SolventRho 0.03333 Dielectric 78.39 $end
The OpticalDielectric keyword in $solvent is always needed. The LR energy is automatically calculated while the CIS/TDDFT calculations are performed with PCM, but it is turned off while the perturbation scheme is employed.
$molecule 0 1 C 0 0 0.0 O 0 0 1.21 $end $rem EXCHANGE B3LYP CIS_N_ROOTS 10 CIS_SINGLETS TRUE CIS_TRIPLETS TRUE RPA TRUE BASIS 6-31+G* SOLVENT_METHOD PCM $end $pcm Theory CPCM Method SWIG Solver Inversion Radii Bondi $end $solvent Dielectric 78.39 OpticalDielectric 1.777849 $end
$molecule 0 1 C 0.000046 -0.000398 1.904953 C 1.210027 0.000379 1.186051 C 1.214640 -0.000065 -0.194515 C 0.000164 -0.000616 -0.933832 C -1.214349 -0.001557 -0.194687 C -1.209753 -0.001846 1.185775 H 2.151949 0.001377 1.722018 H 2.164371 0.000481 -0.709640 H -2.164082 -0.002008 -0.709781 H -2.151763 -0.002287 1.721615 C -0.000227 0.001061 3.325302 N -0.000475 0.002405 4.484321 N 0.000053 -0.000156 -2.297372 C -1.258656 0.001284 -3.036994 H -1.041042 0.001615 -4.102376 H -1.860897 -0.885647 -2.811117 H -1.859247 0.889133 -2.810237 C 1.258563 -0.000660 -3.037285 H 1.860651 0.886208 -2.810755 H 1.859362 -0.888604 -2.811461 H 1.040664 -0.000097 -4.102609 $end $rem EXCHANGE LRC-wPBEPBE OMEGA 260 BASIS 6-31G* CIS_N_ROOTS 10 RPA 2 CIS_SINGLETS 1 CIS_TRIPLETS 0 CIS_DYNAMIC_MEM TRUE CIS_RELAXED_DENSITY TRUE USE_NEW_FUNCTIONAL TRUE SOLVENT_METHOD PCM PCM_PRINT 1 $end $pcm NonEquilibrium Theory IEFPCM StateSpecific Perturb $end $solvent Dielectric 35.688000 ! Acetonitrile OpticalDielectric 1.806874 $end
$molecule 0 1 C -0.189057 -1.215927 -0.000922 H -0.709319 -2.157526 -0.001587 C 1.194584 -1.155381 -0.000067 H 1.762373 -2.070036 -0.000230 C 1.848872 0.069673 0.000936 H 2.923593 0.111621 0.001593 C 1.103041 1.238842 0.001235 H 1.595604 2.196052 0.002078 C -0.283047 1.185547 0.000344 H -0.862269 2.095160 0.000376 C -0.929565 -0.042566 -0.000765 O -2.287040 -0.159171 -0.001759 H -2.663814 0.725029 0.001075 $end $rem EXCHANGE wPBE CORRELATION PBE LRC_DFT 1 OMEGA 300 BASIS 6-31+G* SCF_CONVERGENCE 8 SOLVENT_METHOD PCM PCM_PRINT 1 $end $pcm NonEquilibrium $end $solvent Dielectric 78.39 OpticalDielectric 1.777849 $end @@@ $molecule 1 2 read $end $rem EXCHANGE wPBE CORRELATION PBE LRC_DFT 1 OMEGA 300 BASIS 6-31+G* SCF_CONVERGENCE 8 SOLVENT_METHOD PCM PCM_PRINT 1 SCF_GUESS READ $end $pcm StateSpecific Marcus $end $solvent Dielectric 78.39 OpticalDielectric 1.777849 $end
$molecule 0 1 C -0.000002 -0.000003 1.894447 C 1.230178 0.000003 1.160765 C 1.244862 0.000001 -0.200360 C 0.000003 -0.000013 -0.917901 C -1.244845 -0.000000 -0.200354 C -1.230160 0.000000 1.160753 H 2.171128 0.000013 1.704064 H 2.180670 0.000009 -0.750140 H -2.180657 0.000006 -0.750120 H -2.171112 0.000006 1.704037 C -0.000023 -0.000002 3.293509 N -0.000011 -0.000001 4.459415 N -0.000008 -0.000005 -2.309193 C 0.000001 1.253881 -3.021235 H 0.000005 1.099076 -4.098660 H -0.882908 1.818975 -2.705815 H 0.882911 1.818967 -2.705805 C -0.000009 -1.253886 -3.021246 H 0.882917 -1.818966 -2.705849 H -0.882901 -1.818992 -2.705803 H -0.000038 -1.099071 -4.098670 $end $rem EXCHANGE LRC-wPBEPBE OMEGA 280 BASIS 6-31G* CIS_N_ROOTS 10 RPA 2 CIS_SINGLETS 1 CIS_TRIPLETS 0 CIS_DYNAMIC_MEM TRUE CIS_RELAXED_DENSITY TRUE USE_NEW_FUNCTIONAL TRUE SOLVENT_METHOD PCM PCM_PRINT 1 $end $pcm ChargeSeparation Excited StateSpecific 1 $end $solvent Dielectric 35.688000 ! Acetonitrile OpticalDielectric 1.806874 $end @@@ $molecule READ $end $rem EXCHANGE LRC-wPBEPBE OMEGA 280 BASIS 6-31G* SCF_GUESS READ SCF_CONVERGENCE 8 SOLVENT_METHOD PCM PCM_PRINT 1 $end $pcm StateSpecific Marcus $end $solvent Dielectric 35.688000 ! Acetonitrile OpticalDielectric 1.806874 $end