FYI: difference in chemical shielding calculations between gpu4pyscf and pyscf

[steto123]

Hello,

for your information: i found some small differences between the isotropic shieldings calculated with gpu4pyscf and pyscf (i know the experimental character of this calculations). I have shared a pdf with detailed information because my former issue was not saved here.

Is there a simple reason for this differences. Is there a mistake or error in my code?

https://acrobat.adobe.com/id/urn:aaid:sc:EU:96d4d829-8366-4887-aef5-3e6c411dc82e

Many thanks for your great work.

Steffen

[wxj6000]

@steto123 Hi Steffen, can you share the xyz file? I didn't find the large discrepancy. The difference should be less than 1e-5 ppm. And I added some checks of consistency between PySCF and GPU4PySCF for a simple molecule.

[steto123]

Attached is the xyz file for acetonitrile. The structure comes from the test data of the CSESHIRE data collection.
http://cheshirenmr.info/CalculationFiles.htm
I have calculated all connections for the calibration used there with Orca5 and am in the process of doing the same for pyscf. The deviation in the calculated shielding between orca5 and pyscf lies in the 2nd to 4th decimal place, with deviations for C being larger than for H. At least that has been my observation so far. However, the differences between the values of pyscf and gpu4pyscf are larger, which I have not yet considered and tabulated systematically. I have therefore assumed that it may be due to the way I implemented the calculation based on examples from github. I can also provide you with the jupyter notebooks I used. However, these are commented in German by me.
Feel free to contact me directly as well. I can be reached at the anti-track email [email protected]. I delete such emails from time to time. However, they are forwarded to my actual email, so publishing them here is not a problem.

Acetonitrile.xyz.txt

For security reasons i must rename this to .txt

[wxj6000]

@steto123 Sorry for the delay. The inconsistency between PySCF and GPU4PySCF is due to the inconsistent setup. In your code, mf is generated by mol.RKS().to_gpu().run() which by default is using LDA functional. Here is the code I modified.

import pyscf
import gpu4pyscf
from gpu4pyscf import dft
from gpu4pyscf.properties.shielding import eval_shielding
mol = pyscf.M(atom=xyzin, basis='6-31g*', output=logfile)
mf = dft.RKS(mol)
mf.xc = 'b3lyp'
mf.kernel()
tensor = eval_shielding(mf)
abschirmung = tensor[0].get()+tensor[1].get()
for i in range(mol.natm):
    shift = (abschirmung[i,0,0] + abschirmung[i,1,1] + abschirmung[i,2,2])/3 
    print(shift)

[steto123]

Hello and many thanks for your help and explanations.
Now i have really similar results for molecules wich have only carbon and hydrogen. In the case of molecules with oxygen or nitrogen (like acetone, furan or THF) there are errors like:

---

IndexError Traceback (most recent call last)
Cell In[5], line 36
34 mf.xc = 'b3lyp'
35 mf.kernel()
---> 36 tensor = eval_shielding(mf)
39 abschirmung = tensor[0].get()+tensor[1].get()
40 end_time = time.time()

File ~/miniconda3/lib/python3.12/site-packages/gpu4pyscf/properties/shielding.py:232, in eval_shielding(mf)
229 veff_ai = contract('xav,vi->xai', tmp, mocc)
230 veff_ai -= contract('xai,i->xai', s1ai, mo_energy[idx_occ])
--> 232 mo1 = cphf.solve(fx, mo_energy, mo_occ, veff_ai, max_cycle=20, tol=1e-10)[0]
234 shielding_d = cupy.empty((natom, 3, 3))
235 shielding_p = cupy.empty((natom, 3, 3))

File ~/miniconda3/lib/python3.12/site-packages/gpu4pyscf/scf/cphf.py:43, in solve(fvind, mo_energy, mo_occ, h1, s1, max_cycle, tol, hermi, verbose, level_shift)
33 '''
34 Args:
35 fvind : function
(...)
39 Whether the matrix defined by fvind is Hermitian or not.
40 '''
42 if s1 is None:
---> 43 return solve_nos1(fvind, mo_energy, mo_occ, h1,
44 max_cycle, tol, hermi, verbose)
45 else:
46 return solve_withs1(fvind, mo_energy, mo_occ, h1, s1,
47 max_cycle, tol, hermi, verbose)

File ~/miniconda3/lib/python3.12/site-packages/gpu4pyscf/scf/cphf.py:70, in solve_nos1(fvind, mo_energy, mo_occ, h1, max_cycle, tol, hermi, verbose, level_shift)
68 v = e_ai
69 return v.reshape(-1,nvirnocc)
---> 70 mo1 = krylov(vind_vo, mo1base.reshape(-1,nvir*nocc),
71 tol=tol, max_cycle=max_cycle, hermi=hermi, verbose=log)
72 log.timer('krylov solver in CPHF', *t0)
73 return mo1.reshape(h1.shape), None

File ~/miniconda3/lib/python3.12/site-packages/gpu4pyscf/lib/cupy_helper.py:530, in krylov(aop, b, x0, tol, max_cycle, dot, lindep, callback, hermi, verbose)
528 max_cycle = min(max_cycle, ndim)
529 for cycle in range(max_cycle):
--> 530 axt = aop(x1)
531 if axt.ndim == 1:
532 axt = axt.reshape(1,ndim)

File ~/miniconda3/lib/python3.12/site-packages/gpu4pyscf/scf/cphf.py:65, in solve_nos1..vind_vo(mo1)
64 def vind_vo(mo1):
---> 65 v = fvind(mo1.reshape(-1,nvir,nocc)).reshape(-1,nvir,nocc)
66 if level_shift != 0:
67 v -= mo1 * level_shift

File ~/miniconda3/lib/python3.12/site-packages/gpu4pyscf/properties/shielding.py:55, in gen_vind..fx(mo1)
53 v1 = np.empty((3, nao, nao))
54 for i in range(3):
---> 55 v1[i] = -jk.get_jk(mf.mol, dm1[i].get(), 'ijkl,jk->il')0.5hyb
56 v1 = cupy.array(v1)
57 tmp = contract('nuv,vi->nui', v1, mocc)

File cupy/_core/core.pyx:1527, in cupy._core.core._ndarray_base.getitem()

File cupy/_core/_routines_indexing.pyx:33, in cupy._core._routines_indexing._ndarray_getitem()

File cupy/_core/_routines_indexing.pyx:410, in cupy._core._routines_indexing._view_getitem()

IndexError: Index 1 is out of bounds for axis 0 with size 1

[wxj6000]

@steto123 Can you upload your script and xyz file? I am not able to reproduce the issue with the above code and molecules with oxygen or nitrogen.

[steto123]

Hello,
i have added the extension .txt to this files. The latest test is from now. I have also attached cyclobutane.xyz wich works fine.
i create the .xyz files using coordinates from the csheshire dataset and an script with import functionality from rdkit.

nmr-shift2g.ipynb.txt
Acetone.xyz.txt
Acetone_gpu.log.txt
Cyclobutane.xyz.txt

[steto123]

Today i have checked this on a other machine with older CUDA driver (11.6 i think) with the same result. Running well for cyclobutane and crash for acetone.

[wxj6000]

@steto123 Thank you for providing the molecules and scripts. I am able to reproduce the issue. It has been fixed in #286

[steto123]

I have updated my installation with your fixes manually and now it works. Many Thanks!