Regions¶
Regions are ‘groups’ of atoms within a ChemicalSystem. Regions are used for some of the AMS driver and engine features, but they can also be useful tools for bookkeeping and manipulations (see also the AMS documentation on Regions).
Here is a simple example showing how to create and use regions:
from scm.base import ChemicalSystem
import numpy as np
# Create a ChemicalSystem and define two regions (one for each water molecule)
mol = ChemicalSystem(
"""
System
Atoms
O -2.676678 0.759085 0.370636 region=water_1
H -3.468900 0.415690 0.814339 region=water_1
H -3.005004 1.433129 -0.246320 region=water_1
O 0.039085 0.303872 1.265871 region=water_2
H -0.874303 0.483757 0.975166 region=water_2
H 0.293563 -0.534617 0.849654 region=water_2
End
End
"""
)
print(mol.region_names) # ["water_1", "water_2"]
print(mol.get_atoms_in_region("water_1")) # [0, 1, 2]
print(mol.is_atom_in_region(0, "water_1")) # True
print(mol.get_regions_of_atom(0)) # ["water_1"]
mol.set_atoms_in_region([0, 1], "water_1") # water_1 is now [0, 1]
mol.add_atoms_to_region([2], "water_1") # water_1 is now [0, 1, 2]
mol.remove_atoms_from_region([2], "water_1") # water_1 is now [0, 1]
mol.remove_region("water_1") # water_1 no longer exists
mol.set_atoms_in_region([0, 1, 2], "water_1") # water_1 is now [0, 1, 2]
# Translate all the atoms in the region "water_2"
for i_atom in mol.get_atoms_in_region("water_2"):
mol.atoms[i_atom].coords += np.array([1, 0, 0])
# Region expressions
print(mol.get_atoms_in_region("water_1 + water_2")) # union [0, 1, 2, 3, 4, 5]
print(mol.get_atoms_in_region("water_1 & water_2")) # intersection []
print(mol.get_atoms_in_region("* - water_2")) # [0, 1, 2]
Get information about current regions¶
- property ChemicalSystem.region_names: List[str]
Returns an ASCIIbetically ordered list of the names of all regions in the system.
Note that there is no such thing as an empty region. Every region contains at least one atom.
Update regions or add new regions¶
The following methods can be used to change the assignment of atoms to regions:
- ChemicalSystem.set_atoms_in_region(atom_indices: ArrayLike, region: str) None
Creates or sets an entire region given the indices of the atoms to be part of the region.
If atom_indices is empty and the region previously existed, it will effectively be deleted.
- ChemicalSystem.add_atoms_to_region(atom_indices: ArrayLike, region: str) None
Adds multiple atoms to a region, given their atom indices.
Region information for a specific atom¶
Valid region names¶
Region names are case-sensitive and may not include certain characters.
Whether a string is a valid name for a region can be checked with the is_valid_region_name method:
- classmethod ChemicalSystem.is_valid_region_name(name: str) bool
Checks if a string is a valid region name.
Valid region names do not contain any of the following characters:
, + - * / \ | $ & ^ < > ( ) [ ] { } " 'Valid region names do not have leading or trailing whitespace. Note that region names are case-sensitive, so “A” and “a” correspond to different regions.
Region expressions¶
Almost all methods for querying region information also support so-called region expressions in addition to plain
region names. Region expressions may include region names, the parentheses (), as well as the operators
+ for union, - for set difference, and & for set intersection. Additionally the wildcard * means all
atoms, while $ stands for the set of selected atoms.
# number of atoms that are either in region A or B (or both)
num_at_AB = mol.num_atoms_in_region("A + B")
# loop over all atoms that are not in the intersection of region A and B
for atidx in mol.get_atoms_in_region("* - (A & B)"):
...