Packmol example¶

This example illustrates various ways of using the packmol interface for constructing liquid or gas mixtures or solid/liquid interfaces.

Note: This example requires AMS2023 or later.

Initial imports¶

from scm.plams import *
from ase.visualize.plot import plot_atoms
from ase.build import fcc111, bulk
import matplotlib.pyplot as plt


Helper functions¶

def printsummary(mol, details=None):
if details:
density = details["density"]
else:
density = mol.get_density() * 1e-3
s = f"{len(mol)} atoms, density = {density:.3f} g/cm^3"
if mol.lattice:
s += f", box = {mol.lattice[0][0]:.3f}, {mol.lattice[1][1]:.3f}, {mol.lattice[2][2]:.3f}"
s += f", formula = {mol.get_formula()}"
if details:
s += f'\n#added molecules per species: {details["n_molecules"]}, mole fractions: {details["mole_fractions"]}'
print(s)

def show(mol, figsize=None, **kwargs):
"""Show a molecule in a Jupyter notebook"""
plt.figure(figsize=figsize or (2, 2))
plt.axis("off")
plot_atoms(toASE(mol), **kwargs)


Liquid water (fluid with 1 component)¶

First, create the gasphase molecule:

water = from_smiles("O")
show(water)

print("pure liquid from approximate number of atoms and exact density (in g/cm^3), cubic box with auto-determined size")
out = packmol(water, n_atoms=194, density=1.0)
printsummary(out)
out.write("water-1.xyz")
show(out)

pure liquid from approximate number of atoms and exact density (in g/cm^3), cubic box with auto-determined size
195 atoms, density = 1.000 g/cm^3, box = 12.482, 12.482, 12.482, formula = H130O65

print("pure liquid from approximate density (in g/cm^3) and an orthorhombic box")
out = packmol(water, density=1.0, box_bounds=[0.0, 0.0, 0.0, 8.0, 12.0, 14.0])
printsummary(out)
out.write("water-2.xyz")
show(out)

pure liquid from approximate density (in g/cm^3) and an orthorhombic box
135 atoms, density = 1.002 g/cm^3, box = 8.000, 12.000, 14.000, formula = H90O45

print("pure liquid with explicit number of molecules and exact density")
out = packmol(water, n_molecules=64, density=1.0)
printsummary(out)
out.write("water-3.xyz")
show(out)

pure liquid with explicit number of molecules and exact density
192 atoms, density = 1.000 g/cm^3, box = 12.417, 12.417, 12.417, formula = H128O64

print("pure liquid with explicit number of molecules and box")
out = packmol(water, n_molecules=64, box_bounds=[0.0, 0.0, 0.0, 12.0, 13.0, 14.0])
printsummary(out)
out.write("water-4.xyz")
show(out)

pure liquid with explicit number of molecules and box
192 atoms, density = 0.877 g/cm^3, box = 12.000, 13.000, 14.000, formula = H128O64

print("water-5.xyz: pure liquid in non-orthorhombic box (requires AMS2022 or later)")
# first place the molecules in a cuboid surrounding the desired lattice
# then gradually change into the desired lattice using refine_lattice()
# note that the molecules may become distorted by this procedure
lattice = [[10.0, 2.0, -1.0], [-5.0, 8.0, 0.0], [0.0, -2.0, 11.0]]
temp_out = packmol(
water,
n_molecules=32,
box_bounds=[
0,
0,
0,
max(lattice[i][0] for i in range(3)) - min(lattice[i][0] for i in range(3)),
max(lattice[i][1] for i in range(3)) - min(lattice[i][1] for i in range(3)),
max(lattice[i][2] for i in range(3)) - min(lattice[i][2] for i in range(3)),
],
)
out = refine_lattice(temp_out, lattice=lattice)
if out is not None:
out.write("water-5.xyz")
print(
"Top: system in surrounding orthorhombic box before calling refine_lattice(). Bottom: System in non-orthorhombic box after calling refine_lattice()"
)
show(temp_out)
show(out)

water-5.xyz: pure liquid in non-orthorhombic box (requires AMS2022 or later)
PLAMS working folder: /path/plams_workdir
Top: system in surrounding orthorhombic box before calling refine_lattice(). Bottom: System in non-orthorhombic box after calling refine_lattice()


Water-acetonitrile mixture (fluid with 2 or more components)¶

Let’s also create a single acetonitrile molecule:

