# -*- coding: utf-8 -*-
"""Contains a class for QC-calculations using Molpro."""
from __future__ import print_function
import os
import numpy as np
import xml.etree.ElementTree as et
import pickle
from .iface import iface
from .. import constants as c
[docs]class iface_molpro(iface):
"""
Interface to Molpro for QC-calculations, inherits from `metafalcon.interfaces.iface.iface`.
Parameters
----------
symbols : list of str
atomic symbols of the molecule
coords : np.2darray (shape: (N, 3))
atomic xyz-coordinates of the molecule (bohr)
kwargs :
See below
Keyword Arguments
-----------------
nproc : int
number of processors to be used in parallel calculation
mem : int
amount of memory to be used by molpro (in MB)
Attributes
----------
symbols : list of str
atomic symbols of the molecule
coords : np.2darray (shape: (N, 3))
atomic xyz-coordinates of the molecule (bohr)
energy : float
electronic energy of the molecule calculated most recently (hartrees)
gradient : np.2darray (shape: (N, 3))
electronic energy gradient of the molecule with current coordinates (hartree/bohr)
nproc : int
number of processors to be used in parallel calculation
head : list of str
lines to be used as head for the molpro inputfile
tail : list of str
lines to be used as tail for the molpro inputfile
See Also
--------
iface.set_symbols()
iface.set_coords()
iface.get_energy()
iface.get_gradient()
"""
def __init__(self, symbols, coords, **kwargs):
"""Construct Molpro interface object."""
iface.__init__(self, symbols, coords)
# self.saveorbitals = prm["saveorbitals"]
self.nproc = kwargs["nproc"]
self.mem = kwargs["mem"]
# save temporary wavefunction-file
if os.path.exists("original.wfu"):
os.system("cp original.wfu temp.wfu")
# if self.saveorbitals and not os.path.exists("wfu"):
# os.mkdir("wfu")
# read in head and tail
with open("head") as f, open("tail") as g:
self.head = f.readlines()
self.tail = g.readlines()
[docs] def calculate(self):
"""Calculate energies and gradients for the current set of coordinates."""
def has_error(tree):
return bool(list(tree.iter(ns + "error")))
def get_error():
with open("molpro.xml") as f:
for line in f:
if "ERROR DETECTED" in line:
[f.next() for i in range(2)]
return f.next().strip()
os.system("rm -f molpro.*")
with open("molpro.in", "w") as f:
f.writelines(self.head)
for sym, coord in zip(self.symbols, self.coords*c.a0):
f.write("%s %20.12f %20.12f %20.12f \n" % (sym, coord[0], coord[1], coord[2]))
f.writelines(self.tail)
# if not self.step % 10:
# f.write("\n put,molden,molden_%i.input\n" % self.step)
cwd = os.getcwd()
os.system("molpro -n %d -m %dm -d %s -I %s -W %s molpro.in" % (self.nproc, self.mem, cwd,
cwd, cwd))
ns = "{http://www.molpro.net/schema/molpro-output}"
tree = et.parse("molpro.xml")
if has_error(tree):
error = get_error()
print("found error:\n", error)
if error == "EXCESSIVE GRADIENT IN CI":
pspace = 1
while has_error(tree) and pspace < 100:
print("setting pspace thresh to ", pspace)
os.rename("molpro.in", "molpro.in.bck")
with open("molpro.in.bck") as f, open("molpro.in", "w") as g:
for line in f:
if "pspace" in line:
continue
g.write(line)
if "occ" in line:
g.write("pspace,%d\n" % pspace)
os.system("molpro -n %d -m %dm -d %s -I %s -W %s molpro.in" % (self.nproc,
self.mem, cwd,
cwd, cwd))
tree = et.parse("molpro.xml")
pspace *= 2
else:
if os.path.exists("original.wfu"):
os.system("cp original.wfu temp.wfu")
print("using original orbitals to restart calculation")
print("restart molpro calculation")
os.system("molpro -n %d -m %dm -d %s -I %s -W %s molpro.in" % (self.nproc,
self.mem, cwd, cwd,
cwd))
tree = et.parse("molpro.xml")
if not has_error(tree):
print("problem solved")
else:
print("problem could not be solved")
# # save orbitals every few steps
# if self.saveorbitals and not self.step % self.saveorbitals:
# os.system("cp temp.wfu wfu/%i.wfu" % self.step)
# end save orbitals
energies = [float(prop.attrib["value"]) for prop in
tree.iter(ns + "property") if prop.attrib["name"] == "Energy" and
prop.attrib["method"] == "MCSCF"]
with open("energies.pickle", "w") as f:
pickle.dump(energies, f, protocol=pickle.HIGHEST_PROTOCOL)
gradients = [np.array(gradient.text.split(), dtype=float).reshape(self.coords.shape) for
gradient in tree.iter(ns + "gradient")]
with open("gradients.pickle", "w") as f:
pickle.dump(gradients, f, protocol=pickle.HIGHEST_PROTOCOL)
try:
self.energy = energies[0]
self.gradient = gradients[0]
self.converged = True
except KeyError as e:
self.energy = 0.0
self.gradient = np.zeros_like(self.coords)
self.converged = False
print(e)