Describe the bug
When evaluating coulomb integrals via electron_repulsion_integral, it seems that the function returns different values for bases generated through the from_pyscf wrapper and the make_contractions function. It looks like the integral function is evaluating the bases' shells in different orders since the correct values appear in the matrix.

To Reproduce
basis_bug.pdf

Expected behavior
The integral should return the same value regardless of the method used to generate the basis.

Is your feature request related to a problem? Please describe.

This came up when I was tuning the from_iodata function in gbasis/wrappers.py. At the moment, this function has two return values: basis, which is a list of GeneralizedContractionShell instances and coord_types, which says for each shell whether it is Cartesian or pure. I would suggest to make the elements of coord_types attributes of the corresponding GeneralizedContractionShell instances. (This can be done in various was effectively, not necessarily literally the way I wrote.)

The motivation for this request is that the coord_types part is in practice always tied to the rest of the definition of the basis set. (It is normally not something that is varied during a calculation.) At the moment, when using GBasis, one has to keep track of both variables and pass then into the eval functions every time, which seems redundant and error prone. Forgetting to pass in coord_types would essentially always be a bug.

Is there a specific integral/formula that you would like implemented?

No

Is there a change to the code or algorithm you would like to see?

To be discussed. I can look into options to streamline both variables into one. My first impression is that one could, at the same timel simplify and reduce a lot of the code in the base*.py modules.

Additional context

None

Describe the bug

Big memory usage for density evaluation

To Reproduce

I attach the system I used to test:
iron_hexacarbonyl.tar.gz

mol=load_one("iron_hexacarbonyl.fchk")
dm = mol.one_rdms['scf']
basis, coord_types = from_iodata(mol)
becke = BeckeWeights(order=3)

oned = GaussChebyshevType2(100)
rgrid = BeckeTF(1e-4, 1.5).transform_1d_grid(oned)
grid = MolGrid.horton_molgrid(mol.atcoords, mol.atnums, rgrid, 1030, becke)

evaluate_density(dm, basis, grid.points, coord_type=coord_types)

Expected behaviour

I don't know if this could be improved, I am just reporting the behaviour as an issue.

Screenshots

System Information:

• OS: 18.04.5 LTS
• Python version: 3.6.10
• NumPy version: 1.18.1
• SciPy version: 1.4.1

Additional context

For computing the density it starts with evaluate_density(), then evaluates each orbital by computing the density for the basis functions and then it obtains the total density by using the density matrix. For evaluating the basis functions density, it is done in _deriv.py module inside _eval_deriv_contractions(). To monitor the memory I used several prints with nbytes for arrays and @profile from memory_profiler import profile. To monitor the memory of each contraction I printed the associated memory of matrix_contraction in base one

For a radial grid =100points
angular = 1030
Total = 100 * 1202 * 13atoms = 1562600

Base memory for S-generalized contractions = 12.49MB

The values I obtained profiling the code are shown in the following pictures:

this is the profiling for evaluate density()

You can see evaluates basis increases almost the same memory that it has been calculated in line 78.

But also, there's a spike of memory that it is not shown here that occurs in evaluate_density_using_evaluated_orbs. The profiling for this is shown in the following image.

Here you can see in line 63 where the memory is doubled.

Is your feature request related to a problem? Please describe.
For standard integrals, it would be beneficial to support fast C++ libraries, which will likely remain faster than gbasis regardless of any optimizations we make.

Additional context

There are many possible libraries: libint2, libcint, libmint, libreta, ...

Psi4Numpy uses libmints to access the libint library

Is your feature request related to a problem? Please describe.

I would like a way of specifying the spherical order when evaluating basis-functions or integrals.

For example, IOData specifies the order for D-type orbitals as : ['c0', 'c1', 's1', 'c2', 's2'], whereas PySCF would be ["s2", "s1", "c0", "c1", "c2"]. Similarly for F-type orbitals, IOData does ["s3", "s2", "s1", "c0", "c1", "c2", "c3"] and PySCF is ['c0', 'c1', 's1', 'c2', 's2', 'c3', 's3'], where "c" are cosine-type functions and "s" are sin-type functions.

It appears that all the functionality is there to begin with, the relevant functions are

• The function generate_transform converts from one Cartesian order to a specified Spherical order.
• The default Spherical order in gbasis is the same as PySCF and it is specified by angmom_components_sph inside the class GeneralizedContractionShell.
• When using iodata, you'll have to use the gbasis functoion from_iodata which returns a list of IODataShell, where the function angmom_components_sph is re-written to make it have IOData order.

Is there a change to the code or algorithm you would like to see?

Would like to change the API to better specify the order. Couple of options come to mind:

• specifying the order when evaluating anything, e.g. evaluate_basis, electron_repulsion_integral. The downside of this, you'll have to change all functions that does evaluation/integrals.
• add an option to the constructor of the class GeneralizedContractionShell and wrappers like parse_nwchem, parse_gbs, from_iodata function that specifies common spherical order based on a string: .e.g "gbasis/PySCF", "iodata", "gaussian??", "Psi4?". The downside is that it isn't general.
• add an option to the constructor of the class GeneralizedContractionShell and from_iodata function where the user can specify their own order by re-writting angmom_components_sph function. The downside is that it might be difficult for the user to specify their own order.

Additional context

• One can use IOdata "convert_convention_shell" to create a standard test.

It should have a more welcoming tone. Code of conduct may be moved to a different location.

Describe the bug
The factorial2 function in SciPy returns different values for negative inputs depending on the version. For versions older than 1.11.0, it returns 1, while for versions 1.11.0 and newer, it returns 0.

