{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Populating the interactive namespace from numpy and matplotlib\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "WARNING: pylab import has clobbered these variables: ['nanmean']\n", "`%matplotlib` prevents importing * from pylab and numpy\n" ] } ], "source": [ "import os\n", "import re\n", "import sys\n", "import dataset\n", "from scipy.stats import nanmean\n", "from scipy.optimize import curve_fit\n", "import math\n", "#import numpy\n", "#import matplotlib.pyplot as plt\n", "%pylab inline\n" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def au2kcal(x):\n", " if type(x) is not list: \n", " x = [ x ]\n", " islist=0\n", " if x==[None]:\n", " return None\n", " else:\n", " islist=1\n", " for i in range(len(x)):\n", " x[i]=x[i]*627.50960803\n", " if not islist:\n", " x=x[0]\n", " return x" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def extapolateTau(taus,energies,errors):\n", " def func(x, a,b):\n", " if type(x) is not list: \n", " x = [ x ]\n", " islist=0\n", " else:\n", " islist=1\n", " res=[0]*len(x)\n", " for i in range(len(x)):\n", " res[i]=a*x[i]+b\n", " if islist==0:\n", " res=res[0]\n", " return res\n", " \n", " popt, pcov = curve_fit(func, taus, energies, sigma=errors, absolute_sigma=True, check_finite=True)\n", " perr = np.sqrt(np.diag(pcov))\n", " \n", "# plt.figure()\n", "# plt.errorbar(taus, au2kcal(energies), yerr=au2kcal(errors), fmt='ro')\n", "# t=arange(0.0, 0.01, 0.001)\n", "# s = au2kcal(func(t,popt[0],popt[1]))\n", "# plt.plot(t, s)\n", "# plt.autoscale(enable=True, axis=u'both', tight=True)\n", "# ax=plt.gca().get_ybound()\n", "# plt.ylim([ax[0],ax[0]+10])\n", "# plt.show()\n", " \n", " return[popt[1],perr[1]]" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false }, "outputs": [], "source": [ "moleculeDB='3dMLBE20_JTCT11-2036_2015.db' #/home/kdd/test/3dMLBE20_JTCT11-2036_2015.db'\n", "jobDB='jobDB.db'#/home/kdd/test/jobDB.db'\n", "\n" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": true }, "outputs": [], "source": [ "moleculeDB='3dMLBE20_JTCT11-2036_2015.db' #/home/kdd/test/3dMLBE20_JTCT11-2036_2015.db'\n", "jobDB='jobDB.db'#/home/kdd/test/jobDB.db'\n", "#get list of compounds\n", "with dataset.connect('sqlite:///'+moleculeDB) as DB3d:\n", " compoundList=[]\n", " compoundListMol=[]\n", " for table in ['atoms', 'dimers']:\n", " for molecule in DB3d[table]:\n", " if table=='atoms':\n", " compoundList.append(molecule['atom'])\n", " elif table=='dimers':\n", " compoundList.append(molecule['compound'])\n", " compoundListMol.append(molecule['compound'])" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false }, "outputs": [], "source": [ "dfttypList=['B3LYP']\n", "compoundListShort=['Ti', 'Cr', 'Mn', 'Cl', 'S', 'TiCl', 'CrCl', 'MnS']\n", "compoundListMolShort=['TiCl', 'CrCl', 'MnS']\n", "typList=['EMIN_casino_sJ','DMC_casino_sJ','EMIN_casino_lJ','DMC_casino_lJ', 'EMIN_qwalk_sJ','DMC_qwalk_sJ', 'EMIN_qwalk_lJ','DMC_qwalk_lJ']\n", "dtList=[0.002,0.005,0.008]\n", "energy=numpy.zeros((len(compoundListShort),len(dfttypList)*len(typList)))\n", "error=numpy.zeros((len(compoundListShort),len(dfttypList)*len(typList)))\n", "variance=numpy.zeros((len(compoundListShort),len(dfttypList)*len(typList)))\n", "dissocEnergy=numpy.zeros((len(compoundListMolShort),len(dfttypList)*len(typList)))\n", "dissocError=numpy.zeros((len(compoundListMolShort),len(dfttypList)*len(typList)))\n", "dissocEnergyExp=numpy.zeros(len(compoundListMolShort))\n", "energyExtrap=numpy.zeros((len(compoundListShort),3*len(dfttypList)*len(typList)))\n", "errorExtrap=numpy.zeros((len(compoundListShort),3*len(dfttypList)*len(typList)))\n", "dissocEnergyExtrap=numpy.zeros((len(compoundListMolShort),3*len(dfttypList)*len(typList)))\n", "dissocErrorExtrap=numpy.zeros((len(compoundListMolShort),3*len(dfttypList)*len(typList)))\n" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false }, "outputs": [], "source": [ "#get emin energy results casino small J\n", "typ=0\n", "for compound in range(len(compoundListShort)):\n", " for dfttyp in range(len(dfttypList)):\n", " jobRecord={}\n", " jobRecord['dfttyp']=dfttypList[dfttyp]\n", " jobRecord['scftyp']='UHF'\n", " jobRecord['basisXMLFiles']='[\"BFD_Library_noPseudoH.xml\", \"aug-cc-pVTZ-DK_Diffuse.xml\"]'\n", " jobRecord['basisNames']='[\"vtz\", \"aug-cc-pVTZ-DK_Diffuse\"]'\n", " jobRecord['executingProgram']='casino'\n", " jobRecord['Nu']=4\n", " jobRecord['Nchi']=4\n", " jobRecord['Nf']=2\n", " jobRecord['qmctyp']='vmc'\n", " jobRecord['compound']=compoundListShort[compound]\n", "\n", " with dataset.connect('sqlite:///'+jobDB) as DBdata:\n", " results=DBdata['data'].find(**jobRecord)\n", "\n", " for result in results:\n", " if result['jastrowFrom']==-1:\n", " continue\n", " elif '_true' in result['symmetry']:\n", " continue\n", " else:\n", " with dataset.connect('sqlite:///'+jobDB) as DBdata:\n", " jastrowFrom=DBdata['data'].find_one(id=result['jastrowFrom'])\n", " if jastrowFrom['qmctyp']!='emin':\n", " continue\n", " with dataset.connect('sqlite:///'+jobDB) as DBdata:\n", " jastrowFrom2=DBdata['data'].find_one(id=jastrowFrom['jastrowFrom'])\n", " if jastrowFrom2['qmctyp']!='varmin':\n", " continue \n", "\n", " if energy[compound,dfttyp+len(dfttypList)*typ]!=0:\n", " sys.exit(\"something wrong here: two energies found for one entry.\")\n", " if result['status']=='Done':\n", " energy[compound,dfttyp+len(dfttypList)*typ]=au2kcal(result['energy'])\n", " error[compound,dfttyp+len(dfttypList)*typ]=au2kcal(result['error'])\n", " variance[compound,dfttyp+len(dfttypList)*typ]=result['variance']\n", " else:\n", " energy[compound,dfttyp+len(dfttypList)*typ]=None\n", " error[compound,dfttyp+len(dfttypList)*typ]=0.0\n", " variance[compound,dfttyp+len(dfttypList)*typ]=None" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": true }, "outputs": [], "source": [ "#get emin energy results casino large J\n", "typ=2\n", "for compound in range(len(compoundListShort)):\n", " for dfttyp in range(len(dfttypList)):\n", " jobRecord={}\n", " jobRecord['dfttyp']=dfttypList[dfttyp]\n", " jobRecord['scftyp']='UHF'\n", " jobRecord['basisXMLFiles']='[\"BFD_Library_noPseudoH.xml\", \"aug-cc-pVTZ-DK_Diffuse.xml\"]'\n", " jobRecord['basisNames']='[\"vtz\", \"aug-cc-pVTZ-DK_Diffuse\"]'\n", " jobRecord['executingProgram']='casino'\n", " jobRecord['Nu']=6\n", " jobRecord['Nchi']=6\n", " jobRecord['Nf']=3\n", " jobRecord['qmctyp']='vmc'\n", " jobRecord['compound']=compoundListShort[compound]\n", "\n", " with dataset.connect('sqlite:///'+jobDB) as DBdata:\n", " results=DBdata['data'].find(**jobRecord)\n", "\n", " for result in results:\n", " \n", " if result['jastrowFrom']==-1:\n", " continue\n", " elif '_true' in result['symmetry']:\n", " continue\n", " else:\n", " with dataset.connect('sqlite:///'+jobDB) as DBdata:\n", " jastrowFrom=DBdata['data'].find_one(id=result['jastrowFrom'])\n", " if jastrowFrom['qmctyp']!='emin':\n", " continue\n", " with dataset.connect('sqlite:///'+jobDB) as DBdata:\n", " jastrowFrom2=DBdata['data'].find_one(id=jastrowFrom['jastrowFrom'])\n", " if jastrowFrom2['qmctyp']!='varmin':\n", " continue \n", "\n", " if energy[compound,dfttyp+len(dfttypList)*typ]!=0:\n", " sys.exit(\"something wrong here: two energies found for one entry.\")\n", " if result['status']=='Done':\n", " energy[compound,dfttyp+len(dfttypList)*typ]=au2kcal(result['energy'])\n", " error[compound,dfttyp+len(dfttypList)*typ]=au2kcal(result['error'])\n", " variance[compound,dfttyp+len(dfttypList)*typ]=result['variance']\n", " else:\n", " energy[compound,dfttyp+len(dfttypList)*typ]=None\n", " error[compound,dfttyp+len(dfttypList)*typ]=0.0\n", " variance[compound,dfttyp+len(dfttypList)*typ]=None" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": false }, "outputs": [], "source": [ "#get emin energy results qwalk small J\n", "typ=4\n", "for compound in range(len(compoundListShort)):\n", " for dfttyp in range(len(dfttypList)):\n", " jobRecord={}\n", " jobRecord['dfttyp']=dfttypList[dfttyp]\n", " jobRecord['scftyp']='UHF'\n", " jobRecord['basisXMLFiles']='[\"BFD_Library_noPseudoH.xml\", \"aug-cc-pVTZ-DK_Diffuse.xml\"]'\n", " jobRecord['basisNames']='[\"vtz\", \"aug-cc-pVTZ-DK_Diffuse\"]'\n", " jobRecord['executingProgram']='qwalk'\n", " jobRecord['qmctyp']='vmc'\n", " jobRecord['Nu']=3\n", " jobRecord['Nchi']=4\n", " jobRecord['Nf']=4\n", " jobRecord['compound']=compoundListShort[compound]\n", "\n", " with dataset.connect('sqlite:///'+jobDB) as DBdata:\n", " results=DBdata['data'].find(**jobRecord)\n", "\n", " for result in results:\n", " if result['jastrowFrom']==-1:\n", " continue\n", " elif '_true' in result['symmetry']:\n", " continue\n", " else:\n", " with dataset.connect('sqlite:///'+jobDB) as DBdata:\n", " jastrowFrom=DBdata['data'].find_one(id=result['jastrowFrom'])\n", " if jastrowFrom['qmctyp']!='emin':\n", " continue\n", " with dataset.connect('sqlite:///'+jobDB) as DBdata:\n", " jastrowFrom2=DBdata['data'].find_one(id=jastrowFrom['jastrowFrom'])\n", " if jastrowFrom2['qmctyp']!='varmin':\n", " continue\n", " \n", "\n", " if energy[compound,dfttyp+len(dfttypList)*typ]!=0:\n", " sys.exit(\"something wrong here: two energies found for one entry.\")\n", " if result['status'] in ['Done','LargeErrorbar']:\n", " energy[compound,dfttyp+len(dfttypList)*typ]=au2kcal(result['energy'])\n", " error[compound,dfttyp+len(dfttypList)*typ]=au2kcal(result['error'])\n", " variance[compound,dfttyp+len(dfttypList)*typ]=result['variance']\n", " else:\n", " energy[compound,dfttyp+len(dfttypList)*typ]=None\n", " error[compound,dfttyp+len(dfttypList)*typ]=0.0\n", " variance[compound,dfttyp+len(dfttypList)*typ]=None" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": true }, "outputs": [], "source": [ "#get emin energy results qwalk large J\n", "typ=6\n", "for compound in range(len(compoundListShort)):\n", " for dfttyp in range(len(dfttypList)):\n", " jobRecord={}\n", " jobRecord['dfttyp']=dfttypList[dfttyp]\n", " jobRecord['scftyp']='UHF'\n", " jobRecord['basisXMLFiles']='[\"BFD_Library_noPseudoH.xml\", \"aug-cc-pVTZ-DK_Diffuse.xml\"]'\n", " jobRecord['basisNames']='[\"vtz\", \"aug-cc-pVTZ-DK_Diffuse\"]'\n", " jobRecord['executingProgram']='qwalk'\n", " jobRecord['qmctyp']='vmc'\n", " jobRecord['Nu']=3\n", " jobRecord['Nchi']=4\n", " jobRecord['Nf']=5\n", " jobRecord['compound']=compoundListShort[compound]\n", "\n", " with dataset.connect('sqlite:///'+jobDB) as DBdata:\n", " results=DBdata['data'].find(**jobRecord)\n", "\n", " for result in results:\n", " if result['jastrowFrom']==-1:\n", " continue\n", " elif '_true' in result['symmetry']:\n", " continue\n", " else:\n", " with dataset.connect('sqlite:///'+jobDB) as DBdata:\n", " jastrowFrom=DBdata['data'].find_one(id=result['jastrowFrom'])\n", " if jastrowFrom['qmctyp']!='emin':\n", " continue\n", " with dataset.connect('sqlite:///'+jobDB) as DBdata:\n", " jastrowFrom2=DBdata['data'].find_one(id=jastrowFrom['jastrowFrom'])\n", " if jastrowFrom2['qmctyp']!='varmin':\n", " continue\n", " \n", "\n", " if energy[compound,dfttyp+len(dfttypList)*typ]!=0:\n", " sys.exit(\"something wrong here: two energies found for one entry.\")\n", " if result['status']=='Done':\n", " energy[compound,dfttyp+len(dfttypList)*typ]=au2kcal(result['energy'])\n", " error[compound,dfttyp+len(dfttypList)*typ]=au2kcal(result['error'])\n", " variance[compound,dfttyp+len(dfttypList)*typ]=result['variance']\n", " else:\n", " energy[compound,dfttyp+len(dfttypList)*typ]=None\n", " error[compound,dfttyp+len(dfttypList)*typ]=0.0\n", " variance[compound,dfttyp+len(dfttypList)*typ]=None" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "collapsed": false }, "outputs": [], "source": [ "#get dmc energy casino small J extrapolated to tau=0 and at tau=[0.002,0.005,0.008]\n", "typ=1\n", "for compound in range(len(compoundListShort)):\n", " for dfttyp in range(len(dfttypList)):\n", " jobRecord={}\n", " jobRecord['dfttyp']=dfttypList[dfttyp]\n", " jobRecord['scftyp']='UHF'\n", " jobRecord['basisXMLFiles']='[\"BFD_Library_noPseudoH.xml\", \"aug-cc-pVTZ-DK_Diffuse.xml\"]'\n", " jobRecord['basisNames']='[\"vtz\", \"aug-cc-pVTZ-DK_Diffuse\"]'\n", " jobRecord['executingProgram']='casino'\n", " jobRecord['qmctyp']='vmc_dmc'\n", " jobRecord['compound']=compoundListShort[compound]\n", " jobRecord['Nu']=4\n", " jobRecord['Nchi']=4\n", " jobRecord['Nf']=2\n", "\n", " with dataset.connect('sqlite:///'+jobDB) as DBdata:\n", " results=DBdata['data'].find(**jobRecord)\n", "\n", " resultsExtrap=[]\n", " for result in results:\n", " if result['jastrowFrom']==-1:\n", " continue\n", " elif '_true' in result['symmetry']:\n", " continue\n", " else:\n", " with dataset.connect('sqlite:///'+jobDB) as DBdata:\n", " jastrowFrom=DBdata['data'].find_one(id=result['jastrowFrom'])\n", " if jastrowFrom['qmctyp']!='emin':\n", " continue\n", " with dataset.connect('sqlite:///'+jobDB) as DBdata:\n", " jastrowFrom2=DBdata['data'].find_one(id=jastrowFrom['jastrowFrom'])\n", " if jastrowFrom2['qmctyp']!='varmin':\n", " continue\n", " resultsExtrap.append(result)\n", " \n", " for i in reversed(range(len(resultsExtrap))):\n", " if not resultsExtrap[i]['status'] in [\"Done\",\"LargeErrorbar2\", \"BadReblock\"]:\n", " resultsExtrap.pop(i)\n", " \n", " if len(resultsExtrap)<2:\n", " energy[compound,dfttyp+len(dfttypList)*typ]=None\n", " error[compound,dfttyp+len(dfttypList)*typ]=0.0\n", " else: #extrapolate to zero Time\n", "# print([resultsExtrap[i]['dt'] for i in range(len(resultsExtrap))])\n", "# print([resultsExtrap[i]['energy'] for i in range(len(resultsExtrap))])\n", "# print([resultsExtrap[i]['error'] for i in range(len(resultsExtrap))])\n", " [energy[compound,dfttyp+len(dfttypList)*typ], error[compound,dfttyp+len(dfttypList)*typ]]=\\\n", " au2kcal(extapolateTau([resultsExtrap[i]['dt'] for i in range(len(resultsExtrap))],\n", " [resultsExtrap[i]['energy'] for i in range(len(resultsExtrap))],\n", " [resultsExtrap[i]['error'] for i in range(len(resultsExtrap))]))\n", " #non extrapolated results\n", " energyExtrap[compound,len(dtList)*(dfttyp+len(dfttypList)*typ):len(dtList)*(dfttyp+len(dfttypList)*typ)+len(dtList)]=None\n", " errorExtrap[compound,len(dtList)*(dfttyp+len(dfttypList)*typ):len(dtList)*(dfttyp+len(dfttypList)*typ)+len(dtList)]=None\n", " for i in range(len(resultsExtrap)):\n", " if resultsExtrap[i]['dt']==dtList[0]:\n", " dt=0\n", " elif resultsExtrap[i]['dt']==dtList[1]:\n", " dt=1\n", " elif resultsExtrap[i]['dt']==dtList[2]:\n", " dt=2\n", " energyExtrap[compound,len(dtList)*(dfttyp+len(dfttypList)*typ)+dt]=au2kcal(resultsExtrap[i]['energy'])\n", " errorExtrap[compound,len(dtList)*(dfttyp+len(dfttypList)*typ)+dt]=au2kcal(resultsExtrap[i]['error'])\n", " " ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "collapsed": true }, "outputs": [], "source": [ "#get dmc energy casino large J extrapolated to tau=0 and at tau=[0.002,0.005,0.008]\n", "typ=3\n", "for compound in range(len(compoundListShort)):\n", " for dfttyp in range(len(dfttypList)):\n", " jobRecord={}\n", " jobRecord['dfttyp']=dfttypList[dfttyp]\n", " jobRecord['scftyp']='UHF'\n", " jobRecord['basisXMLFiles']='[\"BFD_Library_noPseudoH.xml\", \"aug-cc-pVTZ-DK_Diffuse.xml\"]'\n", " jobRecord['basisNames']='[\"vtz\", \"aug-cc-pVTZ-DK_Diffuse\"]'\n", " jobRecord['executingProgram']='casino'\n", " jobRecord['qmctyp']='vmc_dmc'\n", " jobRecord['compound']=compoundListShort[compound]\n", " jobRecord['Nu']=6\n", " jobRecord['Nchi']=6\n", " jobRecord['Nf']=3\n", "\n", " with dataset.connect('sqlite:///'+jobDB) as DBdata:\n", " results=DBdata['data'].find(**jobRecord)\n", "\n", " resultsExtrap=[]\n", " for result in results:\n", " if result['jastrowFrom']==-1:\n", " continue\n", " elif '_true' in result['symmetry']:\n", " continue\n", " else:\n", " with dataset.connect('sqlite:///'+jobDB) as DBdata:\n", " jastrowFrom=DBdata['data'].find_one(id=result['jastrowFrom'])\n", " if jastrowFrom['qmctyp']!='emin':\n", " continue\n", " with dataset.connect('sqlite:///'+jobDB) as DBdata:\n", " jastrowFrom2=DBdata['data'].find_one(id=jastrowFrom['jastrowFrom'])\n", " if jastrowFrom2['qmctyp']!='varmin':\n", " continue\n", " resultsExtrap.append(result)\n", " \n", " for i in reversed(range(len(resultsExtrap))):\n", " if not resultsExtrap[i]['status'] in [\"Done\",\"LargeErrorbar2\", \"BadReblock\"]:\n", " resultsExtrap.pop(i)\n", " \n", " if len(resultsExtrap)<2:\n", " energy[compound,dfttyp+len(dfttypList)*typ]=None\n", " error[compound,dfttyp+len(dfttypList)*typ]=0.0\n", " else: #extrapolate to zero Time\n", "# print([resultsExtrap[i]['dt'] for i in range(len(resultsExtrap))])\n", "# print([resultsExtrap[i]['energy'] for i in range(len(resultsExtrap))])\n", "# print([resultsExtrap[i]['error'] for i in range(len(resultsExtrap))])\n", " [energy[compound,dfttyp+len(dfttypList)*typ], error[compound,dfttyp+len(dfttypList)*typ]]=\\\n", " au2kcal(extapolateTau([resultsExtrap[i]['dt'] for i in range(len(resultsExtrap))],\n", " [resultsExtrap[i]['energy'] for i in range(len(resultsExtrap))],\n", " [resultsExtrap[i]['error'] for i in range(len(resultsExtrap))]))\n", " #non extrapolated results\n", " energyExtrap[compound,len(dtList)*(dfttyp+len(dfttypList)*typ):len(dtList)*(dfttyp+len(dfttypList)*typ)+len(dtList)]=None\n", " errorExtrap[compound,len(dtList)*(dfttyp+len(dfttypList)*typ):len(dtList)*(dfttyp+len(dfttypList)*typ)+len(dtList)]=None\n", " for i in range(len(resultsExtrap)):\n", " if resultsExtrap[i]['dt']==dtList[0]:\n", " dt=0\n", " elif resultsExtrap[i]['dt']==dtList[1]:\n", " dt=1\n", " elif resultsExtrap[i]['dt']==dtList[2]:\n", " dt=2\n", " energyExtrap[compound,len(dtList)*(dfttyp+len(dfttypList)*typ)+dt]=au2kcal(resultsExtrap[i]['energy'])\n", " errorExtrap[compound,len(dtList)*(dfttyp+len(dfttypList)*typ)+dt]=au2kcal(resultsExtrap[i]['error'])\n", " " ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "collapsed": false }, "outputs": [], "source": [ "#get dmc energy qwalk small J extrapolated to tau=0 and at tau=[0.002,0.005,0.008]\n", "typ=5\n", "for compound in range(len(compoundListShort)):\n", " for dfttyp in range(len(dfttypList)):\n", " jobRecord={}\n", " jobRecord['dfttyp']=dfttypList[dfttyp]\n", " jobRecord['scftyp']='UHF'\n", " jobRecord['basisXMLFiles']='[\"BFD_Library_noPseudoH.xml\", \"aug-cc-pVTZ-DK_Diffuse.xml\"]'\n", " jobRecord['basisNames']='[\"vtz\", \"aug-cc-pVTZ-DK_Diffuse\"]'\n", " jobRecord['executingProgram']='qwalk'\n", " jobRecord['qmctyp']='vmc_dmc'\n", " jobRecord['compound']=compoundListShort[compound]\n", " jobRecord['Nu']=3\n", " jobRecord['Nchi']=4\n", " jobRecord['Nf']=4\n", "\n", " with dataset.connect('sqlite:///'+jobDB) as DBdata:\n", " results=DBdata['data'].find(**jobRecord)\n", "\n", " resultsExtrap=[]\n", " for result in results:\n", " if result['jastrowFrom']==-1:\n", " continue\n", " elif '_true' in result['symmetry']:\n", " continue\n", " else:\n", " with dataset.connect('sqlite:///'+jobDB) as DBdata:\n", " jastrowFrom=DBdata['data'].find_one(id=result['jastrowFrom'])\n", " if jastrowFrom['qmctyp']!='emin':\n", " continue\n", " with dataset.connect('sqlite:///'+jobDB) as DBdata:\n", " jastrowFrom2=DBdata['data'].find_one(id=jastrowFrom['jastrowFrom'])\n", " if jastrowFrom2['qmctyp']!='varmin':\n", " continue\n", " resultsExtrap.append(result)\n", " \n", " for i in reversed(range(len(resultsExtrap))):\n", " if not resultsExtrap[i]['status'] in [\"Done\",\"LargeErrorbar2\", \"BadReblock\"]:\n", " resultsExtrap.pop(i)\n", " \n", " if len(resultsExtrap)<2:\n", " energy[compound,dfttyp+len(dfttypList)*typ]=None\n", " error[compound,dfttyp+len(dfttypList)*typ]=0.0\n", " else: #extrapolate to zero Time\n", "# print([resultsExtrap[i]['dt'] for i in range(len(resultsExtrap))])\n", "# print([resultsExtrap[i]['energy'] for i in range(len(resultsExtrap))])\n", "# print([resultsExtrap[i]['error'] for i in range(len(resultsExtrap))])\n", " [energy[compound,dfttyp+len(dfttypList)*typ], error[compound,dfttyp+len(dfttypList)*typ]]=\\\n", " au2kcal(extapolateTau([resultsExtrap[i]['dt'] for i in range(len(resultsExtrap))],\n", " [resultsExtrap[i]['energy'] for i in range(len(resultsExtrap))],\n", " [resultsExtrap[i]['error'] for i in range(len(resultsExtrap))]))\n", " #non extrapolated results\n", " energyExtrap[compound,len(dtList)*(dfttyp+len(dfttypList)*typ):len(dtList)*(dfttyp+len(dfttypList)*typ)+len(dtList)]=None\n", " errorExtrap[compound,len(dtList)*(dfttyp+len(dfttypList)*typ):len(dtList)*(dfttyp+len(dfttypList)*typ)+len(dtList)]=None\n", " for i in range(len(resultsExtrap)):\n", " if resultsExtrap[i]['dt']==dtList[0]:\n", " dt=0\n", " elif resultsExtrap[i]['dt']==dtList[1]:\n", " dt=1\n", " elif resultsExtrap[i]['dt']==dtList[2]:\n", " dt=2\n", " energyExtrap[compound,len(dtList)*(dfttyp+len(dfttypList)*typ)+dt]=au2kcal(resultsExtrap[i]['energy'])\n", " errorExtrap[compound,len(dtList)*(dfttyp+len(dfttypList)*typ)+dt]=au2kcal(resultsExtrap[i]['error'])" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "collapsed": true }, "outputs": [], "source": [ "#get dmc energy qwalk large J extrapolated to tau=0 and at tau=[0.002,0.005,0.008]\n", "typ=7\n", "for compound in range(len(compoundListShort)):\n", " for dfttyp in range(len(dfttypList)):\n", " jobRecord={}\n", " jobRecord['dfttyp']=dfttypList[dfttyp]\n", " jobRecord['scftyp']='UHF'\n", " jobRecord['basisXMLFiles']='[\"BFD_Library_noPseudoH.xml\", \"aug-cc-pVTZ-DK_Diffuse.xml\"]'\n", " jobRecord['basisNames']='[\"vtz\", \"aug-cc-pVTZ-DK_Diffuse\"]'\n", " jobRecord['executingProgram']='qwalk'\n", " jobRecord['qmctyp']='vmc_dmc'\n", " jobRecord['compound']=compoundListShort[compound]\n", " jobRecord['Nu']=3\n", " jobRecord['Nchi']=4\n", " jobRecord['Nf']=5\n", "\n", " with dataset.connect('sqlite:///'+jobDB) as DBdata:\n", " results=DBdata['data'].find(**jobRecord)\n", "\n", " resultsExtrap=[]\n", " for result in results:\n", " if result['jastrowFrom']==-1:\n", " continue\n", " elif '_true' in result['symmetry']:\n", " continue\n", " else:\n", " with dataset.connect('sqlite:///'+jobDB) as DBdata:\n", " jastrowFrom=DBdata['data'].find_one(id=result['jastrowFrom'])\n", " if jastrowFrom['qmctyp']!='emin':\n", " continue\n", " with dataset.connect('sqlite:///'+jobDB) as DBdata:\n", " jastrowFrom2=DBdata['data'].find_one(id=jastrowFrom['jastrowFrom'])\n", " if jastrowFrom2['qmctyp']!='varmin':\n", " continue\n", " resultsExtrap.append(result)\n", " \n", " for i in reversed(range(len(resultsExtrap))):\n", " if not resultsExtrap[i]['status'] in [\"Done\",\"LargeErrorbar2\", \"BadReblock\"]:\n", " resultsExtrap.pop(i)\n", " \n", " if len(resultsExtrap)<2:\n", " energy[compound,dfttyp+len(dfttypList)*typ]=None\n", " error[compound,dfttyp+len(dfttypList)*typ]=0.0\n", " else: #extrapolate to zero Time\n", "# print([resultsExtrap[i]['dt'] for i in range(len(resultsExtrap))])\n", "# print([resultsExtrap[i]['energy'] for i in range(len(resultsExtrap))])\n", "# print([resultsExtrap[i]['error'] for i in range(len(resultsExtrap))])\n", " [energy[compound,dfttyp+len(dfttypList)*typ], error[compound,dfttyp+len(dfttypList)*typ]]=\\\n", " au2kcal(extapolateTau([resultsExtrap[i]['dt'] for i in range(len(resultsExtrap))],\n", " [resultsExtrap[i]['energy'] for i in range(len(resultsExtrap))],\n", " [resultsExtrap[i]['error'] for i in range(len(resultsExtrap))]))\n", " #non extrapolated results\n", " energyExtrap[compound,len(dtList)*(dfttyp+len(dfttypList)*typ):len(dtList)*(dfttyp+len(dfttypList)*typ)+len(dtList)]=None\n", " errorExtrap[compound,len(dtList)*(dfttyp+len(dfttypList)*typ):len(dtList)*(dfttyp+len(dfttypList)*typ)+len(dtList)]=None\n", " for i in range(len(resultsExtrap)):\n", " if resultsExtrap[i]['dt']==dtList[0]:\n", " dt=0\n", " elif resultsExtrap[i]['dt']==dtList[1]:\n", " dt=1\n", " elif resultsExtrap[i]['dt']==dtList[2]:\n", " dt=2\n", " energyExtrap[compound,len(dtList)*(dfttyp+len(dfttypList)*typ)+dt]=au2kcal(resultsExtrap[i]['energy'])\n", " errorExtrap[compound,len(dtList)*(dfttyp+len(dfttypList)*typ)+dt]=au2kcal(resultsExtrap[i]['error'])" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\t cas sJ \t cas lJ \t qwalk sJ \t qwalk lJ\n", "Ti\n", "VMC: \t -36486.5(2) \t -36491.3(2) \t -36490.8(2) \t -36491.2(2) \n", "DMC: \t -36516.2(1) \t -36517.1(2) \t -36516.2(2) \t -36515.5(2) \n", "\n", "Cr\n", "VMC: \t -54504.6(2) \t -54511.7(3) \t -54496.7(3) \t -54500.4(2) \n", "DMC: \t -54529.5(2) \t -54530.0(2) \t -54527.9(3) \t -54528.2(2) \n", "\n", "Mn\n", "VMC: \t -65416.6(2) \t -65424.2(3) \t -65413.2(3) \t -65409.2(2) \n", "DMC: \t -65449.4(2) \t -65451.0(2) \t -65447.4(3) \t -65447.6(2) \n", "\n", "Cl\n", "VMC: \t -9386.6(1) \t -9390.3(1) \t -9386.9(1) \t -9387.6(1) \n", "DMC: \t -9394.6(1) \t -9395.1(2) \t -9394.2(1) \t -9394.2(1) \n", "\n", "S\n", "VMC: \t -6351.1(1) \t -6353.3(1) \t -6350.4(1) \t -6350.9(1) \n", "DMC: \t -6357.2(1) \t -6357.8(1) \t -6356.7(1) \t -6356.8(1) \n", "\n", "TiCl\n", "VMC: \t -45961.9(2) \t -45980.3(3) \t -45971.4(3) \t -45974.0(2) \n", "DMC: \t -46007.9(2) \t -46011.8(2) \t -46008.6(2) \t -46009.7(2) \n", "\n", "CrCl\n", "VMC: \t -63960.1(2) \t -63979.7(3) \t -63964.6(3) \t -63966.6(2) \n", "DMC: \t -64011.5(2) \t -64014.8(2) \t -64011.9(3) \t -64011.3(2) \n", "\n", "MnS\n", "VMC: \t -71812.7(2) \t -71811.9(3) \t -71808.3(3) \t -71780.1(2) \n", "DMC: \t -71867.0(2) \t -71868.4(2) \t -71867.5(3) \t -71865.3(2) \n", "\n" ] } ], "source": [ "#total Energies\n", "print(\"\\t cas sJ \\t cas lJ \\t qwalk sJ \\t qwalk lJ\")\n", "for comp in range(len(compoundListShort)):\n", " print(compoundListShort[comp])\n", " print(\"VMC: \\t{:10.1f}({:1.0f}) \\t{:10.1f}({:1.0f}) \\t{:10.1f}({:1.0f}) \\t{:10.1f}({:1.0f}) \".format(energy[comp,0], error[comp,0]*10, energy[comp,2], error[comp,4]*10, energy[comp,4], error[comp,4]*10, energy[comp,6], error[comp,6]*10))\n", " print(\"DMC: \\t{:10.1f}({:1.0f}) \\t{:10.1f}({:1.0f}) \\t{:10.1f}({:1.0f}) \\t{:10.1f}({:1.0f}) \".format(energy[comp,1], error[comp,1]*10, energy[comp,3], error[comp,3]*10, energy[comp,5], error[comp,5]*10, energy[comp,7], error[comp,7]*10))\n", " print(\"\")" ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\t cas sJ \t cas lJ \t qwalk sJ \t qwalk lJ\n", "Ti\n", "VMC: \t 0.0(2) \t -4.8(2) \t -4.3(2) \t -4.7(2) \n", "DMC: \t 0.0(1) \t -1.0(2) \t 0.0(2) \t 0.6(2) \n", "\n", "Cr\n", "VMC: \t 0.0(2) \t -7.1(3) \t 7.9(3) \t 4.1(2) \n", "DMC: \t 0.0(2) \t -0.4(2) \t 1.6(3) \t 1.4(2) \n", "\n", "Mn\n", "VMC: \t 0.0(2) \t -7.6(3) \t 3.4(3) \t 7.4(2) \n", "DMC: \t 0.0(2) \t -1.6(2) \t 2.0(3) \t 1.7(2) \n", "\n", "Cl\n", "VMC: \t 0.0(1) \t -3.7(1) \t -0.2(1) \t -1.0(1) \n", "DMC: \t 0.0(1) \t -0.4(2) \t 0.4(1) \t 0.4(1) \n", "\n", "S\n", "VMC: \t 0.0(1) \t -2.3(1) \t 0.7(1) \t 0.2(1) \n", "DMC: \t 0.0(1) \t -0.6(1) \t 0.4(1) \t 0.4(1) \n", "\n", "TiCl\n", "VMC: \t 0.0(2) \t -18.3(3) \t -9.4(3) \t -12.1(2) \n", "DMC: \t 0.0(2) \t -4.0(2) \t -0.8(2) \t -1.8(2) \n", "\n", "CrCl\n", "VMC: \t 0.0(2) \t -19.6(3) \t -4.5(3) \t -6.5(2) \n", "DMC: \t 0.0(2) \t -3.3(2) \t -0.4(3) \t 0.2(2) \n", "\n", "MnS\n", "VMC: \t 0.0(2) \t 0.8(3) \t 4.4(3) \t 32.6(2) \n", "DMC: \t 0.0(2) \t -1.4(2) \t -0.5(3) \t 1.8(2) \n", "\n" ] } ], "source": [ "#total Energy differences to casino small Jast\n", "print(\"\\t cas sJ \\t cas lJ \\t qwalk sJ \\t qwalk lJ\")\n", "for comp in range(len(compoundListShort)):\n", " print(compoundListShort[comp])\n", " print(\"VMC: \\t{:10.1f}({:1.0f}) \\t{:10.1f}({:1.0f}) \\t{:10.1f}({:1.0f}) \\t{:10.1f}({:1.0f}) \".format(energy[comp,0]-energy[comp,0], error[comp,0]*10,energy[comp,2]-energy[comp,0], error[comp,4]*10, energy[comp,4]-energy[comp,0], error[comp,4]*10, energy[comp,6]-energy[comp,0], error[comp,6]*10))\n", " print(\"DMC: \\t{:10.1f}({:1.0f}) \\t{:10.1f}({:1.0f}) \\t{:10.1f}({:1.0f}) \\t{:10.1f}({:1.0f}) \".format(energy[comp,1]-energy[comp,1], error[comp,1]*10, energy[comp,3]-energy[comp,1], error[comp,3]*10, energy[comp,5]-energy[comp,1], error[comp,5]*10, energy[comp,7]-energy[comp,1], error[comp,7]*10))\n", " print(\"\")" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\t cas sJ \t cas lJ \t qwalk sJ \t qwalk lJ\n", "Ti\n", "VMC: \t 0.0(2) \t 0.0(2) \t -4.3(2) \t 0.1(2) \n", "DMC: \t 0.0(1) \t 0.0(2) \t 0.0(2) \t 1.6(2) \n", "\n", "Cr\n", "VMC: \t 0.0(2) \t 0.0(3) \t 7.9(3) \t 11.2(2) \n", "DMC: \t 0.0(2) \t 0.0(2) \t 1.6(3) \t 1.8(2) \n", "\n", "Mn\n", "VMC: \t 0.0(2) \t 0.0(3) \t 3.4(3) \t 15.0(2) \n", "DMC: \t 0.0(2) \t 0.0(2) \t 2.0(3) \t 3.4(2) \n", "\n", "Cl\n", "VMC: \t 0.0(1) \t 0.0(1) \t -0.2(1) \t 2.7(1) \n", "DMC: \t 0.0(1) \t 0.0(2) \t 0.4(1) \t 0.8(1) \n", "\n", "S\n", "VMC: \t 0.0(1) \t 0.0(1) \t 0.7(1) \t 2.5(1) \n", "DMC: \t 0.0(1) \t 0.0(1) \t 0.4(1) \t 1.0(1) \n", "\n", "TiCl\n", "VMC: \t 0.0(2) \t 0.0(3) \t -9.4(3) \t 6.3(2) \n", "DMC: \t 0.0(2) \t 0.0(2) \t -0.8(2) \t 2.2(2) \n", "\n", "CrCl\n", "VMC: \t 0.0(2) \t 0.0(3) \t -4.5(3) \t 13.1(2) \n", "DMC: \t 0.0(2) \t 0.0(2) \t -0.4(3) \t 3.5(2) \n", "\n", "MnS\n", "VMC: \t 0.0(2) \t 0.0(3) \t 4.4(3) \t 31.7(2) \n", "DMC: \t 0.0(2) \t 0.0(2) \t -0.5(3) \t 3.1(2) \n", "\n" ] } ], "source": [ "#total Energy differences casino vs. qwalk for small and large Jast\n", "print(\"\\t cas sJ \\t cas lJ \\t qwalk sJ \\t qwalk lJ\")\n", "for comp in range(len(compoundListShort)):\n", " print(compoundListShort[comp])\n", " print(\"VMC: \\t{:10.1f}({:1.0f}) \\t{:10.1f}({:1.0f}) \\t{:10.1f}({:1.0f}) \\t{:10.1f}({:1.0f}) \".format(energy[comp,0]-energy[comp,0], error[comp,0]*10,energy[comp,2]-energy[comp,2], error[comp,4]*10, energy[comp,4]-energy[comp,0], error[comp,4]*10, energy[comp,6]-energy[comp,2], error[comp,6]*10))\n", " print(\"DMC: \\t{:10.1f}({:1.0f}) \\t{:10.1f}({:1.0f}) \\t{:10.1f}({:1.0f}) \\t{:10.1f}({:1.0f}) \".format(energy[comp,1]-energy[comp,1], error[comp,1]*10, energy[comp,3]-energy[comp,3], error[comp,3]*10, energy[comp,5]-energy[comp,1], error[comp,5]*10, energy[comp,7]-energy[comp,3], error[comp,7]*10))\n", " print(\"\")" ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\t cas sJ \t cas lJ \t qwalk sJ \t qwalk lJ\n", "Ti\n", "VMC: \t 0.0(2) \t -4.8(2) \t 0.0(2) \t -0.4(2) \n", "DMC: \t 0.0(1) \t -1.0(2) \t 0.0(2) \t 0.6(2) \n", "\n", "Cr\n", "VMC: \t 0.0(2) \t -7.1(3) \t 0.0(3) \t -3.8(2) \n", "DMC: \t 0.0(2) \t -0.4(2) \t 0.0(3) \t -0.2(2) \n", "\n", "Mn\n", "VMC: \t 0.0(2) \t -7.6(3) \t 0.0(3) \t 4.0(2) \n", "DMC: \t 0.0(2) \t -1.6(2) \t 0.0(3) \t -0.2(2) \n", "\n", "Cl\n", "VMC: \t 0.0(1) \t -3.7(1) \t 0.0(1) \t -0.7(1) \n", "DMC: \t 0.0(1) \t -0.4(2) \t 0.0(1) \t 0.0(1) \n", "\n", "S\n", "VMC: \t 0.0(1) \t -2.3(1) \t 0.0(1) \t -0.5(1) \n", "DMC: \t 0.0(1) \t -0.6(1) \t 0.0(1) \t -0.1(1) \n", "\n", "TiCl\n", "VMC: \t 0.0(2) \t -18.3(3) \t 0.0(3) \t -2.7(2) \n", "DMC: \t 0.0(2) \t -4.0(2) \t 0.0(2) \t -1.1(2) \n", "\n", "CrCl\n", "VMC: \t 0.0(2) \t -19.6(3) \t 0.0(3) \t -2.0(2) \n", "DMC: \t 0.0(2) \t -3.3(2) \t 0.0(3) \t 0.7(2) \n", "\n", "MnS\n", "VMC: \t 0.0(2) \t 0.8(3) \t 0.0(3) \t 28.2(2) \n", "DMC: \t 0.0(2) \t -1.4(2) \t 0.0(3) \t 2.2(2) \n", "\n" ] } ], "source": [ "#total Energy differences casino and qwalk seperately\n", "print(\"\\t cas sJ \\t cas lJ \\t qwalk sJ \\t qwalk lJ\")\n", "for comp in range(len(compoundListShort)):\n", " print(compoundListShort[comp])\n", " print(\"VMC: \\t{:10.1f}({:1.0f}) \\t{:10.1f}({:1.0f}) \\t{:10.1f}({:1.0f}) \\t{:10.1f}({:1.0f}) \".format(energy[comp,0]-energy[comp,0], error[comp,0]*10,energy[comp,2]-energy[comp,0], error[comp,4]*10, energy[comp,4]-energy[comp,4], error[comp,4]*10, energy[comp,6]-energy[comp,4], error[comp,6]*10))\n", " print(\"DMC: \\t{:10.1f}({:1.0f}) \\t{:10.1f}({:1.0f}) \\t{:10.1f}({:1.0f}) \\t{:10.1f}({:1.0f}) \".format(energy[comp,1]-energy[comp,1], error[comp,1]*10, energy[comp,3]-energy[comp,1], error[comp,3]*10, energy[comp,5]-energy[comp,5], error[comp,5]*10, energy[comp,7]-energy[comp,5], error[comp,7]*10))\n", " print(\"\")" ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\t cas sJ \t cas lJ \t qwalk sJ \t qwalk lJ\n", "TiCl\n", "VMC: \t 0.0(2) \t 9.9(2) \t 4.9(3) \t 6.4(2) \n", "DMC: \t 0.0(3) \t 2.6(3) \t 1.2(3) \t 2.9(3) \n", "\n", "CrCl\n", "VMC: \t 0.0(2) \t 8.9(2) \t 12.1(4) \t 9.7(2) \n", "DMC: \t 0.0(3) \t 2.4(4) \t 2.4(4) \t 1.6(3) \n", "\n", "MnS\n", "VMC: \t 0.0(2) \t -10.7(2) \t -0.3(5) \t -24.9(2) \n", "DMC: \t 0.0(3) \t -0.9(3) \t 2.9(4) \t 0.4(3) \n", "\n" ] } ], "source": [ "#dissociation energy differences to casino small Jast\n", "print(\"\\t cas sJ \\t cas lJ \\t qwalk sJ \\t qwalk lJ\")\n", "De=np.zeros((len(compoundListMolShort),8))\n", "for comp in range(len(compoundListMolShort)):\n", " atoms=re.findall('[A-Z][^A-Z]*', compoundListMolShort[comp])\n", " a1=compoundListShort.index(atoms[0])\n", " a2=compoundListShort.index(atoms[1])\n", " m=compoundListShort.index(compoundListMolShort[comp])\n", " De[comp,:]=-(energy[m,:]-energy[a1,:]-energy[a2,:])\n", " print(compoundListMolShort[comp])\n", " print(\"VMC: \\t{:10.1f}({:1.0f}) \\t{:10.1f}({:1.0f}) \\t{:10.1f}({:1.0f}) \\t{:10.1f}({:1.0f}) \".format(De[comp,0]-De[comp,0], sqrt(error[m,0]**2+error[a1,0]**2+error[a2,0]**2)*10, De[comp,2]-De[comp,0], sqrt(error[m,2]**2+error[a1,2]**2+error[a2,2]**2)*10, De[comp,4]-De[comp,0], sqrt(error[m,4]**2+error[a1,4]**2+error[a2,4]**2)*10, De[comp,6]-De[comp,0], sqrt(error[m,6]**2+error[a1,6]**2+error[a2,6]**2)*10))\n", " print(\"DMC: \\t{:10.1f}({:1.0f}) \\t{:10.1f}({:1.0f}) \\t{:10.1f}({:1.0f}) \\t{:10.1f}({:1.0f}) \".format(De[comp,1]-De[comp,1], sqrt(error[m,1]**2+error[a1,1]**2+error[a2,1]**2)*10, De[comp,3]-De[comp,1], sqrt(error[m,3]**2+error[a1,3]**2+error[a2,3]**2)*10, De[comp,5]-De[comp,1], sqrt(error[m,5]**2+error[a1,5]**2+error[a2,5]**2)*10, De[comp,7]-De[comp,1], sqrt(error[m,7]**2+error[a1,7]**2+error[a2,7]**2)*10))\n", " print(\"\")\n", " " ] }, { "cell_type": "code", "execution_count": 20, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\t cas sJ \t cas lJ \t qwalk sJ \t qwalk lJ\n", "TiCl\n", "VMC: \t 0.0(2) \t 9.9(2) \t 0.0(3) \t 1.5(2) \n", "DMC: \t 0.0(3) \t 2.6(3) \t 0.0(3) \t 1.7(3) \n", "\n", "CrCl\n", "VMC: \t 0.0(2) \t 8.9(2) \t 0.0(4) \t -2.5(2) \n", "DMC: \t 0.0(3) \t 2.4(4) \t 0.0(4) \t -0.9(3) \n", "\n", "MnS\n", "VMC: \t 0.0(2) \t -10.7(2) \t 0.0(5) \t -24.7(2) \n", "DMC: \t 0.0(3) \t -0.9(3) \t 0.0(4) \t -2.5(3) \n", "\n" ] } ], "source": [ "#dissociation energy differences to casino small Jast and qwalk small Jast seperatly\n", "print(\"\\t cas sJ \\t cas lJ \\t qwalk sJ \\t qwalk lJ\")\n", "De=np.zeros((len(compoundListMolShort),8))\n", "for comp in range(len(compoundListMolShort)):\n", " atoms=re.findall('[A-Z][^A-Z]*', compoundListMolShort[comp])\n", " a1=compoundListShort.index(atoms[0])\n", " a2=compoundListShort.index(atoms[1])\n", " m=compoundListShort.index(compoundListMolShort[comp])\n", " De[comp,:]=-(energy[m,:]-energy[a1,:]-energy[a2,:])\n", " print(compoundListMolShort[comp])\n", " print(\"VMC: \\t{:10.1f}({:1.0f}) \\t{:10.1f}({:1.0f}) \\t{:10.1f}({:1.0f}) \\t{:10.1f}({:1.0f}) \".format(De[comp,0]-De[comp,0], sqrt(error[m,0]**2+error[a1,0]**2+error[a2,0]**2)*10, De[comp,2]-De[comp,0], sqrt(error[m,2]**2+error[a1,2]**2+error[a2,2]**2)*10, De[comp,4]-De[comp,4], sqrt(error[m,4]**2+error[a1,4]**2+error[a2,4]**2)*10, De[comp,6]-De[comp,4], sqrt(error[m,6]**2+error[a1,6]**2+error[a2,6]**2)*10))\n", " print(\"DMC: \\t{:10.1f}({:1.0f}) \\t{:10.1f}({:1.0f}) \\t{:10.1f}({:1.0f}) \\t{:10.1f}({:1.0f}) \".format(De[comp,1]-De[comp,1], sqrt(error[m,1]**2+error[a1,1]**2+error[a2,1]**2)*10, De[comp,3]-De[comp,1], sqrt(error[m,3]**2+error[a1,3]**2+error[a2,3]**2)*10, De[comp,5]-De[comp,5], sqrt(error[m,5]**2+error[a1,5]**2+error[a2,5]**2)*10, De[comp,7]-De[comp,5], sqrt(error[m,7]**2+error[a1,7]**2+error[a2,7]**2)*10))\n", " print(\"\")\n", " " ] }, { "cell_type": "code", "execution_count": 21, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\t cas sJ \t cas lJ \t qwalk sJ \t qwalk lJ\n", "TiCl\n", "VMC: \t -12.9(2) \t -3.0(2) \t -8.0(3) \t -6.5(2) \n", "DMC: \t -4.7(3) \t -2.1(3) \t -3.5(3) \t -1.8(3) \n", "\n", "CrCl\n", "VMC: \t -22.0(2) \t -13.1(2) \t -9.9(4) \t -12.4(2) \n", "DMC: \t -3.6(3) \t -1.1(4) \t -1.1(4) \t -2.0(3) \n", "\n", "MnS\n", "VMC: \t -26.1(2) \t -36.7(2) \t -26.3(5) \t -51.0(2) \n", "DMC: \t -10.6(3) \t -11.5(3) \t -7.7(4) \t -10.3(3) \n", "\n" ] } ], "source": [ "#dissociation energy differences to exp\n", "print(\"\\t cas sJ \\t cas lJ \\t qwalk sJ \\t qwalk lJ\")\n", "De=np.zeros((len(compoundListMolShort),8))\n", "for comp in range(len(compoundListMolShort)):\n", " #get exp value and SO for mol\n", " experimentMol=DB3d['dimers'].find_one(compound=compoundListMolShort[comp])\n", " dissocEnergyExp[comp]=experimentMol['De(exp)[kcal/mol]']\n", " SOm=experimentMol['ESo(fromCAS)[kcal/mol]']\n", " #get atoms and mol index and SO coupling for atoms\n", " atoms=re.findall('[A-Z][^A-Z]*', compoundListMolShort[comp])\n", " a1=compoundListShort.index(atoms[0])\n", " a2=compoundListShort.index(atoms[1])\n", " m=compoundListShort.index(compoundListMolShort[comp])\n", " a1Exp=DB3d['atoms'].find_one(atom=compoundListShort[a1])\n", " SOa1=a1Exp['E_SO [kcal/mol]']\n", " a2Exp=DB3d['atoms'].find_one(atom=compoundListShort[a2])\n", " SOa2=a2Exp['E_SO [kcal/mol]']\n", " \n", " De[comp,:]=-(energy[m,:]+SOm-energy[a1,:]-energy[a2,:]-SOa1-SOa2)\n", " print(compoundListMolShort[comp])\n", " print(\"VMC: \\t{:10.1f}({:1.0f}) \\t{:10.1f}({:1.0f}) \\t{:10.1f}({:1.0f}) \\t{:10.1f}({:1.0f}) \".format(De[comp,0]-dissocEnergyExp[comp], sqrt(error[m,0]**2+error[a1,0]**2+error[a2,0]**2)*10, De[comp,2]-dissocEnergyExp[comp], sqrt(error[m,2]**2+error[a1,2]**2+error[a2,2]**2)*10, De[comp,4]-dissocEnergyExp[comp], sqrt(error[m,4]**2+error[a1,4]**2+error[a2,4]**2)*10, De[comp,6]-dissocEnergyExp[comp], sqrt(error[m,6]**2+error[a1,6]**2+error[a2,6]**2)*10))\n", " print(\"DMC: \\t{:10.1f}({:1.0f}) \\t{:10.1f}({:1.0f}) \\t{:10.1f}({:1.0f}) \\t{:10.1f}({:1.0f}) \".format(De[comp,1]-dissocEnergyExp[comp], sqrt(error[m,1]**2+error[a1,1]**2+error[a2,1]**2)*10, De[comp,3]-dissocEnergyExp[comp], sqrt(error[m,3]**2+error[a1,3]**2+error[a2,3]**2)*10, De[comp,5]-dissocEnergyExp[comp], sqrt(error[m,5]**2+error[a1,5]**2+error[a2,5]**2)*10, De[comp,7]-dissocEnergyExp[comp], sqrt(error[m,7]**2+error[a1,7]**2+error[a2,7]**2)*10))\n", " print(\"\")" ] }, { "cell_type": "code", "execution_count": 22, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "array([[-36486.52577383, -36491.30307117, -36490.81673829, -36491.24815742],\n", " [-54504.56772822, -54511.68267493, -54496.67564263, -54500.4419804 ],\n", " [-65416.59130762, -65424.18532424, -65413.17702942, -65409.17232585],\n", " [ -9386.63516885, -9390.28698662, -9386.87823709, -9387.6169665 ],\n", " [ -6351.0515223 , -6353.34058182, -6350.37067078, -6350.8665728 ],\n", " [-45961.94338972, -45980.2915879 , -45971.35724145, -45974.02701873],\n", " [-63960.09148337, -63979.72980055, -63964.57066003, -63966.60008885],\n", " [-71812.68835535, -71811.87841945, -71808.31770577, -71780.13737879]])" ] }, "execution_count": 22, "metadata": {}, "output_type": "execute_result" } ], "source": [ "energy[:,[0,2,4,6,]]" ] }, { "cell_type": "code", "execution_count": 23, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "array([[ 0.15585474, 0.14507129, 0.15550605, 0.20244914, 0.15110695,\n", " 0.193746 , 0.15847858, 0.18538343],\n", " [ 0.15453558, 0.20321073, 0.15581343, 0.2087834 , 0.25632291,\n", " 0.26122836, 0.15336928, 0.19887026],\n", " [ 0.15550274, 0.21403978, 0.15564591, 0.19829414, 0.28497524,\n", " 0.2812528 , 0.15902891, 0.195217 ],\n", " [ 0.06237133, 0.08531348, 0.05641394, 0.20778377, 0.07550463,\n", " 0.08403533, 0.06172284, 0.07592126],\n", " [ 0.0531228 , 0.07392015, 0.04585663, 0.07495823, 0.06188733,\n", " 0.08806976, 0.05556057, 0.0836905 ],\n", " [ 0.15479762, 0.19891102, 0.15714778, 0.18997519, 0.25320764,\n", " 0.21680709, 0.15143158, 0.18353363],\n", " [ 0.1546164 , 0.18564631, 0.15375001, 0.20627245, 0.33340264,\n", " 0.28034566, 0.15336338, 0.19779689],\n", " [ 0.15421479, 0.15938983, 0.15506915, 0.22263585, 0.34273235,\n", " 0.27107058, 0.15150493, 0.23287521]])" ] }, "execution_count": 23, "metadata": {}, "output_type": "execute_result" } ], "source": [ "error" ] }, { "cell_type": "code", "execution_count": 24, "metadata": { "collapsed": false }, "outputs": [ { "ename": "IndexError", "evalue": "index 3 is out of bounds for axis 0 with size 3", "output_type": "error", "traceback": [ "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[1;31mIndexError\u001b[0m Traceback (most recent call last)", "\u001b[1;32m\u001b[0m in \u001b[0;36m\u001b[1;34m()\u001b[0m\n\u001b[0;32m 6\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 7\u001b[0m \u001b[0mexperimentMol\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mDB3d\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;34m'dimers'\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mfind_one\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mcompound\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mcompoundListMol\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mcompound\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 8\u001b[1;33m \u001b[0mdissocEnergyExp\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mcompound\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mexperimentMol\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;34m'De(exp)[kcal/mol]'\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 9\u001b[0m \u001b[0mexperimentAtom\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 10\u001b[0m \u001b[1;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mlen\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0matoms\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", "\u001b[1;31mIndexError\u001b[0m: index 3 is out of bounds for axis 0 with size 3" ] } ], "source": [ "#get dissocE for DFT\n", "typ=0\n", "for compound in range(len(compoundListMol)):\n", " \n", " atoms=re.findall('[A-Z][^A-Z]*', compoundListMol[compound])\n", " \n", " experimentMol=DB3d['dimers'].find_one(compound=compoundListMol[compound])\n", " dissocEnergyExp[compound]=experimentMol['De(exp)[kcal/mol]']\n", " experimentAtom=[]\n", " for i in range(len(atoms)):\n", " experimentAtom.append(DB3d['atoms'].find_one(atom=atoms[i]))\n", " \n", " \n", " for dfttyp in range(len(dfttypList)):\n", " jobRecord={}\n", " jobRecord['dfttyp']=dfttypList[dfttyp]\n", " jobRecord['scftyp']='UHF'\n", " jobRecord['basisXMLFiles']='[\"BFD_Library_noPseudoH.xml\", \"aug-cc-pVTZ-DK_Diffuse.xml\"]'\n", " jobRecord['basisNames']='[\"vtz\", \"aug-cc-pVTZ-DK_Diffuse\"]'\n", " jobRecord['executingProgram']='gamessCasinoInt'\n", " jobRecord['compound']=compoundListMol[compound] \n", "\n", " with dataset.connect('sqlite:///'+jobDB) as DBdata:\n", " resultMol=DBdata['data'].find_one(**jobRecord)\n", " \n", " resultAtom=[]\n", " for i in range(len(atoms)):\n", " jobRecord['compound']=atoms[i]\n", " with dataset.connect('sqlite:///'+jobDB) as DBdata:\n", " resultAtom.append(DBdata['data'].find_one(**jobRecord))\n", " if resultMol['status']=='Done' and all( [resultAtom[i]['status']=='Done' for i in range(len(atoms))]):\n", " dissocEnergy[compound,dfttyp+len(dfttypList)*typ]=au2kcal(-resultMol['energy']+sum(resultAtom[i]['energy'] for i in range(len(atoms))))\n", " #correct energy for SO coupling:\n", " dissocEnergy[compound,dfttyp+len(dfttypList)*typ]+=(-experimentMol['ESo(fromCAS)[kcal/mol]']+\n", " sum(experimentAtom[i]['E_SO [kcal/mol]'] for i in range(len(atoms))))\n", " dissocError[compound,dfttyp+len(dfttypList)*typ]=0.0\n", " else:\n", " dissocEnergy[compound, dfttyp+len(dfttypList)*typ]=None\n", " dissocError[compound,dfttyp+len(dfttypList)*typ]=0.0\n", "\n", " " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "#get dissocE for emin for casino\n", "typ=1\n", "for compound in range(len(compoundListMol)):\n", " atoms=re.findall('[A-Z][^A-Z]*', compoundListMol[compound])\n", " \n", " experimentMol=DB3d['dimers'].find_one(compound=compoundListMol[compound])\n", " dissocEnergyExp[compound]=experimentMol['De(exp)[kcal/mol]']\n", " experimentAtom=[]\n", " for i in range(len(atoms)):\n", " experimentAtom.append(DB3d['atoms'].find_one(atom=atoms[i]))\n", " \n", " for dfttyp in range(len(dfttypList)):\n", " jobRecord={}\n", " jobRecord['dfttyp']=dfttypList[dfttyp]\n", " jobRecord['scftyp']='UHF'\n", " jobRecord['basisXMLFiles']='[\"BFD_Library_noPseudoH.xml\", \"aug-cc-pVTZ-DK_Diffuse.xml\"]'\n", " jobRecord['basisNames']='[\"vtz\", \"aug-cc-pVTZ-DK_Diffuse\"]'\n", " jobRecord['executingProgram']='casino'\n", " jobRecord['qmctyp']='vmc'\n", " jobRecord['compound']=compoundListMol[compound]\n", "\n", " with dataset.connect('sqlite:///'+jobDB) as DBdata:\n", " results=DBdata['data'].find(**jobRecord)\n", "\n", " for result in results:\n", " if result['jastrowFrom']==-1:\n", " continue\n", " elif '_true' in result['symmetry']:\n", " continue\n", " else:\n", " with dataset.connect('sqlite:///'+jobDB) as DBdata:\n", " jastrowFrom=DBdata['data'].find_one(id=result['jastrowFrom'])\n", " if jastrowFrom['qmctyp']!='emin':\n", " continue\n", " with dataset.connect('sqlite:///'+jobDB) as DBdata:\n", " jastrowFrom2=DBdata['data'].find_one(id=jastrowFrom['jastrowFrom'])\n", " if jastrowFrom2['qmctyp']!='varmin':\n", " continue\n", " resultMol=result\n", " \n", " resultAtom=[]\n", " for i in range(len(atoms)):\n", " jobRecord['compound']=atoms[i]\n", " with dataset.connect('sqlite:///'+jobDB) as DBdata:\n", " results=DBdata['data'].find(**jobRecord)\n", " for result in results:\n", " if result['jastrowFrom']==-1:\n", " continue\n", " elif '_true' in result['symmetry']:\n", " continue\n", " else:\n", " with dataset.connect('sqlite:///'+jobDB) as DBdata:\n", " jastrowFrom=DBdata['data'].find_one(id=result['jastrowFrom'])\n", " if jastrowFrom['qmctyp']!='emin':\n", " continue\n", " with dataset.connect('sqlite:///'+jobDB) as DBdata:\n", " jastrowFrom2=DBdata['data'].find_one(id=jastrowFrom['jastrowFrom'])\n", " if jastrowFrom2['qmctyp']!='varmin':\n", " continue\n", " resultAtom.append(result)\n", " \n", " if resultMol['status']=='Done' and all( [resultAtom[i]['status']=='Done' for i in range(len(atoms))]):\n", " dissocEnergy[compound,dfttyp+len(dfttypList)*typ]=au2kcal(-resultMol['energy']+sum(resultAtom[i]['energy'] for i in range(len(atoms))))\n", " #correct energy for SO coupling:\n", " dissocEnergy[compound,dfttyp+len(dfttypList)*typ]+=(-experimentMol['ESo(fromCAS)[kcal/mol]']+\n", " sum(experimentAtom[i]['E_SO [kcal/mol]'] for i in range(len(atoms))))\n", " dissocError[compound,dfttyp+len(dfttypList)*typ]=au2kcal(sqrt(resultMol['error']*resultMol['error']+sum(resultAtom[i]['error']*resultAtom[i]['error'] for i in range(len(atoms)))))\n", " else:\n", " dissocEnergy[compound, dfttyp+len(dfttypList)*typ]=None\n", " dissocError[compound, dfttyp+len(dfttypList)*typ]=None" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "#get dissocE for emin for qwalk\n", "typ=3\n", "for compound in range(len(compoundListMol)):\n", " print(compound)\n", " atoms=re.findall('[A-Z][^A-Z]*', compoundListMol[compound])\n", " \n", " experimentMol=DB3d['dimers'].find_one(compound=compoundListMol[compound])\n", " dissocEnergyExp[compound]=experimentMol['De(exp)[kcal/mol]']\n", " experimentAtom=[]\n", " for i in range(len(atoms)):\n", " experimentAtom.append(DB3d['atoms'].find_one(atom=atoms[i]))\n", " \n", " for dfttyp in range(len(dfttypList)):\n", " jobRecord={}\n", " jobRecord['dfttyp']=dfttypList[dfttyp]\n", " jobRecord['scftyp']='UHF'\n", " jobRecord['basisXMLFiles']='[\"BFD_Library_noPseudoH.xml\", \"aug-cc-pVTZ-DK_Diffuse.xml\"]'\n", " jobRecord['basisNames']='[\"vtz\", \"aug-cc-pVTZ-DK_Diffuse\"]'\n", " jobRecord['executingProgram']='qwalk'\n", " jobRecord['qmctyp']='vmc'\n", " jobRecord['compound']=compoundListMol[compound]\n", "\n", " with dataset.connect('sqlite:///'+jobDB) as DBdata:\n", " results=DBdata['data'].find(**jobRecord)\n", "\n", " for result in results:\n", " if result['jastrowFrom']==-1:\n", " continue\n", " elif '_true' in result['symmetry']:\n", " continue\n", " else:\n", " with dataset.connect('sqlite:///'+jobDB) as DBdata:\n", " jastrowFrom=DBdata['data'].find_one(id=result['jastrowFrom'])\n", " if jastrowFrom['qmctyp']!='emin':\n", " continue\n", " with dataset.connect('sqlite:///'+jobDB) as DBdata:\n", " jastrowFrom2=DBdata['data'].find_one(id=jastrowFrom['jastrowFrom'])\n", " if jastrowFrom2['qmctyp']!='varmin':\n", " continue\n", " resultMol=result\n", " \n", " resultAtom=[]\n", " for i in range(len(atoms)):\n", " jobRecord['compound']=atoms[i]\n", " with dataset.connect('sqlite:///'+jobDB) as DBdata:\n", " results=DBdata['data'].find(**jobRecord)\n", " for result in results:\n", " if result['jastrowFrom']==-1:\n", " continue\n", " elif '_true' in result['symmetry']:\n", " continue\n", " else:\n", " with dataset.connect('sqlite:///'+jobDB) as DBdata:\n", " jastrowFrom=DBdata['data'].find_one(id=result['jastrowFrom'])\n", " if jastrowFrom['qmctyp']!='emin':\n", " continue\n", " with dataset.connect('sqlite:///'+jobDB) as DBdata:\n", " jastrowFrom2=DBdata['data'].find_one(id=jastrowFrom['jastrowFrom'])\n", " if jastrowFrom2['qmctyp']!='varmin':\n", " continue\n", " resultAtom.append(result)\n", " \n", " if resultMol['status']=='Done' and all( [resultAtom[i]['status']=='Done' for i in range(len(resultAtom))]) and len(resultAtom)==2:\n", " dissocEnergy[compound,dfttyp+len(dfttypList)*typ]=au2kcal(-resultMol['energy']+sum(resultAtom[i]['energy'] for i in range(len(atoms))))\n", " #correct energy for SO coupling:\n", " dissocEnergy[compound,dfttyp+len(dfttypList)*typ]+=(-experimentMol['ESo(fromCAS)[kcal/mol]']+\n", " sum(experimentAtom[i]['E_SO [kcal/mol]'] for i in range(len(atoms))))\n", " dissocError[compound,dfttyp+len(dfttypList)*typ]=au2kcal(sqrt(resultMol['error']*resultMol['error']+sum(resultAtom[i]['error']*resultAtom[i]['error'] for i in range(len(atoms)))))\n", " else:\n", " dissocEnergy[compound, dfttyp+len(dfttypList)*typ]=None\n", " dissocError[compound, dfttyp+len(dfttypList)*typ]=None" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "#get dissocE for dmc (zero tau extrap) for casino\n", "typ=2\n", "for compound in range(len(compoundListMol)):\n", " atoms=re.findall('[A-Z][^A-Z]*', compoundListMol[compound])\n", " \n", " experimentMol=DB3d['dimers'].find_one(compound=compoundListMol[compound])\n", " dissocEnergyExp[compound]=experimentMol['De(exp)[kcal/mol]']\n", " experimentAtom=[]\n", " for i in range(len(atoms)):\n", " experimentAtom.append(DB3d['atoms'].find_one(atom=atoms[i]))\n", " \n", " for dfttyp in range(len(dfttypList)):\n", " jobRecord={}\n", " jobRecord['dfttyp']=dfttypList[dfttyp]\n", " jobRecord['scftyp']='UHF'\n", " jobRecord['basisXMLFiles']='[\"BFD_Library_noPseudoH.xml\", \"aug-cc-pVTZ-DK_Diffuse.xml\"]'\n", " jobRecord['basisNames']='[\"vtz\", \"aug-cc-pVTZ-DK_Diffuse\"]'\n", " jobRecord['executingProgram']='casino'\n", " jobRecord['qmctyp']='vmc_dmc'\n", " jobRecord['compound']=compoundListMol[compound]\n", "\n", " with dataset.connect('sqlite:///'+jobDB) as DBdata:\n", " results=DBdata['data'].find(**jobRecord)\n", "\n", " resultsExtrap=[]\n", " resultMol={}\n", " resultMolQuad={}\n", " for result in results:\n", " if result['jastrowFrom']==-1:\n", " continue\n", " elif '_true' in result['symmetry']:\n", " continue\n", " else:\n", " with dataset.connect('sqlite:///'+jobDB) as DBdata:\n", " jastrowFrom=DBdata['data'].find_one(id=result['jastrowFrom'])\n", " if jastrowFrom['qmctyp']!='emin':\n", " continue\n", " with dataset.connect('sqlite:///'+jobDB) as DBdata:\n", " jastrowFrom2=DBdata['data'].find_one(id=jastrowFrom['jastrowFrom'])\n", " if jastrowFrom2['qmctyp']!='varmin':\n", " continue\n", " resultsExtrap.append(result)\n", " for i in reversed(range(len(resultsExtrap))):\n", " if not resultsExtrap[i]['status'] in [\"Done\",\"LargeErrorbar2\", \"BadReblock\"]:\n", " resultsExtrap.pop(i)\n", "\n", " if len(resultsExtrap)!=3:\n", " resultMol['status']='Error'\n", " else: #extrapolate to zero Time\n", " [resultMol['energy'], resultMol['error']]=\\\n", " extapolateTau([resultsExtrap[i]['dt'] for i in range(len(resultsExtrap))],\n", " [resultsExtrap[i]['energy'] for i in range(len(resultsExtrap))],\n", " [resultsExtrap[i]['error'] for i in range(len(resultsExtrap))])\n", " resultMol['status']='Done'\n", " if len(resultsExtrap)<3:\n", " resultMolQuad['status']='Error'\n", " else:\n", " [resultMolQuad['energy'], resultMolQuad['error']]=\\\n", " extapolateTauQuad([resultsExtrap[i]['dt'] for i in range(len(resultsExtrap))],\n", " [resultsExtrap[i]['energy'] for i in range(len(resultsExtrap))],\n", " [resultsExtrap[i]['error'] for i in range(len(resultsExtrap))])\n", " resultMolQuad['status']='Done'\n", " \n", " resultAtom=[]\n", " resultAtomQuad=[]\n", " for i in range(len(atoms)):\n", " resultsExtrap=[]\n", " resultAtomI={}\n", " resultAtomIQuad={}\n", " jobRecord['compound']=atoms[i]\n", " with dataset.connect('sqlite:///'+jobDB) as DBdata:\n", " results=DBdata['data'].find(**jobRecord)\n", " for result in results:\n", " if result['jastrowFrom']==-1:\n", " continue\n", " elif '_true' in result['symmetry']:\n", " continue\n", " else:\n", " with dataset.connect('sqlite:///'+jobDB) as DBdata:\n", " jastrowFrom=DBdata['data'].find_one(id=result['jastrowFrom'])\n", " if jastrowFrom['qmctyp']!='emin':\n", " continue\n", " with dataset.connect('sqlite:///'+jobDB) as DBdata:\n", " jastrowFrom2=DBdata['data'].find_one(id=jastrowFrom['jastrowFrom'])\n", " if jastrowFrom2['qmctyp']!='varmin':\n", " continue\n", " resultsExtrap.append(result)\n", " for i in reversed(range(len(resultsExtrap))):\n", " if not resultsExtrap[i]['status'] in [\"Done\",\"LargeErrorbar2\", \"BadReblock\"]:\n", " resultsExtrap.pop(i)\n", "\n", " if len(resultsExtrap)!=3:\n", " resultAtomI['status']='Error'\n", " resultAtom.append(resultAtomI)\n", " else: #extrapolate to zero Time\n", " [resultAtomI['energy'], resultAtomI['error']]=\\\n", " extapolateTau([resultsExtrap[i]['dt'] for i in range(len(resultsExtrap))],\n", " [resultsExtrap[i]['energy'] for i in range(len(resultsExtrap))],\n", " [resultsExtrap[i]['error'] for i in range(len(resultsExtrap))])\n", " resultAtomI['status']='Done' \n", " resultAtom.append(resultAtomI)\n", " if len(resultsExtrap)<3:\n", " resultAtomIQuad['status']='Error'\n", " resultAtomQuad.append(resultAtomIQuad)\n", " else:\n", " [resultAtomIQuad['energy'], resultAtomIQuad['error']]=\\\n", " extapolateTauQuad([resultsExtrap[i]['dt'] for i in range(len(resultsExtrap))],\n", " [resultsExtrap[i]['energy'] for i in range(len(resultsExtrap))],\n", " [resultsExtrap[i]['error'] for i in range(len(resultsExtrap))])\n", " resultAtomIQuad['status']='Done' \n", " resultAtomQuad.append(resultAtomIQuad)\n", " #print(compoundListMol[compound])\n", " if resultMol['status']=='Done' and all( [resultAtom[i]['status']=='Done' for i in range(len(atoms))]):\n", " dissocEnergy[compound,dfttyp+len(dfttypList)*typ]=au2kcal(-resultMol['energy']+sum(resultAtom[i]['energy'] for i in range(len(atoms))))\n", " #correct energy for SO coupling:\n", " dissocEnergy[compound,dfttyp+len(dfttypList)*typ]+=(-experimentMol['ESo(fromCAS)[kcal/mol]']+\n", " sum(experimentAtom[i]['E_SO [kcal/mol]'] for i in range(len(atoms))))\n", " dissocError[compound,dfttyp+len(dfttypList)*typ]=au2kcal(sqrt(resultMol['error']*resultMol['error']+sum(resultAtom[i]['error']*resultAtom[i]['error'] for i in range(len(atoms)))))\n", " else:\n", " #print(\"{} {} {}\".format(resultMol['status'], resultAtom[0]['status'], resultAtom[1]['status']))\n", " dissocEnergy[compound, dfttyp+len(dfttypList)*typ]=None\n", " dissocError[compound, dfttyp+len(dfttypList)*typ]=None\n", " if resultMolQuad['status']=='Done' and all( [resultAtomQuad[i]['status']=='Done' for i in range(len(atoms))]): \n", " dissocEnergyQuad[compound,dfttyp+len(dfttypList)*typ]=au2kcal(-resultMolQuad['energy']+sum(resultAtomQuad[i]['energy'] for i in range(len(atoms))))\n", " #correct energy for SO coupling:\n", " dissocEnergyQuad[compound,dfttyp+len(dfttypList)*typ]+=(-experimentMol['ESo(fromCAS)[kcal/mol]']+\n", " sum(experimentAtom[i]['E_SO [kcal/mol]'] for i in range(len(atoms))))\n", " dissocErrorQuad[compound,dfttyp+len(dfttypList)*typ]=au2kcal(sqrt(resultMolQuad['error']*resultMolQuad['error']+sum(resultAtomQuad[i]['error']*resultAtomQuad[i]['error'] for i in range(len(atoms)))))\n", " else:\n", " dissocEnergyQuad[compound,dfttyp+len(dfttypList)*typ]=None\n", " dissocErrorQuad[compound,dfttyp+len(dfttypList)*typ]=None" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "#get dissocE for dmc (zero tau extrap) for qwalk\n", "typ=4\n", "for compound in range(len(compoundListMol)):\n", " atoms=re.findall('[A-Z][^A-Z]*', compoundListMol[compound])\n", " \n", " experimentMol=DB3d['dimers'].find_one(compound=compoundListMol[compound])\n", " dissocEnergyExp[compound]=experimentMol['De(exp)[kcal/mol]']\n", " experimentAtom=[]\n", " for i in range(len(atoms)):\n", " experimentAtom.append(DB3d['atoms'].find_one(atom=atoms[i]))\n", " \n", " for dfttyp in range(len(dfttypList)):\n", " jobRecord={}\n", " jobRecord['dfttyp']=dfttypList[dfttyp]\n", " jobRecord['scftyp']='UHF'\n", " jobRecord['basisXMLFiles']='[\"BFD_Library_noPseudoH.xml\", \"aug-cc-pVTZ-DK_Diffuse.xml\"]'\n", " jobRecord['basisNames']='[\"vtz\", \"aug-cc-pVTZ-DK_Diffuse\"]'\n", " jobRecord['executingProgram']='qwalk'\n", " jobRecord['qmctyp']='vmc_dmc'\n", " jobRecord['compound']=compoundListMol[compound]\n", "\n", " with dataset.connect('sqlite:///'+jobDB) as DBdata:\n", " results=DBdata['data'].find(**jobRecord)\n", "\n", " resultsExtrap=[]\n", " resultMol={}\n", " resultMolQuad={}\n", " for result in results:\n", " if result['jastrowFrom']==-1:\n", " continue\n", " elif '_true' in result['symmetry']:\n", " continue\n", " else:\n", " with dataset.connect('sqlite:///'+jobDB) as DBdata:\n", " jastrowFrom=DBdata['data'].find_one(id=result['jastrowFrom'])\n", " if jastrowFrom['qmctyp']!='emin':\n", " continue\n", " with dataset.connect('sqlite:///'+jobDB) as DBdata:\n", " jastrowFrom2=DBdata['data'].find_one(id=jastrowFrom['jastrowFrom'])\n", " if jastrowFrom2['qmctyp']!='varmin':\n", " continue\n", " resultsExtrap.append(result)\n", " for i in reversed(range(len(resultsExtrap))):\n", " if not resultsExtrap[i]['status'] in [\"Done\",\"LargeErrorbar2\", \"BadReblock\"]:\n", " resultsExtrap.pop(i)\n", "\n", " if len(resultsExtrap)!=3:\n", " resultMol['status']='Error'\n", " else: #extrapolate to zero Time\n", " [resultMol['energy'], resultMol['error']]=\\\n", " extapolateTau([resultsExtrap[i]['dt'] for i in range(len(resultsExtrap))],\n", " [resultsExtrap[i]['energy'] for i in range(len(resultsExtrap))],\n", " [resultsExtrap[i]['error'] for i in range(len(resultsExtrap))])\n", " resultMol['status']='Done'\n", " if len(resultsExtrap)<3:\n", " resultMolQuad['status']='Error'\n", " else:\n", " [resultMolQuad['energy'], resultMolQuad['error']]=\\\n", " extapolateTauQuad([resultsExtrap[i]['dt'] for i in range(len(resultsExtrap))],\n", " [resultsExtrap[i]['energy'] for i in range(len(resultsExtrap))],\n", " [resultsExtrap[i]['error'] for i in range(len(resultsExtrap))])\n", " resultMolQuad['status']='Done'\n", " \n", " resultAtom=[]\n", " resultAtomQuad=[]\n", " for i in range(len(atoms)):\n", " resultsExtrap=[]\n", " resultAtomI={}\n", " resultAtomIQuad={}\n", " jobRecord['compound']=atoms[i]\n", " with dataset.connect('sqlite:///'+jobDB) as DBdata:\n", " results=DBdata['data'].find(**jobRecord)\n", " for result in results:\n", " if result['jastrowFrom']==-1:\n", " continue\n", " elif '_true' in result['symmetry']:\n", " continue\n", " else:\n", " with dataset.connect('sqlite:///'+jobDB) as DBdata:\n", " jastrowFrom=DBdata['data'].find_one(id=result['jastrowFrom'])\n", " if jastrowFrom['qmctyp']!='emin':\n", " continue\n", " with dataset.connect('sqlite:///'+jobDB) as DBdata:\n", " jastrowFrom2=DBdata['data'].find_one(id=jastrowFrom['jastrowFrom'])\n", " if jastrowFrom2['qmctyp']!='varmin':\n", " continue\n", " resultsExtrap.append(result)\n", " for i in reversed(range(len(resultsExtrap))):\n", " if not resultsExtrap[i]['status'] in [\"Done\",\"LargeErrorbar2\", \"BadReblock\"]:\n", " resultsExtrap.pop(i)\n", "\n", " if len(resultsExtrap)!=3:\n", " resultAtomI['status']='Error'\n", " resultAtom.append(resultAtomI)\n", " else: #extrapolate to zero Time\n", " [resultAtomI['energy'], resultAtomI['error']]=\\\n", " extapolateTau([resultsExtrap[i]['dt'] for i in range(len(resultsExtrap))],\n", " [resultsExtrap[i]['energy'] for i in range(len(resultsExtrap))],\n", " [resultsExtrap[i]['error'] for i in range(len(resultsExtrap))])\n", " resultAtomI['status']='Done' \n", " resultAtom.append(resultAtomI)\n", " if len(resultsExtrap)<3:\n", " resultAtomIQuad['status']='Error'\n", " resultAtomQuad.append(resultAtomIQuad)\n", " else:\n", " [resultAtomIQuad['energy'], resultAtomIQuad['error']]=\\\n", " extapolateTauQuad([resultsExtrap[i]['dt'] for i in range(len(resultsExtrap))],\n", " [resultsExtrap[i]['energy'] for i in range(len(resultsExtrap))],\n", " [resultsExtrap[i]['error'] for i in range(len(resultsExtrap))])\n", " resultAtomIQuad['status']='Done' \n", " resultAtomQuad.append(resultAtomIQuad)\n", " #print(compoundListMol[compound])\n", " if resultMol['status']=='Done' and all( [resultAtom[i]['status']=='Done' for i in range(len(atoms))]):\n", " dissocEnergy[compound,dfttyp+len(dfttypList)*typ]=au2kcal(-resultMol['energy']+sum(resultAtom[i]['energy'] for i in range(len(atoms))))\n", " #correct energy for SO coupling:\n", " dissocEnergy[compound,dfttyp+len(dfttypList)*typ]+=(-experimentMol['ESo(fromCAS)[kcal/mol]']+\n", " sum(experimentAtom[i]['E_SO [kcal/mol]'] for i in range(len(atoms))))\n", " dissocError[compound,dfttyp+len(dfttypList)*typ]=au2kcal(sqrt(resultMol['error']*resultMol['error']+sum(resultAtom[i]['error']*resultAtom[i]['error'] for i in range(len(atoms)))))\n", " else:\n", " #print(\"{} {} {}\".format(resultMol['status'], resultAtom[0]['status'], resultAtom[1]['status']))\n", " dissocEnergy[compound, dfttyp+len(dfttypList)*typ]=None\n", " dissocError[compound, dfttyp+len(dfttypList)*typ]=None\n", " if resultMolQuad['status']=='Done' and all( [resultAtomQuad[i]['status']=='Done' for i in range(len(atoms))]): \n", " dissocEnergyQuad[compound,dfttyp+len(dfttypList)*typ]=au2kcal(-resultMolQuad['energy']+sum(resultAtomQuad[i]['energy'] for i in range(len(atoms))))\n", " #correct energy for SO coupling:\n", " dissocEnergyQuad[compound,dfttyp+len(dfttypList)*typ]+=(-experimentMol['ESo(fromCAS)[kcal/mol]']+\n", " sum(experimentAtom[i]['E_SO [kcal/mol]'] for i in range(len(atoms))))\n", " dissocErrorQuad[compound,dfttyp+len(dfttypList)*typ]=au2kcal(sqrt(resultMolQuad['error']*resultMolQuad['error']+sum(resultAtomQuad[i]['error']*resultAtomQuad[i]['error'] for i in range(len(atoms)))))\n", " else:\n", " dissocEnergyQuad[compound,dfttyp+len(dfttypList)*typ]=None\n", " dissocErrorQuad[compound,dfttyp+len(dfttypList)*typ]=None" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "#plot dissoc energy results\n", "typForPlot=[0,2,4]\n", "molPerPlot=min(11, len(compoundListMol))\n", "colors=['r','g','b','y','c','m']\n", "edgecolors=['r','g','b','y','c','m']\n", "opac=[0.3,1,0.5]\n", "error_config = {'ecolor': '0.3'}\n", "barwidth=1/(4*len(typForPlot)*1.2)\n", "#barwidth=0.2\n", "index = numpy.arange(molPerPlot)\n", "\n", "for l in range(int(ceil(float(len(compoundListMol))/molPerPlot))):#range(int(len(compoundListMol)/molPerPlot)+1):\n", " a=l*(molPerPlot-1)\n", " b=min((l+1)*(molPerPlot-1), len(compoundListMol))\n", " fig, ax = plt.subplots(figsize=(13, 6))\n", " for i in range(len(typForPlot)):\n", " for j in range(4):\n", " if i==0 and j==0:\n", " continue\n", " rects=plt.bar(index[0:b-a]+(i*(4)+j)*barwidth+0.1*4*len(typForPlot)*barwidth*i,\n", " (dissocEnergy[a:b,typForPlot[i]*(len(dfttypList))+j]-dissocEnergyExp[a:b]), \n", " barwidth,\n", " label=typList[typForPlot[i]]+\"(\"+dfttypList[j]+\")\",\n", " color=colors[j],\n", " alpha=opac[i],\n", " linewidth=1,\n", " yerr=dissocError[a:b,typForPlot[i]*(len(dfttypList))+j],\n", " error_kw=error_config,)\n", " for k in range(b-a):\n", " plt.axvline(x=k, color='0.85')\n", " plt.xticks(index+0.4, compoundListMol[a:b])\n", " #plt.xlabel('Molecules, Methods') \n", " plt.ylabel(r'$E_\\mathrm{dissoc.}-E_\\mathrm{dissoc.}^\\mathrm{exp}$', fontsize=14) \n", " if l==0:\n", " #plt.title('dissoc. energies grouped by molecules and methods') \n", " plt.legend(loc=2,bbox_to_anchor=(1.05, 1), prop={'size':10})#(loc='best',prop={'size':8})\n", " plt.autoscale( enable=True,axis=u'both', tight=True)\n", " ax.text(real(molPerPlot-4)/2.0, 15, \"overbinding\", fontsize=10,color='green',fontweight='bold')\n", " ax.text(real(molPerPlot-4)/2.0, -15, \"underbinding\", fontsize=10,color='green',fontweight='bold')\n", " ylim([-12,12])\n", " xlim([0,molPerPlot-1])\n", " fig.set_size_inches(12,2.5)\n", " plt.savefig('dissocE'+str(l)+'.svg', format='svg', dpi=300)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "#plot difference in dissoc energy results casino vs. qwalk for DMC(B3LYP)\n", "typForPlot=[2]\n", "RefTypForPlot=[4]\n", "molPerPlot=min(20, len(compoundListMol))\n", "colors=['r','g','b','y','c','m']\n", "edgecolors=['r','g','b','y','c','m']\n", "opac=[0.3,1,0.5]\n", "error_config = {'ecolor': '0.3'}\n", "barwidth=1/(1.2)\n", "#barwidth=0.2\n", "index = numpy.arange(molPerPlot)\n", "\n", "for l in range(int(ceil(float(len(compoundListMol))/molPerPlot))):#range(int(len(compoundListMol)/molPerPlot)+1):\n", " a=l*(molPerPlot-1)\n", " b=min((l+1)*(molPerPlot-1), len(compoundListMol))\n", " fig, ax = plt.subplots(figsize=(13, 6))\n", " for i in range(len(typForPlot)):\n", " for j in [2]:\n", " if i==0 and j==0:\n", " continue\n", " rects=plt.bar(index[0:b-a]+(i*(4))*barwidth+0.1*4*len(typForPlot)*barwidth*i,\n", " (dissocEnergy[a:b,typForPlot[i]*(len(dfttypList))+j]-dissocEnergy[a:b,RefTypForPlot[i]*(len(dfttypList))+j]), \n", " barwidth,\n", " label=typList[typForPlot[i]]+\"(\"+dfttypList[j]+\")\",\n", " color=colors[j],\n", " alpha=opac[i],\n", " linewidth=1,\n", " yerr=(dissocError[a:b,typForPlot[i]*(len(dfttypList))+j]+dissocError[a:b,RefTypForPlot[i]*(len(dfttypList))+j]),\n", " error_kw=error_config,)\n", " for k in range(b-a):\n", " plt.axvline(x=k, color='0.85')\n", " plt.xticks(index+0.4, compoundListMol[a:b])\n", " #plt.xlabel('Molecules, Methods') \n", " plt.ylabel(r'$E_\\mathrm{dissoc.}^\\mathrm{casino}-E_\\mathrm{dissoc.}^\\mathrm{qwalk}$', fontsize=14) \n", " if l==0:\n", " #plt.title('dissoc. energies grouped by molecules and methods') \n", " plt.legend(loc=2,bbox_to_anchor=(1.05, 1), prop={'size':10})#(loc='best',prop={'size':8})\n", " plt.autoscale( enable=True,axis=u'both', tight=True)\n", " ylim([-2.6,2.6])\n", " xlim([0,molPerPlot-1])\n", " fig.set_size_inches(12,2.5)\n", " plt.savefig('dissocE_casinoVsQwalk'+str(l)+'.svg', format='svg', dpi=300)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "#plot dissoc energy results casino vs qwalk for B3LYP\n", "typForPlot=[2,4]\n", "molPerPlot=min(11, len(compoundListMol))\n", "colors=['r','g','b','y','c','m']\n", "edgecolors=['r','g','b','y','c','m']\n", "opac=[1,0.5]\n", "error_config = {'ecolor': '0.3'}\n", "barwidth=1/(len(typForPlot)*1.2)\n", "#barwidth=0.2\n", "index = numpy.arange(molPerPlot)\n", "\n", "for l in range(int(ceil(float(len(compoundListMol))/molPerPlot))):#range(int(len(compoundListMol)/molPerPlot)+1):\n", " a=l*(molPerPlot-1)\n", " b=min((l+1)*(molPerPlot-1), len(compoundListMol))\n", " fig, ax = plt.subplots(figsize=(13, 6))\n", " for i in range(len(typForPlot)):\n", " for j in [2]:\n", " if i==0 and j==0:\n", " continue\n", " rects=plt.bar(index[0:b-a]+(i)*barwidth+0.1*len(typForPlot)*barwidth*i,\n", " (dissocEnergy[a:b,typForPlot[i]*(len(dfttypList))+j]-dissocEnergyExp[a:b]), \n", " barwidth,\n", " label=typList[typForPlot[i]]+\"(\"+dfttypList[j]+\")\",\n", " color=colors[j],\n", " alpha=opac[i],\n", " linewidth=1,\n", " yerr=dissocError[a:b,typForPlot[i]*(len(dfttypList))+j],\n", " error_kw=error_config,)\n", " for k in range(b-a):\n", " plt.axvline(x=k, color='0.85')\n", " plt.xticks(index+0.4, compoundListMol[a:b])\n", " #plt.xlabel('Molecules, Methods') \n", " plt.ylabel(r'$E_\\mathrm{dissoc.}-E_\\mathrm{dissoc.}^\\mathrm{exp}$', fontsize=14) \n", " if l==0:\n", " #plt.title('dissoc. energies grouped by molecules and methods') \n", " plt.legend(loc=2,bbox_to_anchor=(1.05, 1), prop={'size':10})#(loc='best',prop={'size':8})\n", " plt.autoscale( enable=True,axis=u'both', tight=True)\n", " ax.text(real(molPerPlot-4)/2.0, 15, \"overbinding\", fontsize=10,color='green',fontweight='bold')\n", " ax.text(real(molPerPlot-4)/2.0, -15, \"underbinding\", fontsize=10,color='green',fontweight='bold')\n", " ylim([-11,7])\n", " xlim([0,molPerPlot-1])\n", " fig.set_size_inches(12,2.5)\n", " plt.savefig('dissocE'+str(l)+'.svg', format='svg', dpi=300)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "#print the energy differences for DMC for casino\n", "typ=2\n", "for compound in range(len(compoundListMol)):\n", " print(compoundListMol[compound])\n", " print(dissocEnergy[compound,typ*(len(dfttypList)):typ*len(dfttypList)+4]-dissocEnergyExp[compound])\n", " print(dissocError[compound,typ*(len(dfttypList)):typ*len(dfttypList)+4])" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "#print the energy differences for DMC for qwalk\n", "typ=4\n", "for compound in range(len(compoundListMol)):\n", " print(compoundListMol[compound])\n", " print(dissocEnergy[compound,typ*(len(dfttypList)):typ*len(dfttypList)+4]-dissocEnergyExp[compound])\n", " print(dissocError[compound,typ*(len(dfttypList)):typ*len(dfttypList)+4])" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "#print the energy differences for DFT\n", "typ=0\n", "for compound in range(len(compoundListMol)):\n", " print(compoundListMol[compound])\n", " print(dissocEnergy[compound,typ*(len(dfttypList)):typ*len(dfttypList)+4]-dissocEnergyExp[compound])\n", " print(dissocError[compound,typ*(len(dfttypList)):typ*len(dfttypList)+4])" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "#calculate variance in the dissoc energy for DFT and DMC after DFT (without HF):\n", "var=zeros(5)\n", "\n", "for typ in [0,2,4]:\n", " n=0\n", " molNan=0\n", " for compound in range(len(compoundListMol)):\n", " if any(np.isnan(dissocEnergy[compound, 1:4])) or any(np.isnan(dissocEnergy[compound, 11:14])):\n", " molNan+=1\n", " continue\n", " else:\n", " var[typ]+=sum((dissocEnergy[compound, 1+typ*len(dfttypList):4+typ*len(dfttypList)]\n", " -mean(dissocEnergy[compound, 1+typ*len(dfttypList):4+typ*len(dfttypList)]))**2)\n", " n+=1\n", " print(compoundListMol[compound])\n", " print(sqrt(sum((dissocEnergy[compound, 1+typ*len(dfttypList):4+typ*len(dfttypList)]\n", " -mean(dissocEnergy[compound, 1+typ*len(dfttypList):4+typ*len(dfttypList)]))**2)/3))\n", " var[typ]=var[typ]/(n*3-n)\n", "\n", "print(var)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "#calculate variance in the dissoc energy for DFT and DMC after DFT (without HF) and without VO, VH, and FeH, FeCl, CoH, CoCl(bad PBE):\n", "var=zeros(5)\n", "for typ in [0,2,4]:\n", " n=0\n", " molNan=0\n", " for compound in range(len(compoundListMol)):\n", " if any(np.isnan(dissocEnergy[compound, 1:4])) or any(np.isnan(dissocEnergy[compound, 11:14])):\n", " molNan+=1\n", " continue\n", " if 'V' in compoundListMol[compound] or 'Fe' in compoundListMol[compound] or 'Co' in compoundListMol[compound]:\n", " continue\n", " else:\n", " var[typ]+=sum((dissocEnergy[compound, 1+typ*len(dfttypList):4+typ*len(dfttypList)]\n", " -mean(dissocEnergy[compound, 1+typ*len(dfttypList):4+typ*len(dfttypList)]))**2)\n", " n+=1\n", " #print(compoundListMol[compound])\n", " #print(sqrt(sum((dissocEnergy[compound, 1+typ*len(dfttypList):4+typ*len(dfttypList)]\n", " # -mean(dissocEnergy[compound, 1+typ*len(dfttypList):4+typ*len(dfttypList)]))**2)/3))\n", " print(n)\n", " var[typ]=var[typ]/(n*3-n)\n", "\n", "print(var)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "#analyze errorbars of dissocEnergy casino\n", "typ=2\n", "boxplot(dissocError[:,typ*len(dfttypList):4+typ*len(dfttypList)])\n", "plt.xticks([1,2,3,4], dfttypList);" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "#analyze errorbars of dissocEnergy qwalk\n", "typ=4\n", "boxplot(dissocError[:,typ*len(dfttypList):4+typ*len(dfttypList)])\n", "plt.xticks([1,2,3,4], dfttypList);" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "#calculate all mean unsigned errors\n", "for i in [0,2,4]:#range(len(typList)):\n", " for j in range(len(dfttypList)):\n", " print(\"{} {}: {}\".format(typList[i],dfttypList[j],nanmean(abs(dissocEnergy[:,j+i*len(dfttypList)]-dissocEnergyExp))))\n", " print(np.count_nonzero(~np.isnan(dissocEnergy[:,j+i*len(dfttypList)])))\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "#calculate all mean signed errors\n", "for i in [0,2,4]:#range(len(typList)):\n", " for j in range(len(dfttypList)):\n", " print(\"{} {}: {}\".format(typList[i],dfttypList[j],nanmean(dissocEnergy[:,j+i*len(dfttypList)]-dissocEnergyExp)))\n", " print(np.count_nonzero(~np.isnan(dissocEnergy[:,j+i*len(dfttypList)])))\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "#calculate standard deviation of the mean of the errors\n", "for i in [0,2,4]:#range(len(typList)):\n", " for j in range(len(dfttypList)-1):\n", " mue=nanmean(abs(dissocEnergy[:,j+i*len(dfttypList)]-dissocEnergyExp))\n", " n=0\n", " var=0\n", " for k in range(len(compoundListMol)):\n", " if isnan(dissocEnergy[k,j+i*len(dfttypList)]):\n", " continue\n", " n+=1\n", " var+=(abs(dissocEnergy[k,j+i*len(dfttypList)]-dissocEnergyExp[k])-mue)**2\n", " if not n==0:\n", " print(\"{} {}: {}\".format(typList[i],dfttypList[j],sqrt(var/(n-1)/n)))\n", " print(np.count_nonzero(~np.isnan(dissocEnergy[:,j+i*len(dfttypList)])))\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "#calculate all mean errors for the mols which are converged in HF, PBE, B3LY and B97-1 and without V..., Fe..., and Co... (bad PBE)\n", "dissocEnergyAllConv=zeros([len(compoundListMol), 5*len(dfttypList)])\n", "for i in range(len(compoundListMol)):\n", " if any([isnan(dissocEnergy[i, j]) for j in range(4)]) or any([isnan(dissocEnergy[i,2*len(dfttypList)+j:2*len(dfttypList)+j]) for j in range(4)]) or any([isnan(dissocEnergy[i,4*len(dfttypList)+j:4*len(dfttypList)+j]) for j in range(4)]) or 'V' in compoundListMol[i] or 'Fe' in compoundListMol[i] or 'Co' in compoundListMol[i]:\n", " dissocEnergyAllConv[i, 0:4]=None\n", " dissocEnergyAllConv[i,2*len(dfttypList):2*len(dfttypList)+4]=None\n", " dissocEnergyAllConv[i,4*len(dfttypList):4*len(dfttypList)+4]=None\n", " else:\n", " dissocEnergyAllConv[i, 0:4]=dissocEnergy[i, 0:4]\n", " dissocEnergyAllConv[i,2*len(dfttypList):2*len(dfttypList)+4]=dissocEnergy[i,2*len(dfttypList):2*len(dfttypList)+4]\n", " dissocEnergyAllConv[i,4*len(dfttypList):4*len(dfttypList)+4]=dissocEnergy[i,4*len(dfttypList):4*len(dfttypList)+4]\n", "for i in [0,2,4]:#range(len(typList)):\n", " for j in range(4):\n", " print(\"{} {}: {}\".format(typList[i],dfttypList[j],nanmean(abs(dissocEnergyAllConv[:,j+i*len(dfttypList)]-dissocEnergyExp))))\n", " print(np.count_nonzero(~np.isnan(dissocEnergyAllConv[:,j+i*len(dfttypList)])))" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "#calculate standard deviation of the mean of the errors without V..., Fe..., and Co...(bad PBE)\n", "for i in [0,2,4]:#range(len(typList)):\n", " for j in range(len(dfttypList)-1):\n", " mue=nanmean(abs(dissocEnergy[:,j+i*len(dfttypList)]-dissocEnergyExp))\n", " n=0\n", " var=0\n", " for k in range(len(compoundListMol)):\n", " if isnan(dissocEnergy[k,j+i*len(dfttypList)])or any([isnan(dissocEnergy[k,2*len(dfttypList)+j:2*len(dfttypList)+j]) for j in range(4)]) or 'V' in compoundListMol[k] or 'Fe' in compoundListMol[k] or 'Co' in compoundListMol[k]:\n", " continue\n", " n+=1\n", " var+=(abs(dissocEnergy[k,j+i*len(dfttypList)]-dissocEnergyExp[k])-mue)**2\n", " if not n==0:\n", " print(\"{} {}: {}\".format(typList[i],dfttypList[j],sqrt(var/(n-1)/n)))\n", " print(n)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "#calculate mean errors for single Det character mols:\n", "for i in [0,2,4]:#range(len(typList)):\n", " for j in range(len(dfttypList)):\n", " print(\"{} {}: {}\".format(typList[i],dfttypList[j],nanmean([abs(dissocEnergy[k,j+i*len(dfttypList)]-dissocEnergyExp[k]) for k in range(len(compoundListMol)) if compoundListMol[k] in ['CuCl', 'ZnH', 'ZnS', 'ZnCl', 'CrCl', 'MnCl', 'FeCl']])))\n", " print(np.count_nonzero(~np.isnan([dissocEnergy[k,j+i*len(dfttypList)]for k in range(len(compoundListMol)) if compoundListMol[k] in ['CuCl', 'ZnH', 'ZnS', 'ZnCl', 'CrCl', 'MnCl', 'FeCl']])))\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "#calculate mean errors for multi Det character mols:\n", "for i in [0,2,4]:#range(len(typList)):\n", " for j in range(len(dfttypList)):\n", " print(\"{} {}: {}\".format(typList[i],dfttypList[j],nanmean([abs(dissocEnergy[k,j+i*len(dfttypList)]-dissocEnergyExp[k]) for k in range(len(compoundListMol)) if not compoundListMol[k] in ['CuCl', 'ZnH', 'ZnS', 'ZnCl', 'CrCl', 'MnCl', 'FeCl']])))\n", " print(np.count_nonzero(~np.isnan([dissocEnergy[k,j+i*len(dfttypList)]for k in range(len(compoundListMol)) if not compoundListMol[k] in ['CuCl', 'ZnH', 'ZnS', 'ZnCl', 'CrCl', 'MnCl', 'FeCl']])))\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "#calculate mean errors for mols with exp error <1.2kcal/mol (['FeH', 'NiCl', 'CuCl', 'CoH', 'CrO', ZnH', 'ZnO', 'ZnS', 'ZnCl'])\n", "for i in [0,2,4]:#range(len(typList)):\n", " for j in range(len(dfttypList)):\n", " print(\"{} {}: {}\".format(typList[i],dfttypList[j],nanmean([abs(dissocEnergy[k,j+i*len(dfttypList)]-dissocEnergyExp[k]) for k in range(len(compoundListMol)) if not compoundListMol[k] in ['FeH', 'NiCl', 'CuCl', 'CoH', 'CrO','ZnH', 'ZnO', 'ZnS', 'ZnCl']])))\n", " print(np.count_nonzero(~np.isnan([dissocEnergy[k,j+i*len(dfttypList)]for k in range(len(compoundListMol)) if not compoundListMol[k] in ['FeH', 'NiCl', 'CuCl', 'CoH', 'CrO','ZnH', 'ZnO', 'ZnS', 'ZnCl']])))" ] } ], "metadata": { "kernelspec": { "display_name": "Python 2", "language": "python", "name": "python2" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 2 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", "version": "2.7.8" } }, "nbformat": 4, "nbformat_minor": 0 }