SciCodePile/SciCode-Domain-Code
DATA1: Domain-Specific Code Dataset Dataset Overview DATA1 is a large-scale domain-specific code dataset focusing on code samples from interdisciplinary fields such as biology, chemistry, materials science, and related areas. The dataset is collected and organized from GitHub repositories, covering 178 different domain topics with over 1.1 billion lines of code. Dataset Statistics Total Datasets: 178 CSV files Total Data Size: ~115 GB Total Lines… See the full description on the dataset page: https://huggingface.co/datasets/SciCodePile/SciCode-Domain-Code.
42.3k
1"keyword","repo_name","file_path","file_extension","file_size","line_count","content","language"
2"Homo-Lumo","setten/GW100","docs/index.md",".md","40","4","test3 4[name](../pages/compare_all.html)5","Markdown"
6"Homo-Lumo","setten/GW100","scripts/datasetclass.ipynb",".ipynb","3056","126","{7 ""cells"": [8 {9 ""cell_type"": ""code"",10 ""execution_count"": 1,11 ""metadata"": {12 ""collapsed"": true13 },14 ""outputs"": [],15 ""source"": [16 ""import json""17 ]18 },19 {20 ""cell_type"": ""code"",21 ""execution_count"": 2,22 ""metadata"": {23 ""collapsed"": true24 },25 ""outputs"": [],26 ""source"": [27 ""class GW100DataSet(object):\n"",28 "" \n"",29 "" def __init__(self):\n"",30 "" self.data = {'data': {}}\n"",31 "" \n"",32 "" @classmethod\n"",33 "" def from_txt(cls, file_name):\n"",34 "" set = GW100DataSet()\n"",35 "" set.read_txt(file_name)\n"",36 "" return set\n"",37 "" \n"",38 "" def name(self, post_fix=''):\n"",39 "" return '%s_%s_%sv%s_%s%s' % (self.data['calc_type'], self.data['orbital'], self.data['code'][0], \n"",40 "" self.data['code_version'], self.data['basis_name'], post_fix)\n"",41 "" \n"",42 "" def read_txt(self, file_name):\n"",43 "" with open(file_name, 'r') as f:\n"",44 "" lines = f.readlines()\n"",45 "" for line in lines:\n"",46 "" if '=' in line:\n"",47 "" d = line.strip().split('=')\n"",48 "" if d[1].startswith('{'):\n"",49 "" self.data[d[0]] = json.loads(d[1])\n"",50 "" else:\n"",51 "" self.data[d[0]] = d[1]\n"",52 "" else:\n"",53 "" d = line.strip().split( )\n"",54 "" try:\n"",55 "" self.data['data'][d[0]]=float(d[1])\n"",56 "" except (IndexError, ValueError):\n"",57 "" self.data['data'][d[0]]='null'\n"",58 "" \n"",59 "" def dump_json(self, post):\n"",60 "" with open(self.name(post)+'.json','w') as fp:\n"",61 "" json.dump(self.data, fp, indent=4)""62 ]63 },64 {65 ""cell_type"": ""code"",66 ""execution_count"": 8,67 ""metadata"": {},68 ""outputs"": [69 {70 ""data"": {71 ""text/plain"": [72 ""'G0W0@PBE_HOMO_Tv7.0_def2-TQZVP_cbas'""73 ]74 },75 ""execution_count"": 8,76 ""metadata"": {},77 ""output_type"": ""execute_result""78 }79 ],80 ""source"": [81 ""dataset = GW100DataSet.from_txt('test_data/tmcbastqex.txt')\n"",82 ""dataset.name(\""_cbas\"")""83 ]84 },85 {86 ""cell_type"": ""code"",87 ""execution_count"": 9,88 ""metadata"": {89 ""collapsed"": true90 },91 ""outputs"": [],92 ""source"": [93 ""dataset.dump_json(\""_cbas\"")""94 ]95 },96 {97 ""cell_type"": ""code"",98 ""execution_count"": null,99 ""metadata"": {100 ""collapsed"": true101 },102 ""outputs"": [],103 ""source"": [104 ""ll\n""105 ]106 }107 ],108 ""metadata"": {109 ""anaconda-cloud"": {},110 ""kernelspec"": {111 ""display_name"": ""Python [default]"",112 ""language"": ""python"",113 ""name"": ""python2""114 },115 ""language_info"": {116 ""codemirror_mode"": {117 ""name"": ""ipython"",118 ""version"": 2119 },120 ""file_extension"": "".py"",121 ""mimetype"": ""text/x-python"",122 ""name"": ""python"",123 ""nbconvert_exporter"": ""python"",124 ""pygments_lexer"": ""ipython2"",125 ""version"": ""2.7.12""126 }127 },128 ""nbformat"": 4,129 ""nbformat_minor"": 1130}131","Unknown"
132"Homo-Lumo","Mishima-syk/psikit","setup.py",".py","740","25","#!/usr/bin/env python133 134#from distutils.core import setup135from setuptools import setup136 137with open('README.md') as readme_file:138 long_description = readme_file.read()139 140setup(name='Psikit',141 version='0.2.0',142 description='A thin wrapper library for Psi4 and RDKit',143 long_description=long_description,144 long_description_content_type='text/markdown',145 author='Kazufumi Ohkawa, Takayuki Serizawa',146 author_email='kerolinq@gmail.com, seritaka@gmail.com',147 url='https://github.com/Mishima-syk/psikit',148 packages=['psikit'],149 install_requires=['debtcollector'],150 license='MIT',151 classifiers = [152 ""Programming Language :: Python"",153 ""Programming Language :: Python :: 3""]154 )155 156","Python"
157"Homo-Lumo","Mishima-syk/psikit","psikit/pymol_helper.py",".py","3380","80","from glob import glob158import os159 160def run_pymol_server(tmpdir, target='FRONTIER', maprange=0.05):161 '''162 To use the function, user need to install pymol and run the pymol for server mode163 The command is pymol -R164 target ['ESP', 'FRONTIER', 'DUAL']165 '''166 targetlist = ['ESP', 'FRONTIER', 'DUAL']167 if target not in targetlist:168 raise Exception(f'Please set target from ESP, FRONTIER, DENSITY!!')169 import sys170 import xmlrpc.client as xmlrpc171 srv = xmlrpc.ServerProxy('http://localhost:9123')172 srv.do('delete *')173 srv.do('load '+os.path.join(tmpdir, 'target.mol'))174 srv.do('as sticks, target')175 if target == 'FRONTIER':176 homof = glob(os.path.join(tmpdir, 'Psi*_HOMO.cube'))[0]177 lumof = glob(os.path.join(tmpdir, 'Psi*_LUMO.cube'))[0]178 srv.do('load ' + homof + ',HOMO')179 srv.do('load ' + lumof + ',LUMO')180 srv.do(f'isosurface HOMO_A, HOMO, -0.02')181 srv.do(f'isosurface HOMO_B, HOMO, 0.02')182 srv.do(f'isosurface LUMO_A, LUMO, -0.02')183 srv.do(f'isosurface LUMO_B, LUMO, 0.02')184 srv.do('color blue, HOMO_A')185 srv.do('color red, HOMO_B')186 srv.do('color green, LUMO_A')187 srv.do('color yellow, LUMO_B')188 srv.do('set transparency, 0.2')189 srv.do('disable HOMO_A')190 srv.do('disable HOMO_B')191 abb = 'frontier_'192 elif target == 'ESP':193 srv.do('load '+ 'ESP.cube' + ', ESP')194 srv.do('show surface, target')195 srv.do(f'ramp_new cmap, ESP, [-{maprange}, {maprange}]')196 srv.do('color cmap, target')197 srv.do('set transparency, 0.2')198 abb = 'esp_'199 elif target == 'DUAL':200 dualf = glob.glob('DUAL*.cube')[0]201 srv.do('load '+ dualf + ', DUAL_DESC')202 srv.do('show surface, target')203 srv.do(f'ramp_new cmap, DUAL_DESC, [-{maprange}, {maprange}]')204 srv.do('color cmap, target')205 srv.do('set transparency, 0.2')206 abb = 'dual_'207 208 outputpath = abb + 'mo.pse'209 srv.do(f'save {outputpath}')210 print('finished !')211 212 213def save_pyscript(tmpdir, isotype=""isosurface""):214 homof = glob(os.path.join(tmpdir, 'Psi*_HOMO.cube'))[0]215 lumof = glob(os.path.join(tmpdir, 'Psi*_LUMO.cube'))[0]216 with open(""frontier.py"", ""w"") as f:217 f.write('from pymol import *\n')218 f.write('cmd.load(""{0}"", ""HOMO"")\n'.format(homof))219 f.write('cmd.load(""{0}"", ""LUMO"")\n'.format(lumof))220 f.write('cmd.load(""{0}"")\n'.format(os.path.join(tmpdir, ""target.mol"")))221 if isotype == ""isomesh"": 222 f.write('cmd.isomesh(""HOMO_A"", ""HOMO"", -0.02)\n')223 f.write('cmd.isomesh(""HOMO_B"", ""HOMO"", 0.02)\n')224 f.write('cmd.isomesh(""LUMO_A"", ""LUMO"", 0.02)\n')225 f.write('cmd.isomesh(""LUMO_B"", ""LUMO"", -0.02)\n')226 else:227 f.write('cmd.isosurface(""HOMO_A"", ""HOMO"", -0.02)\n')228 f.write('cmd.isosurface(""HOMO_B"", ""HOMO"", 0.02)\n')229 f.write('cmd.isosurface(""LUMO_A"", ""LUMO"", 0.02)\n')230 f.write('cmd.isosurface(""LUMO_B"", ""LUMO"", -0.02)\n') 231 f.write('cmd.color(""blue"", ""HOMO_A"")\n') 232 f.write('cmd.color(""red"", ""HOMO_B"")\n') 233 f.write('cmd.color(""blue"", ""LUMO_A"")\n')234 f.write('cmd.color(""red"", ""LUMO_B"")\n')235 f.write('cmd.disable(""LUMO_A"")\n') 236 f.write('cmd.disable(""LUMO_B"")\n') ","Python"
237"Homo-Lumo","Mishima-syk/psikit","psikit/sapt.py",".py","11353","254","from rdkit import Chem238from psikit import Psikit239import numpy as np240import os241from .util import mol2xyz242 243class Sapt():244 def __init__(self, threads=4, memory=4, debug=False):245 import psi4246 from . import helper_SAPT247 248 self.psi4 = psi4249 self.helper_SAPT = helper_SAPT250 self.psi4.set_memory(""{} GB"".format(memory))251 #self.psi4.set_options({""save_jk"": True}) # for JK calculation252 self.psi4.set_num_threads(threads)253 self.wfn = None254 self.monomer1 = None255 self.monomer2 = None256 self.dimer = None257 self.debug = debug258 if self.debug:259 self.psi4.core.set_output_file(""psikit_out.dat"", True)260 else:261 self.psi4.core.be_quiet()262 263 def monomer1_from_molfile(self, molfile, opt=True, removeHs=False):264 self.monomer1 = Chem.MolFromMolFile(molfile, removeHs=removeHs)265 266 def monomer2_from_molfile(self, molfile, opt=True, removeHs=False):267 self.monomer2 = Chem.MolFromMolFile(molfile, removeHs=removeHs)268 269 def make_dimer(self):270 xyz1 = mol2xyz(self.monomer1)271 xyz2 = mol2xyz(self.monomer2)272 self.dimer = ""{}--\n{}"".format(xyz1, xyz2)273 self.dimer += ""no_reorient\n""274 self.dimer += ""no_com\n""275 self.dimer += ""units angstrom\n""276 277 278 def run_sapt(self, basis='aug-cc-pvdz', e_convergence=1e-8, d_convergence=1e-8, memory=4):279 self.psi4.set_options({'basis':basis,280 'e_convergence':e_convergence,281 'd_convergence':d_convergence})282 dimer = self.psi4.geometry(self.dimer)283 sapt = self.helper_SAPT.helper_SAPT(dimer, memory=memory)284 ### Start E100 Electostatics285 elst_timer = self.helper_SAPT.sapt_timer('electrostatics')286 Elst10 = 4 * np.einsum('abab', sapt.vt('abab'))287 elst_timer.stop()288 ### End E100 Electrostatics289 290 ### Start E100 Exchange291 exch_timer =self.helper_SAPT.sapt_timer('exchange')292 vt_abba = sapt.vt('abba')293 vt_abaa = sapt.vt('abaa')294 vt_abbb = sapt.vt('abbb')295 vt_abab = sapt.vt('abab')296 s_ab = sapt.s('ab')297 298 Exch100 = np.einsum('abba', vt_abba)299 300 tmp = 2 * vt_abaa - vt_abaa.swapaxes(2, 3)301 Exch100 += np.einsum('Ab,abaA', s_ab, tmp)302 303 tmp = 2 * vt_abbb - vt_abbb.swapaxes(2, 3)304 Exch100 += np.einsum('Ba,abBb', s_ab.T, tmp)305 306 Exch100 -= 2 * np.einsum('Ab,BA,abaB', s_ab, s_ab.T, vt_abab)307 Exch100 -= 2 * np.einsum('AB,Ba,abAb', s_ab, s_ab.T, vt_abab)308 Exch100 += np.einsum('Ab,Ba,abAB', s_ab, s_ab.T, vt_abab)309 310 Exch100 *= -2311 exch_timer.stop()312 ### End E100 (S^2) Exchange313 314 ### Start E200 Disp315 disp_timer = self.helper_SAPT.sapt_timer('dispersion')316 v_abrs = sapt.v('abrs')317 v_rsab = sapt.v('rsab')318 e_rsab = 1/(-sapt.eps('r', dim=4) - sapt.eps('s', dim=3) + sapt.eps('a', dim=2) + sapt.eps('b'))319 320 Disp200 = 4 * np.einsum('rsab,rsab,abrs->', e_rsab, v_rsab, v_abrs)321 ### End E200 Disp322 323 ### Start E200 Exchange-Dispersion324 325 # Build t_rsab326 t_rsab = np.einsum('rsab,rsab->rsab', v_rsab, e_rsab)327 328 # Build h_abrs329 vt_abar = sapt.vt('abar')330 vt_abra = sapt.vt('abra')331 vt_absb = sapt.vt('absb')332 vt_abbs = sapt.vt('abbs')333 334 tmp = 2 * vt_abar - vt_abra.swapaxes(2, 3)335 h_abrs = np.einsum('as,AbAr->abrs', sapt.s('as'), tmp)336 337 tmp = 2 * vt_abra - vt_abar.swapaxes(2, 3)338 h_abrs += np.einsum('As,abrA->abrs', sapt.s('as'), tmp)339 340 tmp = 2 * vt_absb - vt_abbs.swapaxes(2, 3)341 h_abrs += np.einsum('br,aBsB->abrs', sapt.s('br'), tmp)342 343 tmp = 2 * vt_abbs - vt_absb.swapaxes(2, 3)344 h_abrs += np.einsum('Br,abBs->abrs', sapt.s('br'), tmp)345 346 # Build q_abrs347 vt_abas = sapt.vt('abas')348 q_abrs = np.einsum('br,AB,aBAs->abrs', sapt.s('br'), sapt.s('ab'), vt_abas)349 q_abrs -= 2 * np.einsum('Br,AB,abAs->abrs', sapt.s('br'), sapt.s('ab'), vt_abas)350 q_abrs -= 2 * np.einsum('br,aB,ABAs->abrs', sapt.s('br'), sapt.s('ab'), vt_abas)351 q_abrs += 4 * np.einsum('Br,aB,AbAs->abrs', sapt.s('br'), sapt.s('ab'), vt_abas)352 353 vt_abrb = sapt.vt('abrb')354 q_abrs -= 2 * np.einsum('as,bA,ABrB->abrs', sapt.s('as'), sapt.s('ba'), vt_abrb)355 q_abrs += 4 * np.einsum('As,bA,aBrB->abrs', sapt.s('as'), sapt.s('ba'), vt_abrb)356 q_abrs += np.einsum('as,BA,AbrB->abrs', sapt.s('as'), sapt.s('ba'), vt_abrb)357 q_abrs -= 2 * np.einsum('As,BA,abrB->abrs', sapt.s('as'), sapt.s('ba'), vt_abrb)358 359 vt_abab = sapt.vt('abab')360 q_abrs += np.einsum('Br,As,abAB->abrs', sapt.s('br'), sapt.s('as'), vt_abab)361 q_abrs -= 2 * np.einsum('br,As,aBAB->abrs', sapt.s('br'), sapt.s('as'), vt_abab)362 q_abrs -= 2 * np.einsum('Br,as,AbAB->abrs', sapt.s('br'), sapt.s('as'), vt_abab)363 364 vt_abrs = sapt.vt('abrs')365 q_abrs += np.einsum('bA,aB,ABrs->abrs', sapt.s('ba'), sapt.s('ab'), vt_abrs)366 q_abrs -= 2 * np.einsum('bA,AB,aBrs->abrs', sapt.s('ba'), sapt.s('ab'), vt_abrs)367 q_abrs -= 2 * np.einsum('BA,aB,Abrs->abrs', sapt.s('ba'), sapt.s('ab'), vt_abrs)368 369 # Sum it all together370 xd_absr = sapt.vt('absr')371 xd_absr += h_abrs.swapaxes(2, 3)372 xd_absr += q_abrs.swapaxes(2, 3)373 ExchDisp20 = -2 * np.einsum('absr,rsab->', xd_absr, t_rsab)374 375 disp_timer.stop()376 ### End E200 Exchange-Dispersion377 378 379 ### Start E200 Induction and Exchange-Induction380 381 # E200Induction and CPHF orbitals382 ind_timer = self.helper_SAPT.sapt_timer('induction')383 384 CPHF_ra, Ind20_ba = sapt.chf('B', ind=True)385 self.helper_SAPT.sapt_printer('Ind20,r (A<-B)', Ind20_ba)386 387 CPHF_sb, Ind20_ab = sapt.chf('A', ind=True)388 self.helper_SAPT.sapt_printer('Ind20,r (A->B)', Ind20_ab)389 390 Ind20r = Ind20_ba + Ind20_ab391 392 # Exchange-Induction393 394 # A <- B395 vt_abra = sapt.vt('abra')396 vt_abar = sapt.vt('abar')397 ExchInd20_ab = np.einsum('ra,abbr', CPHF_ra, sapt.vt('abbr'))398 ExchInd20_ab += 2 * np.einsum('rA,Ab,abar', CPHF_ra, sapt.s('ab'), vt_abar)399 ExchInd20_ab += 2 * np.einsum('ra,Ab,abrA', CPHF_ra, sapt.s('ab'), vt_abra)400 ExchInd20_ab -= np.einsum('rA,Ab,abra', CPHF_ra, sapt.s('ab'), vt_abra)401 402 vt_abbb = sapt.vt('abbb')403 vt_abab = sapt.vt('abab')404 ExchInd20_ab -= np.einsum('ra,Ab,abAr', CPHF_ra, sapt.s('ab'), vt_abar)405 ExchInd20_ab += 2 * np.einsum('ra,Br,abBb', CPHF_ra, sapt.s('br'), vt_abbb)406 ExchInd20_ab -= np.einsum('ra,Br,abbB', CPHF_ra, sapt.s('br'), vt_abbb)407 ExchInd20_ab -= 2 * np.einsum('rA,Ab,Br,abaB', CPHF_ra, sapt.s('ab'), sapt.s('br'), vt_abab)408 409 vt_abrb = sapt.vt('abrb')410 ExchInd20_ab -= 2 * np.einsum('ra,Ab,BA,abrB', CPHF_ra, sapt.s('ab'), sapt.s('ba'), vt_abrb)411 ExchInd20_ab -= 2 * np.einsum('ra,AB,Br,abAb', CPHF_ra, sapt.s('ab'), sapt.s('br'), vt_abab)412 ExchInd20_ab -= 2 * np.einsum('rA,AB,Ba,abrb', CPHF_ra, sapt.s('ab'), sapt.s('ba'), vt_abrb)413 414 ExchInd20_ab += np.einsum('ra,Ab,Br,abAB', CPHF_ra, sapt.s('ab'), sapt.s('br'), vt_abab)415 ExchInd20_ab += np.einsum('rA,Ab,Ba,abrB', CPHF_ra, sapt.s('ab'), sapt.s('ba'), vt_abrb)416 417 ExchInd20_ab *= -2418 self.helper_SAPT.sapt_printer('Exch-Ind20,r (A<-B)', ExchInd20_ab)419 420 # B <- A421 vt_abbs = sapt.vt('abbs')422 vt_absb = sapt.vt('absb')423 ExchInd20_ba = np.einsum('sb,absa', CPHF_sb, sapt.vt('absa'))424 ExchInd20_ba += 2 * np.einsum('sB,Ba,absb', CPHF_sb, sapt.s('ba'), vt_absb)425 ExchInd20_ba += 2 * np.einsum('sb,Ba,abBs', CPHF_sb, sapt.s('ba'), vt_abbs)426 ExchInd20_ba -= np.einsum('sB,Ba,abbs', CPHF_sb, sapt.s('ba'), vt_abbs)427 428 vt_abaa = sapt.vt('abaa')429 vt_abab = sapt.vt('abab')430 ExchInd20_ba -= np.einsum('sb,Ba,absB', CPHF_sb, sapt.s('ba'), vt_absb)431 ExchInd20_ba += 2 * np.einsum('sb,As,abaA', CPHF_sb, sapt.s('as'), vt_abaa)432 ExchInd20_ba -= np.einsum('sb,As,abAa', CPHF_sb, sapt.s('as'), vt_abaa)433 ExchInd20_ba -= 2 * np.einsum('sB,Ba,As,abAb', CPHF_sb, sapt.s('ba'), sapt.s('as'), vt_abab)434 435 vt_abas = sapt.vt('abas')436 ExchInd20_ba -= 2 * np.einsum('sb,Ba,AB,abAs', CPHF_sb, sapt.s('ba'), sapt.s('ab'), vt_abas)437 ExchInd20_ba -= 2 * np.einsum('sb,BA,As,abaB', CPHF_sb, sapt.s('ba'), sapt.s('as'), vt_abab)438 ExchInd20_ba -= 2 * np.einsum('sB,BA,Ab,abas', CPHF_sb, sapt.s('ba'), sapt.s('ab'), vt_abas)439 440 ExchInd20_ba += np.einsum('sb,Ba,As,abAB', CPHF_sb, sapt.s('ba'), sapt.s('as'), vt_abab)441 ExchInd20_ba += np.einsum('sB,Ba,Ab,abAs', CPHF_sb, sapt.s('ba'), sapt.s('ab'), vt_abas)442 443 ExchInd20_ba *= -2444 self.helper_SAPT.sapt_printer('Exch-Ind20,r (A->B)', ExchInd20_ba)445 ExchInd20r = ExchInd20_ba + ExchInd20_ab446 447 ind_timer.stop()448 ### End E200 Induction and Exchange-Induction449 450 print('\nSAPT0 Results')451 print('-' * 70)452 self.helper_SAPT.sapt_printer('Exch10 (S^2)', Exch100)453 self.helper_SAPT.sapt_printer('Elst10', Elst10)454 self.helper_SAPT.sapt_printer('Disp20', Disp200)455 self.helper_SAPT.sapt_printer('Exch-Disp20', ExchDisp20)456 self.helper_SAPT.sapt_printer('Ind20,r', Ind20r)457 self.helper_SAPT.sapt_printer('Exch-Ind20,r', ExchInd20r)458 459 print('-' * 70)460 sapt0 = Exch100 + Elst10 + Disp200 + ExchDisp20 + Ind20r + ExchInd20r461 self.helper_SAPT.sapt_printer('Total SAPT0', sapt0) 462 return sapt0, Exch100, Elst10, Disp200, ExchDisp20, Ind20r, ExchInd20r463 464 def run_fisapt(self, basis='jun-cc-pvdz', scf_type='df', d_convergence=1e-8, memory=4, fisapt_path='fsapt/', return_wfn=False):465 import shutil466 from distutils.dir_util import copy_tree467 from . import fsapt_helper468 self.psi4.set_options({'basis':basis,469 'scf_type':scf_type,470 'd_convergence':d_convergence,471 'FISAPT_FSAPT_FILEPATH':fisapt_path,472 'FISAPT_DO_PLOT': 'true'473 })474 self.psi4.geometry(self.dimer)475 res = self.psi4.energy('fisapt0', return_wfn=return_wfn)476 copy_tree(self.psi4.core.get_datadir()+'/fsapt', fisapt_path)477 #sapt = self.helper_SAPT.helper_SAPT(dimer, memory=memory)478 feats1 = fsapt_helper.make_feat_data(self.monomer1, 1)479 feats2 = fsapt_helper.make_feat_data(self.monomer2, 1 + self.monomer1.GetNumAtoms())480 with open(os.path.join(fisapt_path, 'fA.dat'), 'w') as fA:481 for feat1 in feats1:482 fA.write(' '.join(feat1) + '\n')483 484 with open(os.path.join(fisapt_path,'fB.dat'), 'w') as fB:485 for feat2 in feats2:486 fB.write(' '.join(feat2) + '\n')487 488 return res489 490","Python"
491"Homo-Lumo","Mishima-syk/psikit","psikit/psikit.py",".py","8801","242","# -*- coding: utf-8 -*-492import rdkit493from rdkit import Chem494from rdkit.Chem import AllChem495import numpy as np496import glob497import os498import uuid499import warnings500from collections import defaultdict501from tempfile import mkdtemp502from shutil import rmtree503from debtcollector import moves504from .util import mol2xyz505from .pymol_helper import run_pymol_server, save_pyscript506warnings.simplefilter('ignore')507 508 509class Psikit(object):510 def __init__(self, threads=4, memory=4, debug=False):511 import psi4512 self.psi4 = psi4513 self.psi4.set_memory(""{} GB"".format(memory))514 #self.psi4.set_options({""save_jk"": True}) # for JK calculation515 self.psi4.set_num_threads(threads)516 self.wfn = None517 self.mol = None518 self.debug = debug519 self.tempdir = mkdtemp()520 if self.debug:521 self.psi4.core.set_output_file(""psikit_out.dat"", True)522 else:523 self.psi4.core.be_quiet()524 525 def clean(self):526 rmtree(self.tempdir)527 528 def read_from_smiles(self, smiles_str, opt=True):529 self.mol = Chem.MolFromSmiles(smiles_str)530 if opt:531 self.rdkit_optimize() 532 533 def read_from_molfile(self, molfile, opt=True, removeHs=False):534 self.mol = Chem.MolFromMolFile(molfile, removeHs=removeHs)535 if opt:536 self.rdkit_optimize() 537 538 def rdkit_optimize(self, addHs=True):539 if addHs:540 self.mol = Chem.AddHs(self.mol)541 AllChem.EmbedMolecule(self.mol, useExpTorsionAnglePrefs=True,useBasicKnowledge=True)542 AllChem.UFFOptimizeMolecule(self.mol)543 544 def geometry(self, multiplicity=1):545 xyz = self.mol2xyz(multiplicity=multiplicity)546 self.psi4.geometry(xyz)547 548 def energy(self, basis_sets= ""scf/6-31g**"", return_wfn=True, multiplicity=1):549 self.geometry(multiplicity=multiplicity)550 scf_energy, wfn = self.psi4.energy(basis_sets, return_wfn=return_wfn)551 self.psi4.core.clean()552 self.wfn = wfn553 self.mol = self.xyz2mol()554 return scf_energy555 556 def optimize(self, basis_sets= ""scf/6-31g**"", return_wfn=True, name=None, multiplicity=1, maxiter=50):557 if not name:558 name = uuid.uuid4().hex559 self.psi4.core.IO.set_default_namespace(name)560 self.geometry(multiplicity=multiplicity)561 self.psi4.set_options({'GEOM_MAXITER':maxiter})562 try:563 scf_energy, wfn = self.psi4.optimize(basis_sets, return_wfn=return_wfn)564 self.wfn = wfn565 except self.psi4.OptimizationConvergenceError as cError:566 print('Convergence error caught: {0}'.format(cError))567 self.wfn = cError.wfn568 scf_energy = self.wfn.energy()569 self.psi4.core.clean()570 self.mol = self.xyz2mol()571 572 if not self.debug:573 self.psi4.core.opt_clean() # Seg fault will occured when the function is called before optimize.574 return scf_energy575 576 def set_options(self, **kwargs):577 """"""578 http://www.psicode.org/psi4manual/1.2/psiapi.html579 IV. Analysis of Intermolecular Interactions580 and 581 http://forum.psicode.org/t/how-can-i-change-max-iteration-in-energy-method/1238/2582 """"""583 self.psi4.set_options(kwargs)584 585 def mol2xyz(self, multiplicity=1):586 return mol2xyz(self.mol)587 588 def xyz2mol(self, confId=0):589 natom = self.wfn.molecule().natom()590 mol_array_bohr = self.wfn.molecule().geometry().to_array()591 mol_array = mol_array_bohr * 0.52917721092592 nmol = Chem.Mol(self.mol)593 conf = nmol.GetConformer(confId)594 for i in range(natom):595 conf.SetAtomPosition(i, tuple(mol_array[i]))596 return nmol597 598 599 def clone_mol(self):600 return Chem.Mol(self.mol)601 602 def create_cube_files(self, gridspace=0.3):603 if self.wfn == None:604 print('wfn not found. run optimze()/energy()')605 return None606 else:607 a = self.wfn.nalpha() # HOMO608 b = a + 1 # LUMO609 self.psi4.set_options({""cubeprop_tasks"": ['ESP', 'FRONTIER_ORBITALS', 'Density', 'DUAL_DESCRIPTOR'],610 ""cubic_grid_spacing"": [gridspace, gridspace, gridspace],611 ""cubeprop_filepath"": self.tempdir612 })613 Chem.MolToMolFile(self.mol, os.path.join(self.tempdir, 'target.mol'))614 self.psi4.cubeprop(self.wfn)615 616 getMOview = moves.moved_function(create_cube_files, 'getMOview', __name__)617 618 def view_on_pymol(self, target='FRONTIER', maprange=0.05, gridspace=0.3):619 self.create_cube_files(gridspace=gridspace)620 run_pymol_server(self.tempdir, target=target, maprange=maprange)621 622 def save_frontier(self, gridspace=0.3, isotype=""isosurface""):623 self.create_cube_files(gridspace=gridspace)624 save_pyscript(self.tempdir, isotype=isotype) 625 626 def save_fchk(self, filename=""output.fchk""):627 fchk_writer = self.psi4.core.FCHKWriter(self.wfn)628 fchk_writer.write(filename)629 630 def save_cube(self):631 self.psi4.cubeprop(self.wfn)632 633 def calc_resp_charges(self, constrain_symmetric_atoms=False):634 if self.wfn.molecule() == None:635 print('please run optimze() at first!')636 return None637 try:638 import resp639 except:640 print('please install resp at first')641 print('conda install -c psi4 resp')642 return None643 # https://www.cgl.ucsf.edu/chimerax/docs/user/radii.html644 options = {'VDW_SCALE_FACTORS' : [1.4, 1.6, 1.8, 2.0],645 'VDW_POINT_DENSITY' : 1.0,646 'RESP_A' : 0.0005,647 'RESP_B' : 0.1,648 'RESTRAINT' : True,649 'RADIUS' : {'Br':1.98, 'I':2.09,}650 }651 652 if constrain_symmetric_atoms:653 ranks = Chem.CanonicalRankAtoms(self.mol, breakTies=False)654 groups = defaultdict(list)655 for idx, rank in enumerate(ranks):656 groups[rank].append(idx + 1) # as RESP atoms are 1-indexed but RDKit 0-indexed657 constraint_groups = [constraint_group for constraint_group in groups.values()]658 options['CONSTRAINT_GROUP'] = constraint_groups659 660 charges = resp.resp([self.wfn.molecule()], options)661 #breakpoint()662 663 options['resp_a'] = 0.001664 resp.set_stage2_constraint(self.wfn.molecule(), charges[1], options)665 options['grid']=['%i_%s_grid.dat'%(1, self.wfn.molecule().name())]666 options['esp']=['%i_%s_grid_esp.dat'%(1, self.wfn.molecule().name())]667 668 charges2 = resp.resp([self.wfn.molecule()], options)669 670 atoms = self.mol.GetAtoms()671 for idx, atom in enumerate(atoms):672 atom.SetDoubleProp(""EP"", charges2[0][idx])673 atom.SetDoubleProp(""RESP"", charges2[1][idx])674 return charges2[1]675 676 677 def calc_mulliken_charges(self):678 '''679 Compute Mulliken Charges680 And return the results as numpy array.681 '''682 if self.wfn.molecule() == None:683 print('please run optimze() at first!')684 return None685 self.psi4.oeprop(self.wfn, 'MULLIKEN_CHARGES')686 mulliken_acp = self.wfn.atomic_point_charges()687 atoms = self.mol.GetAtoms()688 for idx, atom in enumerate(atoms):689 atom.SetDoubleProp(""MULLIKEN"", mulliken_acp.np[idx])690 return mulliken_acp.np691 692 def calc_lowdin_charges(self):693 '''694 Compute Lowdin Charges695 And return the results as numpy array.696 '''697 if self.wfn.molecule() == None:698 print('please run optimze() at first!')699 return None700 self.psi4.oeprop(self.wfn, 'LOWDIN_CHARGES')701 lowdin_acp = self.wfn.atomic_point_charges()702 atoms = self.mol.GetAtoms()703 for idx, atom in enumerate(atoms):704 atom.SetDoubleProp(""LOWDIN"", lowdin_acp.np[idx])705 return lowdin_acp.np706 707 708 @property709 def dipolemoment(self, basis_sets=""scf/6-31g**"", return_wfn=True):710 # The three components of the SCF dipole [Debye]711 x = self.psi4.get_variable('SCF DIPOLE X')712 y = self.psi4.get_variable('SCF DIPOLE Y')713 z = self.psi4.get_variable('SCF DIPOLE Z')714 total = np.sqrt(x * x + y * y + z * z)715 return (x, y, z, total)716 717 @property718 def HOMO(self):719 return self.wfn.epsilon_a_subset('AO', 'ALL').np[self.wfn.nalpha()-1]720 721 @property722 def LUMO(self):723 return self.wfn.epsilon_a_subset('AO', 'ALL').np[self.wfn.nalpha()]724 725 @property726 def coulomb_matrix(self):727 return self.wfn.jk().J[0].to_array()728 729 @property730 def exchange_matrix(self):731 return self.wfn.jk().K[0].to_array()732","Python"
733"Homo-Lumo","Mishima-syk/psikit","psikit/__init__.py",".py","75","5","# -*- coding: utf-8 -*-734from .psikit import Psikit735from .sapt import Sapt736 737","Python"
738"Homo-Lumo","Mishima-syk/psikit","psikit/util.py",".py","356","9","from rdkit import Chem739 740def mol2xyz(mol, multiplicity=1):741 charge = Chem.GetFormalCharge(mol)742 xyz_string = ""\n{} {}\n"".format(charge, multiplicity)743 for atom in mol.GetAtoms():744 pos = mol.GetConformer().GetAtomPosition(atom.GetIdx())745 xyz_string += ""{} {} {} {}\n"".format(atom.GetSymbol(), pos.x, pos.y, pos.z)746 return xyz_string","Python"
747"Homo-Lumo","Mishima-syk/psikit","psikit/helper_SAPT.py",".py","16733","442","""""""748Helper classes and functions for the SAPT directory.749 750References:751- Equations and algorithms from [Szalewicz:2005:43], [Jeziorski:1994:1887],752[Szalewicz:2012:254], and [Hohenstein:2012:304]753""""""754 755__authors__ = ""Daniel G. A. Smith""756__credits__ = [""Daniel G. A. Smith""]757 758__copyright__ = ""(c) 2014-2018, The Psi4NumPy Developers""759__license__ = ""BSD-3-Clause""760__date__ = ""2015-12-01""761 762import numpy as np763import time764import psi4765 766class helper_SAPT(object):767 768 def __init__(self, dimer, memory=8, algorithm='MO', reference='RHF'):769 print(""\nInitializing SAPT object...\n"")770 tinit_start = time.time()771 772 # Set a few crucial attributes773 self.alg = algorithm.upper()774 self.reference = reference.upper()775 dimer.reset_point_group('c1')776 dimer.fix_orientation(True)777 dimer.fix_com(True)778 dimer.update_geometry()779 nfrags = dimer.nfragments()780 if nfrags != 2:781 psi4.core.clean()782 raise Exception(""Found %d fragments, must be 2."" % nfrags)783 784 # Grab monomers in DCBS785 monomerA = dimer.extract_subsets(1, 2)786 monomerA.set_name('monomerA')787 monomerB = dimer.extract_subsets(2, 1)788 monomerB.set_name('monomerB')789 self.mult_A = monomerA.multiplicity()790 self.mult_B = monomerB.multiplicity()791 792 # Compute monomer properties793 794 tstart = time.time()795 self.rhfA, self.wfnA = psi4.energy('SCF', return_wfn=True, molecule=monomerA)796 self.V_A = np.asarray(psi4.core.MintsHelper(self.wfnA.basisset()).ao_potential())797 print(""RHF for monomer A finished in %.2f seconds."" % (time.time() - tstart))798 799 tstart = time.time()800 self.rhfB, self.wfnB = psi4.energy('SCF', return_wfn=True, molecule=monomerB)801 self.V_B = np.asarray(psi4.core.MintsHelper(self.wfnB.basisset()).ao_potential())802 print(""RHF for monomer B finished in %.2f seconds."" % (time.time() - tstart))803 804 # Setup a few variables805 self.memory = memory806 self.nmo = self.wfnA.nmo()807 808 # Monomer A809 self.nuc_rep_A = monomerA.nuclear_repulsion_energy()810 self.ndocc_A = self.wfnA.doccpi()[0]811 self.nvirt_A = self.nmo - self.ndocc_A812 if reference == 'ROHF':813 self.idx_A = ['i', 'a', 'r']814 self.nsocc_A = self.wfnA.soccpi()[0]815 occA = self.ndocc_A + self.nsocc_A816 else:817 self.idx_A = ['a', 'r']818 self.nsocc_A = 0819 occA = self.ndocc_A 820 821 self.C_A = np.asarray(self.wfnA.Ca())822 self.Co_A = self.C_A[:, :self.ndocc_A]823 self.Ca_A = self.C_A[:, self.ndocc_A:occA]824 self.Cv_A = self.C_A[:, occA:]825 self.eps_A = np.asarray(self.wfnA.epsilon_a())826 827 # Monomer B828 self.nuc_rep_B = monomerB.nuclear_repulsion_energy()829 self.ndocc_B = self.wfnB.doccpi()[0]830 self.nvirt_B = self.nmo - self.ndocc_B831 if reference == 'ROHF':832 self.idx_B = ['j', 'b', 's']833 self.nsocc_B = self.wfnB.soccpi()[0]834 occB = self.ndocc_B + self.nsocc_B835 else:836 self.idx_B = ['b', 's']837 self.nsocc_B = 0838 occB = self.ndocc_B 839 840 self.C_B = np.asarray(self.wfnB.Ca())841 self.Co_B = self.C_B[:, :self.ndocc_B]842 self.Ca_B = self.C_B[:, self.ndocc_B:occB]843 self.Cv_B = self.C_B[:, occB:]844 self.eps_B = np.asarray(self.wfnB.epsilon_a())845 846 # Dimer847 self.nuc_rep = dimer.nuclear_repulsion_energy() - self.nuc_rep_A - self.nuc_rep_B848 self.vt_nuc_rep = self.nuc_rep / ((2 * self.ndocc_A + self.nsocc_A)849 * (2 * self.ndocc_B + self.nsocc_B))850 851 # Make slice, orbital, and size dictionaries852 if reference == 'ROHF':853 self.slices = {854 'i': slice(0, self.ndocc_A),855 'a': slice(self.ndocc_A, occA),856 'r': slice(occA, None),857 'j': slice(0, self.ndocc_B),858 'b': slice(self.ndocc_B, occB),859 's': slice(occB, None)860 }861 862 self.orbitals = {'i': self.Co_A,863 'a': self.Ca_A,864 'r': self.Cv_A,865 'j': self.Co_B,866 'b': self.Ca_B,867 's': self.Cv_B868 }869 870 self.sizes = {'i': self.ndocc_A,871 'a': self.nsocc_A,872 'r': self.nvirt_A,873 'j': self.ndocc_B,874 'b': self.nsocc_B,875 's': self.nvirt_B}876 877 else:878 self.slices = {879 'a': slice(0, self.ndocc_A),880 'r': slice(occA, None),881 'b': slice(0, self.ndocc_B),882 's': slice(occB, None)883 }884 885 self.orbitals = {'a': self.Co_A,886 'r': self.Cv_A,887 'b': self.Co_B,888 's': self.Cv_B889 }890 891 self.sizes = {'a': self.ndocc_A,892 'r': self.nvirt_A,893 'b': self.ndocc_B,894 's': self.nvirt_B}895 896 # Compute size of ERI tensor in GB897 self.dimer_wfn = psi4.core.Wavefunction.build(dimer, psi4.core.get_global_option('BASIS'))898 mints = psi4.core.MintsHelper(self.dimer_wfn.basisset())899 self.mints = mints900 ERI_Size = (self.nmo ** 4) * 8.e-9901 memory_footprint = ERI_Size * 4902 if memory_footprint > self.memory:903 psi4.core.clean()904 raise Exception(""Estimated memory utilization (%4.2f GB) exceeds numpy_memory \905 limit of %4.2f GB."" % (memory_footprint, self.memory))906 907 # Integral generation from Psi4's MintsHelper908 print('Building ERI tensor...')909 tstart = time.time()910 # Leave ERI as a Psi4 Matrix911 self.I = np.asarray(self.mints.ao_eri()).swapaxes(1,2)912 print('...built ERI tensor in %.3f seconds.' % (time.time() - tstart))913 print(""Size of the ERI tensor is %4.2f GB, %d basis functions."" % (ERI_Size, self.nmo))914 self.S = np.asarray(self.mints.ao_overlap())915 916 # Save additional rank 2 tensors917 self.V_A_BB = np.einsum('ui,vj,uv->ij', self.C_B, self.C_B, self.V_A)918 self.V_A_AB = np.einsum('ui,vj,uv->ij', self.C_A, self.C_B, self.V_A)919 self.V_B_AA = np.einsum('ui,vj,uv->ij', self.C_A, self.C_A, self.V_B)920 self.V_B_AB = np.einsum('ui,vj,uv->ij', self.C_A, self.C_B, self.V_B)921 922 self.S_AB = np.einsum('ui,vj,uv->ij', self.C_A, self.C_B, self.S)923 924 if self.alg == ""AO"":925 tstart = time.time()926 aux_basis = psi4.core.BasisSet.build(self.dimer_wfn.molecule(), ""DF_BASIS_SCF"",927 psi4.core.get_option(""SCF"", ""DF_BASIS_SCF""),928 ""JKFIT"", psi4.core.get_global_option('BASIS'),929 puream=self.dimer_wfn.basisset().has_puream())930 931 self.jk = psi4.core.JK.build(self.dimer_wfn.basisset(), aux_basis)932 self.jk.set_memory(int(memory * 1e9))933 self.jk.initialize()934 print(""\n...initialized JK objects in %5.2f seconds."" % (time.time() - tstart))935 936 print(""\n...finished initializing SAPT object in %5.2f seconds."" % (time.time() - tinit_start))937 938 # Compute MO ERI tensor (v) on the fly939 def v(self, string):940 if len(string) != 4:941 psi4.core.clean()942 raise Exception('v: string %s does not have 4 elements' % string)943 944 # ERI's from mints are of type (11|22) - need <12|12>945 V = np.einsum('pA,pqrs->Aqrs', self.orbitals[string[0]], self.I)946 V = np.einsum('qB,Aqrs->ABrs', self.orbitals[string[1]], V)947 V = np.einsum('rC,ABrs->ABCs', self.orbitals[string[2]], V)948 V = np.einsum('sD,ABCs->ABCD', self.orbitals[string[3]], V)949 return V950 951 # Grab MO overlap matrices952 def s(self, string):953 if len(string) != 2:954 psi4.core.clean()955 raise Exception('S: string %s does not have 2 elements.' % string)956 957 for alpha in 'ijab':958 if (alpha in string) and (self.sizes[alpha] == 0):959 return np.array([0]).reshape(1,1)960 961 s1 = string[0]962 s2 = string[1]963 964 # Compute on the fly965 return (self.orbitals[string[0]].T).dot(self.S).dot(self.orbitals[string[1]])966 #return np.einsum('ui,vj,uv->ij', self.orbitals[string[0]], self.orbitals[string[1]], self.S)967 968 # Grab epsilons, reshape if requested969 def eps(self, string, dim=1):970 if len(string) != 1:971 psi4.core.clean()972 raise Exception('Epsilon: string %s does not have 1 element.' % string)973 974 shape = (-1,) + tuple([1] * (dim - 1))975 976 if (string == 'i') or (string == 'a') or (string == 'r'):977 return self.eps_A[self.slices[string]].reshape(shape)978 elif (string == 'j') or (string == 'b') or (string == 's'):979 return self.eps_B[self.slices[string]].reshape(shape)980 else:981 psi4.core.clean()982 raise Exception('Unknown orbital type in eps: %s.' % string)983 984 # Grab MO potential matrices985 def potential(self, string, side):986 if len(string) != 2:987 psi4.core.clean()988 raise Exception('Potential: string %s does not have 2 elements.' % string)989 990 s1 = string[0]991 s2 = string[1]992 993 # Two separate cases994 if side == 'A':995 # Compute on the fly996 return (self.orbitals[string[0]].T).dot(self.V_A).dot(self.orbitals[string[1]])997 #return np.einsum('ui,vj,uv->ij', self.orbitals[s1], self.orbitals[s2], self.V_A)998 999 elif side == 'B':1000 # Compute on the fly1001 return (self.orbitals[string[0]].T).dot(self.V_B).dot(self.orbitals[string[1]])1002 #return np.einsum('ui,vj,uv->ij', self.orbitals[s1], self.orbitals[s2], self.V_B)1003 else:1004 psi4.core.clean()1005 raise Exception('helper_SAPT.potential side must be either A or B, not %s.' % side)1006 1007 # Compute V tilde, Index as V_{1,2}^{3,4}1008 def vt(self, string):1009 if len(string) != 4:1010 psi4.core.clean()1011 raise Exception('Compute tilde{V}: string %s does not have 4 elements' % string)1012 1013 for alpha in 'ijab':1014 if (alpha in string) and (self.sizes[alpha] == 0):1015 return np.array([0]).reshape(1,1,1,1)1016 1017 # Grab left and right strings1018 s_left = string[0] + string[2]1019 s_right = string[1] + string[3]1020 1021 # ERI term1022 V = self.v(string)1023 # Potential A1024 S_A = self.s(s_left)1025 V_A = self.potential(s_right, 'A') / (2 * self.ndocc_A + self.nsocc_A)1026 V += np.einsum('ik,jl->ijkl', S_A, V_A)1027 1028 # Potential B1029 S_B = self.s(s_right)1030 V_B = self.potential(s_left, 'B') / (2 * self.ndocc_B + self.nsocc_B)1031 #print s_right, np.abs(V_B).sum()1032 V += np.einsum('ik,jl->ijkl', V_B, S_B)1033 1034 # Nuclear1035 V += np.einsum('ik,jl->ijkl', S_A, S_B) * self.vt_nuc_rep1036 1037 return V1038 1039 # Compute CPHF orbitals1040 def chf(self, monomer, ind=False):1041 if monomer not in ['A', 'B']:1042 psi4.core.clean()1043 raise Exception('%s is not a valid monomer for CHF.' % monomer)1044 1045 if self.reference == 'ROHF':1046 psi4.core.clean()1047 raise Exception('CPHF for a ROHF reference not implemented yet.')1048 1049 if monomer == 'A':1050 # Form electrostatic potential1051 w_n = 2 * np.einsum('saba->bs', self.v('saba'))1052 w_n += self.V_A_BB[self.slices['b'], self.slices['s']]1053 eps_ov = (self.eps('b', dim=2) - self.eps('s'))1054 1055 # Set terms1056 v_term1 = 'sbbs'1057 v_term2 = 'sbsb'1058 no, nv = self.ndocc_B, self.nvirt_B1059 1060 if monomer == 'B':1061 w_n = 2 * np.einsum('rbab->ar', self.v('rbab'))1062 w_n += self.V_B_AA[self.slices['a'], self.slices['r']]1063 eps_ov = (self.eps('a', dim=2) - self.eps('r'))1064 v_term1 = 'raar'1065 v_term2 = 'rara'1066 no, nv = self.ndocc_A, self.nvirt_A1067 1068 # Form A matrix (LHS)1069 voov = self.v(v_term1)1070 v_vOoV = 2 * voov - self.v(v_term2).swapaxes(2, 3)1071 v_ooaa = voov.swapaxes(1, 3)1072 v_vVoO = 2 * v_ooaa - v_ooaa.swapaxes(2, 3)1073 A_ovOV = np.einsum('vOoV->ovOV', v_vOoV + v_vVoO.swapaxes(1, 3))1074 1075 # Mangled the indices so badly with strides we need to copy back to C contiguous1076 nov = nv * no1077 A_ovOV = A_ovOV.reshape(nov, nov).copy(order='C')1078 A_ovOV[np.diag_indices_from(A_ovOV)] -= eps_ov.ravel()1079 1080 # Call DGESV, need flat ov array1081 B_ov = -1 * w_n.ravel()1082 t = np.linalg.solve(A_ovOV, B_ov)1083 # Our notation wants vo array1084 t = t.reshape(no, nv).T1085 1086 if ind:1087 # E200 Induction energy is free at the point1088 e20_ind = 2 * np.einsum('vo,ov->', t, w_n)1089 return (t, e20_ind)1090 else:1091 return t1092 1093 def compute_sapt_JK(self, Cleft, Cright, tensor=None):1094 1095 if self.alg != ""AO"":1096 raise Exception(""Attempted a call to JK builder in an MO algorithm"")1097 1098 if self.reference == ""ROHF"":1099 raise Exception(""AO algorithm not yet implemented for ROHF reference."")1100 1101 return_single = False1102 if not isinstance(Cleft, (list, tuple)):1103 Cleft = [Cleft]1104 return_single = True1105 if not isinstance(Cright, (list, tuple)):1106 Cright = [Cright]1107 return_single = True1108 if (not isinstance(tensor, (list, tuple))) and (tensor is not None):1109 tensor = [tensor]1110 return_single = True1111 1112 if len(Cleft) != len(Cright):1113 raise Exception(""Cleft list is not the same length as Cright list"")1114 1115 zero_append = []1116 num_compute = 01117 1118 for num in range(len(Cleft)):1119 Cl = Cleft[num]1120 Cr = Cright[num]1121 1122 if (Cr.shape[1] == 0) or (Cl.shape[1] == 0):1123 zero_append.append(num)1124 continue1125 1126 if tensor is not None:1127 mol = Cl.shape[1]1128 mor = Cr.shape[1]1129 1130 if (tensor[num].shape[0] != mol) or (tensor[num].shape[1] != mor):1131 raise Exception(""compute_sapt_JK: Tensor size does not match Cl (%d) /Cr (%d) : %s"" %1132 (mol, mor, str(tensor[num].shape)))1133 if mol < mor:1134 Cl = np.dot(Cl, tensor[num])1135 else:1136 Cr = np.dot(Cr, tensor[num].T)1137 1138 Cl = psi4.core.Matrix.from_array(Cl)1139 Cr = psi4.core.Matrix.from_array(Cr)1140 1141 self.jk.C_left_add(Cl)1142 self.jk.C_right_add(Cr)1143 num_compute += 11144 1145 self.jk.compute() 1146 1147 J_list = []1148 K_list = []1149 for num in range(num_compute):1150 J_list.append(np.array(self.jk.J()[num])) 1151 K_list.append(np.array(self.jk.K()[num])) 1152 1153 self.jk.C_clear()1154 1155 z = np.zeros((self.nmo, self.nmo))1156 for num in zero_append:1157 J_list.insert(num, z)1158 K_list.insert(num, z)1159 1160 if return_single:1161 return J_list[0], K_list[0]1162 else:1163 return J_list, K_list1164 1165 def chain_dot(self, *dot_list):1166 result = dot_list[0]1167 for x in range(len(dot_list) - 1):1168 result = np.dot(result, dot_list[x + 1])1169 return result1170 1171# End SAPT helper1172 1173class sapt_timer(object):1174 def __init__(self, name):1175 self.name = name1176 self.start = time.time()1177 print('\nStarting %s...' % name)1178 1179 def stop(self):1180 t = time.time() - self.start1181 print('...%s took a total of % .2f seconds.' % (self.name, t))1182 1183 1184def sapt_printer(line, value):1185 spacer = ' ' * (20 - len(line))1186 print(line + spacer + '% 16.8f mH % 16.8f kcal/mol' % (value * 1000, value * 627.509))1187# End SAPT helper1188","Python"
1189"Homo-Lumo","Mishima-syk/psikit","psikit/fsapt_helper.py",".py","1697","47","from rdkit import Chem1190from rdkit.Chem import AllChem1191from rdkit.Chem import Recap1192from rdkit.Chem import rdChemReactions1193import os1194 1195Recap.reactions += tuple([rdChemReactions.ReactionFromSmarts('[c:1]-[X:2]>>[c:1]*.*[X:2]'), rdChemReactions.ReactionFromSmarts('[c:1]-[OH1:2]>>[c:1]*.*[OH1:2]')])1196 1197def get_neighbor_h(atom_idx, mol):1198 atom = mol.GetAtomWithIdx(atom_idx)1199 neis = atom.GetNeighbors()1200 res = []