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"""Löwdin population analysis.""" 

9 

10import random 

11 

12import numpy 

13 

14from cclib.method.population import Population 

15 

16 

17class LPA(Population): 

18 """The Löwdin population analysis""" 

19 def __init__(self, *args): 

20 

21 # Call the __init__ method of the superclass. 

22 super(LPA, self).__init__(logname="LPA", *args) 

23 

24 def __str__(self): 

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

26 return "LPA of %s" % (self.data) 

27 

28 def __repr__(self): 

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

30 return 'LPA("%s")' % (self.data) 

31 

32 def calculate(self, indices=None, x=0.5, fupdate=0.05): 

33 """Perform a calculation of Löwdin population analysis. 

34 

35 Inputs: 

36 indices - list of lists containing atomic orbital indices of fragments 

37 x - overlap matrix exponent in wavefunxtion projection (x=0.5 for Lowdin) 

38 """ 

39 

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

41 nbasis = self.data.nbasis 

42 

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

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

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

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

47 nstep = alpha 

48 

49 if unrestricted: 

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

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

52 nstep += beta 

53 

54 # intialize progress if available 

55 if self.progress: 

56 self.progress.initialize(nstep) 

57 

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

59 S = self.data.aooverlaps 

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

61 S = self.data.fooverlaps 

62 

63 # Get eigenvalues and matrix of eigenvectors for transformation decomposition (U). 

64 # Find roots of diagonal elements, and transform backwards using eigevectors. 

65 # We need two matrices here, one for S^x, another for S^(1-x). 

66 # We don't need to invert U, since S is symmetrical. 

67 eigenvalues, U = numpy.linalg.eig(S) 

68 UI = U.transpose() 

69 Sdiagroot1 = numpy.identity(len(S))*numpy.power(eigenvalues, x) 

70 Sdiagroot2 = numpy.identity(len(S))*numpy.power(eigenvalues, 1-x) 

71 Sroot1 = numpy.dot(U, numpy.dot(Sdiagroot1, UI)) 

72 Sroot2 = numpy.dot(U, numpy.dot(Sdiagroot2, UI)) 

73 

74 step = 0 

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

76 

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

78 

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

80 self.progress.update(step, "Lowdin Population Analysis") 

81 

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

83 

84 temp1 = numpy.dot(ci, Sroot1) 

85 temp2 = numpy.dot(ci, Sroot2) 

86 self.aoresults[spin][i] = numpy.multiply(temp1, temp2).astype("d") 

87 

88 step += 1 

89 

90 if self.progress: 

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

92 

93 retval = super(LPA, self).partition(indices) 

94 

95 if not retval: 

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

97 return False 

98 

99 # Create array for charges. 

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

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

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

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

104 if unrestricted: 

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

106 

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

108 

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

110 

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

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

113 if spin == 0: 

114 alpha = numpy.add(alpha, temp) 

115 elif spin == 1: 

116 beta = numpy.add(beta, temp) 

117 

118 if not unrestricted: 

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

120 else: 

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

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

123 

124 return True