Metadata-Version: 2.4
Name: PyFraME
Version: 0.5.0
Summary: PyFraME: Python framework for Fragment-based Multiscale Embedding
Author-email: PyFraME Developers <pyframe-devs@pm.me>
Maintainer-email: Jógvan Magnus Haugaard Olsen <pyframe-devs@pm.me>
License-Expression: MPL-2.0
Project-URL: Homepage, https://gitlab.com/pyframe-project/pyframe
Project-URL: Download, https://gitlab.com/pyframe-project/pyframe/-/releases
Project-URL: Documentation, https://gitlab.com/pyframe-project/pyframe/-/tree/master/docs
Project-URL: Changelog, https://gitlab.com/pyframe-project/pyframe/-/blob/master/CHANGELOG.md
Project-URL: Zenodo deposit, https://doi.org/10.5281/zenodo.775113
Project-URL: Issue Tracker, https://gitlab.com/pyframe-project/pyframe/-/work_items
Project-URL: Source, https://gitlab.com/pyframe-project/pyframe
Keywords: computational chemistry,multiscale modeling,QM/MM,polarizable embedding,fragment-based methods
Classifier: Development Status :: 4 - Beta
Classifier: Intended Audience :: Science/Research
Classifier: Natural Language :: English
Classifier: Operating System :: MacOS
Classifier: Operating System :: POSIX :: Linux
Classifier: Programming Language :: C++
Classifier: Programming Language :: Python :: 3 :: Only
Classifier: Programming Language :: Python :: 3.11
Classifier: Programming Language :: Python :: 3.12
Classifier: Programming Language :: Python :: 3.13
Classifier: Programming Language :: Python :: 3.14
Classifier: Topic :: Scientific/Engineering
Classifier: Topic :: Scientific/Engineering :: Chemistry
Classifier: Typing :: Typed
Requires-Python: >=3.11
Description-Content-Type: text/markdown
License-File: LICENSE
Requires-Dist: fabric>=2.0
Requires-Dist: h5py>=3.6
Requires-Dist: MDAnalysis>=2.10
Requires-Dist: mpi4py>=4.1
Requires-Dist: numpy>=1.23.2
Requires-Dist: qcelemental>=0.26
Requires-Dist: scikit-learn>=1.0.2
Requires-Dist: scipy>=1.8
Provides-Extra: dalton
Requires-Dist: daltools; extra == "dalton"
Requires-Dist: loprop; extra == "dalton"
Dynamic: license-file

PyFraME: Python framework for Fragment-based Multiscale Embedding
=================================================================

Copyright (C) 2017-2026  PyFraME Contributors

This Source Code Form is subject to the terms of the Mozilla Public
License, v. 2.0. If a copy of the MPL was not distributed with this
file, You can obtain one at https://mozilla.org/MPL/2.0/.

