Skip to content

Polarizable Continuum Model (PCM)

TeraChem has a GPU-accelerated implementation of the Polarizable Continuum Model1 that runs roughly an order of magnitude faster than typical CPU implementations. All implemented variants belong to the conductor-like solvation family (COSMO):

  • C-PCM2 (also known as G-COSMO3) with smooth first-order energy derivatives, using the ISWIG4 or SWIG5 cavity discretization. This is the recommended choice for general-purpose work.
  • The original COSMO by Klamt,6 using a polyhedron grid. Useful for inspection but not recommended for production because its gradients are not smooth.

Both ground-state SCF and excited-state TDDFT/CIS calculations can be coupled to a PCM solvent.

Quick start: ground-state PCM

The minimal PCM input requires the pcm and epsilon keywords:

pcm_ethylene.in
basis        6-31g
method       rhf
coordinates  c2h4.xyz
charge       0
spinmult     1
maxit        100
gpus         1
run          energy

pcm          cosmo
epsilon      78.39
end

pcm cosmo activates PCM, and epsilon is the solvent dielectric (water \(\approx 78.39\) at 298 K). All other PCM parameters take their defaults, so this is a C-PCM ISWIG calculation with Lebedev grid level 17 (~110 points/atom). The run type can be any of energy, gradient, minimize, or md.

For an SCF/TDDFT job at the camB3LYP level with explicit grid level and Bondi radii:

pcm_cytosine.in
basis         6-31+g*
method        camb3lyp
coordinates   cytosin.1.xyz
charge        0
spinmult      1
scf           diis+a
maxit         100
xtol          1e-4
threall       1e-14
convthre      1.0e-8
dftgrid       1
gpus          all
run           energy

pcm             cosmo
pcm_grid        iswig
pcmgrid_heavy   17
pcmgrid_h       17
epsilon         78.39
pcm_radii       bondi
pcm_rad_scale   1.1
solvent_radius  0

# Conjugate-gradient solver
dynamiccg       3
cgprecond       blockjacobi3
cgblocksize     100
end

Cavity choice (pcm_grid)

The pcm_grid keyword is the most important choice — it selects the COSMO variant:

Value Recommended use
iswig (default) Recommended for production. Smooth analytical first derivatives.
swig Alternative smooth ISWIG-family discretization.
lebedev Lebedev grid without switching functions or Gaussian polarization charges. Useful for understanding why ISWIG/SWIG are preferred; not recommended for production.
polyhedron Original Klamt-style polyhedron grid. Gradients are not smooth — diagnostic use only.
sphere Single spherical cavity enclosing the whole system; the cavity center is fixed and does not follow the COM during MD. Only available when spherical boundary conditions are turned on (mdbc spherical).

For ISWIG/SWIG/Lebedev cavities, the surface point density is controlled per atom type by pcmgrid_h (hydrogens) and pcmgrid_heavy (heavy atoms). Both take Lebedev order indices in \([1, 131]\).

Orders 12, 13 and 24–27 are not usable for PCM

The 74-, 230- and 266-point Lebedev–Laikov rules — reached by pcmgrid_* orders 12/13, 24/25 and 26/27 respectively — contain negative weights. PCM reuses the Lebedev weights as segment areas and takes their square root, both for the SWIG Gaussian width \(\zeta_i = g/\sqrt{4\pi w_i}\,/R_A\) and for the COSMO self-interaction diagonal \(A_{ii} = a_\mathrm{diag}/\sqrt{\text{area}_i}\). A negative weight gives \(\sqrt{\text{negative}} = \mathrm{NaN}\), which propagates into the Fock matrix and produces a NaN energy rather than a diagnosable failure.

Use a neighboring order — 11, 17, 23 or 29 — instead. Depending on your version, these orders either fail silently in this way or are rejected with an explicit error at startup.

For the polyhedron grid use nspa and nppa instead — pcmgrid_h and pcmgrid_heavy are not honored.

For the sphere grid, pcm_global_radius sets the radius of the enclosing sphere.

It is often useful to inspect the cavity directly. Setting print_ms 1 writes the discretized surface to a file that can be loaded in a visualizer.

Atomic radii

The Van der Waals radii used to construct the cavity are selected with pcm_radii:

Value Radii set
bondi (default) Bondi radii
klamt Original Klamt radii
read Custom radii from an external file (see below)

The radii are uniformly scaled by pcm_rad_scale (default 1.2); typical values are 1.0-1.2 for ISWIG/SWIG cavities. The solvent probe radius is solvent_radius (Å); a value of 0 (the default for many cases) gives no probe sphere.

Custom atomic radii file

For elements where the built-in tables are unsuitable, supply your own radii:

pcm_radii       read
pcm_radii_file  /absolute/path/to/radii.dat

The path must be absolute. The file lists one element per line with the atomic number in the first column and the radius in Å in the second column, separated by whitespace:

