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