GFN-FF runs into a segfault

[BuggsBuster]

Describe the bug
A clear and concise description of what the bug is.

Using the latest stable XTB version 6.6.1 I run into a segfault when performing a single point calculation on our cluster using GFN-FF for either of the two structures attached.
Running the same jobs on my personal laptop, I just get a warning instead of the segfault:

########################################################################
[WARNING] Runtime exception occurred
-1- gfnff_setup: Could not read topology file.
########################################################################

When I rerun these job on my laptop in the same folder, XTB is able to read the topology file produced in the first iteration.
In the attached files, I piped stdout and stderr into the respective res.out files.

To Reproduce
Steps to reproduce the behaviour:

• I use the precompiled XTB binary that I downloaded from this Github repo.
• The command I used: xtb inp.xyz --gfnff --verbose

1. happens with input (include input files)
2. start xtb with (all the options here)
3. run xtb with your options and the --verbose flag
4. output showing the error

Please provide all input and output file such that we confirm your report.

Expected behaviour
A clear and concise description of what you expected to happen.

A single point calculation without any segfaults nor any warnings.

Additional context
Add any other context about the problem here.

Based on a rough analysis, I assume the bug originates from the following code snippet obtained from the specialTorsList procedure contained in gfnff_ini.f90:

integer :: i,j,k,ii,jj,kk,ll,idx
  logical :: iiok, llok
  idx=0
  do i=1, mol%n
    ! carbon with two neighbors bonded to other carbon* with two neighbors
    if (mol%at(i).eq.6.and.topo%nb(20,i).eq.2) then
      do j=1, 2
        nbi=topo%nb(j,i)
        if (mol%at(nbi).eq.6.and.topo%nb(20,nbi).eq.2) then  ! *other carbon
          ! check carbon triple bond distance
          if (NORM2(mol%xyz(1:3,i)-mol%xyz(1:3,nbi)).le.2.37) then
            ! at this point we know that i and nbi are carbons bonded through triple bond
            ! check C2 and C3
            do k=1, 2  ! C2 is other nb of Ci
              if (topo%nb(k,i).ne.nbi.and.mol%at(topo%nb(k,i)).eq.6) then
               print *, "Setting jj"
                jj=topo%nb(k,i)
              endif
            enddo
            do k=1, 2  ! C3 is other nb of Cnbi
              if (topo%nb(k,nbi).ne.i.and.mol%at(topo%nb(k,nbi)).eq.6) then
                kk=topo%nb(k,nbi)
                print *, "Setting kk"
              endif
            enddo
            ! check C1 through C4 are sp2 carbon
            if (topo%hyb(jj).eq.2.and.topo%hyb(kk).eq.2 &
            &   .and.mol%at(jj).eq.6.and.mol%at(kk).eq.6) then

The variables jj and kk are declared, but not initialized at the beginning. I assume that the two if statements, that are meant to assign values to jj and kk, are never entered for the given examples, leaving the two variables uninitialized. Thus, in the last if statement, jj and kk end up with values beyond the size of the arrays topo%hyb and mol%at, resulting in the segfault.

Can you confirm the issue?

xtb661_issue.zip

[MtoLStoN]

Hi @BuggsBuster,
thank you for reporting this Issue and doing a rough bug analysis. While I cannot replicate this bug at the moment, I believe you are right in your analysis.

xtb/src/gfnff/gfnff_ini.f90

Lines 1924 to 1933 in 3f25dfd

	do k=1, 2 ! C2 is other nb of Ci
	if (topo%nb(k,i).ne.nbi.and.mol%at(topo%nb(k,i)).eq.6) then
	jj=topo%nb(k,i)
	endif
	enddo
	do k=1, 2 ! C3 is other nb of Cnbi
	if (topo%nb(k,nbi).ne.i.and.mol%at(topo%nb(k,nbi)).eq.6) then
	kk=topo%nb(k,nbi)
	endif
	enddo

Looking at the code, in your example of an Acetylene molecule, the if statement would never trigger, leading to uninitialized jj and kk variables, which is most likely the error source, when they are used here:

xtb/src/gfnff/gfnff_ini.f90

Lines 1935 to 1936 in 3f25dfd

	if (topo%hyb(jj).eq.2.and.topo%hyb(kk).eq.2 &
	& .and.mol%at(jj).eq.6.and.mol%at(kk).eq.6) then

We should probably just break the subroutine, when C2 and C3 can not be found, but I will let the Force Field professionals decide ;). @Thomas3R, can you take a look at this?

[Thomas3R]

Thanks for the detailed bug report! Though I could not replicate the bug, I agree with Marcel and think merge #854 should resolve this issue. Please open a new issue if your problem persists.

[BuggsBuster]

@Thomas3R @MtoLStoN
I can confirm that this fixes the bug. Thanks a lot!