radii.dat
1  1.301
8  1.7221

This file is sufficient for a system containing only H and O. The file may contain more elements than the system uses; extras are ignored. However, if any element present in the geometry is missing from the radii file, TeraChem will abort.

Conjugate-gradient solver tuning

The PCM polarization charges are obtained by solving a linear system with a preconditioned conjugate-gradient (CG) solver. The defaults are usually fine, but for large systems performance can be substantially improved by tuning the preconditioner:

Keyword Description
dynamiccg Enables an adaptive CG tolerance schedule (recommended)
cgtoltight Tight convergence threshold used at the final SCF iteration
cgtolscale Number of decimal places by which to loosen the tolerance during the SCF
dynamiccgtol Custom tolerance for the dynamic schedule
cgprecond Preconditioner: jacobi or blockjacobi (alias blockjacobi3). Variants blockjacobi1, blockjacobi2, blockjacobi4, blockjacobi5 are for performance benchmarking only.
cgblocksize Block size for the block-Jacobi family of preconditioners
pcm_matrix Matrix storage mode: full, packed, no, no_sp

For systems beyond ~1000 atoms, cgprecond blockjacobi3 with cgblocksize of order 100 typically gives the best performance.

Excited-state PCM solvation

For CIS or TDDFT excited-state calculations in solvent, TeraChem supports both linear-response and state-specific PCM solvation.

Linear-response PCM (lr_pcm_solvation)

LR-PCM treats the solvent reaction field self-consistently with the linear response equations. Values:

Value Description
eq Equilibrium solvation — the solvent is fully relaxed for the excited state. Appropriate for slow processes such as fluorescence.
neq Non-equilibrium solvation — only the fast electronic component of the solvent responds. Appropriate for fast processes such as vertical absorption.

The non-equilibrium variant requires a second dielectric constant for the optical-frequency response, given by fast_epsilon.

lr_pcm_cis_grad.in
basis           6-31g
method          rhf
coordinates     ground_pcm.xyz
charge          0
spinmult        1
maxit           100
threall         1e-13
convthre        1.0e-7
dftgrid         5

cis             yes
cisnumstates    1
cismaxiter      100
cistarget       1
cisconvtol      1e-6

gpus            all
run             gradient

pcm                 cosmo
pcm_grid            iswig
lr_pcm_solvation    eq
pcmgrid_heavy       17
pcmgrid_h           17
epsilon             78.39
fast_epsilon        78.39
pcm_radii           bondi
pcm_rad_scale       1.2
solvent_radius      0
dynamiccg           2
cgprecond           blockjacobi
cgblocksize         100
end

State-specific PCM (ss_pcm_solvation)

State-specific PCM treats the solvent response self-consistently with the target excited state itself:

Value Description
eq Equilibrium SS-PCM
ground_neq Ground-state non-equilibrium SS-PCM
neq Non-equilibrium SS-PCM

SS-PCM is iterative; convergence is controlled by sspcm_maxit and sspcm_convthre. Like LR-PCM, the non-equilibrium variants need fast_epsilon.

Extreme-Pressure PCM (XP-PCM)

The extreme-pressure PCM7 is an implicit-pressure model for simulating high-pressure conditions on the 1-300 GPa scale. TeraChem implements8 the Cammi-style formulation. Currently only single-point XP-PCM energies are supported — analytical gradients are not yet available.

Required keywords

XP-PCM is activated by pcm xppcm. The grid must be iswig (the only supported choice). Three keywords describe the solvent at the molecular level:

Keyword Meaning
xppcm_nb Solvent valence-electron number
xppcm_mb Solvent molecular weight
xppcm_rhob Solvent number density
xppcm_eta Polarizability scaling parameter (advanced)
xppcm_xi Pauli-repulsion scaling parameter (advanced)

XP-PCM cavities need a denser grid than ground-state PCM — Lebedev level 35 or higher is recommended for pcmgrid_heavy and pcmgrid_h.

To reflect compression of the molecular volume at high pressure, the radii scaling factor pcm_rad_scale is shrunk below the usual value (typically not above 1.2).

Example: argon at 1.1 × Bondi radius

xppcm_ar.in
method         b3lyp
basis          aug-cc-pvdz
coordinates    Ar.xyz
charge         0
spinmult       1
scf            diis
maxit          100
convthre       1e-6
gpus           1
run            energy

pcm            xppcm
pcm_grid       iswig
pcmgrid_heavy  35
pcmgrid_h      35
epsilon        2.0165
pcm_rad_scale  1.1
xppcm_nb       36
xppcm_mb       84.16
xppcm_rhob     0.779
jobname        energy_f1.1_rhob0.779_g35
end

Getting pressure from XP-PCM

An XP-PCM single point gives the free energy of the solute at one cavity volume, not the pressure. To extract pressure, run a series of XP-PCM jobs at several values of pcm_rad_scale (cavity scale factors) and fit the free-energy-vs-volume curve to an equation of state.9

