WS construction in Make_lattice may generate wrong number of real-space points near boundaries

[ALF-Import-Bot]

In GitLab by @pangaopei on Jan 22, 2026, 02:31 UTC:

Recently, Disha encountered a problem very similar to what I reported before in issue #294.

When the lattice type is triangular, honeycomb, or kagome, the subroutine

Libraries/Modules/lattices_v3_mod.F90
  subroutine Make_lattice

(approximately lines 176–240)

uses a Wigner–Seitz (WS) construction to generate the real-space lattice. However, in many cases, this procedure does not correctly handle points on or near the boundary of the WS cell, which leads to an incorrect number of lattice points being generated.

In particular, the number of generated points can be larger than the expected L1*L2.

Example

Using the following input:

&VAR_lattice               !! Parameters defining the specific lattice and base model 
  L1           = 10            ! Length in direction a_1
  L2           = 4             ! Length in direction a_2
  Lattice_type = "triangular"  ! 
/

the code reports an error:

Make_lattice: Error           40          41
lattices_v3_mod.F90 (274):
Generic error
Terminating program

Where the problem seems to come from

The issue seems to originate from the WS construction part in lines ~216–240 of Libraries/Modules/lattices_v3_mod.F90, namely:

           nc = 0
           do i1 = -L,L
              do i2 = -L,L
                 x_p  = dble(i1)*a1_p + dble(i2)*a2_p
                 L_f = 1
                 do i = 1,4
                    if (i.eq.1) a_p =  L2_p
                    if (i.eq.2) a_p =  L1_p
                    if (i.eq.3) a_p =  L2_p - L1_p
                    if (i.eq.4) a_p =  L2_p + L1_p
                    if  (  Iscalar( x_p, a_p )    .le.  xnorm(a_p)**2/2.d0 + Zero   .and.   &
                         & Iscalar( x_p, a_p )    .ge. -xnorm(a_p)**2/2.d0 + Zero    ) then
                       L_f = L_f * 1
                    else
                       L_f = 0
                    endif
                 enddo
                 if (L_f .eq. 1) then
                    nc = nc + 1
                    Latt%list(nc,1) = i1
                    Latt%list(nc,2) = i2
                    Latt%invlist(i1, i2 ) = nc
                 endif
              enddo
           enddo

and also Line 256-257.

More generally

There are many combinations of (L1, L2) that trigger this problem (I have not exhaustively tested them, but a large fraction of combinations seem to fail).

For example, combinations like:
(6,2),(8,2),(10,2)...(10,4),(12,4),(14,4)...(2,6)(14,6)(16,6)(18,6)...(2,8)(18,8)(20,8)...
can also lead to wrong behavior.

(To Disha: I think your case is (20,8)? )

Likely cause

The current WS construction in Make_lattice seems to:

• Include boundary points in an inconsistent way
• Not treat floating-point tolerances (eps) carefully
• As a result, some boundary points are double-counted or incorrectly included

This leads to:

The number of generated real-space lattice points being larger than L1 * L2.

Suggestion

We probably need to:

• Revisit the WS-cell inclusion test
• And treat the epsilon tolerance in a more robust and consistent way to avoid double counting or missing of boundary points.

Migrated from GitLab: https://git.physik.uni-wuerzburg.de/ALF/ALF/-/issues/333

[ALF-Import-Bot]

In GitLab by @pangaopei on Jan 22, 2026, 02:44 UTC:

marked this issue as related to #294

[ALF-Import-Bot]

In GitLab by @fassaad on Jan 22, 2026, 07:44 UTC:

created branch 333-ws-construction-in-make_lattice-may-generate-wrong-number-of-real-space-points-near-boundaries to address this issue

[ALF-Import-Bot]

In GitLab by @DishaHou on Jan 22, 2026, 11:53 UTC:

Hi Gaopei, thanks for reporting this!

For my case, I encountered the issue with input:

&VAR_lattice               !! Parameters defining the specific lattice and base model
L1           = 20          ! Length in direction a_1
L2           = 8          ! Length in direction a_2
Lattice_type = "Triangular"     ! Sets a_1 = (1,0), a_2=(0,1), Norb=1, N_coord=2
Model        = "SSH"    ! Sets the Hubbard model, to be specified in &VAR_Hubbard
/

This leads to the error:

 Make_lattice: Error          160         162
 lattices_v3_mod.F90 (274):
 Generic error                 

which is from the mismatch between nc (number of points in 1st BZ) and Latt%N (number of lattice points in real space) in code file lattices_v3_mod.F90:

272            If (nc.ne.Latt%N) Then
273               write(error_unit,*) 'Make_lattice: Error ', nc, Latt%N
274               Call Terminate_on_error(ERROR_GENERIC,__FILE__,__LINE__)
275            endif

Here, nc is correct, but the code counts Latt%N wrong, which comes from lines 176-208 when it counts the number of lattice points in real-space:

176            ! Count the number of lattice points.
177            L      =   abs(nint ( Iscalar(Latt%BZ1_p,L1_p) / (2.d0*pi) ))
178            L1     =   abs(nint ( Iscalar(Latt%BZ2_p,L1_p) / (2.d0*pi) ))
179            if (L1 .gt. L) L = L1
180            L1     =   abs(nint ( Iscalar(Latt%BZ1_p,L2_p) / (2.d0*pi) ))
181            if (L1 .gt. L) L = L1
182            L1     =   abs(nint ( Iscalar(Latt%BZ2_p,L2_p) / (2.d0*pi) ))
183            if (L1 .gt. L) L = L1
184            nc = 0
185            do i1 = -L,L
186               do i2 = -L,L
187                  x_p  = dble(i1)*a1_p + dble(i2)*a2_p
188                  L_f = 1
189                  do i = 1,4
190                     if (i.eq.1) a_p =  L2_p
191                     if (i.eq.2) a_p =  L1_p
192                     if (i.eq.3) a_p =  L2_p - L1_p
193                     if (i.eq.4) a_p =  L2_p + L1_p
194                     if  (  Iscalar(x_p, a_p)  .le.  xnorm(a_p)**2/2.d0 + Zero   .and.   &
195                          & Iscalar(x_p, a_p)  .ge. -xnorm(a_p)**2/2.d0 + Zero    ) then
196                        L_f = L_f * 1
197                     else
198                        L_f = 0
199                     endif
200                  enddo
201                  if (L_f .eq. 1) then
202                     nc = nc + 1
203                  endif
204               enddo
205            enddo
206            LQ = nc
207            Latt%Ns = LQ
208            Latt%N  = LQ

[ALF-Import-Bot]

In GitLab by @fassaad on Jan 24, 2026, 05:31 UTC:

Dear all, the branch attached to this issue solves the above problems. However this is still in a certain sense a quick fix, since it changes the form of the Brillouin zone. There are many gauge equivalent ways of defining the BZ, and in the master branch we had opted for the Wigner Seitz cell construction. In this branch the BZ is chosen to be a parallelogram. This implies that the set of k-points in the analysis files will differ between this branch and the master branch. As soon as I have more time, I will try to amend this.