<> Navigation: * [[AJMPublic/camcasp|CamCASP]] * [[AJMGrpOnly|AJM main]] = Bugs = == Units in SAPT(DFT) output == The output correctly converts to mH and kcal/mol, but does not convert to kJ/mol. To fix this edit the file {{{ lib/psi4/driver/procrouting/sapt/sapt_util.py }}} to make the change to '''constants.hartree2kJmol''': [[attachment:sapt_util.py (fragment)]] {{{ def print_sapt_var(name, value, short=False, start_spacer=" "): """ Converts the incoming value as hartree to a correctly formatted Psi print format. """ vals = (name, value * 1000, value * constants.hartree2kcalmol, value * constants.hartree2kJmol) if short: return start_spacer + "%-20s % 15.8f [mEh]" % vals[:2] else: return start_spacer + "%-20s % 15.8f [mEh] % 15.8f [kcal/mol] % 15.8f [kJ/mol]" % vals }}} == Cannot use GGAs == The SAPT(DFT) module works correctly when the default, PBE0, functional is used. But when this is changed to a GGA like PBE it fails with the error message: {{{ => Auxiliary Basis Set <= Basis Set: ANONYMOUS400E6A7D Blend: AUG-CC-PVDZ-JKFIT + MB-SET-RI Number of shells: 478 Number of basis function: 1452 Number of Cartesian functions: 1697 Spherical Harmonics?: true Max angular momentum: 4 Traceback (most recent call last): File "/home/alston/Psi4/install/1.2-gcc/bin/psi4", line 259, in exec(content) File "", line 71, in File "/home/alston/Psi4/install/1.2-gcc/lib//psi4/driver/driver.py", line 460, in energy wfn = procedures['energy'][lowername](lowername, molecule=molecule, '''kwargs) File "/home/alston/Psi4/install/1.2-gcc/lib//psi4/driver/procrouting/sapt/sapt_proc.py", line 234, in run_sapt_dft cache = sapt_jk_terms.build_sapt_jk_cache(wfn_A, wfn_B, sapt_jk, True) File "/home/alston/Psi4/install/1.2-gcc/lib//psi4/driver/procrouting/sapt/sapt_jk_terms.py", line 108, in build_sapt_jk_cache cache["K_A"] = jk.K()[0].clone() IndexError: list index out of range '''* Psi4 encountered an error. Buy a developer more coffee! '''* Resources and help at github.com/psi4/psi4. }}} The PBE functional can be selected using: {{{ set { scf_type df sapt_dft_functional PBE sapt_dft_do_dHF False sapt_dft_grac_shift_a 0.12589 sapt_dft_grac_shift_b 0.12589 } }}} The problem seems to be in ''Psi4/install/1.2-gcc/lib/psi4/driver/procrouting/sapt/sapt_jk_terms.py'': [[attachment:sapt_jk_terms.py]] {{{ # J and K matrices jk.C_clear() # Normal J/K for Monomer A jk.C_left_add(wfn_A.Ca_subset("SO", "OCC")) jk.C_right_add(wfn_A.Ca_subset("SO", "OCC")) # Normal J/K for Monomer B jk.C_left_add(wfn_B.Ca_subset("SO", "OCC")) jk.C_right_add(wfn_B.Ca_subset("SO", "OCC")) # K_O J/K C_O_A = core.Matrix.triplet( cache["D_B"], cache["S"], cache["Cocc_A"], False, False, False) jk.C_left_add(C_O_A) jk.C_right_add(cache["Cocc_A"]) jk.compute() # Clone them as the JK object will overwrite. cache["J_A"] = jk.J()[0].clone() cache["K_A"] = jk.K()[0].clone() cache["J_B"] = jk.J()[1].clone() cache["K_B"] = jk.K()[1].clone() }}} If a GGA is used, then the exchange matrices are not created so the clone operations on jk.K() fail. {{{#!wiki important Important 15-03-18 This bug has been fixed by Daniel. Version 1.1-938-gfc51d9a (current master on github) works with GGAs. }}} == SAPT(DFT) not possible without Asymptotic correction == Curiously, SAPT(DFT) calculations are not possible without the asymptotic correction. In ''Psi4/install/1.2-gcc/lib/psi4/driver/procrouting/sapt/sapt_proc.py'' we have the lines: [[attachment:sapt_proc.py]] {{{ if (sapt_dft_functional != "HF") and ((mon_a_shift == 0.0) or (mon_b_shift == 0.0)): raise ValidationError('SAPT(DFT): must set both "SAPT_DFT_GRAC_SHIFT_A" and "B".') if (core.get_option('SCF', 'REFERENCE') != 'RHF'): raise ValidationError('SAPT(DFT) currently only supports restricted references.') }}} This needs to be changed as it should be possible to switch off the AC completely. = Testing Psi4 = == SAPT(DFT) == Currently, Psi4 can perform SAPT(DFT) calculations only with the ALDA kernel. This is the biggest limitation of the code. The following are the features of Psi4 1.2: * ALDA kernel * $E||{(2)}_{\rm ind,exch}$ is computed with the $S||2$ approximation. * GRAC asymptotic correction. * Range-separation also a possibility. * DC or DC+ basis sets. * Second-order terms only. Higher-order with the $\delta||{\rm HF}_{\rm int}$ approx. '''Questions to be answered:''' - What is the accuracy of SAPT(DFT) with the PBE0/GRAC functional with ALDA kernel? - Same for wPBE0(IP-tuned) with ALDA? - How fast is this code? What are the kinds of calculations we may expect to perform? - How well does the code parallelize? === Accuracy: PBE0+(GRAC/LB94) + ALDA === ==== Water dimers ==== Data from ''molecules/water/induction-damping/2017-ISA-ISApol-noS2/saptdft/aTZ-MC+/PBE0-GRLB/'' * Ref(1) : SAPT(DFT) DALTON+CamCASP: PBE0/GRLB : ALDA+CHF kernel : aTZ MC+BS (Mb=3s2p1d1f) * Psi4(1) : SAPT(DFT) PBE0/GRLB : ALDA kernel : aTZ DC+BS (Mb=3s2p1d1f). Here the $S||2$ approx is used to calculate the exch-ind energy. This is quite a poor approx for the 900 geom. Notice how the error is also present in the dHF term. * Psi4(2) : SAPT(DFT) PBE/GRLB : AC-shift = 0.202938 a.u. Other details as with Psi4(1). * '''UNITS''': kJ/mol * water2_900: $R_{O..O} = 2.2 $ Angstrom * water2_907: $R_{O..O} = 2.9924$ Angstrom || || water2_900 |||||| water2_907 |||||| || Energy || Ref(1) || Psi4(1) || Psi4(2) || Ref(1) || Psi4(1) || Psi4(2) || || $E||{(1)}_{\rm elst}$ || -160.995 || -160.459 || -162.198 || -28.934 || -28.814 || -28.317 || || $E||{(1)}_{\rm exch}$ || 378.326 || 376.891 || 392.954 || 25.311 || 24.687 || 26.413 || || $E||{(2)}_{\rm IND}$ || -44.758 || -71.218 || -73.479 || -4.474 || -4.606 || -4.654 || || $E||{(2)}_{\rm ind}$ || -217.284 || -215.744 || -230.843 || -11.054 || -10.782 || -11.348 || || $E||{(2)}_{\rm ind,exch}$ || 172.526 || 144.526 || 157.364 || 6.580 || 6.176 || 6.693 || || $E||{(2)}_{\rm DISP}$ || -44.617 || -36.794 || -38.893 || -8.521 || -7.409 || -8.324 || || $E||{(2)}_{\rm int}$ || 127.955 || 108.419 || 118.383 || -16.619 || -16.143 || -14.882 || || $\delta||{\rm HF}_{\rm int}$ || -56.708 || -35.685 || -35.685 || -2.931 || -2.864 || -2.864 || || $E||{(2)}_{\rm int}+\delta||{\rm HF}_{\rm int}$ || 71.246 || 72.734 || 82.698 || -19.550 || -19.008 || -17.747 || * The differences in $E|{(2)}_{\rm DISP}$ is mainly from the kernel. There is also a minor difference because of the basis (MC+BS vs DC+BS). The ALDA kernel results in too little dispersion. * The differences in $E|{(2)}_{\rm IND}$ is from the $S|2$ approx and the kernel. * The $S|2$ approx is the major source of error for the induction. * '''What is the origin of the differences in the first-order exchange energies?''' * The GRLB asymptotic correction in Psi4 uses the LB94+VWN at long-range. I would have thought they should have retained the PBEc functional entirely and corrected only the PBEx functional. What does DALTON do? '''Q: How does PBE/ALDA fare (Psi4(2))?''' * The second-order energies appear better, but any gains are offset by a too-repulsive first-order exchange energy. * So there is insufficient binding with PBE. * __Almost want to combine PBE0 first-order with PBE second-order...__ === Pyridine dimers === Data from: ''/home/alston/molecules/pyridine/scans/dimers_from_Model3_noS2/Min1-Hb1/aDZMb3-MC+/psi4'' * Ref(1) : PBE0/CS00 '''ALDA+CHF''' aDZ+Mb3 MC+BS * Ref(2) : PBE0/CS00 '''ALDA''' aDZ+Mb3 MC+BS * Psi4-PBE0 : PBE0/GRLB ''' ALDA''' aDZ+Mb3 DC+BS * Psi4-PBE : PBE/GRLB ALDA || || Hb1_5 aDZ+Mb3 |||||||| || Energy || Ref(1) || Ref(2) || Psi4-PBE0 || Psi4-PBE || || $E||{(1)}_{\rm elst}$ || -20.5598 || -20.5598 || -20.9766 || -20.8048 || || $E||{(1)}_{\rm exch}$ || 28.0175 || 28.0175 || 28.5536 || 31.0925 || || $E||{(2)}_{\rm IND}$ || -3.5901 || -3.1628 || -3.7196 || -3.6984 || || $E||{(2)}_{\rm ind}$ || -8.5829 || -7.8752 || -8.9282 || -9.5670 || || $E||{(2)}_{\rm ind,exch}$ || 4.9928 || 4.7123 || 5.2085 || 5.8686 || || $E||{(2)}_{\rm DISP}$ || -16.2333 || -15.3627 || -14.1152 || -15.6091 || || $E||{(2)}_{\rm int}$ || -12.3657 || -11.0679 || -10.2580 || -9.0198 || || $\delta||{\rm HF}_{\rm int}$ || -2.3259 || -2.3259 || -2.3422 || -2.3422 || || $E||{(2)}_{\rm int}+\delta||{\rm HF}_{\rm int}$ || -14.6917 || -13.3939 || -12.6002 || -11.3620 || Data from ''/home/alston/molecules/pyridine/scans/dimers_from_Model3_noS2/Min2-S1/aTZMb3-MC+/psi4'' * Psi4-PBE: Notice how large the first-order exchange is for Psi4-PBE. This seems to be the case in general. * Psi4-PBE: The second-order energies are closer to the reference values. * Psi4-hybrid: first-order from Psi4-PBE0 and second from Psi4-PBE: Eint(2) = -11.7573 kJ/mol which becomes -14.0995 kJ/mol with deltaHF. This is much better. But does this hold for other configurations? {{{ Similar data for the Hb1_5 dimer with the aTZ+Mb3 basis: Cam-PBE0 : CamCASP PBE0/CS00 ALDA+CHF MC+BS basis Cam-PBE : CamCASP PBE/CS00 ALDA MC+BS basis Psi-PBE0 : Psi4 PBE0/GRAC ALDA DC+BS Psi-PBE : Psi4 PBE/GRAC ALDA DC+BS Psi-Hyb : first-order from Psi-PBE0 and second-order from Psi-PBE Psi-Cam-FUNC-KER : Psi4 MOs passed to CamCASP with functional FUNC and kernel KER. Kernels are : ALDA or HYB=ALDA+CHF kJ/mol elst exch IND2 dHF DISP2 Eint CT2 POL2 Cam-PBE0 -20.45658 27.84655 -3.61371 -2.36352 -17.02206 -15.60932 -0.52554 -3.08817 Psi-PBE0 -20.83304 28.52921 -3.68395 -2.36648 -14.53605 -12.89031 Psi-Cam-PBE0-HYB -20.84369 28.56603 -3.61448 -2.36352 -17.22447 -15.48013 -0.50929 -3.10519 Psi-Cam-PBE0-ALDA -20.84369 28.56603 -3.22050 -2.36352 -15.89760 -13.75929 -0.50848 -2.71203 --------------------------------------------------------------------------------------- Psi-PBE0-scaled -20.83304 28.52921 -3.68395 -2.36648 -15.89755 -14.25181 Psi-Cam-PBE0-mix -20.84369 28.56603 -3.61448 -2.36352 -15.89760 -14.15326 --------------------------------------------------------------------------------------- --------------------------------------------------------------------------------------- Cam-PBE -20.29438 31.92628 -3.44943 -2.36352 -18.02445 -12.20549 -0.51709 -2.93234 Psi-PBE -20.69212 31.03235 -3.65223 -2.36648 -16.07874 -11.75723 Psi-PBE-scaled -20.69212 31.03235 -3.65223 -2.36648 -18.06334 -13.74182 Psi-Cam-PBE-ALDA -20.69365 31.03205 -3.56705 -2.36352 -17.92890 -13.52107 -0.54404 -3.02301 --------------------------------------------------------------------------------------- Psi-Hyb -20.83304 28.52921 -3.65223 -2.36648 -16.07874 -14.40128 --------------------------------------------------------------------------------------- CCSD(T) augA-Sadlej+Mb3 DC+BS Psi4: -15.84991 aug-cc-pVDZ+Mb3 DC+BS Psi4: -15.16522 --------------------------------------------------------------------------------------- }}} Points to note: * Cam-PBE0 is very close to CCSD(T). The basis set differences are inconsequential here. I have data for Cam-PBE0 with the augA-Sadlej basis too. * Both Psi4 results are off by a lot. * I would have expected Cam-PBE to be similar to Psi-PBE, but there are substantial differences in the E(1)exch and E(2)DISP. This could be because the CS00 and GRAC schemes are substantially different at this dimer geometry. Both use the LB94 long-range functional, but differ in how this is spliced to the PBE functional. '''TO DO: use the Psi4 orbitals in CamCASP. We should get the same energies!''' '''How does Psi4 calculate the second-order energies?''' * Induction: CHF+ALDA kernel. This is possible as it uses the iterated solver for the coupled equations. S2 approx used for exch-ind. * Dispersion: ALDA kernel. Exch-disp is calculated at un-coupled level. No scaling is applied. We have long recommended that scaling should be applied: $$E||{(2)}_{\rm disp,exch}[C] \approx E||{(2)}_{\rm disp,exch}[UC] \times \frac{E||{(2)}_{\rm disp}[C]}{E||{(2)}_{\rm disp}[UC]}$$ * The '''Psi-FUNC-scaled''' results above use this expression and result in much better total interaction energies. {{{#!wiki important Important '''Are CamCASP and Psi4 consistent?''' I have now got CamCASP to work with Psi4 1.2 MOs. There were changes needed with the interface and a few other issues. From the above discussion, we see that we should compare '''Psi-PBE0-scaled''' energies to CamCASP calculations computed using the first-order energies and the induction and delta-HF energies taken from '''Psi-Cam-PBE0-HYB''' and the second-order dispersion energy from '''Psi-Cam-PBE0-ALDA'''. This resulting theory is termed '''Psi-Cam-PBE0-mix''' above. And indeed, we see good agreement between the results from the two codes. There are differences of a tenth of a kJ/mol, but these are due to density-fitting differences and the S||2 approx made by Psi4. }}} {{{#!wiki important Important '''Recommendation:''' If we must use Psi4 only, then we should compute the interaction energy using PBE0 (as first-order is good) with scaled exchange-dispersion as shown above. * This needs to be tested on different dimers. * The error in the dispersion is still present, but is smaller. }}}