4th order correction D4V breaks symmetry

[ionerrea]

We found a system in which the 4th order correction breaks the symmetry, while the old fortran code does not. There must be a small bug somewhere in the Python version. The funny thing is that the 3rd order does not break it, so it seems to be something that only happens in the 4th order correction in the Python version.

[FlemkMeserath]

The symmetry breaking with the latest SSCHA version seem to appear only when calculating the fourth order corrections in a phonon supercell while in the unitary cell, or while just including the third order corrections, they seem to be correctly applied. Confronting the calculations with the old fortran code, the third order calculations appear to be the exact same, so the issue seem to be related just to the forth order of the newest version.

[mesonepigreco]

Maybe, it could be that I moved one line in a for loop when apply the d4 correction (FORTRAN code), because I thought it was a mistake that was causing a huge slow-down. The code gained a polynomial speedup, but it could be that it spoiled symmetries somehow.

It was in get_odd_straight.f90 (line 131 of the current version):

  do i = 1, n_mode
    maux(:,:) = v3(i,:,:)
    ja = 0 
    do mu = 1, n_mode 
      laux1 = l(mu,:)
      call dgemv ('N',ns,ns,1.0d0,maux,ns,laux1,1,0.0d0,lres1,1)    ! INSERTED HERE
      do nu = 1, n_mode         
        laux2 = l(nu,:)
        ja = ja + 1
!        call dgemv ('T',ns,1,1.0d0,lres1,ns,laux2,1,0.0d0,lsum,1)    ! REMOVED FROM HERE
        v1(i,ja) = dot_product(lres1,laux2) 
        v2(i,ja) = v1(i,ja) * g(mu,nu)
        v1(i,ja) = v1(i,ja) * 0.5d0
!        v2(i,ja) = lsum * g(mu,nu)
!        v1(i,ja) = lsum * 0.5d0
      end do
    end do
  end do

However this should be called only in the d3 bubble computation, not the D4, so I do not see how this can explain the issue.