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

2# 

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

9 

10import random 

11 

12import numpy 

13 

14from cclib.method.population import Population 

15 

16 

17class MPA(Population): 

18 """Mulliken population analysis.""" 

19 

20 def __init__(self, *args): 

21 

22 # Call the __init__ method of the superclass. 

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

24 

25 def __str__(self): 

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

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

28 

29 def __repr__(self): 

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

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

32 

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

34 """Perform a Mulliken population analysis.""" 

35 

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

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

38 nbasis = self.data.nbasis 

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

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

41 nstep = alpha 

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

43 if unrestricted: 

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

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

46 nstep += beta 

47 

48 # Intialize progress if available. 

49 if self.progress: 

50 self.progress.initialize(nstep) 

51 

52 step = 0 

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

54 

55 for i in range(len(self.data.mocoeffs[spin])): 

56 

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

58 self.progress.update(step, "Mulliken Population Analysis") 

59 

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

61 # = c_{ai} \sum_b c_{bi} S_{ab} 

62 # = c_{ai} C(i) \cdot S(a) 

63 # X = C(i) * [C(i) \cdot S] 

64 # C(i) is 1xn and S is nxn, result of matrix mult is 1xn 

65 

66 ci = self.data.mocoeffs[spin][i] 

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

68 temp = numpy.dot(ci, self.data.aooverlaps) 

69 

70 # handle spin-unrestricted beta case 

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

72 temp = numpy.dot(ci, self.data.fooverlaps2) 

73 

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

75 temp = numpy.dot(ci, self.data.fooverlaps) 

76 

77 self.aoresults[spin][i] = numpy.multiply(ci, temp).astype("d") 

78 

79 step += 1 

80 

81 if self.progress: 

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

83 

84 retval = super(MPA, self).partition(indices) 

85 

86 if not retval: 

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

88 return False 

89 

90 # Create array for Mulliken charges. 

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

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

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

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

95 if unrestricted: 

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

97 

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

99 

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

101 

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

103 self.fragcharges = numpy.add(self.fragcharges, temp) 

104 if spin == 0: 

105 alpha = numpy.add(alpha, temp) 

106 elif spin == 1: 

107 beta = numpy.add(beta, temp) 

108 

109 if not unrestricted: 

110 self.fragcharges = numpy.multiply(self.fragcharges, 2) 

111 else: 

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

113 self.fragspins = numpy.subtract(alpha, beta) 

114 

115 return True