"""This module defines an ASE interface to MOPAC.
Set $ASE_MOPAC_COMMAND to something like::
LD_LIBRARY_PATH=/path/to/lib/ \
MOPAC_LICENSE=/path/to/license \
/path/to/MOPAC2012.exe PREFIX.mop 2> /dev/null
"""
import os
import numpy as np
from ase import Atoms
from ase.calculators.calculator import FileIOCalculator, ReadError, Parameters
from ase.units import kcal, mol, Debye
[docs]class MOPAC(FileIOCalculator):
implemented_properties = ['energy', 'forces', 'dipole', 'magmom']
command = 'mopac PREFIX.mop 2> /dev/null'
discard_results_on_any_change = True
default_parameters = dict(
method='PM7',
task='1SCF GRADIENTS',
charge=0,
relscf=0.0001)
methods = ['AM1', 'MNDO', 'MNDOD', 'PM3', 'PM6', 'PM6-D3', 'PM6-DH+',
'PM6-DH2', 'PM6-DH2X', 'PM6-D3H4', 'PM6-D3H4X', 'PMEP', 'PM7',
'PM7-TS', 'RM1']
def __init__(self, restart=None,
ignore_bad_restart_file=FileIOCalculator._deprecated,
label='mopac', atoms=None, **kwargs):
"""Construct MOPAC-calculator object.
Parameters:
label: str
Prefix for filenames (label.mop, label.out, ...)
Examples:
Use default values to do a single SCF calculation and print
the forces (task='1SCF GRADIENTS'):
>>> from ase.build import molecule
>>> from ase.calculators.mopac import MOPAC
>>> atoms = molecule('O2')
>>> atoms.calc = MOPAC(label='O2')
>>> atoms.get_potential_energy()
>>> eigs = atoms.calc.get_eigenvalues()
>>> somos = atoms.calc.get_somo_levels()
>>> homo, lumo = atoms.calc.get_homo_lumo_levels()
Use the internal geometry optimization of Mopac:
>>> atoms = molecule('H2')
>>> atoms.calc = MOPAC(label='H2', task='GRADIENTS')
>>> atoms.get_potential_energy()
Read in and start from output file:
>>> atoms = MOPAC.read_atoms('H2')
>>> atoms.calc.get_homo_lumo_levels()
"""
FileIOCalculator.__init__(self, restart, ignore_bad_restart_file,
label, atoms, **kwargs)
def write_input(self, atoms, properties=None, system_changes=None):
FileIOCalculator.write_input(self, atoms, properties, system_changes)
p = self.parameters
# Build string to hold .mop input file:
s = p.method + ' ' + p.task + ' '
if p.relscf:
s += 'RELSCF={0} '.format(p.relscf)
# Write charge:
if p.charge:
charge = p.charge
else:
charge = atoms.get_initial_charges().sum()
if charge != 0:
s += 'CHARGE={0} '.format(int(round(charge)))
magmom = int(round(abs(atoms.get_initial_magnetic_moments().sum())))
if magmom:
s += (['DOUBLET', 'TRIPLET', 'QUARTET', 'QUINTET'][magmom - 1] +
' UHF ')
s += '\nTitle: ASE calculation\n\n'
# Write coordinates:
for xyz, symbol in zip(atoms.positions, atoms.get_chemical_symbols()):
s += ' {0:2} {1} 1 {2} 1 {3} 1\n'.format(symbol, *xyz)
for v, p in zip(atoms.cell, atoms.pbc):
if p:
s += 'Tv {0} {1} {2}\n'.format(*v)
with open(self.label + '.mop', 'w') as fd:
fd.write(s)
def get_spin_polarized(self):
return self.nspins == 2
def get_index(self, lines, pattern):
for i, line in enumerate(lines):
if line.find(pattern) != -1:
return i
def read(self, label):
FileIOCalculator.read(self, label)
if not os.path.isfile(self.label + '.out'):
raise ReadError
with open(self.label + '.out') as fd:
lines = fd.readlines()
self.parameters = Parameters(task='', method='')
p = self.parameters
parm_line = self.read_parameters_from_file(lines)
for keyword in parm_line.split():
if 'RELSCF' in keyword:
p.relscf = float(keyword.split('=')[-1])
elif keyword in self.methods:
p.method = keyword
else:
p.task += keyword + ' '
p.task.rstrip()
self.atoms = self.read_atoms_from_file(lines)
self.read_results()
def read_atoms_from_file(self, lines):
"""Read the Atoms from the output file stored as list of str in lines.
Parameters:
lines: list of str
"""
# first try to read from final point (last image)
i = self.get_index(lines, 'FINAL POINT AND DERIVATIVES')
if i is None: # XXX should we read it from the input file?
assert 0, 'Not implemented'
lines1 = lines[i:]
i = self.get_index(lines1, 'CARTESIAN COORDINATES')
j = i + 2
symbols = []
positions = []
while not lines1[j].isspace(): # continue until we hit a blank line
l = lines1[j].split()
symbols.append(l[1])
positions.append([float(c) for c in l[2: 2 + 3]])
j += 1
return Atoms(symbols=symbols, positions=positions)
def read_parameters_from_file(self, lines):
"""Find and return the line that defines a Mopac calculation
Parameters:
lines: list of str
"""
for i, line in enumerate(lines):
if line.find('CALCULATION DONE:') != -1:
break
lines1 = lines[i:]
for i, line in enumerate(lines1):
if line.find('****') != -1:
return lines1[i + 1]
def read_results(self):
"""Read the results, such as energy, forces, eigenvalues, etc.
"""
FileIOCalculator.read(self, self.label)
if not os.path.isfile(self.label + '.out'):
raise ReadError
with open(self.label + '.out') as fd:
lines = fd.readlines()
for i, line in enumerate(lines):
if line.find('TOTAL ENERGY') != -1:
self.results['energy'] = float(line.split()[3])
elif line.find('FINAL HEAT OF FORMATION') != -1:
self.final_hof = float(line.split()[5]) * kcal / mol
elif line.find('NO. OF FILLED LEVELS') != -1:
self.nspins = 1
self.no_occ_levels = int(line.split()[-1])
elif line.find('NO. OF ALPHA ELECTRON') != -1:
self.nspins = 2
self.no_alpha_electrons = int(line.split()[-1])
self.no_beta_electrons = int(lines[i+1].split()[-1])
self.results['magmom'] = abs(self.no_alpha_electrons -
self.no_beta_electrons)
elif line.find('FINAL POINT AND DERIVATIVES') != -1:
forces = [-float(line.split()[6])
for line in lines[i + 3:i + 3 + 3 * len(self.atoms)]]
self.results['forces'] = np.array(
forces).reshape((-1, 3)) * kcal / mol
elif line.find('EIGENVALUES') != -1:
if line.find('ALPHA') != -1:
j = i + 1
eigs_alpha = []
while not lines[j].isspace():
eigs_alpha += [float(eps) for eps in lines[j].split()]
j += 1
elif line.find('BETA') != -1:
j = i + 1
eigs_beta = []
while not lines[j].isspace():
eigs_beta += [float(eps) for eps in lines[j].split()]
j += 1
eigs = np.array([eigs_alpha, eigs_beta]).reshape(2, 1, -1)
self.eigenvalues = eigs
else:
eigs = []
j = i + 1
while not lines[j].isspace():
eigs += [float(e) for e in lines[j].split()]
j += 1
self.eigenvalues = np.array(eigs).reshape(1, 1, -1)
elif line.find('DIPOLE ') != -1:
self.results['dipole'] = np.array(
lines[i + 3].split()[1:1 + 3], float) * Debye
def get_eigenvalues(self, kpt=0, spin=0):
return self.eigenvalues[spin, kpt]
def get_homo_lumo_levels(self):
eigs = self.eigenvalues
if self.nspins == 1:
nocc = self.no_occ_levels
return np.array([eigs[0, 0, nocc - 1], eigs[0, 0, nocc]])
else:
na = self.no_alpha_electrons
nb = self.no_beta_electrons
if na == 0:
return None, self.eigenvalues[1, 0, nb - 1]
elif nb == 0:
return self.eigenvalues[0, 0, na - 1], None
else:
eah, eal = eigs[0, 0, na - 1: na + 1]
ebh, ebl = eigs[1, 0, nb - 1: nb + 1]
return np.array([max(eah, ebh), min(eal, ebl)])
def get_somo_levels(self):
assert self.nspins == 2
na, nb = self.no_alpha_electrons, self.no_beta_electrons
if na == 0:
return None, self.eigenvalues[1, 0, nb - 1]
elif nb == 0:
return self.eigenvalues[0, 0, na - 1], None
else:
return np.array([self.eigenvalues[0, 0, na - 1],
self.eigenvalues[1, 0, nb - 1]])
def get_final_heat_of_formation(self):
"""Final heat of formation as reported in the Mopac output file
"""
return self.final_hof