Welcome to PCMSolver’s documentation!

This is the documentation for the PCMSolver application programming interface. PCMSolver is an API for solving the Polarizable Continuum Model electrostatic problem [TMC05]

_images/pcmsolver-scheme.png

With PCMSolver we aim to:

  1. provide a plug-and-play library for adding the PCM functionality to any quantum chemistry program;
  2. create a playground for easily extending the implementation of the model.

PCMSolver is distributed under the terms of the GNU Lesser General Public License. An archive with the currently released source can be found on GitHub or on zenodo.

@misc{PCMSolver,
  note = {\texttt{PCMSolver}, v1.1.9 an Application Programming Interface for the
	Polarizable Continuum Model electrostatic problem, written by R.~Di~Remigio, L.~Frediani and K.~Mozgawa
        with contributions from R.~Bast, J.~Juselius and V.~Weijo
        (see http://pcmsolver.readthedocs.io/)"
}

Contents:

PCMSolver Users’ Manual

Building the module

PCMSolver configuration and build process is managed through CMake.

Prerequisites and dependencies

A number of prerequisites and dependencies are to be satisfied to successfully build the module. It will be here assumed that you want to perform a “full” build, i.e. you want to build the static libraries to be linked to your QM program, the unit test suite and an offline copy of this documentation.

Compilers
  • a C++ compiler, compliant with the 1998 ISO C++ standard plus the 2003 technical corrigendum and some additional defect reports.
  • a C compiler, compliant with the ISO C99 standard.
  • a Fortran compiler, compliant with the Fortran 2003 standard.

The list of primary test environments can be found in the README.md file. It is entirely possible that using other compiler versions you might be able to build the module. In order to ensure that you have a sane build, you will have to run the unit test suite.

Libraries and toolchain programs
  • CMake version 2.8.9 and higher;
  • Git version 1.7.1 and higher;
  • Python interpreter 2.4 and higher;
  • Boost libraries version 1.54.0 and higher;

Note

Version 1.54.0 of Boost libraries is shipped with the module and resides in the cmake/downloaded subdirectory. Unless you want to use another version of Boost, you should not worry about satisfying this dependency.

  • zlib version 1.2 and higher (unit test suite only);
  • Doxygen version 1.7.6 and higher (documentation only)
  • Perl (documentation only)
  • PyYAML Python module (documentation only)
  • Sphinx (documentation only)

PCMSolver relies on the Eigen template libraries version 3.3.0 and higher. Version 3.3.0 of Eigen libraries is shipped with the module and resides in the external subdirectory.

Configuration

Configuration is managed through the front-end script setup residing in the repository main directory. Issuing:

./setup [options] [build path]

will create the build directory in build path and run CMake with the given options. By default, files are configured in the build directory. The -h or –help option will list the available options and their effect. Options can be forwarded directly to CMake by using the –cmake-options flag and listing the -D... options. Usually the following command is sufficient to get the configuration done for a debug build, including compilation of the unit test suite:

./setup --type=debug

The unit tests suite is always compiled in standalone mode, unless the -DENABLE_TESTS=OFF option is forwarded to CMake.

Getting Boost

You can get Boost libraries in two ways:

  • already packaged by your Linux distribution or through MacPorts/Brew;
  • by downloading the archive from http://www.boost.org/ and building it yourself.

In case your distribution packages a version older than 1.54.0 you might chose to either build Boost on your own or to rely on the automated build of the necessary Boost libraries when compiling the module (recommended). Full documentation on how to build Boost on Unix variants is available here. It is here assumed that the user does not have root access to the machine and will install the libraries to a local prefix, a subdirectory of /home/user-name tipically. Once you’ve downloaded and unpacked the archive, run the bootstrap script to configure:

cd path/to/boost
./bootstrap.sh --prefix=/home/user-name/boost

Running ./bootstrap.sh –help will list the available options for the script. To build run:

./b2 install

This might take a while. After a successful build you will find the headers in /home/user-name/boost/include and libraries in /home/user-name/boost/lib Now, you will have Boost in a nonstandard location. Without hints CMake will not be able to find it and configuration of PCMSolver will fail. To avoid this, you will have to pass the location of the headers and libraries to the setup script, either with:

./setup --boost-headers=/home/user-name/boost/include --boost-libs=/home/user-name/boost/lib

or with:

./setup -DBOOST_INCLUDEDIR=/home/user-name/boost/include -DBOOST_LIBRARYDIR=/home/user-name/boost/lib
Advanced configuration options

These options are marked as advanced as it is highly unlikely they will be useful when not programming the library:

  • –exdiag Enable C++ extended diagnostics flags. Disabled by default.
  • –ccache Enable use of ccache for C/C++ compilation caching. Enabled by default, unless ccache is not available.
  • –build-boost Deactivate Boost detection and build on-the-fly. Disabled by default.
  • –eigen Root directory for Eigen3. Search for Eigen3 in the location provided by the user. If search fails, fall back to the version bundled with the library.
  • –static Create only static library. Disabled by default.

Some options can only be tweaked via –cmake-options to the setup script:

  • ENABLE_CXX11_SUPPORT Enable C++11 support. Tries to detect which C++11 features are supported by the compiler and enables use of the new standard. Enabled by default.

    Warning

    This option is always overridden for some compilers that have buggy C++11 support.

  • ENABLE_DOCS Enable build of documentation. This requires a number of additional dependencies. If any of these are not met, documentation is not built. Enabled by default.

  • ENABLE_LOGGER Enable compilation of logger sources. Disabled by default.

    Warning

    The logger is not currently in use in any part of the code.

  • ENABLE_TIMER Enable compilation of timer sources. Enabled by default.

  • BUILD_STANDALONE Enable compilation of standalone run_pcm executable. Enabled by default.

  • ENABLE_FORTRAN_API Enable compilation of the Fortran90 bindings for the API. Disabled by default.

  • ENABLE_GENERIC Enable mostly static linking in shared library. Disabled by default.

  • ENABLE_TESTS Enable compilation of unit tests suite. Enabled by default.

  • SHARED_LIBRARY_ONLY Create only shared library. Opposite of –static.

  • PYMOD_INSTALL_LIBDIR If set, installs python scripts/modules to ${CMAKE_INSTALL_LIBDIR}${PYMOD_INSTALL_LIBDIR}/pcmsolver rather than the default ${CMAKE_INSTALL_BINDIR} (i.e., bin).

  • CMAKE_INSTALL_BINDIR Where to install executables, if not to bin.

  • CMAKE_INSTALL_LIBDIR Where to install executables, if not to bin.

  • CMAKE_INSTALL_INCLUDESDIR Where to install executables, if not to bin.

  • CMAKE_INSTALL_BINDIR Location within CMAKE_INSTALL_PREFIX (--prefix) to which executables are installed (default: bin).

  • CMAKE_INSTALL_LIBDIR Location within CMAKE_INSTALL_PREFIX (--prefix) to which libraries are installed (default: lib).

  • CMAKE_INSTALL_INCLUDEDIR Location within CMAKE_INSTALL_PREFIX (--prefix) to which headers are installed (default: include).

  • PYMOD_INSTALL_LIBDIR If set, location within CMAKE_INSTALL_LIBDIR to which python modules are installed, ${CMAKE_INSTALL_LIBDIR}${PYMOD_INSTALL_LIBDIR}/pcmsolver. If not set, python modules installed to default ${CMAKE_INSTALL_BINDIR} (i.e., bin).

Build and test

To compile and link, just go to the build directory and run:

make -j N

where N is the number of cores you want to use when building.

Note

Building on more than one core can sometimes result in a “race condition” and a crash. If that happens, please report the problem as an issue on our issue tracker on GitHub. Running make on a single core might get you through compilation.

To run the whole test suite:

ctest -j N

You can also use CTest to run a specific test or a set of tests. For example:

ctest -R gepol

will run all the test containing the string “gepol” in their name.

If Doxygen was found, an offline copy of this documentation can be built by:

make doc

and visualized by opening the doc/html/index.html file in your browser.

Input description

PCMSolver needs a number of input parameters at runtime. The API provides two ways of providing them:

  1. by means of an additional input file, parsed by the pcmsolver.py script;
  2. by means of a special section in the host program input.

Method 1 is more flexible: all parameters that can be modified by the user are available. The host program needs only copy the additional input file to the scratch directory before execution. Method 2 just gives access to the core parameters.

In this page, input style and input parameters available in Method 1 will be documented.

Input style

The input for PCMSolver is parsed through the Getkw library written by Jonas Juselius and is organized in sections and keywords. Input reading is case-insensitive. An example input structure is shown below, there are also some working examples in the directory examples. A general input parameter has the following form (Keyword = [Data type]):

Units = [String]
CODATA = [Integer]
Cavity {
       Type = [String]
       NpzFile = [String]
       Area = [Double]
       Scaling = [Bool]
       RadiiSet = [String]
       MinRadius = [Double]
       Mode = [String]
       Atoms = [Array of Integers]
       Radii = [Array of Doubles]
       Spheres = [Array of Doubles]
}
Medium {
       Solvent = [String]
       SolverType = [String]
       MatrixSymm = [Bool]
       Correction = [Double]
       ProbeRadius = [Double]
       Green<GreenTag> {
             Type = [String]
             Der = [String]
             Eps = [Double]
             EpsDyn = [Double]
             Eps1 = [Double]
             EpsDyn1 = [Double]
             Eps2 = [Double]
             EpsDyn2 = [Double]
             Center = [Double]
             Width = [Double]
             InterfaceOrigin = [Array of Doubles]
             MaxL = [Integer]
       }
}

Array-valued keywords will expect the array to be given in comma-separated format and enclosed in square brackets. The purpose of tags is to distinguish between cases in which multiple instances of the same kind of object can be managed by the program. There exist only certain legal tagnames and these are determined in the C++ code. Be aware that the input parsing script does not check the correctness of tags.

Input parameters

Available sections:

  • top section: sets up parameters affecting the module globally;
  • Cavity: sets up all information needed to form the cavity and discretize its surface;
  • Medium: sets up the solver to be used and the properties of the medium, i.e. the Green’s functions inside and outside the cavity;
  • Green, subsection of medium. Sets up the Green’s function inside and outside the cavity.

Warning

Exactly matching results obtained from implementations of IEFPCM and/or CPCM (COSMO) given in other program packages requires careful selection of all the parameters involved. A partial checklist of parameters you should always keep in mind:

  • solvent permittivities (static and optical)
  • atomic radii set
  • scaling of the atomic radii
  • cavity surface
  • cavity partition (tesselation)
  • PCM matrix formation algorithm
  • strategy used to solve the PCM linear equations system.

Top section keywords

Units

Units of measure used in the input file. If Angstrom is given, all relevant input parameters are first converted in au and subsequently parsed.

  • Type: string
  • Valid values: AU | Angstrom
  • Default: No Default
CODATA

Set of fundamental physical constants to be used in the module.

  • Type: integer
  • Valid values: 2010 | 2006 | 2002 | 1998
  • Default: 2010

Cavity section keywords

Type

The type of the cavity. Completely specifies type of molecular surface and its discretization. Only one type is allowed. Restart cavity will read the file specified by NpzFile keyword and create a GePol cavity from that.

  • Type: string
  • Valid values: GePol | Restart
  • Default: none
NpzFile

The name of the .npz file to be used for the GePol cavity restart.

  • Type: string
  • Default: empty string
Area

Average area (weight) of the surface partition for the GePol (TsLess) cavity.

  • Type: double
  • Valid values: \(d \geq 0.01\,\text{a.u.}^2\)
  • Valid for: GePol cavity
  • Default value: \(0.3\,\text{a.u.}^2\)
Scaling

If true, the radii for the spheres will be scaled by 1.2. For finer control on the scaling factor for each sphere, select explicit creation mode.

  • Type: bool
  • Valid for: all cavities except Restart
  • Default value: True
RadiiSet

Select set of atomic radii to be used. Currently Bondi-Mantina [Bondi64][MantinaChamberlinValero+09], UFF [RCC+92] and Allinger’s MM3 [AZB94] sets available, see Available radii.

  • Type: string
  • Valid values: Bondi | UFF | Allinger
  • Valid for: all cavities except Restart
  • Default value: Bondi

Note

Radii in Allinger’s MM3 set are obtained by dividing the value in the original paper by 1.2, as done in the ADF COSMO implementation We advise to turn off scaling of the radii by 1.2 when using this set.

MinRadius

Minimal radius for additional spheres not centered on atoms. An arbitrarily big value is equivalent to switching off the use of added spheres, which is the default.

  • Type: double
  • Valid values: \(d \geq 0.4\,\text{a.u.}\)
  • Valid for: GePol cavity
  • Default value: \(100.0\,\text{a.u.}\)
Mode

How to create the list of spheres for the generation of the molecular surface:

  • in Implicit mode, the atomic coordinates and charges will be obtained from the QM host program. Spheres will be centered on the atoms and the atomic radii, as specified in one the built-in sets, will be used. Scaling by 1.2 will be applied according to the keyword Scaling;
  • in Atoms mode, the atomic coordinates and charges will be obtained from the QM host program. For the atoms specified by the array given in keyword Atoms, the built-in radii will be substituted by the radii provided in the keyword Radii. Scaling by 1.2 will be applied according to the keyword Scaling;
  • in Explicit mode, both centers and radii of the spheres are to be specified in the keyword Spheres. The user has full control over the generation of the list of spheres. Scaling by 1.2 is not applied, regardless of the value of the Scaling keyword.
  • Type: string
  • Valid values: Implicit | Atoms | Explicit
  • Valid for: all cavities except Restart
  • Default value: Implicit
Atoms

Array of atoms whose radius has to be substituted by a custom value.

  • Type: array of integers
  • Valid for: all cavities except Restart
Radii

Array of radii replacing the built-in values for the selected atoms.

  • Type: array of doubles
  • Valid for: all cavities except Restart
Spheres

Array of coordinates and centers for construction of the list of spheres in explicit mode. Format is \([\ldots, x_i, y_i, z_i, R_i, \ldots]\)

  • Type: array of doubles
  • Valid for: all cavities except Restart

Medium section keywords

SolverType

Type of solver to be used. All solvers are based on the Integral Equation Formulation of the Polarizable Continuum Model [CancesMennucci98]

  • IEFPCM. Collocation solver for a general dielectric medium
  • CPCM. Collocation solver for a conductor-like approximation to the dielectric medium
  • Type: string
  • Valid values: IEFPCM | CPCM
  • Default value: IEFPCM
Nonequilibrium

Initializes an additional solver using the dynamic permittivity. To be used in response calculations.

  • Type: bool
  • Valid for: all solvers
  • Default value: False
Solvent

Specification of the dielectric medium outside the cavity. This keyword must always be given a value. If the solvent name given is different from Explicit any other settings in the Green’s function section will be overridden by the built-in values for the solvent specified. See Table Available solvents for details. Solvent = Explicit, triggers parsing of the Green’s function sections.

  • Type: string

  • Valid values:

    • Water , H2O;
    • Propylene Carbonate , C4H6O3;
    • Dimethylsulfoxide , DMSO;
    • Nitromethane , CH3NO2;
    • Acetonitrile , CH3CN;
    • Methanol , CH3OH;
    • Ethanol , CH3CH2OH;
    • Acetone , C2H6CO;
    • 1,2-Dichloroethane , C2H4CL2;
    • Methylenechloride , CH2CL2;
    • Tetrahydrofurane , THF;
    • Aniline , C6H5NH2;
    • Chlorobenzene , C6H5CL;
    • Chloroform , CHCL3;
    • Toluene , C6H5CH3;
    • 1,4-Dioxane , C4H8O2;
    • Benzene , C6H6;
    • Carbon Tetrachloride , CCL4;
    • Cyclohexane , C6H12;
    • N-heptane , C7H16;
    • Explicit.
MatrixSymm

If True, the PCM matrix obtained by the IEFPCM collocation solver is symmetrized \(\mathbf{K} := \frac{\mathbf{K} + \mathbf{K}^\dagger}{2}\)

  • Type: bool
  • Valid for: IEFPCM solver
  • Default: True
Correction

Correction, \(k\) for the apparent surface charge scaling factor in the CPCM solver \(f(\varepsilon) = \frac{\varepsilon - 1}{\varepsilon + k}\)

  • Type: double
  • Valid values: \(k > 0.0\)
  • Valid for: CPCM solver
  • Default: 0.0
DiagonalIntegrator

Type of integrator for the diagonal of the boundary integral operators

  • Type: string
  • Valid values: COLLOCATION
  • Valid for: IEFPCM, CPCM
  • Default: COLLOCATION
  • Notes: in future releases we will add PURISIMA and NUMERICAL as options
DiagonalScaling

Scaling factor for diagonal of collocation matrices

  • Type: double
  • Valid values: \(f > 0.0\)
  • Valid for: IEFPCM, CPCM
  • Default: 1.07
  • Notes: values commonly used in the literature are 1.07 and 1.0694
ProbeRadius

Radius of the spherical probe approximating a solvent molecule. Used for generating the solvent-excluded surface (SES) or an approximation of it. Overridden by the built-in value for the chosen solvent.

  • Type: double
  • Valid values: \(d \in [0.1, 100.0]\,\text{a.u.}\)
  • Valid for: all solvers
  • Default: 1.0

Green section keywords

If Solvent = Explicit, two Green’s functions sections must be specified with tags inside and outside, i.e. Green<inside> and Green<outside>. The Green’s function inside will always be the vacuum, while the Green’s function outside might vary.

Type

Which Green’s function characterizes the medium.

  • Type: string
  • Valid values: Vacuum | UniformDielectric | SphericalDiffuse
  • Default: Vacuum
Der

How to calculate the directional derivatives of the Green’s function:

  • Numerical, perform numerical differentiation debug option;
  • Derivative, use automatic differentiation to get the directional derivative;
  • Gradient, use automatic differentiation to get the full gradient debug option;
  • Hessian, use automatic differentiation to get the full hessian debug option;
  • Type: string
  • Valid values: Numerical | Derivative | Gradient | Hessian
  • Default: Derivative

Note

The spherical diffuse Green’s function always uses numerical differentiation.

Eps

Static dielectric permittivity of the medium

  • Type: double
  • Valid values: \(\varepsilon \geq 1.0\)
  • Default: 1.0
EpsDyn

Dynamic dielectric permittivity of the medium

  • Type: double
  • Valid values: \(\varepsilon \geq 1.0\)
  • Default: 1.0
rofile

Functional form of the dielectric profile

  • Type: string
  • Valid values: Tanh | Erf
  • Valid for: SphericalDiffuse
  • Default: Tanh
Eps1

Static dielectric permittivity inside the interface

  • Type: double
  • Valid values: \(\varepsilon \geq 1.0\)
  • Default: 1.0
EpsDyn1

Dynamic dielectric permittivity inside the interface

  • Type: double
  • Valid values: \(\varepsilon \geq 1.0\)
  • Default: 1.0
Eps2

Static dielectric permittivity outside the interface

  • Type: double
  • Valid values: \(\varepsilon \geq 1.0\)
  • Default: 1.0
EpsDyn2

Dynamic dielectric permittivity outside the interface

  • Type: double
  • Valid values: \(\varepsilon \geq 1.0\)
  • Default: 1.0
Center

Center of the interface layer. This corresponds to the radius of the spherical droplet.

  • Type: double
  • Valid for: SphericalDiffuse
  • Default: 100.0 a.u.
Width

Physical width of the interface layer. This value is divided by 6.0 internally.

  • Type: double
  • Valid for: SphericalDiffuse
  • Default: 5.0 a.u.

Warning

Numerical instabilities may arise if a too small value is selected.

InterfaceOrigin

Center of the spherical droplet

  • Type: array of doubles
  • Valid for: SphericalDiffuse
  • Default: \([0.0, 0.0, 0.0]\)
MaxL

Maximum value of the angular momentum in the expansion of the Green’s function for the spherical diffuse Green’s function

  • Type: integer
  • Valid for: SphericalDiffuse
  • Default: 30

Molecule section keywords

It is possible to run the module standalone and use a classical charge distribution as specified in this section of the input. The run_pcm.x executable has to be compiled for a standalone run with:

python pcmsolver.py -x molecule.inp

where the molecule.inp input file looks like:

units = angstrom
codata = 2002
medium
{
	solvertype = cpcm
	correction = 0.5
    solvent = cyclohexane
}

cavity
{
	type = gepol
	area = 0.6
	radiiset = uff
	mode = implicit
}

molecule
{
    # x, y, z, q
    geometry = [0.000000000, 0.00000000,  0.08729478, 9.0,
                0.000000000, 0.00000000, -1.64558444, 1.0]
}
Geometry

Coordinates and charges of the molecular aggregate. Format is \([\ldots, x_i, y_i, z_i, Q_i, \ldots]\) Charges are always assumed to be in atomic units

  • Type: array of doubles

Available radii

_images/bondi_mantina.png _images/uff.png _images/allingerMM3.png

Available solvents

The macroscopic properties for the built-in list of solvents are:

  • static permittivity, \(\varepsilon_s\)
  • optical permittivity, \(\varepsilon_\infty\)
  • probe radius, \(r_\mathrm{probe}\) in Angstrom.

The following table summarizes the built-in solvents and their properties. Solvents are ordered by decreasing static permittivity.

Name Formula \(\varepsilon_s\) \(\varepsilon_\infty\) \(r_\mathrm{probe}\)
Water H2O 78.39 1.776 1.385
Propylene Carbonate C4H6O3 64.96 2.019 1.385
Dimethylsulfoxide DMSO 46.7 2.179 2.455
Nitromethane CH3NO2 38.20 1.904 2.155
Acetonitrile CH3CN 36.64 1.806 2.155
Methanol CH3OH 32.63 1.758 1.855
Ethanol CH3CH2OH 24.55 1.847 2.180
Acetone C2H6CO 20.7 1.841 2.38
1,2-Dichloroethane C2H4Cl2 10.36 2.085 2.505
Methylenechloride CH2Cl2 8.93 2.020 2.27
Tetrahydrofurane THF 7.58 1.971 2.9
Aniline C6H5NH2 6.89 2.506 2.80
Chlorobenzene C6H5Cl 5.621 2.320 2.805
Chloroform CHCl3 4.90 2.085 2.48
Toluene C6H5CH3 2.379 2.232 2.82
1,4-Dioxane C4H8O2 2.250 2.023 2.630
Benzene C6H6 2.247 2.244 2.630
Carbon tetrachloride CCl4 2.228 2.129 2.685
Cyclohexane C6H12 2.023 2.028 2.815
N-heptane C7H16 1.92 1.918 3.125

Interfacing a QM program and PCMSolver

For the impatients: tl;dr

In these examples, we want to show how every function in the API works. If your program is written in Fortran, head over to Interfacing with a Fortran host If your program is written in C/C++, head over to Interfacing with a C host

How PCMSolver handles potentials and charges: surface functions

Electrostatic potential vectors and the corresponding apparent surface charge vectors are handled internally as surface functions. The actual values are stored into Eigen vectors and saved into a map. The mapping is between the name of the surface function, given by the programmer writing the interface to the library, and the vector holding the values.

What you should care about: API functions

These are the contents of the pcmsolver.h file defining the public API of the PCMSolver library. The Fortran bindings for the API are in the pcmsolver.f90 file. The indexing of symmetry operations and their mapping to a bitstring is explained in the following Table. This is important when passing symmetry information to the pcmsolver_new() function.

Symmetry operations indexing within the module
Index zyx Generator Parity
0 000 E 1.0
1 001 Oyz -1.0
2 010 Oxz -1.0
3 011 C2z 1.0
4 100 Oxy -1.0
5 101 C2y 1.0
6 110 C2x 1.0
7 111 i -1.0

Defines

PCMSolver_EXPORT

PCMSolver, an API for the Polarizable Continuum Model Copyright (C) 2017 Roberto Di Remigio, Luca Frediani and collaborators.

This file is part of PCMSolver.

PCMSolver is free software: you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.

PCMSolver is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details.

You should have received a copy of the GNU Lesser General Public License along with PCMSolver. If not, see http://www.gnu.org/licenses/.

For information on the complete list of contributors to the PCMSolver API, see: http://pcmsolver.readthedocs.io/

pcmsolver_bool_t_DEFINED

Typedefs

typedef bool pcmsolver_bool_t
typedef pcmsolver_context_t

Workaround to have pcmsolver_context_t available to C

typedef HostWriter

Flushes module output to host program

Parameters
  • message: contents of the module output

Enums

enum pcmsolver_reader_t

Input processing strategies.

Values:

PCMSOLVER_READER_OWN

Module reads input on its own

PCMSOLVER_READER_HOST

Module receives input from host

Functions

pcmsolver_context_t *pcmsolver_new(pcmsolver_reader_t input_reading, int nr_nuclei, double charges[], double coordinates[], int symmetry_info[], struct PCMInput *host_input, HostWriter writer)

Creates a new PCM context object.

The molecular point group information is passed as an array of 4 integers: number of generators, first, second and third generator respectively. Generators map to integers as in table :ref:symmetry-ops

Parameters
  • input_reading: input processing strategy
  • nr_nuclei: number of atoms in the molecule
  • charges: atomic charges
  • coordinates: atomic coordinates
  • symmetry_info: molecular point group information
  • host_input: input to the module, as read by the host

void pcmsolver_delete(pcmsolver_context_t *context)

Deletes a PCM context object.

Parameters
  • context: the PCM context object to be deleted

pcmsolver_bool_t pcmsolver_is_compatible_library(void)

Whether the library is compatible with the header file Checks that the compiled library and header file version match. Host should abort when that is not the case.

Warning
This function should be called before instantiating any PCM context objects.

void pcmsolver_print(pcmsolver_context_t *context)

Prints citation and set up information.

Parameters
  • context: the PCM context object

int pcmsolver_get_cavity_size(pcmsolver_context_t *context)

Getter for the number of finite elements composing the molecular cavity.

Return
the size of the cavity
Parameters
  • context: the PCM context object

int pcmsolver_get_irreducible_cavity_size(pcmsolver_context_t *context)

Getter for the number of irreducible finite elements composing the molecular cavity.

Return
the number of irreducible finite elements
Parameters
  • context: the PCM context object

void pcmsolver_get_centers(pcmsolver_context_t *context, double centers[])

Getter for the centers of the finite elements composing the molecular cavity.

Parameters
  • context: the PCM context object
  • centers: array holding the coordinates of the finite elements centers

void pcmsolver_get_center(pcmsolver_context_t *context, int its, double center[])

Getter for the center of the i-th finite element.

Parameters
  • context: the PCM context object
  • its: index of the finite element
  • center: array holding the coordinates of the finite element center

void pcmsolver_get_areas(pcmsolver_context_t *context, double areas[])

Getter for the areas/weights of the finite elements.

Parameters
  • context: the PCM context object
  • areas: array holding the weights/areas of the finite elements

void pcmsolver_compute_asc(pcmsolver_context_t *context, const char *mep_name, const char *asc_name, int irrep)

Computes ASC given a MEP and the desired irreducible representation.

Parameters
  • context: the PCM context object
  • mep_name: label of the MEP surface function
  • asc_name: label of the ASC surface function
  • irrep: index of the desired irreducible representation The module uses the surface function concept to handle potentials and charges. Given labels for each, this function retrieves the MEP and computes the corresponding ASC.

void pcmsolver_compute_response_asc(pcmsolver_context_t *context, const char *mep_name, const char *asc_name, int irrep)

Computes response ASC given a MEP and the desired irreducible representation.

Parameters
  • context: the PCM context object
  • mep_name: label of the MEP surface function
  • asc_name: label of the ASC surface function
  • irrep: index of the desired irreducible representation If Nonequilibrium = True in the input, calculates a response ASC using the dynamic permittivity. Falls back to the solver with static permittivity otherwise.

double pcmsolver_compute_polarization_energy(pcmsolver_context_t *context, const char *mep_name, const char *asc_name)

Computes the polarization energy.

Return
the polarization energy This function calculates the dot product of the given MEP and ASC vectors.
Parameters
  • context: the PCM context object
  • mep_name: label of the MEP surface function
  • asc_name: label of the ASC surface function

double pcmsolver_get_asc_dipole(pcmsolver_context_t *context, const char *asc_name, double dipole[])

Getter for the ASC dipole.

Return
the ASC dipole, i.e. { ^2}
Parameters
  • context: the PCM context object
  • asc_name: label of the ASC surface function
  • dipole: the Cartesian components of the ASC dipole moment

void pcmsolver_get_surface_function(pcmsolver_context_t *context, int size, double values[], const char *name)

Retrieves data wrapped in a given surface function.

Parameters
  • context: the PCM context object
  • size: the size of the surface function
  • values: the values wrapped in the surface function
  • name: label of the surface function

void pcmsolver_set_surface_function(pcmsolver_context_t *context, int size, double values[], const char *name)

Sets a surface function given data and label.

Parameters
  • context: the PCM context object
  • size: the size of the surface function
  • values: the values to be wrapped in the surface function
  • name: label of the surface function

void pcmsolver_print_surface_function(pcmsolver_context_t *context, const char *name)

Prints surface function contents to host output.

Parameters
  • context: the PCM context object
  • name: label of the surface function

void pcmsolver_save_surface_functions(pcmsolver_context_t *context)

Dumps all currently saved surface functions to NumPy arrays in .npy files.

Parameters
  • context: the PCM context object

void pcmsolver_save_surface_function(pcmsolver_context_t *context, const char *name)

Dumps a surface function to NumPy array in .npy file.

Parameters
  • context: the PCM context object
  • name: label of the surface function

void pcmsolver_load_surface_function(pcmsolver_context_t *context, const char *name)

Loads a surface function from a .npy file.

Parameters
  • context: the PCM context object
  • name: label of the surface function

void pcmsolver_write_timings(pcmsolver_context_t *context)

Writes timing results for the API.

Parameters
  • context: the PCM context object

Host input forwarding

struct PCMInput

Data structure for host-API input communication.

Forward-declare PCMInput input wrapping struct

Public Members

char cavity_type[8]

Type of cavity requested.

int patch_level

Wavelet cavity mesh patch level.

double coarsity

Wavelet cavity mesh coarsity.

double area

Average tesserae area.

char radii_set[8]

The built-in radii set to be used.

double min_distance

Minimal distance between sampling points.

int der_order

Derivative order for the switching function.

pcmsolver_bool_t scaling

Whether to scale or not the atomic radii.

char restart_name[20]

Name of the .npz file for GePol cavity restart.

double min_radius

Minimal radius for the added spheres.

char solver_type[7]

Type of solver requested.

double correction

Correction in the CPCM apparent surface charge scaling factor.

char solvent[16]

Name of the solvent.

double probe_radius

Radius of the spherical probe mimicking the solvent.

char equation_type[11]

Type of the integral equation to be used.

char inside_type[7]

Type of Green’s function requested inside the cavity.

double outside_epsilon

Value of the static permittivity outside the cavity.

char outside_type[22]

Type of Green’s function requested outside the cavity.

Internal details of the API

Warning

doxygenclass: Cannot find class “Meddle” in doxygen xml output for project “PCMSolver” from directory: xml

Interfacing with a Fortran host

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
!
! PCMSolver, an API for the Polarizable Continuum Model
! Copyright (C) 2017 Roberto Di Remigio, Luca Frediani and collaborators.
!
! This file is part of PCMSolver.
!
! PCMSolver is free software: you can redistribute it and/or modify
! it under the terms of the GNU Lesser General Public License as published by
! the Free Software Foundation, either version 3 of the License, or
! (at your option) any later version.
!
! PCMSolver is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
! GNU Lesser General Public License for more details.
!
! You should have received a copy of the GNU Lesser General Public License
! along with PCMSolver.  If not, see <http://www.gnu.org/licenses/>.
!
! For information on the complete list of contributors to the
! PCMSolver API, see: <http://pcmsolver.readthedocs.io/>
!

program pcm_fortran_host

      use, intrinsic :: iso_c_binding
      use, intrinsic :: iso_fortran_env, only: output_unit, error_unit
      use pcmsolver
      use utilities
      use testing

      implicit none

      type(c_ptr) :: pcm_context
      integer(c_int) :: nr_nuclei
      real(c_double), allocatable :: charges(:)
      real(c_double), allocatable :: coordinates(:)
      integer(c_int) :: symmetry_info(4)
      type(PCMInput) :: host_input
      logical :: log_open, log_exist
      ! Shows two different, but equivalent ways of defining labels for surface functions
      character(kind=c_char, len=1) :: mep_lbl(7) = (/'N', 'u', 'c', 'M', 'E', 'P', c_null_char/)
      character(kind=c_char, len=1) :: asc_lbl(7) = (/'N', 'u', 'c', 'A', 'S', 'C', c_null_char/)
      character(kind=c_char, len=1) :: asc_B3g_lbl(7) = (/'O', 'I', 'T', 'A', 'S', 'C', c_null_char/)
      character(kind=c_char, len=1) :: asc_neq_B3g_lbl(7) = (/'A', 'S', 'C', 'B', '3', 'g', c_null_char/)
      real(c_double), allocatable :: grid(:), mep(:), asc_Ag(:), asc_B3g(:), asc_neq_B3g(:), areas(:)
      integer(c_int) :: irrep
      integer(c_int) :: grid_size, irr_grid_size
      real(c_double) :: energy
      ! Reference values for scalar quantities
      integer(c_int), parameter :: ref_size = 576, ref_irr_size = 72
      real(c_double), parameter :: ref_energy = -0.437960027982

      if (.not. pcmsolver_is_compatible_library()) then
        write(error_unit, *) 'PCMSolver library not compatible!'
        stop
      end if

      ! Open a file for the output...
      inquire(file = 'fortran_host.log', opened = log_open, &
        exist = log_exist)
      if (log_exist) then
        open(unit = output_unit,                   &
          file = 'fortran_host.log',       &
          status = 'unknown',       &
          form = 'formatted',       &
          access = 'sequential')
        close(unit = output_unit, status = 'delete')
      end if
      open(unit = output_unit,                      &
        file = 'fortran_host.log',          &
        status = 'new',              &
        form = 'formatted',          &
        access = 'sequential')
      rewind(output_unit)
      write(output_unit, *) 'Starting a PCMSolver calculation'

      nr_nuclei = 6_c_int
      allocate(charges(nr_nuclei))
      allocate(coordinates(3*nr_nuclei))

      ! Use C2H4 in D2h symmetry
      charges = (/ 6.0_c_double, 1.0_c_double, 1.0_c_double, &
                   6.0_c_double, 1.0_c_double, 1.0_c_double /)
      coordinates = (/ 0.0_c_double,  0.0_c_double,       1.257892_c_double, &
                      0.0_c_double,  1.745462_c_double,  2.342716_c_double, &
                      0.0_c_double, -1.745462_c_double,  2.342716_c_double, &
                      0.0_c_double,  0.0_c_double,      -1.257892_c_double, &
                      0.0_c_double,  1.745462_c_double, -2.342716_c_double, &
                      0.0_c_double, -1.745462_c_double, -2.342716_c_double /)

      ! This means the molecular point group has three generators:
      ! the Oxy, Oxz and Oyz planes
      symmetry_info = (/ 3, 4, 2, 1 /)

      host_input = pcmsolver_input()

      pcm_context = pcmsolver_new(PCMSOLVER_READER_HOST,           &
                                  nr_nuclei, charges, coordinates, &
                                  symmetry_info, host_input,       &
                                  c_funloc(host_writer))

      call pcmsolver_print(pcm_context)

      grid_size = pcmsolver_get_cavity_size(pcm_context)
      irr_grid_size = pcmsolver_get_irreducible_cavity_size(pcm_context)
      allocate(grid(3*grid_size))
      grid = 0.0_c_double
      call pcmsolver_get_centers(pcm_context, grid)
      allocate(areas(grid_size))
      call pcmsolver_get_areas(pcm_context, areas)

      allocate(mep(grid_size))
      mep = 0.0_c_double
      mep = nuclear_mep(nr_nuclei, charges, reshape(coordinates, (/ 3_c_int, nr_nuclei /)), &
                        grid_size, reshape(grid, (/ 3_c_int, grid_size /)))
      call pcmsolver_set_surface_function(pcm_context, grid_size, mep, mep_lbl)
      ! This is the Ag irreducible representation (totally symmetric)
      irrep = 0_c_int
      call pcmsolver_compute_asc(pcm_context, mep_lbl, asc_lbl, irrep)
      allocate(asc_Ag(grid_size))
      asc_Ag = 0.0_c_double
      call pcmsolver_get_surface_function(pcm_context, grid_size, asc_Ag, asc_lbl)

      energy = pcmsolver_compute_polarization_energy(pcm_context, mep_lbl, asc_lbl)

      write(output_unit, '(A, F20.12)') 'Polarization energy = ', energy

      allocate(asc_neq_B3g(grid_size))
      asc_neq_B3g = 0.0_c_double
      ! This is the B3g irreducible representation
      irrep = 3_c_int
      call pcmsolver_compute_response_asc(pcm_context, mep_lbl, asc_neq_B3g_lbl, irrep)
      call pcmsolver_get_surface_function(pcm_context, grid_size, asc_neq_B3g, asc_neq_B3g_lbl)

      ! Equilibrium ASC in B3g symmetry.
      ! This is an internal check: the relevant segment of the vector
      ! should be the same as the one calculated using pcmsolver_compute_response_asc
      allocate(asc_B3g(grid_size))
      asc_B3g = 0.0_c_double
      ! This is the B3g irreducible representation
      irrep = 3_c_int
      call pcmsolver_compute_asc(pcm_context, mep_lbl, asc_B3g_lbl, irrep)
      call pcmsolver_get_surface_function(pcm_context, grid_size, asc_B3g, asc_B3g_lbl)

      ! Check that everything calculated is OK
      ! Cavity size
      if (grid_size .ne. ref_size) then
        write(error_unit, *) 'Error in the cavity size, please file an issue on: https://github.com/PCMSolver/pcmsolver'
        stop
      else
        write(output_unit, *) 'Test on cavity size: PASSED'
      end if
      ! Irreducible cavity size
      if (irr_grid_size .ne. ref_irr_size) then
        write(error_unit, *) 'Error in the irreducible cavity size, please file an issue on: https://github.com/PCMSolver/pcmsolver'
        stop
      else
        write(output_unit, *) 'Test on irreducible cavity size: PASSED'
      end if
      ! Polarization energy
      if (.not. check_unsigned_error(energy, ref_energy, 1.0e-7_c_double)) then
        write(error_unit, *) 'Error in the polarization energy, please file an issue on: https://github.com/PCMSolver/pcmsolver'
        stop
      else
        write(output_unit, *) 'Test on polarization energy: PASSED'
      end if
      ! Surface functions
      call test_surface_functions(grid_size, mep, asc_Ag, asc_B3g, asc_neq_B3g, areas)

      call pcmsolver_write_timings(pcm_context)

      call pcmsolver_delete(pcm_context)

      deallocate(charges)
      deallocate(coordinates)
      deallocate(grid)
      deallocate(mep)
      deallocate(asc_Ag)
      deallocate(asc_B3g)
      deallocate(asc_neq_B3g)
      deallocate(areas)

      close(output_unit)

end program pcm_fortran_host

Interfacing with a C host

Warning

Multidimensional arrays are handled in column-major ordering (i.e. Fortran ordering) by the module.

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
/**
 * PCMSolver, an API for the Polarizable Continuum Model
 * Copyright (C) 2016 Roberto Di Remigio, Luca Frediani and collaborators.
 *
 * This file is part of PCMSolver.
 *
 * PCMSolver is free software: you can redistribute it and/or modify
 * it under the terms of the GNU Lesser General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * PCMSolver is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU Lesser General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public License
 * along with PCMSolver.  If not, see <http://www.gnu.org/licenses/>.
 *
 * For information on the complete list of contributors to the
 * PCMSolver API, see: <http://pcmsolver.readthedocs.io/>
 */

#include <stddef.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include "pcmsolver.h"
#include "PCMInput.h"

#include "C_host-functions.h"

#define NR_NUCLEI 6

FILE * output;

void host_writer(const char * message) { fprintf(output, "%s\n", message); }

int main() {

  output = fopen("C_host.log", "w+");
  if (!pcmsolver_is_compatible_library()) {
    fprintf(stderr, "%s\n", "PCMSolver library not compatible");
    exit(EXIT_FAILURE);
  }

  fprintf(output, "%s\n", "Starting a PCMSolver calculation");
  // Use C2H4 in D2h symmetry
  double charges[NR_NUCLEI] = {6.0, 1.0, 1.0, 6.0, 1.0, 1.0};
  double coordinates[3 * NR_NUCLEI] = {
      0.0, 0.000000, 1.257892,  0.0, 1.745462, 2.342716,  0.0, -1.745462, 2.342716,
      0.0, 0.000000, -1.257892, 0.0, 1.745462, -2.342716, 0.0, -1.745462, -2.342716};
  // This means the molecular point group has three generators:
  // the Oxy, Oxz and Oyz planes
  int symmetry_info[4] = {3, 4, 2, 1};
  struct PCMInput host_input = pcmsolver_input();

  pcmsolver_context_t * pcm_context =
      pcmsolver_new(PCMSOLVER_READER_HOST, NR_NUCLEI, charges, coordinates,
                    symmetry_info, &host_input, host_writer);

  pcmsolver_print(pcm_context);

  int grid_size = pcmsolver_get_cavity_size(pcm_context);
  int irr_grid_size = pcmsolver_get_irreducible_cavity_size(pcm_context);
  double * grid = (double *)calloc(3 * grid_size, sizeof(double));
  pcmsolver_get_centers(pcm_context, grid);
  double * areas = (double *)calloc(grid_size, sizeof(double));
  pcmsolver_get_areas(pcm_context, areas);

  double * mep = nuclear_mep(NR_NUCLEI, charges, coordinates, grid_size, grid);
  const char * mep_lbl = {"NucMEP"};
  pcmsolver_set_surface_function(pcm_context, grid_size, mep, mep_lbl);
  const char * asc_lbl = {"NucASC"};
  // This is the Ag irreducible representation (totally symmetric)
  int irrep = 0;
  pcmsolver_compute_asc(pcm_context, mep_lbl, asc_lbl, irrep);
  double * asc_Ag = (double *)calloc(grid_size, sizeof(double));
  pcmsolver_get_surface_function(pcm_context, grid_size, asc_Ag, asc_lbl);

  double energy =
      pcmsolver_compute_polarization_energy(pcm_context, mep_lbl, asc_lbl);

  fprintf(output, "Polarization energy: %20.12f\n", energy);

  double * asc_neq_B3g = (double *)calloc(grid_size, sizeof(double));
  const char * asc_neq_B3g_lbl = {"OITASC"};
  // This is the B3g irreducible representation
  irrep = 3;
  pcmsolver_compute_response_asc(pcm_context, mep_lbl, asc_neq_B3g_lbl, irrep);
  pcmsolver_get_surface_function(pcm_context, grid_size, asc_neq_B3g,
                                 asc_neq_B3g_lbl);

  // Equilibrium ASC in B3g symmetry.
  // This is an internal check: the relevant segment of the vector
  // should be the same as the one calculated using pcmsolver_compute_response_asc
  double * asc_B3g = (double *)calloc(grid_size, sizeof(double));
  const char * asc_B3g_lbl = {"ASCB3g"};
  pcmsolver_compute_asc(pcm_context, mep_lbl, asc_B3g_lbl, irrep);
  pcmsolver_get_surface_function(pcm_context, grid_size, asc_B3g, asc_B3g_lbl);

  // Check that everything calculated is OK
  // Cavity size
  const int ref_size = 576;
  if (grid_size != ref_size) {
    fprintf(stderr, "%s\n", "Error in the cavity size, please file an issue on: "
                            "https://github.com/PCMSolver/pcmsolver");
    exit(EXIT_FAILURE);
  } else {
    fprintf(output, "%s\n", "Test on cavity size: PASSED");
  }
  // Irreducible cavity size
  const int ref_irr_size = 72;
  if (irr_grid_size != ref_irr_size) {
    fprintf(stderr, "%s\n", "Error in the irreducible cavity size, please file an "
                            "issue on: https://github.com/PCMSolver/pcmsolver");
    exit(EXIT_FAILURE);
  } else {
    fprintf(output, "%s\n", "Test on irreducible cavity size: PASSED");
  }
  // Polarization energy
  const double ref_energy = -0.437960027982;
  if (!check_unsigned_error(energy, ref_energy, 1.0e-7)) {
    fprintf(stderr, "%s\n", "Error in the polarization energy, please file an issue "
                            "on: https://github.com/PCMSolver/pcmsolver");
    exit(EXIT_FAILURE);
  } else {
    fprintf(output, "%s\n", "Test on polarization energy: PASSED");
  }
  // Surface functions
  test_surface_functions(output, grid_size, mep, asc_Ag, asc_B3g, asc_neq_B3g,
                         areas);

  pcmsolver_write_timings(pcm_context);

  pcmsolver_delete(pcm_context);

  free(grid);
  free(mep);
  free(asc_Ag);
  free(asc_B3g);
  free(asc_neq_B3g);
  free(areas);

  fclose(output);

  return EXIT_SUCCESS;
}

Publications

Theses

Presentations

Posters

Other stuff

PCMSolver Programmers’ Manual

General Structure

_images/structure.png

External libraries:

  • parts of the C++ Boost libraries are used to provide various functionality, like timing and metaprogramming. The source for the 1.54.0 release is shipped with the module’s source code. Some of the libraries used need to be compiled. Boost is released under the terms of the Boost Software License, v1.0 (see also http://www.boost.org/users/license.html) We encourage the use of Boost whenever some functionality has already been coded within those libraries. However, consider carefully the introduction of functionality depending on compiler Boost libraries.
  • the Eigen template library for linear algebra. Almost every operation involving matrices and vectors is performed through Eigen. Eigen provides convenient type definitions for vectors and matrices (of arbitrary dimensions) and the corresponding operations. Have a look here for a quick reference guide to the API and at the getting started guide to get started. Eigen is distributed under the terms of the Mozilla Public License, v2.0
  • the Getkw library by Jonas Juselius is used to manage input. It is distributed under the terms of the GNU General Public License, v2.0
  • the libtaylor library implementing automatic differentiation and available under the terms of the MIT License.

Third-party code snippets:

  • Fortran subroutines dsyevv3, dsyevh3, dsyevj3 for the diagonalization of 3x3 Hermitian matrices. These subroutines were copied verbatim from the source code provided by Joachim Kopp and described in [Kop08] (also available on the arXiv) The diagonalization subroutines are made available under the terms of the GNU Lesser General Public License, v2.1
  • C++ cnpy library for saving arrays in C++ into Numpy arrays. The library is from Carl Rogers under the terms of the MIT License. The version in PCMSolver is slightly different.

Coding standards

General Object-Oriented design principles you should try to follow:
  1. Identify the aspects of your application that vary and separate them from what stays the same;
  2. Program to an interface, not an implementation;
  3. Favor composition over inheritance;
  4. Strive for loosely coupled designs between objects that interact;
  5. Classes should be open for extension, but closed for modification;
  6. Depend upon abstractions. Do not depend upon concrete classes;
  7. Principle of Least Knowledge. Talk only to your immediate friends;

[SA04][CGL98][Cli]

Including header files

Do not include header files unnecessarily. Even if PCMSolver is not a big project, unnecessary include directives and/or forward declarations introduce nasty interdependencies among different parts of the code. This reflects mainly in longer compilation times, but also in uglier looking code (see also the discussion in [Sut99]).

Follow these guidelines to decide whether to include or forward declare:
  1. class A makes no reference to class B. Neither include nor forward declare B;
  2. class A refers to class B as a friend. Neither include nor forward declare B;
  3. class A contains a pointer/reference to a class B object. Forward declare B;
  4. class A contains functions with a class B object (value/pointer/reference) as parameter/return value. Forward declare B;
  5. class A is derived from class B. include B;
  6. class A contains a class B object. include B.
#ifndef MYCLASS_HPP
#define MYCLASS_HPP

//==============================
// Forward declared dependencies
class Foo;
class Bar;

//==============================
// Included dependencies
#include <vector>
#include "Parent.hpp"

//==============================
// The actual class
class MyClass : public Parent // Parent object, so #include "Parent.h"
{
  public:
    std::vector<int> avector; // vector object, so #include <vector>
    Foo * foo;                // Foo pointer, so forward declare
    void Func(Bar & bar);     // Bar reference as parameter, so forward declare

    friend class MyFriend;    // friend declaration is not a dependency
                              //    don't do anything about MyFriend
};

#endif // MYCLASS_HPP

Proper overloading of operator<<

Suppose we have an inheritance hierarchy made of an abstract base class, Base, and two derived classes, Derived1 and Derived2. In the Base class header file we will define a pure virtual private function printObject and provide a public friend overload of operator<<:

#include <iosfwd>

class Base
{
  public:
    // All your other very fancy public members
    friend std::ostream & operator<<(std::ostream & os, Base & base)
    {
            return base.printObject(os);
    }
  protected:
    // All your other very fancy protected members
  private:
    // All your other very fancy private members
    virtual std::ostream & printObject(std::ostream & os) = 0;
}

The printObject method can also be made (impure) virtual, it really depends on your class hierarchy. Derived1 and Derived2 header files will provide a public friend overload of operator<< (friendliness isn’t inherited, transitive or reciprocal) and an override for the printObject method:

#include <iosfwd>

#include "Base.hpp"

class Derived1 : public Base
{
  public:
    // All your other very fancy public members
    friend std::ostream & operator<<(std::ostream & os, Derived1 & derived)
    {
      return derived.printObject(os);
    }
  protected:
    // All your other very fancy protected members
  private:
    // All your other very fancy private members
    virtual std::ostream & printObject(std::ostream & os);
}

class Derived2 : public Base
{
  public:
    // All your other very fancy public members
    friend std::ostream & operator<<(std::ostream & os, Derived2 & derived)
    {
      return derived.printObject(os);
    }
  protected:
    // All your other very fancy protected members
  private:
    // All your other very fancy private members
    virtual std::ostream & printObject(std::ostream & os);
}

Code formatting

We conform to the so-called Linux (aka kernel) formatting style for C/C++ code (see http://en.wikipedia.org/wiki/Indent_style#Kernel_style) with minimal modifications. If uncertain on your code formatting use the Artistic Style command-line utility to rectify it:

astyle --style=linux --max-code-length=85 --indent-namespaces --keep-one-line-blocks filename

Documentation

This documentation is generated using Sphinx and Doxygen The two softwares are bridged by means of the Breathe extension The online version of this documentation is built and served by Read The Docs. The webpage http://pcmsolver.readthedocs.org/ is updated on each push to the public GitHub repository.

How and what to document

Doxygen enables documenting the code in the source code files thus removing a “barrier” for developers. To avoid that the code degenerates into a Big Ball of Mud, it is mandatory to document directly within the source code classes and functions. To document general programming principles, design choices, maintenance etc. you can create a .rst file in the doc directory. Remember to refer the new file inside the index.rst file (it won’t be parsed otherwise). Sphing uses reStructuredText and Markdown. Support for Markdown is not as extensive as for reStructuredText, see these comments Follow the guidelines in [WAB+14] regarding what to document.

Documeting methods in derived classes

Virtual methods should only be documented in the base classes. This avoids unnecessary verbosity and conforms to the principle: “Document _what_, not _how_” [WAB+14] If you feel the _how_ needs to be explicitly documented, add some notes in the appropriate .rst file.

How does this work?

To have an offline version of the documentation just issue make doc in the build directory. The HTML will be stored in doc/html. Open the doc/html/index.html file with your browser to see and browse the documentation.

Warning

Building the documentation requires Doxygen, Sphinx, Perl and the Python modules PyYAML, Breathe and Matplotlib.

CMake usage

This is a brief guide to our CMake infrastructure which is managed via Autocmake

Warning

The minimum required CMake version is 2.8.8

Adding new source subdirectories and/or files

Developers HAVE TO manually list the sources in a given subdirectory of the main source directory src/. In our previous infrastructure this was not necessary, but the developers needed to trigger CMake to regenerate the Makefiles manually.

New subdirectory

First of all, you will have to let CMake know that a new source-containing subdirectory has been added to the source tree. Due to the hierarchical approach CMake is based upon you will need to modify the CMakeLists.txt in the src directory and create a new one in your new subdirectory. For the first step:

1. if your new subdirectory contains header files, add a line like the following to the CMakeLists.txt file contained in the src directory:

${CMAKE_CURRENT_LIST_DIR}/subdir_name

to the command setting the list of directories containing headers. This sets up the list of directories where CMake will look for headers with definitions of classes and functions. If your directory contains Fortran code you can skip this step;

2. add a line like the following to the CMakeLists.txt file contained in the src directory:

add_subdirectory(subdir_name)

This will tell CMake to go look inside subdir_name for a CMakeLists.txt containing more sets of instructions. It is preferable to add these new lines in alphabetic order

Inside your new subdirectory you will need to add a CMakeLists.txt file containing the set of instructions to build your cutting edge code. This is the second step. Run the make_cmake_files.py Python script in the src/ directory:

python make_cmake_files.py --libname=cavity --lang=CXX

to generate a template CMakeLists.txt.try file:

# List of headers
list(APPEND headers_list Cavity.hpp ICavity.hpp Element.hpp GePolCavity.hpp RegisterCavityToFactory.hpp RestartCavity.hpp)

# List of sources
list(APPEND sources_list ICavity.cpp Element.cpp GePolCavity.cpp RestartCavity.cpp)

add_library(cavity OBJECT ${sources_list} ${headers_list})
set_target_properties(cavity PROPERTIES POSITION_INDEPENDENT_CODE 1 )
set_property(GLOBAL APPEND PROPERTY PCMSolver_HEADER_DIRS ${CMAKE_CURRENT_LIST_DIR})
# Sets install directory for all the headers in the list
foreach(_header ${headers_list})
   install(FILES ${_header} DESTINATION include/cavity)
endforeach()

The template might need additional editing. Each source subdirectory is the lowest possible in the CMake hierarchy and it contains set of instructions for:

  1. exporting a list of header files (.h or .hpp) to the upper level in the hierarchy, possibly excluding some of them
  2. define install targets for the files in this subdirectory.

All the source files are compiled into the unique static library libpcm.a and unique dynamic library libpcm.so. This library is the one the host QM program need to link.

Searching for libraries

In general, the use of the find_package macro is to be preferred, as it is standardized and ensured to work on any platform. Use of find_package requires that the package/library you want to use has already a module inside the CMake distribution. If that’s not the case, you should never use the following construct for third-party libraries:

target_link_libraries(myexe -lsomesystemlib)

If the library does not exist, the end result is a cryptic linker error. See also Jussi Pakkanen’s blog You will first need to find the library, using the macro find_library, and then use the target_link_libraries command.

Maintenance

Description and how-to for maintenance operations. Some of the maintenance scripts have been moved to the pcmsolvermeta repository

Bump version

Version numbering follows the guidelines of semantic versioning To update, change the relevant field in the README.md file.

Changelog

We follow the guidelines of Keep a CHANGELOG On all but the release branches, there is an Unreleased section under which new additions should be listed. To simplify perusal of the CHANGELOG.md, use the following subsections:

  1. Added for new features.
  2. Changed for changes in existing functionality.
  3. Deprecated for once-stable features removed in upcoming releases.
  4. Removed for deprecated features removed in this release.
  5. Fixed for any bug fixes.
  6. Security to invite users to upgrade in case of vulnerabilities.

Updating Eigen distribution

The C++ linear algebra library Eigen comes bundled with the module. To update the distributed version one has to:

  1. download the desired version of the library to a scratch location. Eigen’s website is: http://eigen.tuxfamily.org/

  2. unpack the downloaded archive;

  3. go into the newly created directory and create a build directory;

  4. go into the newly created build directory and type the following (remember to substitute @PROJECT_SOURCE_DIR@ with the actual path)

    cmake .. -DCMAKE_INSTALL_PREFIX=@PROJECT_SOURCE_DIR@/external/eigen3
    

Remember to commit and push your modifications.

Release process

We have two repositories one public for the release, hosted on GitHub and one private for the development, hosted on GitLab. At release time the master branch on the private repository is synced to that of the public repository.

Warning

This means that WHATEVER is on master at release time is considered ready for release. Protection of functionality happens EXCLUSIVELY by making use of branches/forks on the private repository.

You need to compile the to-be-released code and run the unit test suite. If compilation works and all unit tests are passing then the code is ready to be released:

git push Origin release

Notice that Origin has been spelled with a capital O the reason being that the release branch gets pushed both to the private and the public repositories (trick explanation) In brief, you need to have a .git/config file that resembles the following:

[remote "origin"]
    url = git@gitlab.com:PCMSolver/pcmsolver.git
    fetch = +refs/heads/*:refs/remotes/origin/*
[remote "GitHub"]
    url = git@github.com:PCMSolver/pcmsolver.git
    fetch = +refs/heads/*:refs/remotes/GitHub/*
[remote "Origin"]
    url = git@gitlab.com:PCMSolver/pcmsolver.git
    url = git@github.com:PCMSolver/pcmsolver.git

Testing

We perform unit testing of our API. The unit testing framework used is Catch The framework provides quite an extensive set of macros to test various data types, it also provides facilities for easily setting up test fixtures. Usage is extremely simple and the documentation is very well written. For a quick primer on how to use Catch refer to: https://github.com/philsquared/Catch/blob/master/docs/tutorial.md The basic idea of unit testing is to test each building block of the code separataly. In our case, the term “building block” is used to mean a class.

To add new tests for your class you have to:

  1. create a new subdirectory inside tests/ and add a line like the following to the CMakeLists.txt

    add_subdirectory(new_subdir)
    
  2. create a CMakeLists.txt inside your new subdirectory. This CMakeLists.txt adds the source for a given unit test to the global UnitTestsSources property and notifies CTest that a test with given name is part of the test suite. The generation of the CMakeLists.txt can be managed by make_cmake_files.py Python script. This will take care of also setting up CTest labels. This helps in further grouping the tests for our convenience. Catch uses tags to index tests and tags are surrounded by square brackets. The Python script inspects the sources and extracts labels from Catch tags. The add_Catch_test CMake macro takes care of the rest.

    We require that each source file containing tests follows the naming convention new_subdir_testname and that testname gives some clue to what is being tested. Depending on the execution of tests in a different subdirectory is bad practice. A possible workaround is to add some kind of input file and create a text fixture that sets up the test environment. Have a look in the tests/input directory for an example

  3. create the .cpp files containing the tests. Use the following template:

      1
      2
      3
      4
      5
      6
      7
      8
      9
     10
     11
     12
     13
     14
     15
     16
     17
     18
     19
     20
     21
     22
     23
     24
     25
     26
     27
     28
     29
     30
     31
     32
     33
     34
     35
     36
     37
     38
     39
     40
     41
     42
     43
     44
     45
     46
     47
     48
     49
     50
     51
     52
     53
     54
     55
     56
     57
     58
     59
     60
     61
     62
     63
     64
     65
     66
     67
     68
     69
     70
     71
     72
     73
     74
     75
     76
     77
     78
     79
     80
     81
     82
     83
     84
     85
     86
     87
     88
     89
     90
     91
     92
     93
     94
     95
     96
     97
     98
     99
    100
    101
    102
    103
    104
    105
    106
    107
    108
    109
    110
    111
    112
    113
    114
    115
    116
    117
    118
    119
    120
    121
    122
    123
    124
    125
    /**
     * PCMSolver, an API for the Polarizable Continuum Model
     * Copyright (C) 2016 Roberto Di Remigio, Luca Frediani and collaborators.
     *
     * This file is part of PCMSolver.
     *
     * PCMSolver is free software: you can redistribute it and/or modify
     * it under the terms of the GNU Lesser General Public License as published by
     * the Free Software Foundation, either version 3 of the License, or
     * (at your option) any later version.
     *
     * PCMSolver is distributed in the hope that it will be useful,
     * but WITHOUT ANY WARRANTY; without even the implied warranty of
     * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
     * GNU Lesser General Public License for more details.
     *
     * You should have received a copy of the GNU Lesser General Public License
     * along with PCMSolver.  If not, see <http://www.gnu.org/licenses/>.
     *
     * For information on the complete list of contributors to the
     * PCMSolver API, see: <http://pcmsolver.readthedocs.io/>
     */
    
    #include "catch.hpp"
    
    #include <vector>
    #include <cmath>
    
    #include <Eigen/Core>
    
    #include "cavity/GePolCavity.hpp"
    #include "TestingMolecules.hpp"
    
    SCENARIO("GePol cavity for a single sphere", "[gepol][gepol_point]") {
      GIVEN("A single sphere") {
        double area = 0.4;
        double probeRadius = 0.0;
        double minRadius = 100.0;
        WHEN("the sphere is obtained from a Molecule object") {
          Molecule point = dummy<0>();
          GePolCavity cavity = GePolCavity(point, area, probeRadius, minRadius, "point");
          cavity.saveCavity("point.npz");
    
          /*! \class GePolCavity
           *  \test \b GePolCavityTest_size tests GePol cavity size for a point charge in
           * C1 symmetry without added spheres
           */
          THEN("the size of the cavity is") {
            int size = 32;
            int actualSize = cavity.size();
            REQUIRE(size == actualSize);
          }
          /*! \class GePolCavity
           *  \test \b GePolCavityTest_area tests GePol cavity surface area for a point
           * charge in C1 symmetry without added spheres
           */
          AND_THEN("the surface area of the cavity is") {
            double area = 4.0 * M_PI * pow(1.0, 2);
            double actualArea = cavity.elementArea().sum();
            REQUIRE(area == Approx(actualArea));
          }
          /*! \class GePolCavity
           *  \test \b GePolCavityTest_volume tests GePol cavity volume for a point
           * charge in C1 symmetry without added spheres
           */
          AND_THEN("the volume of the cavity is") {
            double volume = 4.0 * M_PI * pow(1.0, 3) / 3.0;
            Eigen::Matrix3Xd elementCenter = cavity.elementCenter();
            Eigen::Matrix3Xd elementNormal = cavity.elementNormal();
            double actualVolume = 0;
            for (int i = 0; i < cavity.size(); ++i) {
              actualVolume +=
                  cavity.elementArea(i) * elementCenter.col(i).dot(elementNormal.col(i));
            }
            actualVolume /= 3;
            REQUIRE(volume == Approx(actualVolume));
          }
        }
      }
    
      GIVEN("A single sphere") {
        double area = 0.4;
        double probeRadius = 0.0;
        double minRadius = 100.0;
        WHEN("the sphere is obtained from a Sphere object") {
          Sphere sph(Eigen::Vector3d::Zero(), 1.0);
          GePolCavity cavity = GePolCavity(sph, area, probeRadius, minRadius, "point");
    
          /*! \class GePolCavity
           *  \test \b GePolCavitySphereCTORTest_size tests GePol cavity size for a point
           * charge in C1 symmetry without added spheres
           */
          THEN("the size of the cavity is") {
            int size = 32;
            int actualSize = cavity.size();
            REQUIRE(size == actualSize);
          }
          /*! \class GePolCavity
           *  \test \b GePolCavitySphereCTORTest_area tests GePol cavity surface area for
           * a point charge in C1 symmetry without added spheres
           */
          AND_THEN("the surface area of the cavity is") {
            double area = 4.0 * M_PI * pow(1.0, 2);
            double actualArea = cavity.elementArea().sum();
            REQUIRE(area == Approx(actualArea));
          }
          /*! \class GePolCavity
           *  \test \b GePolCavitySphereCTORTest_volume tests GePol cavity volume for a
           * point charge in C1 symmetry without added spheres
           */
          AND_THEN("the volume of the cavity is") {
            double volume = 4.0 * M_PI * pow(1.0, 3) / 3.0;
            Eigen::Matrix3Xd elementCenter = cavity.elementCenter();
            Eigen::Matrix3Xd elementNormal = cavity.elementNormal();
            double actualVolume = 0;
            for (int i = 0; i < cavity.size(); ++i) {
              actualVolume +=
                  cavity.elementArea(i) * elementCenter.col(i).dot(elementNormal.col(i));
            }
            actualVolume /= 3;
            REQUIRE(volume == Approx(actualVolume));
          }
        }
      }
    }
    

    In this example we are creating a test fixture. The fixture will instatiate a GePolCavity with fixed parameters. The result is then tested against reference values in the various SECTIONs. It is important to add the documentation lines on top of the tests, to help other developers understand which class is being tested and what parameters are being tested. Within Catch fixtures are created behind the curtains, you do not need to worry about those details. This results in somewhat terser test source files.

Timer class

The Timer class enables timing of execution throughout the module. Timer support is enabled by passing -DENABLE_TIMER=ON to the setup.py script. Timing macros are available by inclusion of the Config.hpp header file.

The class is basically a wrapper around an ordered map of strings and cpu timers. To time a code snippet:

TIMER_ON("code-snippet");
// code-snippet
TIMER_OFF("code-snippet");

The timings are printed out to the pcmsolver.timer.dat by a call to the TIMER_DONE macro. This should obviously happen at the very end of the execution!

Defines

TIMER_ON(...)

PCMSolver, an API for the Polarizable Continuum Model Copyright (C) 2017 Roberto Di Remigio, Luca Frediani and collaborators.

This file is part of PCMSolver.

PCMSolver is free software: you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.

PCMSolver is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details.

You should have received a copy of the GNU Lesser General Public License along with PCMSolver. If not, see http://www.gnu.org/licenses/.

For information on the complete list of contributors to the PCMSolver API, see: http://pcmsolver.readthedocs.io/

TIMER_OFF(...)
TIMER_DONE(...)

Classes and functions reference

Cavities

We will here describe the inheritance hierarchy for generating cavities, in order to use and extend it properly. The runtime creation of cavity objects relies on the Factory Method pattern [GHJV94][Ale01], implemented through the generic Factory class.

_images/cavity.png

ICavity

class ICavity

Abstract Base Class for cavities.

This class represents a cavity made of spheres, its surface being discretized in terms of finite elements.

Author
Krzysztof Mozgawa
Date
2011

GePolCavity

class GePolCavity

A class for GePol cavity.

This class is an interface to the Fortran code PEDRA for the generation of cavities according to the GePol algorithm.

Author
Krzysztof Mozgawa, Roberto Di Remigio
Date
2011, 2016

PCMSolver, an API for the Polarizable Continuum Model Copyright (C) 2016 Roberto Di Remigio, Luca Frediani and collaborators.

This file is part of PCMSolver.

PCMSolver is free software: you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.

PCMSolver is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details.

You should have received a copy of the GNU Lesser General Public License along with PCMSolver. If not, see http://www.gnu.org/licenses/.

For information on the complete list of contributors to the PCMSolver API, see: http://pcmsolver.readthedocs.io/

RestartCavity

class RestartCavity

A class for Restart cavity.

Author
Roberto Di Remigio
Date
2014

Green’s Functions

We will here describe the inheritance hierarchy for generating Green’s functions, in order to use and extend it properly. The runtime creation of Green’s functions objects relies on the Factory Method pattern [GHJV94][Ale01], implemented through the generic Factory class.

_images/green.png

IGreensFunction

class IGreensFunction

Interface for Green’s function classes.

We define as Green’s function a function:

\[ G(\mathbf{r}, \mathbf{r}^\prime) : \mathbb{R}^6 \rightarrow \mathbb{R} \]
Green’s functions and their directional derivatives appear as kernels of the \(\mathcal{S}\) and \(\mathcal{D}\) integral operators. Forming the matrix representation of these operators requires performing integrations over surface finite elements. Since these Green’s functions present a Coulombic divergence, the diagonal elements of the operators will diverge unless appropriately formulated. This is possible, but requires explicit access to the subtype of this abstract base object. This justifies the need for the singleLayer and doubleLayer functions. The code uses the Non-Virtual Interface (NVI) idiom.
Author
Luca Frediani and Roberto Di Remigio
Date
2012-2016

GreensFunction

class GreensFunction

Templated interface for Green’s functions.

Author
Luca Frediani and Roberto Di Remigio
Date
2012-2016
Template Parameters
  • DerivativeTraits: evaluation strategy for the function and its derivatives
  • ProfilePolicy: dielectric profile type

Vacuum

class Vacuum

Green’s function for vacuum.

Author
Luca Frediani and Roberto Di Remigio
Date
2012-2016
Template Parameters
  • DerivativeTraits: evaluation strategy for the function and its derivatives

UniformDielectric

class UniformDielectric

Green’s function for uniform dielectric.

Author
Luca Frediani and Roberto Di Remigio
Date
2012-2016
Template Parameters
  • DerivativeTraits: evaluation strategy for the function and its derivatives

IonicLiquid

class IonicLiquid

Green’s functions for ionic liquid, described by the linearized Poisson-Boltzmann equation.

Author
Luca Frediani, Roberto Di Remigio
Date
2013-2016
Template Parameters
  • DerivativeTraits: evaluation strategy for the function and its derivatives

AnisotropicLiquid

class AnisotropicLiquid

Green’s functions for anisotropic liquid, described by a tensorial permittivity.

Author
Roberto Di Remigio
Date
2016
Template Parameters
  • DerivativeTraits: evaluation strategy for the function and its derivatives

SphericalDiffuse

class SphericalDiffuse

Green’s function for a diffuse interface with spherical symmetry.

This class is general, in the sense that no specific dielectric profile has been set in its definition. In principle any profile that can be described by:

  1. a left-side dielectric constant;
  2. a right-side dielectric constant;
  3. an interface layer width;
  4. an interface layer center can be used to define a new diffuse interface with spherical symmetry. The origin of the dielectric sphere can be changed by means of the constructor. The solution of the differential equation defining the Green’s function is always** performed assuming that the dielectric sphere is centered in the origin of the coordinate system. Whenever the public methods are invoked to “sample” the Green’s function at a pair of points, a translation of the sampling points is performed first.
Author
Hui Cao, Ville Weijo, Luca Frediani and Roberto Di Remigio
Date
2010-2015
Template Parameters
  • ProfilePolicy: functional form of the diffuse layer

Solvers

We will here describe the inheritance hierarchy for generating solvers, in order to use and extend it properly. The runtime creation of solver objects relies on the Factory Method pattern [GHJV94][Ale01], implemented through the generic Factory class.

_images/solver.png

ISolver

class ISolver

Abstract Base Class for solvers inheritance hierarchy.

We use the Non-Virtual Interface idiom.

Author
Luca Frediani, Roberto Di Remigio
Date
2011, 2015, 2016

IEFSolver

class IEFSolver

IEFPCM, collocation-based solver.

Author
Luca Frediani, Roberto Di Remigio
Date
2011, 2015, 2016
Note
We store the non-Hermitian, symmetry-blocked T(epsilon) and Rinfinity matrices. The ASC is obtained by multiplying the MEP by Rinfinity and then using a partially pivoted LU decomposition of T(epsilon) on the resulting vector. In case the polarization weights are requested, we use the approach suggested in [2]. First, the adjoint problem is solved:
\[ \mathbf{T}_\varepsilon^\dagger \tilde{v} = v \]
Also in this case a partially pivoted LU decomposition is used. The “transposed” ASC is obtained by the matrix-vector multiplication:
\[ q^* = \mathbf{R}_\infty^\dagger \tilde{v} \]
Eventually, the two sets of charges are summed and divided by 2 This avoids computing and storing the inverse explicitly, at the expense of storing both T(epsilon) and Rinfinity.

CPCMSolver

class CPCMSolver

Solver for conductor-like approximation: C-PCM (COSMO)

Author
Roberto Di Remigio
Date
2013, 2016
Note
We store the scaled, Hermitian, symmetrized S matrix and use a robust Cholesky decomposition to solve for the ASC. This avoids computing and storing the inverse explicitly. The S matrix is already scaled by the dielectric factor entering the definition of the conductor model!

Helper classes and functions

_images/utils.png

Input

Warning

doxygenclass: Cannot find class “Input” in doxygen xml output for project “PCMSolver” from directory: xml

Sphere

struct Sphere

POD describing a sphere.

Author
Roberto Di Remigio
Date
2011, 2016

Atom

struct Atom

A POD describing an atom.

Author
Roberto Di Remigio
Date
2011, 2016

Molecule

Warning

doxygenclass: Cannot find class “Molecule” in doxygen xml output for project “PCMSolver” from directory: xml

Solvent

struct Solvent

POD describing a solvent.

A Solvent object contains all the solvent-related experimental data needed to set up the Green’s functions and the non-electrostatic terms calculations.

Author
Roberto Di Remigio
Date
2011, 2016

Symmetry

class Symmetry

Contains very basic info about symmetry (only Abelian groups)

Just a wrapper around a vector containing the generators of the group

Author
Roberto Di Remigio
Date
2014

Public Functions

Symmetry()

Default constructor sets up C1 point group.

Private Members

int nrGenerators_

Number of generators

int generators_[3]

Generators

int nrIrrep_

Number of irreps

Mathematical utilities

namespace pcm

PCMSolver, an API for the Polarizable Continuum Model Copyright (C) 2017 Roberto Di Remigio, Luca Frediani and collaborators.

This file is part of PCMSolver.

PCMSolver is free software: you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.

PCMSolver is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details.

You should have received a copy of the GNU Lesser General Public License along with PCMSolver. If not, see http://www.gnu.org/licenses/.

For information on the complete list of contributors to the PCMSolver API, see: http://pcmsolver.readthedocs.io/

namespace utils

Functions

template <size_t nBits>
int parity(std::bitset<nBits> bitrep)

Calculate the parity of the bitset as defined by: bitrep[0] XOR bitrep[1] XOR ... XOR bitrep[nBits-1]

Parameters
  • bitrep: a bitset
Template Parameters
  • nBits: lenght of the input bitset

double parity(unsigned int i)

Returns parity of input integer. The parity is defined as the result of using XOR on the bitrep of the given integer. For example: 2 -> 010 -> 0^1^0 = 1 -> -1.0 6 -> 110 -> 1^1^0 = 0 -> 1.0

Parameters
  • i: an integer, usually an index for an irrep or a symmetry operation

It can also be interpreted as the action of a given operation on the Cartesian axes: zyx Parity 0 000 E 1.0 1 001 Oyz -1.0 2 010 Oxz -1.0 3 011 C2z 1.0 4 100 Oxy -1.0 5 101 C2y 1.0 6 110 C2x 1.0 7 111 i -1.0

bool isZero(double value, double threshold)

Returns true if value is less or equal to threshold

Parameters
  • value: the value to be checked
  • threshold: the threshold

bool numericalZero(double value)

Returns true if value is less than 1.0e-14

Parameters
  • value: the value to be checked

void symmetryBlocking(Eigen::MatrixXd &matrix, PCMSolverIndex cavitySize, PCMSolverIndex ntsirr, int nr_irrep)
void symmetryPacking(std::vector<Eigen::MatrixXd> &blockedMatrix, const Eigen::MatrixXd &fullMatrix, int dimBlock, int nrBlocks)

Parameters
  • blockedMatrix: the result of packing fullMatrix
  • fullMatrix: the matrix to be packed
  • dimBlock: the dimension of the square blocks
  • nrBlocks: the number of square blocks

template <typename Derived>
void hermitivitize(Eigen::MatrixBase<Derived> &obj_)

Given obj_ returns 0.5 * (obj_ + obj_^dagger)

Note
We check if a matrix or vector was given, since in the latter case we only want the complex conjugation operation to happen.
Parameters
  • obj_: the Eigen object to be hermitivitized
Template Parameters
  • Derived: the numeric type of obj_ elements

void eulerRotation(Eigen::Matrix3d &R_, const Eigen::Vector3d &eulerAngles_)

Build rotation matrix between two reference frames given the Euler angles.

We assume the convention \( R = Z_3 X_2 Z_1 \) for the ordering of the extrinsic elemental rotations (see http://en.wikipedia.org/wiki/Euler_angles) The Euler angles are given in the order \( \phi, \theta, \psi \). If we write \( c_i, s_i \,\, i = 1, 3 \) for their cosines and sines the rotation matrix will be:

\[\begin{split} R = \begin{pmatrix} c_1c_3 - s_1c_2s_3 & -s_1c_3 - c_1c_2s_3 & s_2s_3 \\ c_1s_3 + s_1c_2c_3 & -s_1s_3 + c_1c_2c_3 & -s_2c_3 \\ s_1s_2 & c_1s_2 & c_2 \end{pmatrix} \end{split}\]
Eigen’s geometry module is used to calculate the rotation matrix
Parameters
  • R_: the rotation matrix
  • eulerAngles_: the Euler angles, in degrees, describing the rotation

double linearInterpolation(const double point, const std::vector<double> &grid, const std::vector<double> &function)

Return value of function defined on grid at an arbitrary point.

This function finds the nearest values for the given point and performs a linear interpolation.

Warning
This function assumes that grid has already been sorted!
Parameters
  • point: where the function has to be evaluated
  • grid: holds points on grid where function is known
  • function: holds known function values

double splineInterpolation(const double point, const std::vector<double> &grid, const std::vector<double> &function)

Return value of function defined on grid at an arbitrary point.

This function finds the nearest values for the given point and performs a cubic spline interpolation.

Warning
This function assumes that grid has already been sorted!
Parameters
  • point: where the function has to be evaluated
  • grid: holds points on grid where function is known
  • function: holds known function values

template <typename Derived>
void print_eigen_matrix(const Eigen::MatrixBase<Derived> &matrix, const std::string &fname)

Prints Eigen object (matrix or vector) to file.

Note
This is for debugging only, the format is in fact rather ugly. Row index Column index Matrix entry 0 0 0.0000
Parameters
  • matrix: Eigen object
  • fname: name of the file
Template Parameters
  • Derived: template parameters of the MatrixBase object

namespace cnpy
namespace custom

Custom overloads for cnpy load and save functions

Functions

template <typename Scalar, int Rows, int Cols>
void npy_save(const std::string &fname, const Eigen::Matrix<Scalar, Rows, Cols> &obj)

Save Eigen object to NumPy array file.

Parameters
  • fname: name of the NumPy array file
  • obj: Eigen object to be saved, either a matrix or a vector
Template Parameters
  • Scalar: the data type of the matrix to be returned. Default is double
  • Rows: number of rows in the Eigen object. Default is dynamic e
  • Cols: number of columns in the Eigen object. Default is dynamic

template <typename Scalar, int Rows, int Cols>
void npz_save(const std::string &fname, const std::string &name, const Eigen::Matrix<Scalar, Rows, Cols> &obj, bool overwrite = false)

Save Eigen object to a compressed NumPy file.

Parameters
  • fname: name of the compressed NumPy file
  • name: tag for the given object in the compressed NumPy file
  • obj: Eigen object to be saved, either a matrix or a vector
  • overwrite: if file exists, overwrite. Appends by default.
Template Parameters
  • Scalar: the data type of the matrix to be returned. Default is double
  • Rows: number of rows in the Eigen object. Default is dynamic
  • Cols: number of columns in the Eigen object. Default is dynamic

template <typename Scalar>
Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> npy_to_eigen(const NpyArray &npy_array)

Load NpyArray object into Eigen object.

Return
An Eigen object (matrix or vector) with the data
Warning
We check that the rank of the object read is not more than 2 Eigen cannot handle general tensors.
Parameters
  • npy_array: the NpyArray object
Template Parameters
  • Scalar: the data type of the matrix to be returned. Default is double

template <typename Scalar>
Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> npy_load(const std::string &fname)

Load NumPy array file into Eigen object.

Return
An Eigen object (matrix or vector) with the data
Parameters
  • fname: name of the NumPy array file
Template Parameters
  • Scalar: the data type of the matrix to be returned. Default is double

Namespaces

We use namespaces to delimit the visibility of functions and classes defined in the various subdirectories of the project. Namespaces provide a convenient layered structure to the project and we use them as a convention to signal which functions and classes are supposed to be used in any given layer. The top-level namespace is called pcm and includes all functions and classes that can be called from the outside world, i.e. a C++ API. Each subdirectory introduces a new namespace of the same name, nested into pcm. Code that can be used outside of a given subdirectory is put directly in the pcm namespace, i.e. the outermost layer. Finally, the namespace detail, at the third level of nesting, is used for functions and classes that are used exclusively within the code in a given subdirectory.

References

[Ale01]Andrei Alexandrescu. Modern C++ design: generic programming and design patterns applied. Addison-Wesley Longman Publishing Co., Inc., 2001. ISBN 0-201-70431-5.
[AZB94]Norman L Allinger, Xuefeng Zhou, and John Bergsma. Molecular mechanics parameters. Journal of Molecular Structure: THEOCHEM, 312(1):69–83, 1~January 1994. doi:10.1016/S0166-1280(09)80008-0.
[Cli]Marshall P. Cline. C++ FAQ. URL: http://www.parashift.com/c++-faq.
[CGL98]Marshall P. Cline, Mike Girou, and Greg Lomow. C++ FAQs. Addison-Wesley Longman Publishing Co., Inc., 1998. ISBN 0201309831.
[GHJV94]Erich Gamma, Richard Helm, Ralph Johnson, and John Vlissides. Design patterns: elements of reusable object-oriented software. Addison-Wesley Longman Publishing Co., Inc., 1994. ISBN 0-201-63361-2.
[Kop08]Joachim Kopp. Efficient numerical diagonalization of hermitian 3x3~matrices. Int. J. Mod. Phys. C, 19(03):523–548, 2008. doi:10.1142/S0129183108012303.
[RCC+92]A. K. Rappe, C. J. Casewit, K. S. Colwell, W. A. Goddard, and W. M. Skiff. UFF, a full periodic table force field for molecular mechanics and molecular dynamics simulations. J. Am. Chem. Soc., 114(25):10024–10035, 1992. URL: http://pubs.acs.org/doi/abs/10.1021/ja00051a040, doi:10.1021/ja00051a040.
[Sut99]Herb Sutter. Exceptional C++: 47 engineering puzzles, programming problems, and solutions. Addison-Wesley Longman Publishing Co., Inc., 1999. ISBN 0-201-61562-2.
[SA04]Herb Sutter and Andrei Alexandrescu. C++ Coding Standards: 101 Rules, Guidelines, and Best Practices (C++ in Depth Series). Addison-Wesley Professional, 2004. ISBN 0321113586.
[TMC05]Jacopo Tomasi, Benedetta Mennucci, and Roberto Cammi. Quantum mechanical continuum solvation models. Chem. Rev., 105(8):2999–3093, 2005. doi:10.1021/cr9904009.
[WAB+14]Greg Wilson, D a Aruliah, C Titus Brown, Neil P Chue Hong, Matt Davis, Richard T Guy, Steven H D Haddock, Kathryn D Huff, Ian M Mitchell, Mark D Plumbley, Ben Waugh, Ethan P White, and Paul Wilson. Best practices for scientific computing. PLoS Biol., 12(1):e1001745, 2014. URL: http://www.pubmedcentral.nih.gov/articlerender.fcgi?artid=3886731&tool=pmcentrez&rendertype=abstract, doi:10.1371/journal.pbio.1001745.
[Bondi64]A. Bondi. van der Waals Volumes and Radii. J. Phys. Chem., 68(3):441–451, 1964. URL: http://pubs.acs.org/doi/pdf/10.1021/j100785a001, doi:10.1021/j100785a001.
[CancesMennucci98]Eric Cancès and Benedetta Mennucci. New Applications of Integral Equations Methods for Solvation Continuum Models: Ionic Solutions and Liquid Crystals. J. Math. Chem., 23:309–326, 1998. doi:10.1023/A:1019133611148.
[MantinaChamberlinValero+09]Manjeera Mantina, Adam C. Chamberlin, Rosendo Valero, Christopher J. Cramer, and Donald G. Truhlar. Consistent van der Waals Radii for the Whole Main Group. J. Phys. Chem. A, 113:5806–5812, 2009. URL: http://pubs.acs.org/doi/pdf/10.1021/jp8111556, doi:10.1021/jp8111556.

Indices and tables