Improvements to continuous space

[PhilipVinc]

Things that come to my mind, by reading @gpescia 's tutorial:

• Constructor of Particle. To avoid having to specify the annoying jnp.inf we could have some default arguments:

# This is the current constructor
Particle(N=N, L=(jnp.inf, jnp.inf, jnp.inf), pbc=False)

# I propose adding the following, which is equivalent to the above
Particle(N=N, D=3)

# Specifying pbc to True will error
Particle(N=N, D=3, pbc=True)
# ValueError: "If you want to specify periodic boundary conditions you should have finite size for the bounding box by using the L argument.

# Specifying both D and L errors as it is redundant
Particle(N=N, L=(jnp.inf, jnp.inf, jnp.inf), D=3, pbc=False)
# ValueError: "You should not specify both D and L, as the dimensionality is inferred by the number of dimensions in the bounding box

WDYT? Internally, there will be still those jnp.infty (I know people don't like them. I'm happy if a pr gets rid of them but dunno how)

Second thing:
It seems to me one often has to define a minimal image distance to use inter particle potentials.
I propose adding a distance function to the Particle Hilbert space that returns the distance according to the minimum distance convention (or not, if the bounding box is infinite)

hi = Particle(N=N, D=3, pbc=True)
hi.distance(hi.random_states(...))

so that we can use that a bit everywhere instead of re-coding it every time, to make life easier for users.

Also we could add a new class TwoBodyPotential that automatically evaluates the potential among particles, to make it easier to implement this kind of potential.
It would be just a wrapper around Potential. We could then add some standard potentials ourselves (Aziz, Coulomb`...)

Mainly I would like feedback from the people that use the continuous space stuff to know if this makes sense. For example, do you always use minimum image or do you sometimes need some different, weird distance functions? Would a TwoBodyPotential function be helpful? Is there some pain point that is missing?

[gcarleo]

I agree with both proposals, they would improve code readability

[gcarleo]

Concerning long-range potentials, they should be implemented as well, I agree, also taking into account that computing the sum over periodic images can be tricky, so we want to have this implemented once for all!

[gpescia]

I also agree with both. Though do you think it makes sense to put the second proposal into nk.utils? I‘m not sure it belongs to hilbert and like this, other functionality can be added as well such as Ewald or radial distribution etc. (or maybe these actually belong in operator…)

[gcarleo]

The main issue is that in continuous space we are not using Graph or similar as we do on the lattice, so this means that all the geometry info (pbc etc) that should be contained in the geometry is instead contained only in hilbert. In theory, I agree with you that the distance shouldn't be in hilbert, it should be in some Geometry class that you pass to Hilbert. Geometry could also take into account specific things such as Ewald, returning reciprocal vectors etc.

For the observables, I am not sure... I think we should just consider them as operators?

[gcarleo]

so a more complete proposal is something like:

geo=Cell(d=3, L=1, pbc=True)

hi = Particle(N=N, geometry=geo)
... 
geo.distance(R1,R2)

notice that I would also add to Particle the possibility of specifying other quantum numbers for each particle, for example the spin:

geo_fermions=Cell(d=3, L=1, pbc=True)
geo_bosons=FreeSpace(d=2)

#mixture of fermions and bosons where bosons are confined in 2D
hi = Particle(N=N1, spin=1/2, geometry=geo_fermions) + Particle(N=N2, spin=0, geometry=geo_bosons) 
... 
geo.distance(R1,R2)

in theory, this allows to use Particle also for the lattice, and have a unified interface:

lattice=Square(L=3)

hi = Particle(N=N,geometry=lattice,spin=1/2)

#equivalent to 
hi = SpinOrbitalFermions(n_orbitals=lattice.n_nodes,n_fermions=N)
 
#similarly for Particle(N=N,geometry=lattice,spin=0) calling Fock(...) internally 

[gcarleo]

@dalin27 you are our best hope here

[gcarleo]

I think implementing Geometry. Cell and Geometry.FreeSpace shouldn't be too complicated to start with , this will give us also an idea on the abstract interface of Geometry objects

[gcarleo]

Hierarchy wise I think that abstract Graph should be inherited from the abstract Geometry or, we should have the constructors above accept both Graphs and Geometries (a bit less elegant), what do you think @PhilipVinc ?

[PhilipVinc]

I think they have not much to do with each other

[gcarleo]

so we make Particle accept both Geometry and Graph? makes sense, it's just that a Graph however is a special form of Geometry ...

[dalin27]

Regarding @PhilipVinc's second point about adding a distance method to the Particle Hilbert space, should the method return a distance function (min_imag_conv) like this:

def distance(self):
    
    pbc = np.array(self._pbc)
    L = np.array(self._extent)

    if np.all(pbc):

        def min_imag_conv(x: Array):
            """Compute distances between particles using the minimum image convention"""
            x = x.reshape(self._N, self._D)  # not batched!
            distances = x[None,:,:] - x[:,None,:]
            distances = distances[jnp.triu_indices(self._N, 1)]  # upper triangular part of the matrix
            distances = jnp.remainder(distances + L/2.0, L) - L/2.0  # shape=(N(N-1)/2,D)
            distances = jnp.linalg.norm(distances, axis=-1)
            return distances
        
        return min_imag_conv

    elif np.all(pbc == False):
        raise ValueError(
            " Not yet implemented. "
        )
    else:
        raise ValueError(
            " Not yet implemented. "
        )

or instead return the evaluated distance values? The code above also requires adding the line self._D = len(L) if D is None else D in the __init__ method

[PhilipVinc]

I think it should be a function

class AbstractGeometry:
    def min_imag_convention(self, x: Array) -> Array:
       ....
       return min image 

An important thing that we should keep in mind is how to handle the composed space of two different geometries. (Aka, imagine you have two sublattices each with a different set of fermions, for example..

Also, Is there anything else that would be useful to compute something?

Also, the AbstractGeometry should define __hash__ and __eq__ so that it plays nice with jax, as Hilbert spaces do.