Complete guide to storing and converting PLAMS Molecules between Python libraries and file formats
**************************************************************************************************

Technical pymatgen ChemicalSystem ASE

**Requires:** AMS2026 or later

-[ Related examples ]-

* Using AMS as an ASE Calculator

* Charged Systems with the AMS ASE Calculator

* Simulating XRD Patterns from a CIF Structure

-[ Related tutorials ]-

* Python Getting Started tutorial

* Building molecules tutorial

-[ Related documentation ]-

* PLAMS Molecule class

* ChemicalSystem file formats

* ChemicalSystem conversion to/from Python libraries

[image]


Initial imports, files, and helper function
===========================================

   import os
   from pathlib import Path
   from scm.plams import view

   # Make sure to source amsbashrc.sh before launching this example so that the AMSHOME environment variable is set.

   AMSHOME = os.environ["AMSHOME"]
   cif_file = f"{AMSHOME}/atomicdata/Molecules/IZA-Zeolites/ABW.cif"
   xyz_file = f"{AMSHOME}/scripting/scm/params/examples/benchmark/ISOL6/e_13.xyz"
   badxyz_file = f"{AMSHOME}/scripting/scm/plams/unit_tests/xyz/reactant2.xyz"

   assert Path(cif_file).exists(), f"{cif_file} does not exist."
   assert Path(xyz_file).exists(), f"{xyz_file} does not exist."

   output_dir = Path("molecule_formats_outputs")
   output_dir.mkdir(exist_ok=True)


   def head(filename, n: int = 4):
       """Print the first ``n`` lines of a file"""
       with open(filename, "r") as f:
           lines = f.readlines()
           lines = lines[: min(n, len(lines))]
       print("".join(lines))


SMILES
======


Load PLAMS Molecule from SMILES string
--------------------------------------

   from scm.plams import from_smiles, Molecule, plot_molecule

   mol = from_smiles("CCCCO")
   print(f"{type(mol)=}")
   view(mol, padding=-0.5, width=300, height=300, picture_path=output_dir / "example_mol.png")

   type(mol)=<class 'scm.plams.mol.molecule.Molecule'>
   [10.04|11:43:57] Starting Xvfb...
   [10.04|11:43:57] Xvfb started

   [image: image generated from notebook][image]


Convert PLAMS Molecule to SMILES string
---------------------------------------

Note: This requires that bonds are defined in the PLAMS Molecule.

   from scm.plams import to_smiles

   smiles = to_smiles(mol)
   print(smiles)

   CCCCO


.xyz files
==========


Load PLAMS Molecule from .xyz file
----------------------------------

   from scm.plams import Molecule, plot_molecule

   mol = Molecule(xyz_file)
   print(f"{type(mol)=}")
   view(mol, guess_bonds=True, padding=-2, width=500, height=300)

   type(mol)=<class 'scm.plams.mol.molecule.Molecule'>

   [image: image generated from notebook][image]


Write PLAMS Molecule to .xyz file
---------------------------------

   mol.properties.comment = "The comment line (2nd line after the number of atoms)"
   mol.write(output_dir / "out.xyz")

   head(output_dir / "out.xyz")

   26
   The comment line (2nd line after the number of atoms)
            C      -1.47627000      -1.15316000      -0.29279600
            C      -2.86892000      -1.15881000      -0.29469200


.cif files
==========


Load PLAMS Molecule from .cif file
----------------------------------

