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/.

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 and Usage example.

Contents

• How to cite
• Contributors
• Requirements
• Installation
• Development
• Testing
• Documentation
• Issues
• Contributing
• Quick start
• Trajectories
• Usage example
• Moving from MolecularSystem
• Parameter sets
• JSON potentials
• Embedding calculations
• 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 represents all versions, and will always resolve to the latest one. A possible BibTeX entry can be found in the CITATION file, and the same entry is printed by pyframe.print_citation(). Alternatively, BibTeX and other formats can be generated by Zenodo.

Contributors

The list of past and current contributors is found here.

Requirements

To use PyFraME you need:

• Python (3.11+)
• NumPy
• SciPy
• scikit-learn
• h5py
• mpi4py
• Fabric
• QCElemental
• MDAnalysis (2.9+)
• daltools and LoProp for Dalton

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 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:

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, which is the default and is used together with LoProp
• Molcas (not tested recently)
• 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:

• pytest
• pytest-timeout
• parse
• PyYAML

Installation

The PyFraME package can be installed from PyPI directly using pip

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.

The entire source including history can be found at GitLab. All releases are also deposited at Zenodo.

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:

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 are needed here as well.

The hooks that pre-commit install sets up run ruff for linting, import sorting and formatting, mypy 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:

pre-commit run --all-files

Testing

The tests are not part of the wheel on PyPI. They come with the source, from GitLab or the sdist, and are run against an editable install as described under Development. To run the full test suite, i.e., the unit tests under pyframe and the integration tests under tests, type

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:

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 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 walks through the embedding library described under Embedding calculations, starting from a UML diagram of it. 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.

For the background of the method and for examples of how PyFraME is used in computational studies, see our tutorial review article.

Issues

Please report issues here.

Contributing

Please take a look at the contribution guide.

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:

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 shows both together.

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

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, 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.

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:

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:

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.

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:

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:

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:

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:

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:

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:

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:

{
  "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:

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:

"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:

{"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:

"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:

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 walks through the library.

License

PyFraME is licensed under the Mozilla Public License, v. 2.0. The full text is in the LICENSE file, and pyframe.print_license() shows it as well.