1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
|
# -*- coding: utf-8 -*-
#
# Copyright (c) 2018, the cclib development team
#
# This file is part of cclib (http://cclib.github.io) and is distributed under
# the terms of the BSD 3-Clause License.
"""Calculation of overlap population analysis based on cclib data."""
import random
import numpy
from cclib.method.calculationmethod import Method
from cclib.method.population import Population
def func(x):
if x==1:
return 1
else:
return x+func(x-1)
class OPA(Population):
"""Overlap population analysis."""
def __init__(self, *args):
# Call the __init__ method of the superclass.
super(OPA, self).__init__(logname="OPA", *args)
def __str__(self):
"""Return a string representation of the object."""
return "OPA of %s" % (self.data)
def __repr__(self):
"""Return a representation of the object."""
return 'OPA("%s")' % (self.data)
def calculate(self, indices=None, fupdate=0.05):
"""Perform an overlap population analysis given the results of a parser"""
if not indices:
# Build list of groups of orbitals in each atom for atomresults.
if hasattr(self.data, "aonames"):
names = self.data.aonames
elif hasattr(self.data, "foonames"):
names = self.data.fonames
atoms = []
indices = []
name = names[0].split('_')[0]
atoms.append(name)
indices.append([0])
for i in range(1, len(names)):
name = names[i].split('_')[0]
try:
index = atoms.index(name)
except ValueError: #not found in atom list
atoms.append(name)
indices.append([i])
else:
indices[index].append(i)
# Determine number of steps, and whether process involves beta orbitals.
nfrag = len(indices) #nfrag
nstep = func(nfrag - 1)
unrestricted = (len(self.data.mocoeffs) == 2)
alpha = len(self.data.mocoeffs[0])
nbasis = self.data.nbasis
self.logger.info("Creating attribute results: array[4]")
results= [ numpy.zeros([nfrag, nfrag, alpha], "d") ]
if unrestricted:
beta = len(self.data.mocoeffs[1])
results.append(numpy.zeros([nfrag, nfrag, beta], "d"))
nstep *= 2
if hasattr(self.data, "aooverlaps"):
overlap = self.data.aooverlaps
elif hasattr(self.data,"fooverlaps"):
overlap = self.data.fooverlaps
#intialize progress if available
if self.progress:
self.progress.initialize(nstep)
size = len(self.data.mocoeffs[0])
step = 0
preresults = []
for spin in range(len(self.data.mocoeffs)):
two = numpy.array([2.0]*len(self.data.mocoeffs[spin]),"d")
# OP_{AB,i} = \sum_{a in A} \sum_{b in B} 2 c_{ai} c_{bi} S_{ab}
for A in range(len(indices)-1):
for B in range(A+1, len(indices)):
if self.progress: #usually only a handful of updates, so remove random part
self.progress.update(step, "Overlap Population Analysis")
for a in indices[A]:
ca = self.data.mocoeffs[spin][:,a]
for b in indices[B]:
cb = self.data.mocoeffs[spin][:,b]
temp = ca * cb * two *overlap[a,b]
results[spin][A,B] = numpy.add(results[spin][A,B],temp)
results[spin][B,A] = numpy.add(results[spin][B,A],temp)
step += 1
temparray2 = numpy.swapaxes(results[0],1,2)
self.results = [ numpy.swapaxes(temparray2,0,1) ]
if unrestricted:
temparray2 = numpy.swapaxes(results[1],1,2)
self.results.append(numpy.swapaxes(temparray2, 0, 1))
if self.progress:
self.progress.update(nstep, "Done")
return True
|