PLAMS cannot natively read .cif files. Instead, go through another
library, for example ASE or pymatgen.

   from ase.io import read
   from scm.plams import fromASE

   mol: Molecule = fromASE(read(cif_file))
   print(f"{type(mol)=}")
   view(mol, fixed_atom_size=False, padding=-2, width=400, height=300)

   type(mol)=<class 'scm.plams.mol.molecule.Molecule'>


   /Users/ormrodmorley/Documents/code/ams/amshome_fix2026/bin.auto/python3.8/lib/python3.8/site-packages/ase/io/cif.py:401: UserWarning: crystal system 'orthorhombic' is not interpreted for space group Spacegroup(74, setting=1). This may result in wrong setting!
     warnings.warn(

   [image: image generated from notebook][image]


Write PLAMS Molecule to .cif file
---------------------------------

PLAMS cannot natively export to .cif files. Instead, go through
another library, for example ASE or pymatgen.

ASE can be used to write many file formats. See
https://wiki.fysik.dtu.dk/ase/ase/io/io.html

   from scm.plams import toASE

   toASE(mol).write(output_dir / "out.cif")
   head(output_dir / "out.cif")

   data_image0
   _chemical_formula_structural       O16Si8
   _chemical_formula_sum              "O16 Si8"
   _cell_length_a       9.873


AMS .in system block format
===========================


Write PLAMS Molecule to AMS .in system file
-------------------------------------------

   mol.write(output_dir / "ams_system_block.in")
   head(output_dir / "ams_system_block.in")

   System
     Atoms
                 O       4.9365000000       1.3135000000       3.0896710000
                 O       0.0000000000       3.9405000000       7.4746710000


Load PLAMS Molecule from AMS .in system file
--------------------------------------------

   from scm.plams import Molecule

   mol = Molecule(output_dir / "ams_system_block.in")
   view(mol, fixed_atom_size=False, padding=-2, width=400, height=300)

   [image: image generated from notebook][image]


POSCAR/CONTCAR (VASP input format)
==================================


Write PLAMS Molecule to POSCAR/CONTCAR (VASP input format)
----------------------------------------------------------

ASE can be used to write many file formats. See
https://wiki.fysik.dtu.dk/ase/ase/io/io.html

   from scm.plams import toASE

   toASE(mol).write(output_dir / "POSCAR")
   head(output_dir / "POSCAR", 10)

    O Si
    1.0000000000000000
        9.8729999999999993    0.0000000000000000    0.0000000000000000
        0.0000000000000000    5.2539999999999996    0.0000000000000000
        0.0000000000000000    0.0000000000000000    8.7699999999999996
   ... output trimmed ....
     16   8
   Cartesian
     4.9364999999999997  1.3134999999999999  3.0896710000000001
     0.0000000000000000  3.9405000000000001  7.4746709999999998


Load PLAMS Molecule from POSCAR/CONTCAR (VASP input format)
-----------------------------------------------------------

   from scm.plams import fromASE
   from ase.io import read

   mol: Molecule = fromASE(read(output_dir / "POSCAR"))

   print(f"{type(mol)=}")
   view(mol, fixed_atom_size=False, padding=-2, width=400, height=300)

   type(mol)=<class 'scm.plams.mol.molecule.Molecule'>

   [image: image generated from notebook][image]


ASE Atoms Python class
======================


Convert PLAMS Molecule to ASE Atoms
-----------------------------------

   from scm.plams import toASE
   from ase import Atoms
   from ase.visualize.plot import plot_atoms
   import matplotlib.pyplot as plt

   print(f"{type(mol)=}")
   print(f"{mol.get_formula()=}")

   ase_atoms: Atoms = toASE(mol)
   print(f"{type(ase_atoms)=}")
   print(f"{ase_atoms.get_chemical_formula()=}")

   _, ax = plt.subplots(figsize=(2, 2))
   plot_atoms(ase_atoms, rotation="-85x,5y,0z", ax=ax);

   type(mol)=<class 'scm.plams.mol.molecule.Molecule'>
   mol.get_formula()='O16Si8'
   type(ase_atoms)=<class 'ase.atoms.Atoms'>
   ase_atoms.get_chemical_formula()='O16Si8'

   [image: image generated from notebook][image]


Convert ASE Atoms to PLAMS Molecule
-----------------------------------

   from scm.plams import fromASE, plot_molecule, Molecule

   mol: Molecule = fromASE(ase_atoms)
   print(f"{type(mol)=}")
   view(
       mol,
       direction="tilt_z",
       show_lattice_vectors=True,
       padding=-1,
       fixed_atom_size=False,
       width=400,
       height=400,
   )

   type(mol)=<class 'scm.plams.mol.molecule.Molecule'>

   [image: image generated from notebook][image]


RDKit Mol Python class
======================


Convert PLAMS Molecule to RDKit Mol
-----------------------------------

   from scm.plams import to_rdmol, Molecule
   from rdkit.Chem import Draw
   from rdkit.Chem.Draw import IPythonConsole

   IPythonConsole.ipython_useSVG = True
   IPythonConsole.molSize = 250, 250

   plams_mol = Molecule(xyz_file)
   # guess bonds, the bonds will be included in the RDKit molecule
   plams_mol.guess_bonds()

   rdkit_mol = to_rdmol(plams_mol)
   print(f"{type(rdkit_mol)=}")
   rdkit_mol

   type(rdkit_mol)=<class 'rdkit.Chem.rdchem.Mol'>

   [image: image generated from notebook][image]


Convert RDKit Mol to PLAMS Molecule
-----------------------------------

   from scm.plams import from_rdmol, plot_molecule, Molecule

   mol: Molecule = from_rdmol(rdkit_mol)

   print(f"{type(rdkit_mol)=}")
   print(f"{type(mol)=}")
   view(mol, padding=-2, width=500, height=300)

   type(rdkit_mol)=<class 'rdkit.Chem.rdchem.Mol'>
   type(mol)=<class 'scm.plams.mol.molecule.Molecule'>

   [image: image generated from notebook][image]


Convert problematic PLAMS Molecule to RDKit Mol
-----------------------------------------------

   mol = Molecule(badxyz_file)
   mol.guess_bonds()
   view(mol, padding=-0.5, width=300, height=300)

   [image: image generated from notebook][image]

This molecule will fail to convert to an RDKit Mol object, because
RDKit does not like the AMS assignment of double bonds.

   try:
       rdkit_mol = to_rdmol(mol)
   except ValueError as exc:
       print("Failed to convert")

   [10.04|11:44:05] RDKit Sanitization Error.
   [10.04|11:44:05] Most likely this is a problem with the assigned bond orders: Use chemical insight to adjust them.
   [10.04|11:44:05] Note that the atom indices below start at zero, while the AMS-GUI indices start at 1.
   Failed to convert


   RDKit ERROR: [11:44:05] Can't kekulize mol.  Unkekulized atoms: 10 11 12 13 14
   [11:44:05] Can't kekulize mol.  Unkekulized atoms: 10 11 12 13 14

   RDKit ERROR:

The problem can be fixed by passing the argument "presanitize" to the
"to_rdmol" function.

   rdkit_mol = to_rdmol(mol, presanitize=True)
   rdkit_mol

   RDKit ERROR: [11:44:05] Can't kekulize mol.  Unkekulized atoms: 10 11 12 13 14
   RDKit ERROR:
   [11:44:05] Can't kekulize mol.  Unkekulized atoms: 10 11 12 13 14

   [image: image generated from notebook][image]


SCM Base ChemicalSystem Python class
====================================


Convert PLAMS Molecule to ChemicalSystem
----------------------------------------

   from scm.utils.conversions import plams_molecule_to_chemsys, chemsys_to_plams_molecule
   from scm.plams import Molecule
   from scm.base import ChemicalSystem

   mol = Molecule(xyz_file)
   chemsys = plams_molecule_to_chemsys(mol)
   print(f"{type(chemsys)=}")
   print(chemsys)

   type(chemsys)=<class 'scm.base._internal.ChemicalSystem'>
   System
      Atoms
         C -1.47627 -1.15316 -0.292796
         C -2.86892 -1.15881 -0.294692
   ... output trimmed ....
         H 5.434540000000001 0.731266 -0.466878
         H -5.434540000000001 -0.731266 -0.466878
         H -5.43534 0.887836 -0.024637
      End
   End


Convert ChemicalSystem to PLAMS Molecule
----------------------------------------

   from scm.utils.conversions import plams_molecule_to_chemsys, chemsys_to_plams_molecule
   from scm.plams import Molecule
   from scm.base import ChemicalSystem

   mol = chemsys_to_plams_molecule(chemsys)
   print(f"{type(chemsys)=}")
   print(f"{type(mol)=}")
   view(mol, guess_bonds=True, padding=-2, width=500, height=300)

   type(chemsys)=<class 'scm.base._internal.ChemicalSystem'>
   type(mol)=<class 'scm.plams.mol.molecule.Molecule'>

   [image: image generated from notebook][image]


pymatgen Structure and Molecule Python classes
==============================================


Convert PLAMS Molecule to pymatgen Structure (periodic)
-------------------------------------------------------

There is no builtin converter between PLAMS Molecule and pymatgen
Structure (periodic crystal). Instead, you need to go through the ASE
interface to both packages:

   from pymatgen.core.structure import Structure
   from pymatgen.io.ase import AseAtomsAdaptor
   import scm.plams
   from scm.plams import fromASE, toASE, Molecule
   from ase.io import read


   def convert_plams_molecule_to_pymatgen_structure(mol: Molecule) -> Structure:
       return AseAtomsAdaptor().get_structure(toASE(mol))


   mol: scm.plams.Molecule = fromASE(read(cif_file))

   pymatgen_structure: Structure = convert_plams_molecule_to_pymatgen_structure(mol)

   print(f"{type(mol)=}")
   print(f"{type(pymatgen_structure)=}")
   print(pymatgen_structure)

   type(mol)=<class 'scm.plams.mol.molecule.Molecule'>
   type(pymatgen_structure)=<class 'pymatgen.core.structure.Structure'>
   Full Formula (Si8 O16)
   Reduced Formula: SiO2
   abc   :   9.873000   5.254000   8.770000
   ... output trimmed ....
    19  Si    0.1572  0.75  0.9007
    20  Si    0.6572  0.75  0.5993
    21  Si    0.1572  0.25  0.0993
    22  Si    0.3428  0.75  0.5993
    23  Si    0.8428  0.25  0.0993


   /Users/ormrodmorley/Documents/code/ams/amshome_fix2026/bin.auto/python3.8/lib/python3.8/site-packages/ase/io/cif.py:401: UserWarning: crystal system 'orthorhombic' is not interpreted for space group Spacegroup(74, setting=1). This may result in wrong setting!
     warnings.warn(


Convert pymatgen Structure (periodic) to PLAMS Molecule
-------------------------------------------------------

Go through the ASE interface:

   from pymatgen.io.ase import AseAtomsAdaptor
   from pymatgen.core.structure import Structure
   from scm.plams import fromASE
   from scm.plams import Molecule


   def pymatgen_structure_to_plams_molecule(pymatgen_structure: Structure) -> Molecule:
       return fromASE(AseAtomsAdaptor().get_atoms(pymatgen_structure))


   print(f"{type(pymatgen_structure)=}")

   mol = pymatgen_structure_to_plams_molecule(pymatgen_structure)
   print(f"{type(mol)=}")

   type(pymatgen_structure)=<class 'pymatgen.core.structure.Structure'>
   type(mol)=<class 'scm.plams.mol.molecule.Molecule'>


Convert PLAMS Molecule to pymatgen Molecule (non-periodic)
----------------------------------------------------------

pymatgen has a special "Molecule" class for non-periodic systems. In
PLAMS, the "Molecule" class is used for both periodic and non-periodic
systems.

   import pymatgen.core.structure
   import scm.plams
   from pymatgen.io.ase import AseAtomsAdaptor
   from scm.plams import toASE


   def convert_plams_molecule_to_pymatgen_molecule(
       mol: scm.plams.Molecule,
   ) -> pymatgen.core.structure.Molecule:
       return AseAtomsAdaptor().get_molecule(toASE(mol))


   plams_molecule = scm.plams.Molecule(xyz_file)

   pymatgen_molecule: pymatgen.core.structure.Molecule = convert_plams_molecule_to_pymatgen_molecule(
       plams_molecule
   )

   print(f"{type(plams_molecule)=}")
   print(f"{type(pymatgen_molecule)=}")
   print(pymatgen_molecule)

   type(plams_molecule)=<class 'scm.plams.mol.molecule.Molecule'>
   type(pymatgen_molecule)=<class 'pymatgen.core.structure.Molecule'>
   Full Formula (H12 C12 N2)
   Reduced Formula: H6C6N
   Charge = 0, Spin Mult = 1
   ... output trimmed ....
   21 H     3.407760     2.066590    -0.565539
   22 H     5.435340    -0.887836    -0.024637
   23 H     5.434540     0.731266    -0.466878
   24 H    -5.434540    -0.731266    -0.466878
   25 H    -5.435340     0.887836    -0.024637


Convert pymatgen Molecule (non-periodic) to PLAMS Molecule
----------------------------------------------------------

   from pymatgen.io.ase import AseAtomsAdaptor
   import pymatgen.core.structure
   from scm.plams import fromASE
   from scm.plams import Molecule


   def pymatgen_molecule_to_plams_molecule(
       pymatgen_molecule: pymatgen.core.structure.Molecule,
   ) -> scm.plams.Molecule:
       return fromASE(AseAtomsAdaptor().get_atoms(pymatgen_molecule))


   print(f"{type(pymatgen_molecule)=}")

   mol = pymatgen_molecule_to_plams_molecule(pymatgen_molecule)
   print(f"{type(mol)=}")
   view(
       mol,
       guess_bonds=True,
       padding=-2,
       width=500,
       height=300,
       picture_path=output_dir / "example_mol.png",
   )

   type(pymatgen_molecule)=<class 'pymatgen.core.structure.Molecule'>
   type(mol)=<class 'scm.plams.mol.molecule.Molecule'>

   [image: image generated from notebook][image]


See also
========


Related examples
----------------

* Using AMS as an ASE Calculator

* Charged Systems with the AMS ASE Calculator

* Simulating XRD Patterns from a CIF Structure


Related tutorials
-----------------

* Python Getting Started tutorial

* Building molecules tutorial


Related documentation
---------------------

* PLAMS Molecule class

* ChemicalSystem file formats

* ChemicalSystem conversion to/from Python libraries


Python Script
=============

   #!/usr/bin/env python
   # coding: utf-8

   # ## Initial imports, files, and helper function

   import os
   from pathlib import Path
   from scm.plams import view

   # Make sure to source amsbashrc.sh before launching this example so that the AMSHOME environment variable is set.

   AMSHOME = os.environ["AMSHOME"]
   cif_file = f"{AMSHOME}/atomicdata/Molecules/IZA-Zeolites/ABW.cif"
   xyz_file = f"{AMSHOME}/scripting/scm/params/examples/benchmark/ISOL6/e_13.xyz"
   badxyz_file = f"{AMSHOME}/scripting/scm/plams/unit_tests/xyz/reactant2.xyz"

   assert Path(cif_file).exists(), f"{cif_file} does not exist."
   assert Path(xyz_file).exists(), f"{xyz_file} does not exist."

   output_dir = Path("molecule_formats_outputs")
   output_dir.mkdir(exist_ok=True)


   def head(filename, n: int = 4):
       """Print the first ``n`` lines of a file"""
       with open(filename, "r") as f:
           lines = f.readlines()
           lines = lines[: min(n, len(lines))]
       print("".join(lines))


   # ## SMILES

   # ### Load PLAMS Molecule from SMILES string

   from scm.plams import from_smiles, Molecule, plot_molecule

   mol = from_smiles("CCCCO")
   print(f"{type(mol)=}")
   view(mol, padding=-0.5, width=300, height=300, picture_path=output_dir / "example_mol.png")


   # ### Convert PLAMS Molecule to SMILES string
   #
   # Note: This requires that bonds are defined in the PLAMS Molecule.

   from scm.plams import to_smiles

   smiles = to_smiles(mol)
   print(smiles)


   # ## .xyz files

   # ### Load PLAMS Molecule from .xyz file

   from scm.plams import Molecule, plot_molecule

   mol = Molecule(xyz_file)
   print(f"{type(mol)=}")
   view(mol, guess_bonds=True, padding=-2, width=500, height=300, picture_path="picture1.png")


   # ### Write PLAMS Molecule to .xyz file

   mol.properties.comment = "The comment line (2nd line after the number of atoms)"
   mol.write(output_dir / "out.xyz")


   head(output_dir / "out.xyz")


   # ## .cif files

   # ### Load PLAMS Molecule from .cif file
   #
   # PLAMS cannot natively read .cif files. Instead, go through another library, for example ASE or pymatgen.

   from ase.io import read
   from scm.plams import fromASE

   mol: Molecule = fromASE(read(cif_file))
   print(f"{type(mol)=}")
   view(mol, fixed_atom_size=False, padding=-2, width=400, height=300, picture_path="picture2.png")


   # ### Write PLAMS Molecule to .cif file
   #
   # PLAMS cannot natively export to .cif files. Instead, go through another library, for example ASE or pymatgen.
   #
   # ASE can be used to write many file formats. See https://wiki.fysik.dtu.dk/ase/ase/io/io.html

   from scm.plams import toASE

   toASE(mol).write(output_dir / "out.cif")
   head(output_dir / "out.cif")


   # ## AMS .in system block format
   #
   # ### Write PLAMS Molecule to AMS .in system file

   mol.write(output_dir / "ams_system_block.in")
   head(output_dir / "ams_system_block.in")


   # ### Load PLAMS Molecule from AMS .in system file

   from scm.plams import Molecule

   mol = Molecule(output_dir / "ams_system_block.in")
   view(mol, fixed_atom_size=False, padding=-2, width=400, height=300, picture_path="picture3.png")


   # ## POSCAR/CONTCAR (VASP input format)

   # ### Write PLAMS Molecule to POSCAR/CONTCAR (VASP input format)
   #
   # ASE can be used to write many file formats. See https://wiki.fysik.dtu.dk/ase/ase/io/io.html

   from scm.plams import toASE

   toASE(mol).write(output_dir / "POSCAR")
   head(output_dir / "POSCAR", 10)


   # ### Load PLAMS Molecule from POSCAR/CONTCAR (VASP input format)

   from scm.plams import fromASE
   from ase.io import read

   mol: Molecule = fromASE(read(output_dir / "POSCAR"))

   print(f"{type(mol)=}")
   view(mol, fixed_atom_size=False, padding=-2, width=400, height=300, picture_path="picture4.png")


   # ## ASE Atoms Python class

   # ### Convert PLAMS Molecule to ASE Atoms

   from scm.plams import toASE
   from ase import Atoms
   from ase.visualize.plot import plot_atoms
   import matplotlib.pyplot as plt

   print(f"{type(mol)=}")
   print(f"{mol.get_formula()=}")

   ase_atoms: Atoms = toASE(mol)
   print(f"{type(ase_atoms)=}")
   print(f"{ase_atoms.get_chemical_formula()=}")

   _, ax = plt.subplots(figsize=(2, 2))
   plot_atoms(ase_atoms, rotation="-85x,5y,0z", ax=ax)
   # ### Convert ASE Atoms to PLAMS Molecule

   from scm.plams import fromASE, plot_molecule, Molecule

   mol: Molecule = fromASE(ase_atoms)
   print(f"{type(mol)=}")
   view(
       mol,
       direction="tilt_z",
       show_lattice_vectors=True,
       padding=-1,
       fixed_atom_size=False,
       width=400,
       height=400,
       picture_path="picture5.png",
   )


   # ## RDKit Mol Python class

   # ### Convert PLAMS Molecule to RDKit Mol

   from scm.plams import to_rdmol, Molecule
   from rdkit.Chem import Draw
   from rdkit.Chem.Draw import IPythonConsole

   IPythonConsole.ipython_useSVG = True
   IPythonConsole.molSize = 250, 250

   plams_mol = Molecule(xyz_file)
   # guess bonds, the bonds will be included in the RDKit molecule
   plams_mol.guess_bonds()

   rdkit_mol = to_rdmol(plams_mol)
   print(f"{type(rdkit_mol)=}")
   rdkit_mol


   # ### Convert RDKit Mol to PLAMS Molecule

   from scm.plams import from_rdmol, plot_molecule, Molecule

   mol: Molecule = from_rdmol(rdkit_mol)

   print(f"{type(rdkit_mol)=}")
   print(f"{type(mol)=}")
   view(mol, padding=-2, width=500, height=300, picture_path="picture6.png")


   # ### Convert problematic PLAMS Molecule to RDKit Mol

   mol = Molecule(badxyz_file)
   mol.guess_bonds()
   view(mol, padding=-0.5, width=300, height=300, picture_path="picture7.png")


   # This molecule will fail to convert to an RDKit Mol object, because RDKit does not like the AMS assignment of double bonds.

   try:
       rdkit_mol = to_rdmol(mol)
   except ValueError as exc:
       print("Failed to convert")


   # The problem can be fixed by passing the argument `presanitize` to the `to_rdmol` function.

   rdkit_mol = to_rdmol(mol, presanitize=True)
   rdkit_mol


   # ## SCM Base ChemicalSystem Python class
   #
   # ### Convert PLAMS Molecule to ChemicalSystem

   from scm.utils.conversions import plams_molecule_to_chemsys, chemsys_to_plams_molecule
   from scm.plams import Molecule
   from scm.base import ChemicalSystem

   mol = Molecule(xyz_file)
   chemsys = plams_molecule_to_chemsys(mol)
   print(f"{type(chemsys)=}")
   print(chemsys)


   # ### Convert ChemicalSystem to PLAMS Molecule

   from scm.utils.conversions import plams_molecule_to_chemsys, chemsys_to_plams_molecule
   from scm.plams import Molecule
   from scm.base import ChemicalSystem

   mol = chemsys_to_plams_molecule(chemsys)
   print(f"{type(chemsys)=}")
   print(f"{type(mol)=}")
   view(mol, guess_bonds=True, padding=-2, width=500, height=300, picture_path="picture8.png")


   # ## pymatgen Structure and Molecule Python classes

   # ### Convert PLAMS Molecule to pymatgen Structure (periodic)
   #
   # There is no builtin converter between PLAMS Molecule and pymatgen Structure (periodic crystal). Instead, you need to go through the ASE interface to both packages:

   from pymatgen.core.structure import Structure
   from pymatgen.io.ase import AseAtomsAdaptor
   import scm.plams
   from scm.plams import fromASE, toASE, Molecule
   from ase.io import read


   def convert_plams_molecule_to_pymatgen_structure(mol: Molecule) -> Structure:
       return AseAtomsAdaptor().get_structure(toASE(mol))


   mol: scm.plams.Molecule = fromASE(read(cif_file))

   pymatgen_structure: Structure = convert_plams_molecule_to_pymatgen_structure(mol)

   print(f"{type(mol)=}")
   print(f"{type(pymatgen_structure)=}")
   print(pymatgen_structure)


   # ### Convert pymatgen Structure (periodic) to PLAMS Molecule
   #
   # Go through the ASE interface:

   from pymatgen.io.ase import AseAtomsAdaptor
   from pymatgen.core.structure import Structure
   from scm.plams import fromASE
   from scm.plams import Molecule


   def pymatgen_structure_to_plams_molecule(pymatgen_structure: Structure) -> Molecule:
       return fromASE(AseAtomsAdaptor().get_atoms(pymatgen_structure))


   print(f"{type(pymatgen_structure)=}")

   mol = pymatgen_structure_to_plams_molecule(pymatgen_structure)
   print(f"{type(mol)=}")


   # ### Convert PLAMS Molecule to pymatgen Molecule (non-periodic)
   #
   # pymatgen has a special ``Molecule`` class for non-periodic systems. In PLAMS, the ``Molecule`` class is used for both periodic and non-periodic systems.

   import pymatgen.core.structure
   import scm.plams
   from pymatgen.io.ase import AseAtomsAdaptor
   from scm.plams import toASE


   def convert_plams_molecule_to_pymatgen_molecule(
       mol: scm.plams.Molecule,
   ) -> pymatgen.core.structure.Molecule:
       return AseAtomsAdaptor().get_molecule(toASE(mol))


   plams_molecule = scm.plams.Molecule(xyz_file)

   pymatgen_molecule: pymatgen.core.structure.Molecule = convert_plams_molecule_to_pymatgen_molecule(plams_molecule)

   print(f"{type(plams_molecule)=}")
   print(f"{type(pymatgen_molecule)=}")
   print(pymatgen_molecule)


   # ### Convert pymatgen Molecule (non-periodic) to PLAMS Molecule

   from pymatgen.io.ase import AseAtomsAdaptor
   import pymatgen.core.structure
   from scm.plams import fromASE
   from scm.plams import Molecule


   def pymatgen_molecule_to_plams_molecule(
       pymatgen_molecule: pymatgen.core.structure.Molecule,
   ) -> scm.plams.Molecule:
       return fromASE(AseAtomsAdaptor().get_atoms(pymatgen_molecule))


   print(f"{type(pymatgen_molecule)=}")

   mol = pymatgen_molecule_to_plams_molecule(pymatgen_molecule)
   print(f"{type(mol)=}")
   view(
       mol,
       guess_bonds=True,
       padding=-2,
       width=500,
       height=300,
       picture_path=output_dir / "example_mol.png",
   )
