# -*- 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 pickle
from .cv import cv
from .. import constants as c
[docs]class cv_energygap(cv):
"""
Energy gap between two states 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
-----------------
state1 : int
index of first electronic state (ground state = 0)
state2 : int
index of second electronic state (ground state = 0)
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., 1.0, 200)
self.state1 = self.get_kwargs("state1") # lower state index, ground state = 0
self.state2 = self.get_kwargs("state2") # upper state index, ground state = 0
self.do_multistate = self.get_kwargs("do_multistate", default=False)
self.width *= c.from_energyunit(self.get_kwargs("width_unit", default="eV"))
[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("deltaE_meta.dat") and self.restart):
with open("deltaE_BO.dat", "w") as f:
f.write("# energy gap / hartree\n")
with open("deltaE_meta.dat", "w") as f:
f.write("# energy gap / hartree\n")
# with open("ddeltaE_meta.dat", "w") as f:
# f.write("derivative of gap\n")
[docs] def set_s_and_ds(self, coords):
"""
Calculate energy gap and its gradient for the current set of coordinates.
Parameters
----------
coords : np.2darray
cartesian coordinates of the molecule (shape: (N, 3))
"""
# if self.iface == "om2":
# try:
# reader = om2reader("fort.15")
# reader.getEnergies()
# en1 = reader.getEnergy(self.state1)*kcalmol2au
# en2 = reader.getEnergy(self.state2)*kcalmol2au
# reader.getGradients()
# grad1 = reader.getGradient(self.state1)*kcalmol2au*au2ang
# grad2 = reader.getGradient(self.state2)*kcalmol2au*au2ang
# except AttributeError:
# beststep = np.argmin(np.loadtxt("deltaE_meta.dat"))
# print "step with lowest energy gap: ", beststep
# print "saving coords to ci_%d.npy" % (self.ci_index+1)
# s, bestcoords, v = initial.get_initial_from_log(beststep)
## sqdm = get_squared_distancematrix(bestcoords, swap = self.swap, ref = self.ref)
# np.save("ci_%d.npy" % (self.ci_index+1), bestcoords)
# return 0.0, np.zeros_like(coords)
if self.iface == "gaussian": # only CASSCF calculations
nat = coords.shape[0]
en = []
grad = []
def read_energy(i, en):
with open(str(i)+"/out", "r") as f:
for line in f:
if "EIGENVALUE " in line:
tmp = line.split()
if tmp[1] == str(i+1)+")" or tmp[0] == "("+str(i+1)+")":
en.append(float(tmp[-1]))
return en
for i in range(2):
en = read_energy(i, en)
if len(en) < i+1:
os.system("cp %d/inp.chk %d/inp.chk" % ((i+1) % 2, i))
os.chdir(str(i))
os.system("g09 < inp > out")
os.chdir("..")
en = read_energy(i, en)
grad.append([])
with open(str(i)+"/Test.FChk") as f:
for line in f:
if "Cartesian Gradient" in line:
while len(grad[i]) < 3*nat:
grad[i] += next(f).split()
break
grad[i] = np.array(grad[i], dtype=float).reshape((nat, 3))
en1 = en[self.state1]
en2 = en[self.state2]
grad1 = grad[self.state1]
grad2 = grad[self.state2]
elif self.iface == "turbomole":
# from ..interfaces.turbomole import iface_turbomole
# iface = iface_turbomole(self.symbols, coords, directory="tm_state2")
# iface.calculate()
# en2 = iface.get_energy()
# grad2 = iface.get_gradient()
with open("tm/energies.pickle") as f, open("tm/gradients.pickle") as g:
en1 = pickle.load(f)
grad1 = pickle.load(g)
with open("tm_state2/energies.pickle") as f, open("tm_state2/gradients.pickle") as g:
en2 = pickle.load(f)
grad2 = pickle.load(g)
if en1 > en2:
en = en1
grad = grad1
en1 = en2
grad1 = grad2
en2 = en
grad2 = grad
print("ENERGIES: ", en1, en2)
else:
with open("energies.pickle") as f, open("gradients.pickle") as g:
energies = pickle.load(f)
gradients = pickle.load(g)
if self.iface == "molpro":
# index1 = str(self.state1+1)
# index2 = str(self.state2+1)
index1 = self.state1
index2 = self.state2
en1 = energies[index1]
en2 = energies[index2]
grad1 = gradients[index1]
grad2 = gradients[index2]
fmt = "%8.2f %16.8f\n"
if self.do_multistate:
# print("energies before coupling: \n", en2, en1)
# oldgap = en2-en1
with open("deltaE_BO.dat", "a") as f:
f.write(fmt % (self.step*self.tstep, en2-en1))
# f.write("%16.8f\n" % (en2-en1))
# en1, en2, grad1, grad2 = self.couple_energies_and_grads(en1, en2, grad1, grad2, coords)
en1, en2, grad1, grad2 = self.offdiag.get_coupled_energies_grads(coords, en1, en2,
grad1, grad2,
step=self.step)
with open("deltaE_meta.dat", "a") as f:
f.write(fmt % (self.step*self.tstep, en2-en1))
# print("gap after coupling: \n", en2, en1)
# if en2-en1 <= self.ci_thresh:
## sqdm = get_squared_distancematrix(coords, swap = self.swap, ref = self.ref)
# np.save("ci_%d.npy" % (self.ci_index+1), coords)
# newgap = en2-en1
# with open("coupling.dat", "a") as f:
# f.write("%12.6f\n" % (newgap-oldgap))
gap = en2-en1
dgap = grad2-grad1
# for s, r in zip(self.steps, result):
# print("append gap to deltaE_meta.dat: ", gap)
# f.write("%16.8f\n" % gap)
# with open("ddeltaE_meta.dat", "a") as f:
# f.write("%16.8f\n" % np.linalg.norm(dgap))
self.s = gap
self.ds_dr = dgap