1# -*- coding: utf-8 -*-
2# MolMod is a collection of molecular modelling tools for python.
3# Copyright (C) 2007 - 2019 Toon Verstraelen <Toon.Verstraelen@UGent.be>, Center
4# for Molecular Modeling (CMM), Ghent University, Ghent, Belgium; all rights
5# reserved unless otherwise stated.
6#
7# This file is part of MolMod.
8#
9# MolMod is free software; you can redistribute it and/or
10# modify it under the terms of the GNU General Public License
11# as published by the Free Software Foundation; either version 3
12# of the License, or (at your option) any later version.
13#
14# MolMod is distributed in the hope that it will be useful,
15# but WITHOUT ANY WARRANTY; without even the implied warranty of
16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
17# GNU General Public License for more details.
18#
19# You should have received a copy of the GNU General Public License
20# along with this program; if not, see <http://www.gnu.org/licenses/>
21#
22# --
23"""Evaluation of internal coordinates with first and second order derivatives
24
25This implementation is pure python and sacrifices computational efficiency on
26the altar of programming flexibility. It is really easy to implement new types
27of internal coordinates since one only has to enter the formula that evaluates
28the internal coordinate. First and second order derivatives towards Cartesian
29coordinates require only a minimum of extra work.
30
31Two auxiliary classes Scalar and Vector3 support most of the mathematical
32operations required to compute the internal coordinates. Additionally they also
33know the chain rule for each operation and can therefore evaluate the
34derivatives simultaneously.
35"""
36
37
38from __future__ import division
39
40import numpy as np
41
42
43__all__ = [
44    "Scalar", "Vector3", "dot", "cross",
45    "bond_length", "pair_distance",
46    "bend_cos", "bend_angle",
47    "dihed_cos", "dihed_angle",
48    "opbend_cos", "opbend_angle", "opbend_dist",
49    "opbend_mcos", "opbend_mangle",
50]
51
52
53#
54# Auxiliary classes
55#
56
57class Scalar(object):
58    """A scalar object with optional first and second order derivates
59
60       Each input value to which the derivative is computed has its own index.
61       The numerical value of the derivatives are stored in arrays self.d and
62       self.dd. The value of the scalar itself if self.v
63    """
64
65    def __init__(self, size, deriv=0, value=0, index=None):
66        """
67           Arguments:
68            | ``size`` -- The number of inputs on which this ic depends. e.g. a
69                          distance depends on 6 Cartesian coordinates.
70            | ``deriv`` -- Consider up to deriv order derivatives. (max=2)
71            | ``value`` -- The initial value.
72            | ``index`` -- If this scalar is one of the input variables, this is
73                           its index.
74
75           The scalar object supports several in place modifications.
76        """
77        self.deriv = deriv
78        self.size = size
79        self.v = value
80        if deriv > 0:
81            self.d = np.zeros(size, float)
82            if index is not None:
83                self.d[index] = 1
84        if deriv > 1:
85            self.dd = np.zeros((size, size), float)
86        if deriv > 2:
87            raise ValueError("This implementation (only) supports up to second order derivatives.")
88
89    def copy(self):
90        """Return a deep copy"""
91        result = Scalar(self.size, self.deriv)
92        result.v = self.v
93        if self.deriv > 0: result.d[:] = self.d[:]
94        if self.deriv > 1: result.dd[:] = self.dd[:]
95        return result
96
97    def results(self):
98        """Return the value and optionally derivative and second order derivative"""
99        if self.deriv == 0:
100            return self.v,
101        if self.deriv == 1:
102            return self.v, self.d
103        if self.deriv == 2:
104            return self.v, self.d, self.dd
105
106    def __iadd__(self, other):
107        if self.deriv > 1: self.dd += other.dd
108        if self.deriv > 0: self.d += other.d
109        self.v += other.v
110        return self
111
112    def __add__(self, other):
113        result = self.copy()
114        result += other
115        return result
116
117    def __isub__(self, other):
118        if self.deriv > 1: self.dd -= other.dd
119        if self.deriv > 0: self.d -= other.d
120        self.v -= other.v
121        return self
122
123    def __sub__(self, other):
124        result = self.copy()
125        result -= other
126        return result
127
128    def __imul__(self, other):
129        if isinstance(other, int) or isinstance(other, float):
130            self.v *= other
131            if self.deriv > 0:
132                self.d *= other
133            if self.deriv > 1:
134                self.dd *= other
135        elif isinstance(other, Scalar):
136            # trying to avoid temporaries as much as possible
137            if self.deriv > 1:
138                self.dd *= other.v
139                self.dd += self.v*other.dd
140                tmp = np.outer(self.d, other.d)
141                self.dd += tmp
142                self.dd += tmp.transpose()
143            if self.deriv > 0:
144                self.d *= other.v
145                self.d += self.v*other.d
146            self.v *= other.v
147        else:
148            raise TypeError("Second argument must be float, int or Scalar")
149        return self
150
151    def __mul__(self, other):
152        result = self.copy()
153        result *= other
154        return result
155
156    def __itruediv__(self, other):
157        if isinstance(other, int) or isinstance(other, float):
158            self.v /= other
159            if self.deriv > 0:
160                self.d /= other
161            if self.deriv > 1:
162                self.dd /= other
163        elif isinstance(other, Scalar):
164            # trying to avoid temporaries as much as possible
165            self.v /= other.v
166            if self.deriv > 0:
167                self.d -= self.v*other.d
168                self.d /= other.v
169            if self.deriv > 1:
170                self.dd -= self.v*other.dd
171                tmp = np.outer(self.d, other.d)
172                self.dd -= tmp
173                self.dd -= tmp.transpose()
174                self.dd /= other.v
175        else:
176            raise TypeError("Second argument must be float, int or Scalar")
177        return self
178
179    def inv(self):
180        """In place invert"""
181        self.v = 1/self.v
182        tmp = self.v**2
183        if self.deriv > 1:
184            self.dd[:] = tmp*(2*self.v*np.outer(self.d, self.d) - self.dd)
185        if self.deriv > 0:
186            self.d[:] = -tmp*self.d[:]
187
188
189class Vector3(object):
190    """A Three dimensional vector with optional first and second order derivatives.
191
192       This object is nothing more than a tier for three Scalar objects.
193    """
194
195    def __init__(self, size, deriv=0, values=(0, 0, 0), indexes=(None, None, None)):
196        """
197           Arguments:
198            | ``size`` -- The number of inputs on which this ic depends. e.g. a
199                          distance depends on 6 Cartesian coordinates.
200            | ``deriv`` -- Consider up to deriv order derivatives. (max=2)
201            | ``values`` -- The initial values.
202            | ``indexes`` -- If this vector is one of the input variables, these
203                             are the indexes of the components.
204        """
205        self.deriv = deriv
206        self.size = size
207        self.x = Scalar(size, deriv, values[0], indexes[0])
208        self.y = Scalar(size, deriv, values[1], indexes[1])
209        self.z = Scalar(size, deriv, values[2], indexes[2])
210
211    def copy(self):
212        """Return a deep copy"""
213        result = Vector3(self.size, self.deriv)
214        result.x.v = self.x.v
215        result.y.v = self.y.v
216        result.z.v = self.z.v
217        if self.deriv > 0:
218            result.x.d[:] = self.x.d
219            result.y.d[:] = self.y.d
220            result.z.d[:] = self.z.d
221        if self.deriv > 1:
222            result.x.dd[:] = self.x.dd
223            result.y.dd[:] = self.y.dd
224            result.z.dd[:] = self.z.dd
225        return result
226
227    def __iadd__(self, other):
228        self.x += other.x
229        self.y += other.y
230        self.z += other.z
231        return self
232
233    def __isub__(self, other):
234        self.x -= other.x
235        self.y -= other.y
236        self.z -= other.z
237        return self
238
239    def __imul__(self, other):
240        self.x *= other
241        self.y *= other
242        self.z *= other
243        return self
244
245    def __itruediv__(self, other):
246        self.x /= other
247        self.y /= other
248        self.z /= other
249        return self
250
251    def norm(self):
252        """Return a Scalar object with the norm of this vector"""
253        result = Scalar(self.size, self.deriv)
254        result.v = np.sqrt(self.x.v**2 + self.y.v**2 + self.z.v**2)
255        if self.deriv > 0:
256            result.d += self.x.v*self.x.d
257            result.d += self.y.v*self.y.d
258            result.d += self.z.v*self.z.d
259            result.d /= result.v
260        if self.deriv > 1:
261            result.dd += self.x.v*self.x.dd
262            result.dd += self.y.v*self.y.dd
263            result.dd += self.z.v*self.z.dd
264            denom = result.v**2
265            result.dd += (1 - self.x.v**2/denom)*np.outer(self.x.d, self.x.d)
266            result.dd += (1 - self.y.v**2/denom)*np.outer(self.y.d, self.y.d)
267            result.dd += (1 - self.z.v**2/denom)*np.outer(self.z.d, self.z.d)
268            tmp = -self.x.v*self.y.v/denom*np.outer(self.x.d, self.y.d)
269            result.dd += tmp+tmp.transpose()
270            tmp = -self.y.v*self.z.v/denom*np.outer(self.y.d, self.z.d)
271            result.dd += tmp+tmp.transpose()
272            tmp = -self.z.v*self.x.v/denom*np.outer(self.z.d, self.x.d)
273            result.dd += tmp+tmp.transpose()
274            result.dd /= result.v
275        return result
276
277
278#
279# Auxiliary functions
280#
281
282
283def dot(r1, r2):
284    """Compute the dot product
285
286       Arguments:
287        | ``r1``, ``r2``  -- two :class:`Vector3` objects
288
289       (Returns a Scalar)
290    """
291    if r1.size != r2.size:
292        raise ValueError("Both arguments must have the same input size.")
293    if r1.deriv != r2.deriv:
294        raise ValueError("Both arguments must have the same deriv.")
295    return r1.x*r2.x + r1.y*r2.y + r1.z*r2.z
296
297
298def cross(r1, r2):
299    """Compute the cross product
300
301       Arguments:
302        | ``r1``, ``r2``  -- two :class:`Vector3` objects
303
304       (Returns a Vector3)
305    """
306    if r1.size != r2.size:
307        raise ValueError("Both arguments must have the same input size.")
308    if r1.deriv != r2.deriv:
309        raise ValueError("Both arguments must have the same deriv.")
310    result = Vector3(r1.size, r1.deriv)
311    result.x = r1.y*r2.z - r1.z*r2.y
312    result.y = r1.z*r2.x - r1.x*r2.z
313    result.z = r1.x*r2.y - r1.y*r2.x
314    return result
315
316
317#
318# Internal coordinate functions
319#
320
321def bond_length(rs, deriv=0):
322    """Compute the distance between the two points rs[0] and rs[1]
323
324       Arguments:
325        | ``rs``  --  two numpy array with three elements
326        | ``deriv``  --  the derivatives to be computed: 0, 1 or 2 [default=0]
327
328       When derivatives are computed a tuple with a single result is returned
329    """
330    return _bond_transform(rs, _bond_length_low, deriv)
331
332pair_distance = bond_length
333
334
335def bend_cos(rs, deriv=0):
336    """Compute the cosine of the angle between the vectors rs[0]-rs[1] and rs[2]-rs[1]
337
338       Arguments:
339        | ``rs``  --  three numpy array with three elements
340        | ``deriv``  --  the derivatives to be computed: 0, 1 or 2 [default=0]
341
342       When derivatives are computed a tuple with a single result is returned
343    """
344    return _bend_transform(rs, _bend_cos_low, deriv)
345
346
347def bend_angle(rs, deriv=0):
348    """Compute the angle between the vectors rs[0]-rs[1] and rs[2]-rs[1]
349
350       Arguments:
351        | ``rs``  --  three numpy array with three elements
352        | ``deriv``  --  the derivatives to be computed: 0, 1 or 2 [default=0]
353
354       When derivatives are computed a tuple with a single result is returned
355    """
356    return _bend_transform(rs, _bend_angle_low, deriv)
357
358
359def dihed_cos(rs, deriv=0):
360    """Compute the cosine of the angle between the planes rs[0], rs[1], rs[2] and rs[1], rs[2], rs[3]
361
362       Arguments:
363        | ``rs``  --  four numpy array with three elements
364        | ``deriv``  --  the derivatives to be computed: 0, 1 or 2 [default=0]
365    """
366    return _dihed_transform(rs, _dihed_cos_low, deriv)
367
368
369def dihed_angle(rs, deriv=0):
370    """Compute the angle between the planes rs[0], rs[1], rs[2] and rs[1], rs[2], rs[3]
371
372       The sign convention corresponds to the IUPAC definition of the torsion
373       angle: http://dx.doi.org/10.1351/goldbook.T06406
374
375       Arguments:
376        | ``rs``  --  four numpy array with three elements
377        | ``deriv``  --  the derivatives to be computed: 0, 1 or 2 [default=0]
378
379       When derivatives are computed a tuple with a single result is returned
380    """
381    return _dihed_transform(rs, _dihed_angle_low, deriv)
382
383
384def opbend_dist(rs, deriv=0):
385    """Compute the out-of-plane distance, i.e. the distance between atom rs[3] and plane rs[0],rs[1],rs[2]
386
387       Arguments:
388        | ``rs``  --  four numpy array with three elements
389        | ``deriv``  --  the derivatives to be computed: 0, 1 or 2 [default=0]
390    """
391    return _opbend_transform(rs, _opdist_low, deriv)
392
393
394def opbend_cos(rs, deriv=0):
395    """Compute the cosine of the angle between the vector (rs[0],rs[3]) and plane rs[0],rs[1],rs[2]
396
397       Arguments:
398        | ``rs``  --  four numpy array with three elements
399        | ``deriv``  --  the derivatives to be computed: 0, 1 or 2 [default=0]
400    """
401    return _opbend_transform(rs, _opbend_cos_low, deriv)
402
403
404def opbend_angle(rs, deriv=0):
405    """Compute the angle between the vector rs[0], rs[3] and the plane rs[0], rs[1], rs[2]
406
407       The sign convention is as follows: positive if rs[3] lies in the space
408       above plane rs[0], rs[1], rs[2] and negative if rs[3] lies below. Above
409       is defined by right hand rule from rs[0]-rs[1] to rs[0]-rs[2].
410
411       Arguments:
412        | ``rs``  --  four numpy array with three elements
413        | ``deriv``  --  the derivatives to be computed: 0, 1 or 2 [default=0]
414
415       When no derivatives are computed a tuple with a single result is returned.
416    """
417    return _opbend_transform(rs, _opbend_angle_low, deriv)
418
419
420def opbend_mangle(rs, deriv=0):
421    """Compute the mean value of the 3 opbend_angles
422    """
423    return _opbend_transform_mean(rs, _opbend_angle_low, deriv)
424
425
426def opbend_mcos(rs, deriv=0):
427    """Compute the mean cos of the 3 opbend_angles
428    """
429    return _opbend_transform_mean(rs, _opbend_cos_low, deriv)
430
431
432#
433# Transformers
434#
435
436
437def _bond_transform(rs, fn_low, deriv):
438    r = rs[0] - rs[1]
439    result = fn_low(r, deriv)
440    v = result[0]
441    if deriv == 0:
442        return v,
443    d = np.zeros((2, 3), float)
444    d[0] = result[1]
445    d[1] = -result[1]
446    if deriv == 1:
447        return v, d
448    dd = np.zeros((2, 3, 2, 3), float)
449    dd[0, :, 0, :] = result[2]
450    dd[1, :, 1, :] = result[2]
451    dd[0, :, 1, :] = -result[2]
452    dd[1, :, 0, :] = -result[2]
453    if deriv == 2:
454        return v, d, dd
455    raise ValueError("deriv must be 0, 1 or 2.")
456
457
458def _bend_transform(rs, fn_low, deriv):
459    a = rs[0] - rs[1]
460    b = rs[2] - rs[1]
461    result = fn_low(a, b, deriv)
462    v = result[0]
463    if deriv == 0:
464        return v,
465    d = np.zeros((3, 3), float)
466    d[0] = result[1][:3]
467    d[1] = -result[1][:3]-result[1][3:]
468    d[2] = result[1][3:]
469    if deriv == 1:
470        return v, d
471    dd = np.zeros((3, 3, 3, 3), float)
472    aa = result[2][:3, :3]
473    ab = result[2][:3, 3:]
474    ba = result[2][3:, :3]
475    bb = result[2][3:, 3:]
476    dd[0, :, 0, :] =   aa
477    dd[0, :, 1, :] = - aa - ab
478    dd[0, :, 2, :] =   ab
479    dd[1, :, 0, :] = - aa - ba
480    dd[1, :, 1, :] =   aa + ba + ab + bb
481    dd[1, :, 2, :] = - ab - bb
482    dd[2, :, 0, :] =   ba
483    dd[2, :, 1, :] = - ba - bb
484    dd[2, :, 2, :] =   bb
485    if deriv == 2:
486        return v, d, dd
487    raise ValueError("deriv must be 0, 1 or 2.")
488
489
490def _dihed_transform(rs, fn_low, deriv):
491    a = rs[0] - rs[1]
492    b = rs[2] - rs[1]
493    c = rs[3] - rs[2]
494    result = fn_low(a, b, c, deriv)
495    v = result[0]
496    if deriv == 0:
497        return v,
498    d = np.zeros((4, 3), float)
499    d[0] = result[1][:3]
500    d[1] = -result[1][:3]-result[1][3:6]
501    d[2] = result[1][3:6]-result[1][6:]
502    d[3] = result[1][6:]
503    if deriv == 1:
504        return v, d
505    dd = np.zeros((4, 3, 4, 3), float)
506    aa = result[2][:3, :3]
507    ab = result[2][:3, 3:6]
508    ac = result[2][:3, 6:]
509    ba = result[2][3:6, :3]
510    bb = result[2][3:6, 3:6]
511    bc = result[2][3:6, 6:]
512    ca = result[2][6:, :3]
513    cb = result[2][6:, 3:6]
514    cc = result[2][6:, 6:]
515
516    dd[0, :, 0, :] =   aa
517    dd[0, :, 1, :] = - aa - ab
518    dd[0, :, 2, :] =   ab - ac
519    dd[0, :, 3, :] =   ac
520
521    dd[1, :, 0, :] = - aa - ba
522    dd[1, :, 1, :] =   aa + ba + ab + bb
523    dd[1, :, 2, :] = - ab - bb + ac + bc
524    dd[1, :, 3, :] = - ac - bc
525
526    dd[2, :, 0, :] =   ba - ca
527    dd[2, :, 1, :] = - ba + ca - bb + cb
528    dd[2, :, 2, :] =   bb - cb - bc + cc
529    dd[2, :, 3, :] =   bc - cc
530
531    dd[3, :, 0, :] =   ca
532    dd[3, :, 1, :] = - ca - cb
533    dd[3, :, 2, :] =   cb - cc
534    dd[3, :, 3, :] =   cc
535    if deriv == 2:
536        return v, d, dd
537    raise ValueError("deriv must be 0, 1 or 2.")
538
539
540def _opbend_transform(rs, fn_low, deriv):
541    a = rs[1] - rs[0]
542    b = rs[2] - rs[0]
543    c = rs[3] - rs[0]
544    result = fn_low(a, b, c, deriv)
545    v = result[0]
546    if deriv == 0:
547        return v,
548    d = np.zeros((4, 3), float)
549    d[0] = -result[1][:3]-result[1][3:6]-result[1][6:]
550    d[1] = result[1][:3]
551    d[2] = result[1][3:6]
552    d[3] = result[1][6:]
553    if deriv == 1:
554        return v, d
555    dd = np.zeros((4, 3, 4, 3), float)
556    aa = result[2][:3, :3]
557    ab = result[2][:3, 3:6]
558    ac = result[2][:3, 6:]
559    ba = result[2][3:6, :3]
560    bb = result[2][3:6, 3:6]
561    bc = result[2][3:6, 6:]
562    ca = result[2][6:, :3]
563    cb = result[2][6:, 3:6]
564    cc = result[2][6:, 6:]
565
566    dd[0, :, 0, :] =   aa + ab + ac + ba + bb + bc + ca + cb + cc
567    dd[0, :, 1, :] = - aa - ba - ca
568    dd[0, :, 2, :] = - ab - bb - cb
569    dd[0, :, 3, :] = - ac - bc - cc
570
571    dd[1, :, 0, :] = - aa - ab - ac
572    dd[1, :, 1, :] =   aa
573    dd[1, :, 2, :] =   ab
574    dd[1, :, 3, :] =   ac
575
576    dd[2, :, 0, :] = - ba - bb - bc
577    dd[2, :, 1, :] =   ba
578    dd[2, :, 2, :] =   bb
579    dd[2, :, 3, :] =   bc
580
581    dd[3, :, 0, :] = - ca - cb - cc
582    dd[3, :, 1, :] =   ca
583    dd[3, :, 2, :] =   cb
584    dd[3, :, 3, :] =   cc
585    if deriv == 2:
586        return v, d, dd
587    raise ValueError("deriv must be 0, 1 or 2.")
588
589
590def _opbend_transform_mean(rs, fn_low, deriv=0):
591    """Compute the mean of the 3 opbends
592    """
593    v = 0.0
594    d = np.zeros((4,3), float)
595    dd = np.zeros((4,3,4,3), float)
596    #loop over the 3 cyclic permutations
597    for p in np.array([[0,1,2], [2,0,1], [1,2,0]]):
598        opbend = _opbend_transform([rs[p[0]], rs[p[1]], rs[p[2]], rs[3]], fn_low, deriv)
599        v += opbend[0]/3
600        index0 = np.where(p==0)[0][0] #index0 is the index of the 0th atom (rs[0])
601        index1 = np.where(p==1)[0][0]
602        index2 = np.where(p==2)[0][0]
603        index3 = 3
604        if deriv>0:
605            d[0] += opbend[1][index0]/3
606            d[1] += opbend[1][index1]/3
607            d[2] += opbend[1][index2]/3
608            d[3] += opbend[1][index3]/3
609        if deriv>1:
610            dd[0, :, 0, :] += opbend[2][index0, :, index0, :]/3
611            dd[0, :, 1, :] += opbend[2][index0, :, index1, :]/3
612            dd[0, :, 2, :] += opbend[2][index0, :, index2, :]/3
613            dd[0, :, 3, :] += opbend[2][index0, :, index3, :]/3
614
615            dd[1, :, 0, :] += opbend[2][index1, :, index0, :]/3
616            dd[1, :, 1, :] += opbend[2][index1, :, index1, :]/3
617            dd[1, :, 2, :] += opbend[2][index1, :, index2, :]/3
618            dd[1, :, 3, :] += opbend[2][index1, :, index3, :]/3
619
620            dd[2, :, 0, :] += opbend[2][index2, :, index0, :]/3
621            dd[2, :, 1, :] += opbend[2][index2, :, index1, :]/3
622            dd[2, :, 2, :] += opbend[2][index2, :, index2, :]/3
623            dd[2, :, 3, :] += opbend[2][index2, :, index3, :]/3
624
625            dd[3, :, 0, :] += opbend[2][index3, :, index0, :]/3
626            dd[3, :, 1, :] += opbend[2][index3, :, index1, :]/3
627            dd[3, :, 2, :] += opbend[2][index3, :, index2, :]/3
628            dd[3, :, 3, :] += opbend[2][index3, :, index3, :]/3
629    if deriv==0:
630        return v,
631    elif deriv==1:
632        return v, d
633    elif deriv==2:
634        return v, d, dd
635    else:
636        raise ValueError("deriv must be 0, 1 or 2.")
637
638
639#
640# Low level Internal coordinate functions
641#
642
643
644def _bond_length_low(r, deriv):
645    """Similar to bond_length, but with a relative vector"""
646    r = Vector3(3, deriv, r, (0, 1, 2))
647    d = r.norm()
648    return d.results()
649
650
651def _bend_cos_low(a, b, deriv):
652    """Similar to bend_cos, but with relative vectors"""
653    a = Vector3(6, deriv, a, (0, 1, 2))
654    b = Vector3(6, deriv, b, (3, 4, 5))
655    a /= a.norm()
656    b /= b.norm()
657    return dot(a, b).results()
658
659
660def _bend_angle_low(a, b, deriv):
661    """Similar to bend_angle, but with relative vectors"""
662    result = _bend_cos_low(a, b, deriv)
663    return _cos_to_angle(result, deriv)
664
665
666def _dihed_cos_low(a, b, c, deriv):
667    """Similar to dihed_cos, but with relative vectors"""
668    a = Vector3(9, deriv, a, (0, 1, 2))
669    b = Vector3(9, deriv, b, (3, 4, 5))
670    c = Vector3(9, deriv, c, (6, 7, 8))
671    b /= b.norm()
672    tmp = b.copy()
673    tmp *= dot(a, b)
674    a -= tmp
675    tmp = b.copy()
676    tmp *= dot(c, b)
677    c -= tmp
678    a /= a.norm()
679    c /= c.norm()
680    return dot(a, c).results()
681
682
683def _dihed_angle_low(av, bv, cv, deriv):
684    """Similar to dihed_cos, but with relative vectors"""
685    a = Vector3(9, deriv, av, (0, 1, 2))
686    b = Vector3(9, deriv, bv, (3, 4, 5))
687    c = Vector3(9, deriv, cv, (6, 7, 8))
688    b /= b.norm()
689    tmp = b.copy()
690    tmp *= dot(a, b)
691    a -= tmp
692    tmp = b.copy()
693    tmp *= dot(c, b)
694    c -= tmp
695    a /= a.norm()
696    c /= c.norm()
697    result = dot(a, c).results()
698    # avoid trobles with the gradients by either using arccos or arcsin
699    if abs(result[0]) < 0.5:
700        # if the cosine is far away for -1 or +1, it is safe to take the arccos
701        # and fix the sign of the angle.
702        sign = 1-(np.linalg.det([av, bv, cv]) > 0)*2
703        return _cos_to_angle(result, deriv, sign)
704    else:
705        # if the cosine is close to -1 or +1, it is better to compute the sine,
706        # take the arcsin and fix the sign of the angle
707        d = cross(b, a)
708        side = (result[0] > 0)*2-1 # +1 means angle in range [-pi/2,pi/2]
709        result = dot(d, c).results()
710        return _sin_to_angle(result, deriv, side)
711
712
713def _opdist_low(av, bv, cv, deriv):
714    """Similar to opdist, but with relative vectors"""
715    a = Vector3(9, deriv, av, (0, 1, 2))
716    b = Vector3(9, deriv, bv, (3, 4, 5))
717    c = Vector3(9, deriv, cv, (6, 7, 8))
718    n  = cross(a, b)
719    n /= n.norm()
720    dist = dot(c, n)
721    return dist.results()
722
723
724def _opbend_cos_low(a, b, c, deriv):
725    """Similar to opbend_cos, but with relative vectors"""
726    a = Vector3(9, deriv, a, (0, 1, 2))
727    b = Vector3(9, deriv, b, (3, 4, 5))
728    c = Vector3(9, deriv, c, (6, 7, 8))
729    n  = cross(a,b)
730    n /= n.norm()
731    c /= c.norm()
732    temp = dot(n,c)
733    result = temp.copy()
734    result.v = np.sqrt(1.0-temp.v**2)
735    if result.deriv > 0:
736        result.d *= -temp.v
737        result.d /= result.v
738    if result.deriv > 1:
739        result.dd *= -temp.v
740        result.dd /= result.v
741        temp2 = np.array([temp.d]).transpose()*temp.d
742        temp2 /= result.v**3
743        result.dd -= temp2
744    return result.results()
745
746
747def _opbend_angle_low(a, b, c, deriv=0):
748    """Similar to opbend_angle, but with relative vectors"""
749    result = _opbend_cos_low(a, b, c, deriv)
750    sign = np.sign(np.linalg.det([a, b, c]))
751    return _cos_to_angle(result, deriv, sign)
752
753
754#
755# Cosine and sine to angle conversion
756#
757
758
759def _cos_to_angle(result, deriv, sign=1):
760    """Convert a cosine and its derivatives to an angle and its derivatives"""
761    v = np.arccos(np.clip(result[0], -1, 1))
762    if deriv == 0:
763        return v*sign,
764    if abs(result[0]) >= 1:
765        factor1 = 0
766    else:
767        factor1 = -1.0/np.sqrt(1-result[0]**2)
768    d = factor1*result[1]
769    if deriv == 1:
770        return v*sign, d*sign
771    factor2 = result[0]*factor1**3
772    dd = factor2*np.outer(result[1], result[1]) + factor1*result[2]
773    if deriv == 2:
774        return v*sign, d*sign, dd*sign
775    raise ValueError("deriv must be 0, 1 or 2.")
776
777
778def _sin_to_angle(result, deriv, side=1):
779    """Convert a sine and its derivatives to an angle and its derivatives"""
780    v = np.arcsin(np.clip(result[0], -1, 1))
781    sign = side
782    if sign == -1:
783        if v < 0:
784            offset = -np.pi
785        else:
786            offset = np.pi
787    else:
788        offset = 0.0
789    if deriv == 0:
790        return v*sign + offset,
791    if abs(result[0]) >= 1:
792        factor1 = 0
793    else:
794        factor1 = 1.0/np.sqrt(1-result[0]**2)
795    d = factor1*result[1]
796    if deriv == 1:
797        return v*sign + offset, d*sign
798    factor2 = result[0]*factor1**3
799    dd = factor2*np.outer(result[1], result[1]) + factor1*result[2]
800    if deriv == 2:
801        return v*sign + offset, d*sign, dd*sign
802    raise ValueError("deriv must be 0, 1 or 2.")
803