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"""Fragment analysis based on parsed ADF data.""" 

9 

10import logging 

11import random 

12 

13import numpy 

14numpy.inv = numpy.linalg.inv 

15 

16from cclib.method.calculationmethod import Method 

17 

18 

19class FragmentAnalysis(Method): 

20 """Convert a molecule's basis functions from atomic-based to fragment MO-based""" 

21 def __init__(self, data, progress=None, loglevel=logging.INFO, 

22 logname="FragmentAnalysis of"): 

23 

24 # Call the __init__ method of the superclass. 

25 super(FragmentAnalysis, self).__init__(data, progress, loglevel, logname) 

26 self.parsed = False 

27 

28 def __str__(self): 

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

30 return "Fragment molecule basis of %s" % (self.data) 

31 

32 def __repr__(self): 

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

34 return 'Fragment molecular basis("%s")' % (self.data) 

35 

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

37 

38 nFragBasis = 0 

39 nFragAlpha = 0 

40 nFragBeta = 0 

41 self.fonames = [] 

42 

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

44 

45 self.logger.info("Creating attribute fonames[]") 

46 

47 # Collect basis info on the fragments. 

48 for j in range(len(fragments)): 

49 nFragBasis += fragments[j].nbasis 

50 nFragAlpha += fragments[j].homos[0] + 1 

51 if unrestricted and len(fragments[j].homos) == 1: 

52 nFragBeta += fragments[j].homos[0] + 1 #assume restricted fragment 

53 elif unrestricted and len(fragments[j].homos) == 2: 

54 nFragBeta += fragments[j].homos[1] + 1 #assume unrestricted fragment 

55 

56 #assign fonames based on fragment name and MO number 

57 for i in range(fragments[j].nbasis): 

58 if hasattr(fragments[j],"name"): 

59 self.fonames.append("%s_%i"%(fragments[j].name,i+1)) 

60 else: 

61 self.fonames.append("noname%i_%i"%(j,i+1)) 

62 

63 nBasis = self.data.nbasis 

64 nAlpha = self.data.homos[0] + 1 

65 if unrestricted: 

66 nBeta = self.data.homos[1] + 1 

67 

68 # Check to make sure calcs have the right properties. 

69 if nBasis != nFragBasis: 

70 self.logger.error("Basis functions don't match") 

71 return False 

72 

73 if nAlpha != nFragAlpha: 

74 self.logger.error("Alpha electrons don't match") 

75 return False 

76 

77 if unrestricted and nBeta != nFragBeta: 

78 self.logger.error("Beta electrons don't match") 

79 return False 

80 

81 if len(self.data.atomcoords) != 1: 

82 self.logger.warning("Molecule calc appears to be an optimization") 

83 

84 for frag in fragments: 

85 if len(frag.atomcoords) != 1: 

86 msg = "One or more fragment appears to be an optimization" 

87 self.logger.warning(msg) 

88 break 

89 

90 last = 0 

91 for frag in fragments: 

92 size = frag.natom 

93 if self.data.atomcoords[0][last:last+size].tolist() != \ 

94 frag.atomcoords[0].tolist(): 

95 self.logger.error("Atom coordinates aren't aligned") 

96 return False 

97 if self.data.atomnos[last:last+size].tolist() != \ 

98 frag.atomnos.tolist(): 

99 self.logger.error("Elements don't match") 

100 return False 

101 

102 last += size 

103 

104 # And let's begin! 

105 self.mocoeffs = [] 

106 self.logger.info("Creating mocoeffs in new fragment MO basis: mocoeffs[]") 

107 

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

109 blockMatrix = numpy.zeros((nBasis,nBasis), "d") 

110 pos = 0 

111 

112 # Build up block-diagonal matrix from fragment mocoeffs. 

113 # Need to switch ordering from [mo,ao] to [ao,mo]. 

114 for i in range(len(fragments)): 

115 size = fragments[i].nbasis 

116 if len(fragments[i].mocoeffs) == 1: 

117 temp = numpy.transpose(fragments[i].mocoeffs[0]) 

118 blockMatrix[pos:pos+size, pos:pos+size] = temp 

119 else: 

120 temp = numpy.transpose(fragments[i].mocoeffs[spin]) 

121 blockMatrix[pos:pos+size, pos:pos+size] = temp 

122 pos += size 

123 

124 # Invert and mutliply to result in fragment MOs as basis. 

125 iBlockMatrix = numpy.inv(blockMatrix) 

126 temp = numpy.transpose(self.data.mocoeffs[spin]) 

127 results = numpy.transpose(numpy.dot(iBlockMatrix, temp)) 

128 

129 self.mocoeffs.append(results) 

130 

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

132 tempMatrix = numpy.dot(self.data.aooverlaps, blockMatrix) 

133 tBlockMatrix = numpy.transpose(blockMatrix) 

134 if spin == 0: 

135 self.fooverlaps = numpy.dot(tBlockMatrix, tempMatrix) 

136 self.logger.info("Creating fooverlaps: array[x,y]") 

137 elif spin == 1: 

138 self.fooverlaps2 = numpy.dot(tBlockMatrix, tempMatrix) 

139 self.logger.info("Creating fooverlaps (beta): array[x,y]") 

140 else: 

141 self.logger.warning("Overlap matrix missing") 

142 

143 self.parsed = True 

144 self.nbasis = nBasis 

145 self.homos = self.data.homos 

146 

147 return True