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