#!/usr/bin/env python
"""Contains a class for calculation of Mulliken charges from Turbomole oder DFTBaby outputs."""
from __future__ import print_function
import numpy as np
import os
try:
from DFTB import utils, DFTB2
except ImportError:
pass
[docs]class mulliken():
def __init__(self, symbols=None):
self.z = {"H": 1, "C": 6, "N": 7, "O": 8, "S": 16, "Cl": 17, "Br": 35}
self.basis = {"H": (), "C": (), "N": (), "O": (), "S": (), "Cl": (), "Br": ()}
if symbols:
self.symbols = symbols
else:
self.set_symbols()
self.is_openshell = False
self.set_S()
if self.is_openshell:
self.Da, self.Db = self.get_D_os()
self.D = self.Da+self.Db
else:
self.D = self.get_D()
self.set_order()
[docs] def set_symbols(self):
""" get atom symbols """
with open("coord") as f:
for line in f:
if "$coord" in line:
self.symbols = []
l = f.next().split()
while len(l) == 4:
self.symbols.append(l[3].capitalize())
l = f.next().split()
break
[docs] def set_S(self):
""" get overlap matrix S from tm.out """
def parse_basis(bas_str):
ns, np, nd, nf, ng = 0, 0, 0, 0, 0
if "s" in bas_str:
ns, bas_str = bas_str.split("s")
if "p" in bas_str:
np, bas_str = bas_str.split("p")
if "d" in bas_str:
nd, bas_str = bas_str.split("d")
if "f" in bas_str:
nf, bas_str = bas_str.split("f")
if "g" in bas_str:
nf, bas_str = bas_str.split("g")
return (int(ns), 3*int(np), 5*int(nd), 7*int(nf), 9*int(ng))
with open("tm.out") as f:
for line in f:
if "UHF" in line:
self.is_openshell = True
elif "basis set information" in line:
tmp = f.next()
while "type" not in tmp:
tmp = f.next()
f.next()
tmp = f.next().split()
while "----" not in tmp[0]:
self.basis[tmp[0].capitalize()] = parse_basis(tmp[-1].split("|")[0][1:])
tmp = f.next().split()
elif "mo occupation" in line:
f.next()
self.norb, self.nocc = [int(n) for n in f.next().split()[1:]]
elif "OVERLAP" in line:
overlap = []
l = f.next()
l = f.next()
while not "----" in l:
overlap += l.split()
l = f.next()
overlap = np.array(overlap, dtype=float)
break
S_tri = np.zeros((self.norb, self.norb))
S_tri[np.tril_indices(self.norb)] = overlap
self.S = S_tri.T+S_tri
np.fill_diagonal(self.S, np.diag(S_tri))
[docs] def get_coefficients(self, mofile):
coefficients = []
with open(mofile) as f:
for line in f:
if "scfmo" in line or "uhfmo" in line:
ncol, ndig = [int(n) for n in line.split()[2].split(".")[0][7:].split("d")]
elif "eigenvalue" in line:
ci = []
while len(ci) < self.norb:
tmp = f.next().replace("D", "E")
if len(tmp) > ndig+1:
ci += [tmp[i:i+ndig] for i in range(0, ncol*ndig, ndig)]
else:
ci.append(tmp)
coefficients.append(ci)
if len(coefficients) >= self.nocc:
break
coefficients = [c[:self.norb] for c in coefficients]
return np.array(coefficients, dtype=float)
[docs] def get_D(self):
""" take MO coefficients from mos file and calculate density matrix D"""
coeffs = self.get_coefficients("mos")
D = np.zeros((self.norb, self.norb))
for orb in range(self.nocc):
D += 2*np.prod(np.dstack(np.meshgrid(coeffs[orb], coeffs[orb])), axis=2)
return D
[docs] def get_D_os(self):
""" take MO coefficients from mos file and calculate density matrix D"""
coeffs_a = self.get_coefficients("alpha")
coeffs_b = self.get_coefficients("beta")
Da = np.zeros((self.norb, self.norb))
Db = np.zeros((self.norb, self.norb))
for orb in range(self.nocc):
Da += np.prod(np.dstack(np.meshgrid(coeffs_a[orb], coeffs_a[orb])), axis=2)
if orb < self.nocc-1:
Db += np.prod(np.dstack(np.meshgrid(coeffs_b[orb], coeffs_b[orb])), axis=2)
return Da, Db
[docs] def set_order(self):
""" set order of basis functions """
self.order = []
for i, s in enumerate(self.symbols):
self.order += [i]*sum(self.basis[s])
return self.order
[docs] def get_mulliken_charges(self):
""" calculate mulliken charges from D and S """
P = np.dot(self.D, self.S)
q = np.zeros(len(self.symbols))
q0 = np.array([self.z[s] for s in self.symbols])
for i in range(self.norb):
q[self.order[i]] += P[i, i]
return q0-q
[docs]def do_dft_calculation(symbols, coord):
with open("coord","w") as f:
f.write("$coord\n")
for s, c in zip(symbols, coord):
f.write("%20.12f %20.12f %20.12f %s \n" %(c[0],c[1],c[2],s))
f.write("$end")
ridft = False
with open("control") as f:
for line in f:
if "$rij" in line:
ridft = True
break
if ridft:
os.system("ridft > tm.out")
else:
os.system("dscf > tm.out")
[docs]def get_charges(symbols=None):
"""
Calculate Mulliken charges from Turbomole ridft calculations.
Files mos (alpha, beta) and tm.out are needed in the execution directory. Please include
the keyword "$intsdebug cao" in the control file to print out the overlap matrix. "$scfconv"
should be increased to raise the accuracy of MO coefficients. With the default $scfconv=6,
charges can be obtained much faster directly from Turbomole in similar accuracy.
Parameters
----------
symbols: list of strings (optional)
if symbols are not passed, they are extracted from the coord file
Returns
-------
charges: numpy array
Mulliken charges, shape (nat)
Note
----
This program does not work with Turbomole 6.4, because the overlap matrix is not printed.
"""
os.chdir("tm")
mul = mulliken(symbols)
os.chdir("..")
return mul.get_mulliken_charges()
[docs]def get_new_charges_tm(args):
i, symbols, coords = args
os.system("cp -r tm tm_"+str(i))
os.chdir("tm_"+str(i))
do_dft_calculation(symbols, coords)
mul = mulliken(symbols)
charges = mul.get_mulliken_charges()
os.chdir("..")
os.system("rm -r tm_%i" % i)
return charges
[docs]def get_new_charges_dftb(args):
i, symbols, coords = args
elements={"H": 1, "C": 6, "N": 7, "O": 8, "Cl": 17, "Br": 35, "Ru": 44}
atomlist = [(elements[sym.capitalize()], tuple(coord)) for sym, coord in zip(symbols, coords)]
parser = utils.OptionParserFuncWrapper([DFTB2.DFTB2.__init__], "")
(options, args) = parser.parse_args()
dftb = DFTB2.DFTB2(atomlist, **options)
dftb.setGeometry(atomlist)
options = {}
dftb.getEnergy(**options)
return dftb.dq
if __name__ == "__main__":
print(get_charges())