acetonitrile = from_smiles("CC#N")
show(acetonitrile)


Set the desired mole fractions and density. Here, the density is calculated as the weighted average of water (1.0 g/cm^3) and acetonitrile (0.76 g/cm^3) densities, but you could use any other density.

# MIXTURES
x_water = 0.666  # mole fraction
x_acetonitrile = 1 - x_water  # mole fraction
density = (x_water * 1.0 + x_acetonitrile * 0.76) / (
x_water + x_acetonitrile
)  # weighted average of pure component densities

print("MIXTURES")
print(f"x_water = {x_water:.3f}")
print(f"x_acetonitrile = {x_acetonitrile:.3f}")
print(f"target density = {density:.3f} g/cm^3")

MIXTURES
x_water = 0.666
x_acetonitrile = 0.334
target density = 0.920 g/cm^3


By setting return_details=True, you can get information about the mole fractions of the returned system. They may not exactly match the mole fractions you put in.

print(
"2-1 water-acetonitrile from approximate number of atoms and exact density (in g/cm^3), cubic box with auto-determined size"
)
out, details = packmol(
molecules=[water, acetonitrile],
mole_fractions=[x_water, x_acetonitrile],
n_atoms=200,
density=density,
return_details=True,
)
printsummary(out, details)
out.write("water-acetonitrile-1.xyz")
show(out)

2-1 water-acetonitrile from approximate number of atoms and exact density (in g/cm^3), cubic box with auto-determined size
201 atoms, density = 0.920 g/cm^3, box = 13.263, 13.263, 13.263, formula = C34H117N17O33
#added molecules per species: [33, 17], mole fractions: [0.66, 0.34]


The details is a dictionary as follows:

for k, v in details.items():
print(f"{k}: {v}")

