The Python-based Simulations of Chemistry Framework (PySCF)
Qiming Sun Timothy C. Berkelbach Nick S. Blunt George H. Booth Sheng Guo Zhendong Li Junzi Liu James McClain Elvira R. Sayfutyarova Sandeep Sharma Sebastian Wouters Garnet Kin-Lic Chan


1]Qiming Sunosirpt.sun@gmail.com
2]Timothy C. Berkelbach
3,4]Nick S. Blunt
5]George H. Booth
1,6]Sheng Guo
1]Zhendong Li
7]Junzi Liu
1,6]James D. McClain
1,6]Elvira R. Sayfutyarova
8]Sandeep Sharma
9]Sebastian Wouters
1]Garnet Kin-Lic Changkc1000@gmail.com
[1]Division of Chemistry and Chemical Engineering, California Institute of Technology, Pasadena CA 91125, USA
[2]Department of Chemistry and James Franck Institute, University of Chicago, Chicago, Illinois 60637, USA
[3]Chemical Science Division, Lawrence Berkeley National Laboratory, Berkeley, California 94720, USA
[4]Department of Chemistry, University of California, Berkeley, California 94720, USA
[5]Department of Physics, King's College London, Strand, London WC2R 2LS, United Kingdom
[6]Department of Chemistry, Princeton University, Princeton, New Jersey 08544, USA
[7]Institute of Chemistry Chinese Academy of Sciences, Beijing 100190, P. R. China
[8]Department of Chemistry and Biochemistry, University of Colorado Boulder, Boulder, CO 80302, USA
[9]Brantsandpatents, Pauline van Pottelsberghelaan 24, 9051 Sint-Denijs-Westrem, Belgium