Code that utilizes the factorial2 function may encounter potential issues, such as those found in the contractions.py file:

    @property
    def norm_prim_cart(self):
        r"""Return the normalization constants of the Cartesian Gaussian primitives.

        For a Cartesian primitive with exponent :math:`\alpha_i`, the normalization constant is:

        .. math::
           N(\alpha_i, \vec{a}) = \sqrt {
           \left(\frac{2\alpha_i}{\pi}\right)^\frac{3}{2}
           \frac{(4\alpha_i)^{a_x + a_y + a_z}}{(2a_x - 1)!! (2a_y - 1)!! (2a_z - 1)!!}}

        Returns
        -------
        norm_prim_cart : np.ndarray(L, K)
            The normalization constants of the Cartesian Gaussian primitives.
            `L` is the number of contracted Cartesian Gaussian functions for the given angular
            momentum, i.e. :math:`(\ell + 1) * (\ell + 2) / 2`
            `K` is the number of exponents (i.e. primitives).

        """
        exponents = self.exps[np.newaxis, :]
        angmom_components_cart = self.angmom_components_cart[:, :, np.newaxis]

        return (
            (2 * exponents / np.pi) ** (3 / 4)
            * ((4 * exponents) ** (self.angmom / 2))
            / np.sqrt(np.prod(factorial2(2 * angmom_components_cart - 1), axis=1))
        )

An unwanted value will be returned when angmom_components_cart is 0, i.e., $s$ orbitals.

To Reproduce
Use SciPy with version $\geq 1.11.0$

Describe the bug
I was integrating density^(1/3) when I noticed that the electron density computed by evaluate_density gives a negative, but very small value, for faraway points (e.g. -2.77765871e-128).

To Reproduce
Helium FCHK file: helium.fchk.tar.gz

import numpy as np

from iodata import load_one

from grid.molgrid import MolGrid
from grid.becke import BeckeWeights
from grid.rtransform import BeckeTF
from grid.onedgrid import GaussChebyshevType2

from gbasis.wrappers import from_iodata
from gbasis.evals.density import evaluate_density

# load molecule
mol = load_one("helium.fchk")
dm = mol.one_rdms['scf']
print("molecular charge   = ", mol.charge)
print("atomic numbers     = ", mol.atnums)
print("atomic coordinates = ", mol.atcoords)

# make grid
oned = GaussChebyshevType2(150)
rgrid = BeckeTF(1.0e-6, 1.5).transform_1d_grid(oned)
grid = MolGrid.horton_molgrid(mol.atcoords, mol.atnums, rgrid, 350, BeckeWeights(order=3))

# evaluate density
basis, coord_types = from_iodata(mol)
dens = evaluate_density(dm, basis, grid.points, coord_type=coord_types)
print("min & max density = ", np.min(dens), np.max(dens))

index = np.where(dens < 0.0)
print("density     at points with negative densiry = ", dens[index])
print("coordinates of points with negative density = ", grid.points[index])

The output is:

molecular charge   =  0.0
atomic numbers     =  [2]
atomic coordinates =  [[0. 0. 0.]]
min & max density =  -2.77765870781763e-128 2.930734188084798
density     at points with negative densiry =  [-2.77765871e-128]
coordinates of points with negative density =  [[ 10.15753086 -15.57090172 -10.15753086]]

Expected behaviour
The electron density is expected to be positive at every point.

System Information:

• OS: macOS
• Python version: 3.6.10
• NumPy version: 1.18.1
• SciPy version: 1.4.1

Additional context
I am using the snippet below as a workaround, until this get fixed.
dens[np.where((dens < 0.0) & (abs(dens) < 1.0e-15))] = 0.0

Given 4 shells and centers compute integrals.

Your function electron_repulsion_integral agrees with PySCF:

import pyscf 
import gbasis

mol        = pyscf.gto.M(atom="H 0 0 0; H 0 0 1; ")
pyscf_eri  = mol.intor("int2e")
gbasis_eri = gbasis.electron_repulsion.electron_repulsion_integral(gbasis.wrappers.from_pyscf(mol))

assert np.allclose( pyscf_eri, gbasis_eri ) # this is True

I need the derivative of the two electron integral which in PySCF is:

pyscf_eri_deriv = mol.intor("int2e_ip1")

I've rewritten your electron_repulsion_integral from numpy to jax.numpy so I can get the Jacobian in one line of code:

import jax
J_eri = jax.jacfwd( jax_electron_repulsion_integral ) ( atom_positions )

The problem: pyscf_eri_deriv and J_eri almost perfectly agree, but there are two entries that are different.

Feature Request: An implementation of the derivative of electron_repulsion_integral such that the output is np.allclose to mol.intor("int2e_ip1") from PySCF.

Note: Using jax.jacfwd I managed to get jax.numpy implementations of your code to agree with pyscf (which just uses libcint) for

• gbasis.overlap.overlap_integral
• gbasis.kinetic_energy.kinetic_energy.integral
• gbasis.point_charge.point_charge_integral

Check: https://github.com/theochem/cgbasis/blob/master/python/gbasis/gobasis.py

This issue is to request a set of Jupyter notebooks to show the functionality/features of Gbasis package. This follows recent decisions to incorporate notebooks for qc-devs packages that are easy to use and modify so all the features can be accessed readily.

Add functions for evaluation on a grid using obasis, points, & dm:

