backflow issue when C=0

[Vladimir_Konjkov]

Dear Casino developer!

Quite randomly, I generate test data for backflow calculations I found a combination of parameters that issues SIGSEGV (files in the attachments).

As far as I understand, the error occurs when executing SUBROUTINE construct_C or somewhere nearby, I specially chose C=0, to follow the corresponding branches in the code, and for all combinations N_eN and N_ee less than five with the exception of (5, 4) everything is fine.

Unfortunately, I myself cannot figure out what the error is.

Best, Vladimir.

PS:

may be in highloghted line of code
if(k>1.and.m<N_phi_eN(set))C(n,p+nphi-1*inc_k+1*inc_m)=-real(m,dp)
m should be compared with N_phi_ee(set) because m loop is between do m=1,N_phi_ee(set)?

Code: Select all

! Do irrotational bit of the constraint matrix.
 if(set_phi_is_ae(set))then
  offset_irrot=offset_eN3b+NN_2
 else
  offset_irrot=offset_eN1a+NN_1
 endif

 if(set_phi_is_irrot(set))then
  p=0 ; n=offset_irrot
  inc_k=1
  inc_l=inc_k*N_phi_eN(set)
  inc_m=inc_l*N_phi_eN(set)
  nphi=inc_m*N_phi_ee(set)
  do m=1,N_phi_ee(set)
   do l=1,N_phi_eN(set)
    do k=1,N_phi_eN(set)
     p=p+1 ; n=n+1
     if(C_trunc>0)then
      if(m>1)then
       C(n,p-1*inc_m)=real(C_trunc+k-1,dp)
       if(k<N_phi_eN(set))C(n,p+1*inc_k-1*inc_m)=-Lc*real(k,dp)
      endif
      if(m<N_phi_ee(set))then
       if(k>2)C(n,p+nphi-2*inc_k+1*inc_m)=-real(m,dp)
       if(k>1)C(n,p+nphi-1*inc_k+1*inc_m)=Lc*real(m,dp)
      endif
     else ! C=0
      if(m>1.and.k<N_phi_eN(set))C(n,p+1*inc_k-1*inc_m)=real(k,dp)
      if(k>1.and.m<N_phi_eN(set))C(n,p+nphi-1*inc_k+1*inc_m)=-real(m,dp)
     endif ! C>0 or not
    enddo ! k
   enddo ! l
  enddo ! m

[Neil Drummond]

Dear Vladimir,

Thanks very much for this. As you suggest, replacing

Code: Select all

if(k>1.and.m<N_phi_eN(set))C(n,p+nphi-1*inc_k+1*inc_m)=-real(m,dp)

with

Code: Select all

if(k>1.and.m<N_phi_ee(set))C(n,p+nphi-1*inc_k+1*inc_m)=-real(m,dp)

in pbackflow.f90 should be correct. (I've checked with Pablo.)

Sorry about the problem.

Best wishes,

Neil.

[Mike Towler]

Now fixed in the published distribution.

Thanks Vladimir!

Cheers,
Mike

[Vladimir_Konjkov]

Hello Mike.

Nice to read you.

I am writing a program in python, numpy and numba that reads Casino input files and performs calculations given in INPUT and now I came to BACKFLOW implementation...everything is so terribly complicated there. It seems to me that the ETA term is calculated incorrectly, it looks like displacement of first (o last) electrons is always zero. May be in the formula the summation is performed on i>j (or i<j) only therefore, for the first (last) electron there are no terms.

Screenshot_20210415_113922.png (4.71 KiB) Viewed 138085 times

If you tell me how to print out backflow displacements (for ETA-term) I'll be very grateful or there might be an error in my code then it's even worse...

Best, Vladimir.

[Pablo_Lopez_Rios]

Does the output of `plot_backflow` help with your debugging? (See `casinohelp plot_backflow` for info.)

[Vladimir_Konjkov]

Does the output of `plot_backflow` help with your debugging? (See `casinohelp plot_backflow` for info.)

Thank you Pablo, I'll try to check right now. My doubts arose from my implementation of eta-term as:

Code: Select all

    def eta_term(self, e_vectors, e_powers):
        """
        :param e_vectors: e-e vectors
        :param e_powers: powers of e-e distances
        :return: displacements of electrons - array(nelec, 3)
        """
        res = np.zeros((self.neu + self.ned, 3))
        if not self.eta_cutoff.any():
            return res

        C = self.trunc
        parameters = self.eta_parameters
        for i in range(1, self.neu + self.ned):
            for j in range(i):
                r_vec = e_vectors[i, j]
                r = e_powers[i, j, 1]
                eta_set = (int(i >= self.neu) + int(j >= self.neu)) % parameters.shape[1]
                # I don't think it works fast if NO SPIN-DEP
                L = self.eta_cutoff[eta_set] or self.eta_cutoff[0]
                if r < L:
                    poly = 0
                    for k in range(parameters.shape[0]):
                        poly += parameters[k, eta_set] * e_powers[i, j, k]
                    res[i] += (1 - r/L) ** C * poly * r_vec
                    # res[j] -= ... commented line
        return res

this implementation gives the same energy as CASINO, calculation of gradient and Laplacian I can make numerically or analiticaly - it doesn't change the result.

[Vladimir_Konjkov]

It was a terrible mistake in my code in the function calculating smooth cutoffs around AE atoms.
Thanks everyone for the help.

[Pablo_Lopez_Rios]

Glad you found the problem!

Best,
Pablo