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 electric multipole moments based on data parsed by cclib.""" 

9 

10import sys 

11from collections.abc import Iterable 

12 

13import numpy 

14 

15from cclib.parser.utils import convertor 

16from cclib.method.calculationmethod import Method 

17 

18 

19class Moments(Method): 

20 """A class used to calculate electric multipole moments. 

21 

22 The obtained results are stored in `results` attribute as a 

23 dictionary whose keys denote the used charge population scheme. 

24 """ 

25 def __init__(self, data): 

26 self.required_attrs = ('atomcoords', 'atomcharges') 

27 self.results = {} 

28 

29 super(Moments, self).__init__(data) 

30 

31 def __str__(self): 

32 """Returns a string representation of the object.""" 

33 return "Multipole moments of %s" % (self.data) 

34 

35 def __repr__(self): 

36 """Returns a representation of the object.""" 

37 return 'Moments("%s")' % (self.data) 

38 

39 def _calculate_dipole(self, charges, coords, origin): 

40 """Calculate the dipole moment from the given atomic charges 

41 and their coordinates with respect to the origin. 

42 """ 

43 transl_coords_au = convertor(coords - origin, 'Angstrom', 'bohr') 

44 dipole = numpy.dot(charges, transl_coords_au) 

45 

46 return convertor(dipole, 'ebohr', 'Debye') 

47 

48 def _calculate_quadrupole(self, charges, coords, origin): 

49 """Calculate the traceless quadrupole moment from the given 

50 atomic charges and their coordinates with respect to the origin. 

51 """ 

52 transl_coords_au = convertor(coords - origin, 'Angstrom', 'bohr') 

53 

54 delta = numpy.eye(3) 

55 Q = numpy.zeros([3, 3]) 

56 for i in range(3): 

57 for j in range(3): 

58 for q, r in zip(charges, transl_coords_au): 

59 Q[i,j] += 1/2 * q * (3 * r[i] * r[j] - \ 

60 numpy.linalg.norm(r)**2 * delta[i,j]) 

61 

62 triu_idxs = numpy.triu_indices_from(Q) 

63 raveled_idxs = numpy.ravel_multi_index(triu_idxs, Q.shape) 

64 quadrupole = numpy.take(Q.flatten(), raveled_idxs) 

65 

66 return convertor(quadrupole, 'ebohr2', 'Buckingham') 

67 

68 def calculate(self, origin='nuccharge', population='mulliken', 

69 masses=None): 

70 """Calculate electric dipole and quadrupole moments using parsed 

71 partial atomic charges. 

72 

73 Inputs: 

74 origin - a choice of the origin of coordinate system. Can be 

75 either a three-element iterable or a string. If 

76 iterable, then it explicitly defines the origin (in 

77 Angstrom). If string, then the value can be any one of 

78 the following and it describes what is used as the 

79 origin: 

80 * 'nuccharge' -- center of positive nuclear charge 

81 * 'mass' -- center of mass 

82 population - a type of population analysis used to extract 

83 corresponding atomic charges from the output file. 

84 masses - if None, then use default atomic masses. Otherwise, 

85 the user-provided will be used. 

86 

87 Returns: 

88 A list where the first element is the origin of coordinates, 

89 while other elements are dipole and quadrupole moments 

90 expressed in terms of Debye and Buckingham units 

91 respectively. 

92 

93 Raises: 

94 ValueError when an argument with incorrect value or of 

95 inappropriate type is passed to a method. 

96 

97 Notes: 

98 To calculate the quadrupole moment the Buckingham definition 

99 [1]_ is chosen. Hirschfelder et al. [2]_ define it two times 

100 as much. 

101 

102 References: 

103 .. [1] Buckingham, A. D. (1959). Molecular quadrupole moments. 

104 Quarterly Reviews, Chemical Society, 13(3), 183. 

105 https://doi.org:10.1039/qr9591300183. 

106 .. [2] Hirschfelder J. O., Curtiss C. F. and Bird R. B. (1954). 

107 The Molecular Theory of Gases and Liquids. New York: Wiley. 

108 """ 

109 coords = self.data.atomcoords[-1] 

110 try: 

111 charges = self.data.atomcharges[population] 

112 except KeyError as e: 

113 msg = ("charges coming from requested population analysis" 

114 "scheme are not parsed") 

115 raise ValueError(msg, e) 

116 

117 if isinstance(origin, Iterable) and not isinstance(origin, str): 

118 origin_pos = numpy.asarray(origin) 

119 elif origin == 'nuccharge': 

120 origin_pos = numpy.average(coords, weights=self.data.atomnos, axis=0) 

121 elif origin == 'mass': 

122 if masses: 

123 atommasses = numpy.asarray(masses) 

124 else: 

125 try: 

126 atommasses = self.data.atommasses 

127 except AttributeError as e: 

128 msg = ("atomic masses were not parsed, consider provide " 

129 "'masses' argument instead") 

130 raise ValueError(msg, e) 

131 origin_pos = numpy.average(coords, weights=atommasses, axis=0) 

132 else: 

133 raise ValueError("{} is invalid value for 'origin'".format(origin)) 

134 

135 dipole = self._calculate_dipole(charges, coords, origin_pos) 

136 quadrupole = self._calculate_quadrupole(charges, coords, origin_pos) 

137 

138 rv = [origin_pos, dipole, quadrupole] 

139 self.results.update({population: rv}) 

140 

141 return rv