from ase import Atoms
import ase
from ase.build import add_adsorbate
from ase.visualize import view
from ase.constraints import FixAtoms
from ase.optimize import LBFGS, QuasiNewton, FIRE
from ase.io import write,read
from gpaw import GPAW, PW
from ase.calculators.dftd3 import DFTD3
import numpy as np
from ase.neb import NEB
import matplotlib.pyplot as plt
from ase.optimize.sciopt import SciPyFminCG

slab = read('Ni.cif') # from the article
view(slab )
## Find the initial and final states for the reaction.
Mg2S2 = Atoms('Mg2S2',[[0,0,0],[3.7,0,0],[1.85,1.85,0],[1.85,-1.85,0]]) # Mg2S2 from Materials Project 
## slab.set_pbc((1, 1, 1)) # Set periodic boundary conditions. (1,1,0) means periodic in x and y, but not in z.
add_adsorbate(slab,Mg2S2,height=4.0,position=(2.383, 5.455)) # height means the distance between the adsorbate and the surface. position means the x-y position of the adsorbate, 'fcc' means the position of the (0,0,0) with shift to fcc postion
## slab.center(vacuum=10, axis=2) # set a 20 A vacuum along the z-axis
view(slab)

slab.calc = GPAW(xc='PBE', mode=PW(340))
relax = QuasiNewton(slab) # Structure optimization using the QuasiNewton method (BFGSLineSearch)
relax.run()

view(slab)
write('NiRelaxed.cif', slab)
write('NiRelaxed.traj', slab)