n_molecules: [33, 17]
mole_fractions: [0.66, 0.34]
n_atoms: 201
molecule_type_indices: [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
molecule_indices: [0, 0, 0, 1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4, 5, 5, 5, 6, 6, 6, 7, 7, 7, 8, 8, 8, 9, 9, 9, 10, 10, 10, 11, 11, 11, 12, 12, 12, 13, 13, 13, 14, 14, 14, 15, 15, 15, 16, 16, 16, 17, 17, 17, 18, 18, 18, 19, 19, 19, 20, 20, 20, 21, 21, 21, 22, 22, 22, 23, 23, 23, 24, 24, 24, 25, 25, 25, 26, 26, 26, 27, 27, 27, 28, 28, 28, 29, 29, 29, 30, 30, 30, 31, 31, 31, 32, 32, 32, 33, 33, 33, 33, 33, 33, 34, 34, 34, 34, 34, 34, 35, 35, 35, 35, 35, 35, 36, 36, 36, 36, 36, 36, 37, 37, 37, 37, 37, 37, 38, 38, 38, 38, 38, 38, 39, 39, 39, 39, 39, 39, 40, 40, 40, 40, 40, 40, 41, 41, 41, 41, 41, 41, 42, 42, 42, 42, 42, 42, 43, 43, 43, 43, 43, 43, 44, 44, 44, 44, 44, 44, 45, 45, 45, 45, 45, 45, 46, 46, 46, 46, 46, 46, 47, 47, 47, 47, 47, 47, 48, 48, 48, 48, 48, 48, 49, 49, 49, 49, 49, 49]
atom_indices_in_molecule: [0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5]
volume: 2333.0853879652004
density: 0.9198400000000004

print("2-1 water-acetonitrile from approximate density (in g/cm^3) and box bounds")
out, details = packmol(
molecules=[water, acetonitrile],
mole_fractions=[x_water, x_acetonitrile],
box_bounds=[0, 0, 0, 13.2, 13.2, 13.2],
density=density,
return_details=True,
)
printsummary(out, details)
out.write("water-acetonitrile-2.xyz")
show(out)

2-1 water-acetonitrile from approximate density (in g/cm^3) and box bounds
201 atoms, density = 0.933 g/cm^3, box = 13.200, 13.200, 13.200, formula = C34H117N17O33
#added molecules per species: [33, 17], mole fractions: [0.66, 0.34]

print("2-1 water-acetonitrile from explicit number of molecules and density, cubic box with auto-determined size")
out, details = packmol(molecules=[water, acetonitrile], n_molecules=[32, 16], density=density, return_details=True)
printsummary(out, details)
out.write("water-acetonitrile-3.xyz")
show(out)

2-1 water-acetonitrile from explicit number of molecules and density, cubic box with auto-determined size
192 atoms, density = 0.920 g/cm^3, box = 13.058, 13.058, 13.058, formula = C32H112N16O32
#added molecules per species: [32, 16], mole fractions: [0.6666666666666666, 0.3333333333333333]

print("2-1 water-acetonitrile from explicit number of molecules and box")
out = packmol(molecules=[water, acetonitrile], n_molecules=[32, 16], box_bounds=[0, 0, 0, 13.2, 13.2, 13.2])
printsummary(out)
out.write("water-acetonitrile-4.xyz")
show(out)

2-1 water-acetonitrile from explicit number of molecules and box
192 atoms, density = 0.890 g/cm^3, box = 13.200, 13.200, 13.200, formula = C32H112N16O32


Pack inside sphere¶

Set sphere=True to pack in a sphere (non-periodic) instead of in a periodic box. The sphere will be centered near the origin.

print("water in a sphere from exact density and number of molecules")
out, details = packmol(molecules=[water], n_molecules=[100], density=1.0, return_details=True, sphere=True)
printsummary(out, details)
print(f"Center of mass xyz (ang): {out.get_center_of_mass()}")
out.write("water-sphere.xyz")
show(out)

water in a sphere from exact density and number of molecules
300 atoms, density = 1.000 g/cm^3, formula = H200O100
#added molecules per species: [100], mole fractions: [1.0]
Center of mass xyz (ang): (0.35956557054572336, 0.23551764976716527, -0.8914888983730765)

print(
"2-1 water-acetonitrile in a sphere from exact density (in g/cm^3) and approximate number of atoms and mole fractions"
)
out, details = packmol(
molecules=[water, acetonitrile],
mole_fractions=[x_water, x_acetonitrile],
n_atoms=500,
density=density,
return_details=True,
sphere=True,
)
printsummary(out, details)
out.write("water-acetonitrile-sphere.xyz")
show(out)

2-1 water-acetonitrile in a sphere from exact density (in g/cm^3) and approximate number of atoms and mole fractions
501 atoms, density = 0.920 g/cm^3, formula = C84H292N42O83
#added molecules per species: [83, 42], mole fractions: [0.664, 0.336]


Packing ions, total system charge¶

The total system charge will be sum of the charges of the constituent molecules.

In PLAMS, molecule.properties.charge specifies the charge:

ammonium = from_smiles("[NH4+]")  # ammonia.properties.charge == +1
chloride = from_smiles("[Cl-]")  # chloride.properties.charge == -1
print("3 water molecules, 3 ammonium, 1 chloride (non-periodic)")
print("Initial charges:")
print(f"Water: {water.properties.get('charge', 0)}")
print(f"Ammonia: {ammonium.properties.get('charge', 0)}")
print(f"Chloride: {chloride.properties.get('charge', 0)}")
out = packmol(molecules=[water, ammonium, chloride], n_molecules=[3, 3, 1], density=0.4, sphere=True)
tot_charge = out.properties.get("charge", 0)
print(f"Total charge of packmol-generated system: {tot_charge}")
out.write("water-ammonium-chloride.xyz")
show(out)

3 water molecules, 3 ammonium, 1 chloride (non-periodic)
Initial charges:
Water: 0
Ammonia: 1
Chloride: -1
Total charge of packmol-generated system: 2


Microsolvation¶

packmol_microsolvation can create a microsolvation sphere around a solute.

out = packmol_microsolvation(solute=acetonitrile, solvent=water, density=1.5, threshold=4.0)
# for microsolvation it's a good idea to have a higher density than normal to get enough solvent molecules
print(f"Microsolvated structure: {len(out)} atoms.")
out.write("acetonitrile-microsolvated.xyz")

figsize = (3, 3)
show(out, figsize=figsize)

Microsolvated structure: 81 atoms.


Solid-liquid or solid-gas interfaces¶

First, create a slab using the ASE fcc111 function

rotation = "90x,0y,0z"  # sideview of slab
slab = fromASE(fcc111("Al", size=(4, 6, 3), vacuum=15.0, orthogonal=True, periodic=True))
show(slab, figsize=figsize, rotation=rotation)

print("water surrounding an Al slab, from an approximate density")
out = packmol_on_slab(slab, water, density=1.0)
printsummary(out)
out.write("al-water-pure.xyz")
show(out, figsize=figsize, rotation=rotation)

water surrounding an Al slab, from an approximate density
534 atoms, density = 1.325 g/cm^3, box = 11.455, 14.881, 34.677, formula = Al72H308O154

print("2-1 water-acetonitrile mixture surrounding an Al slab, from mole fractions and an approximate density")
out = packmol_on_slab(slab, [water, acetonitrile], mole_fractions=[x_water, x_acetonitrile], density=density)
printsummary(out)
out.write("al-water-acetonitrile.xyz")
show(out, figsize=figsize, rotation=rotation)

2-1 water-acetonitrile mixture surrounding an Al slab, from mole fractions and an approximate density
468 atoms, density = 1.260 g/cm^3, box = 11.455, 14.881, 34.677, formula = C66H231Al72N33O66


Pack inside voids in crystals¶

Use the packmol_in_void function. You can decrease tolerance if you need to pack very tightly. The default value for tolerance is 2.0.

bulk_Al = fromASE(bulk("Al", cubic=True).repeat((3, 3, 3)))
rotation = "90x,5y,5z"

out = packmol_in_void(
host=bulk_Al, molecules=[from_smiles("[H]"), from_smiles("[He]")], n_molecules=[50, 20], tolerance=1.5
)
printsummary(out)
out.write("al-bulk-with-h-he.xyz")

178 atoms, density = 2.819 g/cm^3, box = 12.150, 12.150, 12.150, formula = Al108H50He20


Bonds, atom properties (force field types, regions, …)¶

The packmol() function accepts the arguments keep_bonds and keep_atom_properties. These options will keep the bonds defined for the constitutent molecules, as well as any atomic properties.

The bonds and atom properties are easiest to see by printing the System block for an AMS job:

water = from_smiles("O")
n2 = from_smiles("N#N")

# delete properties coming from from_smiles
for at in water:
at.properties = Settings()
for at in n2:
at.properties = Settings()

water[1].properties.region = "oxygen_atom"
water[2].properties.mass = 2.014  # deuterium
water.delete_bond(water[1, 2])  # delete bond between atoms 1 and 2 (O and H)

out = packmol([water, n2], n_molecules=[2, 1], density=0.5)
print(AMSJob(molecule=out).get_input())

system
Atoms
O       3.0728760000       3.9143770000       1.9903040000 region=mol0,oxygen_atom
H       3.9160850000       3.5184940000       1.6850930000 mass=2.014 region=mol0
H       2.7876040000       4.6565520000       1.4140990000 region=mol0
O       4.9258210000       3.8909400000       3.9982150000 region=mol0,oxygen_atom
H       4.9810380000       3.6502800000       4.9468530000 mass=2.014 region=mol0
H       5.0008460000       4.8604790000       3.8619060000 region=mol0
N       1.1338120000       1.0294860000       0.9890770000 region=mol1
N       0.9243670000       1.6667980000       1.8734330000 region=mol1
End
BondOrders
1 3 1.0
4 6 1.0
7 8 3.0
End
Lattice
5.9692549746     0.0000000000     0.0000000000
0.0000000000     5.9692549746     0.0000000000
0.0000000000     0.0000000000     5.9692549746
End
End


By default, the packmol() function assigns regions called mol0, mol1, etc. to the different added molecules. The region_names option lets you set custom names.

out = packmol([water, n2], n_molecules=[2, 1], density=0.5, region_names=["water", "nitrogen_molecule"])
print(AMSJob(molecule=out).get_input())

system
Atoms
O       3.0728760000       3.9143770000       1.9903040000 region=oxygen_atom,water
H       3.9160850000       3.5184940000       1.6850930000 mass=2.014 region=water
H       2.7876040000       4.6565520000       1.4140990000 region=water
O       4.9258210000       3.8909400000       3.9982150000 region=oxygen_atom,water
H       4.9810380000       3.6502800000       4.9468530000 mass=2.014 region=water
H       5.0008460000       4.8604790000       3.8619060000 region=water
N       1.1338120000       1.0294860000       0.9890770000 region=nitrogen_molecule
N       0.9243670000       1.6667980000       1.8734330000 region=nitrogen_molecule
End
BondOrders
1 3 1.0
4 6 1.0
7 8 3.0
End
Lattice
5.9692549746     0.0000000000     0.0000000000
0.0000000000     5.9692549746     0.0000000000
0.0000000000     0.0000000000     5.9692549746
End
End


Below, we also set keep_atom_properties=False, this will remove the previous regions (in this example “oxygen_atom”) and mass.

out = packmol([water, n2], n_molecules=[2, 1], density=0.5, keep_atom_properties=False)
print(AMSJob(molecule=out).get_input())

system
Atoms
O       3.0728760000       3.9143770000       1.9903040000 region=mol0
H       3.9160850000       3.5184940000       1.6850930000 region=mol0
H       2.7876040000       4.6565520000       1.4140990000 region=mol0
O       4.9258210000       3.8909400000       3.9982150000 region=mol0
H       4.9810380000       3.6502800000       4.9468530000 region=mol0
H       5.0008460000       4.8604790000       3.8619060000 region=mol0
N       1.1338120000       1.0294860000       0.9890770000 region=mol1
N       0.9243670000       1.6667980000       1.8734330000 region=mol1
End
BondOrders
1 3 1.0
4 6 1.0
7 8 3.0
End
Lattice
5.9692549746     0.0000000000     0.0000000000
0.0000000000     5.9692549746     0.0000000000
0.0000000000     0.0000000000     5.9692549746
End
End


keep_bonds=False will additionally ignore any defined bonds:

out = packmol(
[water, n2],
n_molecules=[2, 1],
density=0.5,
region_names=["water", "nitrogen_molecule"],
keep_bonds=False,
keep_atom_properties=False,
)
print(AMSJob(molecule=out).get_input())

system
Atoms
O       3.0728760000       3.9143770000       1.9903040000 region=water
H       3.9160850000       3.5184940000       1.6850930000 region=water
H       2.7876040000       4.6565520000       1.4140990000 region=water
O       4.9258210000       3.8909400000       3.9982150000 region=water
H       4.9810380000       3.6502800000       4.9468530000 region=water
H       5.0008460000       4.8604790000       3.8619060000 region=water
N       1.1338120000       1.0294860000       0.9890770000 region=nitrogen_molecule
N       0.9243670000       1.6667980000       1.8734330000 region=nitrogen_molecule
End
Lattice
5.9692549746     0.0000000000     0.0000000000
0.0000000000     5.9692549746     0.0000000000
0.0000000000     0.0000000000     5.9692549746
End
End


Complete Python code¶

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

# ## Initial imports

from scm.plams import *
from ase.visualize.plot import plot_atoms
from ase.build import fcc111, bulk
import matplotlib.pyplot as plt

# ## Helper functions

def printsummary(mol, details=None):
if details:
density = details["density"]
else:
density = mol.get_density() * 1e-3
s = f"{len(mol)} atoms, density = {density:.3f} g/cm^3"
if mol.lattice:
s += f", box = {mol.lattice[0][0]:.3f}, {mol.lattice[1][1]:.3f}, {mol.lattice[2][2]:.3f}"
s += f", formula = {mol.get_formula()}"
if details:
s += f'\n#added molecules per species: {details["n_molecules"]}, mole fractions: {details["mole_fractions"]}'
print(s)

def show(mol, figsize=None, **kwargs):
"""Show a molecule in a Jupyter notebook"""
plt.figure(figsize=figsize or (2, 2))
plt.axis("off")
plot_atoms(toASE(mol), **kwargs)

# ## Liquid water (fluid with 1 component)
# First, create the gasphase molecule:

water = from_smiles("O")
show(water)

print("pure liquid from approximate number of atoms and exact density (in g/cm^3), cubic box with auto-determined size")
out = packmol(water, n_atoms=194, density=1.0)
printsummary(out)
out.write("water-1.xyz")
show(out)

print("pure liquid from approximate density (in g/cm^3) and an orthorhombic box")
out = packmol(water, density=1.0, box_bounds=[0.0, 0.0, 0.0, 8.0, 12.0, 14.0])
printsummary(out)
out.write("water-2.xyz")
show(out)

print("pure liquid with explicit number of molecules and exact density")
out = packmol(water, n_molecules=64, density=1.0)
printsummary(out)
out.write("water-3.xyz")
show(out)

print("pure liquid with explicit number of molecules and box")
out = packmol(water, n_molecules=64, box_bounds=[0.0, 0.0, 0.0, 12.0, 13.0, 14.0])
printsummary(out)
out.write("water-4.xyz")
show(out)

print("water-5.xyz: pure liquid in non-orthorhombic box (requires AMS2022 or later)")
# first place the molecules in a cuboid surrounding the desired lattice
# then gradually change into the desired lattice using refine_lattice()
# note that the molecules may become distorted by this procedure
lattice = [[10.0, 2.0, -1.0], [-5.0, 8.0, 0.0], [0.0, -2.0, 11.0]]
temp_out = packmol(
water,
n_molecules=32,
box_bounds=[
0,
0,
0,
max(lattice[i][0] for i in range(3)) - min(lattice[i][0] for i in range(3)),
max(lattice[i][1] for i in range(3)) - min(lattice[i][1] for i in range(3)),
max(lattice[i][2] for i in range(3)) - min(lattice[i][2] for i in range(3)),
],
)
out = refine_lattice(temp_out, lattice=lattice)
if out is not None:
out.write("water-5.xyz")
print(
"Top: system in surrounding orthorhombic box before calling refine_lattice(). Bottom: System in non-orthorhombic box after calling refine_lattice()"
)
show(temp_out)
show(out)

# ## Water-acetonitrile mixture (fluid with 2 or more components)
# Let's also create a single acetonitrile molecule:

acetonitrile = from_smiles("CC#N")
show(acetonitrile)

# Set the desired mole fractions and density. Here, the density is calculated as the weighted average of water (1.0 g/cm^3) and acetonitrile (0.76 g/cm^3) densities, but you could use any other density.

# MIXTURES
x_water = 0.666  # mole fraction
x_acetonitrile = 1 - x_water  # mole fraction
density = (x_water * 1.0 + x_acetonitrile * 0.76) / (
x_water + x_acetonitrile
)  # weighted average of pure component densities

print("MIXTURES")
print(f"x_water = {x_water:.3f}")
print(f"x_acetonitrile = {x_acetonitrile:.3f}")
print(f"target density = {density:.3f} g/cm^3")

# By setting return_details=True, you can get information about the mole fractions of the returned system. They may not exactly match the mole fractions you put in.

print(
"2-1 water-acetonitrile from approximate number of atoms and exact density (in g/cm^3), cubic box with auto-determined size"
)
out, details = packmol(
molecules=[water, acetonitrile],
mole_fractions=[x_water, x_acetonitrile],
n_atoms=200,
density=density,
return_details=True,
)
printsummary(out, details)
out.write("water-acetonitrile-1.xyz")
show(out)

# The details is a dictionary as follows:

for k, v in details.items():
print(f"{k}: {v}")

print("2-1 water-acetonitrile from approximate density (in g/cm^3) and box bounds")
out, details = packmol(
molecules=[water, acetonitrile],
mole_fractions=[x_water, x_acetonitrile],
box_bounds=[0, 0, 0, 13.2, 13.2, 13.2],
density=density,
return_details=True,
)
printsummary(out, details)
out.write("water-acetonitrile-2.xyz")
show(out)

print("2-1 water-acetonitrile from explicit number of molecules and density, cubic box with auto-determined size")
out, details = packmol(molecules=[water, acetonitrile], n_molecules=[32, 16], density=density, return_details=True)
printsummary(out, details)
out.write("water-acetonitrile-3.xyz")
show(out)

print("2-1 water-acetonitrile from explicit number of molecules and box")
out = packmol(molecules=[water, acetonitrile], n_molecules=[32, 16], box_bounds=[0, 0, 0, 13.2, 13.2, 13.2])
printsummary(out)
out.write("water-acetonitrile-4.xyz")
show(out)

# ## Pack inside sphere
#
# Set sphere=True to pack in a sphere (non-periodic) instead of in a periodic box. The sphere will be centered near the origin.

print("water in a sphere from exact density and number of molecules")
out, details = packmol(molecules=[water], n_molecules=[100], density=1.0, return_details=True, sphere=True)
printsummary(out, details)
print(f"Center of mass xyz (ang): {out.get_center_of_mass()}")
out.write("water-sphere.xyz")
show(out)

print(
"2-1 water-acetonitrile in a sphere from exact density (in g/cm^3) and approximate number of atoms and mole fractions"
)
out, details = packmol(
molecules=[water, acetonitrile],
mole_fractions=[x_water, x_acetonitrile],
n_atoms=500,
density=density,
return_details=True,
sphere=True,
)
printsummary(out, details)
out.write("water-acetonitrile-sphere.xyz")
show(out)

# ## Packing ions, total system charge
#
# The total system charge will be sum of the charges of the constituent molecules.
#
# In PLAMS, molecule.properties.charge specifies the charge:

ammonium = from_smiles("[NH4+]")  # ammonia.properties.charge == +1
chloride = from_smiles("[Cl-]")  # chloride.properties.charge == -1
print("3 water molecules, 3 ammonium, 1 chloride (non-periodic)")
print("Initial charges:")
print(f"Water: {water.properties.get('charge', 0)}")
print(f"Ammonia: {ammonium.properties.get('charge', 0)}")
print(f"Chloride: {chloride.properties.get('charge', 0)}")
out = packmol(molecules=[water, ammonium, chloride], n_molecules=[3, 3, 1], density=0.4, sphere=True)
tot_charge = out.properties.get("charge", 0)
print(f"Total charge of packmol-generated system: {tot_charge}")
out.write("water-ammonium-chloride.xyz")
show(out)

# ## Microsolvation
# packmol_microsolvation can create a microsolvation sphere around a solute.

out = packmol_microsolvation(solute=acetonitrile, solvent=water, density=1.5, threshold=4.0)
# for microsolvation it's a good idea to have a higher density than normal to get enough solvent molecules
print(f"Microsolvated structure: {len(out)} atoms.")
out.write("acetonitrile-microsolvated.xyz")

figsize = (3, 3)
show(out, figsize=figsize)

# ## Solid-liquid or solid-gas interfaces
# First, create a slab using the ASE fcc111 function

rotation = "90x,0y,0z"  # sideview of slab
slab = fromASE(fcc111("Al", size=(4, 6, 3), vacuum=15.0, orthogonal=True, periodic=True))
show(slab, figsize=figsize, rotation=rotation)

print("water surrounding an Al slab, from an approximate density")
out = packmol_on_slab(slab, water, density=1.0)
printsummary(out)
out.write("al-water-pure.xyz")
show(out, figsize=figsize, rotation=rotation)

print("2-1 water-acetonitrile mixture surrounding an Al slab, from mole fractions and an approximate density")
out = packmol_on_slab(slab, [water, acetonitrile], mole_fractions=[x_water, x_acetonitrile], density=density)
printsummary(out)
out.write("al-water-acetonitrile.xyz")
show(out, figsize=figsize, rotation=rotation)

# ## Pack inside voids in crystals
#
# Use the packmol_in_void function. You can decrease tolerance if you need to pack very tightly. The default value for tolerance is 2.0.

bulk_Al = fromASE(bulk("Al", cubic=True).repeat((3, 3, 3)))
rotation = "90x,5y,5z"

out = packmol_in_void(
host=bulk_Al, molecules=[from_smiles("[H]"), from_smiles("[He]")], n_molecules=[50, 20], tolerance=1.5
)
printsummary(out)
out.write("al-bulk-with-h-he.xyz")

# ## Bonds, atom properties (force field types, regions, ...)
#
# The packmol() function accepts the arguments keep_bonds and keep_atom_properties. These options will keep the bonds defined for the constitutent molecules, as well as any atomic properties.
#
# The bonds and atom properties are easiest to see by printing the System block for an AMS job:

water = from_smiles("O")
n2 = from_smiles("N#N")

# delete properties coming from from_smiles
for at in water:
at.properties = Settings()
for at in n2:
at.properties = Settings()

water[1].properties.region = "oxygen_atom"
water[2].properties.mass = 2.014  # deuterium
water.delete_bond(water[1, 2])  # delete bond between atoms 1 and 2 (O and H)

out = packmol([water, n2], n_molecules=[2, 1], density=0.5)
print(AMSJob(molecule=out).get_input())

# By default, the packmol() function assigns regions called mol0, mol1, etc. to the different added molecules. The region_names option lets you set custom names.

out = packmol([water, n2], n_molecules=[2, 1], density=0.5, region_names=["water", "nitrogen_molecule"])
print(AMSJob(molecule=out).get_input())

# Below, we also set keep_atom_properties=False, this will remove the previous regions (in this example "oxygen_atom") and mass.

out = packmol([water, n2], n_molecules=[2, 1], density=0.5, keep_atom_properties=False)
print(AMSJob(molecule=out).get_input())

# keep_bonds=False will additionally ignore any defined bonds:

out = packmol(
[water, n2],
n_molecules=[2, 1],
density=0.5,
region_names=["water", "nitrogen_molecule"],
keep_bonds=False,
keep_atom_properties=False,
)
print(AMSJob(molecule=out).get_input())