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"""Charge Decomposition Analysis (CDA)""" 

9 

10import random 

11 

12import numpy 

13 

14from cclib.method.fragments import FragmentAnalysis 

15 

16 

17class CDA(FragmentAnalysis): 

18 """Charge Decomposition Analysis (CDA)""" 

19 

20 def __init__(self, *args): 

21 

22 # Call the __init__ method of the superclass. 

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

24 

25 def __str__(self): 

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

27 return "CDA of %s" % (self.data) 

28 

29 def __repr__(self): 

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

31 return 'CDA("%s")' % (self.data) 

32 

33 def calculate(self, fragments, cupdate=0.05): 

34 """Perform the charge decomposition analysis. 

35 

36 Inputs: 

37 fragments - list of ccData data objects 

38 """ 

39 

40 

41 retval = super(CDA, self).calculate(fragments, cupdate) 

42 if not retval: 

43 return False 

44 

45 # At this point, there should be a mocoeffs and fooverlaps 

46 # in analogy to a ccData object. 

47 

48 donations = [] 

49 bdonations = [] 

50 repulsions = [] 

51 residuals = [] 

52 

53 if len(self.mocoeffs) == 2: 

54 occs = 1 

55 else: 

56 occs = 2 

57 

58 # Intialize progress if available. 

59 nstep = self.data.homos[0] 

60 if len(self.data.homos) == 2: 

61 nstep += self.data.homos[1] 

62 if self.progress: 

63 self.progress.initialize(nstep) 

64 

65 # Begin the actual method. 

66 step = 0 

67 for spin in range(len(self.mocoeffs)): 

68 

69 size = len(self.mocoeffs[spin]) 

70 homo = self.data.homos[spin] 

71 

72 if len(fragments[0].homos) == 2: 

73 homoa = fragments[0].homos[spin] 

74 else: 

75 homoa = fragments[0].homos[0] 

76 

77 if len(fragments[1].homos) == 2: 

78 homob = fragments[1].homos[spin] 

79 else: 

80 homob = fragments[1].homos[0] 

81 

82 self.logger.info("handling spin unrestricted") 

83 if spin == 0: 

84 fooverlaps = self.fooverlaps 

85 elif spin == 1 and hasattr(self, "fooverlaps2"): 

86 fooverlaps = self.fooverlaps2 

87 

88 offset = fragments[0].nbasis 

89 

90 self.logger.info("Creating donations, bdonations, and repulsions: array[]") 

91 donations.append(numpy.zeros(size, "d")) 

92 bdonations.append(numpy.zeros(size, "d")) 

93 repulsions.append(numpy.zeros(size, "d")) 

94 residuals.append(numpy.zeros(size, "d")) 

95 

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

97 

98 # Calculate donation for each MO. 

99 for k in range(0, homoa + 1): 

100 for n in range(offset + homob + 1, self.data.nbasis): 

101 donations[spin][i] += 2 * occs * self.mocoeffs[spin][i,k] \ 

102 * self.mocoeffs[spin][i,n] * fooverlaps[k][n] 

103 

104 for l in range(offset, offset + homob + 1): 

105 for m in range(homoa + 1, offset): 

106 bdonations[spin][i] += 2 * occs * self.mocoeffs[spin][i,l] \ 

107 * self.mocoeffs[spin][i,m] * fooverlaps[l][m] 

108 

109 for k in range(0, homoa + 1): 

110 for m in range(offset, offset+homob + 1): 

111 repulsions[spin][i] += 2 * occs * self.mocoeffs[spin][i,k] \ 

112 * self.mocoeffs[spin][i, m] * fooverlaps[k][m] 

113 

114 for m in range(homoa + 1, offset): 

115 for n in range(offset + homob + 1, self.data.nbasis): 

116 residuals[spin][i] += 2 * occs * self.mocoeffs[spin][i,m] \ 

117 * self.mocoeffs[spin][i, n] * fooverlaps[m][n] 

118 

119 step += 1 

120 if self.progress and random.random() < cupdate: 

121 self.progress.update(step, "Charge Decomposition Analysis...") 

122 

123 if self.progress: 

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

125 

126 self.donations = donations 

127 self.bdonations = bdonations 

128 self.repulsions = repulsions 

129 self.residuals = residuals 

130 

131 return True