Defining molecular systems

Before diving into any computational chemistry simulation, you need to answer two key questions: What are you studying, and what do you want to do with it? In other words, you have to define the system and the type of calculation that you need to run.

To help with the first part, describing the system, the spycci.systems submodule has you covered.

This submodule provides a set of classes designed to organize and store information about the chemical system you’re working with. Currently the following objects have been implemented:

You can bring these classes into your project using the syntax:

from spycci.systems import System
from spycci.systems import Ensemble

In this section of the guide, we’ll show you how to define molecular systems using these classes, their features and the most useful mode of operation.

Single systems

A System object represents a molecular system containing one or more individual molecules. A System object, as a real molecular system, is defined by a series of attributes:

  • name: the name of the system.

  • geometry: the 3D molecular geometry of the system encoding the position of each atom.

  • charge: the total charge of the system.

  • spin: the spin multiplicity of the system as \(2S+1\), where \(S\) is the total spin quantum number.

  • box_size: for periodic systems the size of the cubic periodic boundary conditions box.

Furthermore each system is associated with a set of properties that are usually set and accessed during computational chemistry calculations and workflows. We will cover these in subsequent chapters of this user guide.

Loading data from an .xyz file

A System object can be initialized in various ways the most common of which is that of providing a molecular structure via a .xyz file. This can be done using the from_xyz() classmethod. In the following example a system object, containing the structure of the water molecule, is loaded from the water.xyz file and some of its properies are printed:

from spycci.systems import System

# Create a System object for a neutral water molecule from a `.xyz` file
mol = System.from_xyz("../example_files/water.xyz")

# Print some attributes
print(f"Name: {mol.name}")
print(f"Charge: {mol.charge}")
print(f"Spin multiplicity: {mol.spin}")
print(f"Number of atoms: {mol.geometry.atomcount}")
print(f"Atoms list: {mol.geometry.atoms}")

Name: water
Charge: 0
Spin multiplicity: 1
Number of atoms: 3
Atoms list: ['O', 'H', 'H']

Please notice how the name of the molecule has been set using the path basename stripping the .xyz file extension and the unspecified attributes have been set to default values:

  • charge: total charge of the system. Default: 0(not charged)

  • spin: total spin multiplicity of the system. Default: 1 (singlet)

  • box_side: length of the cubic simulation box (in \(\mathrm{Å}\)). Default: None (non periodic)

These can of course be specified explicilty in the constructor when needed. For example, a periodic system of side \(\mathrm{18.27Å}\) containing a radical cation can be defined according to:

from spycci.systems import System
my_mol = System.from_xyz(
    "path/to/xyz/water.xyz", 
    charge=1, 
    spin=2,
    box_side=18.27,
)

Creating a System object form a SMILES string

Another convenient way to initialize a System object is that of using Simplified Molecular Input Line Entry System (SMILES) strings. This can be done using the form_smiles class method providing both a name for the system and the SMILES string. As an example the System object representing a ethanol molecule can be derived from is smile CCO according to:

from spycci.systems import System

mol = System.from_smiles("ethanol", "CCO")

# Print some attributes
print(f"Name: {mol.name}")
print(f"Charge: {mol.charge}")
print(f"Spin multiplicity: {mol.spin}")
print(f"Number of atoms: {mol.geometry.atomcount}")
print(f"Atoms list: {mol.geometry.atoms}")

Name: ethanol
Charge: 0
Spin multiplicity: 1
Number of atoms: 9
Atoms list: ['C', 'C', 'O', 'H', 'H', 'H', 'H', 'H', 'H']

The class constructor automatically adds implicit hydrogens to the structure, embeds it in the 3D space and runs a rough geometry optimization using either the MMFFs or the UFF force fields. Additional optional keywords that can be used to configure the SMILES translation process can be found in the System object API.

Saving an loading a System to and form a .json file

