# Butanol QM/MM optimisation with electrostatic embedding
# ORCA/GULP VERSION

# Initialise ChemShell
from chemsh import *

# Read in the butanol geometry from disk
my_mol = Fragment(coords="butanol.cjson")

# Define GULP MM forcefield for butanol
# Here we are following AMBER conventions (e.g. vdw and coulomb scaling)
ff = '''#
# molmec removes Coulomb terms between bonded atoms
# and between atoms bonded to a common third atom
keyword molmec
# atom types
species
C_T core
H_C core
O_H core
H_O core
# Harmonic bond terms
harmonic bond kcal
C_T C_T  620.0   1.526
C_T H_C  680.0   1.090
C_T O_H  640.0   1.410
H_O O_H 1106.0   0.960
# Angle terms
three bond kcal
C_T C_T C_T    80.0  109.50
C_T C_T H_C   100.0  109.50
C_T C_T O_H   100.0  109.50
C_T H_C H_C    70.0  109.50
O_H C_T H_O   110.0  108.50
C_T H_C O_H    70.0  109.50
O_H H_O H_O    94.0  104.50
# Torsional terms
torsion bond kcal
X C_T C_T X   0.1555555  3  0.0
X C_T O_H X   0.1666666  3  0.0
# VDW terms
# scale vdw by 0.5
# last number in entries is a cutoff
lennard 12 6 kcal x13 0.5
 H_O   H_O  0.0  0.0   99.9
 H_O   H_C  0.0  0.0   99.9
 H_O   O_H  0.0  0.0   99.9
 H_O   C_T  0.0  0.0   99.9
 H_C   H_C     7516.0     21.73    99.9
 H_C   O_H    68280.0    125.3     99.9
 H_C   C_T    97170.0    126.9     99.9
 O_H   O_H   581800.0    699.7     99.9
 O_H   C_T   791500.0    693.1     99.9
 C_T   C_T  1043000.0    675.6     99.9
# scale 1-4 Coulomb by 0.5 by subtraction
coulomb o14
 C_T C_T 0.5 0.0 99.9
 C_T H_C 0.5 0.0 99.9
 C_T O_H 0.5 0.0 99.9
 C_T H_O 0.5 0.0 99.9
 H_C H_C 0.5 0.0 99.9
 H_C O_H 0.5 0.0 99.9
 H_C H_O 0.5 0.0 99.9
 O_H O_H 0.5 0.0 99.9
 O_H H_O 0.5 0.0 99.9
 H_O H_O 0.5 0.0 99.9
'''

# Define QM and MM parts of the calculation
my_qm = ORCA(method="dft", functional="blyp", basis="3-21g", convergence="TightSCF")
my_mm = GULP(ff=ff, molecule=True)

# Define the combined QM/MM calculation
my_qmmm = QMMM(frag=my_mol, qm_region=range(5), 
            qm=my_qm, mm=my_mm, embedding="electrostatic", coupling="covalent")

# Request an energy minimisation
my_opt = Opt(theory=my_qmmm,
             algorithm="lbfgs", tolerance=0.0003,
             maxcycle=300, maxene=400, trust_radius="energy")
my_opt.run()