• molecular orbitals (doesn't need dm)
• electron density
• gradient of electron density
• hessian of electron density
• higher order derivatives of electron density (i.e., arbitrary order derivative of Gaussian basis)
• laplacian of electron density
• positive definite kinetic energy density
• derivative of positive definite kinetic energy density
• electrostatic potential

API:

compute_grid_mo(points, obasis, mo_coeffs)
compute_grid_density(points, obasis, dm)

Is there a specific integral/formula that you would like implemented?

We currently have the Ehrenfest Force/Hessian implemented using the virial theorem (via the stress tensor). It would be nice to have a direct implementation, using (fancier) integrals. This would have better numerical characteristics, as the virial-theorem-based definition has spurious nodes resulting from the fourth (force) and fifth (Hessian) derivatives of GTOs.

Additional context
James Anderson has notes on this, so the formulas are explicit (and worked out).

Is there a change to the code or algorithm you would like to see?
Make the functions/methods consistent between the original gbasis implementation and the lincint interface. For example, CBasis.overlap should change to CBasis.overlap_integral. @msricher can you please fix these to add _integral to all methods?

Describe the bug

The method construct_array_mix in base_two_symm.py, base_two_asymm.py and base_four.py constructs transform_* variables when they are not needed. When the wrapper to generate the basis set does not define the transformation, this will result in an unwarranted error message. This typically happens when loading wavefunction files with IOData from file formats that always work with Cartesian shells, but there are other cases too, see example below.

A similar problem was fixed in #104 for base_one.py.

To Reproduce

The following script (with IOData installed) generates the exception:

from importlib.resources import path
from iodata import load_one
from gbasis.wrappers import from_iodata
from gbasis.integrals.overlap import overlap_integral

with path("iodata.test.data", "he2_ghost_psi4_1.0.molden") as fn_full:
        iodata = load_one(str(fn_full))
basis, coord_types = from_iodata(iodata)
olp = overlap_integral(basis, coord_type=coord_types)
print(olp)

Exception:

Traceback (most recent call last):
  File "example.py", line 9, in <module>
    olp = overlap_integral(basis, coord_type=coord_types)
  File "/home/toon/miniconda3/envs/horton3/lib/python3.7/site-packages/gbasis/integrals/overlap.py", line 145, in overlap_integral
    return Overlap(basis).construct_array_mix(coord_type)
  File "/home/toon/miniconda3/envs/horton3/lib/python3.7/site-packages/gbasis/base_two_symm.py", line 304, in construct_array_mix
    cont_one.angmom_components_sph,
  File "/home/toon/miniconda3/envs/horton3/lib/python3.7/site-packages/gbasis/wrappers.py", line 117, in angmom_components_sph
    "momentum {0}".format(self.angmom)
ValueError: Given convention does not support spherical contractions for the angular momentum 0

Expected behaviour

[[1.00000146e+00 5.95216385e-01 2.59436827e-04 4.41546274e-02]
 [5.95216385e-01 1.00000000e+00 4.41546274e-02 2.14669756e-01]
 [2.59436827e-04 4.41546274e-02 1.00000146e+00 5.95216385e-01]
 [4.41546274e-02 2.14669756e-01 5.95216385e-01 1.00000000e+00]]

Screenshots

System Information:

• OS: Fedora 32
• Python version: 3.7.7
• NumPy version: 1.18.5
• SciPy version: 1.5.0

Additional context

The fix for this problem should be relatively easy, essentially replicating the changes to base_one.py (and its unit tests) made in #104 in all the other base_*.py. The unit tests do not need to make use of iodata.

Convert the notes.pdf to a scientific documentation on the website, excluding the code-snippet examples.

Is your feature request related to a problem? Please describe.
In the docs, we mention that the user can use the transform key to perform linear transformations of basis functions. Later on in the notes, we show that it's possible to use AO to MO and MO to AO transformations for all of our high-level functions. It would be useful if we actually provided the functionality to generate these transformations, instead of leaving that up to the user to generate.

Currently we don't support ECPs in gbasis. It would be nice to do that. Given the (new) strategy of falling back on a non-Python implementation where available as necessary, it may be easier to do that. There is an open-source library for ECP integrals (C++).
https://github.com/robashaw/libecpint

For now I think this is low(ish) priority but this is a placeholder so we don't forget. It's something that I think can be supported relatively easily when/if it is needed.

Describe the bug

As I was making the tutorial notebook that I've proposed in PR #142, I believe I came across a bug that comes about when trying to compute the electrostatic potential for a system using the density matrix expressed with respect to molecular orbitals. I passed a density matrix expressed w.r.t. molecular orbitals of shape (14, 14), whereas the density matrix expressed w.r.t. atomic orbitals would've been of shape (18, 18). Even though I provided a transform argument (a (14, 18) matrix) to tell this higher level function that I want to work in terms of the MO basis, it threw an error saying that my one-electron density matrix was the wrong shape (it was expected to be (18, 18)). I believe the problem is that the presence or absence of the transform matrix when calling electrostatic_potential is irrelevant when it comes to checking the shape of the one-electron density matrix.

I don't think there was a test in gbasis/tests to check this behaviour, despite it being apparently possible to carry out this operation in the notes.

To Reproduce

*see screenshot

Expected behaviour

I would expect no exception of this sort to be raised when the density matrix expressed w.r.t. molecular orbitals is passed to electrostatic_potential along with a proper transform argument. This can likely be (at least partially) fixed by modifying the conditional statements 1. here, 2. here, and 3. here to include an and transform is None. If transform is None, then the shape of the one-electron density matrix should be (# AOs, # AOs). However, when transform is passed, I believe it would have to be of shape (# MOs, # MOs). This is confirmed in the documentation for electrostatic_potential.py. Right now, it doesn't look like it matters if transform is passed or not; the one-electron density matrix is always expected to be of shape (# AOs, # AOs) according to the sum(cont.num_sph * cont.num_seg_cont for cont in basis) != one_density_matrix.shape[0] clause.

Screenshots

System Information:

• OS:
• Python version:
• NumPy version:
• SciPy version:

Additional context

Describe the bug
The generated transformation matrix for cartesian to spherical contractions/primitives is different from the one from HORTON (https://github.com/theochem/horton/blob/master/horton/gbasis/cartpure.cpp) for angular momentum 4 and 5 (and probably all of the higher angular momentums).

In the case of angular momentum 4, gbasis's transformation matrix contains an extra element at position (3, 3), which corresponds to C_{42} (or m=2) and X(xxyy) (or (2, 2, 0)).

In the case of angular momentum 5, gbasis's transformation matrix contains extra elements at positions (3, 7), which correspond to C_{52} (m=2) and X(xxyyz) ((2, 2, 1)) and (6,6) which corresponds to S_{53} (m=-3) and X(xxyyy) ((2, 3, 0)). horton's transformation matrix contains an extra element at position (5, 3) which corresponds to C_{53} (m=3) and X(xxxyy) ((3, 2, 0)).

To Reproduce
Follow code was used to create the HORTON's transformation matrix:

import itertools as it

from gbasis.spherical import generate_transformation
import numpy as np


# taken from HORTON's gbasis/cartpure.cpp
horton_transform = [
    [0,  0, 0.375],
    [0,  3, 0.21957751641341996535],
    [0,  5, -0.87831006565367986142],
    [0, 10, 0.375],
    [0, 12, -0.87831006565367986142],
    [0, 14, 1.0],
    [1,  2, -0.89642145700079522998],
    [1,  7, -0.40089186286863657703],
    [1,  9, 1.19522860933439364],
    [2,  4, -0.40089186286863657703],
    [2, 11, -0.89642145700079522998],
    [2, 13, 1.19522860933439364],
    [3,  0, -0.5590169943749474241],
    [3,  5, 0.9819805060619657157],
    [3, 10, 0.5590169943749474241],
    [3, 12, -0.9819805060619657157],
    [4,  1, -0.42257712736425828875],
    [4,  6, -0.42257712736425828875],
    [4,  8, 1.1338934190276816816],
    [5,  2, 0.790569415042094833],
    [5,  7, -1.0606601717798212866],
    [6,  4, 1.0606601717798212866],
    [6, 11, -0.790569415042094833],
    [7,  0, 0.73950997288745200532],
    [7,  3, -1.2990381056766579701],
    [7, 10, 0.73950997288745200532],
    [8,  1, 1.1180339887498948482],
    [8,  6, -1.1180339887498948482],
]
answer = np.zeros((9, 15))
for i in horton_transform:
    answer[i[0], i[1]] = i[2]
test = (
    np.abs(
        generate_transformation(
            4,
            np.array(
                [
                    (i.count(0), i.count(1), i.count(2))
                    for i in it.combinations_with_replacement(range(3), 4)
                ]
            ),
            (0, 1, -1, 2, -2, 3, -3, 4, -4),
            "left",
        ),
    )
    - np.abs(answer)
)
print(
    np.where((test > 1e-9))
)
print(
    np.where((test < -1e-9))
)

# taken from HORTON's gbasis/cartpure.cpp
horton_transform = [
    [0,  2, 0.625],
    [0,  7, 0.36596252735569994226],
    [0,  9, -1.0910894511799619063],
    [0, 16, 0.625],
    [0, 18, -1.0910894511799619063],
    [0, 20, 1.0],
    [1,  0, 0.48412291827592711065],
    [1,  3, 0.21128856368212914438],
    [1,  5, -1.2677313820927748663],
    [1, 10, 0.16137430609197570355],
    [1, 12, -0.56694670951384084082],
    [1, 14, 1.2909944487358056284],
    [2,  1, 0.16137430609197570355],
    [2,  6, 0.21128856368212914438],
    [2,  8, -0.56694670951384084082],
    [2, 15, 0.48412291827592711065],
    [2, 17, -1.2677313820927748663],
    [2, 19, 1.2909944487358056284],
    [3,  2, -0.85391256382996653193],
    [3,  9, 1.1180339887498948482],
    [3, 16, 0.85391256382996653193],
    [3, 18, -1.1180339887498948482],
    [4,  4, -0.6454972243679028142],
    [4, 11, -0.6454972243679028142],
    [4, 13, 1.2909944487358056284],
    [5,  0, -0.52291251658379721749],
    [5,  3, 0.22821773229381921394],
    [5,  5, 0.91287092917527685576],
    [5, 10, 0.52291251658379721749],
    [5, 12, -1.2247448713915890491],
    [6,  1, -0.52291251658379721749],
    [6,  6, -0.22821773229381921394],
    [6,  8, 1.2247448713915890491],
    [6, 15, 0.52291251658379721749],
    [6, 17, -0.91287092917527685576],
    [7,  2, 0.73950997288745200532],
    [7,  7, -1.2990381056766579701],
    [7, 16, 0.73950997288745200532],
    [8,  4, 1.1180339887498948482],
    [8, 11, -1.1180339887498948482],
    [9,  0, 0.7015607600201140098],
    [9,  3, -1.5309310892394863114],
    [9, 10, 1.169267933366856683],
    [10,  1, 1.169267933366856683],
    [10,  6, -1.5309310892394863114],
    [10, 15, 0.7015607600201140098],
]
answer = np.zeros((11, 21))
for i in horton_transform:
    answer[i[0], i[1]] = i[2]

test = (
    np.abs(
        generate_transformation(
            5,
            np.array(
                [
                    (i.count(0), i.count(1), i.count(2))
                    for i in it.combinations_with_replacement(range(3), 5)
                ]
            ),
            (0, 1, -1, 2, -2, 3, -3, 4, -4, 5, -5),
            "left",
        ),
    )
    - np.abs(answer)
)
print(
    np.where((test > 1e-9))
)
print(
    np.where((test < -1e-9))
)

Expected behaviour
They are expected to be the same.

System Information:

• OS: Fedora 26
• Python version: Python 3.6.8
• NumPy version: 1.16.4
• SciPy version: 1.3.0

Additional context
I'm looking into this, but could @wilhadams check the gbasis.spherical part of the code for bugs and @tovrstra verify that the HORTON's transformation matrix is correct? It'd also be great if we are able to compare against another reference.

The num_contr property states that it returns the number of contractions in the ContractedCartesianGaussians class, however, this seems ambiguously worded. While the data class contains the information for every contracted GTO of a given angular momentum (l ), for a total of (l + 1)*(l + 2) / 2 CGTOs (which is the number given as num_contr), each CGTO is a contraction of a given number of Cartesian primitive GTOs, specified by the number of coefficients/exponents given. It seems like we should be careful to be consistent with our usage in the docstring.

Is your feature request related to a problem? Please describe.
The goal is to be able to evaluate the spherically-averaged density, its derivatives, and the kinetic-energy density for an atom.

Is there a specific integral/formula that you would like implemented?
There are several possible ways to do this, including using fractional occupation numbers or (probably easier) neglecting the spherical harmonic contribution when evaluating the density and its derivatives.

Is there a change to the code or algorithm you would like to see?
This might be a different set of routines for spherically averaging or maybe it is a flag for existing routines.

API:
compute_overlap(center_1, shell_1, convention_1, center_2, shell_2, convention_2)
compute_overlap_multicenter(shells, centers, convention)

Describe the bug

The function evaluate_density from module gbasis.evals.density fails to evaluate $\rho(r)$ in a grid of points for open-shell systems.

To Reproduce

# Load open shell system from fchk
#--------------------------------------------------
from iodata import load_one

# fchk of oxygen atom fetched from :
# https://github.com/theochem/chemtools/blob/master/chemtools/data/atom_08_O_N08_M3_ub3lyp_ccpvtz_g09.fchk
mol = load_one("atom_08_O_N08_M3_ub3lyp_ccpvtz_g09.fchk")

# Importing grid dependencies, and creating grid
#----------------------------------------------------------------------
# Make Becke-Lebedev molecular grid (using preset grid)
from grid.becke import BeckeWeights
from grid.molgrid import MolGrid
from grid.onedgrid import GaussChebyshev
from grid.rtransform import BeckeRTransform

oned = GaussChebyshev(100)
rgrid = BeckeRTransform(1e-4, 1.5).transform_1d_grid(oned)
grid = MolGrid.from_preset(mol.atnums, mol.atcoords, rgrid, "coarse", BeckeWeights())

# Evaluate the electron density for the grid points
#----------------------------------------------------------------------

from gbasis.wrappers import from_iodata
from gbasis.evals.density import evaluate_density

one_rdm = np.dot(mol.mo.coeffs * mol.mo.occs, mol.mo.coeffs.T)
basis = from_iodata(mol)
density = evaluate_density(one_rdm, basis[0], grid.points, coord_type=basis[1])

Expected behaviour

density should be a numpy.ndarray with the values of the electron density for each point of the grid.

Screenshots

System Information:

• OS: "openSUSE Leap 15.5"
• Python version: 3.11.4
• NumPy version: 1.25.2
• SciPy version: 1.11.1

Additional context

I am using your project in one of my packages.
I will be publishing it soon. I was wondering what's the best way of citing this work.

Thanks! :)

Avoid using *, +, and ().

Describe the bug

While I was fixing up the linter complaints. It mentioned that two tests didn't have assertion statements in them. If I did add them I get an error where the two arrays are not close at all. I haven't looked into but thought I put an issue for it.

As you can see below, the two tests don't have the obvious assert np.allclose(a, b) statement to them. If I add this fix, the shape is correct but the error is very large.

The first failed test is test_compute_differential_operator_integrals_multiarray in test_diff_operator_int.py

The second test is test_compute_multipole_moment_integrals_multiarray in test_moment_int.py:

To Reproduce

pytest -v .

Additional context

Is your feature request related to a problem? Please describe.

When incorporating gbaiss in chemtools, I noticed that gradient and hessian computations are much slower than anticipated when compared to horton2 which is currently used in chemtools. In gbasis, evaluating density seems to be pretty fast and I was hoping that gradient and hessian evaluations take 3 and 6 times as much. So, I am wondering whether there is room for improving the implementations.

Is there a specific integral/formula that you would like implemented?

@PaulWAyers has included a direct implementation of gradient and hessian in his notes.

Is there a change to the code or algorithm you would like to see?

If possible, it would be good to have a faster implementation of gradient and hessian computation, by either hard-coding these derivative and/or taking advantage of Hessian symmetries or using cython.

Additional context

Here are the code snippet and files to time the computations for H2O using cc-pVTZ basis on a cubic grid with 105,840 points (please unzip the attached files.tar.gz to get the h2o.fchk and h2o.npz). On my system, using time module, computing density, gradient, and hessian take ~0.8, ~10, ~55 seconds, respectively.

#!/usr/bin/env python

import time
import numpy as np

from iodata import load_one

from gbasis.wrappers import from_iodata
from gbasis.evals.density import evaluate_density
from gbasis.evals.density import evaluate_density_gradient, evaluate_density_hessian

# load points & density matrix
data = np.load("h2o.npz", allow_pickle=True)
points, dm = data["points"], data["dm"]
print("points shape = ", points.shape)

# build basis & compute
mol = load_one("h2o.fchk")
basis = from_iodata(mol)

start = time.time()
# dens = evaluate_density(dm, basis, points)
# grad = evaluate_density_gradient(dm, basis, points)
hess = evaluate_density_hessian(dm, basis, points)
end = time.time()
print("computation time = ", end - start)

files.tar.gz

Is your feature request related to a problem? Please describe.

At the moment, with the exception of the data class ContractedCartesianGaussians, everything is a function. We're starting to see some duplication of the code especially for input check and linear combinations of constructed arrays. We have composite function calls to retrieves the appropriate array to manipulate. At the moment, we don't have enough features for this to be too unwieldy, but it will be nice to build a class structure that pieces everything together nicely.

I can think of two (not necessarily exclusive) structures that seem to make sense:

1. Each class is associated with a mathematical entity. We already have ContractedCartesianGaussians for the Cartesian contractions, but there are also spherical AO's, MO's, other linear combinations of orbitals, and even densities that are all manipulated to construct other structures.
2. Each class is associated with the number of indices (associated with the Gaussian related object) needed to describe a given array. For example, the evaluations of the orbitals and their derivatives would be a one-index, and overlap of orbitals would be two-index. I guess in this sense, density and reduced density matrices would be zero-index. This lets us store develop a consistent API, for example for linearly transforming these indices.
3. Some combination of the two.

There are some questions that come to mind:

• Should the overlaps and other integrals computed via methods of the orbitals (contraction, AO, MO)?
• What about the evaluation of orbitals and their derivatives?
• Should the methods for producing these arrays be fed in orbitals instead?
• Which orbitals should be supported?
• Should the orbitals be represented with composite structure? i.e. need contractions to instantiate AO's, need AO's to instantiate MO's.
• Should arrays be represented with composite structure? e.g. need orbital evaluations to instantiate density evaluations

Though I like the idea of having objects that I can manipulate just as I would on paper, I think that these classes will get quite unwieldy and we will end up with a Russian doll situation whenever we want to retrieve a piece of low level information (e.g. exponents). We can store the low-level information directly in the object, but instantiating the object will be less intuitive. Seeing as how different orbitals have the same structure (they're just linear combinations of one another), we will have repeated information and structure that will make things difficult to read.

Of course, we can ignore these issues and simply say "we'll just make some porcelain (probably iodata) handle it" but the line between porcelain and plumbing is not too obvious if we expect our users to also be developers. If we want people to develop and contribute in a sustainable manner, we need to build a stable enough structure. Otherwise, we will need a gatekeeper to pore through each contribution.

I think this is quite hard, especially since our goal is quite fuzzy. If I'm not mistaken, the current goal is to

• Make a python only code that is not too slow
• Have enough features for the release of chemtools
• Be extendible enough to facilitate future scientific works
• Be well tested and documented
• Be stable enough to not have an active maintainer
  which isn't easy to accomplish within a short time period. At the very least, this will require someone who's invested enough in accomplishing these goals to spend a good amount of time learning about the science, thinking about the code structure, and be willing to rip apart the entire structure a couple of times. Without this person, we will likely end up with a messy mixture that will make it even harder to refactor in the future.

For now (next few days?), I might create classes according to the number of indices with specific structures based on my whims.

Describe the bug

I was trying to compute kinetic energy atomic integrals with gbasis.integrals.kinetic_energy. kinetic_energy_integral() with a fchk file where there are mix cartesian and spherical shells (
neon.fchk.tar.gz ) with the code in Reproduce section and got the next error:

File "issue_code.py", line 11, in <module>
    kinetic_integrals = kinetic_energy_integral(basis, coord_type = coord)
  File "/home/leila/gbasis/gbasis/gbasis/integrals/kinetic_energy.py", line 159, in kinetic_energy_integral
    return KineticEnergyIntegral(basis).construct_array_mix(coord_type)
  File "/home/leila/gbasis/gbasis/gbasis/base_two_symm.py", line 304, in construct_array_mix
    cont_one.angmom_components_sph,
  File "/home/leila/gbasis/gbasis/gbasis/wrappers.py", line 117, in angmom_components_sph
    "momentum {0}".format(self.angmom)
ValueError: Given convention does not support spherical contractions for the angular momentum 0

If I more or less understood the code correctly in gbasis, this happens because when the code in kinetic_energy.py tries to construct the contractions array with KineticEnergyIntegral(basis).construct_array_mix(coord_type) (coord_type = coord in the reproduce code and coord being the list of cartesian/spherical values for each shell) it goes to the base class module in base_two_symm.py, to use construct_array_mix() function defined there and in this case it builds the transformation variable for all the basis functions

In contrast to this behaviour for BaseTwoIndexSymmetric, for BaseOneIndex in base_one.py this sorting is done before building the transform variable and avoiding this mistake.

To Reproduce

#! usr/bin/env python
from iodata import load_one

from gbasis.wrappers import from_iodata
from gbasis.integrals.kinetic_energy import kinetic_energy_integral, KineticEnergyIntegral
import numpy as np

mol = load_one("neon.fchk")
basis, coord = from_iodata(mol)

kinetic_integrals = kinetic_energy_integral(basis, coord_type = coord)

Expected behaviour

Screenshots

System Information:

• OS: Ubuntu 18.04.5 LTS
• Python version: 3.6.10
• NumPy version: 1.18.1
• SciPy version: 1.4.1

Additional context

There's a difference in some positive and negative signs in the outputs.

tests/test_libcint.py::test_integral[STO-6G-Be-Cartesian-AngularMomentum] 
GBASIS
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00-1.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+1.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00-1.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00-1.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+1.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00-1.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
CBASIS
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+1.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00-1.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00-1.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00+1.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+1.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00-1.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j, 0.000000000e+00+0.000000000e+00j
FAILED

Originally posted by @msricher in #137 (comment)

Attached image is what the Gbasis notes say angular momentum is, which matches up with the code as far as I can tell.

Libcint doesn't give as much documentation, but it says it's r cross Nabla; does this exclusively imply the formula in the image attached, or could there be a difference of convention? I feel like a difference of convention could explain why there's a factor of -1 difference in some elements, but not others.

Is there anyway I can evaluate this array numerically using Grid, to verify whether GBasis or Libcint is correct?

Check: theochem/iodata#42

Hello Everyone,

I will try to summarize my finding as far as I could.
is_overlap_included: Return True or False to indicate that these two contractions' overlaps cannot be ignored. It comes in response to whether overlap screening is required by the basis. if required, we verify the possibility by checking whether the shells are far away, and can be screened or not.

If the first contraction overlap screening is not required, we return True. Otherwise, we calculate the distance between the centers of the two shells, and then the cutoff distance (or distance limit) which depends on the shells' minimum exponents and tolerance.
If the distance between the two shells is greater than the cutoff, we return False, otherwise True.

Now, I have two questions in mind regarding overlap and asymmetric overlap in particular.

1. For Asymmetric overlap /gbasis/integrals/overlap.py Line 185:
   is checking for only the first contraction's ovr_screen sufficient?

if not contractions_one.ovr_screen:
   return True

While it does not pose any problems with Symmetric Overlap, but I don't know if the same logic can be applied with Asymmetric Overlap between two contractions from two different basis set. As one of them can be False and the other True?

2. in overlap_integral, I need to ask what does this format in Line 145

coord_type = [ct for ct in [shell.coord_type for shell in basis]]

serve rather than just using

coord_type = [shell.coord_type for shell in basis]

Describe
Hello folks,
Thanks for developing gbasis! :)
I am using gbasis in a package I am writing. I want to make part of my CI, but the pypi installation has a version from 2017. Would it be possible to get it updated to the most recent commit? Unless there has not been any relevant changes in the meantime.

Thanks a lot!

Arbitrary order multipoles and derivative (a special case is kinetic energy density)

Describe the bug
It is not clear what are the objects returned by the from_iodata interface. If the returned lists contains the contracted shells, then the overlap and other integrals should return square matrix with the number of rows/columns
equal to the returned value from from_iodata.

To Reproduce

from iodata import load_one
from gbasis.wrappers import from_iodata
import numpy as np
from gbasis.integrals.overlap import overlap_integral

mol = load_one("./data/c2h4_q0.fchk")
ao_basis = from_iodata(mol)

print("Number of contracted basis functions: ", len(ao_basis))

# compute overlap of atomic orbitals (AO)
int1e_overlap_ao = overlap_integral(ao_basis)
print("Overlap matrix (S) of atomic orbitals (AO) shape: ", int1e_overlap_ao.shape)

returns:

Number of contracted basis functions:  64
Overlap matrix (S) of atomic orbitals (AO) shape:  (184, 184)

Expected behaviour
The overlap matrix should be of dimension NxN where N is the number of contracted shells that from_iodata returns

Screenshots

System Information:

• OS:
• Python version:
• NumPy version:
• SciPy version:

Additional context

When #106 was pulled in, I didn't notice that we should update the notes.pdf file as well, because the wrappers return both the basis and, newly, coord_types. I'll try to update the docs soon.

Since functions inside gbasis.wrappers depend on extra dependencies, we would need to install them in order to test for coverage. There are a couple of options I think:

1. Make sure all of the test environments install the extra dependencies. This would ensure that every possible feature of the code gets tested fully. However, these are optional dependencies and it would be nice to distinguish the testing environment for the core features and the optional features. And we don't want the tests to depend on another (optional) package.
2. Make a separate test environment for testing coverage. Maybe this is the way to go.

Is there a specific integral/formula that you would like implemented?
Current Gbasis code performs arbitrary order derivatives by using specific formula involving Hermite polynomials. Based on later talks, it would be interesting to prototype an test a version of Gbasis that could use automatic differentiation. This would give access to nuclear derivatives as well, which is a missing feature.

Is there a change to the code or algorithm you would like to see?
Right now to compute a specific order derivative you would use funct evaluate_deriv_reduced_density_matrix and specify the orders for first and second grid points. This ultimately calls _eval_deriv_contraction. The goal would be to use Automatic differentiation to avoid the explicit use of this function.

Describe the bug

While I'm trying to set up github actions to use pytest. I've came across an error on ubuntu systems for only python 3.9, python 3.10, python=3.11 and for windows system, python=3.9,3.10,3.11. It passes on my own system and python=3.7

The test is test_compute_differential_operator_integrals_multiarray in test_diff_operator_int.

It seems that they both return nans. Note this requires to fix the test as done in #151
Python=3.7, I'll get the following value: order, angmom_a, angmom_b:
[3 0 0] [2 1 0] [2 0 0] -0.0007958223786576885
Python-3.10, I'm getting this array to be all infinity and nans

The error comes from teh normalization factor calculated in

which comes from the factorial2

To Reproduce

conda create -n py310 python=3.10
conda activate py310
pip install .
pytest -v .

Expected behaviour

Screenshots

See the run: https://github.com/theochem/gbasis/actions/runs/7491475696/job/20392775547?pr=151

         for i, single_order in enumerate(orders_diff):
            for j, angmom_a in enumerate(angmoms_a):
                for k, angmom_b in enumerate(angmoms_b):
>                   assert np.allclose(
                        _compute_differential_operator_integrals(
                            np.array([single_order]),
                            coord_a,
                            np.array([angmom_a]),
                            exps_a,
                            coeffs_a,
                            norm_a,
                            coord_b,
                            np.array([angmom_b]),
                            exps_b,
                            coeffs_b,
                            norm_b,
                        )[0, 0, j, 0, k],
                        test[i, 0, j, 0, k]
                    )
E                   assert False
E                    +  where False = <function allclose at 0x7fe17171d4b0>(nan, nan)
E                    +    where <function allclose at 0x7fe17171d4b0> = np.allclose

System Information:

• OS:
• Python version:
• NumPy version:
• SciPy version:

Additional context

Is your feature request related to a problem? Please describe.

We need to be very clear about units. If I'm not mistaken, everything in GBasis is in atomic units. While we are specific about atomic coordinates being in a.u., we should be specific about point-sets passed in, basis functions, and the properties/integrals too, as it isn't entirely obvious.  

One could document the most important functions explicitly and put a note at the bottom of the main web page/Readme to cover cases where we don't explicitly indicate units.

Description

Use screening based on the overlap between basis functions to improve performance.

📚 Package Description and Impact

For large molecules, textbook expressions for quantities expanded in Gaussian basis functions (e.g., the electron density) or integrals based on Gaussian basis functions (e.g., the kinetic-energy integral) typically include many negligible terms. By screening out these terms, and only evaluating terms that are nonnegligible, the performance of GBasis can be greatly enhanced.

👷 What will you do?

In GBasis we provide a utility for screening these terms using their overlap, is_overlap_included. When this expression is small, one can also neglect other one-electron integrals. A generalization of this approach allows fast evaluation of spatial quantities and 2-electron integrals. Your main goal would be to screen 1-electron integrals and the evaluation of quantities at (grid) points using overlap screening and its generalization.

🏁 Expected Outcomes

1. Adapt is_overlap_included to screen other one-electron integrals.
2. Extend is_overlap_included to three functions, which allows screening spatial evaluations.
3. Write tests to ensure correctness and assess performance.
4. 🏆 An ambitious stretch goal is to implement screening of 2-electron integrals.

🙋 Mentors

Could I get admin access to this repo?

The module iodata is supported, technically, but not very well. iodata for which the interface was created is not available on neither pypi nor conda. It was designed for the latest iodata in the master branch of the repository https://github.com/theochem/iodata. There are some major API breaking changes from the previous version which makes the current wrapper not compatible with the ones in pypi or conda.

The setup.py references that the iodata must be greater than 0.1.7 (current version in pypi and conda), but due to reasons the latest iodata has not been pushed or versioned yet.

The fix, at the moment, is to clone the repository and install it manually

https://github.com/theochem/iodata
cd iodata
pip install --user -e .

and the pip install gbasis[iodata] is broken until iodata is available (This command wouldn't work anyways because the current version of gbasis isn't up on pypi or conda, but that'll happen soon).

See PR #69 for some discussion

Describe the bug

While navigating through the notes.pdf and later at parser.py, I think I found some unintended words in place.

To Reproduce

1. /docs/notes.pdf Section 6.2 Multipole moment intergral:
   Instead of: "To compute the overlap of a set of molecular orbitals"
   It should be: "To compute the moment of a set of molecular orbitals"
2. /gbasis/parser.py Line 217:
   Instead of: "Tolerance must be provided as True or False."
   It should be: "Overlap must be provided as True or False."

Expected behaviour

Minimal Effect on users and developers' experience.

Screenshots

2.

System Information:

• OS:
• Python version:
• NumPy version:
• SciPy version:

Additional context

compute_external_potential(shell1, center1, shell2, center2, conventions, point_charges, point_coords)

Compute the potential due to a product of 2 shells of Gaussians with an array of point charges and their centers. Gives the matrix of external potential integrals

• generalize to compute_external_potential_type() to compute this for different types of potentials (not necessarily point charges). The only change is that a different Boys function is used, so there should be a low-level function of the form
  compute_extpot(shell1, center1, shell2, center2, conventions, Boys_function)
  (or something like that).

Algorithm could be McMurchie-Davidson or Obara-Saika but probably Obara-Saika is a bit easier?

Apparently, there is an old website in the web: https://gbasis.readthedocs.io/en/latest/index.html
@PaulWAyers @FarnazH can we somehow crush it?
I think here is an instruction on how to do that https://groups.google.com/g/read-the-docs/c/xJyqg7QSX9Y

Gbasis Documentation website

This issue is to request the creation of a Gbasis Documentation website similar to https://grid.qcdevs.org/. The most complete set of notes about Gbasis functionality and what is already implemented can be found within the repo here. Example tutorials are being developed and updates kept track in issue #133