Once created, a System object, together with all its computed properties and attributes, can be saved to a standardized .json file. This can be done using the method save_json providing the path to the destination file. The syntax of the command is simple and can be examined in the following example:

from spycci.systems import System

mol = System.from_xyz("water.xyz")

# ... SOME CALCULATION CAN BE DONE HERE ...

mol.save_json("water.json")

The saved .json file contains all the information of the System object and can be used to rebuild the same object at a later time. This can be done using the from_json() classmethod according to:

from spycci.systems import System

mol = System.from_json("water.json")
print(mol)

=========================================================
SYSTEM: water
=========================================================

Number of atoms: 3
Charge: 0
Spin multiplicity: 1

********************** GEOMETRY *************************

Total system mass: 18.0153 amu

----------------------------------------------
 index  atom    x (Å)      y (Å)      z (Å)   
----------------------------------------------
 0       O     0.00000   -0.37937    0.00000  
 1       H     0.77221    0.18968   -0.00000  
 2       H    -0.77221    0.18969    0.00000  
----------------------------------------------

********************** PROPERTIES *************************

Geometry level of theory: XtbInput || method: gfn2 | solvent: None
Electronic level of theory: XtbInput || method: gfn2 | solvent: None
Vibronic level of theory: XtbInput || method: gfn2 | solvent: None

Electronic energy: -5.070544446692 Eh
Gibbs free energy correction G-E(el): 0.002508973529 Eh
Gibbs free energy: -5.0680354731629995 Eh
pKa: pKa object status is NOT SET


MULLIKEN ANALYSIS
----------------------------------------------
 index  atom   charge    spin
----------------------------------------------
 0       O    -0.56500  0.00000  
 1       H    0.28200   0.00000  
 2       H    0.28200   0.00000

Creating a System manually from coordinates

The most basic approach to create a System object is that of manually entering the coordinates of each atom in the molecule. This approach is rarely used by the computational chemistry practitioner but, not involving the use of external files, finds large applications in developing new library features. The workflow is quite simple an involves the creation of a MolecularGeometry object that can be directly fed to the System constructor. As an example the following script can be used to generate a water molecule system directly from coordinates.

from spycci.systems import System
from spycci.core.geometry import MolecularGeometry

WATER = [
    ["O", -5.02534, 1.26595, 0.01097],
    ["H", -4.05210, 1.22164, -0.01263],
    ["H", -5.30240, 0.44124, -0.42809],
]

# Creating a MolecularGeometry onject
water = MolecularGeometry()
for l in WATER:
    water.append(l[0], l[1::])

# Creating the `System` object from the generated geometry
mol = System("water", geometry=water)

# Print some attributes
print(f"Name: {mol.name}")
print(f"Charge: {mol.charge}")
print(f"Spin multiplicity: {mol.spin}")
print(f"Number of atoms: {mol.geometry.atomcount}")
print(f"Atoms list: {mol.geometry.atoms}")

Name: water
Charge: 0
Spin multiplicity: 1
Number of atoms: 3
Atoms list: ['O', 'H', 'H']

Please notice how this approach is the only one relaying directly on the __init__ method of the System class. In fact, all the classmethods explored so far simply generate a MolecularGeometry object that is then fed to the __init__ method.

Ensembles

An Ensemble object is defined as a collection of related Systems, for example the collection of lowest-energy conformers generated by CREST. To initialise an Ensemble object you need to provide an iterable (such as a list) yielding System objects:

from spycci.systems import System
from spycci.systems import Ensemble

my_mol1 = System.from_xyz("path/to/xyz/my_mol1.xyz")
my_mol2 = System.from_xyz("path/to/xyz/my_mol2.xyz")
my_mol3 = System.from_xyz("path/to/xyz/my_mol3.xyz")

my_list = [my_mol1, my_mol2, my_mol3]

my_ens = Ensemble(my_list)

To get a list of all the available properties, please refer to the API.