Fix division by zero at geoopt (#1170)

* Avoid division by zero at first optimization step Signed-off-by: Igor S. Gerasimov <[email protected]> * Avoid division by zero during computing dPhi/dR for torsions Signed-off-by: Igor S. Gerasimov <[email protected]> --------- Signed-off-by: Igor S. Gerasimov <[email protected]>
grimme-lab · Feb 10, 2025 · 40c757e · 40c757e
1 parent 945d7e6
commit 40c757e
Show file tree

Hide file tree

Showing 2 changed files with 27 additions and 17 deletions.
diff --git a/src/constr.f90 b/src/constr.f90
@@ -613,7 +613,11 @@ Subroutine dphidrPBC(mode,nat,xyz,i,j,k,l,vTrR,vTrB,vTrC,phi,&
       dphidrj=0
       dphidrk=0
       dphidrl=0
-      onenner=1.0d0/(nan*nbn)
+      if (abs(nan*nbn).gt.eps) then
+         onenner=1.0d0/(nan*nbn)
+      else
+         onenner=0.0d0
+      endif
    else
       onenner=1.d0/nenner
    endif
@@ -630,21 +634,27 @@ Subroutine dphidrPBC(mode,nat,xyz,i,j,k,l,vTrR,vTrB,vTrC,phi,&
    call crossprod(rbpc,na,rbpca)
    call crossprod(rbpc,nb,rbpcb)
 
-   ! ... dphidri
-   do ic=1,3
-      dphidri(ic)=onenner*(cosphi*nbn/nan*rab(ic)-rbb(ic))
-
-      ! ... dphidrj
-      dphidrj(ic)=onenner*(cosphi*(nbn/nan*rapba(ic)&
-         &                                +nan/nbn*rbc(ic))&
-         &                        -(rac(ic)+rapbb(ic)))
-      ! ... dphidrk
-      dphidrk(ic)=onenner*(cosphi*(nbn/nan*raa(ic)&
-         &                             +nan/nbn*rbpcb(ic))&
-         &                        -(rba(ic)+rbpca(ic)))
-      ! ... dphidrl
-      dphidrl(ic)=onenner*(cosphi*nan/nbn*rbb(ic)-rab(ic))
-   end do
+   if (abs(onenner).gt.eps) then
+      do ic=1,3
+         ! ... dphidri
+         dphidri(ic)=onenner*(cosphi*nbn/nan*rab(ic)-rbb(ic))
+         ! ... dphidrj
+         dphidrj(ic)=onenner*(cosphi*(nbn/nan*rapba(ic)&
+            &                                +nan/nbn*rbc(ic))&
+            &                        -(rac(ic)+rapbb(ic)))
+         ! ... dphidrk
+         dphidrk(ic)=onenner*(cosphi*(nbn/nan*raa(ic)&
+            &                             +nan/nbn*rbpcb(ic))&
+            &                        -(rba(ic)+rbpca(ic)))
+         ! ... dphidrl
+         dphidrl(ic)=onenner*(cosphi*nan/nbn*rbb(ic)-rab(ic))
+      end do
+   else
+      dphidri=0.0d0
+      dphidrj=0.0d0
+      dphidrk=0.0d0
+      dphidrl=0.0d0
+   endif
 
 End subroutine
 

diff --git a/src/optimizer.f90 b/src/optimizer.f90
@@ -850,7 +850,7 @@ subroutine relax(env,iter,mol,anc,restart,maxcycle,maxdispl,ethr,gthr, &
       write(env%unit,'(5x,"change   ",e18.7,1x,"Eh")')                      echng
       write(env%unit,'(3x,"gradient norm :",f14.7,1x,"Eh/α")',advance='no') gnorm
       write(env%unit,'(3x,"predicted",e18.7)',advance='no')                 depred
-      write(env%unit,'(1x,"("f7.2"%)")')         (depred-echng)/echng*100
+      write(env%unit,'(1x,"("f7.2"%)")')         (depred-echng)/(echng+1e-34_wp)*100
    endif
 
    ! check 0 energy case !