1# -*- coding: utf-8 -*- 

2# 

3# Copyright (c) 2017, the cclib development team 

4# 

5# This file is part of cclib (http://cclib.github.io) and is distributed under 

6# the terms of the BSD 3-Clause License. 

7 

8"""Calculation of Mayer's bond orders based on data parsed by cclib.""" 

9 

10import random 

11 

12import numpy 

13 

14from cclib.method.density import Density 

15 

16 

17class MBO(Density): 

18 """Mayer's bond orders""" 

19 

20 def __init__(self, *args): 

21 

22 # Call the __init__ method of the superclass. 

23 super(MBO, self).__init__(logname="MBO", *args) 

24 

25 def __str__(self): 

26 """Return a string representation of the object.""" 

27 return "Mayer's bond order of %s" % (self.data) 

28 

29 def __repr__(self): 

30 """Return a representation of the object.""" 

31 return 'Mayer\'s bond order("%s")' % (self.data) 

32 

33 def calculate(self, indices=None, fupdate=0.05): 

34 """Calculate Mayer's bond orders.""" 

35 

36 retval = super(MBO, self).calculate(fupdate) 

37 if not retval: #making density didn't work 

38 return False 

39 

40 # Do we have the needed info in the ccData object? 

41 if not (hasattr(self.data, "aooverlaps") 

42 or hasattr(self.data, "fooverlaps")): 

43 self.logger.error("Missing overlap matrix") 

44 return False #let the caller of function know we didn't finish 

45 

46 if not indices: 

47 

48 # Build list of groups of orbitals in each atom for atomresults. 

49 if hasattr(self.data, "aonames"): 

50 names = self.data.aonames 

51 overlaps = self.data.aooverlaps 

52 elif hasattr(self.data, "fonames"): 

53 names = self.data.fonames 

54 overlaps = self.data.fooverlaps 

55 else: 

56 self.logger.error("Missing aonames or fonames") 

57 return False 

58 

59 atoms = [] 

60 indices = [] 

61 

62 name = names[0].split('_')[0] 

63 atoms.append(name) 

64 indices.append([0]) 

65 

66 for i in range(1, len(names)): 

67 name = names[i].split('_')[0] 

68 try: 

69 index = atoms.index(name) 

70 except ValueError: #not found in atom list 

71 atoms.append(name) 

72 indices.append([i]) 

73 else: 

74 indices[index].append(i) 

75 

76 self.logger.info("Creating attribute fragresults: array[3]") 

77 size = len(indices) 

78 

79 # Determine number of steps, and whether process involves beta orbitals. 

80 PS = [] 

81 PS.append(numpy.dot(self.density[0], overlaps)) 

82 nstep = size**2 #approximately quadratic in size 

83 unrestricted = (len(self.data.mocoeffs) == 2) 

84 if unrestricted: 

85 self.fragresults = numpy.zeros([2, size, size], "d") 

86 PS.append(numpy.dot(self.density[1], overlaps)) 

87 else: 

88 self.fragresults = numpy.zeros([1, size, size], "d") 

89 

90 # Intialize progress if available. 

91 if self.progress: 

92 self.progress.initialize(nstep) 

93 

94 step = 0 

95 for i in range(len(indices)): 

96 

97 if self.progress and random.random() < fupdate: 

98 self.progress.update(step, "Mayer's Bond Order") 

99 

100 for j in range(i+1, len(indices)): 

101 

102 tempsumA = 0 

103 tempsumB = 0 

104 

105 for a in indices[i]: 

106 

107 for b in indices[j]: 

108 

109 if unrestricted: 

110 tempsumA += 2 * PS[0][a][b] * PS[0][b][a] 

111 tempsumB += 2 * PS[1][a][b] * PS[1][b][a] 

112 else: 

113 tempsumA += PS[0][a][b] * PS[0][b][a] 

114 

115 self.fragresults[0][i, j] = tempsumA 

116 self.fragresults[0][j, i] = tempsumA 

117 

118 if unrestricted: 

119 self.fragresults[1][i, j] = tempsumB 

120 self.fragresults[1][j, i] = tempsumB 

121 

122 if self.progress: 

123 self.progress.update(nstep, "Done") 

124 

125 return True