How to model adsorption using Quacc.
Placing the Adsorbate
imports:
from ase import Atoms
from ase.io import read, write
import numpy as np
Aquire a structure object of the thing you're adsorbing to (ie .cif/POSCAR/CONTCAR/.xyz )
load in the structure object using ase (syntax may be slightly different depending on where you're getting your structure from).
structure = read("CONTCAR")
Identify where you want to adsorb your molecule (example: H2)
This can be done by finding the coordinates of the nearest atom to where you want your molecule to be in VESTA (load in structure object of your original file and double click on atom to get coordinates)
alternatively, from ase.visualize import view and then view(structure) has the same functionality (although this must be done on your local machine)
click here for more information on ase gui on the clusters Then, alter the coordinates slightly (may involve some trial and error) until your molecule is in the correct location.
declare the molecule center using those coordinates you found
h2_center = np.array([3.0, 5.0, 4.0])
H2 is a diatomic molecule (w/ a bond length of .74 A) so we set each individual H atom's location around the H2_center
in this case the y axis (aka b in VESTA) is where we want the H2 molecule to extend along, but this can change depending on what orientation suits your system.
you can do an analogous procedure if you have a different molecule
bond_length = 0.74
h2_positions = [
h2_center - np.array([0, bond_length/2, 0]),
h2_center + np.array([0, bond_length/2, 0])
]
Now, place the molecule and combine the two systems
h2_atoms = Atoms('H2', positions=h2_positions)
combined_system = mewdid + h2_atoms
Finally, write out a POSCAR which will serve as an input to your DFT calculation
write('POSCAR', combined_system, format='vasp')
Setting up a Relaxation
- go about this in the usual fashion except set
isif = 2 to avoid cell volume changes
here's a sample:
from ase.io import read
from quacc.calculators.vasp import Vasp
atoms = read("POSCAR")
calc = Vasp(atoms, preset="RosenSetPBE", isif=2, nsw=200)
atoms.calc = calc
atoms.get_potential_energy()