# -*- coding: utf-8 -*-
"""Contains a class for using the energy gap between two electronic states as a CV."""
from __future__ import print_function
import numpy as np
import os
import multiprocessing as mp
import xml.etree.ElementTree as et
from .cv import cv
from .. import constants as c
from .. import atoms as at
[docs]class cv_noongap(cv):
"""
NOON gap as collective variable, inherits from `metafalcon.cvs.cv.cv`.
Parameters
----------
idx : int
index of collective variable in the input file
symbols : list of str
atomic symbols of the molecule
**kwargs :
See below
Keyword Arguments
-----------------
no_index1 : int
index of first NOON (close to 0)
no_index2 : int
index of second NOON (close to 2)
nstates : int
number of electronic states -(default: 1)
state : int
electronic state index (default: 0)
nproc : int
number of parallel processes for the calculation of numerical derivatives (default: 1)
do_multistate : bool
whether to use this CV for multistate metadynamics (default: False)
"""
def __init__(self, idx, symbols, **kwargs):
"""Construct energy gap CV object."""
cv.__init__(self, idx, symbols, **kwargs)
self.x = np.linspace(0., 2.0, 200)
self.ticks = 0.2
self.no_index1 = self.get_kwargs("no_index1") # lower NOON
self.no_index2 = self.get_kwargs("no_index2") # upper NOON
self.nstates = self.get_kwargs("nstates", default=1)
self.state = self.get_kwargs("state", default=0)
self.nproc = self.get_kwargs("nproc", default=1)
self.iface = self.get_kwargs("iface", default="molpro")
self.h = self.get_kwargs("h", default=1e-2)
if self.iface == "molpro":
self.norbitals = get_norbitals_molpro()
self.do_multistate = self.get_kwargs("do_multistate", default=False)
[docs] def initialize(self):
cv.initialize(self)
if self.do_multistate:
from ..meta import metadynamics
self.offdiag = metadynamics(symbols=self.symbols, file_extension="multistate",
restart=self.restart)
self.offdiag.read_input()
self.offdiag.initialize_vg() # overwrites exisiting files even during reconstruct!!
if not (os.path.exists("deltan_meta.dat") and self.restart):
with open("deltan_org.dat", "w") as f:
f.write("# NOON gap\n")
with open("deltan_meta.dat", "w") as f:
f.write("# NOON gap\n")
with open("noons_all.dat", "w") as f:
f.write("# Natural Orbital Occupation Numbers (NOONs) \n")
[docs] def get_noon(self, index, qcout=None):
if qcout == None and self.iface == "molpro":
qcout = "molpro.out"
all_ev = getNatPop(self.nstates, self.norbitals, qcout, iface=self.iface)
if qcout == "job.out" or qcout == "molpro.out":
time = "%8.2f " % (self.step*self.tstep)
noons = " ".join(["%16.9e" % noon for noon in all_ev[self.state][0]])
with open("noons_all.dat", "a") as f:
f.write(time + noons + "\n")
return all_ev[self.state][0][self.norbitals-index-1]
[docs] def get_dnoon_numerical(self, index, s_current, coords, h=1e-2):
h_exp = int(-np.log10(h))
if not [dirname for dirname in os.listdir(".") if
dirname.startswith("molpro_%i_" % self.step)]:
os.system("rm -rf molpro_*_*" % h_exp)
for i in range(coords.size):
if self.iface == "molpro":
os.mkdir("molpro_%i_%i" % (self.step, i))
os.system("cp head_noon tail_noon temp.wfu molpro_%i_%i" % (self.step, i))
pool = mp.Pool(processes=self.nproc)
coords_h = []
for i in range(coords.size):
coords_h.append(np.copy(coords))
coords_h[i][i/3, i%3] += h
if self.iface == "molpro":
step_i_sym_coords_h = [(self.step, i, self.symbols, coord) for i, coord in
enumerate(coords_h)]
pool.map(run_molpro, step_i_sym_coords_h)
pool.close()
pool.join()
if self.iface == "molpro":
s_h = np.array([self.get_noon(index, qcout="molpro_%i_%i/molpro.out" % (self.step, i))
for i in range(coords.size)]).reshape(coords.shape)
return (s_h-s_current)/h
[docs] def set_s_and_ds(self, coords):
"""
Calculate NOON gap and its gradient for the current set of coordinates.
Parameters
----------
coords : np.2darray
cartesian coordinates of the molecule (shape: (N, 3))
"""
noon1 = self.get_noon(self.no_index1)
noon2 = self.get_noon(self.no_index2)
dnoon1 = self.get_dnoon_numerical(self.no_index1, noon1, coords, self.h)
dnoon2 = self.get_dnoon_numerical(self.no_index2, noon2, coords, self.h)
fmt = "%8.2f %16.9e\n"
if self.do_multistate:
with open("deltan_org.dat", "a") as f:
f.write(fmt % (self.step*self.tstep, noon2-noon1))
noon1, noon2, dnoon1, dnoon2 = self.offdiag.get_coupled_energies_grads(coords, noon1, noon2,
dnoon1, dnoon2,
step=self.step)
with open("deltan_meta.dat", "a") as f:
f.write(fmt % (self.step*self.tstep, noon2-noon1))
self.s = noon2-noon1
self.ds_dr = dnoon2-dnoon1
[docs]def get_norbitals_molpro():
with open("molpro.out") as f:
for line in f:
if "NUMBER OF CONTRACTIONS" in line:
return int(line.split()[3])
[docs]def get_coefficients_molpro(norbs, molproout):
def split_coefficients(line, ncols):
tmp = line.split()
while len(tmp) < ncols:
for k, t in enumerate(tmp):
if len(t) > 10:
tmp = tmp[:k] + [t[:-10]] + [t[-10:]] + tmp[k+1:]
break
return tmp
ncols = 10
coefficients = np.zeros((norbs, norbs))
with open(molproout) as f:
for line in f:
if "Orbital matrix ORB" in line:
while not " 1.1 " in line:
line = f.next()
for j in range(norbs):
if norbs > ncols:
coefficients[j, :ncols] = split_coefficients(line, ncols+1)[1:]
else:
coefficients[j, :norbs] = split_coefficients(line, norbs+1)[1:]
for i in range((norbs-ncols) / ncols):
coefficients[j, (i+1)*ncols:(i+2)*ncols] = split_coefficients(f.next(),
ncols)
if norbs % ncols and norbs > ncols:
coefficients[j, norbs/ncols*ncols:] = split_coefficients(f.next(),
norbs % ncols)
line = f.next()
line = f.next()
coefficients[j]
return coefficients.T
[docs]def get_overlap_molpro(norb, molproout):
with open(molproout) as f:
for line in f:
if "MATRIX S" in line:
f.next()
f.next()
# if norb <= 99:
# S = np.array([f.next().split() for i in range(norb)], dtype=float)
# else:
ncol = 99
nlines = norb / ncol
S = np.zeros((norb, norb))
for i in range(norb):
for j in range(nlines):
S[i][j*ncol:(j+1)*ncol] = f.next().split()
if norb % ncol:
S[i][nlines*ncol:] = f.next().split()
return S
[docs]def get_coefficients_ortho(C, S):
ev, evec = np.linalg.eigh(S)
s_sqrt = np.diag(ev**(-0.5))
U = evec
X = np.dot(U, s_sqrt)
return np.dot(np.linalg.inv(X), C)
[docs]def get_population_molpro(molproout):
with open(molproout) as f:
pop = np.array([line.split()[7] for line in f if "DUMP ORBITAL" in line], dtype=float)
return pop
[docs]def get_orbitals_molproxml(xml):
ns = "{http://www.molpro.net/schema/molpro-output}"
tree = et.parse(xml)
orbitals = [prop for prop in tree.iter(ns + "orbital")]
occupation = np.array([orbital.attrib["occupation"] for orbital in orbitals], dtype=float)
# energies = np.array([orbital.attrib["energy"] for orbital in orbitals], dtype=float)
coefficients = np.array([orbital.text.split() for orbital in orbitals], dtype=float).T
return occupation, coefficients
[docs]def getNatPop(nstates, norb, qcout, iface="molpro"):
""" Get NOONs and NO coefficients. """
if iface == "molpro":
# get natural orbital populations from Molpro output
if qcout == None or qcout == "molpro.out":
xml = "natural.xml"
else:
xml = qcout.split("/")[0] + "/natural.xml"
if os.path.exists(xml):
pop, C = get_orbitals_molproxml(xml)
else:
C = get_coefficients_molpro(norb, qcout)
pop = get_population_molpro(qcout)
all_ev = [(pop[::-1], C)]
return all_ev
[docs]def run_molpro((step, i, symbols, coords)):
os.chdir("molpro_%i_%i" % (step, i))
with open("molpro.in", "w") as f:
with open("head_noon") as g:
head = g.readlines()
f.writelines(head)
for sym, coord in zip(symbols, coords*c.a0):
f.write("%s %20.12f %20.12f %20.12f \n" % (sym, coord[0], coord[1], coord[2]))
with open("tail_noon") as g:
tail = g.readlines()
f.writelines(tail)
cwd = os.getcwd()
os.system("molpro -n 1 -d %s -I %s -W %s molpro.in" % (cwd, cwd, cwd))
os.chdir("..")