1#!/usr/bin/env python 2# Copyright 2014-2020 The PySCF Developers. All Rights Reserved. 3# 4# Licensed under the Apache License, Version 2.0 (the "License"); 5# you may not use this file except in compliance with the License. 6# You may obtain a copy of the License at 7# 8# http://www.apache.org/licenses/LICENSE-2.0 9# 10# Unless required by applicable law or agreed to in writing, software 11# distributed under the License is distributed on an "AS IS" BASIS, 12# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 13# See the License for the specific language governing permissions and 14# limitations under the License. 15# 16# Author: Qiming Sun <osirpt.sun@gmail.com> 17# 18 19''' 20Mole class and helper functions to handle paramters and attributes for GTO 21integrals. This module serves the interface to the integral library libcint. 22''' 23 24import os, sys 25import types 26import re 27import platform 28import gc 29import time 30 31import json 32import ctypes 33import numpy 34import numpy as np 35import h5py 36import scipy.special 37import scipy.linalg 38import contextlib 39import threading 40from pyscf import lib 41from pyscf.lib import param 42from pyscf.data import elements 43from pyscf.lib import logger 44from pyscf.gto import cmd_args 45from pyscf.gto import basis 46from pyscf.gto import moleintor 47from pyscf.gto.eval_gto import eval_gto 48from pyscf.gto.ecp import core_configuration 49from pyscf import __config__ 50 51from pyscf.data.elements import ELEMENTS, ELEMENTS_PROTON, \ 52 _rm_digit, charge, _symbol, _std_symbol, _atom_symbol, is_ghost_atom, \ 53 _std_symbol_without_ghost 54 55from pyscf.lib.exceptions import BasisNotFoundError, PointGroupSymmetryError 56import warnings 57 58# For code compatibility in python-2 and python-3 59if sys.version_info >= (3,): 60 unicode = str 61 62 63# for _atm, _bas, _env 64CHARGE_OF = 0 65PTR_COORD = 1 66NUC_MOD_OF = 2 67PTR_ZETA = PTR_FRAC_CHARGE = 3 68ATM_SLOTS = 6 69ATOM_OF = 0 70ANG_OF = 1 71NPRIM_OF = 2 72NCTR_OF = 3 73RADI_POWER = 3 # for ECP 74KAPPA_OF = 4 75SO_TYPE_OF = 4 # for ECP 76PTR_EXP = 5 77PTR_COEFF = 6 78BAS_SLOTS = 8 79# pointer to env 80PTR_EXPCUTOFF = 0 81PTR_COMMON_ORIG = 1 82PTR_RINV_ORIG = 4 83PTR_RINV_ZETA = 7 84PTR_RANGE_OMEGA = 8 85PTR_F12_ZETA = 9 86PTR_GTG_ZETA = 10 87NGRIDS = 11 88PTR_GRIDS = 12 89AS_RINV_ORIG_ATOM = 17 90AS_ECPBAS_OFFSET = 18 91AS_NECPBAS = 19 92PTR_ENV_START = 20 93# parameters from libcint 94NUC_POINT = 1 95NUC_GAUSS = 2 96# nucleus with fractional charges. It can be used to mimic MM particles 97NUC_FRAC_CHARGE = 3 98NUC_ECP = 4 # atoms with pseudo potential 99 100BASE = getattr(__config__, 'BASE', 0) 101NORMALIZE_GTO = getattr(__config__, 'NORMALIZE_GTO', True) 102DISABLE_EVAL = getattr(__config__, 'DISABLE_EVAL', False) 103# Whether to disable the explicit call to gc.collect(). gc.collect() may cause 104# non-negligible overhead (https://github.com/pyscf/pyscf/issues/1038). 105DISABLE_GC = getattr(__config__, 'DISABLE_GC', False) 106 107def M(**kwargs): 108 r'''This is a shortcut to build up Mole object. 109 110 Args: Same to :func:`Mole.build` 111 112 Examples: 113 114 >>> from pyscf import gto 115 >>> mol = gto.M(atom='H 0 0 0; F 0 0 1', basis='6-31g') 116 ''' 117 mol = Mole() 118 mol.build(**kwargs) 119 return mol 120 121def gaussian_int(n, alpha): 122 r'''int_0^inf x^n exp(-alpha x^2) dx''' 123 n1 = (n + 1) * .5 124 return scipy.special.gamma(n1) / (2. * alpha**n1) 125 126def gto_norm(l, expnt): 127 r'''Normalized factor for GTO radial part :math:`g=r^l e^{-\alpha r^2}` 128 129 .. math:: 130 131 \frac{1}{\sqrt{\int g^2 r^2 dr}} 132 = \sqrt{\frac{2^{2l+3} (l+1)! (2a)^{l+1.5}}{(2l+2)!\sqrt{\pi}}} 133 134 Ref: H. B. Schlegel and M. J. Frisch, Int. J. Quant. Chem., 54(1995), 83-87. 135 136 Args: 137 l (int): 138 angular momentum 139 expnt : 140 exponent :math:`\alpha` 141 142 Returns: 143 normalization factor 144 145 Examples: 146 147 >>> print(gto_norm(0, 1)) 148 2.5264751109842591 149 ''' 150 if l >= 0: 151 #f = 2**(2*l+3) * math.factorial(l+1) * (2*expnt)**(l+1.5) \ 152 # / (math.factorial(2*l+2) * math.sqrt(math.pi)) 153 #return math.sqrt(f) 154 return 1/numpy.sqrt(gaussian_int(l*2+2, 2*expnt)) 155 else: 156 raise ValueError('l should be > 0') 157 158def cart2sph(l, c_tensor=None, normalized=None): 159 ''' 160 Cartesian to real spherical transformation matrix 161 162 Kwargs: 163 normalized : 164 How the Cartesian GTOs are normalized. 'sp' means the s and p 165 functions are normalized (this is the convention used by libcint 166 library). 167 ''' 168 nf = (l+1)*(l+2)//2 169 if c_tensor is None: 170 c_tensor = numpy.eye(nf) 171 else: 172 c_tensor = numpy.asarray(c_tensor, order='F').reshape(-1,nf) 173 if l == 0 or l == 1: 174 if normalized == 'sp': 175 return c_tensor 176 elif l == 0: 177 return c_tensor * 0.282094791773878143 178 else: 179 return c_tensor * 0.488602511902919921 180 else: 181 assert(l <= 12) 182 nd = l * 2 + 1 183 ngrid = c_tensor.shape[0] 184 c2sph = numpy.zeros((ngrid,nd), order='F') 185 fn = moleintor.libcgto.CINTc2s_ket_sph 186 fn(c2sph.ctypes.data_as(ctypes.c_void_p), ctypes.c_int(ngrid), 187 c_tensor.ctypes.data_as(ctypes.c_void_p), ctypes.c_int(l)) 188 return c2sph 189 190def cart2spinor_kappa(kappa, l=None, normalized=None): 191 '''Cartesian to spinor transformation matrix for kappa 192 193 Kwargs: 194 normalized : 195 How the Cartesian GTOs are normalized. 'sp' means the s and p 196 functions are normalized (this is the convention used by libcint 197 library). 198 ''' 199 if kappa < 0: 200 l = -kappa - 1 201 nd = l * 2 + 2 202 elif kappa > 0: 203 l = kappa 204 nd = l * 2 205 else: 206 assert(l is not None) 207 assert(l <= 12) 208 nd = l * 4 + 2 209 nf = (l+1)*(l+2)//2 210 c2smat = numpy.zeros((nf*2,nd), order='F', dtype=numpy.complex128) 211 cmat = numpy.eye(nf) 212 fn = moleintor.libcgto.CINTc2s_ket_spinor_sf1 213 fn(c2smat.ctypes.data_as(ctypes.c_void_p), 214 c2smat[nf:].ctypes.data_as(ctypes.c_void_p), 215 cmat.ctypes.data_as(ctypes.c_void_p), 216 ctypes.c_int(nf*2), ctypes.c_int(nf), 217 ctypes.c_int(1), ctypes.c_int(kappa), ctypes.c_int(l)) 218 if normalized != 'sp': 219 if l == 0: 220 c2smat *= 0.282094791773878143 221 elif l == 1: 222 c2smat *= 0.488602511902919921 223 # c2smat[0] is the transformation for spin up 224 # c2smat[1] is the transformation for spin down 225 c2smat = c2smat.reshape(2,nf,nd) 226 return c2smat 227cart2j_kappa = cart2spinor_kappa 228 229def cart2spinor_l(l, normalized=None): 230 '''Cartesian to spinor transformation matrix for angular moment l 231 232 Kwargs: 233 normalized : 234 How the Cartesian GTOs are normalized. 'sp' means the s and p 235 functions are normalized (this is the convention used by libcint 236 library). 237 ''' 238 return cart2spinor_kappa(0, l, normalized) 239cart2j_l = cart2spinor_l 240 241def sph2spinor_kappa(kappa, l=None): 242 '''Real spherical to spinor transformation matrix for kappa''' 243 from pyscf.symm.sph import sph2spinor 244 ua, ub = sph2spinor(l) 245 if kappa < 0: 246 l = -kappa - 1 247 ua = ua[:,l*2:] 248 ub = ub[:,l*2:] 249 elif kappa > 0: 250 l = kappa 251 ua = ua[:,:l*2] 252 ub = ub[:,:l*2] 253 else: 254 assert(l is not None) 255 assert(l <= 12) 256 return ua, ub 257 258def sph2spinor_l(l): 259 '''Real spherical to spinor transformation matrix for angular moment l''' 260 return sph2spinor_kappa(0, l) 261 262def atom_types(atoms, basis=None): 263 '''symmetry inequivalent atoms''' 264 atmgroup = {} 265 for ia, a in enumerate(atoms): 266 if 'GHOST' in a[0].upper(): 267 a = ['X'+a[0][5:]] + list(a[1:]) 268 if a[0] in atmgroup: 269 atmgroup[a[0]].append(ia) 270 elif basis is None: 271 atmgroup[a[0]] = [ia] 272 else: 273 stdsymb = _std_symbol(a[0]) 274 if a[0] in basis: 275 if stdsymb in basis and basis[a[0]] == basis[stdsymb]: 276 if stdsymb in atmgroup: 277 atmgroup[stdsymb].append(ia) 278 else: 279 atmgroup[stdsymb] = [ia] 280 else: 281 atmgroup[a[0]] = [ia] 282 elif stdsymb in atmgroup: 283 atmgroup[stdsymb].append(ia) 284 else: 285 atmgroup[stdsymb] = [ia] 286 return atmgroup 287 288 289def format_atom(atoms, origin=0, axes=None, 290 unit=getattr(__config__, 'UNIT', 'Ang')): 291 '''Convert the input :attr:`Mole.atom` to the internal data format. 292 Including, changing the nuclear charge to atom symbol, converting the 293 coordinates to AU, rotate and shift the molecule. 294 If the :attr:`~Mole.atom` is a string, it takes ";" and "\\n" 295 for the mark to separate atoms; "," and arbitrary length of blank space 296 to spearate the individual terms for an atom. Blank line will be ignored. 297 298 Args: 299 atoms : list or str 300 the same to :attr:`Mole.atom` 301 302 Kwargs: 303 origin : ndarray 304 new axis origin. 305 axes : ndarray 306 (new_x, new_y, new_z), new coordinates 307 unit : str or number 308 If unit is one of strings (B, b, Bohr, bohr, AU, au), the coordinates 309 of the input atoms are the atomic unit; If unit is one of strings 310 (A, a, Angstrom, angstrom, Ang, ang), the coordinates are in the 311 unit of angstrom; If a number is given, the number are considered 312 as the Bohr value (in angstrom), which should be around 0.53. 313 Set unit=1 if wishing to preserve the unit of the coordinates. 314 315 Returns: 316 "atoms" in the internal format. The internal format is 317 | atom = [[atom1, (x, y, z)], 318 | [atom2, (x, y, z)], 319 | ... 320 | [atomN, (x, y, z)]] 321 322 Examples: 323 324 >>> gto.format_atom('9,0,0,0; h@1 0 0 1', origin=(1,1,1)) 325 [['F', [-1.0, -1.0, -1.0]], ['H@1', [-1.0, -1.0, 0.0]]] 326 >>> gto.format_atom(['9,0,0,0', (1, (0, 0, 1))], origin=(1,1,1)) 327 [['F', [-1.0, -1.0, -1.0]], ['H', [-1, -1, 0]]] 328 ''' 329 def str2atm(line): 330 dat = line.split() 331 try: 332 coords = [float(x) for x in dat[1:4]] 333 except ValueError: 334 if DISABLE_EVAL: 335 raise ValueError('Failed to parse geometry %s' % line) 336 else: 337 coords = list(eval(','.join(dat[1:4]))) 338 if len(coords) != 3: 339 raise ValueError('Coordinates error in %s' % line) 340 return [_atom_symbol(dat[0]), coords] 341 342 if isinstance(atoms, (str, unicode)): 343 # The input atoms points to a geometry file 344 if os.path.isfile(atoms): 345 try: 346 atoms = fromfile(atoms) 347 except ValueError: 348 sys.stderr.write('\nFailed to parse geometry file %s\n\n' % atoms) 349 raise 350 351 atoms = str(atoms.replace(';','\n').replace(',',' ').replace('\t',' ')) 352 fmt_atoms = [] 353 for dat in atoms.split('\n'): 354 dat = dat.strip() 355 if dat and dat[0] != '#': 356 fmt_atoms.append(dat) 357 358 if len(fmt_atoms[0].split()) < 4: 359 fmt_atoms = from_zmatrix('\n'.join(fmt_atoms)) 360 else: 361 fmt_atoms = [str2atm(line) for line in fmt_atoms] 362 else: 363 fmt_atoms = [] 364 for atom in atoms: 365 if isinstance(atom, (str, unicode)): 366 if atom.lstrip()[0] != '#': 367 fmt_atoms.append(str2atm(atom.replace(',',' '))) 368 else: 369 if isinstance(atom[1], (int, float)): 370 fmt_atoms.append([_atom_symbol(atom[0]), atom[1:4]]) 371 else: 372 fmt_atoms.append([_atom_symbol(atom[0]), atom[1]]) 373 374 if len(fmt_atoms) == 0: 375 return [] 376 377 if axes is None: 378 axes = numpy.eye(3) 379 380 if isinstance(unit, (str, unicode)): 381 if unit.upper().startswith(('B', 'AU')): 382 unit = 1. 383 else: #unit[:3].upper() == 'ANG': 384 unit = 1./param.BOHR 385 else: 386 unit = 1./unit 387 388 c = numpy.array([a[1] for a in fmt_atoms], dtype=numpy.double) 389 c = numpy.einsum('ix,kx->ki', axes * unit, c - origin) 390 z = [a[0] for a in fmt_atoms] 391 return list(zip(z, c.tolist())) 392 393#TODO: sort exponents 394def format_basis(basis_tab): 395 '''Convert the input :attr:`Mole.basis` to the internal data format. 396 397 ``{ atom: [(l, ((-exp, c_1, c_2, ..), 398 (-exp, c_1, c_2, ..))), 399 (l, ((-exp, c_1, c_2, ..), 400 (-exp, c_1, c_2, ..)))], ... }`` 401 402 Args: 403 basis_tab : dict 404 Similar to :attr:`Mole.basis`, it **cannot** be a str 405 406 Returns: 407 Formated :attr:`~Mole.basis` 408 409 Examples: 410 411 >>> gto.format_basis({'H':'sto-3g', 'H^2': '3-21g'}) 412 {'H': [[0, 413 [3.4252509099999999, 0.15432897000000001], 414 [0.62391373000000006, 0.53532813999999995], 415 [0.16885539999999999, 0.44463454000000002]]], 416 'H^2': [[0, 417 [5.4471780000000001, 0.15628500000000001], 418 [0.82454700000000003, 0.90469100000000002]], 419 [0, [0.18319199999999999, 1.0]]]} 420 ''' 421 basis_converter = _generate_basis_converter() 422 fmt_basis = {} 423 for atom, atom_basis in basis_tab.items(): 424 symb = _atom_symbol(atom) 425 fmt_basis[symb] = basis_converter(symb, atom_basis) 426 427 if len(fmt_basis[symb]) == 0: 428 raise BasisNotFoundError('Basis not found for %s' % symb) 429 return fmt_basis 430 431def _generate_basis_converter(): 432 def nparray_to_list(item): 433 val = [] 434 for x in item: 435 if isinstance(x, (tuple, list)): 436 val.append(nparray_to_list(x)) 437 elif isinstance(x, numpy.ndarray): 438 val.append(x.tolist()) 439 else: 440 val.append(x) 441 return val 442 443 def load(basis_name, symb): 444 if basis_name.lower().startswith('unc'): 445 return uncontract(basis.load(basis_name[3:], symb)) 446 else: 447 return basis.load(basis_name, symb) 448 449 def converter(symb, raw_basis): 450 if isinstance(raw_basis, (str, unicode)): 451 bset = load(str(raw_basis), _std_symbol_without_ghost(symb)) 452 elif (any(isinstance(x, (str, unicode)) for x in raw_basis) 453 # The first element is the basis of internal format 454 or not isinstance(raw_basis[0][0], int)): 455 stdsymb = _std_symbol_without_ghost(symb) 456 bset = [] 457 for rawb in raw_basis: 458 if isinstance(rawb, (str, unicode)): 459 bset += load(str(rawb), stdsymb) 460 else: 461 bset += nparray_to_list(rawb) 462 else: 463 bset = nparray_to_list(raw_basis) 464 return bset 465 return converter 466 467def uncontracted_basis(_basis): 468 '''Uncontract internal format _basis 469 470 Examples: 471 472 >>> gto.uncontract(gto.load('sto3g', 'He')) 473 [[0, [6.36242139, 1]], [0, [1.158923, 1]], [0, [0.31364979, 1]]] 474 ''' 475 MAXL = 10 476 ubasis_raw = [[] for l in range(MAXL)] 477 ubasis_exp = [[] for l in range(MAXL)] 478 for b in _basis: 479 angl = b[0] 480 kappa = b[1] 481 if isinstance(kappa, int): 482 coeffs = b[2:] 483 else: 484 coeffs = b[1:] 485 486 if isinstance(kappa, int) and kappa != 0: 487 warnings.warn('For basis with kappa != 0, the uncontract basis might be wrong. ' 488 'Please double check the resultant attribute mol._basis') 489 for p in coeffs: 490 ubasis_raw[angl].append([angl, kappa, [p[0], 1]]) 491 ubasis_exp[angl].append(p[0]) 492 else: 493 for p in coeffs: 494 ubasis_raw[angl].append([angl, [p[0], 1]]) 495 ubasis_exp[angl].append(p[0]) 496 497 # Check linear dependency 498 ubasis = [] 499 for l in range(MAXL): 500 basis_l = ubasis_raw[l] 501 if basis_l: 502 es = numpy.array(ubasis_exp[l]) 503 # Remove duplicated primitive basis functions 504 es, e_idx = numpy.unique(es.round(9), True) 505 # from large exponent to small exponent 506 for i in reversed(e_idx): 507 ubasis.append(basis_l[i]) 508 return ubasis 509uncontract = uncontracted_basis 510contract = contracted_basis = basis.to_general_contraction 511 512def to_uncontracted_cartesian_basis(mol): 513 '''Decontract the basis of a Mole or a Cell. Returns a Mole (Cell) object 514 with the uncontracted basis environment and a list of coefficients that 515 transform the uncontracted cartesian basis to the original basis. Each 516 element in the list corresponds to one shell of the original Mole (Cell). 517 518 Examples: 519 520 >>> mol = gto.M(atom='Ne', basis='ccpvdz') 521 >>> pmol, ctr_coeff = mol.to_uncontracted_cartesian_basis() 522 >>> c = scipy.linalg.block_diag(*ctr_coeff) 523 >>> s = reduce(numpy.dot, (c.T, pmol.intor('int1e_ovlp'), c)) 524 >>> abs(s-mol.intor('int1e_ovlp')).max() 525 0.0 526 ''' 527 import copy 528 lmax = mol._bas[:,ANG_OF].max() 529 if mol.cart: 530 c2s = [numpy.eye((l+1)*(l+2)//2) for l in range(lmax+1)] 531 else: 532 c2s = [cart2sph(l, normalized='sp') for l in range(lmax+1)] 533 534 pmol = copy.copy(mol) 535 pmol.cart = True 536 _bas = [] 537 _env = mol._env.copy() 538 contr_coeff = [] 539 for ib in range(mol.nbas): 540 l = mol._bas[ib,ANG_OF] 541 ncart = (l+1)*(l+2)//2 542 es = mol.bas_exp(ib) 543 cs = mol._libcint_ctr_coeff(ib) 544 nprim = cs.shape[0] 545 norm = gto_norm(l, es) 546 c = numpy.einsum('pi,p,xm->pxim', cs, 1./norm, c2s[l]) 547 contr_coeff.append(c.reshape(nprim*ncart,-1)) 548 549 pexp = mol._bas[ib,PTR_EXP] 550 pc = mol._bas[ib,PTR_COEFF] 551 bs = numpy.empty((nprim,8), dtype=numpy.int32) 552 bs[:] = mol._bas[ib] 553 bs[:,NCTR_OF] = bs[:,NPRIM_OF] = 1 554 bs[:,PTR_EXP] = numpy.arange(pexp, pexp+nprim) 555 bs[:,PTR_COEFF] = numpy.arange(pc, pc+nprim) 556 _env[pc:pc+nprim] = norm 557 _bas.append(bs) 558 559 pmol._bas = numpy.asarray(numpy.vstack(_bas), dtype=numpy.int32) 560 pmol._env = _env 561 return pmol, contr_coeff 562 563def format_ecp(ecp_tab): 564 r'''Convert the input :attr:`ecp` (dict) to the internal data format:: 565 566 { atom: (nelec, # core electrons 567 ((l, # l=-1 for UL, l>=0 for Ul to indicate |l><l| 568 (((exp_1, c_1), # for r^0 569 (exp_2, c_2), 570 ...), 571 ((exp_1, c_1), # for r^1 572 (exp_2, c_2), 573 ...), 574 ((exp_1, c_1), # for r^2 575 ...))))), 576 ...} 577 ''' 578 fmt_ecp = {} 579 for atom, atom_ecp in ecp_tab.items(): 580 symb = _atom_symbol(atom) 581 582 if isinstance(atom_ecp, (str, unicode)): 583 stdsymb = _std_symbol_without_ghost(symb) 584 ecp_dat = basis.load_ecp(str(atom_ecp), stdsymb) 585 if ecp_dat is None or len(ecp_dat) == 0: 586 #raise BasisNotFoundError('ECP not found for %s' % symb) 587 sys.stderr.write('ECP %s not found for %s\n' % 588 (atom_ecp, symb)) 589 else: 590 fmt_ecp[symb] = ecp_dat 591 else: 592 fmt_ecp[symb] = atom_ecp 593 return fmt_ecp 594 595# transform etb to basis format 596def expand_etb(l, n, alpha, beta): 597 r'''Generate the exponents of even tempered basis for :attr:`Mole.basis`. 598 .. math:: 599 600 e = e^{-\alpha * \beta^{i-1}} for i = 1 .. n 601 602 Args: 603 l : int 604 Angular momentum 605 n : int 606 Number of GTOs 607 608 Returns: 609 Formated :attr:`~Mole.basis` 610 611 Examples: 612 613 >>> gto.expand_etb(1, 3, 1.5, 2) 614 [[1, [6.0, 1]], [1, [3.0, 1]], [1, [1.5, 1]]] 615 ''' 616 return [[l, [alpha*beta**i, 1]] for i in reversed(range(n))] 617def expand_etbs(etbs): 618 r'''Generate even tempered basis. See also :func:`expand_etb` 619 620 Args: 621 etbs = [(l, n, alpha, beta), (l, n, alpha, beta),...] 622 623 Returns: 624 Formated :attr:`~Mole.basis` 625 626 Examples: 627 628 >>> gto.expand_etbs([(0, 2, 1.5, 2.), (1, 2, 1, 2.)]) 629 [[0, [6.0, 1]], [0, [3.0, 1]], [1, [1., 1]], [1, [2., 1]]] 630 ''' 631 return lib.flatten([expand_etb(*etb) for etb in etbs]) 632etbs = expand_etbs 633 634# concatenate two mol 635def conc_env(atm1, bas1, env1, atm2, bas2, env2): 636 r'''Concatenate two Mole's integral parameters. This function can be used 637 to construct the environment for cross integrals like 638 639 .. math:: 640 641 \langle \mu | \nu \rangle, \mu \in mol1, \nu \in mol2 642 643 Returns: 644 Concatenated atm, bas, env 645 646 Examples: 647 Compute the overlap between H2 molecule and O atom 648 649 >>> mol1 = gto.M(atom='H 0 1 0; H 0 0 1', basis='sto3g') 650 >>> mol2 = gto.M(atom='O 0 0 0', basis='sto3g') 651 >>> atm3, bas3, env3 = gto.conc_env(mol1._atm, mol1._bas, mol1._env, 652 ... mol2._atm, mol2._bas, mol2._env) 653 >>> gto.moleintor.getints('int1e_ovlp_sph', atm3, bas3, env3, range(2), range(2,5)) 654 [[ 0.04875181 0.44714688 0. 0.37820346 0. ] 655 [ 0.04875181 0.44714688 0. 0. 0.37820346]] 656 ''' 657 off = len(env1) 658 natm_off = len(atm1) 659 atm2 = numpy.copy(atm2) 660 bas2 = numpy.copy(bas2) 661 atm2[:,PTR_COORD] += off 662 atm2[:,PTR_ZETA ] += off 663 bas2[:,ATOM_OF ] += natm_off 664 bas2[:,PTR_EXP ] += off 665 bas2[:,PTR_COEFF] += off 666 return (numpy.asarray(numpy.vstack((atm1,atm2)), dtype=numpy.int32), 667 numpy.asarray(numpy.vstack((bas1,bas2)), dtype=numpy.int32), 668 numpy.hstack((env1,env2))) 669 670def conc_mol(mol1, mol2): 671 '''Concatenate two Mole objects. 672 ''' 673 if not mol1._built: 674 logger.warn(mol1, 'Warning: object %s not initialized. Initializing %s', 675 mol1, mol1) 676 mol1.build() 677 if not mol2._built: 678 logger.warn(mol2, 'Warning: object %s not initialized. Initializing %s', 679 mol2, mol2) 680 mol2.build() 681 682 mol3 = Mole() 683 mol3._built = True 684 685 mol3._atm, mol3._bas, mol3._env = \ 686 conc_env(mol1._atm, mol1._bas, mol1._env, 687 mol2._atm, mol2._bas, mol2._env) 688 off = len(mol1._env) 689 natm_off = len(mol1._atm) 690 if len(mol2._ecpbas) == 0: 691 mol3._ecpbas = mol1._ecpbas 692 else: 693 ecpbas2 = numpy.copy(mol2._ecpbas) 694 ecpbas2[:,ATOM_OF ] += natm_off 695 ecpbas2[:,PTR_EXP ] += off 696 ecpbas2[:,PTR_COEFF] += off 697 if len(mol1._ecpbas) == 0: 698 mol3._ecpbas = ecpbas2 699 else: 700 mol3._ecpbas = numpy.hstack((mol1._ecpbas, ecpbas2)) 701 702 mol3.verbose = mol1.verbose 703 mol3.output = mol1.output 704 mol3.max_memory = mol1.max_memory 705 mol3.charge = mol1.charge + mol2.charge 706 mol3.spin = abs(mol1.spin - mol2.spin) 707 mol3.symmetry = False 708 mol3.symmetry_subgroup = None 709 mol3.cart = mol1.cart and mol2.cart 710 711 mol3._atom = mol1._atom + mol2._atom 712 mol3.atom = mol3._atom 713 mol3.unit = 'Bohr' 714 715 mol3._basis.update(mol2._basis) 716 mol3._basis.update(mol1._basis) 717 mol3._pseudo.update(mol2._pseudo) 718 mol3._pseudo.update(mol1._pseudo) 719 mol3._ecp.update(mol2._ecp) 720 mol3._ecp.update(mol1._ecp) 721 mol3.basis = mol3._basis 722 mol3.ecp = mol3._ecp 723 724 mol3.nucprop.update(mol1.nucprop) 725 mol3.nucprop.update(mol2.nucprop) 726 return mol3 727 728# <bas-of-mol1|intor|bas-of-mol2> 729def intor_cross(intor, mol1, mol2, comp=None, grids=None): 730 r'''1-electron integrals from two molecules like 731 732 .. math:: 733 734 \langle \mu | intor | \nu \rangle, \mu \in mol1, \nu \in mol2 735 736 Args: 737 intor : str 738 Name of the 1-electron integral, such as int1e_ovlp_sph (spherical overlap), 739 int1e_nuc_cart (cartesian nuclear attraction), int1e_ipovlp_spinor 740 (spinor overlap gradients), etc. Ref to :func:`getints` for the 741 full list of available 1-electron integral names 742 mol1, mol2: 743 :class:`Mole` objects 744 745 Kwargs: 746 comp : int 747 Components of the integrals, e.g. int1e_ipovlp_sph has 3 components 748 grids : ndarray 749 Coordinates of grids for the int1e_grids integrals 750 751 Returns: 752 ndarray of 1-electron integrals, can be either 2-dim or 3-dim, depending on comp 753 754 Examples: 755 Compute the overlap between H2 molecule and O atom 756 757 >>> mol1 = gto.M(atom='H 0 1 0; H 0 0 1', basis='sto3g') 758 >>> mol2 = gto.M(atom='O 0 0 0', basis='sto3g') 759 >>> gto.intor_cross('int1e_ovlp_sph', mol1, mol2) 760 [[ 0.04875181 0.44714688 0. 0.37820346 0. ] 761 [ 0.04875181 0.44714688 0. 0. 0.37820346]] 762 ''' 763 nbas1 = len(mol1._bas) 764 nbas2 = len(mol2._bas) 765 atmc, basc, envc = conc_env(mol1._atm, mol1._bas, mol1._env, 766 mol2._atm, mol2._bas, mol2._env) 767 if '_grids' in intor: 768 assert grids is not None 769 envc = numpy.append(envc, grids.ravel()) 770 envc[NGRIDS] = grids.shape[0] 771 envc[PTR_GRIDS] = envc.size - grids.size 772 773 shls_slice = (0, nbas1, nbas1, nbas1+nbas2) 774 775 if (intor.endswith('_sph') or intor.startswith('cint') or 776 intor.endswith('_spinor') or intor.endswith('_cart')): 777 return moleintor.getints(intor, atmc, basc, envc, shls_slice, comp, 0) 778 elif mol1.cart == mol2.cart: 779 intor = mol1._add_suffix(intor) 780 return moleintor.getints(intor, atmc, basc, envc, shls_slice, comp, 0) 781 elif mol1.cart: 782 mat = moleintor.getints(intor+'_cart', atmc, basc, envc, shls_slice, comp, 0) 783 return numpy.dot(mat, mol2.cart2sph_coeff()) 784 else: 785 mat = moleintor.getints(intor+'_cart', atmc, basc, envc, shls_slice, comp, 0) 786 return numpy.dot(mol1.cart2sph_coeff().T, mat) 787 788# append (charge, pointer to coordinates, nuc_mod) to _atm 789def make_atm_env(atom, ptr=0, nuclear_model=NUC_POINT, nucprop={}): 790 '''Convert the internal format :attr:`Mole._atom` to the format required 791 by ``libcint`` integrals 792 ''' 793 nuc_charge = charge(atom[0]) 794 if nuclear_model == NUC_POINT: 795 zeta = 0 796 elif nuclear_model == NUC_GAUSS: 797 zeta = dyall_nuc_mod(nuc_charge, nucprop) 798 else: # callable(nuclear_model) 799 zeta = nuclear_model(nuc_charge, nucprop) 800 nuclear_model = NUC_GAUSS 801 802 _env = numpy.hstack((atom[1], zeta)) 803 _atm = numpy.zeros(6, dtype=numpy.int32) 804 _atm[CHARGE_OF] = nuc_charge 805 _atm[PTR_COORD] = ptr 806 _atm[NUC_MOD_OF] = nuclear_model 807 _atm[PTR_ZETA ] = ptr + 3 808 return _atm, _env 809 810# append (atom, l, nprim, nctr, kappa, ptr_exp, ptr_coeff, 0) to bas 811# absorb normalization into GTO contraction coefficients 812def make_bas_env(basis_add, atom_id=0, ptr=0): 813 '''Convert :attr:`Mole.basis` to the argument ``bas`` for ``libcint`` integrals 814 ''' 815 _bas = [] 816 _env = [] 817 for b in basis_add: 818 if not b: # == [] 819 continue 820 angl = b[0] 821 if angl > 14: 822 sys.stderr.write('Warning: integral library does not support basis ' 823 'with angular momentum > 14\n') 824 825 if isinstance(b[1], int): 826 kappa = b[1] 827 b_coeff = numpy.array(sorted(list(b[2:]), reverse=True)) 828 else: 829 kappa = 0 830 b_coeff = numpy.array(sorted(list(b[1:]), reverse=True)) 831 es = b_coeff[:,0] 832 cs = b_coeff[:,1:] 833 nprim, nctr = cs.shape 834 cs = numpy.einsum('pi,p->pi', cs, gto_norm(angl, es)) 835 if NORMALIZE_GTO: 836 cs = _nomalize_contracted_ao(angl, es, cs) 837 838 _env.append(es) 839 _env.append(cs.T.reshape(-1)) 840 ptr_exp = ptr 841 ptr_coeff = ptr_exp + nprim 842 ptr = ptr_coeff + nprim * nctr 843 _bas.append([atom_id, angl, nprim, nctr, kappa, ptr_exp, ptr_coeff, 0]) 844 _env = lib.flatten(_env) # flatten nested lists 845 return (numpy.array(_bas, numpy.int32).reshape(-1,BAS_SLOTS), 846 numpy.array(_env, numpy.double)) 847 848def _nomalize_contracted_ao(l, es, cs): 849 #ee = numpy.empty((nprim,nprim)) 850 #for i in range(nprim): 851 # for j in range(i+1): 852 # ee[i,j] = ee[j,i] = gaussian_int(angl*2+2, es[i]+es[j]) 853 #s1 = 1/numpy.sqrt(numpy.einsum('pi,pq,qi->i', cs, ee, cs)) 854 ee = es.reshape(-1,1) + es.reshape(1,-1) 855 ee = gaussian_int(l*2+2, ee) 856 s1 = 1. / numpy.sqrt(numpy.einsum('pi,pq,qi->i', cs, ee, cs)) 857 return numpy.einsum('pi,i->pi', cs, s1) 858 859def make_env(atoms, basis, pre_env=[], nucmod={}, nucprop={}): 860 '''Generate the input arguments for ``libcint`` library based on internal 861 format :attr:`Mole._atom` and :attr:`Mole._basis` 862 ''' 863 _atm = [] 864 _bas = [] 865 _env = [] 866 ptr_env = len(pre_env) 867 868 for ia, atom in enumerate(atoms): 869 symb = atom[0] 870 stdsymb = _rm_digit(symb) 871 872 if ia+1 in nucprop: 873 prop = nucprop[ia+1] 874 elif symb in nucprop: 875 prop = nucprop[symb] 876 else: 877 prop = nucprop.get(stdsymb, {}) 878 879 nuclear_model = NUC_POINT 880 if nucmod: 881 if nucmod is None: 882 nuclear_model = NUC_POINT 883 elif isinstance(nucmod, (int, str, unicode, types.FunctionType)): 884 nuclear_model = _parse_nuc_mod(nucmod) 885 elif ia+1 in nucmod: 886 nuclear_model = _parse_nuc_mod(nucmod[ia+1]) 887 elif symb in nucmod: 888 nuclear_model = _parse_nuc_mod(nucmod[symb]) 889 elif stdsymb in nucmod: 890 nuclear_model = _parse_nuc_mod(nucmod[stdsymb]) 891 atm0, env0 = make_atm_env(atom, ptr_env, nuclear_model, prop) 892 ptr_env = ptr_env + len(env0) 893 _atm.append(atm0) 894 _env.append(env0) 895 896 _basdic = {} 897 for symb, basis_add in basis.items(): 898 bas0, env0 = make_bas_env(basis_add, 0, ptr_env) 899 ptr_env = ptr_env + len(env0) 900 _basdic[symb] = bas0 901 _env.append(env0) 902 903 for ia, atom in enumerate(atoms): 904 symb = atom[0] 905 puresymb = _rm_digit(symb) 906 if symb in _basdic: 907 b = _basdic[symb].copy() 908 elif puresymb in _basdic: 909 b = _basdic[puresymb].copy() 910 else: 911 if symb[:2].upper() == 'X-': 912 symb = symb[2:] 913 elif symb[:6].upper() == 'GHOST-': 914 symb = symb[6:] 915 puresymb = _rm_digit(symb) 916 if symb in _basdic: 917 b = _basdic[symb].copy() 918 elif puresymb in _basdic: 919 b = _basdic[puresymb].copy() 920 else: 921 sys.stderr.write('Warning: Basis not found for atom %d %s\n' % (ia, symb)) 922 continue 923 b[:,ATOM_OF] = ia 924 _bas.append(b) 925 926 if _atm: 927 _atm = numpy.asarray(numpy.vstack(_atm), numpy.int32).reshape(-1, ATM_SLOTS) 928 else: 929 _atm = numpy.zeros((0,ATM_SLOTS), numpy.int32) 930 if _bas: 931 _bas = numpy.asarray(numpy.vstack(_bas), numpy.int32).reshape(-1, BAS_SLOTS) 932 else: 933 _bas = numpy.zeros((0,BAS_SLOTS), numpy.int32) 934 if _env: 935 _env = numpy.hstack((pre_env,numpy.hstack(_env))) 936 else: 937 _env = numpy.array(pre_env, copy=False) 938 return _atm, _bas, _env 939 940def make_ecp_env(mol, _atm, ecp, pre_env=[]): 941 '''Generate the input arguments _ecpbas for ECP integrals 942 ''' 943 _env = [] 944 ptr_env = len(pre_env) 945 946 _ecpdic = {} 947 for symb, ecp_add in ecp.items(): 948 ecp0 = [] 949 nelec = ecp_add[0] 950 for lb in ecp_add[1]: 951 for rorder, bi in enumerate(lb[1]): 952 if len(bi) > 0: 953 ec = numpy.array(sorted(bi, reverse=True)) 954 nexp, ncol = ec.shape 955 _env.append(ec[:,0]) 956 _env.append(ec[:,1]) 957 ptr_exp, ptr_coeff = ptr_env, ptr_env + nexp 958 ecp0.append([0, lb[0], nexp, rorder, 0, 959 ptr_exp, ptr_coeff, 0]) 960 ptr_env += nexp * 2 961 962 if ncol == 3: # Has SO-ECP 963 _env.append(ec[:,2]) 964 ptr_coeff, ptr_env = ptr_env, ptr_env + nexp 965 ecp0.append([0, lb[0], nexp, rorder, 1, 966 ptr_exp, ptr_coeff, 0]) 967 968 _ecpdic[symb] = (nelec, numpy.asarray(ecp0, dtype=numpy.int32)) 969 970 _ecpbas = [] 971 if _ecpdic: 972 _atm = _atm.copy() 973 for ia, atom in enumerate(mol._atom): 974 symb = atom[0] 975 if symb in _ecpdic: 976 ecp0 = _ecpdic[symb] 977 elif _rm_digit(symb) in _ecpdic: 978 ecp0 = _ecpdic[_rm_digit(symb)] 979 else: 980 ecp0 = None 981 if ecp0 is not None: 982 _atm[ia,CHARGE_OF ] = charge(symb) - ecp0[0] 983 _atm[ia,NUC_MOD_OF] = NUC_ECP 984 b = ecp0[1].copy().reshape(-1,BAS_SLOTS) 985 b[:,ATOM_OF] = ia 986 _ecpbas.append(b) 987 if _ecpbas: 988 _ecpbas = numpy.asarray(numpy.vstack(_ecpbas), numpy.int32) 989 _env = numpy.hstack([pre_env] + _env) 990 else: 991 _ecpbas = numpy.zeros((0,BAS_SLOTS), numpy.int32) 992 _env = pre_env 993 return _atm, _ecpbas, _env 994 995def tot_electrons(mol): 996 '''Total number of electrons for the given molecule 997 998 Returns: 999 electron number in integer 1000 1001 Examples: 1002 1003 >>> mol = gto.M(atom='H 0 1 0; C 0 0 1', charge=1) 1004 >>> mol.tot_electrons() 1005 6 1006 ''' 1007 if mol._atm.size != 0: 1008 nelectron = mol.atom_charges().sum() 1009 elif mol._atom: 1010 nelectron = sum(charge(a[0]) for a in mol._atom) 1011 else: 1012 nelectron = sum(charge(a[0]) for a in format_atom(mol.atom)) 1013 nelectron -= mol.charge 1014 nelectron_int = int(round(nelectron)) 1015 1016 if abs(nelectron - nelectron_int) > 1e-4: 1017 logger.warn(mol, 'Found fractional number of electrons %f. Round it to %d', 1018 nelectron, nelectron_int) 1019 return nelectron_int 1020 1021def copy(mol): 1022 '''Deepcopy of the given :class:`Mole` object 1023 ''' 1024 import copy 1025 newmol = copy.copy(mol) 1026 newmol._atm = numpy.copy(mol._atm) 1027 newmol._bas = numpy.copy(mol._bas) 1028 newmol._env = numpy.copy(mol._env) 1029 newmol._ecpbas = numpy.copy(mol._ecpbas) 1030 1031 newmol.atom = copy.deepcopy(mol.atom) 1032 newmol._atom = copy.deepcopy(mol._atom) 1033 newmol.basis = copy.deepcopy(mol.basis) 1034 newmol._basis = copy.deepcopy(mol._basis) 1035 newmol.ecp = copy.deepcopy(mol.ecp) 1036 newmol._ecp = copy.deepcopy(mol._ecp) 1037 return newmol 1038 1039def pack(mol): 1040 '''Pack the input args of :class:`Mole` to a dict. 1041 1042 Note this function only pack the input arguments (not the entire Mole 1043 class). Modifications to mol._atm, mol._bas, mol._env are not tracked. 1044 Use :func:`dumps` to serialize the entire Mole object. 1045 ''' 1046 mdic = {'atom' : mol.atom, 1047 'unit' : mol.unit, 1048 'basis' : mol.basis, 1049 'charge' : mol.charge, 1050 'spin' : mol.spin, 1051 'symmetry': mol.symmetry, 1052 'nucmod' : mol.nucmod, 1053 'nucprop' : mol.nucprop, 1054 'ecp' : mol.ecp, 1055 '_nelectron': mol._nelectron, 1056 'verbose' : mol.verbose} 1057 return mdic 1058def unpack(moldic): 1059 '''Unpack a dict which is packed by :func:`pack`, to generate the input 1060 arguments for :class:`Mole` object. 1061 ''' 1062 mol = Mole() 1063 mol.__dict__.update(moldic) 1064 return mol 1065 1066 1067def dumps(mol): 1068 '''Serialize Mole object to a JSON formatted str. 1069 ''' 1070 exclude_keys = set(('output', 'stdout', '_keys', 1071 # Constructing in function loads 1072 'symm_orb', 'irrep_id', 'irrep_name')) 1073 nparray_keys = set(('_atm', '_bas', '_env', '_ecpbas', 1074 '_symm_orig', '_symm_axes')) 1075 1076 moldic = dict(mol.__dict__) 1077 for k in exclude_keys: 1078 if k in moldic: 1079 del(moldic[k]) 1080 for k in nparray_keys: 1081 if isinstance(moldic[k], numpy.ndarray): 1082 moldic[k] = moldic[k].tolist() 1083 moldic['atom'] = repr(mol.atom) 1084 moldic['basis']= repr(mol.basis) 1085 moldic['ecp' ] = repr(mol.ecp) 1086 1087 try: 1088 return json.dumps(moldic) 1089 except TypeError: 1090 def skip_value(dic): 1091 dic1 = {} 1092 for k,v in dic.items(): 1093 if (v is None or 1094 isinstance(v, (str, unicode, bool, int, float))): 1095 dic1[k] = v 1096 elif isinstance(v, (list, tuple)): 1097 dic1[k] = v # Should I recursively skip_vaule? 1098 elif isinstance(v, set): 1099 dic1[k] = list(v) 1100 elif isinstance(v, dict): 1101 dic1[k] = skip_value(v) 1102 else: 1103 msg =('Function mol.dumps drops attribute %s because ' 1104 'it is not JSON-serializable' % k) 1105 warnings.warn(msg) 1106 return dic1 1107 return json.dumps(skip_value(moldic), skipkeys=True) 1108 1109def loads(molstr): 1110 '''Deserialize a str containing a JSON document to a Mole object. 1111 ''' 1112 # the numpy function array is used by eval function 1113 from numpy import array # noqa 1114 moldic = json.loads(molstr) 1115 if sys.version_info < (3,): 1116 # Convert to utf8 because JSON loads fucntion returns unicode. 1117 def byteify(inp): 1118 if isinstance(inp, dict): 1119 return dict([(byteify(k), byteify(v)) for k, v in inp.iteritems()]) 1120 elif isinstance(inp, (tuple, list)): 1121 return [byteify(x) for x in inp] 1122 elif isinstance(inp, unicode): 1123 return inp.encode('utf-8') 1124 else: 1125 return inp 1126 moldic = byteify(moldic) 1127 mol = Mole() 1128 mol.__dict__.update(moldic) 1129 mol.atom = eval(mol.atom) 1130 mol.basis= eval(mol.basis) 1131 mol.ecp = eval(mol.ecp) 1132 mol._atm = numpy.array(mol._atm, dtype=numpy.int32) 1133 mol._bas = numpy.array(mol._bas, dtype=numpy.int32) 1134 mol._env = numpy.array(mol._env, dtype=numpy.double) 1135 mol._ecpbas = numpy.array(mol._ecpbas, dtype=numpy.int32) 1136 1137 if mol.symmetry and mol._symm_orig is not None: 1138 from pyscf import symm 1139 mol._symm_orig = numpy.array(mol._symm_orig) 1140 mol._symm_axes = numpy.array(mol._symm_axes) 1141 mol.symm_orb, mol.irrep_id = \ 1142 symm.symm_adapted_basis(mol, mol.groupname, 1143 mol._symm_orig, mol._symm_axes) 1144 mol.irrep_name = [symm.irrep_id2name(mol.groupname, ir) 1145 for ir in mol.irrep_id] 1146 1147 elif mol.symmetry and mol.symm_orb is not None: 1148 # Backward compatibility. To load symm_orb from chkfile of pyscf-1.6 1149 # and earlier. 1150 symm_orb = [] 1151 1152 # decompress symm_orb 1153 for val, x, y, shape in mol.symm_orb: 1154 if isinstance(val[0], list): 1155 # backward compatibility for chkfile of pyscf-1.4 in which val 1156 # is an array of real floats. In pyscf-1.5, val can be a list 1157 # of list, to include the imaginary part 1158 val_real, val_imag = val 1159 else: 1160 val_real, val_imag = val, None 1161 if val_imag is None: 1162 c = numpy.zeros(shape) 1163 c[numpy.array(x),numpy.array(y)] = numpy.array(val_real) 1164 else: 1165 c = numpy.zeros(shape, dtype=numpy.complex128) 1166 val = numpy.array(val_real) + numpy.array(val_imag) * 1j 1167 c[numpy.array(x),numpy.array(y)] = val 1168 symm_orb.append(c) 1169 mol.symm_orb = symm_orb 1170 1171 return mol 1172 1173 1174def len_spinor(l, kappa): 1175 '''The number of spinor associated with given angular momentum and kappa. If kappa is 0, 1176 return 4l+2 1177 ''' 1178 if kappa == 0: 1179 n = (l * 4 + 2) 1180 elif kappa < 0: 1181 n = (l * 2 + 2) 1182 else: 1183 n = (l * 2) 1184 return n 1185 1186def len_cart(l): 1187 '''The number of Cartesian function associated with given angular momentum. 1188 ''' 1189 return (l + 1) * (l + 2) // 2 1190 1191def npgto_nr(mol, cart=None): 1192 '''Total number of primitive spherical GTOs for the given :class:`Mole` object''' 1193 if cart is None: 1194 cart = mol.cart 1195 l = mol._bas[:,ANG_OF] 1196 if cart: 1197 return ((l+1)*(l+2)//2 * mol._bas[:,NPRIM_OF]).sum() 1198 else: 1199 return ((l*2+1) * mol._bas[:,NPRIM_OF]).sum() 1200def nao_nr(mol, cart=None): 1201 '''Total number of contracted GTOs for the given :class:`Mole` object''' 1202 if cart is None: 1203 cart = mol.cart 1204 if cart: 1205 return nao_cart(mol) 1206 else: 1207 return ((mol._bas[:,ANG_OF]*2+1) * mol._bas[:,NCTR_OF]).sum() 1208def nao_cart(mol): 1209 '''Total number of contracted cartesian GTOs for the given :class:`Mole` object''' 1210 l = mol._bas[:,ANG_OF] 1211 return ((l+1)*(l+2)//2 * mol._bas[:,NCTR_OF]).sum() 1212 1213# nao_id0:nao_id1 corresponding to bas_id0:bas_id1 1214def nao_nr_range(mol, bas_id0, bas_id1): 1215 '''Lower and upper boundary of contracted spherical basis functions associated 1216 with the given shell range 1217 1218 Args: 1219 mol : 1220 :class:`Mole` object 1221 bas_id0 : int 1222 start shell id 1223 bas_id1 : int 1224 stop shell id 1225 1226 Returns: 1227 tupel of start basis function id and the stop function id 1228 1229 Examples: 1230 1231 >>> mol = gto.M(atom='O 0 0 0; C 0 0 1', basis='6-31g') 1232 >>> gto.nao_nr_range(mol, 2, 4) 1233 (2, 6) 1234 ''' 1235 ao_loc = moleintor.make_loc(mol._bas[:bas_id1], 'sph') 1236 nao_id0 = ao_loc[bas_id0] 1237 nao_id1 = ao_loc[-1] 1238 return nao_id0, nao_id1 1239 1240def nao_2c(mol): 1241 '''Total number of contracted spinor GTOs for the given :class:`Mole` object''' 1242 l = mol._bas[:,ANG_OF] 1243 kappa = mol._bas[:,KAPPA_OF] 1244 dims = (l*4+2) * mol._bas[:,NCTR_OF] 1245 dims[kappa<0] = l[kappa<0] * 2 + 2 1246 dims[kappa>0] = l[kappa>0] * 2 1247 return dims.sum() 1248 1249# nao_id0:nao_id1 corresponding to bas_id0:bas_id1 1250def nao_2c_range(mol, bas_id0, bas_id1): 1251 '''Lower and upper boundary of contracted spinor basis functions associated 1252 with the given shell range 1253 1254 Args: 1255 mol : 1256 :class:`Mole` object 1257 bas_id0 : int 1258 start shell id, 0-based 1259 bas_id1 : int 1260 stop shell id, 0-based 1261 1262 Returns: 1263 tupel of start basis function id and the stop function id 1264 1265 Examples: 1266 1267 >>> mol = gto.M(atom='O 0 0 0; C 0 0 1', basis='6-31g') 1268 >>> gto.nao_2c_range(mol, 2, 4) 1269 (4, 12) 1270 ''' 1271 ao_loc = moleintor.make_loc(mol._bas[:bas_id1], '') 1272 nao_id0 = ao_loc[bas_id0] 1273 nao_id1 = ao_loc[-1] 1274 return nao_id0, nao_id1 1275 1276def ao_loc_nr(mol, cart=None): 1277 '''Offset of every shell in the spherical basis function spectrum 1278 1279 Returns: 1280 list, each entry is the corresponding start basis function id 1281 1282 Examples: 1283 1284 >>> mol = gto.M(atom='O 0 0 0; C 0 0 1', basis='6-31g') 1285 >>> gto.ao_loc_nr(mol) 1286 [0, 1, 2, 3, 6, 9, 10, 11, 12, 15, 18] 1287 ''' 1288 if cart is None: 1289 cart = mol.cart 1290 if cart: 1291 return moleintor.make_loc(mol._bas, 'cart') 1292 else: 1293 return moleintor.make_loc(mol._bas, 'sph') 1294 1295def ao_loc_2c(mol): 1296 '''Offset of every shell in the spinor basis function spectrum 1297 1298 Returns: 1299 list, each entry is the corresponding start id of spinor function 1300 1301 Examples: 1302 1303 >>> mol = gto.M(atom='O 0 0 0; C 0 0 1', basis='6-31g') 1304 >>> gto.ao_loc_2c(mol) 1305 [0, 2, 4, 6, 12, 18, 20, 22, 24, 30, 36] 1306 ''' 1307 return moleintor.make_loc(mol._bas, 'spinor') 1308 1309def time_reversal_map(mol): 1310 r'''The index to map the spinor functions and its time reversal counterpart. 1311 The returned indices have postive or negative values. For the i-th basis function, 1312 if the returned j = idx[i] < 0, it means :math:`T|i\rangle = -|j\rangle`, 1313 otherwise :math:`T|i\rangle = |j\rangle` 1314 ''' 1315 tao = [] 1316 i = 0 1317 for b in mol._bas: 1318 l = b[ANG_OF] 1319 if b[KAPPA_OF] == 0: 1320 djs = (l * 2, l * 2 + 2) 1321 elif b[KAPPA_OF] > 0: 1322 djs = (l * 2,) 1323 else: 1324 djs = (l * 2 + 2,) 1325 if l % 2 == 0: 1326 for n in range(b[NCTR_OF]): 1327 for dj in djs: 1328 for m in range(0, dj, 2): 1329 tao.append(-(i + dj - m)) 1330 tao.append( i + dj - m - 1) 1331 i += dj 1332 else: 1333 for n in range(b[NCTR_OF]): 1334 for dj in djs: 1335 for m in range(0, dj, 2): 1336 tao.append( i + dj - m) 1337 tao.append(-(i + dj - m - 1)) 1338 i += dj 1339 return numpy.asarray(tao, dtype=numpy.int32) 1340 1341 1342CHECK_GEOM = getattr(__config__, 'gto_mole_check_geom', True) 1343 1344def energy_nuc(mol, charges=None, coords=None): 1345 '''Compute nuclear repulsion energy (AU) or static Coulomb energy 1346 1347 Returns 1348 float 1349 ''' 1350 if charges is None: charges = mol.atom_charges() 1351 if len(charges) <= 1: 1352 return 0 1353 rr = inter_distance(mol, coords) 1354 rr[numpy.diag_indices_from(rr)] = 1e200 1355 if CHECK_GEOM and numpy.any(rr < 1e-5): 1356 for atm_idx in numpy.argwhere(rr<1e-5): 1357 logger.warn(mol, 'Atoms %s have the same coordinates', atm_idx) 1358 raise RuntimeError('Ill geometry') 1359 e = numpy.einsum('i,ij,j->', charges, 1./rr, charges) * .5 1360 return e 1361 1362def inter_distance(mol, coords=None): 1363 ''' 1364 Inter-particle distance array 1365 ''' 1366 if coords is None: coords = mol.atom_coords() 1367 rr = numpy.linalg.norm(coords.reshape(-1,1,3) - coords, axis=2) 1368 rr[numpy.diag_indices_from(rr)] = 0 1369 return rr 1370 1371def sph_labels(mol, fmt=True, base=BASE): 1372 '''Labels for spherical GTO functions 1373 1374 Kwargs: 1375 fmt : str or bool 1376 if fmt is boolean, it controls whether to format the labels and the 1377 default format is "%d%3s %s%-4s". if fmt is string, the string will 1378 be used as the print format. 1379 1380 Returns: 1381 List of [(atom-id, symbol-str, nl-str, str-of-real-spherical-notation] 1382 or formatted strings based on the argument "fmt" 1383 1384 Examples: 1385 1386 >>> mol = gto.M(atom='H 0 0 0; Cl 0 0 1', basis='sto-3g') 1387 >>> gto.sph_labels(mol) 1388 [(0, 'H', '1s', ''), (1, 'Cl', '1s', ''), (1, 'Cl', '2s', ''), (1, 'Cl', '3s', ''), 1389 (1, 'Cl', '2p', 'x'), (1, 'Cl', '2p', 'y'), (1, 'Cl', '2p', 'z'), (1, 'Cl', '3p', 'x'), 1390 (1, 'Cl', '3p', 'y'), (1, 'Cl', '3p', 'z')] 1391 ''' 1392 count = numpy.zeros((mol.natm, 9), dtype=int) 1393 label = [] 1394 for ib in range(mol.nbas): 1395 ia = mol.bas_atom(ib) 1396 l = mol.bas_angular(ib) 1397 strl = param.ANGULAR[l] 1398 nc = mol.bas_nctr(ib) 1399 symb = mol.atom_symbol(ia) 1400 nelec_ecp = mol.atom_nelec_core(ia) 1401 if nelec_ecp == 0 or l > 3: 1402 shl_start = count[ia,l]+l+1 1403 else: 1404 coreshl = core_configuration(nelec_ecp) 1405 shl_start = coreshl[l]+count[ia,l]+l+1 1406 count[ia,l] += nc 1407 for n in range(shl_start, shl_start+nc): 1408 for m in range(-l, l+1): 1409 label.append((ia+base, symb, '%d%s' % (n, strl), 1410 str(param.REAL_SPHERIC[l][l+m]))) 1411 1412 if isinstance(fmt, (str, unicode)): 1413 return [(fmt % x) for x in label] 1414 elif fmt: 1415 return ['%d %s %s%-4s' % x for x in label] 1416 else: 1417 return label 1418spheric_labels = spherical_labels = sph_labels 1419 1420def cart_labels(mol, fmt=True, base=BASE): 1421 '''Labels of Cartesian GTO functions 1422 1423 Kwargs: 1424 fmt : str or bool 1425 if fmt is boolean, it controls whether to format the labels and the 1426 default format is "%d%3s %s%-4s". if fmt is string, the string will 1427 be used as the print format. 1428 1429 Returns: 1430 List of [(atom-id, symbol-str, nl-str, str-of-xyz-notation)] 1431 or formatted strings based on the argument "fmt" 1432 ''' 1433 cartxyz = [] 1434 for l in range(max(mol._bas[:,ANG_OF])+1): 1435 xyz = [] 1436 for x in range(l, -1, -1): 1437 for y in range(l-x, -1, -1): 1438 z = l-x-y 1439 xyz.append('x'*x + 'y'*y + 'z'*z) 1440 cartxyz.append(xyz) 1441 1442 count = numpy.zeros((mol.natm, 9), dtype=int) 1443 label = [] 1444 for ib in range(len(mol._bas)): 1445 ia = mol.bas_atom(ib) 1446 l = mol.bas_angular(ib) 1447 strl = param.ANGULAR[l] 1448 nc = mol.bas_nctr(ib) 1449 symb = mol.atom_symbol(ia) 1450 nelec_ecp = mol.atom_nelec_core(ia) 1451 if nelec_ecp == 0 or l > 3: 1452 shl_start = count[ia,l]+l+1 1453 else: 1454 coreshl = core_configuration(nelec_ecp) 1455 shl_start = coreshl[l]+count[ia,l]+l+1 1456 count[ia,l] += nc 1457 ncart = (l + 1) * (l + 2) // 2 1458 for n in range(shl_start, shl_start+nc): 1459 for m in range(ncart): 1460 label.append((ia+base, symb, '%d%s' % (n, strl), cartxyz[l][m])) 1461 1462 if isinstance(fmt, (str, unicode)): 1463 return [(fmt % x) for x in label] 1464 elif fmt: 1465 return ['%d%3s %s%-4s' % x for x in label] 1466 else: 1467 return label 1468 1469def ao_labels(mol, fmt=True, base=BASE): 1470 '''Labels of AO basis functions 1471 1472 Kwargs: 1473 fmt : str or bool 1474 if fmt is boolean, it controls whether to format the labels and the 1475 default format is "%d%3s %s%-4s". if fmt is string, the string will 1476 be used as the print format. 1477 1478 Returns: 1479 List of [(atom-id, symbol-str, nl-str, str-of-AO-notation)] 1480 or formatted strings based on the argument "fmt" 1481 ''' 1482 if mol.cart: 1483 return mol.cart_labels(fmt, base) 1484 else: 1485 return mol.sph_labels(fmt, base) 1486 1487def spinor_labels(mol, fmt=True, base=BASE): 1488 ''' 1489 Labels of spinor GTO functions 1490 ''' 1491 count = numpy.zeros((mol.natm, 9), dtype=int) 1492 label = [] 1493 for ib in range(mol.nbas): 1494 ia = mol.bas_atom(ib) 1495 l = mol.bas_angular(ib) 1496 kappa = mol.bas_kappa(ib) 1497 strl = param.ANGULAR[l] 1498 nc = mol.bas_nctr(ib) 1499 symb = mol.atom_symbol(ia) 1500 nelec_ecp = mol.atom_nelec_core(ia) 1501 if nelec_ecp == 0 or l > 3: 1502 shl_start = count[ia,l]+l+1 1503 else: 1504 coreshl = core_configuration(nelec_ecp) 1505 shl_start = coreshl[l]+count[ia,l]+l+1 1506 count[ia,l] += nc 1507 for n in range(shl_start, shl_start+nc): 1508 if kappa >= 0: 1509 for m in range(-l*2+1, l*2, 2): 1510 label.append((ia+base, symb, '%d%s%d/2' % (n, strl, l*2-1), 1511 '%d/2'%m)) 1512 if kappa <= 0: 1513 for m in range(-l*2-1, l*2+2, 2): 1514 label.append((ia+base, symb, '%d%s%d/2' % (n, strl, l*2+1), 1515 '%d/2'%m)) 1516 1517 if isinstance(fmt, (str, unicode)): 1518 return [(fmt % x) for x in label] 1519 elif fmt: 1520 return ['%d %s %s,%-5s' % x for x in label] 1521 else: 1522 return label 1523 1524def search_ao_label(mol, label): 1525 '''Find the index of the AO basis function based on the given ao_label 1526 1527 Args: 1528 ao_label : string or a list of strings 1529 The regular expression pattern to match the orbital labels 1530 returned by mol.ao_labels() 1531 1532 Returns: 1533 A list of index for the AOs that matches the given ao_label RE pattern 1534 1535 Examples: 1536 1537 >>> mol = gto.M(atom='H 0 0 0; Cl 0 0 1', basis='ccpvtz') 1538 >>> mol.search_ao_label('Cl.*p') 1539 [19 20 21 22 23 24 25 26 27 28 29 30] 1540 >>> mol.search_ao_label('Cl 2p') 1541 [19 20 21] 1542 >>> mol.search_ao_label(['Cl.*d', 'Cl 4p']) 1543 [25 26 27 31 32 33 34 35 36 37 38 39 40] 1544 ''' 1545 return _aolabels2baslst(mol, label) 1546 1547def _aolabels2baslst(mol, aolabels_or_baslst, base=BASE): 1548 if callable(aolabels_or_baslst): 1549 baslst = [i for i,x in enumerate(mol.ao_labels(base=base)) 1550 if aolabels_or_baslst(x)] 1551 elif isinstance(aolabels_or_baslst, str): 1552 aolabels = re.sub(' +', ' ', aolabels_or_baslst.strip(), count=1) 1553 aolabels = re.compile(aolabels) 1554 baslst = [i for i,s in enumerate(mol.ao_labels(base=base)) 1555 if re.search(aolabels, s)] 1556 elif len(aolabels_or_baslst) > 0 and isinstance(aolabels_or_baslst[0], str): 1557 aolabels = [re.compile(re.sub(' +', ' ', x.strip(), count=1)) 1558 for x in aolabels_or_baslst] 1559 baslst = [i for i,t in enumerate(mol.ao_labels(base=base)) 1560 if any(re.search(x, t) for x in aolabels)] 1561 else: 1562 baslst = [i-base for i in aolabels_or_baslst] 1563 return numpy.asarray(baslst, dtype=int) 1564 1565def search_shell_id(mol, atm_id, l): 1566 '''Search the first basis/shell id (**not** the basis function id) which 1567 matches the given atom-id and angular momentum 1568 1569 Args: 1570 atm_id : int 1571 atom id, 0-based 1572 l : int 1573 angular momentum 1574 1575 Returns: 1576 basis id, 0-based. If not found, return None 1577 1578 Examples: 1579 1580 >>> mol = gto.M(atom='H 0 0 0; Cl 0 0 1', basis='sto-3g') 1581 >>> mol.search_shell_id(1, 1) # Cl p shell 1582 4 1583 >>> mol.search_shell_id(1, 2) # Cl d shell 1584 None 1585 ''' 1586 for ib in range(len(mol._bas)): 1587 ia = mol.bas_atom(ib) 1588 l1 = mol.bas_angular(ib) 1589 if ia == atm_id and l1 == l: 1590 return ib 1591 1592def search_ao_nr(mol, atm_id, l, m, atmshell): 1593 '''Search the first basis function id (**not** the shell id) which matches 1594 the given atom-id, angular momentum magnetic angular momentum, principal shell. 1595 1596 Args: 1597 atm_id : int 1598 atom id, 0-based 1599 l : int 1600 angular momentum 1601 m : int 1602 magnetic angular momentum 1603 atmshell : int 1604 principal quantum number 1605 1606 Returns: 1607 basis function id, 0-based. If not found, return None 1608 1609 Examples: 1610 1611 >>> mol = gto.M(atom='H 0 0 0; Cl 0 0 1', basis='sto-3g') 1612 >>> mol.search_ao_nr(1, 1, -1, 3) # Cl 3px 1613 7 1614 ''' 1615 ibf = 0 1616 for ib in range(len(mol._bas)): 1617 ia = mol.bas_atom(ib) 1618 l1 = mol.bas_angular(ib) 1619 if mol.cart: 1620 degen = (l1 + 1) * (l1 + 2) // 2 1621 else: 1622 degen = l1 * 2 + 1 1623 nc = mol.bas_nctr(ib) 1624 if ia == atm_id and l1 == l: 1625 if atmshell > nc+l1: 1626 atmshell = atmshell - nc 1627 else: 1628 return ibf + (atmshell-l1-1)*degen + (l1+m) 1629 ibf += degen * nc 1630 raise RuntimeError('Required AO not found') 1631 1632def search_ao_r(mol, atm_id, l, j, m, atmshell): 1633 raise RuntimeError('TODO') 1634#TODO: ibf = 0 1635#TODO: for ib in range(len(mol._bas)): 1636#TODO: ia = mol.bas_atom(ib) 1637#TODO: l1 = mol.bas_angular(ib) 1638#TODO: nc = mol.bas_nctr(ib) 1639#TODO: k = mol.bas_kappa(bas_id) 1640#TODO: degen = len_spinor(l1, k) 1641#TODO: if ia == atm_id and l1 == l and k == kappa: 1642#TODO: if atmshell > nc+l1: 1643#TODO: atmshell = atmshell - nc 1644#TODO: else: 1645#TODO: return ibf + (atmshell-l1-1)*degen + (degen+m) 1646#TODO: ibf += degen 1647 1648def offset_2c_by_atom(mol): 1649 '''2-component AO offset for each atom. Return a list, each item 1650 of the list gives (start-shell-id, stop-shell-id, start-AO-id, stop-AO-id) 1651 ''' 1652 return aoslice_by_atom(mol, mol.ao_loc_2c()) 1653 1654def aoslice_by_atom(mol, ao_loc=None): 1655 '''AO offsets for each atom. Return a list, each item of the list gives 1656 (start-shell-id, stop-shell-id, start-AO-id, stop-AO-id) 1657 ''' 1658 if ao_loc is None: 1659 ao_loc = mol.ao_loc_nr() 1660 1661 aorange = numpy.empty((mol.natm,4), dtype=int) 1662 1663 if mol.natm == 0: 1664 return aorange 1665 1666 bas_atom = mol._bas[:,ATOM_OF] 1667 delimiter = numpy.where(bas_atom[0:-1] != bas_atom[1:])[0] + 1 1668 1669 if mol.natm == len(delimiter) + 1: 1670 aorange[:,0] = shell_start = numpy.append(0, delimiter) 1671 aorange[:,1] = shell_end = numpy.append(delimiter, mol.nbas) 1672 1673 else: # Some atoms miss basis 1674 shell_start = numpy.empty(mol.natm, dtype=int) 1675 shell_start[:] = -1 1676 shell_start[0] = 0 1677 shell_start[bas_atom[0]] = 0 1678 shell_start[bas_atom[delimiter]] = delimiter 1679 1680 shell_end = numpy.empty(mol.natm, dtype=int) 1681 shell_end[0] = 0 1682 shell_end[bas_atom[delimiter-1]] = delimiter 1683 shell_end[bas_atom[-1]] = mol.nbas 1684 1685 for i in range(1, mol.natm): 1686 if shell_start[i] == -1: 1687 shell_start[i] = shell_end[i] = shell_end[i-1] 1688 1689 aorange[:,0] = shell_start 1690 aorange[:,1] = shell_end 1691 aorange[:,2] = ao_loc[shell_start] 1692 aorange[:,3] = ao_loc[shell_end] 1693 return aorange 1694offset_nr_by_atom = aoslice_by_atom 1695 1696def same_basis_set(mol1, mol2): 1697 '''Check whether two molecules use the same basis sets. 1698 The two molecules can have different geometry. 1699 ''' 1700 atomtypes1 = atom_types(mol1._atom, mol1._basis) 1701 atomtypes2 = atom_types(mol2._atom, mol2._basis) 1702 if set(atomtypes1.keys()) != set(atomtypes2.keys()): 1703 return False 1704 for k in atomtypes1: 1705 if len(atomtypes1[k]) != len(atomtypes2[k]): 1706 return False 1707 elif mol1._basis.get(k, None) != mol2._basis.get(k, None): 1708 return False 1709 return True 1710 1711def same_mol(mol1, mol2, tol=1e-5, cmp_basis=True, ignore_chiral=False): 1712 '''Compare the two molecules whether they have the same structure. 1713 1714 Kwargs: 1715 tol : float 1716 In Bohr 1717 cmp_basis : bool 1718 Whether to compare basis functions for the two molecules 1719 ''' 1720 from pyscf import symm 1721 1722 if mol1._atom.__len__() != mol2._atom.__len__(): 1723 return False 1724 1725 chg1 = mol1._atm[:,CHARGE_OF] 1726 chg2 = mol2._atm[:,CHARGE_OF] 1727 if not numpy.all(numpy.sort(chg1) == numpy.sort(chg2)): 1728 return False 1729 1730 if cmp_basis and not same_basis_set(mol1, mol2): 1731 return False 1732 1733 def finger(mol, chgs, coord): 1734 center = numpy.einsum('z,zr->r', chgs, coord) / chgs.sum() 1735 im = inertia_moment(mol, chgs, coord) 1736 # Divid im by chgs.sum(), to normalize im. Otherwise the input tol may 1737 # not reflect the actual deviation. 1738 im /= chgs.sum() 1739 1740 w, v = scipy.linalg.eigh(im) 1741 axes = v.T 1742 if numpy.linalg.det(axes) < 0: 1743 axes *= -1 1744 r = numpy.dot(coord-center, axes.T) 1745 return w, r 1746 1747 coord1 = mol1.atom_coords() 1748 coord2 = mol2.atom_coords() 1749 w1, r1 = finger(mol1, chg1, coord1) 1750 w2, r2 = finger(mol2, chg2, coord2) 1751 if not (numpy.allclose(w1, w2, atol=tol)): 1752 return False 1753 1754 rotate_xy = numpy.array([[-1., 0., 0.], 1755 [ 0.,-1., 0.], 1756 [ 0., 0., 1.]]) 1757 rotate_yz = numpy.array([[ 1., 0., 0.], 1758 [ 0.,-1., 0.], 1759 [ 0., 0.,-1.]]) 1760 rotate_zx = numpy.array([[-1., 0., 0.], 1761 [ 0., 1., 0.], 1762 [ 0., 0.,-1.]]) 1763 1764 def inspect(z1, r1, z2, r2): 1765 place = int(-numpy.log10(tol)) - 1 1766 idx = symm.argsort_coords(r2, place) 1767 z2 = z2[idx] 1768 r2 = r2[idx] 1769 for rot in (1, rotate_xy, rotate_yz, rotate_zx): 1770 r1new = numpy.dot(r1, rot) 1771 idx = symm.argsort_coords(r1new, place) 1772 if (numpy.all(z1[idx] == z2) and 1773 numpy.allclose(r1new[idx], r2, atol=tol)): 1774 return True 1775 return False 1776 1777 return (inspect(chg1, r1, chg2, r2) or 1778 (ignore_chiral and inspect(chg1, r1, chg2, -r2))) 1779is_same_mol = same_mol 1780 1781def chiral_mol(mol1, mol2=None): 1782 '''Detect whether the given molelcule is chiral molecule or two molecules 1783 are chiral isomers. 1784 ''' 1785 if mol2 is None: 1786 mol2 = mol1.copy() 1787 ptr_coord = mol2._atm[:,PTR_COORD] 1788 mol2._env[ptr_coord ] *= -1 1789 mol2._env[ptr_coord+1] *= -1 1790 mol2._env[ptr_coord+2] *= -1 1791 return (not same_mol(mol1, mol2, ignore_chiral=False) and 1792 same_mol(mol1, mol2, ignore_chiral=True)) 1793 1794def inertia_moment(mol, mass=None, coords=None): 1795 if mass is None: 1796 mass = mol.atom_mass_list() 1797 if coords is None: 1798 coords = mol.atom_coords() 1799 mass_center = numpy.einsum('i,ij->j', mass, coords)/mass.sum() 1800 coords = coords - mass_center 1801 im = numpy.einsum('i,ij,ik->jk', mass, coords, coords) 1802 im = numpy.eye(3) * im.trace() - im 1803 return im 1804 1805def atom_mass_list(mol, isotope_avg=False): 1806 '''A list of mass for all atoms in the molecule 1807 1808 Kwargs: 1809 isotope_avg : boolean 1810 Whether to use the isotope average mass as the atomic mass 1811 ''' 1812 if isotope_avg: 1813 mass_table = elements.MASSES 1814 else: 1815 mass_table = elements.ISOTOPE_MAIN 1816 1817 nucprop = mol.nucprop 1818 if nucprop: 1819 mass = [] 1820 for ia in range(mol.natm): 1821 z = mol.atom_charge(ia) 1822 symb = mol.atom_symbol(ia) 1823 stdsymb = _std_symbol(symb) 1824 if ia+1 in nucprop: 1825 prop = nucprop[ia+1] 1826 elif symb in nucprop: 1827 prop = nucprop[symb] 1828 else: 1829 prop = nucprop.get(stdsymb, {}) 1830 1831 mass.append(prop.get('mass', mass_table[z])) 1832 else: 1833 #mass = [mass_table[z] for z in mol.atom_charges()] 1834 mass = [] 1835 for ia in range(mol.natm): 1836 z = charge(mol.atom_symbol(ia)) 1837 mass.append(mass_table[z]) 1838 1839 return numpy.array(mass) 1840 1841def condense_to_shell(mol, mat, compressor='max'): 1842 '''The given matrix is first partitioned to blocks, based on AO shell as 1843 delimiter. Then call compressor function to abstract each block. 1844 ''' 1845 ao_loc = mol.ao_loc_nr() 1846 if callable(compressor): 1847 abstract = numpy.empty((mol.nbas, mol.nbas)) 1848 for i, i0 in enumerate(ao_loc[:mol.nbas]): 1849 for j, j0 in enumerate(ao_loc[:mol.nbas]): 1850 abstract[i,j] = compressor(mat[i0:ao_loc[i+1],j0:ao_loc[j+1]]) 1851 else: 1852 abstract = lib.condense(compressor, mat, ao_loc) 1853 return abstract 1854 1855 1856def tostring(mol, format='raw'): 1857 '''Convert molecular geometry to a string of the required format. 1858 1859 Supported output formats: 1860 | raw: Each line is <symobl> <x> <y> <z> 1861 | xyz: XYZ cartesian coordinates format 1862 | zmat: Z-matrix format 1863 ''' 1864 format = format.lower() 1865 if format == 'xyz' or format == 'raw': 1866 coords = mol.atom_coords() * param.BOHR 1867 output = [] 1868 if format == 'xyz': 1869 output.append('%d' % mol.natm) 1870 output.append('XYZ from PySCF') 1871 1872 for i in range(mol.natm): 1873 symb = mol.atom_pure_symbol(i) 1874 x, y, z = coords[i] 1875 output.append('%-4s %14.5f %14.5f %14.5f' % 1876 (symb, x, y, z)) 1877 return '\n'.join(output) 1878 elif format == 'zmat': 1879 coords = mol.atom_coords() * param.BOHR 1880 zmat = cart2zmat(coords).splitlines() 1881 output = [] 1882 for i, line in enumerate(zmat): 1883 symb = mol.atom_pure_symbol(i) 1884 output.append('%-4s %s' % (symb, line)) 1885 return '\n'.join(output) 1886 else: 1887 raise NotImplementedError 1888 1889def tofile(mol, filename, format=None): 1890 '''Write molecular geometry to a file of the required format. 1891 1892 Supported output formats: 1893 | raw: Each line is <symobl> <x> <y> <z> 1894 | xyz: XYZ cartesian coordinates format 1895 | zmat: Z-matrix format 1896 ''' 1897 if format is None: # Guess format based on filename 1898 format = os.path.splitext(filename)[1][1:] 1899 string = tostring(mol, format) 1900 with open(filename, 'w') as f: 1901 f.write(string) 1902 f.write('\n') 1903 return string 1904 1905def fromfile(filename, format=None): 1906 '''Read molecular geometry from a file 1907 (in testing) 1908 1909 Supported formats: 1910 | raw: Each line is <symobl> <x> <y> <z> 1911 | xyz: XYZ cartesian coordinates format 1912 | zmat: Z-matrix format 1913 ''' 1914 if format is None: # Guess format based on filename 1915 format = os.path.splitext(filename)[1][1:].lower() 1916 if format not in ('xyz', 'zmat', 'sdf'): 1917 format = 'raw' 1918 with open(filename, 'r') as f: 1919 return fromstring(f.read(), format) 1920 1921 1922def fromstring(string, format='xyz'): 1923 '''Convert the string of the specified format to internal format 1924 (in testing) 1925 1926 Supported formats: 1927 | raw: Each line is <symobl> <x> <y> <z> 1928 | xyz: XYZ cartesian coordinates format 1929 | zmat: Z-matrix format 1930 ''' 1931 format = format.lower() 1932 if format == 'zmat': 1933 return from_zmatrix(string) 1934 elif format == 'xyz': 1935 dat = string.splitlines() 1936 natm = int(dat[0]) 1937 return '\n'.join(dat[2:natm+2]) 1938 elif format == 'sdf': 1939 raw = string.splitlines() 1940 natoms, nbonds = raw[3].split()[:2] 1941 atoms = [] 1942 for line in raw[4:4+int(natoms)]: 1943 d = line.split() 1944 atoms.append('%s %s %s %s' % (d[3], d[0], d[1], d[2])) 1945 return '\n'.join(atoms) 1946 elif format == 'raw': 1947 return string 1948 else: 1949 raise NotImplementedError 1950 1951 1952# 1953# Mole class handles three layers: input, internal format, libcint arguments. 1954# The relationship of the three layers are, eg 1955# .atom (input) <=> ._atom (for python) <=> ._atm (for libcint) 1956# .basis (input) <=> ._basis (for python) <=> ._bas (for libcint) 1957# input layer does not talk to libcint directly. Data are held in python 1958# internal fomrat layer. Most of methods defined in this class only operates 1959# on the internal format. Exceptions are make_env, make_atm_env, make_bas_env, 1960# set_common_orig_, set_rinv_orig_ which are used to manipulate the libcint arguments. 1961# 1962class Mole(lib.StreamObject): 1963 '''Basic class to hold molecular structure and global options 1964 1965 Attributes: 1966 verbose : int 1967 Print level 1968 output : str or None 1969 Output file, default is None which dumps msg to sys.stdout 1970 max_memory : int, float 1971 Allowed memory in MB 1972 charge : int 1973 Charge of molecule. It affects the electron numbers 1974 spin : int or None 1975 2S, num. alpha electrons - num. beta electrons to control 1976 multiplicity. If spin = None is set, multiplicity will be guessed 1977 based on the neutral molecule. 1978 symmetry : bool or str 1979 Whether to use symmetry. When this variable is set to True, the 1980 molecule will be rotated and the highest rotation axis will be 1981 placed z-axis. 1982 If a string is given as the name of point group, the given point 1983 group symmetry will be used. Note that the input molecular 1984 coordinates will not be changed in this case. 1985 symmetry_subgroup : str 1986 subgroup 1987 1988 atom : list or str 1989 To define molecluar structure. The internal format is 1990 1991 | atom = [[atom1, (x, y, z)], 1992 | [atom2, (x, y, z)], 1993 | ... 1994 | [atomN, (x, y, z)]] 1995 1996 unit : str 1997 Angstrom or Bohr 1998 basis : dict or str 1999 To define basis set. 2000 nucmod : dict or str or [function(nuc_charge, nucprop) => zeta] 2001 Nuclear model. 0 or None means point nuclear model. Other 2002 values will enable Gaussian nuclear model. If a function is 2003 assigned to this attribute, the function will be called to 2004 generate the nuclear charge distribution value "zeta" and the 2005 relevant nuclear model will be set to Gaussian model. 2006 Default is point nuclear model. 2007 nucprop : dict 2008 Nuclear properties (like g-factor 'g', quadrupole moments 'Q'). 2009 It is needed by pyscf.prop module and submodules. 2010 cart : boolean 2011 Using Cartesian GTO basis and integrals (6d,10f,15g) 2012 2013 ** Following attributes are generated by :func:`Mole.build` ** 2014 2015 stdout : file object 2016 Default is sys.stdout if :attr:`Mole.output` is not set 2017 topgroup : str 2018 Point group of the system. 2019 groupname : str 2020 The supported subgroup of the point group. It can be one of SO3, 2021 Dooh, Coov, D2h, C2h, C2v, D2, Cs, Ci, C2, C1 2022 nelectron : int 2023 sum of nuclear charges - :attr:`Mole.charge` 2024 symm_orb : a list of numpy.ndarray 2025 Symmetry adapted basis. Each element is a set of symm-adapted orbitals 2026 for one irreducible representation. The list index does **not** correspond 2027 to the id of irreducible representation. 2028 irrep_id : a list of int 2029 Each element is one irreducible representation id associated with the basis 2030 stored in symm_orb. One irrep id stands for one irreducible representation 2031 symbol. The irrep symbol and the relevant id are defined in 2032 :attr:`symm.param.IRREP_ID_TABLE` 2033 irrep_name : a list of str 2034 Each element is one irreducible representation symbol associated with the basis 2035 stored in symm_orb. The irrep symbols are defined in 2036 :attr:`symm.param.IRREP_ID_TABLE` 2037 _built : bool 2038 To label whether :func:`Mole.build` has been called. It is to 2039 ensure certain functions being initialized only once. 2040 _basis : dict 2041 like :attr:`Mole.basis`, the internal format which is returned from the 2042 parser :func:`format_basis` 2043 _keys : a set of str 2044 Store the keys appeared in the module. It is used to check misinput attributes 2045 2046 ** Following attributes are arguments used by ``libcint`` library ** 2047 2048 _atm : 2049 :code:`[[charge, ptr-of-coord, nuc-model, ptr-zeta, 0, 0], [...]]` 2050 each element reperesents one atom 2051 natm : 2052 number of atoms 2053 _bas : 2054 :code:`[[atom-id, angular-momentum, num-primitive-GTO, num-contracted-GTO, 0, ptr-of-exps, ptr-of-contract-coeff, 0], [...]]` 2055 each element reperesents one shell 2056 nbas : 2057 number of shells 2058 _env : 2059 list of floats to store the coordinates, GTO exponents, contract-coefficients 2060 2061 Examples: 2062 2063 >>> mol = Mole(atom='H^2 0 0 0; H 0 0 1.1', basis='sto3g').build() 2064 >>> print(mol.atom_symbol(0)) 2065 H^2 2066 >>> print(mol.atom_pure_symbol(0)) 2067 H 2068 >>> print(mol.nao_nr()) 2069 2 2070 >>> print(mol.intor('int1e_ovlp_sph')) 2071 [[ 0.99999999 0.43958641] 2072 [ 0.43958641 0.99999999]] 2073 >>> mol.charge = 1 2074 >>> mol.build() 2075 <class 'pyscf.gto.mole.Mole'> has no attributes Charge 2076 ''' # noqa: E501 2077 2078 verbose = getattr(__config__, 'VERBOSE', logger.NOTE) 2079 2080 # the unit (angstrom/bohr) of the coordinates defined by the input self.atom 2081 unit = getattr(__config__, 'UNIT', 'angstrom') 2082 2083 # Whether to hold everything in memory 2084 incore_anyway = getattr(__config__, 'INCORE_ANYWAY', False) 2085 2086 # Using cartesian GTO (6d,10f,15g) 2087 cart = getattr(__config__, 'gto_mole_Mole_cart', False) 2088 2089 def __init__(self, **kwargs): 2090 self.output = None 2091 self.max_memory = param.MAX_MEMORY 2092 2093 self.charge = 0 2094 self.spin = 0 # 2j == nelec_alpha - nelec_beta 2095 self.symmetry = False 2096 self.symmetry_subgroup = None 2097 self.cart = False 2098 2099# Save inputs 2100# self.atom = [(symb/nuc_charge, (coord(Angstrom):0.,0.,0.)), ...] 2101 self.atom = [] 2102# self.basis = {atom_type/nuc_charge: [l, kappa, (expnt, c_1, c_2,..),..]} 2103 self.basis = 'sto-3g' 2104# self.nucmod = {atom_symbol: nuclear_model, atom_id: nuc_mod}, atom_id is 1-based 2105 self.nucmod = {} 2106# self.ecp = {atom_symbol: [[l, (r_order, expnt, c),...]]} 2107 self.ecp = {} 2108# Nuclear property. self.nucprop = {atom_symbol: {key: value}} 2109 self.nucprop = {} 2110################################################## 2111# don't modify the following private variables, they are not input options 2112 self._atm = numpy.zeros((0,6), dtype=numpy.int32) 2113 self._bas = numpy.zeros((0,8), dtype=numpy.int32) 2114 self._env = numpy.zeros(PTR_ENV_START) 2115 self._ecpbas = numpy.zeros((0,8), dtype=numpy.int32) 2116 2117 self.stdout = sys.stdout 2118 self.groupname = 'C1' 2119 self.topgroup = 'C1' 2120 self.symm_orb = None 2121 self.irrep_id = None 2122 self.irrep_name = None 2123 self._symm_orig = None 2124 self._symm_axes = None 2125 self._nelectron = None 2126 self._nao = None 2127 self._enuc = None 2128 self._atom = [] 2129 self._basis = {} 2130 self._ecp = {} 2131 self._built = False 2132 2133 # _pseudo is created to make the mol object consistenet with the mol 2134 # object converted from Cell.to_mol(). It is initialized in the 2135 # Cell.build() method only. Assigning _pseudo to mol object basically 2136 # has no effects. Mole.build() method does not have code to access the 2137 # contents of _pseudo. 2138 self._pseudo = {} 2139 2140 keys = set(('verbose', 'unit', 'cart', 'incore_anyway')) 2141 self._keys = set(self.__dict__.keys()).union(keys) 2142 self.__dict__.update(kwargs) 2143 2144 @property 2145 def natm(self): 2146 return len(self._atm) 2147 @property 2148 def nbas(self): 2149 return len(self._bas) 2150 2151 @property 2152 def nelec(self): 2153 ne = self.nelectron 2154 nalpha = (ne + self.spin) // 2 2155 nbeta = nalpha - self.spin 2156 assert(nalpha >= 0 and nbeta >= 0) 2157 if nalpha + nbeta != ne: 2158 raise RuntimeError('Electron number %d and spin %d are not consistent\n' 2159 'Note mol.spin = 2S = Nalpha - Nbeta, not 2S+1' % 2160 (ne, self.spin)) 2161 return nalpha, nbeta 2162 @nelec.setter 2163 def nelec(self, neleca_nelecb): 2164 neleca, nelecb = neleca_nelecb 2165 self._nelectron = neleca + nelecb 2166 self.spin = neleca - nelecb 2167 2168 @property 2169 def nelectron(self): 2170 if self._nelectron is None: 2171 return self.tot_electrons() 2172 else: 2173 return self._nelectron 2174 @nelectron.setter 2175 def nelectron(self, n): 2176 self._nelectron = n 2177 2178 @property 2179 def multiplicity(self): 2180 return self.spin + 1 2181 @multiplicity.setter 2182 def multiplicity(self, x): 2183 if x is None: 2184 self.spin = None 2185 else: 2186 self.spin = x - 1 2187 2188 @property 2189 def ms(self): 2190 '''Spin quantum number. multiplicity = ms*2+1''' 2191 if self.spin % 2 == 0: 2192 return self.spin // 2 2193 else: 2194 return self.spin * .5 2195 @ms.setter 2196 def ms(self, x): 2197 if x is None: 2198 self.spin = None 2199 else: 2200 self.spin = int(round(2*x, 4)) 2201 2202 def __getattr__(self, key): 2203 '''To support accessing methods (mol.HF, mol.KS, mol.CCSD, mol.CASSCF, ...) 2204 from Mole object. 2205 ''' 2206 if key[:2] == '__': # Skip Python builtins 2207 raise AttributeError('Mole object has no attribute %s' % key) 2208 elif key in ('_ipython_canary_method_should_not_exist_', 2209 '_repr_mimebundle_'): 2210 # https://github.com/mewwts/addict/issues/26 2211 # https://github.com/jupyter/notebook/issues/2014 2212 raise AttributeError 2213 2214 # Import all available modules. Some methods are registered to other 2215 # classes/modules when importing modules in __all__. 2216 from pyscf import __all__ # noqa 2217 from pyscf import scf, dft 2218 for mod in (scf, dft): 2219 method = getattr(mod, key, None) 2220 if callable(method): 2221 return method(self) 2222 2223 if 'TD' in key[:3]: 2224 if key in ('TDHF', 'TDA'): 2225 mf = scf.HF(self) 2226 else: 2227 mf = dft.KS(self) 2228 xc = key.split('TD', 1)[1] 2229 if xc in dft.XC: 2230 mf.xc = xc 2231 key = 'TDDFT' 2232 else: 2233 mf = scf.HF(self) 2234 2235 method = getattr(mf, key, None) 2236 if method is None: 2237 raise AttributeError('Mole object has no attribute %s' % key) 2238 2239 # Initialize SCF object for post-SCF methods if applicable 2240 if self.nelectron != 0: 2241 mf.run() 2242 return method 2243 2244# need "deepcopy" here because in shallow copy, _env may get new elements but 2245# with ptr_env unchanged 2246# def __copy__(self): 2247# cls = self.__class__ 2248# newmol = cls.__new__(cls) 2249# newmol = ... 2250# do not use __copy__ to aovid iteratively call copy.copy 2251 copy = copy 2252 2253 pack = pack 2254 2255 @classmethod 2256 @lib.with_doc(unpack.__doc__) 2257 def unpack(cls, moldic): 2258 return unpack(moldic) 2259 2260 @lib.with_doc(unpack.__doc__) 2261 def unpack_(self, moldic): 2262 self.__dict__.update(moldic) 2263 return self 2264 2265 dumps = dumps 2266 2267 @classmethod 2268 @lib.with_doc(loads.__doc__) 2269 def loads(cls, molstr): 2270 return loads(molstr) 2271 2272 @lib.with_doc(loads.__doc__) 2273 def loads_(self, molstr): 2274 self.__dict__.update(loads(molstr).__dict__) 2275 return self 2276 2277 def build(self, dump_input=True, parse_arg=True, 2278 verbose=None, output=None, max_memory=None, 2279 atom=None, basis=None, unit=None, nucmod=None, ecp=None, 2280 charge=None, spin=0, symmetry=None, symmetry_subgroup=None, 2281 cart=None): 2282 '''Setup moleclue and initialize some control parameters. Whenever you 2283 change the value of the attributes of :class:`Mole`, you need call 2284 this function to refresh the internal data of Mole. 2285 2286 Kwargs: 2287 dump_input : bool 2288 whether to dump the contents of input file in the output file 2289 parse_arg : bool 2290 whether to read the sys.argv and overwrite the relevant parameters 2291 verbose : int 2292 Print level. If given, overwrite :attr:`Mole.verbose` 2293 output : str or None 2294 Output file. If given, overwrite :attr:`Mole.output` 2295 max_memory : int, float 2296 Allowd memory in MB. If given, overwrite :attr:`Mole.max_memory` 2297 atom : list or str 2298 To define molecluar structure. 2299 basis : dict or str 2300 To define basis set. 2301 nucmod : dict or str 2302 Nuclear model. If given, overwrite :attr:`Mole.nucmod` 2303 charge : int 2304 Charge of molecule. It affects the electron numbers 2305 If given, overwrite :attr:`Mole.charge` 2306 spin : int 2307 2S, num. alpha electrons - num. beta electrons to control 2308 multiplicity. If setting spin = None , multiplicity will be 2309 guessed based on the neutral molecule. 2310 If given, overwrite :attr:`Mole.spin` 2311 symmetry : bool or str 2312 Whether to use symmetry. If given a string of point group 2313 name, the given point group symmetry will be used. 2314 2315 ''' 2316 if not DISABLE_GC: 2317 gc.collect() # To release circular referred objects 2318 2319 if isinstance(dump_input, (str, unicode)): 2320 sys.stderr.write('Assigning the first argument %s to mol.atom\n' % 2321 dump_input) 2322 dump_input, atom = True, dump_input 2323 2324 if verbose is not None: self.verbose = verbose 2325 if output is not None: self.output = output 2326 if max_memory is not None: self.max_memory = max_memory 2327 if atom is not None: self.atom = atom 2328 if basis is not None: self.basis = basis 2329 if unit is not None: self.unit = unit 2330 if nucmod is not None: self.nucmod = nucmod 2331 if ecp is not None: self.ecp = ecp 2332 if charge is not None: self.charge = charge 2333 if spin != 0: self.spin = spin 2334 if symmetry is not None: self.symmetry = symmetry 2335 if symmetry_subgroup is not None: self.symmetry_subgroup = symmetry_subgroup 2336 if cart is not None: self.cart = cart 2337 2338 if parse_arg: 2339 _update_from_cmdargs_(self) 2340 2341 # avoid opening output file twice 2342 if (self.output is not None 2343 # StringIO() does not have attribute 'name' 2344 and getattr(self.stdout, 'name', None) != self.output): 2345 2346 if self.verbose > logger.QUIET: 2347 if os.path.isfile(self.output): 2348 print('overwrite output file: %s' % self.output) 2349 else: 2350 print('output file: %s' % self.output) 2351 2352 if self.output == '/dev/null': 2353 self.stdout = open(os.devnull, 'w') 2354 else: 2355 self.stdout = open(self.output, 'w') 2356 2357 if self.verbose >= logger.WARN: 2358 self.check_sanity() 2359 2360 self._atom = self.format_atom(self.atom, unit=self.unit) 2361 uniq_atoms = set([a[0] for a in self._atom]) 2362 2363 if isinstance(self.basis, (str, unicode, tuple, list)): 2364 # specify global basis for whole molecule 2365 _basis = dict(((a, self.basis) for a in uniq_atoms)) 2366 elif 'default' in self.basis: 2367 default_basis = self.basis['default'] 2368 _basis = dict(((a, default_basis) for a in uniq_atoms)) 2369 _basis.update(self.basis) 2370 del(_basis['default']) 2371 else: 2372 _basis = self.basis 2373 self._basis = self.format_basis(_basis) 2374 2375# TODO: Consider ECP info in point group symmetry initialization 2376 if self.ecp: 2377 # Unless explicitly input, ECP should not be assigned to ghost atoms 2378 if isinstance(self.ecp, (str, unicode)): 2379 _ecp = dict([(a, str(self.ecp)) 2380 for a in uniq_atoms if not is_ghost_atom(a)]) 2381 elif 'default' in self.ecp: 2382 default_ecp = self.ecp['default'] 2383 _ecp = dict(((a, default_ecp) 2384 for a in uniq_atoms if not is_ghost_atom(a))) 2385 _ecp.update(self.ecp) 2386 del(_ecp['default']) 2387 else: 2388 _ecp = self.ecp 2389 self._ecp = self.format_ecp(_ecp) 2390 2391 env = self._env[:PTR_ENV_START] 2392 self._atm, self._bas, self._env = \ 2393 self.make_env(self._atom, self._basis, env, self.nucmod, 2394 self.nucprop) 2395 self._atm, self._ecpbas, self._env = \ 2396 self.make_ecp_env(self._atm, self._ecp, self._env) 2397 2398 if self.spin is None: 2399 self.spin = self.nelectron % 2 2400 else: 2401 # Access self.nelec in which the code checks whether the spin and 2402 # number of electrons are consistent. 2403 self.nelec 2404 2405 if self.symmetry: 2406 self._build_symmetry() 2407 2408 if dump_input and not self._built and self.verbose > logger.NOTE: 2409 self.dump_input() 2410 2411 logger.debug3(self, 'arg.atm = %s', str(self._atm)) 2412 logger.debug3(self, 'arg.bas = %s', str(self._bas)) 2413 logger.debug3(self, 'arg.env = %s', str(self._env)) 2414 logger.debug3(self, 'ecpbas = %s', str(self._ecpbas)) 2415 2416 self._built = True 2417 return self 2418 kernel = build 2419 2420 def _build_symmetry(self): 2421 ''' 2422 Update symmetry related attributes: topgroup, groupname, _symm_orig, 2423 _symm_axes, irrep_id, irrep_name, symm_orb 2424 ''' 2425 from pyscf import symm 2426 self.topgroup, orig, axes = symm.detect_symm(self._atom, self._basis) 2427 2428 if isinstance(self.symmetry, (str, unicode)): 2429 self.symmetry = str(symm.std_symb(self.symmetry)) 2430 try: 2431 self.groupname, axes = symm.as_subgroup(self.topgroup, axes, 2432 self.symmetry) 2433 except PointGroupSymmetryError: 2434 raise PointGroupSymmetryError( 2435 'Unable to identify input symmetry %s. Try symmetry="%s"' % 2436 (self.symmetry, self.topgroup)) 2437 else: 2438 self.groupname, axes = symm.as_subgroup(self.topgroup, axes, 2439 self.symmetry_subgroup) 2440 self._symm_orig = orig 2441 self._symm_axes = axes 2442 2443 if self.cart and self.groupname in ('Dooh', 'Coov', 'SO3'): 2444 if self.groupname == 'Coov': 2445 self.groupname, lgroup = 'C2v', self.groupname 2446 else: 2447 self.groupname, lgroup = 'D2h', self.groupname 2448 logger.warn(self, 'This version does not support symmetry %s ' 2449 'for cartesian GTO basis. Its subgroup %s is used', 2450 lgroup, self.groupname) 2451 2452 self.symm_orb, self.irrep_id = \ 2453 symm.symm_adapted_basis(self, self.groupname, orig, axes) 2454 self.irrep_name = [symm.irrep_id2name(self.groupname, ir) 2455 for ir in self.irrep_id] 2456 return self 2457 2458 @lib.with_doc(format_atom.__doc__) 2459 def format_atom(self, atom, origin=0, axes=None, unit='Ang'): 2460 return format_atom(atom, origin, axes, unit) 2461 2462 @lib.with_doc(format_basis.__doc__) 2463 def format_basis(self, basis_tab): 2464 return format_basis(basis_tab) 2465 2466 @lib.with_doc(format_ecp.__doc__) 2467 def format_ecp(self, ecp_tab): 2468 return format_ecp(ecp_tab) 2469 2470 @lib.with_doc(expand_etb.__doc__) 2471 def expand_etb(self, l, n, alpha, beta): 2472 return expand_etb(l, n, alpha, beta) 2473 2474 @lib.with_doc(expand_etbs.__doc__) 2475 def expand_etbs(self, etbs): 2476 return expand_etbs(etbs) 2477 etbs = expand_etbs 2478 2479 @lib.with_doc(make_env.__doc__) 2480 def make_env(self, atoms, basis, pre_env=[], nucmod={}, nucprop=None): 2481 if nucprop is None: 2482 nucprop = self.nucprop 2483 return make_env(atoms, basis, pre_env, nucmod, nucprop) 2484 2485 @lib.with_doc(make_atm_env.__doc__) 2486 def make_atm_env(self, atom, ptr=0, nucmod=NUC_POINT, nucprop=None): 2487 if nucprop is None: 2488 nucprop = self.nucprop.get(atom[0], {}) 2489 return make_atm_env(atom, ptr, nucmod, nucprop) 2490 2491 @lib.with_doc(make_bas_env.__doc__) 2492 def make_bas_env(self, basis_add, atom_id=0, ptr=0): 2493 return make_bas_env(basis_add, atom_id, ptr) 2494 2495 @lib.with_doc(make_ecp_env.__doc__) 2496 def make_ecp_env(self, _atm, _ecp, pre_env=[]): 2497 if _ecp: 2498 _atm, _ecpbas, _env = make_ecp_env(self, _atm, _ecp, pre_env) 2499 else: 2500 _atm, _ecpbas, _env = _atm, numpy.zeros((0,BAS_SLOTS)), pre_env 2501 return _atm, _ecpbas, _env 2502 2503 tot_electrons = tot_electrons 2504 2505 @lib.with_doc(gto_norm.__doc__) 2506 def gto_norm(self, l, expnt): 2507 return gto_norm(l, expnt) 2508 2509 2510 def dump_input(self): 2511 import __main__ 2512 if hasattr(__main__, '__file__'): 2513 try: 2514 filename = os.path.abspath(__main__.__file__) 2515 finput = open(filename, 'r') 2516 self.stdout.write('#INFO: **** input file is %s ****\n' % filename) 2517 self.stdout.write(finput.read()) 2518 self.stdout.write('#INFO: ******************** input file end ********************\n') 2519 self.stdout.write('\n') 2520 self.stdout.write('\n') 2521 finput.close() 2522 except IOError: 2523 logger.warn(self, 'input file does not exist') 2524 2525 self.stdout.write('System: %s Threads %s\n' % 2526 (str(platform.uname()), lib.num_threads())) 2527 self.stdout.write('Python %s\n' % sys.version) 2528 self.stdout.write('numpy %s scipy %s\n' % 2529 (numpy.__version__, scipy.__version__)) 2530 self.stdout.write('Date: %s\n' % time.ctime()) 2531 import pyscf 2532 self.stdout.write('PySCF version %s\n' % pyscf.__version__) 2533 info = lib.repo_info(os.path.join(__file__, '..', '..')) 2534 self.stdout.write('PySCF path %s\n' % info['path']) 2535 if 'git' in info: 2536 self.stdout.write(info['git'] + '\n') 2537 2538 self.stdout.write('\n') 2539 for key in os.environ: 2540 if 'PYSCF' in key: 2541 self.stdout.write('[ENV] %s %s\n' % (key, os.environ[key])) 2542 if self.verbose >= logger.DEBUG2: 2543 for key in dir(__config__): 2544 if key[:2] != '__': 2545 self.stdout.write('[CONFIG] %s = %s\n' % 2546 (key, getattr(__config__, key))) 2547 else: 2548 conf_file = getattr(__config__, 'conf_file', None) 2549 self.stdout.write('[CONFIG] conf_file %s\n' % conf_file) 2550 2551 self.stdout.write('[INPUT] verbose = %d\n' % self.verbose) 2552 if self.verbose >= logger.DEBUG: 2553 self.stdout.write('[INPUT] max_memory = %s \n' % self.max_memory) 2554 self.stdout.write('[INPUT] num. atoms = %d\n' % self.natm) 2555 self.stdout.write('[INPUT] num. electrons = %d\n' % self.nelectron) 2556 self.stdout.write('[INPUT] charge = %d\n' % self.charge) 2557 self.stdout.write('[INPUT] spin (= nelec alpha-beta = 2S) = %d\n' % self.spin) 2558 self.stdout.write('[INPUT] symmetry %s subgroup %s\n' % 2559 (self.symmetry, self.symmetry_subgroup)) 2560 self.stdout.write('[INPUT] Mole.unit = %s\n' % self.unit) 2561 if self.cart: 2562 self.stdout.write('[INPUT] Cartesian GTO integrals (6d 10f)\n') 2563 2564 for ia,atom in enumerate(self._atom): 2565 coorda = tuple([x * param.BOHR for x in atom[1]]) 2566 coordb = tuple([x for x in atom[1]]) 2567 self.stdout.write('[INPUT]%3d %-4s %16.12f %16.12f %16.12f AA ' 2568 '%16.12f %16.12f %16.12f Bohr\n' 2569 % ((ia+1, _symbol(atom[0])) + coorda + coordb)) 2570 if self.nucmod: 2571 if isinstance(self.nucmod, (int, str, unicode, 2572 types.FunctionType)): 2573 nucatms = [_symbol(atom[0]) for atom in self._atom] 2574 else: 2575 nucatms = self.nucmod.keys() 2576 self.stdout.write('[INPUT] Gaussian nuclear model for atoms %s\n' % 2577 nucatms) 2578 2579 if self.nucprop: 2580 self.stdout.write('[INPUT] nucprop %s\n' % self.nucprop) 2581 2582 if self.verbose >= logger.DEBUG: 2583 self.stdout.write('[INPUT] ---------------- BASIS SET ---------------- \n') 2584 self.stdout.write('[INPUT] l, kappa, [nprim/nctr], ' 2585 'expnt, c_1 c_2 ...\n') 2586 for atom, basis_set in self._basis.items(): 2587 self.stdout.write('[INPUT] %s\n' % atom) 2588 for b in basis_set: 2589 if isinstance(b[1], int): 2590 kappa = b[1] 2591 b_coeff = b[2:] 2592 else: 2593 kappa = 0 2594 b_coeff = b[1:] 2595 nprim = len(b_coeff) 2596 nctr = len(b_coeff[0])-1 2597 if nprim < nctr: 2598 logger.warn(self, 'num. primitives smaller than num. contracted basis') 2599 self.stdout.write('[INPUT] %d %2d [%-5d/%-4d] ' 2600 % (b[0], kappa, nprim, nctr)) 2601 for k, x in enumerate(b_coeff): 2602 if k == 0: 2603 self.stdout.write('%-15.12g ' % x[0]) 2604 else: 2605 self.stdout.write(' '*32+'%-15.12g ' % x[0]) 2606 for c in x[1:]: 2607 self.stdout.write(' %4.12g' % c) 2608 self.stdout.write('\n') 2609 2610 if self.verbose >= logger.INFO: 2611 self.stdout.write('\n') 2612 logger.info(self, 'nuclear repulsion = %.15g', self.energy_nuc()) 2613 if self.symmetry: 2614 if self.topgroup == self.groupname: 2615 logger.info(self, 'point group symmetry = %s', self.topgroup) 2616 else: 2617 logger.info(self, 'point group symmetry = %s, use subgroup %s', 2618 self.topgroup, self.groupname) 2619 for ir in range(self.symm_orb.__len__()): 2620 logger.info(self, 'num. orbitals of irrep %s = %d', 2621 self.irrep_name[ir], self.symm_orb[ir].shape[1]) 2622 logger.info(self, 'number of shells = %d', self.nbas) 2623 logger.info(self, 'number of NR pGTOs = %d', self.npgto_nr()) 2624 logger.info(self, 'number of NR cGTOs = %d', self.nao_nr()) 2625 logger.info(self, 'basis = %s', self.basis) 2626 logger.info(self, 'ecp = %s', self.ecp) 2627 if self.verbose >= logger.DEBUG2: 2628 for i in range(len(self._bas)): 2629 exps = self.bas_exp(i) 2630 logger.debug1(self, 'bas %d, expnt(s) = %s', i, str(exps)) 2631 2632 logger.info(self, 'CPU time: %12.2f', logger.process_clock()) 2633 return self 2634 2635 def set_common_origin(self, coord): 2636 '''Update common origin for integrals of dipole, rxp etc. 2637 **Note** the unit of the coordinates needs to be Bohr 2638 2639 Examples: 2640 2641 >>> mol.set_common_origin(0) 2642 >>> mol.set_common_origin((1,0,0)) 2643 ''' 2644 self._env[PTR_COMMON_ORIG:PTR_COMMON_ORIG+3] = coord 2645 return self 2646 set_common_orig = set_common_origin 2647 set_common_orig_ = set_common_orig # for backward compatibility 2648 set_common_origin_ = set_common_orig # for backward compatibility 2649 2650 def with_common_origin(self, coord): 2651 '''Retuen a temporary mol context which has the rquired common origin. 2652 The required common origin has no effects out of the temporary context. 2653 See also :func:`mol.set_common_origin` 2654 2655 Examples: 2656 2657 >>> with mol.with_common_origin((1,0,0)): 2658 ... mol.intor('int1e_r', comp=3) 2659 ''' 2660 coord0 = self._env[PTR_COMMON_ORIG:PTR_COMMON_ORIG+3].copy() 2661 return self._TemporaryMoleContext(self.set_common_origin, (coord,), (coord0,)) 2662 with_common_orig = with_common_origin 2663 2664 def set_rinv_origin(self, coord): 2665 r'''Update origin for operator :math:`\frac{1}{|r-R_O|}`. 2666 **Note** the unit is Bohr 2667 2668 Examples: 2669 2670 >>> mol.set_rinv_origin(0) 2671 >>> mol.set_rinv_origin((0,1,0)) 2672 ''' 2673 self._env[PTR_RINV_ORIG:PTR_RINV_ORIG+3] = coord[:3] 2674 return self 2675 set_rinv_orig = set_rinv_origin 2676 set_rinv_orig_ = set_rinv_orig # for backward compatibility 2677 set_rinv_origin_ = set_rinv_orig # for backward compatibility 2678 2679 def with_rinv_origin(self, coord): 2680 '''Retuen a temporary mol context which has the rquired origin of 1/r 2681 operator. The required origin has no effects out of the temporary 2682 context. See also :func:`mol.set_rinv_origin` 2683 2684 Examples: 2685 2686 >>> with mol.with_rinv_origin((1,0,0)): 2687 ... mol.intor('int1e_rinv') 2688 ''' 2689 coord0 = self._env[PTR_RINV_ORIG:PTR_RINV_ORIG+3].copy() 2690 return self._TemporaryMoleContext(self.set_rinv_origin, (coord,), (coord0,)) 2691 with_rinv_orig = with_rinv_origin 2692 2693 def set_range_coulomb(self, omega): 2694 '''Switch on range-separated Coulomb operator for **all** 2e integrals 2695 2696 Args: 2697 omega : double 2698 2699 | = 0 : Regular electron repulsion integral 2700 | > 0 : Long-range operator erf(omega r12) / r12 2701 | < 0 : Short-range operator erfc(omega r12) /r12 2702 ''' 2703 if omega is None: 2704 self._env[PTR_RANGE_OMEGA] = 0 2705 else: 2706 self._env[PTR_RANGE_OMEGA] = omega 2707 set_range_coulomb_ = set_range_coulomb # for backward compatibility 2708 2709 @property 2710 def omega(self): 2711 return self._env[PTR_RANGE_OMEGA] 2712 omega = omega.setter(set_range_coulomb) 2713 2714 def with_range_coulomb(self, omega): 2715 '''Retuen a temporary mol context which sets the required parameter 2716 omega for range-separated Coulomb operator. 2717 If omega = None, return the context for regular Coulomb integrals. 2718 See also :func:`mol.set_range_coulomb` 2719 2720 Examples: 2721 2722 >>> with mol.with_range_coulomb(omega=1.5): 2723 ... mol.intor('int2e') 2724 ''' 2725 omega0 = self._env[PTR_RANGE_OMEGA].copy() 2726 return self._TemporaryMoleContext(self.set_range_coulomb, (omega,), (omega0,)) 2727 2728 def with_long_range_coulomb(self, omega): 2729 '''Retuen a temporary mol context for long-range part of 2730 range-separated Coulomb operator. 2731 ''' 2732 return self.with_range_coulomb(abs(omega)) 2733 2734 def with_short_range_coulomb(self, omega): 2735 '''Retuen a temporary mol context for short-range part of 2736 range-separated Coulomb operator. 2737 ''' 2738 return self.with_range_coulomb(-abs(omega)) 2739 2740 def set_f12_zeta(self, zeta): 2741 '''Set zeta for YP exp(-zeta r12)/r12 or STG exp(-zeta r12) type integrals 2742 ''' 2743 self._env[PTR_F12_ZETA] = zeta 2744 2745 def set_nuc_mod(self, atm_id, zeta): 2746 '''Change the nuclear charge distribution of the given atom ID. The charge 2747 distribution is defined as: rho(r) = nuc_charge * Norm * exp(-zeta * r^2). 2748 This function can **only** be called after .build() method is executed. 2749 2750 Examples: 2751 2752 >>> for ia in range(mol.natm): 2753 ... zeta = gto.filatov_nuc_mod(mol.atom_charge(ia)) 2754 ... mol.set_nuc_mod(ia, zeta) 2755 ''' 2756 ptr = self._atm[atm_id,PTR_ZETA] 2757 self._env[ptr] = zeta 2758 if zeta == 0: 2759 self._atm[atm_id,NUC_MOD_OF] = NUC_POINT 2760 else: 2761 self._atm[atm_id,NUC_MOD_OF] = NUC_GAUSS 2762 return self 2763 set_nuc_mod_ = set_nuc_mod # for backward compatibility 2764 2765 def set_rinv_zeta(self, zeta): 2766 '''Assume the charge distribution on the "rinv_origin". zeta is the parameter 2767 to control the charge distribution: rho(r) = Norm * exp(-zeta * r^2). 2768 **Be careful** when call this function. It affects the behavior of 2769 int1e_rinv_* functions. Make sure to set it back to 0 after using it! 2770 ''' 2771 self._env[PTR_RINV_ZETA] = zeta 2772 return self 2773 set_rinv_zeta_ = set_rinv_zeta # for backward compatibility 2774 2775 def with_rinv_zeta(self, zeta): 2776 '''Retuen a temporary mol context which has the rquired Gaussian charge 2777 distribution placed at "rinv_origin": rho(r) = Norm * exp(-zeta * r^2). 2778 See also :func:`mol.set_rinv_zeta` 2779 2780 Examples: 2781 2782 >>> with mol.with_rinv_zeta(zeta=1.5), mol.with_rinv_origin((1.,0,0)): 2783 ... mol.intor('int1e_rinv') 2784 ''' 2785 zeta0 = self._env[PTR_RINV_ZETA].copy() 2786 return self._TemporaryMoleContext(self.set_rinv_zeta, (zeta,), (zeta0,)) 2787 2788 def with_rinv_at_nucleus(self, atm_id): 2789 '''Retuen a temporary mol context in which the rinv operator (1/r) is 2790 treated like the Coulomb potential of a Gaussian charge distribution 2791 rho(r) = Norm * exp(-zeta * r^2) at the place of the input atm_id. 2792 2793 Examples: 2794 2795 >>> with mol.with_rinv_at_nucleus(3): 2796 ... mol.intor('int1e_rinv') 2797 ''' 2798 zeta = self._env[self._atm[atm_id,PTR_ZETA]] 2799 rinv = self.atom_coord(atm_id) 2800 if zeta == 0: 2801 self._env[AS_RINV_ORIG_ATOM] = atm_id # required by ecp gradients 2802 return self.with_rinv_origin(rinv) 2803 else: 2804 self._env[AS_RINV_ORIG_ATOM] = atm_id # required by ecp gradients 2805 rinv0 = self._env[PTR_RINV_ORIG:PTR_RINV_ORIG+3].copy() 2806 zeta0 = self._env[PTR_RINV_ZETA].copy() 2807 2808 def set_rinv(z, r): 2809 self._env[PTR_RINV_ZETA] = z 2810 self._env[PTR_RINV_ORIG:PTR_RINV_ORIG+3] = r 2811 return self._TemporaryMoleContext(set_rinv, (zeta,rinv), (zeta0,rinv0)) 2812 with_rinv_as_nucleus = with_rinv_at_nucleus # For backward compatibility 2813 2814 def with_integral_screen(self, threshold): 2815 '''Retuen a temporary mol context which has the rquired integral 2816 screen threshold 2817 ''' 2818 expcutoff0 = self._env[PTR_EXPCUTOFF] 2819 expcutoff = abs(numpy.log(threshold)) 2820 def set_cutoff(cut): 2821 self._env[PTR_EXPCUTOFF] = cut 2822 return self._TemporaryMoleContext(set_cutoff, (expcutoff,), (expcutoff0,)) 2823 2824 def set_geom_(self, atoms_or_coords, unit=None, symmetry=None, 2825 inplace=True): 2826 '''Update geometry 2827 ''' 2828 import copy 2829 if inplace: 2830 mol = self 2831 else: 2832 mol = copy.copy(self) 2833 mol._env = mol._env.copy() 2834 if unit is None: 2835 unit = mol.unit 2836 else: 2837 mol.unit = unit 2838 if symmetry is None: 2839 symmetry = mol.symmetry 2840 2841 if isinstance(atoms_or_coords, numpy.ndarray): 2842 mol.atom = list(zip([x[0] for x in mol._atom], 2843 atoms_or_coords.tolist())) 2844 else: 2845 mol.atom = atoms_or_coords 2846 2847 if isinstance(atoms_or_coords, numpy.ndarray) and not symmetry: 2848 mol._atom = mol.atom 2849 2850 if isinstance(unit, (str, unicode)): 2851 if unit.upper().startswith(('B', 'AU')): 2852 unit = 1. 2853 else: #unit[:3].upper() == 'ANG': 2854 unit = 1./param.BOHR 2855 else: 2856 unit = 1./unit 2857 ptr = mol._atm[:,PTR_COORD] 2858 mol._env[ptr+0] = unit * atoms_or_coords[:,0] 2859 mol._env[ptr+1] = unit * atoms_or_coords[:,1] 2860 mol._env[ptr+2] = unit * atoms_or_coords[:,2] 2861 else: 2862 mol.symmetry = symmetry 2863 mol.build(False, False) 2864 2865 if mol.verbose >= logger.INFO: 2866 logger.info(mol, 'New geometry') 2867 for ia, atom in enumerate(mol._atom): 2868 coorda = tuple([x * param.BOHR for x in atom[1]]) 2869 coordb = tuple([x for x in atom[1]]) 2870 coords = coorda + coordb 2871 logger.info(mol, ' %3d %-4s %16.12f %16.12f %16.12f AA ' 2872 '%16.12f %16.12f %16.12f Bohr\n', 2873 ia+1, mol.atom_symbol(ia), *coords) 2874 return mol 2875 2876 def update(self, chkfile): 2877 return self.update_from_chk(chkfile) 2878 def update_from_chk(self, chkfile): 2879 with h5py.File(chkfile, 'r') as fh5: 2880 mol = loads(fh5['mol'][()]) 2881 self.__dict__.update(mol.__dict__) 2882 return self 2883 2884 def has_ecp(self): 2885 '''Whether pseudo potential is used in the system.''' 2886 return len(self._ecpbas) > 0 or self._pseudo 2887 2888 def has_ecp_soc(self): 2889 '''Whether spin-orbit coupling is enabled in ECP.''' 2890 return (len(self._ecpbas) > 0 and 2891 numpy.any(self._ecpbas[:,SO_TYPE_OF] == 1)) 2892 2893 2894####################################################### 2895#NOTE: atm_id or bas_id start from 0 2896 def atom_symbol(self, atm_id): 2897 r'''For the given atom id, return the input symbol (without striping special characters) 2898 2899 Args: 2900 atm_id : int 2901 0-based 2902 2903 Examples: 2904 2905 >>> mol.build(atom='H^2 0 0 0; H 0 0 1.1') 2906 >>> mol.atom_symbol(0) 2907 H^2 2908 ''' 2909 return _symbol(self._atom[atm_id][0]) 2910 2911 def atom_pure_symbol(self, atm_id): 2912 r'''For the given atom id, return the standard symbol (striping special characters) 2913 2914 Args: 2915 atm_id : int 2916 0-based 2917 2918 Examples: 2919 2920 >>> mol.build(atom='H^2 0 0 0; H 0 0 1.1') 2921 >>> mol.atom_symbol(0) 2922 H 2923 ''' 2924 return _std_symbol(self._atom[atm_id][0]) 2925 2926 @property 2927 def elements(self): 2928 '''A list of elements in the molecule''' 2929 return [self.atom_pure_symbol(i) for i in range(self.natm)] 2930 2931 def atom_charge(self, atm_id): 2932 r'''Nuclear effective charge of the given atom id 2933 Note "atom_charge /= charge(atom_symbol)" when ECP is enabled. 2934 Number of electrons screened by ECP can be obtained by charge(atom_symbol)-atom_charge 2935 2936 Args: 2937 atm_id : int 2938 0-based 2939 2940 Examples: 2941 2942 >>> mol.build(atom='H 0 0 0; Cl 0 0 1.1') 2943 >>> mol.atom_charge(1) 2944 17 2945 ''' 2946 if self._atm[atm_id,NUC_MOD_OF] != NUC_FRAC_CHARGE: 2947 # regular QM atoms 2948 return self._atm[atm_id,CHARGE_OF] 2949 else: 2950 # MM atoms with fractional charges 2951 return self._env[self._atm[atm_id,PTR_FRAC_CHARGE]] 2952 2953 def atom_charges(self): 2954 '''np.asarray([mol.atom_charge(i) for i in range(mol.natm)])''' 2955 z = self._atm[:,CHARGE_OF] 2956 if numpy.any(self._atm[:,NUC_MOD_OF] == NUC_FRAC_CHARGE): 2957 # Create the integer nuclear charges first then replace the MM 2958 # particles with the MM charges that saved in _env[PTR_FRAC_CHARGE] 2959 z = numpy.array(z, dtype=numpy.double) 2960 idx = self._atm[:,NUC_MOD_OF] == NUC_FRAC_CHARGE 2961 # MM fractional charges can be positive or negative 2962 z[idx] = self._env[self._atm[idx,PTR_FRAC_CHARGE]] 2963 return z 2964 2965 def atom_nelec_core(self, atm_id): 2966 '''Number of core electrons for pseudo potential. 2967 ''' 2968 return charge(self.atom_symbol(atm_id)) - self.atom_charge(atm_id) 2969 2970 def atom_coord(self, atm_id, unit='Bohr'): 2971 r'''Coordinates (ndarray) of the given atom id 2972 2973 Args: 2974 atm_id : int 2975 0-based 2976 2977 Examples: 2978 2979 >>> mol.build(atom='H 0 0 0; Cl 0 0 1.1') 2980 >>> mol.atom_coord(1) 2981 [ 0. 0. 2.07869874] 2982 ''' 2983 ptr = self._atm[atm_id,PTR_COORD] 2984 if unit[:3].upper() == 'ANG': 2985 return self._env[ptr:ptr+3] * param.BOHR 2986 else: 2987 return self._env[ptr:ptr+3].copy() 2988 2989 def atom_coords(self, unit='Bohr'): 2990 '''np.asarray([mol.atom_coords(i) for i in range(mol.natm)])''' 2991 ptr = self._atm[:,PTR_COORD] 2992 c = self._env[numpy.vstack((ptr,ptr+1,ptr+2)).T].copy() 2993 if unit[:3].upper() == 'ANG': 2994 c *= param.BOHR 2995 return c 2996 2997 atom_mass_list = atom_mass_list 2998 2999 def atom_nshells(self, atm_id): 3000 r'''Number of basis/shells of the given atom 3001 3002 Args: 3003 atm_id : int 3004 0-based 3005 3006 Examples: 3007 3008 >>> mol.build(atom='H 0 0 0; Cl 0 0 1.1') 3009 >>> mol.atom_nshells(1) 3010 5 3011 ''' 3012 return (self._bas[:,ATOM_OF] == atm_id).sum() 3013 3014 def atom_shell_ids(self, atm_id): 3015 r'''A list of the shell-ids of the given atom 3016 3017 Args: 3018 atm_id : int 3019 0-based 3020 3021 Examples: 3022 3023 >>> mol.build(atom='H 0 0 0; Cl 0 0 1.1', basis='cc-pvdz') 3024 >>> mol.atom_shell_ids(1) 3025 [3, 4, 5, 6, 7] 3026 ''' 3027 return numpy.where(self._bas[:,ATOM_OF] == atm_id)[0] 3028 3029 def bas_coord(self, bas_id): 3030 r'''Coordinates (ndarray) associated with the given basis id 3031 3032 Args: 3033 bas_id : int 3034 0-based 3035 3036 Examples: 3037 3038 >>> mol.build(atom='H 0 0 0; Cl 0 0 1.1') 3039 >>> mol.bas_coord(1) 3040 [ 0. 0. 2.07869874] 3041 ''' 3042 atm_id = self.bas_atom(bas_id) 3043 ptr = self._atm[atm_id,PTR_COORD] 3044 return self._env[ptr:ptr+3].copy() 3045 3046 def bas_atom(self, bas_id): 3047 r'''The atom (0-based id) that the given basis sits on 3048 3049 Args: 3050 bas_id : int 3051 0-based 3052 3053 Examples: 3054 3055 >>> mol.build(atom='H 0 0 0; Cl 0 0 1.1', basis='cc-pvdz') 3056 >>> mol.bas_atom(7) 3057 1 3058 ''' 3059 return self._bas[bas_id,ATOM_OF].copy() 3060 3061 def bas_angular(self, bas_id): 3062 r'''The angular momentum associated with the given basis 3063 3064 Args: 3065 bas_id : int 3066 0-based 3067 3068 Examples: 3069 3070 >>> mol.build(atom='H 0 0 0; Cl 0 0 1.1', basis='cc-pvdz') 3071 >>> mol.bas_atom(7) 3072 2 3073 ''' 3074 return self._bas[bas_id,ANG_OF].copy() 3075 3076 def bas_nctr(self, bas_id): 3077 r'''The number of contracted GTOs for the given shell 3078 3079 Args: 3080 bas_id : int 3081 0-based 3082 3083 Examples: 3084 3085 >>> mol.build(atom='H 0 0 0; Cl 0 0 1.1', basis='cc-pvdz') 3086 >>> mol.bas_atom(3) 3087 3 3088 ''' 3089 return self._bas[bas_id,NCTR_OF].copy() 3090 3091 def bas_nprim(self, bas_id): 3092 r'''The number of primitive GTOs for the given shell 3093 3094 Args: 3095 bas_id : int 3096 0-based 3097 3098 Examples: 3099 3100 >>> mol.build(atom='H 0 0 0; Cl 0 0 1.1', basis='cc-pvdz') 3101 >>> mol.bas_atom(3) 3102 11 3103 ''' 3104 return self._bas[bas_id,NPRIM_OF].copy() 3105 3106 def bas_kappa(self, bas_id): 3107 r'''Kappa (if l < j, -l-1, else l) of the given shell 3108 3109 Args: 3110 bas_id : int 3111 0-based 3112 3113 Examples: 3114 3115 >>> mol.build(atom='H 0 0 0; Cl 0 0 1.1', basis='cc-pvdz') 3116 >>> mol.bas_kappa(3) 3117 0 3118 ''' 3119 return self._bas[bas_id,KAPPA_OF].copy() 3120 3121 def bas_exp(self, bas_id): 3122 r'''exponents (ndarray) of the given shell 3123 3124 Args: 3125 bas_id : int 3126 0-based 3127 3128 Examples: 3129 3130 >>> mol.build(atom='H 0 0 0; Cl 0 0 1.1', basis='cc-pvdz') 3131 >>> mol.bas_kappa(0) 3132 [ 13.01 1.962 0.4446] 3133 ''' 3134 nprim = self.bas_nprim(bas_id) 3135 ptr = self._bas[bas_id,PTR_EXP] 3136 return self._env[ptr:ptr+nprim].copy() 3137 3138 def _libcint_ctr_coeff(self, bas_id): 3139 nprim = self.bas_nprim(bas_id) 3140 nctr = self.bas_nctr(bas_id) 3141 ptr = self._bas[bas_id,PTR_COEFF] 3142 return self._env[ptr:ptr+nprim*nctr].reshape(nctr,nprim).T 3143 3144 def bas_ctr_coeff(self, bas_id): 3145 r'''Contract coefficients (ndarray) of the given shell 3146 3147 Args: 3148 bas_id : int 3149 0-based 3150 3151 Examples: 3152 3153 >>> mol.M(atom='H 0 0 0; Cl 0 0 1.1', basis='cc-pvdz') 3154 >>> mol.bas_ctr_coeff(0) 3155 [[ 10.03400444] 3156 [ 4.1188704 ] 3157 [ 1.53971186]] 3158 ''' 3159 l = self.bas_angular(bas_id) 3160 es = self.bas_exp(bas_id) 3161 cs = self._libcint_ctr_coeff(bas_id) 3162 cs = numpy.einsum('pi,p->pi', cs, 1/gto_norm(l, es)) 3163 return cs 3164 3165 def bas_len_spinor(self, bas_id): 3166 '''The number of spinor associated with given basis 3167 If kappa is 0, return 4l+2 3168 ''' 3169 l = self.bas_angular(bas_id) 3170 k = self.bas_kappa(bas_id) 3171 return len_spinor(l, k) 3172 3173 def bas_len_cart(self, bas_id): 3174 '''The number of Cartesian function associated with given basis 3175 ''' 3176 return len_cart(self._bas[bas_id,ANG_OF]) 3177 3178 3179 npgto_nr = npgto_nr 3180 3181 nao_nr = nao_nr 3182 nao_2c = nao_2c 3183 nao_cart = nao_cart 3184 3185 nao_nr_range = nao_nr_range 3186 nao_2c_range = nao_2c_range 3187 3188 ao_loc_nr = ao_loc_nr 3189 ao_loc_2c = ao_loc_2c 3190 3191 @property 3192 def nao(self): 3193 if self._nao is None: 3194 return self.nao_nr() 3195 else: 3196 return self._nao 3197 @nao.setter 3198 def nao(self, x): 3199 self._nao = x 3200 3201 ao_loc = property(ao_loc_nr) 3202 3203 tmap = time_reversal_map = time_reversal_map 3204 3205 inertia_moment = inertia_moment 3206 3207 tostring = tostring 3208 tofile = tofile 3209 3210 def fromstring(self, string, format='xyz'): 3211 '''Update the Mole object based on the input geometry string''' 3212 atom = self.format_atom(fromstring(string, format), unit=1) 3213 self.set_geom_(atom, unit='Angstrom', inplace=True) 3214 if format == 'sdf' and 'M CHG' in string: 3215 raise NotImplementedError 3216 #FIXME self.charge = 0 3217 return self 3218 3219 def fromfile(self, filename, format=None): 3220 '''Update the Mole object based on the input geometry file''' 3221 atom = self.format_atom(fromfile(filename, format), unit=1) 3222 self.set_geom_(atom, unit='Angstrom', inplace=True) 3223 if format == 'sdf': 3224 raise NotImplementedError 3225 return self 3226 3227 def intor(self, intor, comp=None, hermi=0, aosym='s1', out=None, 3228 shls_slice=None, grids=None): 3229 '''Integral generator. 3230 3231 Args: 3232 intor : str 3233 Name of the 1e or 2e AO integrals. Ref to :func:`getints` for the 3234 complete list of available 1-electron integral names 3235 3236 Kwargs: 3237 comp : int 3238 Components of the integrals, e.g. int1e_ipovlp_sph has 3 components. 3239 hermi : int 3240 Symmetry of the integrals 3241 3242 | 0 : no symmetry assumed (default) 3243 | 1 : hermitian 3244 | 2 : anti-hermitian 3245 3246 grids : ndarray 3247 Coordinates of grids for the int1e_grids integrals 3248 3249 Returns: 3250 ndarray of 1-electron integrals, can be either 2-dim or 3-dim, depending on comp 3251 3252 Examples: 3253 3254 >>> mol.build(atom='H 0 0 0; H 0 0 1.1', basis='sto-3g') 3255 >>> mol.intor('int1e_ipnuc_sph', comp=3) # <nabla i | V_nuc | j> 3256 [[[ 0. 0. ] 3257 [ 0. 0. ]] 3258 [[ 0. 0. ] 3259 [ 0. 0. ]] 3260 [[ 0.10289944 0.48176097] 3261 [-0.48176097 -0.10289944]]] 3262 >>> mol.intor('int1e_nuc_spinor') 3263 [[-1.69771092+0.j 0.00000000+0.j -0.67146312+0.j 0.00000000+0.j] 3264 [ 0.00000000+0.j -1.69771092+0.j 0.00000000+0.j -0.67146312+0.j] 3265 [-0.67146312+0.j 0.00000000+0.j -1.69771092+0.j 0.00000000+0.j] 3266 [ 0.00000000+0.j -0.67146312+0.j 0.00000000+0.j -1.69771092+0.j]] 3267 ''' 3268 if not self._built: 3269 logger.warn(self, 'Warning: intor envs of %s not initialized.', self) 3270 # FIXME: Whether to check _built and call build? ._bas and .basis 3271 # may not be consistent. calling .build() may leads to wrong intor env. 3272 #self.build(False, False) 3273 intor = self._add_suffix(intor) 3274 bas = self._bas 3275 env = self._env 3276 if 'ECP' in intor: 3277 assert(self._ecp is not None) 3278 bas = numpy.vstack((self._bas, self._ecpbas)) 3279 env[AS_ECPBAS_OFFSET] = len(self._bas) 3280 env[AS_NECPBAS] = len(self._ecpbas) 3281 if shls_slice is None: 3282 shls_slice = (0, self.nbas, 0, self.nbas) 3283 elif '_grids' in intor: 3284 assert grids is not None 3285 env = numpy.append(env, grids.ravel()) 3286 env[NGRIDS] = grids.shape[0] 3287 env[PTR_GRIDS] = env.size - grids.size 3288 return moleintor.getints(intor, self._atm, bas, env, 3289 shls_slice, comp, hermi, aosym, out=out) 3290 3291 def _add_suffix(self, intor, cart=None): 3292 if not (intor[:4] == 'cint' or 3293 intor.endswith(('_sph', '_cart', '_spinor', '_ssc'))): 3294 if cart is None: 3295 cart = self.cart 3296 if cart: 3297 intor = intor + '_cart' 3298 else: 3299 intor = intor + '_sph' 3300 return intor 3301 3302 def intor_symmetric(self, intor, comp=None, grids=None): 3303 '''One-electron integral generator. The integrals are assumed to be hermitian 3304 3305 Args: 3306 intor : str 3307 Name of the 1-electron integral. Ref to :func:`getints` for the 3308 complete list of available 1-electron integral names 3309 3310 Kwargs: 3311 comp : int 3312 Components of the integrals, e.g. int1e_ipovlp_sph has 3 components. 3313 grids : ndarray 3314 Coordinates of grids for the int1e_grids integrals 3315 3316 Returns: 3317 ndarray of 1-electron integrals, can be either 2-dim or 3-dim, depending on comp 3318 3319 Examples: 3320 3321 >>> mol.build(atom='H 0 0 0; H 0 0 1.1', basis='sto-3g') 3322 >>> mol.intor_symmetric('int1e_nuc_spinor') 3323 [[-1.69771092+0.j 0.00000000+0.j -0.67146312+0.j 0.00000000+0.j] 3324 [ 0.00000000+0.j -1.69771092+0.j 0.00000000+0.j -0.67146312+0.j] 3325 [-0.67146312+0.j 0.00000000+0.j -1.69771092+0.j 0.00000000+0.j] 3326 [ 0.00000000+0.j -0.67146312+0.j 0.00000000+0.j -1.69771092+0.j]] 3327 ''' 3328 return self.intor(intor, comp, 1, aosym='s4', grids=grids) 3329 3330 def intor_asymmetric(self, intor, comp=None, grids=None): 3331 '''One-electron integral generator. The integrals are assumed to be anti-hermitian 3332 3333 Args: 3334 intor : str 3335 Name of the 1-electron integral. Ref to :func:`getints` for the 3336 complete list of available 1-electron integral names 3337 3338 Kwargs: 3339 comp : int 3340 Components of the integrals, e.g. int1e_ipovlp has 3 components. 3341 grids : ndarray 3342 Coordinates of grids for the int1e_grids integrals 3343 3344 Returns: 3345 ndarray of 1-electron integrals, can be either 2-dim or 3-dim, depending on comp 3346 3347 Examples: 3348 3349 >>> mol.build(atom='H 0 0 0; H 0 0 1.1', basis='sto-3g') 3350 >>> mol.intor_asymmetric('int1e_nuc_spinor') 3351 [[-1.69771092+0.j 0.00000000+0.j 0.67146312+0.j 0.00000000+0.j] 3352 [ 0.00000000+0.j -1.69771092+0.j 0.00000000+0.j 0.67146312+0.j] 3353 [-0.67146312+0.j 0.00000000+0.j -1.69771092+0.j 0.00000000+0.j] 3354 [ 0.00000000+0.j -0.67146312+0.j 0.00000000+0.j -1.69771092+0.j]] 3355 ''' 3356 return self.intor(intor, comp, 2, aosym='a4', grids=grids) 3357 3358 @lib.with_doc(moleintor.getints_by_shell.__doc__) 3359 def intor_by_shell(self, intor, shells, comp=None, grids=None): 3360 intor = self._add_suffix(intor) 3361 if 'ECP' in intor: 3362 assert(self._ecp is not None) 3363 bas = numpy.vstack((self._bas, self._ecpbas)) 3364 self._env[AS_ECPBAS_OFFSET] = len(self._bas) 3365 self._env[AS_NECPBAS] = len(self._ecpbas) 3366 else: 3367 bas = self._bas 3368 return moleintor.getints_by_shell(intor, shells, self._atm, bas, 3369 self._env, comp) 3370 3371 eval_ao = eval_gto = eval_gto 3372 3373 energy_nuc = energy_nuc 3374 def get_enuc(self): 3375 return self.energy_nuc() 3376 3377 sph_labels = spheric_labels = sph_labels 3378 cart_labels = cart_labels 3379 ao_labels = ao_labels 3380 3381 spinor_labels = spinor_labels 3382 3383 search_ao_label = search_ao_label 3384 3385 def search_shell_id(self, atm_id, l): 3386 return search_shell_id(self, atm_id, l) 3387 3388 search_ao_nr = search_ao_nr 3389 search_ao_r = search_ao_r 3390 3391 aoslice_by_atom = aoslice_nr_by_atom = offset_ao_by_atom = offset_nr_by_atom = aoslice_by_atom 3392 aoslice_2c_by_atom = offset_2c_by_atom = offset_2c_by_atom 3393 3394 condense_to_shell = condense_to_shell 3395 3396 to_uncontracted_cartesian_basis = to_uncontracted_cartesian_basis 3397 3398 __add__ = conc_mol 3399 3400 def cart2sph_coeff(self, normalized='sp'): 3401 '''Transformation matrix that transforms Cartesian GTOs to spherical 3402 GTOs for all basis functions 3403 3404 Kwargs: 3405 normalized : string or boolean 3406 How the Cartesian GTOs are normalized. Except s and p functions, 3407 Cartesian GTOs do not have the universal normalization coefficients 3408 for the different components of the same shell. The value of this 3409 argument can be one of 'sp', 'all', None. 'sp' means the Cartesian s 3410 and p basis are normalized. 'all' means all Cartesian functions are 3411 normalized. None means none of the Cartesian functions are normalized. 3412 The default value 'sp' is the convention used by libcint library. 3413 3414 Examples: 3415 3416 >>> mol = gto.M(atom='H 0 0 0; F 0 0 1', basis='ccpvtz') 3417 >>> c = mol.cart2sph_coeff() 3418 >>> s0 = mol.intor('int1e_ovlp_sph') 3419 >>> s1 = c.T.dot(mol.intor('int1e_ovlp_cart')).dot(c) 3420 >>> print(abs(s1-s0).sum()) 3421 >>> 4.58676826646e-15 3422 ''' 3423 c2s_l = [cart2sph(l, normalized=normalized) for l in range(12)] 3424 c2s = [] 3425 for ib in range(self.nbas): 3426 l = self.bas_angular(ib) 3427 for n in range(self.bas_nctr(ib)): 3428 c2s.append(c2s_l[l]) 3429 return scipy.linalg.block_diag(*c2s) 3430 3431 def sph2spinor_coeff(self): 3432 '''Transformation matrix that transforms real-spherical GTOs to spinor 3433 GTOs for all basis functions 3434 3435 Examples: 3436 3437 >>> mol = gto.M(atom='H 0 0 0; F 0 0 1', basis='ccpvtz') 3438 >>> ca, cb = mol.sph2spinor_coeff() 3439 >>> s0 = mol.intor('int1e_ovlp_spinor') 3440 >>> s1 = ca.conj().T.dot(mol.intor('int1e_ovlp_sph')).dot(ca) 3441 >>> s1+= cb.conj().T.dot(mol.intor('int1e_ovlp_sph')).dot(cb) 3442 >>> print(abs(s1-s0).max()) 3443 >>> 6.66133814775e-16 3444 ''' 3445 from pyscf.symm import sph 3446 return sph.sph2spinor_coeff(self) 3447 3448 3449 def apply(self, fn, *args, **kwargs): 3450 if callable(fn): 3451 return lib.StreamObject.apply(self, fn, *args, **kwargs) 3452 elif isinstance(fn, (str, unicode)): 3453 method = getattr(self, fn.upper()) 3454 return method(*args, **kwargs) 3455 else: 3456 raise TypeError('First argument of .apply method must be a ' 3457 'function/class or a name (string) of a method.') 3458 3459 def ao2mo(self, mo_coeffs, erifile=None, dataname='eri_mo', intor='int2e', 3460 **kwargs): 3461 '''Integral transformation for arbitrary orbitals and arbitrary 3462 integrals. See more detalied documentation in func:`ao2mo.kernel`. 3463 3464 Args: 3465 mo_coeffs (an np array or a list of arrays) : A matrix of orbital 3466 coefficients if it is a numpy ndarray, or four sets of orbital 3467 coefficients, corresponding to the four indices of (ij|kl). 3468 3469 Kwargs: 3470 erifile (str or h5py File or h5py Group object) : The file/object 3471 to store the transformed integrals. If not given, the return 3472 value is an array (in memory) of the transformed integrals. 3473 dataname : str 3474 *Note* this argument is effective if erifile is given. 3475 The dataset name in the erifile (ref the hierarchy of HDF5 format 3476 http://www.hdfgroup.org/HDF5/doc1.6/UG/09_Groups.html). By assigning 3477 different dataname, the existed integral file can be reused. If 3478 the erifile contains the specified dataname, the old integrals 3479 will be replaced by the new one under the key dataname. 3480 intor (str) : integral name Name of the 2-electron integral. Ref 3481 to :func:`getints_by_shell` 3482 for the complete list of available 2-electron integral names 3483 3484 Returns: 3485 An array of transformed integrals if erifile is not given. 3486 Otherwise, return the file/fileobject if erifile is assigned. 3487 3488 3489 Examples: 3490 3491 >>> import pyscf 3492 >>> mol = pyscf.M(atom='O 0 0 0; H 0 1 0; H 0 0 1', basis='sto3g') 3493 >>> mo1 = numpy.random.random((mol.nao_nr(), 10)) 3494 >>> mo2 = numpy.random.random((mol.nao_nr(), 8)) 3495 3496 >>> eri1 = mol.ao2mo(mo1) 3497 >>> print(eri1.shape) 3498 (55, 55) 3499 3500 >>> eri1 = mol.ao2mo(mo1, compact=False) 3501 >>> print(eri1.shape) 3502 (100, 100) 3503 3504 >>> eri1 = mol.ao2mo(eri, (mo1,mo2,mo2,mo2)) 3505 >>> print(eri1.shape) 3506 (80, 36) 3507 3508 >>> eri1 = mol.ao2mo(eri, (mo1,mo2,mo2,mo2), erifile='water.h5') 3509 ''' 3510 from pyscf import ao2mo 3511 return ao2mo.kernel(self, mo_coeffs, erifile, dataname, intor, **kwargs) 3512 3513 @contextlib.contextmanager 3514 def _TemporaryMoleContext(self, method, args, args_bak): 3515 '''Almost every method depends on the Mole environment. Ensure the 3516 modification in temporary environment being thread safe 3517 ''' 3518 haslock = hasattr(self, '_lock') 3519 if not haslock: 3520 self._lock = threading.RLock() 3521 3522 with self._lock: 3523 method(*args) 3524 try: 3525 yield 3526 finally: 3527 method(*args_bak) 3528 if not haslock: 3529 del self._lock 3530 3531def _parse_nuc_mod(str_or_int_or_fn): 3532 nucmod = NUC_POINT 3533 if callable(str_or_int_or_fn): 3534 nucmod = str_or_int_or_fn 3535 elif (isinstance(str_or_int_or_fn, (str, unicode)) and 3536 str_or_int_or_fn[0].upper() == 'G'): # 'gauss_nuc' 3537 nucmod = NUC_GAUSS 3538 elif str_or_int_or_fn != 0: 3539 nucmod = NUC_GAUSS 3540 return nucmod 3541 3542def _update_from_cmdargs_(mol): 3543 try: 3544 # Detect whether in Ipython shell 3545 __IPYTHON__ # noqa: 3546 return 3547 except Exception: 3548 pass 3549 3550 if not mol._built: # parse cmdline args only once 3551 opts = cmd_args.cmd_args() 3552 3553 if opts.verbose: 3554 mol.verbose = opts.verbose 3555 if opts.max_memory: 3556 mol.max_memory = opts.max_memory 3557 3558 if opts.output: 3559 mol.output = opts.output 3560 3561 3562def from_zmatrix(atomstr): 3563 '''>>> a = """H 3564 H 1 2.67247631453057 3565 H 1 4.22555607338457 2 50.7684795164077 3566 H 1 2.90305235726773 2 79.3904651036893 3 6.20854462618583""" 3567 >>> for x in zmat2cart(a): print(x) 3568 ['H', array([ 0., 0., 0.])] 3569 ['H', array([ 2.67247631, 0. , 0. ])] 3570 ['H', array([ 2.67247631, 0. , 3.27310166])] 3571 ['H', array([ 0.53449526, 0.30859098, 2.83668811])] 3572 ''' 3573 from pyscf.symm import rotation_mat 3574 atomstr = atomstr.replace(';','\n').replace(',',' ') 3575 symbols = [] 3576 coord = [] 3577 min_items_per_line = 1 3578 for line_id, line in enumerate(atomstr.splitlines()): 3579 line = line.strip() 3580 if line and line[0] != '#': 3581 rawd = line.split() 3582 if len(rawd) < min_items_per_line: 3583 raise ValueError('Zmatrix format error at L%d %s' % (line_id, line)) 3584 3585 symbols.append(rawd[0]) 3586 if len(rawd) < 3: 3587 coord.append(numpy.zeros(3)) 3588 min_items_per_line = 3 3589 elif len(rawd) == 3: 3590 if DISABLE_EVAL: 3591 coord.append(numpy.array((float(rawd[2]), 0, 0))) 3592 else: 3593 coord.append(numpy.array((eval(rawd[2]), 0, 0))) 3594 min_items_per_line = 5 3595 elif len(rawd) == 5: 3596 if DISABLE_EVAL: 3597 vals = rawd[1:] 3598 else: 3599 vals = eval(','.join(rawd[1:])) 3600 bonda = int(vals[0]) - 1 3601 bond = float(vals[1]) 3602 anga = int(vals[2]) - 1 3603 ang = float(vals[3])/180*numpy.pi 3604 assert(ang >= 0) 3605 v1 = coord[anga] - coord[bonda] 3606 if not numpy.allclose(v1[:2], 0): 3607 vecn = numpy.cross(v1, numpy.array((0.,0.,1.))) 3608 else: # on z 3609 vecn = numpy.array((0.,0.,1.)) 3610 rmat = rotation_mat(vecn, ang) 3611 c = numpy.dot(rmat, v1) * (bond/numpy.linalg.norm(v1)) 3612 coord.append(coord[bonda]+c) 3613 min_items_per_line = 7 3614 else: 3615 if DISABLE_EVAL: 3616 vals = rawd[1:] 3617 else: 3618 vals = eval(','.join(rawd[1:])) 3619 bonda = int(vals[0]) - 1 3620 bond = float(vals[1]) 3621 anga = int(vals[2]) - 1 3622 ang = float(vals[3])/180*numpy.pi 3623 assert(ang >= 0 and ang <= numpy.pi) 3624 v1 = coord[anga] - coord[bonda] 3625 v1 /= numpy.linalg.norm(v1) 3626 if ang < 1e-7: 3627 c = v1 * bond 3628 elif numpy.pi-ang < 1e-7: 3629 c = -v1 * bond 3630 else: 3631 diha = int(vals[4]) - 1 3632 dih = float(vals[5])/180*numpy.pi 3633 v2 = coord[diha] - coord[anga] 3634 vecn = numpy.cross(v2, -v1) 3635 vecn_norm = numpy.linalg.norm(vecn) 3636 if vecn_norm < 1e-7: 3637 if not numpy.allclose(v1[:2], 0): 3638 vecn = numpy.cross(v1, numpy.array((0.,0.,1.))) 3639 else: # on z 3640 vecn = numpy.array((0.,0.,1.)) 3641 rmat = rotation_mat(vecn, ang) 3642 c = numpy.dot(rmat, v1) * bond 3643 else: 3644 rmat = rotation_mat(v1, -dih) 3645 vecn = numpy.dot(rmat, vecn) / vecn_norm 3646 rmat = rotation_mat(vecn, ang) 3647 c = numpy.dot(rmat, v1) * bond 3648 coord.append(coord[bonda]+c) 3649 atoms = list(zip([_atom_symbol(x) for x in symbols], coord)) 3650 return atoms 3651zmat2cart = zmat = from_zmatrix 3652 3653def cart2zmat(coord): 3654 '''>>> c = numpy.array(( 3655 (0.000000000000, 1.889726124565, 0.000000000000), 3656 (0.000000000000, 0.000000000000, -1.889726124565), 3657 (1.889726124565, -1.889726124565, 0.000000000000), 3658 (1.889726124565, 0.000000000000, 1.133835674739))) 3659 >>> print(cart2zmat(c)) 3660 1 3661 1 2.67247631453057 3662 1 4.22555607338457 2 50.7684795164077 3663 1 2.90305235726773 2 79.3904651036893 3 6.20854462618583 3664 ''' 3665 zstr = [] 3666 zstr.append('1') 3667 if len(coord) > 1: 3668 r1 = coord[1] - coord[0] 3669 nr1 = numpy.linalg.norm(r1) 3670 zstr.append('1 %.15g' % nr1) 3671 if len(coord) > 2: 3672 r2 = coord[2] - coord[0] 3673 nr2 = numpy.linalg.norm(r2) 3674 a = numpy.arccos(numpy.dot(r1,r2)/(nr1*nr2)) 3675 zstr.append('1 %.15g 2 %.15g' % (nr2, a*180/numpy.pi)) 3676 if len(coord) > 3: 3677 o0, o1, o2 = coord[:3] 3678 p0, p1, p2 = 1, 2, 3 3679 for k, c in enumerate(coord[3:]): 3680 r0 = c - o0 3681 nr0 = numpy.linalg.norm(r0) 3682 r1 = o1 - o0 3683 nr1 = numpy.linalg.norm(r1) 3684 a1 = numpy.arccos(numpy.dot(r0,r1)/(nr0*nr1)) 3685 b0 = numpy.cross(r0, r1) 3686 nb0 = numpy.linalg.norm(b0) 3687 3688 if abs(nb0) < 1e-7: # o0, o1, c in line 3689 a2 = 0 3690 zstr.append('%d %.15g %d %.15g %d %.15g' % 3691 (p0, nr0, p1, a1*180/numpy.pi, p2, a2)) 3692 else: 3693 b1 = numpy.cross(o2-o0, r1) 3694 nb1 = numpy.linalg.norm(b1) 3695 3696 if abs(nb1) < 1e-7: # o0 o1 o2 in line 3697 a2 = 0 3698 zstr.append('%d %.15g %d %.15g %d %.15g' % 3699 (p0, nr0, p1, a1*180/numpy.pi, p2, a2)) 3700 o2 = c 3701 p2 = 4 + k 3702 else: 3703 if numpy.dot(numpy.cross(b1, b0), r1) < 0: 3704 a2 = numpy.arccos(numpy.dot(b1, b0) / (nb0*nb1)) 3705 else: 3706 a2 =-numpy.arccos(numpy.dot(b1, b0) / (nb0*nb1)) 3707 zstr.append('%d %.15g %d %.15g %d %.15g' % 3708 (p0, nr0, p1, a1*180/numpy.pi, p2, a2*180/numpy.pi)) 3709 3710 return '\n'.join(zstr) 3711 3712def dyall_nuc_mod(nuc_charge, nucprop={}): 3713 ''' Generate the nuclear charge distribution parameter zeta 3714 rho(r) = nuc_charge * Norm * exp(-zeta * r^2) 3715 3716 Ref. L. Visscher and K. Dyall, At. Data Nucl. Data Tables, 67, 207 (1997) 3717 ''' 3718 mass = nucprop.get('mass', elements.ISOTOPE_MAIN[nuc_charge]) 3719 r = (0.836 * mass**(1./3) + 0.570) / 52917.7249 3720 zeta = 1.5 / (r**2) 3721 return zeta 3722 3723def filatov_nuc_mod(nuc_charge, nucprop={}): 3724 ''' Generate the nuclear charge distribution parameter zeta 3725 rho(r) = nuc_charge * Norm * exp(-zeta * r^2) 3726 3727 Ref. M. Filatov and D. Cremer, Theor. Chem. Acc. 108, 168 (2002) 3728 M. Filatov and D. Cremer, Chem. Phys. Lett. 351, 259 (2002) 3729 ''' 3730 c = param.LIGHT_SPEED 3731 nuc_charge = charge(nuc_charge) 3732 r = (-0.263188*nuc_charge + 106.016974 + 138.985999/nuc_charge) / c**2 3733 zeta = 1 / (r**2) 3734 return zeta 3735 3736def fakemol_for_charges(coords, expnt=1e16): 3737 '''Construct a fake Mole object that holds the charges on the given 3738 coordinates (coords). The shape of the charge can be a normal 3739 distribution with the Gaussian exponent (expnt). 3740 ''' 3741 nbas = coords.shape[0] 3742 fakeatm = numpy.zeros((nbas,ATM_SLOTS), dtype=numpy.int32) 3743 fakebas = numpy.zeros((nbas,BAS_SLOTS), dtype=numpy.int32) 3744 fakeenv = [0] * PTR_ENV_START 3745 ptr = PTR_ENV_START 3746 fakeatm[:,PTR_COORD] = numpy.arange(ptr, ptr+nbas*3, 3) 3747 fakeenv.append(coords.ravel()) 3748 ptr += nbas*3 3749 fakebas[:,ATOM_OF] = numpy.arange(nbas) 3750 fakebas[:,NPRIM_OF] = 1 3751 fakebas[:,NCTR_OF] = 1 3752# approximate point charge with gaussian distribution exp(-1e16*r^2) 3753 fakebas[:,PTR_EXP] = ptr 3754 fakebas[:,PTR_COEFF] = ptr+1 3755 fakeenv.append([expnt, 1/(2*numpy.sqrt(numpy.pi)*gaussian_int(2,expnt))]) 3756 ptr += 2 3757 fakemol = Mole() 3758 fakemol._atm = fakeatm 3759 fakemol._bas = fakebas 3760 fakemol._env = numpy.hstack(fakeenv) 3761 fakemol._built = True 3762 return fakemol 3763 3764del(BASE) 3765