# Vibrational analysis of adsorbed CO2

# Initialise ChemShell
from chemsh import *
from chemsh.utils import testutils
from numpy import append as np_append

# Import optimised MgO/CO2 geometry
mgo = Fragment(coords='mgoco2end.pun', connect_mode=None)

# QM region contains MgO cluster region 1 and CO2
co2_atoms = [2370, 2371, 2372]
qm_region = np_append(mgo.getRegion(1), co2_atoms)
print("QM region is: ", qm_region)

# Get QM charge from input fragment
qm_charge = round(mgo.totalcharge)
print("QM charge is: ", qm_charge)

# Only CO2 is active for vibrational analysis
print("Active region is: ", co2_atoms)

# QM/MM settings
qm_theory = NWChem(method='dft',
                   functional='b3lyp',
                   basis='mgoco2_nwchem.basis',
                   ecp='mgo_nwchem.ecp',
                   charge=qm_charge,
                   path="mpirun -np 4 nwchem")

mm_theory = GULP(ff='mgo_2body.ff',
                 conjugate=True)

qmmm = QMMM(qm=qm_theory,
            mm=mm_theory,
            frag=mgo,
            qm_region=qm_region,
            coupling='ionic',
            shl_maxcycles=20)

# Run vibrational analysis
opt = Opt(theory=qmmm, thermal=True, active=co2_atoms)
opt.run()