Summary of keywords

Master switches and cavity geometry

Keyword Type Default Description
pcm string not set cosmo (C-PCM/COSMO) or xppcm (XP-PCM). Master switch.
epsilon float 80.0 Solvent dielectric constant
pcm_grid string iswig Cavity discretization: iswig, swig, lebedev, polyhedron, or sphere
pcmgrid_heavy int 17 Lebedev grid order for heavy atoms (used with iswig/swig/lebedev)
pcmgrid_h int 17 Lebedev grid order for hydrogens
pcmgrid_heavy_mm int — Lebedev grid order for heavy atoms in the MM region of a QM/MM calculation
pcmgrid_h_mm int — Lebedev grid order for hydrogens in the MM region
pcmgrid_ptchrg int — Lebedev grid order for point-charge atoms
pcm_global_radius float — Radius of the enclosing sphere when pcm_grid sphere
nspa int — Polyhedron grid: surface points per atom (only with pcm_grid polyhedron)
nppa int — Polyhedron grid: points per Angstrom² (only with pcm_grid polyhedron)
cosmodsc float — Polyhedron grid discretization parameter (ignored by ISWIG/SWIG)

Atomic radii

Keyword Type Default Description
pcm_radii string bondi Radii set: bondi, klamt, or read
pcm_radii_file filename — Absolute path to a custom radii file (used with pcm_radii read)
pcm_rad_scale float 1.2 Uniform scaling of the radii
solvent_radius float 0.0 Probe sphere radius in Å
pcm_ptchrg_radius float — Cavity radius for point charges in QM/MM PCM jobs

One-electron integral options

Keyword Type Default Description
pcm_pair_driven bool no Use pair-driven evaluation of PCM one-electron integrals
pcm_threoe float — Threshold for one-electron integrals in PCM

CG solver

Keyword Type Default Description
pcm_matrix string full Matrix storage: full, packed, no, no_sp
dynamiccg int 0 Adaptive CG tolerance schedule (recommended values: 2 or 3)
dynamiccgtol float — Custom dynamic-tolerance value
cgtoltight float — Final-step CG convergence threshold
cgtolscale int — Loosen the CG threshold by this many decimal places during SCF
cgprecond string jacobi jacobi or blockjacobi (default is blockjacobi3)
cgblocksize int — Block size for block-Jacobi preconditioners
cg_random_block string yes Use a randomized block partition (no to disable)

Excited-state solvation

Keyword Type Default Description
lr_pcm_solvation string not set eq or neq — linear-response PCM for CIS/TDDFT
ss_pcm_solvation string not set eq, ground_neq, or neq — state-specific PCM for CIS/TDDFT
sspcm_maxit int — Maximum SS-PCM iterations
sspcm_convthre float — SS-PCM convergence threshold
fast_epsilon float — Optical (fast) dielectric constant for non-equilibrium solvation

XP-PCM extension

Keyword Type Default Description
xppcm_nb float — Solvent valence-electron number
xppcm_mb float — Solvent molecular weight (g/mol)
xppcm_rhob float — Solvent number density
xppcm_eta float — Polarizability scaling parameter
xppcm_xi float — Pauli-repulsion scaling parameter

I/O and diagnostics

Keyword Type Default Description
print_ms int 0 If nonzero, write the discretized molecular surface for visualization
pcm_print string yes verbose, yes, or no
pcm_scale int 0 Internal cavity-scaling flag (0 or 1)
pcm_write filename — Save the PCM reaction-field potential to this file
pcm_read filename — Read the PCM reaction-field potential from this file
split_pcm_energy int 0 Print the PCM contribution to the energy as a separate term

  1. F. Liu, N. Luehr, J. K. Kulik, and T. J. Martinez, J. Chem. Theory Comput. 11, 3131-3144 (2015). ↩

  2. V. Barone and M. Cossi, J. Phys. Chem. A 102, 1995-2001 (1998). ↩

  3. T. N. Truong and E. V. Stefanovich, Chem. Phys. Lett. 240, 253-260 (1995). ↩

  4. A. W. Lange and J. M. Herbert, J. Chem. Phys. 133, 244111 (2010). ↩

  5. D. M. York and M. Karplus, J. Phys. Chem. A 103, 11060-11079 (1999). ↩

  6. A. Klamt and G. Schuurmann, J. Chem. Soc. Perkin Trans. 2, 799-805 (1993). ↩

  7. R. Cammi, V. Verdolino, B. Mennucci, and J. Tomasi, Chem. Phys. 344, 135 (2008). ↩

  8. A. Gale, E. Hruska, and F. Liu, J. Chem. Phys. 154, 244103 (2021). ↩

  9. R. Cammi, J. Comp. Chem. 36, 2246-2259 (2015). ↩