problem with i<k|(r-Rc)Xp|n> while interfacing to stda program #50
Comments
|
Dear Marc,
It might be due to the ordering of Cartesian basis. The ordering in libcint code
follows
for nx in l ... 0:
for ny in (l-nx) .. 0:
nz = l - nx - ny
nx * 'x' + ny * 'y' + nz * 'z'
I'm not quite sure the comparison in your example due to the ordering of basis
functions. Here are the numbers for d-type basis in your first system
```
array([[[ 0. , 0. , 0. , 0. , -0. , 0. ],
[ 0. , 0. , 1. , 0. , 0. , 0. ],
[ 0. , -1. , 0. , 0. , 0. , 0. ],
[ 0. , 0. , 0. , 0. , 1.155, 0. ],
[-0. , 0. , 0. , -1.155, 0. , 1.155],
[ 0. , 0. , 0. , 0. , -1.155, 0. ]],
[[ 0. , 0. , -1.155, 0. , 0. , 0. ],
[ 0. , 0. , 0. , 0. , -1. , 0. ],
[ 1.155, 0. , 0. , 0. , 0. , -1.155],
[ 0. , 0. , 0. , 0. , 0. , 0. ],
[ 0. , 1. , 0. , 0. , 0. , 0. ],
[ 0. , 0. , 1.155, 0. , 0. , 0. ]],
[[ 0. , 1.155, 0. , 0. , 0. , 0. ],
[-1.155, 0. , 0. , 1.155, 0. , 0. ],
[ 0. , 0. , 0. , 0. , 1. , 0. ],
[ 0. , -1.155, 0. , 0. , 0. , 0. ],
[ 0. , 0. , -1. , 0. , 0. , 0. ],
[ 0. , 0. , 0. , 0. , 0. , 0. ]]])
```
The three blocks corresponds to three components. In each block, the basis are ordered as
```
['0 H 3dxx ',
'0 H 3dxy ',
'0 H 3dxz ',
'0 H 3dyy ',
'0 H 3dyz ',
'0 H 3dzz ']
```
If you have pyscf installed, you can reproduce the results using the script below
```
mol = pyscf.M(atom='H', basis=[[2, [2.062, 1]]], spin=None, cart=True)
norm = mol.intor('int1e_ovlp').diagonal()**-.5
print(np.einsum('xij,i,j->xij', mol.intor('int1e_cg_irxp'), norm, norm).round(3))
print(mol.ao_labels())
```
…On Tue, Jan 12, 2021 at 02:34:25AM -0800, Marc de Wergifosse wrote:
I am trying to modernize the [stda](https://github.com/grimme-lab/stda) program by removing its old integral deck, replacing it by libcint, also allowing new types of integrals that were not available. This is working perfectly fine for the overlap, dipole moment, and <i|nabla|j> integrals for which I can reproduce values obtained by our integral deck (see https://github.com/grimme-lab/stda/blob/master/intslvm.f and https://github.com/grimme-lab/stda/blob/master/intpack.f90).
But I have some issues with i<k|(r-Rc)Xp|n> type of integrals while using exactly the same interface, just replacing the name of the subroutine by cint1e_cg_irxp_cart. Note that everything is done with cartesian basis functions and that I have well defined Rc in env.
For integrals on the same atom, few elements have the reverse sign while others are ok for this skew-symmetric matrix.
For <d| |d> integrals, <dyy| |dxy>, <dzz| |dxz>, and <dzz| |dyz> signs are swapped with their counterpart.
For <f| |f> integrals, I have the same issue with <fzzz| |fxxy>, <fyyy| |fxxz>, <fyyy| |fyyx>, <fzzz| |fyyx>, <fyyy| |fxzz>, <fzzz| |fxzz>, <fzzz| |fyzz>, <fyyz| |fxzz>, <fyyz| |fxyz>, <fxzz| |fxyz>, and <fyzz| |fxyz>. All others are ok.
In addition to this, for integrals between two different atoms, I obtained different values than with our integral deck for <s| |fxxx>, <s| |fyyy>, <s| |fzzz>, <p| |fxxx>, <p| |fyyy>, <p| |fzzz>, <dxx| |dxx>, <dyy| |dyy>, <dzz| |dzz>, and almost all <f| |f>. Note that for <dxx| |dxx>, <dyy| |dyy>, and <dzz| |dzz> the sign is also wrong.
Those are what I was able to identify with small models.
Here are some examples:
1° case
For H- with one set of d orbitals:
```
[Atoms] Angs
H 1 1 0.00000000 0.00000000 0.00000000
[GTO]
1 0
d 1 1.000000
2.06200000E+00 1.00000000E+00
```
```
-------------- 1 --------------magnetic_moment
old 0.00000 0.00000 0.00000
new 0.00000 0.00000 0.00000
-------------- 2 --------------magnetic_moment
old 0.00000 0.00000 0.00000
new 0.00000 0.00000 0.00000
-------------- 3 --------------magnetic_moment
old 0.00000 0.00000 0.00000
new 0.00000 0.00000 0.00000
-------------- 4 --------------magnetic_moment
old 0.00000 0.00000 0.00000
new 0.00000 0.00000 0.00000
-------------- 5 --------------magnetic_moment
old 0.00000 0.00000 0.00000
new 0.00000 0.00000 0.00000
-------------- 6 --------------magnetic_moment
old 0.00000 0.00000 0.00000
new 0.00000 0.00000 0.00000
-------------- 7 --------------magnetic_moment
old 0.00000 0.00000 -1.15470
new 0.00000 0.00000 -1.15470
-------------- 8 --------------magnetic_moment
old 0.00000 0.00000 1.15470
new 0.00000 0.00000 -1.15470
**********************
-------------- 9 --------------magnetic_moment
old 0.00000 0.00000 0.00000
new 0.00000 0.00000 0.00000
-------------- 10 --------------magnetic_moment
old 0.00000 0.00000 0.00000
new 0.00000 0.00000 0.00000
-------------- 11 --------------magnetic_moment
old 0.00000 1.15470 0.00000
new 0.00000 1.15470 0.00000
-------------- 12 --------------magnetic_moment
old 0.00000 0.00000 0.00000
new 0.00000 0.00000 0.00000
-------------- 13 --------------magnetic_moment
old 0.00000 -1.15470 0.00000
new 0.00000 1.15470 0.00000
**********************
-------------- 14 --------------magnetic_moment
old -1.00000 0.00000 0.00000
new -1.00000 0.00000 0.00000
-------------- 15 --------------magnetic_moment
old 0.00000 0.00000 0.00000
new 0.00000 0.00000 0.00000
-------------- 16 --------------magnetic_moment
old 0.00000 0.00000 0.00000
new 0.00000 0.00000 0.00000
-------------- 17 --------------magnetic_moment
old -1.15470 0.00000 0.00000
new -1.15470 0.00000 0.00000
-------------- 18 --------------magnetic_moment
old 1.15470 0.00000 0.00000
new -1.15470 0.00000 0.00000
**********************
-------------- 19 --------------magnetic_moment
old 0.00000 1.00000 0.00000
new 0.00000 1.00000 0.00000
-------------- 20 --------------magnetic_moment
old 0.00000 0.00000 -1.00000
new 0.00000 0.00000 -1.00000
-------------- 21 --------------magnetic_moment
old 0.00000 0.00000 0.00000
new 0.00000 0.00000 0.00000
```
These are the upper triangle of the skew-symmetric matrix where d orbitals are ordered as dxx dyy dzz dxy dxz dyz. Thus, 8, 13, and 18 are <dyy| |dxy>, <dzz| |dxz>, and <dzz| |dyz>, respectively.
2° case
Same system with one set of f functions.
```
[Atoms] Angs
H 1 1 0.00000000 0.00000000 0.00000000
[GTO]
1 0
f 1 1.000000
1.39700000E+00 1.00000000E+00
```
```
-------------- 1 --------------magnetic_moment
old 0.00000 0.00000 0.00000
new 0.00000 0.00000 0.00000
-------------- 2 --------------magnetic_moment
old 0.00000 0.00000 -0.60000
new 0.00000 0.00000 -0.60000
-------------- 3 --------------magnetic_moment
old 0.00000 0.00000 0.00000
new 0.00000 0.00000 0.00000
-------------- 4 --------------magnetic_moment
old 0.00000 0.60000 0.00000
new 0.00000 0.60000 0.00000
-------------- 5 --------------magnetic_moment
old -0.60000 0.00000 0.00000
new -0.60000 0.00000 0.00000
-------------- 6 --------------magnetic_moment
old 0.00000 0.00000 0.00000
new 0.00000 0.00000 0.00000
-------------- 7 --------------magnetic_moment
old 0.00000 0.00000 -1.34164
new 0.00000 0.00000 -1.34164
-------------- 8 --------------magnetic_moment
old 0.00000 0.00000 0.00000
new 0.00000 0.00000 0.00000
-------------- 9 --------------magnetic_moment
old 0.44721 0.00000 0.00000
new -0.44721 0.00000 0.00000
**********************
-------------- 10 --------------magnetic_moment
old 0.00000 0.00000 0.00000
new 0.00000 0.00000 0.00000
-------------- 11 --------------magnetic_moment
old 0.00000 1.34164 0.00000
new 0.00000 1.34164 0.00000
-------------- 12 --------------magnetic_moment
old -0.44721 0.00000 0.00000
new 0.44721 0.00000 0.00000
**********************
-------------- 13 --------------magnetic_moment
old 0.00000 0.00000 0.00000
new 0.00000 0.00000 0.00000
-------------- 14 --------------magnetic_moment
old -1.00000 0.00000 0.00000
new -1.00000 0.00000 0.00000
-------------- 15 --------------magnetic_moment
old 0.00000 0.00000 0.00000
new 0.00000 0.00000 0.00000
-------------- 16 --------------magnetic_moment
old 0.00000 0.00000 0.00000
new 0.00000 0.00000 0.00000
-------------- 17 --------------magnetic_moment
old 0.00000 0.00000 1.34164
new 0.00000 0.00000 -1.34164
**********************
-------------- 18 --------------magnetic_moment
old 0.00000 -0.44721 0.00000
new 0.00000 0.44721 0.00000
**********************
-------------- 19 --------------magnetic_moment
old 0.00000 0.00000 -1.00000
new 0.00000 0.00000 -1.00000
-------------- 20 --------------magnetic_moment
old 0.00000 -0.33333 0.00000
new 0.00000 -0.33333 0.00000
-------------- 21 --------------magnetic_moment
old 0.00000 0.00000 0.00000
new 0.00000 0.00000 0.00000
-------------- 22 --------------magnetic_moment
old 0.00000 0.44721 0.00000
new 0.00000 0.44721 0.00000
-------------- 23 --------------magnetic_moment
old -1.34164 0.00000 0.00000
new -1.34164 0.00000 0.00000
-------------- 24 --------------magnetic_moment
old 0.00000 0.00000 0.00000
new 0.00000 0.00000 0.00000
-------------- 25 --------------magnetic_moment
old -0.33333 0.00000 0.00000
new -0.33333 0.00000 0.00000
-------------- 26 --------------magnetic_moment
old 0.00000 0.00000 0.00000
new 0.00000 0.00000 0.00000
-------------- 27 --------------magnetic_moment
old 0.00000 1.00000 0.00000
new 0.00000 1.00000 0.00000
-------------- 28 --------------magnetic_moment
old 0.00000 0.00000 0.00000
new 0.00000 0.00000 0.00000
-------------- 29 --------------magnetic_moment
old 0.00000 0.00000 0.00000
new 0.00000 0.00000 0.00000
-------------- 30 --------------magnetic_moment
old 0.00000 0.00000 0.44721
new 0.00000 0.00000 -0.44721
**********************
-------------- 31 --------------magnetic_moment
old 0.00000 -1.34164 0.00000
new 0.00000 1.34164 0.00000
**********************
-------------- 32 --------------magnetic_moment
old 0.00000 0.00000 0.33333
new 0.00000 0.00000 0.33333
-------------- 33 --------------magnetic_moment
old 0.00000 1.00000 0.00000
new 0.00000 1.00000 0.00000
-------------- 34 --------------magnetic_moment
old 0.00000 0.00000 0.00000
new 0.00000 0.00000 0.00000
-------------- 35 --------------magnetic_moment
old 0.00000 -0.33333 0.00000
new 0.00000 0.33333 0.00000
**********************
-------------- 36 --------------magnetic_moment
old 0.00000 0.00000 0.00000
new 0.00000 0.00000 0.00000
-------------- 37 --------------magnetic_moment
old 0.00000 0.00000 -0.44721
new 0.00000 0.00000 -0.44721
-------------- 38 --------------magnetic_moment
old 0.00000 0.00000 0.00000
new 0.00000 0.00000 0.00000
-------------- 39 --------------magnetic_moment
old 1.34164 0.00000 0.00000
new -1.34164 0.00000 0.00000
**********************
-------------- 40 --------------magnetic_moment
old 0.00000 0.00000 0.00000
new 0.00000 0.00000 0.00000
-------------- 41 --------------magnetic_moment
old 0.33333 0.00000 0.00000
new 0.33333 0.00000 0.00000
-------------- 42 --------------magnetic_moment
old 0.00000 0.00000 -0.33333
new 0.00000 0.00000 -0.33333
-------------- 43 --------------magnetic_moment
old -1.00000 0.00000 0.00000
new -1.00000 0.00000 0.00000
-------------- 44 --------------magnetic_moment
old 0.00000 0.00000 -1.00000
new 0.00000 0.00000 -1.00000
-------------- 45 --------------magnetic_moment
old 0.00000 0.00000 0.00000
new 0.00000 0.00000 0.00000
-------------- 46 --------------magnetic_moment
old 0.00000 0.00000 0.00000
new 0.00000 0.00000 0.00000
-------------- 47 --------------magnetic_moment
old 0.00000 0.00000 0.00000
new 0.00000 0.00000 0.00000
-------------- 48 --------------magnetic_moment
old 0.00000 0.00000 0.00000
new 0.00000 0.00000 0.00000
-------------- 49 --------------magnetic_moment
old 0.00000 1.15470 0.00000
new 0.00000 1.15470 0.00000
-------------- 50 --------------magnetic_moment
old 0.00000 0.00000 -1.15470
new 0.00000 0.00000 -1.15470
-------------- 51 --------------magnetic_moment
old -1.15470 0.00000 0.00000
new -1.15470 0.00000 0.00000
-------------- 52 --------------magnetic_moment
old 0.00000 0.00000 1.15470
new 0.00000 0.00000 -1.15470
**********************
-------------- 53 --------------magnetic_moment
old 1.15470 0.00000 0.00000
new -1.15470 0.00000 0.00000
**********************
-------------- 54 --------------magnetic_moment
old 0.00000 -1.15470 0.00000
new 0.00000 1.15470 0.00000
**********************
-------------- 55 --------------magnetic_moment
old 0.00000 0.00000 0.00000
new 0.00000 0.00000 0.00000
```
The f functions are ordered as fxxx, fyyy, fzzz, fxxy, fxxz, fyyx, fyyz, fxzz, fyzz, and fxyz. As for the case above, the upper diagonal is printed and problematic elements are <fzzz| |fxxy>, <fyyy| |fxxz>, <fyyy| |fyyx>, <fzzz| |fyyx>, <fyyy| |fxzz>, <fzzz| |fxzz>, <fzzz| |fyzz>, <fyyz| |fxzz>, <fyyz| |fxyz>, <fxzz| |fxyz>, and <fyzz| |fxyz>.
3° case
For CH4 with s functions on H and one set of f functions for the carbon.
```
[Atoms] Angs
C 1 6 -0.00000000 0.00000000 0.00000000
H 2 1 -0.61776480 0.61776480 0.61776480
H 3 1 0.61776480 -0.61776480 0.61776480
H 4 1 0.61776480 0.61776480 -0.61776480
H 5 1 -0.61776480 -0.61776480 -0.61776480
[GTO]
1 0
f 1 1.000000
1.41900000E+00 1.00000000E+00
2 0
s 1 1.000000
7.97700000E-01 1.00000000E+00
3 0
s 1 1.000000
7.97700000E-01 1.00000000E+00
4 0
s 1 1.000000
7.97700000E-01 1.00000000E+00
5 0
s 1 1.000000
7.97700000E-01 1.00000000E+00
```
```
-------------- 1 --------------magnetic_moment
old 0.00000 0.00000 0.00000
new 0.00000 0.00000 0.00000
-------------- 2 --------------magnetic_moment
old 0.00000 0.00000 -0.60000
new 0.00000 0.00000 -0.60000
-------------- 3 --------------magnetic_moment
old 0.00000 0.00000 0.00000
new 0.00000 0.00000 0.00000
-------------- 4 --------------magnetic_moment
old 0.00000 0.60000 0.00000
new 0.00000 0.60000 0.00000
-------------- 5 --------------magnetic_moment
old -0.60000 0.00000 0.00000
new -0.60000 0.00000 0.00000
-------------- 6 --------------magnetic_moment
old 0.00000 0.00000 0.00000
new 0.00000 0.00000 0.00000
-------------- 7 --------------magnetic_moment
old 0.00000 0.00000 -1.34164
new 0.00000 0.00000 -1.34164
-------------- 8 --------------magnetic_moment
old 0.00000 0.00000 0.00000
new 0.00000 0.00000 0.00000
-------------- 9 --------------magnetic_moment
old 0.44721 0.00000 0.00000
new -0.44721 0.00000 0.00000
**********************
-------------- 10 --------------magnetic_moment
old 0.00000 0.00000 0.00000
new 0.00000 0.00000 0.00000
-------------- 11 --------------magnetic_moment
old 0.00000 1.34164 0.00000
new 0.00000 1.34164 0.00000
-------------- 12 --------------magnetic_moment
old -0.44721 0.00000 0.00000
new 0.44721 0.00000 0.00000
**********************
-------------- 13 --------------magnetic_moment
old 0.00000 0.00000 0.00000
new 0.00000 0.00000 0.00000
-------------- 14 --------------magnetic_moment
old -1.00000 0.00000 0.00000
new -1.00000 0.00000 0.00000
-------------- 15 --------------magnetic_moment
old 0.00000 0.00000 0.00000
new 0.00000 0.00000 0.00000
-------------- 16 --------------magnetic_moment
old 0.00000 0.00000 0.00000
new 0.00000 0.00000 0.00000
-------------- 17 --------------magnetic_moment
old 0.00000 0.00000 1.34164
new 0.00000 0.00000 -1.34164
**********************
-------------- 18 --------------magnetic_moment
old 0.00000 -0.44721 0.00000
new 0.00000 0.44721 0.00000
**********************
-------------- 19 --------------magnetic_moment
old 0.00000 0.00000 -1.00000
new 0.00000 0.00000 -1.00000
-------------- 20 --------------magnetic_moment
old 0.00000 -0.33333 0.00000
new 0.00000 -0.33333 0.00000
-------------- 21 --------------magnetic_moment
old 0.00000 0.00000 0.00000
new 0.00000 0.00000 0.00000
-------------- 22 --------------magnetic_moment
old 0.00000 0.44721 0.00000
new 0.00000 0.44721 0.00000
-------------- 23 --------------magnetic_moment
old -1.34164 0.00000 0.00000
new -1.34164 0.00000 0.00000
-------------- 24 --------------magnetic_moment
old 0.00000 0.00000 0.00000
new 0.00000 0.00000 0.00000
-------------- 25 --------------magnetic_moment
old -0.33333 0.00000 0.00000
new -0.33333 0.00000 0.00000
-------------- 26 --------------magnetic_moment
old 0.00000 0.00000 0.00000
new 0.00000 0.00000 0.00000
-------------- 27 --------------magnetic_moment
old 0.00000 1.00000 0.00000
new 0.00000 1.00000 0.00000
-------------- 28 --------------magnetic_moment
old 0.00000 0.00000 0.00000
new 0.00000 0.00000 0.00000
-------------- 29 --------------magnetic_moment
old 0.00000 0.00000 0.00000
new 0.00000 0.00000 0.00000
-------------- 30 --------------magnetic_moment
old 0.00000 0.00000 0.44721
new 0.00000 0.00000 -0.44721
**********************
-------------- 31 --------------magnetic_moment
old 0.00000 -1.34164 0.00000
new 0.00000 1.34164 0.00000
**********************
-------------- 32 --------------magnetic_moment
old 0.00000 0.00000 0.33333
new 0.00000 0.00000 0.33333
-------------- 33 --------------magnetic_moment
old 0.00000 1.00000 0.00000
new 0.00000 1.00000 0.00000
-------------- 34 --------------magnetic_moment
old 0.00000 0.00000 0.00000
new 0.00000 0.00000 0.00000
-------------- 35 --------------magnetic_moment
old 0.00000 -0.33333 0.00000
new 0.00000 0.33333 0.00000
**********************
-------------- 36 --------------magnetic_moment
old 0.00000 0.00000 0.00000
new 0.00000 0.00000 0.00000
-------------- 37 --------------magnetic_moment
old 0.00000 0.00000 -0.44721
new 0.00000 0.00000 -0.44721
-------------- 38 --------------magnetic_moment
old 0.00000 0.00000 0.00000
new 0.00000 0.00000 0.00000
-------------- 39 --------------magnetic_moment
old 1.34164 0.00000 0.00000
new -1.34164 0.00000 0.00000
**********************
-------------- 40 --------------magnetic_moment
old 0.00000 0.00000 0.00000
new 0.00000 0.00000 0.00000
-------------- 41 --------------magnetic_moment
old 0.33333 0.00000 0.00000
new 0.33333 0.00000 0.00000
-------------- 42 --------------magnetic_moment
old 0.00000 0.00000 -0.33333
new 0.00000 0.00000 -0.33333
-------------- 43 --------------magnetic_moment
old -1.00000 0.00000 0.00000
new -1.00000 0.00000 0.00000
-------------- 44 --------------magnetic_moment
old 0.00000 0.00000 -1.00000
new 0.00000 0.00000 -1.00000
-------------- 45 --------------magnetic_moment
old 0.00000 0.00000 0.00000
new 0.00000 0.00000 0.00000
-------------- 46 --------------magnetic_moment
old 0.00000 0.00000 0.00000
new 0.00000 0.00000 0.00000
-------------- 47 --------------magnetic_moment
old 0.00000 0.00000 0.00000
new 0.00000 0.00000 0.00000
-------------- 48 --------------magnetic_moment
old 0.00000 0.00000 0.00000
new 0.00000 0.00000 0.00000
-------------- 49 --------------magnetic_moment
old 0.00000 1.15470 0.00000
new 0.00000 1.15470 0.00000
-------------- 50 --------------magnetic_moment
old 0.00000 0.00000 -1.15470
new 0.00000 0.00000 -1.15470
-------------- 51 --------------magnetic_moment
old -1.15470 0.00000 0.00000
new -1.15470 0.00000 0.00000
-------------- 52 --------------magnetic_moment
old 0.00000 0.00000 1.15470
new 0.00000 0.00000 -1.15470
**********************
-------------- 53 --------------magnetic_moment
old 1.15470 0.00000 0.00000
new -1.15470 0.00000 0.00000
**********************
-------------- 54 --------------magnetic_moment
old 0.00000 -1.15470 0.00000
new 0.00000 1.15470 0.00000
**********************
-------------- 55 --------------magnetic_moment
old 0.00000 0.00000 0.00000
new 0.00000 0.00000 0.00000
**-------------- 56 --------------magnetic_moment
old 0.00000 0.21382 -0.21382
new 0.00000 0.20626 -0.20626
-------------- 57 --------------magnetic_moment
old -0.21382 0.00000 -0.21382
new -0.20626 -0.00000 -0.20626
-------------- 58 --------------magnetic_moment
old 0.21382 0.21382 0.00000
new 0.20626 0.20626 0.00000**
-------------- 59 --------------magnetic_moment
old -0.15374 -0.13497 -0.01877
new -0.15374 -0.13497 -0.01877
-------------- 60 --------------magnetic_moment
old 0.15374 0.01877 0.13497
new 0.15374 0.01877 0.13497
-------------- 61 --------------magnetic_moment
old 0.13497 0.15374 -0.01877
new 0.13497 0.15374 -0.01877
-------------- 62 --------------magnetic_moment
old 0.01877 0.15374 -0.13497
new 0.01877 0.15374 -0.13497
-------------- 63 --------------magnetic_moment
old -0.13497 0.01877 -0.15374
new -0.13497 0.01877 -0.15374
-------------- 64 --------------magnetic_moment
old -0.01877 0.13497 -0.15374
new -0.01877 0.13497 -0.15374
-------------- 65 --------------magnetic_moment
old 0.00000 0.00000 0.00000
new -0.00000 -0.00000 0.00000
-------------- 66 --------------magnetic_moment
old 0.00000 0.00000 0.00000
new 0.00000 0.00000 0.00000
**-------------- 67 --------------magnetic_moment
old 0.00000 0.21382 0.21382
new 0.00000 0.20626 0.20626
-------------- 68 --------------magnetic_moment
old -0.21382 0.00000 0.21382
new -0.20626 -0.00000 0.20626
-------------- 69 --------------magnetic_moment
old -0.21382 -0.21382 0.00000
new -0.20626 -0.20626 0.00000**
-------------- 70 --------------magnetic_moment
old -0.15374 -0.13497 0.01877
new -0.15374 -0.13497 0.01877
-------------- 71 --------------magnetic_moment
old -0.15374 -0.01877 0.13497
new -0.15374 -0.01877 0.13497
-------------- 72 --------------magnetic_moment
old 0.13497 0.15374 0.01877
new 0.13497 0.15374 0.01877
-------------- 73 --------------magnetic_moment
old -0.01877 -0.15374 -0.13497
new -0.01877 -0.15374 -0.13497
-------------- 74 --------------magnetic_moment
old -0.13497 0.01877 0.15374
new -0.13497 0.01877 0.15374
-------------- 75 --------------magnetic_moment
old -0.01877 0.13497 0.15374
new -0.01877 0.13497 0.15374
-------------- 76 --------------magnetic_moment
old 0.00000 0.00000 0.00000
new 0.00000 0.00000 0.00000
-------------- 77 --------------magnetic_moment
old -0.02810 -0.02810 0.00000
new -0.02810 -0.02810 0.00000
-------------- 78 --------------magnetic_moment
old 0.00000 0.00000 0.00000
new 0.00000 0.00000 0.00000
**-------------- 79 --------------magnetic_moment
old 0.00000 -0.21382 -0.21382
new 0.00000 -0.20626 -0.20626
-------------- 80 --------------magnetic_moment
old 0.21382 0.00000 0.21382
new 0.20626 -0.00000 0.20626
-------------- 81 --------------magnetic_moment
old 0.21382 -0.21382 0.00000
new 0.20626 -0.20626 0.00000**
-------------- 82 --------------magnetic_moment
old 0.15374 -0.13497 0.01877
new 0.15374 -0.13497 0.01877
-------------- 83 --------------magnetic_moment
old 0.15374 -0.01877 0.13497
new 0.15374 -0.01877 0.13497
-------------- 84 --------------magnetic_moment
old 0.13497 -0.15374 -0.01877
new 0.13497 -0.15374 -0.01877
-------------- 85 --------------magnetic_moment
old 0.01877 -0.15374 -0.13497
new 0.01877 -0.15374 -0.13497
-------------- 86 --------------magnetic_moment
old -0.13497 -0.01877 -0.15374
new -0.13497 -0.01877 -0.15374
-------------- 87 --------------magnetic_moment
old 0.01877 0.13497 0.15374
new 0.01877 0.13497 0.15374
-------------- 88 --------------magnetic_moment
old 0.00000 0.00000 0.00000
new 0.00000 -0.00000 0.00000
-------------- 89 --------------magnetic_moment
old 0.02810 0.00000 0.02810
new 0.02810 0.00000 0.02810
-------------- 90 --------------magnetic_moment
old 0.00000 -0.02810 -0.02810
new 0.00000 -0.02810 -0.02810
-------------- 91 --------------magnetic_moment
old 0.00000 0.00000 0.00000
new 0.00000 0.00000 0.00000
**-------------- 92 --------------magnetic_moment
old 0.00000 -0.21382 0.21382
new 0.00000 -0.20626 0.20626
-------------- 93 --------------magnetic_moment
old 0.21382 0.00000 -0.21382
new 0.20626 -0.00000 -0.20626
-------------- 94 --------------magnetic_moment
old -0.21382 0.21382 0.00000
new -0.20626 0.20626 0.00000**
-------------- 95 --------------magnetic_moment
old 0.15374 -0.13497 -0.01877
new 0.15374 -0.13497 -0.01877
-------------- 96 --------------magnetic_moment
old -0.15374 0.01877 0.13497
new -0.15374 0.01877 0.13497
-------------- 97 --------------magnetic_moment
old 0.13497 -0.15374 0.01877
new 0.13497 -0.15374 0.01877
-------------- 98 --------------magnetic_moment
old -0.01877 0.15374 -0.13497
new -0.01877 0.15374 -0.13497
-------------- 99 --------------magnetic_moment
old -0.13497 -0.01877 0.15374
new -0.13497 -0.01877 0.15374
-------------- 100 --------------magnetic_moment
old 0.01877 0.13497 -0.15374
new 0.01877 0.13497 -0.15374
-------------- 101 --------------magnetic_moment
old 0.00000 0.00000 0.00000
new -0.00000 0.00000 0.00000
-------------- 102 --------------magnetic_moment
old 0.00000 0.02810 -0.02810
new 0.00000 0.02810 -0.02810
-------------- 103 --------------magnetic_moment
old -0.02810 0.00000 0.02810
new -0.02810 0.00000 0.02810
-------------- 104 --------------magnetic_moment
old 0.02810 -0.02810 0.00000
new 0.02810 -0.02810 0.00000
-------------- 105 --------------magnetic_moment
old 0.00000 0.00000 0.00000
new 0.00000 0.00000 0.00000
```
Same problem for f functions as above and in addition we have issues with <s| |fxxx>, <s| |fyyy>, and <s| |fzzz> integrals.
Other issues can be revealed in the same way. Can you help me with this?
Many thanks,
Marc
--
You are receiving this because you are subscribed to this thread.
Reply to this email directly or view it on GitHub:
#50
--
|
|
Thank you for your answer but this is not a problem of ordering because I have already taken care to convert libcint orbital ordering to stda one. This is also why overlap, dipole moment in both length and velocity representations are computed correctly and that I am able to compute oscillator strengths and polarizabilities also correctly. Here, if I take the first example for the x component. The matrix should be : |
|
There might be certain convention difference. The following code also gives the
x component of int1e_cg_irxp. Hope it explains what has been computed.
```
mol = pyscf.M(atom='H', basis=[[2,[2.062,1]]], spin=None, cart=True)
norm = mol.intor('int1e_ovlp').diagonal()**-.5
rp = -1j * mol.intor('int1e_irp').reshape(3, 3, 6, 6)
y_pz = rp[1,2]
z_py = rp[2,1]
py_z = z_py.conj().T
rxp_x = y_pz - py_z
print(np.einsum('ij,i,j->ij', (rxp_x * 1j).real, norm, norm).round(3))
```
…On Tue, Jan 12, 2021 at 02:29:30PM -0800, Marc de Wergifosse wrote:
Thank you for your answer but this is not a problem of ordering because I have already taken care to convert libcint orbital ordering to stda one. This is also why overlap, dipole moment in both length and velocity representations are computed correctly and that I am able to compute oscillator strengths and polarizabilities also correctly.
Here, if I take the first example for the x component. The matrix should be :
[[ 0. , 0. , 0. , 0. , -0. , 0. ],
[ 0. , 0. , 1. , 0. , 0. , 0. ],
[ 0. , -1. , 0. , 0. , 0. , 0. ],
[ 0. , 0. , 0. , 0. , 1.155, 0. ],
[-0. , 0. , 0. , -1.155, 0. , **-1.155**],
[ 0. , 0. , 0. , 0. , **1.155**, 0. ]]
where both sign are reversed with respect to what you obtained.
This is what our integral subroutine are giving.
--
You are receiving this because you commented.
Reply to this email directly or view it on GitHub:
#50 (comment)
--
|
|
Thanks a lot for helping me with this! I have checked this and this is not at the level of the vector product that the difference occurs. 0.000 | 0.000 | 0.000 | 0.000 | 0.167 | 0.000 and for z_py 0.000 | 0.000 | 0.000 | 0.000 | 0.167 | 0.000 But with stda subroutines, elements in bold have the sign reversed. |
|
I don't have ideas for other possibilities of the sign difference. Another check
for integral int1e_irp can be the integral of int1e_ipovlp (<\nabla i|j>)
between f_yyz and d_zz. For libcint, it's like
```
mol = pyscf.M(atom='H', basis=[[2,[2.062,1]],[3,[2.062,1]]],spin=None, cart=True)
y_pz = -1j*mol.intor_by_shell('int1e_irp', shells=(0,0)).reshape(3,3,6,6)[1,2]
pz = 1j*mol.intor_by_shell('int1e_ipovlp', shells=(1,0)).reshape(3,10,6)[2]
print(y_pz[4,5], pz[6,5])
```
The conventions for operators and basis are consistent if you compute all
integrals with libcint. I guess it may cause errors if you mix the two integral
libraries in the same program.
…On Wed, Jan 13, 2021 at 07:30:25AM -0800, Marc de Wergifosse wrote:
Thanks a lot for helping me with this! I have checked this and this is not at the level of the vector product that the difference occurs.
The difference is already there when computing i<k|rp|n>.
Still with the same example for the x component, with libcint, y_pz is
0.000 | 0.000 | 0.000 | 0.000 | 0.167 | 0.000
0.000 | 0.000 | 0.167 | 0.000 | 0.000 | 0.000
0.000 | -0.167 | 0.000 | 0.000 | 0.000 | 0.000
0.000 | 0.000 | 0.000 | 0.000 | 0.500 | 0.000
-0.167 | 0.000 | 0.000 | -0.500 | 0.000 | **0.167**
0.000 | 0.000 | 0.000 | 0.000 | **-0.167** | 0.000
and for z_py
0.000 | 0.000 | 0.000 | 0.000 | 0.167 | 0.000
0.000 | 0.000 | -0.167 | 0.000 | 0.000 | 0.000
0.000 | 0.167 | 0.000 | 0.000 | 0.000 | 0.000
0.000 | 0.000 | 0.000 | 0.000 | -0.167 | 0.000
-0.167 | 0.000 | 0.000 | 0.167 | 0.000 | **0.5**
0.000 | 0.000 | 0.000 | 0.000 | **-0.5** | 0.000
But with stda subroutines, elements in bold have the sign reversed.
--
You are receiving this because you commented.
Reply to this email directly or view it on GitHub:
#50 (comment)
--
|
|
Unfortunately, I cannot compute <\nabla i|j> with our more than 20 years old stda integral deck (already used before stda), at least in a straightforward manner. I am just comparing results between both integral decks not mixing them. My motivation for this is because when I use only libcint only, optical rotations and rotary strengths are wrong. Meaning that e.g. CH4 has a non-zero optical rotation value. The interfacing is pretty simple in a way since I only need to feed data to the library and extract integrals with values in the right order. Then it transforms them to the MO basis to uses it. Since it is working well for the overlap and the dipole moments (for real systems, with r and nabla operators, both length and velocity representations), there is no reason it should not work for the magnetic moment. Here, the difference is that the magnetic moment is skew-symmetric and not symmetric. But this should not be affected by the interfacing. Yesterday, I did a lot of testing and printing, again. In the int1e_irp subroutine, the sign is actually given by the nabla part, g1 computed by G1E_D_J. I thought a bit about this and in fact, to be sure we are right, we can reformulate the problem for the first example as this: while with stda, I have In other words, libcint gives i<dyy|(r-R)Xp|dyz>=-<dzz| |dyz>=<dyz| |dzz>=-<dyz| |dyy> while stda gives i<dyy|(r-R)Xp|dyz>=<dzz| |dyz>=-<dyz| |dzz>=-<dyz| |dyy> So, if I look to orbital symmetries, they are all in the same yz plan where going from <dyy| |dyz> to <dzz| |dyz>, dyy has rotated of 90° to result to dzz. For libcint, doing so changes the integral sign while for stda, this is not the case. So the question is which one is right and for what reason the other one gives the reverse sign? Of course, this is not answering other examples for which some of the values are different. |
|
The differences in optical properties may be caused by other interface problems
I guess. One-electron integral can be symoblically calculated. The x component
of i<dzz|rxp|dyz> using maxima
```
assume(a > 0);
dyz: y * z * exp(-a*(x^2+y^2+z^2));
dzz: z * z * exp(-a*(x^2+y^2+z^2));
integrate(
integrate(
integrate(dzz * (y * diff(dyz, z) - z * diff(dyz, y)),
x, minf, inf),
y, minf, inf),
z, minf, inf);
```
The result is -pi^1.5/2^4.5/a^3.5
…On Fri, Jan 15, 2021 at 12:17:39AM -0800, Marc de Wergifosse wrote:
Unfortunately, I cannot compute <\nabla i|j> with our more than 20 years old stda integral deck (already used before stda), at least in a straightforward manner. I am just comparing results between both integral decks not mixing them. My motivation for this is because when I use only libcint only, optical rotations and rotary strengths are wrong. Meaning that e.g. CH4 has a non-zero optical rotation value. The interfacing is pretty simple in a way since I only need to feed data to the library and extract integrals with values in the right order. Then it transforms them to the MO basis to uses it. Since it is working well for the overlap and the dipole moments (for real systems, with r and nabla operators, both length and velocity representations), there is no reason it should not work for the magnetic moment.
Here, the difference is that the magnetic moment is skew-symmetric and not symmetric. But this should not be affected by the interfacing. Yesterday, I did a lot of testing and printing, again. In the int1e_irp subroutine, the sign is actually given by the nabla part, g1 computed by G1E_D_J.
I thought a bit about this and in fact, to be sure we are right, we can reformulate the problem for the first example as this:
With libcint, I obtain this matrix for the x component:
[[ 0. , 0. , 0. , 0. , -0. , 0. ],
[ 0. , 0. , 1. , 0. , 0. , 0. ],
[ 0. , -1. , 0. , 0. , 0. , 0. ],
[ 0. , 0. , 0. , 0. , **1.155**, 0. ],
[-0. , 0. , 0. , **-1.155**, 0. , **1.155**],
[ 0. , 0. , 0. , 0. , **-1.155**, 0. ]]
while with stda, I have
[[ 0. , 0. , 0. , 0. , -0. , 0. ],
[ 0. , 0. , 1. , 0. , 0. , 0. ],
[ 0. , -1. , 0. , 0. , 0. , 0. ],
[ 0. , 0. , 0. , 0. , **1.155**, 0. ],
[-0. , 0. , 0. , **-1.155**, 0. , **-1.155**],
[ 0. , 0. , 0. , 0. , **1.155**, 0. ]]
In other words, libcint gives
i<dyy|(r-R)Xp|dyz>=-<dzz| |dyz>=<dyz| |dzz>=-<dyz| |dyy>
while stda gives
i<dyy|(r-R)Xp|dyz>=<dzz| |dyz>=-<dyz| |dzz>=-<dyz| |dyy>
So, if I look to orbital symmetries, they are all in the same yz plan where going from <dyy| |dyz> to <dzz| |dyz>, dyy has rotated of 90° to result to dzz. For libcint, doing so changes the integral sign while for stda, this is not the case. So the question is which one is right and for what reason the other one gives the reverse sign?
Of course, this is not answering other examples for which some of the values are different.
--
You are receiving this because you commented.
Reply to this email directly or view it on GitHub:
#50 (comment)
--
|
|
Yes you are right i<dzz|rxp|dyz> is positive while i<dyy|rxp|dyz> is negative. |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
I am trying to modernize the stda program by removing its old integral deck, replacing it by libcint, also allowing new types of integrals that were not available. This is working perfectly fine for the overlap, dipole moment, and <i|nabla|j> integrals for which I can reproduce values obtained by our integral deck (see https://github.com/grimme-lab/stda/blob/master/intslvm.f and https://github.com/grimme-lab/stda/blob/master/intpack.f90).
But I have some issues with i<k|(r-Rc)Xp|n> type of integrals while using exactly the same interface, just replacing the name of the subroutine by cint1e_cg_irxp_cart. Note that everything is done with cartesian basis functions and that I have well defined Rc in env.
For integrals on the same atom, few elements have the reverse sign while others are ok for this skew-symmetric matrix.
For <d| |d> integrals, <dyy| |dxy>, <dzz| |dxz>, and <dzz| |dyz> signs are swapped with their counterpart.
For <f| |f> integrals, I have the same issue with <fzzz| |fxxy>, <fyyy| |fxxz>, <fyyy| |fyyx>, <fzzz| |fyyx>, <fyyy| |fxzz>, <fzzz| |fxzz>, <fzzz| |fyzz>, <fyyz| |fxzz>, <fyyz| |fxyz>, <fxzz| |fxyz>, and <fyzz| |fxyz>. All others are ok.
In addition to this, for integrals between two different atoms, I obtained different values than with our integral deck for <s| |fxxx>, <s| |fyyy>, <s| |fzzz>, <p| |fxxx>, <p| |fyyy>, <p| |fzzz>, <dxx| |dxx>, <dyy| |dyy>, <dzz| |dzz>, and almost all <f| |f>. Note that for <dxx| |dxx>, <dyy| |dyy>, and <dzz| |dzz> the sign is also wrong.
Those are what I was able to identify with small models.
Here are some examples:
1° case
For H- with one set of d orbitals:
These are the upper triangle of the skew-symmetric matrix where d orbitals are ordered as dxx dyy dzz dxy dxz dyz. Thus, 8, 13, and 18 are <dyy| |dxy>, <dzz| |dxz>, and <dzz| |dyz>, respectively.
2° case
Same system with one set of f functions.
The f functions are ordered as fxxx, fyyy, fzzz, fxxy, fxxz, fyyx, fyyz, fxzz, fyzz, and fxyz. As for the case above, the upper diagonal is printed and problematic elements are <fzzz| |fxxy>, <fyyy| |fxxz>, <fyyy| |fyyx>, <fzzz| |fyyx>, <fyyy| |fxzz>, <fzzz| |fxzz>, <fzzz| |fyzz>, <fyyz| |fxzz>, <fyyz| |fxyz>, <fxzz| |fxyz>, and <fyzz| |fxyz>.
3° case
For CH4 with s functions on H and one set of f functions for the carbon.
Same problem for f functions as above and in addition we have issues with <s| |fxxx>, <s| |fyyy>, and <s| |fzzz> integrals.
Other issues can be revealed in the same way. Can you help me with this?
Many thanks,
Marc
The text was updated successfully, but these errors were encountered: