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

2# 

3# Copyright (c) 2020, 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 Bickelhaupt population analysis based on data parsed by cclib.""" 

9 

10import random 

11 

12import numpy 

13 

14from cclib.method.population import Population 

15 

16 

17class Bickelhaupt(Population): 

18 """Bickelhaupt population analysis.""" 

19 

20 def __init__(self, *args): 

21 

22 # Call the __init__ method of the superclass. 

23 super(Bickelhaupt, self).__init__(logname="Bickelhaupt Population Analysis", *args) 

24 

25 def __str__(self): 

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

27 return "Bickelhaupt charges of {}".format(self.data) 

28 

29 def __repr__(self): 

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

31 return "Bickelhaupt({})".format(self.data) 

32 

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

34 """Perform a Bickelhaupt population analysis.""" 

35 # Bickelhaupt population analysis uses the relative magnitude of the diagonal terms to 

36 # partition off-diagonal terms. The weights are therefore calculated using the following 

37 # formula: 

38 # W_{ab} = 2 * \sum_k c_{ak}^2 / (\sum_i c_{ai}^2 + \sum_j c_{bj}^2) 

39 # The dimension of W is (nbasis)x(nbasis). 

40 # Then, each term is calculated similarly to Mulliken charges, in a form of 

41 # X_{ai} = \sum_b W_{ab} c_{ai} c_{bi} S_{ab} 

42 # This agrees with quation 11 on 10.1021/om950966x, with the charge (q) terms substituted 

43 # with MO coefficients. 

44 

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

46 self.logger.info("Creating attribute aoresults: [array[2]]") 

47 nbasis = self.data.nbasis 

48 alpha = len(self.data.mocoeffs[0]) 

49 self.aoresults = [numpy.zeros([alpha, nbasis], "d")] 

50 w = [numpy.zeros([self.data.nbasis, self.data.nbasis], "d")] 

51 nstep = alpha + nbasis 

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

53 if unrestricted: 

54 beta = len(self.data.mocoeffs[1]) 

55 self.aoresults.append(numpy.zeros([beta, nbasis], "d")) 

56 w.append(numpy.zeros([self.data.nbasis, self.data.nbasis], "d")) 

57 nstep += beta + nbasis 

58 

59 # Intialize progress if available. 

60 if self.progress: 

61 self.progress.initialize(nstep) 

62 

63 step = 0 

64 

65 # First determine the weight matrix 

66 for spin in range(len(self.data.mocoeffs)): 

67 for i in range(self.data.nbasis): 

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

69 self.progress.update(step, "Bickelhaupt Population Analysis") 

70 for j in range(self.data.nbasis): 

71 aTerm = numpy.dot( 

72 self.data.mocoeffs[spin][0 : (self.data.homos[spin] + 1), i], 

73 self.data.mocoeffs[spin][0 : (self.data.homos[spin] + 1), i], 

74 ) 

75 bTerm = numpy.dot( 

76 self.data.mocoeffs[spin][0 : (self.data.homos[spin] + 1), j], 

77 self.data.mocoeffs[spin][0 : (self.data.homos[spin] + 1), j], 

78 ) 

79 w[spin][i][j] = 2 * aTerm / (aTerm + bTerm) 

80 step += 1 

81 

82 for spin in range(len(self.data.mocoeffs)): 

83 for i in range(len(self.data.mocoeffs[spin])): # for each mo 

84 for a in range(self.data.nbasis): # for each basis 

85 if hasattr(self.data, "aooverlaps"): 

86 overlaps = self.data.aooverlaps[a] 

87 # handle spin-unrestricted beta case 

88 elif hasattr(self.data, "fooverlaps2") and spin == 1: 

89 overlaps = self.data.fooverlaps2[a] 

90 

91 elif hasattr(self.data, "fooverlaps"): 

92 overlaps = self.data.fooverlaps[a] 

93 

94 self.aoresults[spin][i][a] = numpy.sum( 

95 ( 

96 (w[spin][a] * self.data.mocoeffs[spin][i][a]) 

97 * self.data.mocoeffs[spin][i] 

98 ) 

99 * overlaps 

100 ) 

101 

102 step += 1 

103 

104 if self.progress: 

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

106 

107 retval = super(Bickelhaupt, self).partition(indices) 

108 

109 if not retval: 

110 self.logger.error("Error in partitioning results") 

111 return False 

112 

113 # Create array for Bickelhaupt charges. 

114 self.logger.info("Creating fragcharges: array[1]") 

115 size = len(self.fragresults[0][0]) 

116 self.fragcharges = numpy.zeros([size], "d") 

117 alpha = numpy.zeros([size], "d") 

118 if unrestricted: 

119 beta = numpy.zeros([size], "d") 

120 

121 for spin in range(len(self.fragresults)): 

122 

123 for i in range(self.data.homos[spin] + 1): 

124 

125 temp = numpy.reshape(self.fragresults[spin][i], (size,)) 

126 self.fragcharges = self.fragcharges + temp 

127 if spin == 0: 

128 alpha += temp 

129 elif spin == 1: 

130 beta += temp 

131 

132 if not unrestricted: 

133 self.fragcharges = self.fragcharges * 2 

134 else: 

135 self.logger.info("Creating fragspins: array[1]") 

136 self.fragspins = alpha - beta 

137 

138 return True