1################################################################## 2## (c) Copyright 2015- by Jaron T. Krogel ## 3################################################################## 4 5 6import os 7from numpy import array,sqrt,arange,linspace,zeros,exp 8from generic import obj 9from periodic_table import is_element 10from developer import DevBase,error,to_str 11from fileio import TextFile 12from plotting import * 13from debug import * 14 15 16 17def show_plots(): 18 show() 19#end def show_plots 20 21# container class for available basis set files 22class BasisSets(DevBase): 23 def __init__(self,*basissets): 24 if len(basissets)==1 and isinstance(basissets[0],list): 25 basissets = basissets[0] 26 #end if 27 bsfiles = [] 28 bss = [] 29 errors = False 30 for bs in basissets: 31 if isinstance(bs,BasisFile): 32 bss.append(bs) 33 elif isinstance(bs,str): 34 bsfiles.append(bs) 35 else: 36 self.error('expected BasisFile type or filepath, got '+str(type(bs)),exit=False) 37 errors = True 38 #end if 39 #end for 40 if errors: 41 self.error('cannot create Basissets object') 42 #end if 43 44 if len(bss)>0: 45 self.addbs(bss) 46 #end if 47 if len(bsfiles)>0: 48 self.readbs(bsfiles) 49 #end if 50 #end def __init__ 51 52 53 def addbs(self,*basissets): 54 if len(basissets)==1 and isinstance(basissets[0],list): 55 basissets = basissets[0] 56 #end if 57 for bs in basissets: 58 self[bs.filename] = bs 59 #end for 60 #end def addbs 61 62 63 def readbs(self,*bsfiles): 64 if len(bsfiles)==1 and isinstance(bsfiles[0],list): 65 bsfiles = bsfiles[0] 66 #end if 67 bss = [] 68 self.log('') 69 self.log(' Basissets') 70 for filepath in bsfiles: 71 self.log(' reading basis: '+filepath) 72 ext = filepath.split('.')[-1].lower() 73 if ext=='gms_bas' or ext=='bas': 74 bs = gamessBasisFile(filepath) 75 else: 76 bs = BasisFile(filepath) 77 #end if 78 bss.append(bs) 79 #end for 80 self.log('') 81 self.addbs(bss) 82 #end def readbs 83 84 85 # test needed 86 def bases_by_atom(self,*bsfiles): 87 bss = obj() 88 for bsfile in bsfiles: 89 if bsfile in self: 90 bs = self[bsfile] 91 bss[bs.element_label] = bs 92 else: 93 self.error('basis file not found\nmissing file: {0}'.format(bsfile)) 94 #end if 95 #end for 96 return bss 97 #end def bases_by_atom 98#end class BasisSets 99 100 101class BasisFile(DevBase): 102 def __init__(self,filepath=None): 103 self.element = None 104 self.element_label = None 105 self.filename = None 106 self.location = None 107 if filepath!=None: 108 self.filename = os.path.basename(filepath) 109 self.location = os.path.abspath(filepath) 110 elem_label = self.filename.split('.')[0] 111 is_elem,symbol = is_element(elem_label,symbol=True) 112 if not is_elem: 113 self.error('cannot determine element for basis file: {0}\nbasis file names must be prefixed by an atomic symbol or label\n(e.g. Si, Si1, etc)'.format(filepath)) 114 #end if 115 self.element = symbol 116 self.element_label = elem_label 117 #end if 118 #end def __init__ 119 120 def cleaned_text(self): 121 self.not_implemented() 122 #end def cleaned_text 123#end class BasisFile 124 125 126 127 128class gaussBasisFile(BasisFile): 129 angular_terms = 'spdfghiklmn' 130 131 def __init__(self,filepath=None): 132 BasisFile.__init__(self,filepath) 133 self.text = None 134 if filepath!=None: 135 self.read(filepath) 136 #end if 137 #end def __init__ 138 139 def cleaned_text(self): 140 if self.text is None: 141 self.error('text requested prior to read\nfile: {0}'.format(self.location)) 142 #end if 143 return self.text 144 #end def cleaned_text 145 146 def read(self,filepath=None): 147 if filepath is None: 148 filepath = self.location 149 #end if 150 if not os.path.exists(filepath): 151 self.error('file does not exist: {0}'.format(filepath)) 152 #end if 153 file = TextFile(filepath) 154 self.read_file(file) 155 #end def read 156 157 def read_file(self,file): 158 self.not_implemented() 159 #end def read_file 160#end class gaussBasisFile 161 162 163 164class gamessBasisFile(gaussBasisFile): 165 def read_file(self,file): 166 dstart = file.find('$DATA') 167 estart = file.find('$END') 168 if dstart!=-1: 169 file.seek(dstart) 170 file.readline() 171 cur = file.tell() 172 tokens = file.readtokens() 173 if len(tokens)==2 and tokens[0].lower() in self.angular_terms and tokens[1].isdigit(): 174 file.seek(cur) # skip atom header line 175 #end if 176 if estart==-1: 177 estart = file.size() 178 #end if 179 self.text = to_str(file[file.tell():estart].strip()) 180 else: 181 lines = to_str(file[:]).splitlines() 182 blines = [] 183 for line in lines: 184 ls = line.strip() 185 if len(ls)>0 and ls[0]!='!' and ls[0]!='#': 186 blines.append(line) 187 #end if 188 #end for 189 if len(blines)>0: 190 tokens = blines[0].split() 191 if not (len(tokens)==2 and tokens[0].lower() in self.angular_terms and tokens[1].isdigit()): 192 blines = blines[1:] # skip atom header line 193 #end if 194 #end if 195 text = '' 196 for line in blines: 197 text += line + '\n' 198 #end for 199 self.text = text.strip() 200 #end if 201 #end def read_file 202#end class gamessBasisFile 203 204 205 206 207 208 209 210def process_gaussian_text(text,format,pp=True,basis=True,preserve_spacing=False): 211 if format=='gamess' or format=='gaussian' or format=='atomscf': 212 rawlines = text.splitlines() 213 sections = [] 214 last_empty = True 215 for rline in rawlines: 216 line = rline.strip() 217 if (not line.startswith('!')) and (not line.startswith('#')) and (not line.startswith('*')) and len(line)>0 and line!='$DATA' and line!='$END': 218 if last_empty: 219 lines = [] 220 sections.append(lines) 221 #end if 222 if preserve_spacing: 223 lines.append(rline) 224 else: 225 lines.append(line) 226 #end if 227 last_empty = False 228 else: 229 last_empty = True 230 #end if 231 #end for 232 del lines 233 if len(sections)==2: 234 basis_lines = sections[0] 235 pp_lines = sections[1] 236 elif pp: 237 basis_lines = None 238 pp_lines = sections[0] 239 elif basis: 240 basis_lines = sections[0] 241 pp_lines = None 242 #end if 243 elif format=='crystal': 244 rawlines = text.splitlines() 245 pp_lines = [] 246 basis_lines = [] 247 foundpp = False 248 for line in rawlines: 249 if not foundpp: 250 foundpp = len(line.split())==5 251 #end if 252 if not foundpp: 253 pp_lines.append(line) 254 else: 255 basis_lines.append(line) 256 #end if 257 #end for 258 if len(pp_lines)==0: 259 pp_lines = None 260 #end if 261 if len(basis_lines)==0: 262 basis_lines = None 263 #end if 264 else: 265 error('{0} format is unknown'.format(format),'process_gaussian_text') 266 #end if 267 if pp and basis: 268 return pp_lines,basis_lines 269 elif pp: 270 return pp_lines 271 elif basis: 272 return basis_lines 273 else: 274 error('must request pp or basis') 275 #end if 276#end def process_gaussian_text 277 278 279 280 281class GaussianBasisSet(DevBase): 282 lset_full = tuple('spdfghijk') 283 lstyles = obj(s='g-',p='r-',d='b-',f='m-',g='c-',h='k-',i='g-.',j='r-.',k='b-.') 284 formats = 'gaussian gamess'.split() 285 286 crystal_lmap = {0:'s',1:'sp',2:'p',3:'d',4:'f'} 287 crystal_lmap_reverse = dict(s=0,sp=1,p=2,d=3,f=4) 288 289 @staticmethod 290 def process_float(s): 291 return float(s.replace('D','e').replace('d','e')) 292 #end def process_float 293 294 295 def __init__(self,filepath=None,format=None): 296 self.name = None 297 self.basis = obj() 298 if filepath!=None: 299 self.read(filepath,format) 300 #end if 301 #end def __init__ 302 303 304 def read(self,filepath,format=None): 305 if format is None: 306 self.error('format keyword must be specified to read file {0}\nvalid options are: {1}'.format(filepath,self.formats)) 307 elif not format in self.formats: 308 self.error('incorrect format requested: {0}\nvalid options are: {1}'.format(format,self.formats)) 309 #end if 310 if not os.path.exists(filepath): 311 self.error('cannot read {0}, file does not exist'.format(filepath)) 312 #end if 313 #self.name = split_delims(os.path.split(filepath)[1])[0] 314 self.name = os.path.split(filepath)[1].split('.')[0] 315 text = open(filepath,'r').read() 316 self.read_text(text,format) 317 #end def read 318 319 320 def write(self,filepath=None,format=None): 321 if format is None: 322 self.error('format keyword must be specified to write file {0}\nvalid options are: {1}'.format(filepath,self.formats)) 323 elif not format in self.formats: 324 self.error('incorrect format requested: {0}\nvalid options are: {1}'.format(format,self.formats)) 325 #end if 326 text = self.write_text(format) 327 if filepath!=None: 328 open(filepath,'w').write(text) 329 #end if 330 return text 331 #end def write 332 333 334 def read_text(self,text,format=None): 335 basis_lines = process_gaussian_text(text,format,pp=False) 336 self.read_lines(basis_lines,format) 337 #end def read_text 338 339 340 def read_lines(self,basis_lines,format=None): 341 basis = self.basis 342 basis.clear() 343 if format=='gamess': 344 i=1 345 while i<len(basis_lines): 346 tokens = basis_lines[i].split(); i+=1 347 ltext = tokens[0].lower() 348 ngauss = int(tokens[1]) 349 scale = array(tokens[2:],dtype=float) 350 bterms = obj() 351 for j in range(ngauss): 352 index,expon,coeff = basis_lines[i].split(); i+=1 353 expon = GaussianBasisSet.process_float(expon) 354 coeff = GaussianBasisSet.process_float(coeff) 355 bterms.append(obj(expon=expon,coeff=coeff)) 356 #end for 357 basis.append(obj(l=ltext,scale=scale,terms=bterms)) 358 #end while 359 #end if 360 elif format=='gaussian': 361 i=1 362 while i<len(basis_lines): 363 tokens = basis_lines[i].split(); i+=1 364 ltext = tokens[0].lower() 365 ngauss = int(tokens[1]) 366 scale = array(tokens[2:],dtype=float) 367 bterms = obj() 368 for j in range(ngauss): 369 expon,coeff = basis_lines[i].split(); i+=1 370 expon = GaussianBasisSet.process_float(expon) 371 coeff = GaussianBasisSet.process_float(coeff) 372 bterms.append(obj(expon=expon,coeff=coeff)) 373 #end for 374 basis.append(obj(l=ltext,scale=scale,terms=bterms)) 375 #end while 376 elif format=='crystal': 377 i=0 378 while i<len(basis_lines): 379 tokens = basis_lines[i].split(); i+=1 380 if len(tokens)!=5: 381 self.error('could not parse crystal basisset, input may be misformatted') 382 #end if 383 basis_type = int(tokens[0]) 384 l_type = int(tokens[1]) 385 ngauss = int(tokens[2]) 386 formal_charge = float(tokens[3]) 387 scale = array([tokens[4]],dtype=float) 388 ltext = GaussianBasisSet.crystal_lmap[l_type] 389 if ltext!='sp': 390 bterms = obj() 391 for j in range(ngauss): 392 expon,coeff = basis_lines[i].split(); i+=1 393 expon = GaussianBasisSet.process_float(expon) 394 coeff = GaussianBasisSet.process_float(coeff) 395 bterms.append(obj(expon=expon,coeff=coeff)) 396 #end for 397 basis.append(obj(l=ltext,scale=scale,terms=bterms)) 398 else: # sp has shared exponent for s and p, split them now 399 sterms = obj() 400 pterms = obj() 401 for j in range(ngauss): 402 expon,scoeff,pcoeff = basis_lines[i].split(); i+=1 403 expon = GaussianBasisSet.process_float(expon) 404 scoeff = GaussianBasisSet.process_float(scoeff) 405 pcoeff = GaussianBasisSet.process_float(pcoeff) 406 sterms.append(obj(expon=expon,coeff=scoeff)) 407 pterms.append(obj(expon=expon,coeff=pcoeff)) 408 #end for 409 basis.append(obj(l='s',scale=scale,terms=sterms)) 410 basis.append(obj(l='p',scale=scale,terms=pterms)) 411 #end if 412 #end while 413 else: 414 self.error('ability to read file format {0} has not been implemented'.format(format)) 415 #end if 416 # sort the basis in s,p,d,f,... order 417 self.lsort() 418 #end def read_lines 419 420 421 422 def write_text(self,format=None,occ=None): 423 text = '' 424 format = format.lower() 425 if format=='gamess': 426 #text += '{0} {1} 0. 0. 0.\n'.format(self.element,self.Zcore+self.Zval) 427 for ib in range(len(self.basis)): 428 b = self.basis[ib] 429 line = '{0} {1}'.format(b.l,len(b.terms)) 430 for s in b.scale: 431 line += ' {0}'.format(s) 432 #end for 433 text += line + '\n' 434 for it in range(len(b.terms)): 435 t = b.terms[it] 436 text += '{0:<4} {1:12.8E} {2: 12.8E}\n'.format(it+1,t.expon,t.coeff) 437 #end for 438 #end for 439 elif format=='gaussian': 440 #text += '{0} 0\n'.format(self.element) 441 for ib in range(len(self.basis)): 442 b = self.basis[ib] 443 line = '{0} {1}'.format(b.l,len(b.terms)) 444 for s in b.scale: 445 line += ' {0}'.format(s) 446 #end for 447 text += line + '\n' 448 for it in range(len(b.terms)): 449 t = b.terms[it] 450 text += '{0:12.8E} {1: 12.8E}\n'.format(t.expon,t.coeff) 451 #end for 452 #end for 453 elif format=='crystal': 454 if occ is not None: 455 lcounts = dict(s=0,p=0,d=0,f=0) 456 #end if 457 for ib in range(len(self.basis)): 458 b = self.basis[ib] 459 if b.l not in self.crystal_lmap_reverse: 460 self.error('{0} channels cannot be handled by crystal'.format(b.l)) 461 #end if 462 Zf = 0 463 if occ is not None and b.l in occ and lcounts[b.l]<len(occ[b.l]): 464 Zf = occ[b.l][lcounts[b.l]] 465 lcounts[b.l]+=1 466 #end if 467 lnum = self.crystal_lmap_reverse[b.l] 468 line = '0 {0} {1} {2} {3}'.format(lnum,len(b.terms),Zf,b.scale[0]) 469 text += line + '\n' 470 for it in range(len(b.terms)): 471 t = b.terms[it] 472 text += '{0:12.8E} {1: 12.8E}\n'.format(t.expon,t.coeff) 473 #end for 474 #end for 475 else: 476 self.error('ability to write file format {0} has not been implemented'.format(format)) 477 #end if 478 return text 479 #end def write_text 480 481 482 # test needed 483 def size(self): 484 return len(self.basis) 485 #end def size 486 487 488 # test needed 489 def lset(self): 490 lset = set() 491 for bf in self.basis: 492 lset.add(bf.l) 493 #end for 494 return lset 495 #end def lset 496 497 498 # test needed 499 def lcount(self): 500 return len(self.lset()) 501 #end def lcount 502 503 504 # test needed 505 def lbasis(self): 506 lbasis = obj() 507 for n in range(len(self.basis)): 508 bf = self.basis[n] 509 l = bf.l 510 if l not in lbasis: 511 lbasis[l] = obj() 512 #end if 513 lbasis[l].append(bf) 514 #end for 515 return lbasis 516 #end def lbasis 517 518 519 # test needed 520 def lsort(self): 521 lbasis = self.lbasis() 522 self.basis.clear() 523 for l in self.lset_full: 524 if l in lbasis: 525 lbas = lbasis[l] 526 for n in range(len(lbas)): 527 bf = lbas[n] 528 self.basis.append(bf) 529 #end for 530 #end if 531 #end for 532 #end def lsort 533 534 535 # test needed 536 def uncontracted(self): 537 all_uncon = True 538 for bf in self.basis: 539 all_uncon &= len(bf.terms)==1 540 #end for 541 return all_uncon 542 #end def uncontracted 543 544 545 # test needed 546 def contracted(self): 547 return not self.uncontracted() 548 #end def contracted 549 550 551 # test needed 552 def uncontract(self,tol=1e-3): 553 if self.uncontracted(): 554 return 555 #end if 556 lbasis = self.lbasis() 557 self.basis.clear() 558 for l in self.lset_full: 559 if l in lbasis: 560 exponents = [] 561 lbas = lbasis[l] 562 for n in range(len(lbas)): 563 uterms = lbas[n].terms 564 for i in range(len(uterms)): 565 expon = uterms[i].expon 566 if len(exponents)==0: 567 exponents = array([expon],dtype=float) 568 elif abs(exponents-expon).min()>tol: 569 exponents = array(list(exponents)+[expon],dtype=float) 570 #end if 571 #end for 572 #end for 573 for expon in exponents: 574 cterms = obj() 575 cterms.append(obj(expon=expon,coeff=1.0)) 576 bf = obj(l=l,scale=array([1.0]),terms=cterms) 577 self.basis.append(bf) 578 #end for 579 #end if 580 #end for 581 #end def uncontract 582 583 584 # test needed 585 def contracted_basis_size(self): 586 bcount = obj() 587 for bf in self.basis: 588 l = bf.l 589 if l not in bcount: 590 bcount[l]=0 591 #end if 592 bcount[l] += 1 593 #end for 594 bs = '' 595 for l in self.lset_full: 596 if l in bcount: 597 bs += str(bcount[l])+l 598 #end if 599 #end for 600 return bs 601 #end def contracted_basis_size 602 603 604 # test needed 605 def uncontracted_basis_size(self): 606 if self.uncontracted(): 607 return self.contracted_basis_size() 608 #end if 609 uc = self.copy() 610 uc.uncontract() 611 return uc.contracted_basis_size() 612 #end def uncontracted_basis_size 613 614 615 # test needed 616 def basis_size(self): 617 us = self.uncontracted_basis_size() 618 cs = self.contracted_basis_size() 619 return '({0})/[{1}]'.format(us,cs) 620 #end def basis_size 621 622 623 # test needed 624 def prim_expons(self): 625 if self.contracted(): 626 self.error('cannot find primitive gaussian expons because basis is contracted') 627 #end if 628 lbasis = self.lbasis() 629 gexpon = obj() 630 for l,lbas in lbasis.items(): 631 e = [] 632 for n in range(len(lbas)): 633 e.append(lbas[n].terms[0].expon) 634 #end for 635 gexpon[l] = array(e,dtype=float) 636 #end for 637 return gexpon 638 #end def prim_expons 639 640 641 # test needed 642 def prim_widths(self): 643 if self.contracted(): 644 self.error('cannot find primitive gaussian widths because basis is contracted') 645 #end if 646 lbasis = self.lbasis() 647 gwidth = obj() 648 for l,lbas in lbasis.items(): 649 w = [] 650 for n in range(len(lbas)): 651 w.append(1./sqrt(2.*lbas[n].terms[0].expon)) 652 #end for 653 gwidth[l] = array(w,dtype=float) 654 #end for 655 return gwidth 656 #end def prim_widths 657 658 659 # test needed 660 def remove_prims(self,comp=None,keep=None,**lselectors): 661 lbasis = self.lbasis() 662 if comp!=None: 663 gwidths = self.prim_widths() 664 #end if 665 for l,lsel in lselectors.items(): 666 if l not in lbasis: 667 self.error('cannot remove basis functions from channel {0}, channel not present'.format(l)) 668 #end if 669 lbas = lbasis[l] 670 if isinstance(lsel,float): 671 rcut = lsel 672 less = False 673 if comp=='<': 674 less = True 675 elif comp=='>': 676 less = False 677 elif comp is None: 678 self.error('comp argument must be provided (< or >)') 679 else: 680 self.error('comp must be < or >, you provided: {0}'.format(comp)) 681 #end if 682 gw = gwidths[l] 683 iw = arange(len(gw)) 684 nkeep = 0 685 if keep!=None and l in keep: 686 nkeep = keep[l] 687 #end if 688 if less: 689 rem = iw[gw<rcut] 690 for i in range(len(rem)-nkeep): 691 del lbas[rem[i]] 692 #end for 693 else: 694 rem = iw[gw>rcut] 695 for i in range(nkeep,len(rem)): 696 del lbas[rem[i]] 697 #end for 698 #end if 699 elif isinstance(lsel,int): 700 if comp=='<': 701 if lsel>len(lbas): 702 self.error('cannot remove {0} basis functions from channel {1} as it only has {2}'.format(lsel,l,len(lbas))) 703 #end if 704 for i in range(lsel): 705 del lbas[i] 706 #end for 707 elif comp=='>': 708 if lsel>len(lbas): 709 self.error('cannot remove {0} basis functions from channel {1} as it only has {2}'.format(lsel,l,len(lbas))) 710 #end if 711 for i in range(len(lbas)-lsel,len(lbas)): 712 del lbas[i] 713 #end for 714 else: 715 if lsel>=len(lbas): 716 self.error('cannot remove basis function {0} from channel {1} as it only has {2}'.format(lsel,l,len(lbas))) 717 #end if 718 del lbas[lsel] 719 #end if 720 else: 721 for ind in lsel: 722 del lbas[ind] 723 #end for 724 #end if 725 #end for 726 self.basis.clear() 727 for l in self.lset_full: 728 if l in lbasis: 729 lbas = lbasis[l] 730 for k in sorted(lbas.keys()): 731 self.basis.append(lbas[k]) 732 #end for 733 #end if 734 #end for 735 #end def remove_prims 736 737 738 # test needed 739 def remove_small_prims(self,**keep): 740 lsel = obj() 741 for l,lbas in self.lbasis().items(): 742 if l in keep: 743 lsel[l] = len(lbas)-keep[l] 744 #end if 745 #end for 746 self.remove_prims(comp='<',**lsel) 747 #end def remove_small_prims 748 749 750 # test needed 751 def remove_large_prims(self,**keep): 752 lsel = obj() 753 for l,lbas in self.lbasis().items(): 754 if l in keep: 755 lsel[l] = len(lbas)-keep[l] 756 #end if 757 #end for 758 self.remove_prims(comp='>',**lsel) 759 #end def remove_large_prims 760 761 762 # test needed 763 def remove_small_prims_rel(self,other,**keep): 764 gwidths = other.prim_widths() 765 lsel = obj() 766 for l,gw in gwidths.items(): 767 lsel[l] = gw.min() 768 #end for 769 self.remove_prims(comp='<',keep=keep,**lsel) 770 #end def remove_small_prims_rel 771 772 773 # test needed 774 def remove_large_prims_rel(self,other,**keep): 775 gwidths = other.prim_widths() 776 lsel = obj() 777 for l,gw in gwidths.items(): 778 lsel[l] = gw.max() 779 #end for 780 self.remove_prims(comp='>',keep=keep,**lsel) 781 #end def remove_large_prims_rel 782 783 784 # test needed 785 def remove_channels(self,llist): 786 lbasis = self.lbasis() 787 for l in llist: 788 if l in lbasis: 789 del lbasis[l] 790 #end if 791 #end for 792 self.basis.clear() 793 for l in self.lset_full: 794 if l in lbasis: 795 lbas = lbasis[l] 796 for n in range(len(lbas)): 797 bf = lbas[n] 798 self.basis.append(bf) 799 #end for 800 #end if 801 #end for 802 #end def remove_channels 803 804 805 # test needed 806 def incorporate(self,other,tol=1e-3,unique=False): 807 uncontracted = self.uncontracted() and other.uncontracted() 808 lbasis = self.lbasis() 809 lbasis_other = other.lbasis() 810 if uncontracted: 811 gwidths = self.prim_widths() 812 gwidths_other = other.prim_widths() 813 #end if 814 self.basis.clear() 815 if not unique: # simple, direct merge of basis sets 816 for l in self.lset_full: 817 if l in lbasis: 818 lbas = lbasis[l] 819 for n in range(len(lbas)): 820 bf = lbas[n] 821 self.basis.append(bf) 822 #end for 823 #end if 824 if l in lbasis_other: 825 lbas = lbasis_other[l] 826 for n in range(len(lbas)): 827 bf = lbas[n] 828 self.basis.append(bf) 829 #end for 830 #end if 831 #end for 832 else: # merge uncontracted basis sets preserving order 833 for l in self.lset_full: 834 primitives = [] 835 widths = [] 836 orig_widths = array([]) 837 if l in lbasis: 838 primitives.extend(lbasis[l].list()) 839 widths.extend(gwidths[l]) 840 orig_widths = gwidths[l] 841 #end if 842 if l in lbasis_other: 843 prims = lbasis_other[l].list() 844 owidths = gwidths_other[l] 845 for n in range(len(prims)): 846 w = owidths[n] 847 if len(orig_widths)==0 or abs(orig_widths-w).min()>tol: 848 primitives.append(prims[n]) 849 widths.append(w) 850 #end if 851 #end if 852 #end if 853 primitives = array(primitives,dtype=object)[array(widths).argsort()] 854 for bf in primitives: 855 self.basis.append(bf) 856 #end for 857 #end for 858 #end if 859 #end def incorporate 860 861 862 def plot(self,r=None,rmin=0.01,rmax=8.0,show=True,fig=True,sep=False,prim=False,style=None,fmt=None,nsub=None): 863 if r is None: 864 r = linspace(rmin,rmax,1000) 865 #end if 866 if not prim: 867 ptitle = '{0} {1} basis'.format(self.name,self.basis_size()) 868 else: 869 ptitle = '{0} {1} primitives'.format(self.name,self.basis_size()) 870 #end if 871 if fig: 872 figure() 873 #end if 874 r2 = r**2 875 lcount = self.lcount() 876 if nsub!=None: 877 lcount = max(lcount,nsub) 878 #end if 879 lbasis = self.lbasis() 880 lc = 0 881 for l in self.lset_full: 882 if l in lbasis: 883 lc+=1 884 if sep: 885 subplot(lcount,1,lc) 886 ylabel(l) 887 if lc==1: 888 title(ptitle) 889 #end if 890 #end if 891 lstyle=self.lstyles[l] 892 if style!=None: 893 lstyle = lstyle[0]+style 894 #end if 895 if fmt!=None: 896 lstyle=fmt 897 #end if 898 lbas = lbasis[l] 899 for n in range(len(lbas)): 900 bf = lbas[n] 901 br = zeros(r.shape) 902 s = bf.scale[0] 903 for pf in bf.terms: 904 c = pf.coeff 905 a = -pf.expon*s**2 906 pr = exp(a*r2) 907 if not prim: 908 br += c*pr 909 else: 910 plot(r,pr,lstyle,label=l) 911 #end if 912 #end for 913 if not prim: 914 plot(r,br,lstyle,label=l) 915 #end if 916 #end for 917 #end if 918 #end for 919 if fig: 920 if not sep: 921 if self.name!=None: 922 title(ptitle) 923 #end if 924 ylabel('basis functions') 925 legend() 926 #end if 927 xlabel('r') 928 if show: 929 show_plots() 930 #end if 931 #end if 932 #end def plot 933 934 935 def plot_primitives(self): 936 None 937 #end def plot_primitives 938 939 940 def plot_prim_widths(self,show=True,fig=True,sep=False,style='o',fmt=None,nsub=None,semilog=True,label=True): 941 if self.contracted(): 942 self.error('cannot plot primitive gaussian widths because basis is contracted') 943 #end if 944 ptitle = '{0} {1} primitive widths'.format(self.name,self.basis_size()) 945 if fig: 946 figure() 947 #end if 948 pwidths = self.prim_widths() 949 lcount = self.lcount() 950 if nsub!=None: 951 lcount = max(lcount,nsub) 952 #end if 953 lbasis = self.lbasis() 954 lc = 0 955 for l in self.lset_full: 956 if l in lbasis: 957 lwidths = pwidths[l] 958 lc+=1 959 if sep: 960 subplot(lcount,1,lc) 961 ylabel(l) 962 if lc==1: 963 title(ptitle) 964 #end if 965 #end if 966 lstyle=self.lstyles[l] 967 if style!=None: 968 lstyle = lstyle[0]+style 969 #end if 970 if fmt!=None: 971 lstyle=fmt 972 #end if 973 if semilog: 974 plotter = semilogy 975 else: 976 plotter = plot 977 #end if 978 plotter(list(range(len(lwidths))),lwidths,lstyle,label=l) 979 #end if 980 #end for 981 if label: 982 if not sep: 983 if self.name!=None: 984 title(ptitle) 985 #end if 986 ylabel('primitive widths') 987 legend() 988 #end if 989 xlabel('primitive index') 990 #end if 991 if show: 992 show_plots() 993 #end if 994 #end def plot_prim_widths 995#end class GaussianBasisSet 996