[![pipeline status](https://gitlab.com/pyframe-project/pyframe/badges/master/pipeline.svg)](https://gitlab.com/pyframe-project/pyframe/-/pipelines)
[![coverage report](https://gitlab.com/pyframe-project/pyframe/badges/master/coverage.svg)](https://gitlab.com/pyframe-project/pyframe/-/pipelines)
[![PyPI version](https://img.shields.io/pypi/v/PyFraME.svg)](https://pypi.org/project/PyFraME/)
[![Python versions](https://img.shields.io/pypi/pyversions/PyFraME.svg)](https://pypi.org/project/PyFraME/)
[![License: MPL 2.0](https://img.shields.io/badge/license-MPL--2.0-blue.svg)](https://mozilla.org/MPL/2.0/)
[![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.775113.svg)](https://doi.org/10.5281/zenodo.775113)


Description
-----------

PyFraME is a Python package providing a framework for managing fragment-based
multiscale embedding calculations. In such calculations, a molecular system is
divided into two principal domains: a central core and its environment. The
core part is treated at the highest level of theory while the effects from the
environment, which can be further subdivided into several differently treated
domains, are included effectively through an embedding potential. Using
PyFraME the user can automatize the workflow starting from an initial structure
to the final embedding potential. It can be used to build a multilayer
description of the molecular environment. Each layer can be described either by
a standard embedding potential, i.e., using a predefined set of parameters, or
by deriving the embedding-potential parameters based on first-principles
calculations. For the latter, a fragmentation method is used to subdivide large
molecular structures into smaller computationally manageable fragments. The
number of layers, as well as the composition and level of theory used for each
layer, can be fully customized.

The basic workflow consists of three main steps. First, a molecular structure
is given as an input. PyFraME reads structures, and the topology and trajectory
files of a molecular dynamics simulation, in any of the formats MDAnalysis
reads, and builds a potential for each of the frames taken from them, a single
structure being one frame. The structure defines the basic units of the system,
i.e., fragments, one per residue. Small molecules typically constitute a
fragment on their own, but larger molecules need to be broken down into smaller
fragments. For example, for proteins, a fragment would usually consist of an
amino-acid residue. The molecular system to be used for the embedding
calculation is then built by selecting fragments with MDAnalysis atom
selections, e.g. by name, chain ID, distance, or a combination thereof, and
placing them into separate regions. As mentioned above, any number of regions
may be added, and each can be fully customized. Once the system has been built,
the final step is the derivation of the embedding potential. Depending on the
specifics, it may involve a large number of separate calculations on the
individual fragments in order to compute the embedding-potential parameters.
For large molecules, where the parameters cannot be computed directly, PyFraME
uses a fragmentation method based on the molecular fractionation with conjugate
caps (MFCC) approach to derive the parameters. The individual fragment
calculations are typically performed by Dalton and the LoProp Python package
but this can be customized. The fragmentation of the system, fragment
calculations, and subsequent joining of parameters to build the embedding
potential are fully automatized and can make full use of large-scale HPC
resources.

For an example showing how PyFraME can be used, see
[Quick start](#quick-start) and [Usage example](#usage-example).


Contents
--------

* [How to cite](#how-to-cite)
* [Contributors](#contributors)
* [Requirements](#requirements)
* [Installation](#installation)
* [Development](#development)
* [Testing](#testing)
* [Documentation](#documentation)
* [Issues](#issues)
* [Contributing](#contributing)
* [Quick start](#quick-start)
* [Trajectories](#trajectories)
* [Usage example](#usage-example)
* [Moving from MolecularSystem](#moving-from-molecularsystem)
* [Parameter sets](#parameter-sets)
* [JSON potentials](#json-potentials)
* [Embedding calculations](#embedding-calculations)
* [License](#license)


How to cite
-----------

To cite PyFraME please use a format similar to the following:

J. M. H. Olsen, P. Reinholdt, J. Vester, T. N. Heegaard, and G. Berger, *PyFraME: Python framework for Fragment-based
Multiscale Embedding (version 0.5.0)*, **2026**. DOI: 10.5281/zenodo.22996201.
See https://gitlab.com/pyframe-project/pyframe.

where the version and DOI should correspond to the actual version that was
used. Note that the DOI
[10.5281/zenodo.775113](https://doi.org/10.5281/zenodo.775113)
represents all versions, and will always resolve to the latest one. A possible
BibTeX entry can be found in the
[CITATION](https://gitlab.com/pyframe-project/pyframe/-/blob/master/CITATION)
file, and the same entry is printed by `pyframe.print_citation()`.
Alternatively, BibTeX and other formats can be generated by
[Zenodo](https://doi.org/10.5281/zenodo.775113).


Contributors
------------

The list of past and current contributors is found
[here](https://gitlab.com/pyframe-project/pyframe/-/graphs/master).


Requirements
------------

To use PyFraME you need:

* [Python (3.11+)](https://www.python.org/)
* [NumPy](https://www.numpy.org/)
* [SciPy](https://www.scipy.org/)
* [scikit-learn](https://scikit-learn.org/)
* [h5py](https://www.h5py.org/)
* [mpi4py](https://mpi4py.readthedocs.io/)
* [Fabric](https://www.fabfile.org/)
* [QCElemental](https://molssi.github.io/QCElemental/)
* [MDAnalysis (2.9+)](https://www.mdanalysis.org/)
* [daltools](https://github.com/vahtras/daltools) and
  [LoProp for Dalton](https://github.com/vahtras/loprop)

These are all available on PyPI and are installed along with the package,
except daltools and LoProp, which `pip install pyframe[dalton]` adds.
MDAnalysis publishes no wheels for Linux on aarch64, so installing PyFraME
there builds MDAnalysis from source, which needs a C compiler. A parallel
calculation needs an MPI library as well, which mpi4py loads when it is first
given a communicator; a serial one needs none. Building PyFraME from source needs a C++ compiler with OpenMP support and the
[Eigen](https://eigen.tuxfamily.org/) headers, which the compiled embedding
engine is built against.

The headers usually need not be pointed at: the build asks `pkg-config` for
them, and then looks where the distribution, Homebrew, and MacPorts packages
(`libeigen3-dev`, `eigen3-devel`, `eigen`) put them, so installing Eigen through
the package manager is all that is needed. Where it sits somewhere else, or where
a particular copy of it is wanted, point `EIGEN_INCLUDE_DIR` at the directory
that contains `Eigen/Dense`, which overrides the search:

```bash
EIGEN_INCLUDE_DIR=/path/to/eigen3 python -m pip install PyFraME
```

The fragment calculations themselves are run by a quantum chemistry program,
which is not installed by the above:

* [Dalton](https://www.daltonprogram.org/), which is the default and is used
  together with LoProp
* [Molcas](https://www.molcas.org/) (not tested recently)
* [OpenMolcas](https://gitlab.com/Molcas/OpenMolcas) (same)

To run the test suite you need the `test` dependency group of
`pyproject.toml`, which is part of the `dev` group described under
[Development](#development):

* [pytest](https://pytest.org)
* [pytest-timeout](https://github.com/pytest-dev/pytest-timeout)
* [parse](https://github.com/r1chardj0n3s/parse)
* [PyYAML](https://pyyaml.org)


Installation
------------

The PyFraME package can be installed from
[PyPI](https://pypi.org/project/PyFraME/) directly using
[pip](https://pip.pypa.io/en/stable/)

```bash
python -m pip install PyFraME
```

This will also install required dependencies (see above) that are available on
PyPI, i.e., not Dalton, Molcas, etc.

Binary wheels are provided for Linux on x86_64 and aarch64 and for macOS 15 or
newer on Apple silicon, for CPython 3.11 to 3.14. Anywhere else pip builds the
package from source, which needs the compiler, OpenMP, and Eigen described under
[Requirements](#requirements).

The entire source including history can be found at
[GitLab](https://gitlab.com/pyframe-project/pyframe).
All releases are also deposited at
[Zenodo](https://doi.org/10.5281/zenodo.775113).


Development
-----------

To work on PyFraME, clone the repository and install it in editable mode with
the `dev` dependency group, which brings pytest, mypy, and pre-commit along with
the rest. Dependency groups need pip 25.1 or newer:

```bash
git clone https://gitlab.com/pyframe-project/pyframe.git
cd pyframe
python -m pip install --editable . --group dev
pre-commit install
```

The editable install compiles the embedding engine, so the compiler, OpenMP and
Eigen described under [Requirements](#requirements) are needed here as well.

The hooks that `pre-commit install` sets up run
[ruff](https://docs.astral.sh/ruff/) for linting, import sorting and
formatting, [mypy](https://mypy-lang.org/) for type checking,
clang-format over the C++ sources of the engine, and codespell over
everything. They are what the pipeline gates on, so running them locally is the
quickest way to find out whether a change will pass. To run them over the whole
tree rather than over what is staged:

```bash
pre-commit run --all-files
```


Testing
-------

The tests are not part of the wheel on PyPI. They come with the source, from
[GitLab](https://gitlab.com/pyframe-project/pyframe) or the sdist, and are run
against an editable install as described under [Development](#development). To
run the full test suite, i.e., the unit tests under `pyframe` and the
integration tests under `tests`, type

```bash
python -m pytest
```

in the PyFraME root directory, with no path: `python -m pytest pyframe` leaves
the integration tests out.

The tests of the task farming are marked `cluster`. They start workers over SSH
and therefore need a host that can connect to itself. Where it cannot, deselect
them, as the pipeline does:

```bash
python -m pytest -m 'not cluster'
```

The tests that need an MPI library skip themselves where mpi4py cannot load
one, and those that launch several ranks skip themselves where there is no
`mpirun` or `mpiexec` on the `PATH`, so the suite also runs without an MPI. The
pipeline uses MPICH. Open MPI works as well, provided the machine has at least
five cores, since it refuses to start more ranks than there are cores and the
parallel tests start five. On a machine that keeps its MPI in an environment
module, load the module in the same shell before running the tests, or they
skip.


Documentation
-------------

The API is documented in the docstrings of the package, which
`python -m pydoc pyframe` and `help()` read, and the sections below cover the
workflow from a structure or a trajectory to a potential.

Three notebooks go further than the sections below.
[docs/embedding.ipynb](https://gitlab.com/pyframe-project/pyframe/-/blob/master/docs/embedding.ipynb)
runs the whole workflow on a small system from the test suite, from the
structure to the potential and on into an embedding calculation, and needs no
quantum chemistry program to do it.
[docs/electrostatic_embedding.ipynb](https://gitlab.com/pyframe-project/pyframe/-/blob/master/docs/electrostatic_embedding.ipynb)
walks through the embedding library described under
[Embedding calculations](#embedding-calculations), starting from a UML diagram
of it.
[docs/host_program.ipynb](https://gitlab.com/pyframe-project/pyframe/-/blob/master/docs/host_program.ipynb)
embeds a host program in full, with PySCF as the host: an integral driver of
three methods, an embedding model, and the few lines that attach the model to an
SCF. It needs PySCF to run.

What has changed between versions, including the changes that are not backwards
compatible, is recorded in the
[changelog](https://gitlab.com/pyframe-project/pyframe/-/blob/master/CHANGELOG.md).

For the background of the method and for examples of how PyFraME is used in
computational studies, see our
[tutorial review](https://doi.org/10.1002/qua.25717) article.


Issues
------

Please report issues [here](https://gitlab.com/pyframe-project/pyframe/-/work_items).


Contributing
------------

Please take a look at the
[contribution guide](https://gitlab.com/pyframe-project/pyframe/-/blob/master/CONTRIBUTING.md).


Quick start
-----------

The shortest way from a structure to an embedding potential is to treat one
molecule quantum mechanically and to take the parameters of everything around
it from a parameter set, so that no fragment calculation has to be run at
all:

```python
import pyframe

# Read the structure. A single structure is a trajectory of one frame, read
# through MDAnalysis, and a residue of the file is a fragment.
trajectory = pyframe.Trajectory('solvated_solute.pdb')

# Put the solute in the core region, i.e. the part treated quantum
# mechanically, and give it the basis set to be used for it. A region is given
# as an MDAnalysis atom selection.
trajectory.set_core_region('resname LIG', basis='pcseg-2')

# Describe the water within 8 angstrom of it by the TIP3P parameter set.
trajectory.add_region('solvent', 'resname WAT',
                      distance=8.0, reference='resname LIG',
                      use_parameter_set=True,
                      parameter_set='tip3p')

# Assemble the potential and write the two files it takes to run the embedding
# calculation: the core region as a molecule file and the environment as a
# potential file. The structure is the one snapshot of the trajectory.
snapshot = trajectory.snapshots[0]
project = pyframe.Project(work_dir='/path/to/work')
project.create_embedding_potential(snapshot)
project.write_core(snapshot)
project.write_potential(snapshot)
```

Regions whose parameters are computed rather than tabulated are set up the same
way, with the level of theory of the fragment calculations given to
`add_region` instead; the [usage example](#usage-example) shows both together.

The same, for every frame of a molecular dynamics trajectory rather than for a
single structure:

```python
import pyframe

# Read the frames through MDAnalysis, with the solute unwrapped and centered in
# the box in every frame.
trajectory = pyframe.Trajectory(topology_file='topology.tpr',
                                trajectory_file='trajectory.xtc',
                                center_selection='resname LIG')

# The selections are resolved again in every frame, so the solvent region
# follows what is within 8 angstrom of the solute as the frames go by.
trajectory.set_core_region('resname LIG', basis='pcseg-2')
trajectory.add_region('solvent', 'resname SOL',
                      distance=8.0, reference='resname LIG',
                      use_parameter_set=True,
                      parameter_set='tip3p')

# One potential per frame, in files named after the frame they came from.
trajectory.create_potential(pyframe.Project(work_dir='/path/to/work'))
trajectory.write_core()
trajectory.write_potential()
```

The next section covers this in full.


Trajectories
------------

A `Trajectory` reads a topology and a trajectory file through
[MDAnalysis](https://www.mdanalysis.org/), so any pair of formats it supports
can be used, or a single file carrying both, e.g. a PDB, which is a trajectory
of one frame per model. A single structure is a trajectory of one frame. Every
selected frame becomes a `Snapshot`, carrying the fragments, the regions, and
the simulation box of that frame, and a residue of the topology is a fragment.

A topology that carries no elements or partial charges, e.g. a plain PDB, is
read as well. The element is guessed from the atom name, which is warned about
once, and the charge is taken from the charge column of a PDB, so that an ion
keeps its charge. A PQR is read with its fields separated by whitespace, or in
the columns of the format where they run together. The coordinates and box are
read at the precision of the file for the common formats, and a warning says so
where they are not.

```python
import pyframe

# Take every tenth of the first thousand frames. The core region is unwrapped
# and centered in the box, and everything else is wrapped around it, so that a
# molecule split across the periodic boundary is whole in the snapshots.
trajectory = pyframe.Trajectory(topology_file='/path/to/topology.tpr',
                                trajectory_file='/path/to/trajectory.xtc',
                                start=0, stop=1000, step=10,
                                center_selection='resname LIG')

# Regions are given as MDAnalysis atom selections. A fragment belongs to the
# first region that selects it, so a later selection takes only what is left,
# and 'all' takes the rest.
trajectory.set_core_region('resname LIG', basis='pcseg-2')
trajectory.add_region('solvent', 'resname SOL',
                      distance=8.0, reference='resname LIG',
                      use_parameter_set=True,
                      parameter_set='tip3p')

# check_selection reports which fragments a selection covers, which is worth
# checking before anything is run.
print(trajectory.check_selection('resname LIG'))

# The settings of the fragment calculations are those of a Project, which every
# snapshot is handed to in turn. The files are named after the frame they came
# from, e.g. "trajectory_frame-10.pot".
project = pyframe.Project(work_dir='/path/to/work')
trajectory.create_potential(project)
trajectory.write_core()
trajectory.write_potential()
```

A selection is resolved again in every frame by default, so that a region
defined by a distance follows what is within that distance as the frames go
by. With `updating=False` it is resolved once, in the first frame taken, and
the region holds the same fragments throughout. A region that would hold no
fragment in any frame is refused, and the frames in which one holds none are
reported in the log.

The `distance` of a region is measured to the fragments of its `reference`
selection, which may be in a region already, between centers of mass or, with
`use_center_of_mass=False`, between nearest atoms. With `protect_molecules`,
which is the default, every fragment bonded to one within the distance is taken
as well, so that no molecule is cut. The distance is measured as the coordinates
lie, which suits frames centered with `center_selection`, and to the nearest
periodic image with `periodic=True`. A distance written into the selection
itself, e.g. `'resname SOL and around 8.0 resname LIG'`, is MDAnalysis's: always
between atoms, and to the nearest image where the frame has a box.

A fragment can be split into named parts, e.g. a lipid into its head group and
tails, before any region is set. An atom goes to the first part with a pattern
that its name matches, the patterns being regular expressions, and the parts
are residues of their own afterwards, which a selection can name:

```python
trajectory.split_fragments('resname POPC',
                           {'HEAD': ['N', 'C1[1-5]', 'H1[1-5][A-C]', 'P', 'O1[1-4]'],
                            'TAIL': ['.*']})
trajectory.add_region('lipid', 'resname HEAD TAIL', use_mfcc=True,
                      use_multipoles=True, use_polarizabilities=True)
```

`reset()` drops every region, potential, and split and starts the setup over.
A single structure is read the same way, and its snapshot and files are named
after it alone, e.g. `structure.pot`:

```python
trajectory = pyframe.Trajectory('/path/to/structure.pdb')
snapshot = trajectory.snapshots[0]
```

A universe that is already there, e.g. one read with options of MDAnalysis's
own or built in memory, is taken with `Trajectory.from_universe(universe,
name='run')`, which works on a copy of it. The frame each snapshot was built
from is in `Trajectory.frames`, and the MDAnalysis universe is in
`Trajectory.universe`. The frames are selected once, when the trajectory is
built, so a run over a different set of frames is a new `Trajectory`.


Usage example
-------------

The following commented example is based on a molecular system consisting of a
channelrhodopsin protein dimer embedded in a lipid membrane.

```python
import pyframe

# Read the structure. Any format MDAnalysis reads will do, and a residue of the
# file is a fragment.
trajectory = pyframe.Trajectory('/path/to/input/file.pdb')

# By default fragments are defined by the input but fragments can be split into
# parts as shown here, which affects every residue the selection touches.
trajectory.split_fragments(
        'resname RETK',
        {'LYSB': ['N', 'H', 'CA', 'HA', 'C', 'O'],
         'LYSS': ['CB', 'HB1', 'HB2', 'CG', 'HG1', 'HG2', 'CD',
                  'HD1', 'HD2', 'CE', 'HE1', 'HE2'],
         'RET': ['.*']})
trajectory.split_fragments(
        'resname POPE',
        {'POP1': ['N', 'HN1', 'HN2', 'HN3', 'C12', 'H12A', 'H12B', 'C11', 'H11A',
                  'H11B', 'P', 'O13', 'O14', 'O11', 'O12', 'C1', 'HA', 'HB', 'C2',
                  'HS', 'O21', 'C3', 'HX', 'HY', 'O31'],
         'POP2': ['C21', 'O22', 'C22', 'H2R', 'H2S', 'C23', 'H3R', 'H3S', 'C24',
                  'H4R', 'H4S', 'C25', 'H5R', 'H5S', 'C26', 'H6R', 'H6S', 'C27',
                  'H7R', 'H7S', 'C28', 'H8R', 'H8S', 'C29', 'H91'],
         'POP3': ['0C21', '1H10', '1C21', 'H11R', 'H11S', '2C21', 'H12R', 'H12S',
                  '3C21', 'H13R', 'H13S', '4C21', 'H14R', 'H14S', '5C21', 'H15R',
                  'H15S', '6C21', 'H16R', 'H16S', '7C21', 'H17R', 'H17S', '8C21',
                  'H18R', 'H18S', 'H18T'],
         'POP4': ['C31', 'O32', 'C32', 'H2X', 'H2Y', 'C33', 'H3X', 'H3Y', 'C34',
                  'H4X', 'H4Y', 'C35', 'H5X', 'H5Y', 'C36', 'H6X', 'H6Y', 'C37',
                  'H7X', 'H7Y', 'C38', 'H8X', 'H8Y', 'C39', 'H9X', 'H9Y'],
         'POP5': ['0C31', 'H10X', 'H10Y', '1C31', 'H11X', 'H11Y', '2C31', 'H12X',
                  'H12Y', '3C31', 'H13X', 'H13Y', '4C31', 'H14X', 'H14Y', '5C31',
                  'H15X', 'H15Y', '6C31', 'H16X', 'H16Y', 'H16Z']})

# The retinal and whatever is within 3 angstrom of it, by the nearest atoms
# and without taking whole molecules along, make up the core region.
trajectory.set_core_region('all', distance=3.0,
                           reference='resid 248 and chainID A and resname RET',
                           use_center_of_mass=False, protect_molecules=False,
                           basis='pcseg-2')

# Add a region containing the protein (here chain A holds all of it). Note that
# each of these settings have defaults and that there are more than those shown
# here.
trajectory.add_region('protein', 'chainID A', use_mfcc=True,
                      mfcc_order=2, use_multipoles=True, multipole_order=2,
                      use_polarizabilities=True, basis='loprop-6-31+G*')

# Here we repeat for lipids, ions, and solvent within 8 angstrom of the protein,
# which may be in a region already to be measured to.
trajectory.add_region('lipid', 'resname POP1 POP2 POP3 POP4 POP5',
                      distance=8.0, reference='chainID A',
                      use_mfcc=True, mfcc_order=2,
                      use_multipoles=True, multipole_order=2,
                      use_polarizabilities=True, basis='loprop-6-31+G*')
trajectory.add_region('ion', 'resname NA CL', distance=8.0, reference='chainID A',
                      use_multipoles=True, multipole_order=0,
                      use_polarizabilities=True, basis='6-31+G*')
trajectory.add_region('solvent', 'resname SOL', distance=8.0,
                      reference='chainID A', use_multipoles=True,
                      multipole_order=2, use_polarizabilities=True,
                      basis='loprop-6-31+G*')

# Create Project() object that is used to create embedding potentials and write
# input files. Every setting shown below has a default, and there are more of
# them than are shown here.
project = pyframe.Project()

# Set path to scratch directory. This will be used by the auxiliary programs,
# e.g. Dalton or Molcas.
project.scratch_dir = '/path/to/scratch'

# Set path to working directory (it will be created if it does not exist).
# This directory will contain the final output files from PyFraME (e.g. Dalton
# mol and pot files), and the output from the auxiliary program. In addition,
# during execution it will contain temporary directories for each fragment.
project.work_dir = '/path/to/work'

# Specifies the number of jobs that will be run on each node. A fragment may
# require one or more calculations run by an auxiliary program. Each of these
# counts as a job.
project.jobs_per_node = 2

# Specifies memory per job. Note that this amount will be shared by MPI processes.
project.memory_per_job = 2048 * 12

# Number of MPI processes per job.
project.mpi_procs_per_job = 12

# You can manually specify the name of nodes that should be used to run jobs.
# PyFraME will attempt to autodetect nodes from SLURM and PBS/Torque queuing
# systems. For example, with os imported:
# project.node_list = [os.environ['HOSTNAME']]

# This will start the fragment calculations using the auxiliary
# programs and settings defined when creating the regions. The structure is the
# one snapshot of the trajectory.
snapshot = trajectory.snapshots[0]
project.create_embedding_potential(snapshot)

# Write potential file containing all parameters of the embedding potential.
project.write_potential(snapshot)

# Write molecule file containing the core quantum region.
project.write_core(snapshot)
```

The calculations of a snapshot are run in a directory of its own in the work
directory, and are named after the index of the residue of their fragment, e.g.
`fragment-9_dalton_multipoles.out`, or of the two fragments a concap caps, e.g.
`concap-8-9_dalton_multipoles.out`. `fragments.txt` beside them says which
fragment each is of, as messages about a fragment do, e.g. `ILE 10, chain A,
resindex 9, atoms index 120:138`, where `resindex 9` selects it in the
universe of the trajectory. A calculation whose output is there already is not
run again, so one that failed is run again by deleting its output and creating
the potential again.

The fragment calculations are run by workers that PyFraME starts on each node,
over SSH on any node other than the one it runs on, and each worker and each job
runs in a login shell, i.e. `bash -l`. Each job is given the `PATH` and
`LD_LIBRARY_PATH` of the process that started the run. Where that is not enough
for the programs to be found, e.g. where Dalton is kept in an environment
module, two environment variables give commands to run first:
`PyFraME_WORKER_COMMANDS` before each worker is started, and
`PyFraME_JOB_COMMANDS` before each job:

```bash
export PyFraME_JOB_COMMANDS='module load dalton'
```

A job runs only if its commands succeed. A third variable,
`PyFraME_MPI_COMMAND`, sets the command that a fragment calculation is launched
with where the program is run through MPI, `mpirun -map-by core -np` by
default. It is read when the `Project` is created, an `mpi_command` given to the
`Project` takes precedence over it, and it has to end with the option that takes
the number of processes, which PyFraME adds:

```bash
export PyFraME_MPI_COMMAND='srun -n'
```

The uppercase spellings, `PYFRAME_WORKER_COMMANDS`, `PYFRAME_JOB_COMMANDS`, and
`PYFRAME_MPI_COMMAND`, are accepted as well. All three are for setting up the
environment, and never for a password, a token, or any other secret. The worker
and job commands become part of a command line that other users on a node can
read in a process listing, and of the message logged when a worker fails to
start, and the MPI command is written into the run script of every calculation
that uses it and logged.


Moving from MolecularSystem
---------------------------

`MolecularSystem` was removed in 0.5.0. A structure is read as a `Trajectory`
of one frame, and the regions are given as MDAnalysis atom selections rather
than as fragments taken from the system one call at a time. What a script of
0.4 did is done as follows:

| 0.4 | 0.5 |
|---|---|
| `system = pyframe.MolecularSystem(input_file='file.pdb')` | `trajectory = pyframe.Trajectory('file.pdb')`, whose one snapshot is `trajectory.snapshots[0]` |
| `MolecularSystem(..., input_reader=reader)` | read the file into an MDAnalysis universe, and `pyframe.Trajectory.from_universe(universe, name='system')` |
| `get_fragments_by_name(['WAT'])` | `'resname WAT'` |
| `get_fragments_by_number([1, *range(3, 11)])` | `'resid 1 3:10'` |
| `get_fragments_by_chain_id(['A'])` | `'chainID A'` |
| `get_fragments_by_identifier(['248_A_RET'])` | `'resid 248 and chainID A and resname RET'` |
| `get_fragments_by_charge([1, -1])` | no counterpart; select the charged fragments by name, e.g. `'resname NA CL'` |
| `get_fragments_by_distance(8.0, reference=core)` | `distance=8.0, reference='resname LIG'` given to `add_region` or `set_core_region`, with the selection it narrows down, e.g. `'all'` |
| `get_fragments_by_distance_and_name(8.0, ['WAT'], reference=core)` | `add_region('solvent', 'resname WAT', distance=8.0, reference='resname LIG')` |
| `get_remaining_fragments()` | `'all'`, which takes what no earlier region took |
| `system.add_region(name, fragments, ...)` | `trajectory.add_region(name, selection, ...)` |
| `system.set_core_region(fragments, ...)` | `trajectory.set_core_region(selection, ...)` |
| `split_fragment_by_name('POPC', names, definitions)` | `trajectory.split_fragments('resname POPC', dict(zip(names, definitions)))` |
| `split_fragment_by_identifier('1_A_VAL', names, definitions)` | `trajectory.split_fragments('resid 1 and chainID A', dict(zip(names, definitions)))` |
| `system.reset()` | `trajectory.reset()` |
| `project.create_embedding_potential(system)` | `project.create_embedding_potential(snapshot)`, or `trajectory.create_potential(project)` for every snapshot |
| `project.write_core(system)`, `project.write_potential(system)`, `system.to_subsystems()` | the same, given the snapshot |

The options of a selection by distance keep their names and defaults:
`use_center_of_mass=True` and `protect_molecules=True`. A selection is an atom
selection, so a fragment is taken if any of its atoms is selected.

Standard potentials are called parameter sets, and the settings of a region
that take its parameters from one are renamed: `use_standard_potentials`,
`standard_potential_model`, and `standard_potential_exclusion_type` are
`use_parameter_set`, `parameter_set`, and `parameter_set_exclusion_type`. An
old name is refused with an error that gives the new one. A region that names
no parameter set takes `pyframe-2026` rather than `sep`.

A fragment is keyed by the index of its residue rather than by a name such as
`10_A_ILE`, and the files of its calculations are named after that index, e.g.
`fragment-9_...`, so the calculations of a work directory of 0.4 are run
again. Where a structure gives no elements, they are guessed by MDAnalysis,
which reads some atom names better, e.g. those of a protonated histidine,
`HIP`.


Parameter sets
--------------

A region can take its parameters from a parameter set, i.e., a table of
parameters per fragment type and atom name, rather than having them computed:

```python
trajectory.add_region('solvent', 'resname SOL',
                      use_parameter_set=True,
                      parameter_set='tip3p')
```

The parameter sets shipped with PyFraME are `alep`, `cp3`, `ff03`, `ff94`,
`lipid14`, `lipid17`, `pyframe-2026`, `sep`, and `tip3p`. A region that names
none uses `pyframe-2026`, which merges three of them, each the parameters of one
publication: `cp3` for the amino acids, with their terminal variants, `sep` for
water, common organic solvents, and ions, and `alep` for the lipids DMPC, POPC,
and POPS. A protein in a membrane, with its solvent and ions, can therefore be
a single region. The phenol of `sep` is named `PHN`, so that `PHE` is
phenylalanine. A parameter set merged from several tables takes each fragment
whole from the first table that has it.

`alep` gives the parameters of the head group and the tails of a lipid, which
carry fractional charges that only add up to that of the whole lipid. A lipid
is selected and given as a whole, e.g. `'resname POPC'`, and the parameter set
splits it into its parts itself. The sites of a lipid are excluded as the
parameters were made to be used, i.e. each part with caps of three atoms into
the parts it is bonded to, whatever the region says, so the far ends of a lipid
polarize each other. A lipid split by hand with `split_fragments` is treated the
same way, and its parts are corrected to the charge of the lipid as a whole.

A program that offers a choice of them, e.g. a host program, can list the
parameter sets and ask what each provides, i.e. its description, its tables,
the fragments it covers, and its parameters. The description names the
publication the parameters are from, with its DOI, for `alep`, `cp3`, `sep`,
and `pyframe-2026`:

```python
for name in pyframe.available_parameter_sets():
    info = pyframe.describe_parameter_set(name)
    print(name, info.multipole_order, info.polarizabilities, info.lennard_jones, info.description)
```

A file named after a parameter set with a `.csv` extension in the working
directory takes precedence over the one shipped with the package, so a
parameter set of your own is used by placing it there and naming it in
`parameter_set`. It is not listed by `available_parameter_sets`, but it is
described. A `pyframe-2026.csv` there replaces all of `pyframe-2026`, and a
`sep.csv` replaces `sep`, both by itself and as the part of `pyframe-2026` that
it is.

A parameter set is a CSV file whose first two columns hold the fragment
name and the atom name, followed by one column per parameter component. Every
column states the unit of its values in brackets after its label. Lines above
the header that begin with `#` are ignored, except that one beginning with
`# description:` gives the description that `describe_parameter_set`
reports. For example, a parameter set of charges and isotropic dipole-dipole
polarizabilities:

```
name,atom,M_000[e],P1_iso[bohr^3]
WAT,O,-0.67444,5.73935
WAT,H,0.33722,2.30839
```

Two more lines above the header say how the fragments of a parameter set are
used. `# exclusion: mfcc 3` gives the exclusion type and cap order of the sites
of its fragments, which then apply in place of the `parameter_set_exclusion_type`
and `mfcc_order` of the region. `# parts: POPC = POCH POCO POCP` makes a
fragment of the parts listed, each atom taking the parameters of the part that
lists its name, so that a structure can name the whole fragment. These are the
lines `alep` begins with:

```
# exclusion: mfcc 3
# parts: DMPC = DMCH DMCM
# parts: POPC = POCH POCO POCP
# parts: POPS = POSH POSO POSP
```

A component of a multipole is labeled `M_` followed by its multi-index, i.e.,
by how many times x, y, and z occur among its Cartesian indices. `M_000` is the
charge, `M_100`, `M_010`, and `M_001` are the x, y, and z components of the
dipole, and `M_200`, `M_110`, `M_101`, `M_020`, `M_011`, and `M_002` are the six
components of the quadrupole. Multipoles up to order six are supported.

A component of a polarizability carries one multi-index per operator, since the
two operators need not be the same: `P_100_010` is the xy component of the
dipole-dipole polarizability and `P_000_100` the x component of the
charge-dipole one. An isotropic polarizability, which is defined only where the
two orders are equal, is labeled by its order instead, so that `P1_iso` is an
isotropic dipole-dipole polarizability and `P2_iso` an isotropic
quadrupole-quadrupole one.

The Lennard-Jones parameters of a site are labeled `sigma` and `epsilon`. For
example, the TIP3P oxygen:

```
name,atom,M_000[e],epsilon[kJ/mol],sigma[nm]
HOH,O,-0.834,0.636386,0.315061
```

Every value is converted to the atomic units a site holds it in. The units each
column may be given in are:

| Column | Units |
|---|---|
| charge, `M_000` | `e` |
| dipole, `M_100` etc. | `e*bohr`, `e*angstrom`, `debye` |
| quadrupole, `M_200` etc. | `e*bohr^2`, `e*angstrom^2`, `debye*angstrom` |
| multipole of order k | `e*bohr^k`, `e*angstrom^k` |
| polarizability of orders k and l, `P_100_010` or `P1_iso` etc. | `bohr^n`, `angstrom^n`, with n = k + l + 1 |
| `sigma` | `bohr`, `angstrom`, `nm` |
| `epsilon` | `hartree`, `kJ/mol`, `kcal/mol`, `eV`, `K` |

and any column may be given in `au`. A polarizability is given as a
polarizability volume, i.e., the atomic unit of polarizability divided by
4 pi epsilon_0, so that a dipole-dipole polarizability is in `bohr^3` or
`angstrom^3` and the value in `bohr^3` is the value in atomic units. An epsilon
in `K` is the well depth divided by the Boltzmann constant.

A column without a unit, or in a unit its parameter cannot be given in, is
refused, since the same number is a different potential in each unit.

The header is checked before any parameter is read. A multipole or a
polarizability must have a column for each of its components, none of them
twice, and a P(k,l) polarizability with k different from l must have the
columns of its P(l,k) polarizability as well, the two being transposes of each
other only where the two operators are the same. A polarizability given both as
an isotropic value and component by component is refused, as is a label that is
neither of the above, e.g. the `M0` and `P11` of the format that took all
components of a parameter in a single column, or a column holding anything other
than a multipole, a polarizability, `sigma`, or `epsilon`.


JSON potentials
---------------

A finished potential is written either as a PElib potential file or as a JSON
file, i.e. `snapshot.write_potential(filetype='json')`, and the JSON file is the
one that `pyframe.embedding` reads back through `read_input.reader`. It carries
what a potential file has no place for: the simulation box, the core region and
the Lennard-Jones parameters of every site.

A complete file, for a hydrogen fluoride molecule described quantum
mechanically with a sodium and a chloride ion around it:

```json
{
  "generated_by": "Generated by PyFraME (version 0.5.0)",
  "units": {
    "length": "bohr",
    "energy": "hartree",
    "angle": "degree",
    "multipoles": "au",
    "polarizabilities": "au"
  },
  "quantum_subsystems": [
    {
      "name": "core region",
      "charge": 0,
      "basis": [
        "6-31+G*",
        "6-31+G*"
      ],
      "quantum_fragments": [
        {
          "index": 0,
          "name": "HF 20, resindex 0, atoms index 0:1",
          "nuclei": [
            {
              "index": 0,
              "element": "F",
              "charge": 9,
              "coordinate": [
                0.0,
                0.0,
                0.0
              ]
            },
            {
              "index": 1,
              "element": "H",
              "charge": 1,
              "coordinate": [
                1.7385480346557087,
                0.0,
                0.0
              ]
            }
          ]
        }
      ]
    }
  ],
  "classical_subsystems": [
    {
      "classical_fragments": [
        {
          "name": "NA 10, resindex 1, atom index 2",
          "index": 0,
          "atoms": [
            {
              "index": 0,
              "element": "Na",
              "coordinate": [
                5.669178373877311,
                0.0,
                0.0
              ],
              "exclusions": [
                0
              ],
              "multipoles": {
                "elements": [
                  1.0
                ],
                "order": 0
              }
            }
          ]
        },
        {
          "name": "CL 11, resindex 2, atom index 3",
          "index": 1,
          "atoms": [
            {
              "repulsion": {
                "method": "LJ",
                "parameters": {
                  "lj_sigma": 4.4,
                  "lj_epsilon": 0.1
                }
              },
              "dispersion": {
                "method": "LJ",
                "parameters": {
                  "lj_sigma": 4.4,
                  "lj_epsilon": 0.1
                }
              },
              "index": 1,
              "element": "Cl",
              "coordinate": [
                -5.669178373877311,
                0.0,
                0.0
              ],
              "exclusions": [
                1
              ],
              "multipoles": {
                "elements": [
                  -1.0
                ],
                "order": 0
              }
            }
          ]
        }
      ]
    }
  ]
}
```

The file states the units of what it holds in its `units` block, and one
without the block is refused when it is read. The block gives one unit for each
kind of quantity:

| Key | Gives the unit of | Units |
|---|---|---|
| `length` | coordinates, box lengths, Lennard-Jones sigma | `bohr`, `angstrom`, `nm` |
| `energy` | Lennard-Jones epsilon | `hartree`, `kJ/mol`, `kcal/mol`, `eV`, `K` |
| `angle` | box angles | `degree`, `radian` |
| `multipoles` | multipoles, a unit for every order | `au`, `e*bohr^k`, `e*angstrom^k` |
| `polarizabilities` | polarizabilities, a unit for every pair of orders | `au`, `bohr^(k+l+1)`, `angstrom^(k+l+1)` |

and every value is converted to bohr, hartree, degrees, and atomic units as it is
read. A file PyFraME writes is in those already: a system carries its
coordinates in angstrom everywhere else in PyFraME, and the writer converts them
on the way out, while the parameters of a site are held in atomic units. A
multipole or polarizability in other than atomic units is converted order by
order, so its components have to be every component up to the order the
multipoles are of, and the `order` of polarizabilities has to give the two
orders of every block, as the writer gives them.

Everything is numbered from zero: the nuclei, the fragments of either
subsystem, and the atoms of a classical subsystem. An index is an identity
rather than a position, i.e. it is matched by value where an exclusion list
names it, so a file numbered from one is read as well, provided it is numbered
that way throughout. The potential file numbers its sites from one instead,
because its exclusion lists are padded to a common width and a zero is what
marks the padding; a JSON file has no such need, its exclusion lists being
arrays of their own length.

`generated_by` records the version that wrote the file. It is ignored when the
file is read, as is any other key the reader does not know at the top level.

A system read from a structure that carries a simulation box states it under
`system`, in the units of the block, and one that carries none leaves the key out
altogether:

```json
"system": [{"simulation_box": [{"lengths": [56.69178373877310,
                                            56.69178373877310,
                                            56.69178373877310],
                                "angles": [90.0, 90.0, 90.0]}]}]
```

The quantum subsystem is the core region: its formal charge under `charge`, its
basis set under `basis`, and its fragments under `quantum_fragments`. The basis
is one label per nucleus, in the order the nuclei come in, whichever of the
three forms the region was given it in, i.e. one label for the whole region,
one per element or one per atom, so that a program reading the file needs
neither the element table nor the atom order to resolve a label; a list that
does not cover the nuclei one for one is refused rather than read. The nuclei
live in the fragments and nowhere else, and are numbered across the whole
subsystem rather than restarting in each fragment, since a nucleus index and a
basis label are positional in the subsystem. Each nucleus states its `kind`:
`atom` for an atom of the core region, `link` for a hydrogen link atom and
`cap` for a capping atom, where a file that gives none is read as `atom`. A
file may give them in a flat
`nuclei` list beside the fragments instead, which is how they were written
before and is still read; a file that gives both is read from the flat list.

The classical subsystem is the environment, as one classical fragment per
fragment of the system, each named by the description of the fragment it is,
which is how a person finds it, and under `region` the name of the region it
belongs to. Its
atoms are the sites of the potential, numbered across the whole subsystem
rather than within a fragment, i.e. the second fragment of the file above
carries atom 1 rather than a second atom 0, since that is the numbering
`exclusions` refers to. The first entry of an exclusion list is the atom
itself. A water and a sodium ion, as two fragments of one subsystem, therefore
exclude as follows, with the rest of each atom left out:

```json
{"name": "WAT 1, resindex 0, atoms index 0:2", "index": 0,
 "atoms": [{"index": 0, "exclusions": [0, 1, 2]},
           {"index": 1, "exclusions": [1, 0, 2]},
           {"index": 2, "exclusions": [2, 0, 1]}]}
{"name": "NA 2, resindex 1, atom index 3", "index": 1,
 "atoms": [{"index": 3, "exclusions": [3]}]}
```

An exclusion that names no atom of the subsystem is refused when the file is
read, since it excludes nothing and is what an input whose atoms and exclusion
lists are numbered differently comes to. A potential that was not built from
the regions of the system, e.g. one read back from a potential file, does not
line up with them and is written as a single fragment named `all sites`.

The parameters of a site are its multipoles, its polarizabilities, and its
Lennard-Jones parameters. Multipoles and polarizabilities are written as one
flat array of components per site, padded with zeros to the highest orders
found anywhere in the system so that every site has an array of the same shape,
with `order` naming what the array holds: the highest multipole order, and the
two orders of every polarizability block in the order the blocks come in. A
site carrying a charge and an isotropic dipole-dipole polarizability, in a
system whose highest polarizability is that one, is therefore written as:

```json
"multipoles": {"elements": [-0.834], "order": 0},
"polarizabilities": {"elements": [0.0, 0.0, 0.0, 0.0,
                                  9.5, 0.0, 0.0, 9.5, 0.0, 9.5],
                     "order": [[0, 0], [0, 1], [1, 1]]}
```

where the four leading zeros are the (0,0) and (0,1) blocks that no site of
that system carries. Where nothing in the system is polarizable, no site
carries the key at all.

Repulsion and dispersion are kept apart, under `repulsion` and `dispersion`,
because they need not agree: a file written by hand may give a site different
parameters for the two, although a site of a PyFraME potential holds one sigma
and one epsilon and so gives both the same. Either may carry an `lj_exponent`
beside `lj_sigma` and `lj_epsilon`, which defaults to 12 for repulsion and 6
for dispersion. A site that has no Lennard-Jones parameters carries neither
key.


Embedding calculations
----------------------

Where the top-level package builds an embedding potential, the
`pyframe.embedding` subpackage consumes one. It evaluates the interactions
between the quantum subsystem and the classical surroundings it is embedded in,
i.e. their energies, the contributions they make to the Fock matrix, and their
geometric derivatives.

The evaluation is split into a module per interaction type, i.e.
`electrostatic_interactions`, `induction_interactions`,
`repulsion_interactions` and `dispersion_interactions`, which operate on the
subsystem and particle objects that `read_input` builds out of a potential. The
numerically intensive parts are delegated to a compiled engine, which is the
C++ extension module that is built when the package is installed.

A snapshot whose embedding potential has been created gives its subsystems
directly, without a JSON file in between:

```python
project.create_embedding_potential(snapshot)
simulation_box, (quantum_subsystem,), (classical_subsystem,) = snapshot.to_subsystems()
```

The subsystems are those `read_input.reader` builds from the JSON file of the
snapshot, and a communicator is passed on as `snapshot.to_subsystems(comm=comm)`. The embedding module is
imported only then, so `import pyframe` does not load the compiled engine.

The notebook
[docs/electrostatic_embedding.ipynb](https://gitlab.com/pyframe-project/pyframe/-/blob/master/docs/electrostatic_embedding.ipynb)
walks through the library.


License
-------

PyFraME is licensed under the Mozilla Public License, v. 2.0. The full text is
in the
[LICENSE](https://gitlab.com/pyframe-project/pyframe/-/blob/master/LICENSE)
file, and `pyframe.print_license()` shows it as well.
