From 7eaab9aa7d4ba57fd2269ce7b09b581adc849f24 Mon Sep 17 00:00:00 2001
From: Marcel Stahn <70513124+MtoLStoN@users.noreply.github.com>
Date: Tue, 1 Aug 2023 13:44:22 +0200
Subject: [PATCH] Fixes GFN-FF bug for specific pi systems (#836)

* Fixes GFN-FF bug for specific pi systems

Signed-off-by: MtoLStoN <70513124+MtoLStoN@users.noreply.github.com>

---------

Signed-off-by: MtoLStoN <70513124+MtoLStoN@users.noreply.github.com>
Signed-off-by: Johannes Gorges <58849467+gorges97@users.noreply.github.com>
---
 src/gfnff/gfnff_ini.f90  | 30 ++++++++++++++++++++----------
 test/unit/test_gfnff.f90 |  4 ++--
 2 files changed, 22 insertions(+), 12 deletions(-)

diff --git a/src/gfnff/gfnff_ini.f90 b/src/gfnff/gfnff_ini.f90
index 41f605f67..d88b3e6b0 100644
--- a/src/gfnff/gfnff_ini.f90
+++ b/src/gfnff/gfnff_ini.f90
@@ -918,16 +918,26 @@ integer function itabrow6(i)
       if(pr)then
          write(env%unit,'(''Hueckel system :'',i3,'' charge : '',i3,'' ndim/Nel :'',2i5, &
      &         3x, ''eps(HOMO/LUMO)'',2f12.6)')pis,ipis(pis),npi,nelpi,pisip(pis),pisea(pis)
-         if(pisip(pis).gt.0.40) then
-            write(env%unit,*)'WARNING: probably wrong pi occupation. Second attempt with Nel=Nel-1!'
-            do i=1,mol%n
-               if(piadr4(i).ne.0) write(env%unit,*) 'at,nb,topo%hyb,Npiel:', i,mol%sym(i),topo%nb(20,i),topo%hyb(i),piel(i)
-            enddo
-            nelpi=nelpi-1
-            Api = Apisave
-            call gfnffqmsolve(.false.,Api,S,.false.,300.0d0,npi,0,nelpi,dum,occ,eps)  !diagonalize
-            call PREIG(6,occ,1.0d0,eps,1,npi)
-         endif
+      end if
+      if(pisip(pis).gt.0.40) then
+         write(env%unit,'(a,i0,a)')'WARNING: probably wrong pi occupation for system ',pis,'. Second attempt with Nel=Nel-1!'
+         do i=1,mol%n
+            if(piadr4(i).ne.0) write(env%unit,*) 'at,nb,topo%hyb,Npiel:', i,mol%sym(i),topo%nb(20,i),topo%hyb(i),piel(i)
+         enddo
+         nelpi=nelpi-1
+         Api = Apisave
+         call gfnffqmsolve(.false.,Api,S,.false.,4000.0d0,npi,0,nelpi,dum,occ,eps)  !diagonalize
+         call PREIG(6,occ,1.0d0,eps,1,npi)
+         do i=1,npi  ! save IP/EA
+            if(occ(i).gt.0.5) then
+               pisip(pis)=eps(i)   ! IP
+               if(i+1.lt.npi)pisea(pis)=eps(i+1) ! EA
+            endif
+         enddo
+      if(pr)then
+         write(env%unit,'(''Hueckel system :'',i3,'' charge : '',i3,'' ndim/Nel :'',2i5, &
+     &         3x, ''eps(HOMO/LUMO)'',2f12.6)')pis,ipis(pis),npi,nelpi,pisip(pis),pisea(pis)
+      end if
       endif
 ! save BO
       do i=1,topo%nbond
diff --git a/test/unit/test_gfnff.f90 b/test/unit/test_gfnff.f90
index 524c12772..bea7b1daf 100644
--- a/test/unit/test_gfnff.f90
+++ b/test/unit/test_gfnff.f90
@@ -487,10 +487,10 @@ subroutine test_gfnff_scaleup(error)
    character(len=*), parameter :: solvents(5) = [character(len=20) ::&
       & "h2o", "acetonitrile", "toluene", "ether", "dmso"]
    real(wp), parameter :: ref_energies(5) = &
-      &[-4.6919926039901_wp, -8.8807760138817_wp, -13.311107073100_wp, &
+      &[-4.6919926039901_wp, -8.8807760138817_wp, -13.40715596616603_wp, &
       & -13.822994859730_wp, -20.577952329212_wp]
    real(wp), parameter :: ref_gnorms(5) = &
-      &[0.05947676640487_wp, 0.09522104624089_wp, 0.16195989065950_wp, &
+      &[0.05947676640487_wp, 0.09522104624089_wp, 0.1753769495334539_wp, &
       & 0.12496592222660_wp, 0.19366599743810_wp]
 
    call init(env)