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 142