# -*- coding: utf-8 -*-
"""Contains a class for using a Natural Orbital Occupation Number (NOON) 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 atoms as at
from .. import constants as c
[docs]class cv_noon(cv):
"""
NOON 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_index : int
index of natural orbital (sorted by rising occupation number (0 = lowest occupation)
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)
iface : str
interface to QC code
"""
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_index = self.get_kwargs("no_index")
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()
[docs] def initialize(self):
cv.initialize(self)
if not (os.path.exists("noon_%i.dat" % self.no_index) and self.restart):
with open("noon_%i.dat" % self.no_index, "w") as f:
f.write("# Natural Orbital Occupation Number (NOON) \n")
with open("noons_all.dat", "w") as f:
f.write("# Natural Orbital Occupation Numbers (NOONs) \n")
[docs] def get_noon(self, 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" and self.idx == 0:
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-self.no_index-1]
[docs] def get_dnoon_numerical(self, 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(qcout="molpro_%i_%i/molpro.out" % (self.step, i)) for i
in range(coords.size)]).reshape(coords.shape)
return (s_h-self.s)/h
[docs] def set_s_and_ds(self, coords):
"""
Calculate noon and its gradient for the current set of coordinates.
Parameters
----------
coords : np.2darray
cartesian coordinates of the molecule (shape: (N, 3))
"""
self.s = self.get_noon()
self.ds_dr = self.get_dnoon_numerical(coords, self.h)
fmt = "%8.2f %16.9e\n"
with open("noon_%i.dat" % self.no_index, "a") as f:
f.write(fmt % (self.step*self.tstep, self.s))
[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_coefficients_deortho(C_prime, 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(X, C_prime)
[docs]def get_population_molpro(molproout):
### attention: be careful that only one molden dump records is included in molpro.out
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)
# S = get_overlap_molpro(norb, qcout)
# C = get_coefficients_ortho(C, S)
# rho = [get_density_from_coeffs(C, pop)]
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("..")