Dear Are users,
I’m using the excellent ASE program to combine different calculator using the QM/MM subtractive scheme. In particular, I’m using for the QM part the ORCA calculator and for the MM part the ML potential M3GNET. The optimization are running quite fast and good, I report here the code that I’m using if anyone is interested:
import os
import ase
import m3gnet
from m3gnet.models import *
from ase.optimize.fire import FIRE
from ase.io.trajectory import Trajectory
from ase.calculators.orca import ORCA
from ase.atoms import Atoms
from ase.calculators.qmmm import SimpleQMMM
from ase.constraints import FixAtoms, FixBondLengths
from ase.calculators.mixing import SumCalculator
from dftd3.ase import DFTD3
geom=ase.io.read('HF3c_CO_unfix_grain_fix.xyz')
#take the name of the .py file without extension
output_filename = os.path.splitext(os.path.basename(__file__))[0]
#energy = atoms.get_potential_energy()
#forces = atoms.get_forces()
qmindex = [ 9, 10, 11, 138, 139, 140, 198, 199, 200, 219, 220, 221, 261, 262, 263, 399, 400, 401, 408, 409, 410, 435, 436, 437, 447, 448, 449, 465, 466, 467, 555, 556, 557, 588, 589, 590, 600, 601]
mmindex=[]
for element in geom:
if element.index not in qmindex:
mmindex.append(element.index)
geom.constraints=FixAtoms(mmindex)
#defining the calculator for the QM and the QM1-MM part, ase allows to differentiate between the method used for the low level real system and the model, in this example are set identical
QM_calc = ORCA(label= output_filename + '_ORCA', orcasimpleinput='B97-3c', orcablocks='%pal nprocs 48 end \n %maxcore 1000')
MM1_calc = SumCalculator([DFTD3(method="PBE", damping="d3bj"), M3GNetCalculator(potential=Potential(M3GNet.load()), stress_weight=1.0)])
MM2_calc = MM1_calc
#setting the ase subtractive schem calculator
geom.calc= SimpleQMMM(qmindex, QM_calc, MM1_calc, MM2_calc, vacuum=None)
opt = FIRE(geom, trajectory= output_filename + '.traj')
opt.run(fmax=0.01)
traj = Trajectory(output_filename + '.traj')
nsteps = len(traj)
for i in range(0, nsteps):
atoms= traj[i]
ase.io.write(output_filename + '_traj.xyz', atoms, format="xyz", append=True)
At this point however, the vibration calculation is demanding too much time with respect to the vibration directly computed in the ORCA program. This is of course due to the fact that ORCA divide each displacement for each processor of the node. What I was wondering was, is there a way to do something similar with ase?
Thank you in advance