[NO \title GIVEN]
    [
    
=================








  §.§.§ Abstract








PySCF is a general-purpose electronic structure platform designed from the ground
up to emphasize code simplicity, so as to facilitate new method development and enable
flexible computational workflows.  The package provides a wide range of tools to
support simulations of finite-size systems, extended systems with periodic boundary
conditions, low-dimensional periodic systems, and custom Hamiltonians, using mean-field
and post-mean-field methods with standard Gaussian basis functions.  To ensure ease of
extensibility, PySCF uses the Python language to implement almost all of its
features, while computationally critical paths are implemented with heavily optimized C
routines.  Using this combined Python/C implementation, the package is as efficient as the
best existing C or Fortran based quantum chemistry programs.  In this paper we document
the capabilities and design philosophy of the current version of the PySCF
package.








iblabel[1]#1.
apsrev









§ INTRODUCTION


The Python programming language is playing an increasingly important role in scientific
computing.  As a high level language, Python supports rapid development practices and easy
program maintenance.  While programming productivity is hard to measure, it is commonly
thought that it is more efficient to prototype new ideas in Python, rather than in
traditional low-level compiled languages such as Fortran or C/C++. 
Further, through the use of the many high-quality numerical libraries available in Python
– such as NumPy<cit.>, SciPy<cit.>, and MPI4Py<cit.>
– Python programs can perform at competitive levels with optimized Fortran and C/C++
programs, including on large-scale computing architectures.






There have been several efforts in the past to incorporate Python into electronic structure
programs.
Python has been widely adopted in a scripting role: 
the Psi4<cit.> quantum chemistry package uses a custom “Psithon” dialect to
drive the underlying C++ implementation, while general simulation environments
such as ASE<cit.> and PyMatGen<cit.> provide Python frontends to multiple
quantum chemistry and electronic structure packages, to organize complex workflows<cit.>. Python
has also proved popular for implementing symbolic second-quantized algebra and code
generation tools, such as the Tensor Contraction Engine<cit.> and the
SecondQuantizationAlgebra library<cit.>.

In the above cases, Python has been used as a supporting language, with the underlying
quantum chemistry algorithms implemented in a compiled language.  However, Python has also
seen some use as a primary implementation language for electronic structure methods.
PyQuante<cit.> was an early attempt to implement a Gaussian-based
quantum chemistry code in Python, although it did not achieve speed or functionality
competitive with typical packages.  Another early effort was the GPAW<cit.>
code, which implements the projector augmented wave formalism for density functional
theory, and which is still under active development in multiple groups.  Nonetheless, it
is probably fair to say that using Python as an implementation language, rather than a
supporting language, remains the exception rather than the rule in modern quantum
chemistry and electronic structure software efforts.



















With the aim of developing a new highly functional, high-performance computing toolbox
for the quantum chemistry of molecules and materials implemented primarily in the
Python language, we started the open-source project “Python-based Simulations of Chemistry Framework” (PySCF) in 2014.
The program was initially ported from our quantum chemistry density matrix embedding 
theory (DMET) project<cit.> and contained only the Gaussian integral interface, a
basic Hartree-Fock solver, and a few post-Hartree-Fock components required by DMET.
In the next 18 months, multi-configurational self-consistent-field (MCSCF), density
functional theory and coupled cluster theory, as well as relevant modules for
molecular properties, were added into the package.
In 2015, we released the first stable version, PySCF 1.0, wherein we codified our primary
goals for further code development: to produce a package that emphasizes simplicity,
generality, and efficiency, in that order.

As a result of this choice, most functions in PySCF are written purely in Python, with a very 
limited amount of C code only for the most time-critical parts.
The various features and APIs are designed and implemented in the simplest and most
straightforward manner, so that users can easily modify the source code to meet their own
scientific needs and workflow.
However, although we have favored algorithm accessibility and extensibility over
performance, we have found that the efficient use of numerical Python libraries allows
PySCF to perform at least as fast as the best existing quantum chemistry
implementations. 
In this article, we highlight the current capabilities and design philosophy of the
PySCF package.




§ CAPABILITIES




Molecular electronic structure methods are typically the main focus of quantum chemistry
packages. We have put significant effort towards the production of a stable, feature-rich
and efficient molecular simulation environment in PySCF.
In addition to molecular quantum chemistry methods, PySCF also provides a wide
range of quantum chemistry methods for extended systems with periodic boundary conditions
(PBC).  Table <ref> lists the main electronic structure
methods available in the PySCF package.
More detailed descriptions are presented in Section <ref> - Section <ref>.
Although not listed in the table, many auxiliary tools for method development are also part of the package.
They are briefly documented in Section <ref> - Section <ref>.




 §.§ Self-consistent field methods



Self-consistent field (SCF) methods are the starting point for most electronic structure
calculations. In PySCF, the SCF module includes implementations of Hartree-Fock
(HF) and density functional theory (DFT) for restricted, unrestricted, closed-shell and
open-shell Slater determinant references.  A wide range of predefined exchange-correlation (XC) functionals
are supported through a general interface to the Libxc<cit.> and
Xcfun<cit.> functional libraries.  Using the interface, as shown
in Figure <ref>, one can easily customize the XC functionals in DFT calculations.
PySCF uses the Libcint<cit.> Gaussian
integral library, written by one of us (QS) as its integral engine.
In its current implementation, the SCF program can handle over 5000 basis functions on a
single symmetric multiprocessing (SMP) node without any approximations to the integrals.
To obtain rapid convergence in the SCF iterations, we have also developed a second order
co-iterative augmented Hessian (CIAH) algorithm for orbital optimization <cit.>.
Using the direct SCF technique with the CIAH algorithm, we are able to converge
a Hartree-Fock calculation for the open-shell molecule Fe(II)-porphine (2997
AOs) on a 16-core node in one day.




 §.§ Post SCF methods



Single-reference correlation methods can be used on top of the HF or DFT references,
including Møller-Plesset second-order perturbation theory (MP2), configuration
interaction, and coupled cluster theory. 

Canonical single-reference coupled cluster theory has been implemented with 
single and double excitations (CCSD)<cit.> and with perturbative triples
[CCSD(T)].
The associated derivative routines include CCSD and CCSD(T) density matrices, CCSD and
CCSD(T) analytic gradients, and equation-of-motion CCSD for the ionization potentials,
electron affinities, and excitation energies
(EOM-IP/EA/EE-CCSD)<cit.>.
The package contains two complementary implementations of each of these methods.
The first set are straightforward spin-orbital and spatial-orbital implementations, which
are optimized for readability and written in pure Python using syntax of the 
NumPy  function (which can use either the default Numpy implementation
or a custom -based version) for tensor contraction.
These implementations are easy for the user to modify.  A second
spatial-orbital implementation has been intensively optimized to minimize
dataflow and uses asynchronous I/O and a threaded  function for efficient
tensor contractions.  For a system of 25 occupied orbitals and 1500 virtual orbitals 
(H_50 with cc-pVQZ basis), the latter CCSD implementation takes less than 3 hours to
finish one iteration using 28 CPU cores. 

The configuration interaction code implements two solvers: a solver for configuration
interaction with single and double excitations (CISD), and a determinant-based full
configuration interaction (FCI) solver<cit.> for fermion, boson or coupled
fermion-boson Hamiltonians.  The CISD solver has a similar program layout
to the CCSD solver.  The FCI solver additionally implements the spin-squared operator,
second quantized creation and annihilation operators (from which arbitrary second
quantized algebra can be implemented), functions to evaluate the density matrices and
transition density matrices (up to fourth order), as well as a function to evaluate the
overlap of two FCI wavefunctions in different orbital bases.  The FCI solver is
intensively optimized for multi-threaded performance.  It can perform one matrix-vector
operation for 16 electrons and 16 orbitals using 16 CPU cores in 30 seconds.




 §.§ Multireference methods









For multireference problems, the PySCF package provides the complete active space self
consistent field (CASSCF) method<cit.> and N-electron
valence perturbation theory (NEVPT2)<cit.>.
When the size of the active space exceeds the capabilities of the

conventional FCI solver, one can switch to  external variational solvers
such as a density matrix renormalization group (DMRG) program<cit.>
or a full configuration interaction quantum Monte Carlo (FCIQMC) program<cit.>.
Incorporating  external solvers into the CASSCF optimizer widens
the range of possible applications, while raising new challenges for an efficient CASSCF







algorithm.
One challenge is the communication between the external solver and the orbital
optimization driver; communication must be limited to quantities
that are easy to obtain from the external solver.





A second challenge is the cost of handling quantities associated with the active space;
for example, as the active space becomes large, the memory required
to hold integrals involving active labels can easily exceed available memory.
Finally, any approximations introduced in the context of the above two challenges
should not interfere with the quality of convergence of the CASSCF optimizer.


To address these challenges,
we have implemented a general AO-driven CASSCF optimizer<cit.> that
provides second order convergence and which may easily be combined with a wide
variety of external variational solvers, including DMRG, FCIQMC and their
state-averaged solvers.
Only the 2-particle density matrix and Hamiltonian integrals
are communicated between the CASSCF driver and the external CI solver.
Further, the AO-driven algorithm has a low memory and I/O footprint. The current
implementation supports
calculations with 3000 basis functions and 30–50 active orbitals
on a single SMP node with 128 GB memory, without any approximations to the AO integrals.



A simple interface is provided to use an external solver in 
multiconfigurational calculations.
Figure <ref> shows how to perform a
DMRG-CASSCF calculation by replacing the  attribute of the CASSCF
method.

DMRG-SC-NEVPT2<cit.>, and ic-MPS-PT2 and
ic-MPS-LCC<cit.> methods are also available through the interface to the DMRG
program package Block<cit.>, and the ic-MPS-LCC program of





Sharma <cit.>.




 §.§ Molecular properties


At the present stage, the program can compute molecular properties 
such as analytic nuclear gradients, analytic nuclear Hessians, and NMR
shielding parameters at the SCF level.  The CCSD and CCSD(T) modules include solvers for
the Λ-equations.  As a result, we also provide one-particle and two-particle density
matrices, as well as the analytic nuclear gradients, for the CCSD and CCSD(T)
methods<cit.>.

For excited states, 
time-dependent HF (TDHF) and time-dependent DFT (TDDFT)
are implemented on top of the SCF module.  The relevant analytic nuclear
gradients are also programmed<cit.>.
The CCSD module offers another option to obtain
excited states using the EOM-IP/EA/EE-CCSD methods.  The third option
to obtain excited states
is through the multi-root CASCI/CASSCF solvers, optionally followed by the MRPT tool chain.  
Starting from the multi-root CASCI/CASSCF solutions, the program can compute
the density matrices of all the states and the transition density matrices
between any two states.  One can contract these density matrices with
specific AO integrals to obtain different first order molecular properties.




 §.§ Relativistic effects


Many different relativistic treatments are available in PySCF.
Scalar relativistic effects can be added to all SCF and post-SCF methods through relativistic
effective core potentials (ECP)<cit.> or the all-electron spin-free X2C<cit.>
relativistic correction.
For a more advanced treatment, PySCF also provides 4-component relativistic
Hartree-Fock and no-pair MP2 methods with Dirac-Coulomb, Dirac-Coulomb-Gaunt, and
Dirac-Coulomb-Breit Hamiltonians.
Although not programmed as a standalone module, the no-pair
CCSD electron correlation energy can also be computed with the straightforward
spin-orbital version of the CCSD program.
Using the 4-component Hamiltonian, molecular properties including analytic
nuclear gradients and NMR shielding parameters are available at the mean-field level<cit.>.






 §.§ Orbital localizer and result analysis


Two classes of orbital localization methods are available in the package.
The first emphasizes  the atomic character of the basis functions.
The relevant localization functions can generate intrinsic atomic orbitals
(IAO)<cit.>, natural atomic orbitals (NAO)<cit.>,
and meta-Löwdin orbitals<cit.> based on orbital projection and
orthogonalization.
With these AO-based local orbitals, charge distributions can be properly
assigned to atoms in  population analysis<cit.>.
In the PySCF population analysis code, meta-Löwdin orbitals are the
default choice.

The second class, represented by Boys-Foster, Edmiston-Ruedenberg, and
Pipek-Mezey localization, require minimizing (or maximizing) the dipole, the
Coulomb self-energy, or the atomic charges, to obtain the optimal localized orbitals.
The localization routines can take arbitrary orthogonal input orbitals and call
the CIAH algorithm to rapidly converge the solution.
For example, using 16 CPU cores, it takes 3 minutes to localize 1620 HF unoccupied
orbitals for the C_60 molecule using Boys localization.

A common task when analysing the results of an electronic structure calculation
is to visualize the orbitals.
Although PySCF does not have a visualization tool itself, it provides a module
to convert the given orbital coefficients to the <cit.> format which can
be read and visualized by other software, e.g. Jmol<cit.>.
Figure <ref> is an example to run Boys localization for the
C_60 HF occupied orbitals and to generate the orbital surfaces of the localized
σ-bond orbital in a single Python script.




 §.§ Extended systems with periodic boundary conditions


PBC implementations typically use either plane
waves<cit.> or local atomic
functions<cit.> as the underlying orbital basis.








The PBC implementation in PySCF uses the local basis formulation, specifically
crystalline orbital Gaussian basis functions ϕ, expanded in terms of a lattice sum
over local Gaussians χ

    ϕ_k,χ(r) = ∑_T
            e^ik·Tχ(r - T)

where k is a vector in the first Brillouin zone and T is a lattice
translational vector.
We use a pure Gaussian basis in our PBC implementation for two reasons: to 
simplify the development of post-mean-field methods for extended systems and to
have a seamless interface and direct comparability to
finite-sized quantum chemistry calculations.  Local bases are favourable for
post-mean-field methods because they are generally quite compact, resulting in small
virtual spaces <cit.>, and further allow locality to be exploited.  
Due to the use of local bases, various boundary conditions can be easily
applied in the PBC module, from  zero-dimensional systems (molecules) to
extended one-, two- and three-dimensional periodic systems.

The PBC module supports both all-electron and pseudopotential calculations.  Both
separable pseudopotentials (e.g. Goedecker-Teter-Hutter (GTH)
pseudopotentials <cit.>) and non-separable pseudopotentials (quantum chemistry ECPs and
Burkatzi-Filippi-Dolg pseudopotentials<cit.>) can be
used.  
In the separable pseudopotential implementation, the associated orbitals and densities are
guaranteed to be smooth, allowing a grid-based treatment that uses discrete fast Fourier
transforms <cit.>.
In both the pseudopotential and all-electron PBC calculations, 
Coulomb-based integrals are handled via density fitting as described in Section <ref>.


























The PBC implementation is organized in direct correspondence to the molecular implementation.
We implemented the same function interfaces as in the molecular code,
with analogous module and function names.
Consequently, methods defined in the molecular part of the code can be seamlessly mixed 
with the PBC functions without modification, especially in Γ-point calculations 
where the PBC wave functions are real.
Thus, starting from PBC Γ-point mean-field orbitals, one can,
for example, carry out CCSD, CASSCF, TDDFT, etc. calculations using
the molecular implementations.
We also introduce specializations of the PBC methods to support
k-point (Brillouin zone) sampling. 
The k-point methods slightly
modify the Γ-point data structures, but inherit from and reuse
almost all of the Γ-point functionality.
Explicit k-point sampling is supported at the HF and DFT level, and on top of this we
have also implemented k-point MP2, CCSD, CCSD(T) and EOM-CCSD methods<cit.>, with
optimizations to carefully distribute work and data across cores.





On 100 computational cores, mean-field
simulations including unit cells with over 100 atoms, or k-point CCSD calculations with
over 3000 orbitals, can be executed without difficulty.




 §.§ General AO integral evaluator and J/K builds


Integral evaluation forms the foundation of Gaussian-based electronic
structure simulation.

The general integral evaluator library Libcint supports a wide range of
GTO integrals, and PySCF exposes simple APIs to access the Libcint integral
functions.
As the examples in Figure <ref> show, 
the PySCF integral API allows users to access AO integrals either in a
giant array or in individual shells with a single line of Python code. The
integrals provided include,

  
  * integrals in the basis of Cartesian, real-spherical and j-adapted spinor GTOs;
  
  * arbitrary integral expressions built from r, p, and
    σ polynomials;



  
  * 2-center, 3-center and 4-center 2-electron integrals for the Coulomb
    operator 1/r_12,  range-separated Coulomb operator
    erf(ω r_12)/r_12,
    Gaunt interaction, and Breit interaction.
    
    


Using the general AO integral evaluator, the package provides a general
AO-driven J/K contraction function.
J/K-matrix construction involves a contraction over a high order tensor
(e.g. 4-index 2-electron integrals (ij|kl)) and a low order tensor (e.g. the 2-index
density matrix γ)

    J_ij = ∑_kl (ij|kl) γ_kl
    
      K_il = ∑_jk (ij|kl) γ_jk

When both tensors can be held in memory, the Numpy package offers a convenient
tensor contraction function  
to quickly construct J/K matrices. 
However, it is common for the high order tensor to be too large to fit into the
available memory.
Using the Einstein summation notation of the Numpy  function,
our AO-driven J/K contraction implementation offers the capability to contract the
high order tensor (e.g. 2-electron integrals or their high order derivatives) with
multiple density matrices, with a small memory footprint.
The J/K contraction function also supports subsystem contraction, in
which the 4 indices of the 2-electron integrals are distributed over different segments
of the system which may or may not overlap with each other.
This subsystem contraction is particularly useful in two scenarios:
in fragment-based methods, where the evaluation of Coulomb or exchange energies
involves integral contraction over different fragments, and
in parallel algorithms, where one partitions the J/K contraction into
small segments and distributes them to different computing nodes.











 §.§ General integral transformations


Integral transformations are another fundamental operation found in quantum
chemistry programs. A common kind of integral transformation is to transform the 4 indices of the
2-electron integrals by 4 sets of different orbitals.
To satisfy this need, we designed a general integral transformation function to
handle the arbitrary AO integrals provided by the Libcint library and arbitrary kinds of orbitals.

To reduce  disk usage, we use permutation symmetry over i and j, k and l in
(ij|kl) whenever  possible for real integrals.

Integral transformations involve high computational and I/O costs.  
A standard approach to reduce these costs involves precomputation to reduce integral costs and data
compression to increase I/O throughput.  However, we have not adopted such an
optimization strategy in our implementation because it is against the objective of simplicity
for the PySCF package.  In our implementation, initialization is not 
required for the general integral transformation function.  Similarly to the AO integral API, the
integral transformation can thus be launched with one line of Python code.  In the
integral data structure, we store the transformed integrals by chunks in the HDF5
format without compression.  This choice has two advantages.  First, it allows
for fast indexing and hyperslab selection for subblocks of the integral array.
Second, the integral data can be easily accessed by other program packages
without any overhead for parsing the integral storage protocol.













 §.§ Density fitting


The density fitting (DF) technique is implemented for both  finite-sized
systems and crystalline systems with periodic boundary conditions.

In  finite-sized systems, one can use DF to approximate the J/K matrix and the MO
integrals for the HF, DFT and MP2 methods.  To improve the performance of
the CIAH algorithm, one can use the DF orbital Hessian in the CIAH
orbital optimization for Edmiston-Ruedenberg localization and for the HF, DFT and
CASSCF algorithms.

In the PBC module, the 2-electron integrals are represented as the product of
two 3-index tensors which are treated as DF objects.
Based on the requirements of the system being modelled, we have developed various DF 
representations.
When the calculation involves only smooth bases (typically with pseudopotentials),
plane waves are used as the auxiliary fitting functions and the DF 3-index
tensor is computed within a grid-based treatment using discrete fast Fourier transforms <cit.>.
When high accuracy in all-electron calculations is required,
a mixed density fitting technique is invoked in which the fitting functions are
Gaussian functions plus plane waves.
Besides the choice of fitting basis, different metrics (e.g. overlap, kinetic, or Coulomb) can be used
in the fitting to balance performance against computational accuracy.



The 3-index DF tensor is stored as a giant array in the HDF5 format without compression.
With this design, it is straightforward to access the 2-electron integrals
with the functions of the PySCF package.
Moreover, it allows us to supply 2-electron integrals to calculations by
overloading the DF object in cases where direct storage of the 4-index integrals in memory or on disk is infeasible
(see discussion in Section <ref>).




 §.§ Custom Hamiltonians


Most quantum chemistry approximations are not tied to the details of the ab initio molecular or periodic Hamiltonian.
This means that they can also be used with arbitrary model Hamiltonians,
which is of interest for semi-empirical quantum chemistry calculations
as well as condensed-matter  model studies.







In PySCF, overwriting the predefined Hamiltonian is straightforward.
The Hamiltonian is an attribute of the mean-field calculation object.
Once the 1-particle and 2-particle integral attributes of the mean-field object are
defined, they are used by the mean-field calculation and all subsequent
post-Hartree-Fock correlation treatments.
Users can thus carry out correlated calculations with model Hamiltonians in exactly the same way as
with standard ab initio Hamiltonians.
Figure <ref> displays an example of how to input a model Hamiltonian.





 §.§ Interfaces to external programs


PySCF can be used either as the driver to execute external programs or as
an independent solver to use as part of a computational workflow involving other software.
In PySCF, the DMRG programs Block<cit.> and
CheMPS2<cit.> and the FCIQMC program NECI<cit.> can be








used as a replacement for the FCI routine for large active spaces in the
CASCI/CASSCF solver.





In the QM/MM interface, by supplying the charges and the positions of the MM
atoms, one can compute the HF, DFT, MP2, CC, CI and MCSCF energies and their
analytic nuclear gradients.

To communicate with other quantum chemistry programs,
we provide utility functions to read and write Hamiltonians in the
Molpro<cit.>  format, and arbitrary orbitals in the
<cit.> format.
The program also supports to write SCF solution and CI wavefunction in the
GAMESS<cit.>  format and to read orbitals from
Molpro XML output.
The real space electron density can be output on cubic grids in the
Gaussian<cit.>  format.




 §.§ Numerical tools


Although the Numpy and Scipy libraries provide a wide range of
numerical tools for scientific computing, there are some numerical components
commonly found in quantum chemistry algorithms that are not provided by these
libraries.  For example, the direct inversion of the iterative
space (DIIS) method<cit.> is one of the most commonly used 
tools in quantum chemistry
to speed up optimizations when a second order algorithm is not
available.
In  PySCF we provide a general DIIS handler for an object array of
arbitrary size and arbitrary data type.  In the current implementation,
it supports  DIIS optimization both with or without supplying the error vectors.
For the latter case, the differences between the arrays of adjacent iterations are minimized.
Large scale eigenvalue and linear equation solvers are also
common components of many quantum chemistry methods.
The Davidson diagonalization algorithm and Arnoldi/Krylov subspace solver
are accessible in PySCF through simple APIs.




§ DESIGN AND IMPLEMENTATION OF PYSCF


While we have tried to provide rich functionality for quantum chemical simulations
with the built-in functions of the PySCF package, it will
nonetheless often be the case that a user's needs are not covered by the built-in functionality.
A major design goal has been to implement PySCF in a sufficiently flexible way so
that users can easily extend its functionality.  




To provide robust components for complex problems and non-trivial workflows,
we have made the following general design choices in
PySCF:

  
  * Language: Mostly Python, with a little C.  We believe that it is easiest to
  develop and test new functionality in Python.  For this reason, most functions in
  PySCF are written in pure Python.  Only a few computational hot spots have been
  rewritten and optimized in C.
  
  * Style: Mostly functional, with a little object-oriented programming (OOP).
  Although OOP is a successful and widely used programming paradigm, we feel that it is
  hard for users to customize typical OOP programs without learning details of the
  object hierarchy and interfaces.  We have adopted a functional programming style, where
  most functions are pure, and thus can be invoked alone and independently of each other.
  This allows users to mix functionality with a minimal knowledge of the
  PySCF internals.





We elaborate on these choices below.






 §.§ Input language


Almost every  quantum chemistry package today uses its own custom input language.
This is a burden to the user, who must become
familiar with a new domain-specific language for every new package.
In contrast, PySCF does not have an input language.  Rather, the functionality
is simply called from an input script written in the host Python language. 
This choice has clear benefits:

  
  * There is no need to learn a domain-specific language.
    Python, as a general programming language, is already widely used for
    numerical computing, and is taught in modern computer science courses.
    For novices, the language is easy to learn and help is readily available from the
    large Python community.
  
  * One can use all Python language features in the input script.
    This allows the input script to implement complex logic and computational
    workflows, and to carry out tasks (e.g. data processing and plotting) in the same
    script as the electronic structure simulation (see Figure <ref> for an
    example).
  
  * The computational environment is easily extended beyond that provided by the PySCF package.
    The PySCF package is a regular Python module
    which can be mixed and matched with other Python modules to build a personalized
     computing environment.
  
  * Computing can be carried out interactively.  Simulations can be tested, debugged, and executed step by step within
    the Python interpreter shell.












 §.§ Enabling interactive computing


As discussed above, a strength of the PySCF package is that its functionality
can be invoked from the interactive Python shell. However, maximizing its usability
in this interactive mode entails additional design optimizations.
There are three critical considerations to facilitate such interactive computations:

  
  * The functions and data need to be easy to access;
  
  * Functions should be insensitive to execution order (when and how many times
    a function is called should not affect the result);
  
  * Computations should not cause (significant) halts in the interactive shell.

To address these requirements, we have enforced the following design rules 
wherever possible in the package:

  
  * Functions are pure (i.e. state free). This ensures that they are insensitive to execution order;
  
  * Method objects (classes) only hold results and control parameters;
  
  * There is no initialization of functions, or at most a short initialization chain;
  
  * Methods are placed at both the module level and within classes so that
    the methods and their documentation can be easily accessed by the
    interactive shell (see Figure <ref>).

A practical solution to eliminate halting of the interactive shell is to overlap the REPL
(read-eval-print-loop) and task execution.  Such task parallelism requires
the underlying tasks to be independent of each other.  Although certain dependence between methods
is inevitable, the above design rules greatly reduce function call dependence.  Most functions in
PySCF can be safely placed in the background using the standard Python   and  libraries.




 §.§ Methods as plugins


Ease-of-use is the primary design objective of the PySCF package.
However, function simplicity and versatility are difficult to balance in
the same software framework.
To balance readability and complexity,
we have implemented only the basic algorithmic features in the main methods, and placed advanced
features in additional “plugins”.
For instance, the main mean-field module implements only the basic self-consistent loop.
Corrections (such as for relativistic effects) are implemented in an independent
plugin module, which can be activated by reassigning the mean-field 1-electron Hamiltonian
method at runtime.
Although this design increases the complexity of implementation of the plugin functions,
the core methods retain a clear structure  and are easy to comprehend.
Further, this approach decreases the coupling between different features: for example,
independent features can be modified and tested independently and combined in calculations.
In the package, this plugin design has been widely used, for example, to enable molecular point group symmetry,
relativistic corrections, solvation effects, density fitting approximations,
the use of second-order orbital optimization, different variational active space solvers, and many
other features (Figure <ref>).




 §.§ Seamless MPI functionality

The Message Passing Interface (MPI) is the most popular parallel protocol in the
field of high performance computing.  Although MPI provides high efficiency
for parallel programming, it is a challenge to develop a simple and efficient MPI
program.  In compiled languages, the program must explicitly control 
data communication according to the MPI communication protocol.
The most common design is to activate MPI communication from the beginning and
to update the status of the MPI communicator throughout the program.
When developing new methods, this often leads to extra effort in code
development and debugging.
To sustain the simplicity of the PySCF package, we have designed a
different mechanism to execute parallel code with MPI.
We use MPI to start the Python interpreter as a daemon to receive both the
functions and data on the remote nodes.
When a parallel session is activated, the master process sends to the remote Python
daemons both the functions and the data.  The function is decoded remotely and then
executed.
This design allows one to develop code mainly in serial mode and to switch
to the MPI mode only when high performance is required.  Figure <ref> shows an
example to perform a periodic calculation with and without a parallel session.
Comparing to the serial mode invocation, we see that the user only has to change the
density fitting object to acquire parallel functionality.




§ CONCLUSIONS


Python and its large collection of third party libraries are helping to revolutionize how
we carry out and implement numerical simulations.  It is potentially much more productive
to solve computational problems within the Python ecosystem because it frees researchers
to work at the highest level of abstraction without worrying about the details of complex
software implementation.  To bring all the benefits of the Python ecosystem to quantum
chemistry and  electronic structure simulations, we have
started the open-source PySCF project.

PySCF is a simple, lightweight, and efficient computational chemistry program
package, which supports ab initio calculations for both molecular and extended systems.


The package serves as an extensible electronic structure toolbox, providing a large number
of fundamental operations with simple APIs to manipulate methods, integrals, and wave
functions.  We have invested significant effort to ensure simplicity of use and
implementation while preserving competitive functionality and performance. We believe that
this package represents a new style of program and library design that will be
representative of future software developments in the field.





§ ACKNOWLEDGMENTS


QS would like to thank Junbo Lu and Alexander Sokolov for testing functionality and for
useful suggestions for the program package.  The development of different components of
the  PySCF package has been generously supported by several sources.  Most of the
molecular quantum chemistry software infrastructure was developed with support from the US
National Science Foundation, through grants CHE-1650436 and ACI-1657286. The
periodic mean-field infrastructure was developed with support from ACI-1657286. The
excited-state periodic coupled cluster methods were developed with support from the US
Department of Energy, Office of Science, through the grants DE-SC0010530 and DE-SC0008624.
Additional support for the extended-system methods has been provided by the Simons
Foundation through the Simons Collaboration on the Many Electron Problem, a Simons
Investigatorship in Theoretical Physics, the Princeton Center for Theoretical Science, and
startup funds from Princeton University and the California Institute of Technology.




























































