AMS Settings: Chemical System (Molecule)¶
This example shows how to convert between a PLAMS Molecule object and the text input in the AMS System block.
See also
PLAMS documentation: Molecule handling
To follow along, either
Download
ams_settings_system.py(run as$AMSBIN/amspython ams_settings_system.py).Download
ams_settings_system.ipynb(see also: how to install Jupyterlab in AMS)
Initial imports¶
from scm.plams import *
Elements, coordinates, lattice vectors, and charge¶
Manual molecule definition¶
molecule = Molecule()
molecule.add_atom(Atom(symbol='O', coords=(0,0,0)))
molecule.add_atom(Atom(symbol='H', coords=(1,0,0)))
molecule.add_atom(Atom(symbol='H', coords=(0,1,0)))
To see the input that will be passed to AMS, create an AMSJob and print the input:
print(AMSJob(molecule=molecule).get_input())
system
Atoms
O 0.0000000000 0.0000000000 0.0000000000
H 1.0000000000 0.0000000000 0.0000000000
H 0.0000000000 1.0000000000 0.0000000000
End
End
Lattice vectors: 1D-periodic¶
For periodic systems in 1 dimension, the lattice vector must be along the x direction (with 0 components along y and z)
molecule.lattice = [
[10, 0, 0]
]
print(AMSJob(molecule=molecule).get_input())
system
Atoms
O 0.0000000000 0.0000000000 0.0000000000
H 1.0000000000 0.0000000000 0.0000000000
H 0.0000000000 1.0000000000 0.0000000000
End
Lattice
10.0000000000 0.0000000000 0.0000000000
End
End
Lattice vectors: 2D-periodic¶
For 2 dimensions, the two lattice vectors must lie in the xy plane (with 0 component along z).
molecule.lattice = [
[10, 0, 0],
[0, 11, 0],
]
print(AMSJob(molecule=molecule).get_input())
system
Atoms
O 0.0000000000 0.0000000000 0.0000000000
H 1.0000000000 0.0000000000 0.0000000000
H 0.0000000000 1.0000000000 0.0000000000
End
Lattice
10.0000000000 0.0000000000 0.0000000000
0.0000000000 11.0000000000 0.0000000000
End
End
Lattice vectors: 3D-periodic¶
molecule.lattice = [
[10, 0, 0],
[0, 11, 0],
[-1, 0, 12]
]
print(AMSJob(molecule=molecule).get_input())
system
Atoms
O 0.0000000000 0.0000000000 0.0000000000
H 1.0000000000 0.0000000000 0.0000000000
H 0.0000000000 1.0000000000 0.0000000000
End
Lattice
10.0000000000 0.0000000000 0.0000000000
0.0000000000 11.0000000000 0.0000000000
-1.0000000000 0.0000000000 12.0000000000
End
End
Delete lattice vectors¶
molecule.lattice = []
print(AMSJob(molecule=molecule).get_input())
system
Atoms
O 0.0000000000 0.0000000000 0.0000000000
H 1.0000000000 0.0000000000 0.0000000000
H 0.0000000000 1.0000000000 0.0000000000
End
End
Charge¶
molecule.properties.charge = -1
print(AMSJob(molecule=molecule).get_input())
system
Atoms
O 0.0000000000 0.0000000000 0.0000000000
H 1.0000000000 0.0000000000 0.0000000000
H 0.0000000000 1.0000000000 0.0000000000
End
Charge -1
End
To get the charge of a molecule, use
molecule.properties.get("charge", 0). If the charge is not defined
you will then get 0 as the charge.
my_charge = molecule.properties.get("charge", 0)
print(f"The charge is {my_charge}")
The charge is -1
Unset the charge:
if 'charge' in molecule.properties:
del molecule.properties.charge
my_charge = molecule.properties.get("charge", 0)
print(f"The charge is {my_charge}")
The charge is 0
Atomic properties: masses, regions, force field types …¶
In the AMS system block most atomic properties are given as a suffix at the end of the line.
To access an individual atom, use for example molecule[1], which
corresponds to the first atom. Note that the indexing starts with 1,
unlike normal Python lists that start with 0!
Isotopes (atomic masses)¶
molecule[2].properties.mass = 2.014
print(AMSJob(molecule=molecule).get_input())
system
Atoms
O 0.0000000000 0.0000000000 0.0000000000
H 1.0000000000 0.0000000000 0.0000000000 mass=2.014
H 0.0000000000 1.0000000000 0.0000000000
End
End
Regions¶
Regions are used for example to
set special basis sets on a subset of atoms, or
apply a thermostat in molecular dynamics to only a subset of atoms,
visualize atoms easily in the AMS GUI,
and much more!
Use Python sets to specify regions. In this way, one atom can belong to multiple regions.
molecule[1].properties.region = {"region1"}
molecule[2].properties.region = {"region1"}
molecule[3].properties.region = {"region1", "region2"}
print(AMSJob(molecule=molecule).get_input())
system
Atoms
O 0.0000000000 0.0000000000 0.0000000000 region=region1
H 1.0000000000 0.0000000000 0.0000000000 mass=2.014 region=region1
H 0.0000000000 1.0000000000 0.0000000000 region=region1,region2
End
End
Force field types¶
Some force fields need to know the specific atom type and not just the
chemical element. Use ForceField.Type for this when you use the
ForceField engine:
molecule[1].properties.ForceField.Type = "OW" # these types would depend on what type of force field you use!
molecule[2].properties.ForceField.Type = "HW"
molecule[3].properties.ForceField.Type = "HW"
print(AMSJob(molecule=molecule).get_input())
system
Atoms
O 0.0000000000 0.0000000000 0.0000000000 ForceField.Type=OW region=region1
H 1.0000000000 0.0000000000 0.0000000000 ForceField.Type=HW mass=2.014 region=region1
H 0.0000000000 1.0000000000 0.0000000000 ForceField.Type=HW region=region1,region2
End
End
Delete all atom-specific options¶
Loop over the atoms and set atom.properties to an empty
Settings():
for at in molecule:
at.properties = Settings()
print(AMSJob(molecule=molecule).get_input())
system
Atoms
O 0.0000000000 0.0000000000 0.0000000000
H 1.0000000000 0.0000000000 0.0000000000
H 0.0000000000 1.0000000000 0.0000000000
End
End
Bonds¶
Most methods (DFT, DFTB, ML Potential, ReaxFF) ignore any specified bonds.
When using force fields, you sometimes need to specify the bonds that connect atoms. Some force fields (UFF, GAFF) can automatically guess the correct types of bonds.
So most of the time you do not manually need to specify bonds.
If you need to specify bonds, it is easiest
to handle in the AMS GUI: use File -> Export Coordinates -> .in, and then load the file with
molecule = Molecule("my_file.in")to use the
from_smilesfunction to generate a molecule from SMILES code, for examplemolecule = from_smiles("O")for water.
If you need to add bonds manually in PLAMS you can do it as follows:
molecule.add_bond(molecule[1], molecule[2], order=1.0)
molecule.add_bond(molecule[1], molecule[3], order=1.0)
print(AMSJob(molecule=molecule).get_input())
system
Atoms
O 0.0000000000 0.0000000000 0.0000000000
H 1.0000000000 0.0000000000 0.0000000000
H 0.0000000000 1.0000000000 0.0000000000
End
BondOrders
1 2 1.0
1 2 1.0
1 3 1.0
End
End
Multiple systems¶
Some tasks like NEB (nudged elastic band) require more than 1 system in the input file. This can be accomplished by using a Python dictionary.
In AMS,
the “main system” has no name. It should have the key
""(empty string) in the dictionary.every additional system needs to have a name, that is used as the key in the dictionary.
Let’s first define two Molecule in the normal way:
molecule1 = Molecule()
molecule1.add_atom(Atom(symbol='O', coords=(0,0,0)))
molecule1.add_atom(Atom(symbol='H', coords=(1,0,0)))
molecule1.add_atom(Atom(symbol='H', coords=(0,1,0)))
molecule2 = Molecule()
molecule2.add_atom(Atom(symbol='O', coords=(0,0,0)))
molecule2.add_atom(Atom(symbol='H', coords=(3.33333,0,0)))
molecule2.add_atom(Atom(symbol='H', coords=(0,5.555555,0)))
Then create the mol_dict dictionary:
mol_dict = {
"": molecule1, # main system, empty key (no name)
"final": molecule2 # for NEB, use "final" as the name for the other endpoint
}
Pass the mol_dict as the molecule argument to AMSJob:
print(AMSJob(molecule=mol_dict).get_input())
system
Atoms
O 0.0000000000 0.0000000000 0.0000000000
H 1.0000000000 0.0000000000 0.0000000000
H 0.0000000000 1.0000000000 0.0000000000
End
End
system final
Atoms
O 0.0000000000 0.0000000000 0.0000000000
H 3.3333300000 0.0000000000 0.0000000000
H 0.0000000000 5.5555550000 0.0000000000
End
End
Above we see that the main system is printed just as before. A second
system block “system final” is also added with molecule2.
Complete Python code¶
#!/usr/bin/env amspython
# coding: utf-8
# ## Initial imports
from scm.plams import *
# ## Elements, coordinates, lattice vectors, and charge
# ### Manual molecule definition
molecule = Molecule()
molecule.add_atom(Atom(symbol="O", coords=(0, 0, 0)))
molecule.add_atom(Atom(symbol="H", coords=(1, 0, 0)))
molecule.add_atom(Atom(symbol="H", coords=(0, 1, 0)))
# To see the input that will be passed to AMS, create an AMSJob and print the input:
print(AMSJob(molecule=molecule).get_input())
# ### Lattice vectors: 1D-periodic
#
# For periodic systems in 1 dimension, the lattice vector must be along the x direction (with 0 components along y and z)
molecule.lattice = [[10, 0, 0]]
print(AMSJob(molecule=molecule).get_input())
# ### Lattice vectors: 2D-periodic
#
# For 2 dimensions, the two lattice vectors must lie in the xy plane (with 0 component along z).
molecule.lattice = [
[10, 0, 0],
[0, 11, 0],
]
print(AMSJob(molecule=molecule).get_input())
# ### Lattice vectors: 3D-periodic
molecule.lattice = [[10, 0, 0], [0, 11, 0], [-1, 0, 12]]
print(AMSJob(molecule=molecule).get_input())
# ### Delete lattice vectors
molecule.lattice = []
print(AMSJob(molecule=molecule).get_input())
# ### Charge
molecule.properties.charge = -1
print(AMSJob(molecule=molecule).get_input())
# To get the charge of a molecule, use ``molecule.properties.get("charge", 0)``. If the charge is not defined you will then get 0 as the charge.
my_charge = molecule.properties.get("charge", 0)
print(f"The charge is {my_charge}")
# Unset the charge:
if "charge" in molecule.properties:
del molecule.properties.charge
my_charge = molecule.properties.get("charge", 0)
print(f"The charge is {my_charge}")
# ## Atomic properties: masses, regions, force field types ...
#
# In the AMS system block most atomic properties are given as a suffix at the end of the line.
#
# To access an individual atom, use for example ``molecule[1]``, which corresponds to the first atom. **Note that the indexing starts with 1**, unlike normal Python lists that start with 0!
# ### Isotopes (atomic masses)
molecule[2].properties.mass = 2.014
print(AMSJob(molecule=molecule).get_input())
# ### Regions
#
# Regions are used for example to
#
# * set special basis sets on a subset of atoms, or
# * apply a thermostat in molecular dynamics to only a subset of atoms,
# * visualize atoms easily in the AMS GUI,
# * and much more!
#
# Use Python sets to specify regions. In this way, one atom can belong to multiple regions.
molecule[1].properties.region = {"region1"}
molecule[2].properties.region = {"region1"}
molecule[3].properties.region = {"region1", "region2"}
print(AMSJob(molecule=molecule).get_input())
# ### Force field types
#
# Some force fields need to know the specific atom type and not just the chemical element. Use ``ForceField.Type`` for this when you use the ForceField engine:
molecule[1].properties.ForceField.Type = "OW" # these types would depend on what type of force field you use!
molecule[2].properties.ForceField.Type = "HW"
molecule[3].properties.ForceField.Type = "HW"
print(AMSJob(molecule=molecule).get_input())
# ### Delete all atom-specific options
# Loop over the atoms and set ``atom.properties`` to an empty ``Settings()``:
for at in molecule:
at.properties = Settings()
print(AMSJob(molecule=molecule).get_input())
# ## Bonds
#
# Most methods (DFT, DFTB, ML Potential, ReaxFF) ignore any specified bonds.
#
# When using force fields, you sometimes need to specify the bonds that connect atoms. Some force fields (UFF, GAFF) can automatically guess the correct types of bonds.
#
# So **most of the time you do not manually need to specify bonds**.
#
# If you **need** to specify bonds, it is easiest
#
# * to handle in the AMS GUI: use File -> Export Coordinates -> .in, and then load the file with ``molecule = Molecule("my_file.in")``
# * to use the ``from_smiles`` function to generate a molecule from SMILES code, for example ``molecule = from_smiles("O")`` for water.
#
# If you need to add bonds manually in PLAMS you can do it as follows:
molecule.add_bond(molecule[1], molecule[2], order=1.0)
molecule.add_bond(molecule[1], molecule[3], order=1.0)
print(AMSJob(molecule=molecule).get_input())
# ## Multiple systems
#
# Some tasks like NEB (nudged elastic band) require more than 1 system in the input file. This can be accomplished by using a Python dictionary.
#
# In AMS,
#
# * the "main system" has no name. It should have the key ``""`` (empty string) in the dictionary.
#
# * every additional system needs to have a name, that is used as the key in the dictionary.
#
# Let's first define two ``Molecule`` in the normal way:
molecule1 = Molecule()
molecule1.add_atom(Atom(symbol="O", coords=(0, 0, 0)))
molecule1.add_atom(Atom(symbol="H", coords=(1, 0, 0)))
molecule1.add_atom(Atom(symbol="H", coords=(0, 1, 0)))
molecule2 = Molecule()
molecule2.add_atom(Atom(symbol="O", coords=(0, 0, 0)))
molecule2.add_atom(Atom(symbol="H", coords=(3.33333, 0, 0)))
molecule2.add_atom(Atom(symbol="H", coords=(0, 5.555555, 0)))
# Then create the ``mol_dict`` dictionary:
mol_dict = {
"": molecule1, # main system, empty key (no name)
"final": molecule2, # for NEB, use "final" as the name for the other endpoint
}
# Pass the ``mol_dict`` as the ``molecule`` argument to ``AMSJob``:
print(AMSJob(molecule=mol_dict).get_input())
# Above we see that the main system is printed just as before. A second system block "system final" is also added with ``molecule2``.