Don't fail on hessian calculations for atoms (#587)

src/eeq_model.f90

@@ -350,17 +350,17 @@ subroutine goedecker_chrgeq(n,at,xyz,chrg,cn,dcndr,q,dqdr,energy,gradient,&
    real(wp),intent(in)    :: xyz(3,n)      ! geometry
    real(wp),intent(in)    :: chrg          ! total charge
    real(wp),intent(in)    :: cn(n)         ! erf-CN
-   real(wp),intent(in)    :: dcndr(3,n,n)  ! derivative of erf-CN
+   real(wp),intent(in),optional :: dcndr(3,n,n)  ! derivative of erf-CN
    logical, intent(in)    :: lverbose      ! toggles printout
    logical, intent(in)    :: lgrad         ! flag for gradient calculation
    logical, intent(in)    :: lcpq          ! do partial charge derivative
-!! ------------------------------------------------------------------------
+!! --------
 !  Output
 !! ------------------------------------------------------------------------
-   real(wp),intent(out)   :: q(n)          ! partial charges
-   real(wp),intent(out)   :: dqdr(3,n,n) ! derivative of partial charges
-   real(wp),intent(inout) :: energy        ! electrostatic energy
-   real(wp),intent(inout) :: gradient(3,n) ! molecular gradient of IES
+   real(wp),intent(out), optional   :: q(n)          ! partial charges
+   real(wp),intent(out), optional   :: dqdr(3,n,n) ! derivative of partial charges
+   real(wp),intent(inout), optional :: energy        ! electrostatic energy
+   real(wp),intent(inout), optional :: gradient(3,n) ! molecular gradient of IES
 !
 !! ------------------------------------------------------------------------
 !  charge model

src/hessian.F90

@@ -402,30 +402,32 @@ subroutine numhess( &
    end if
 
    ! sort such that rot/trans are modes 1:6, H/isqm are scratch
-   kend=6
-   if(res%linear)then
-      kend=5
-      do i=1,kend
-         izero(i)=i
-      enddo
-      res%freq(1:5)=0
-   endif
-   do k=1,kend
-      h(1:n3,k)=res%hess(1:n3,izero(k))
-      isqm(  k)=res%freq(izero(k))
-   enddo
-   j=kend
-   do k=1,n3
-      if(abs(res%freq(k)).gt.0.01_wp)then
-         j=j+1
-         if(j.gt.n3) then
-            call env%error('internal error while sorting hessian', source)
-            return
-         end if
-         h(1:n3,j)=res%hess(1:n3,k)
-         isqm(  j)=res%freq(   k)
+   if (mol%n > 1) then
+      kend=6
+      if(res%linear)then
+         kend=5
+         do i=1,kend
+            izero(i)=i
+         enddo
+         res%freq(1:5)=0
       endif
-   enddo
+      do k=1,kend
+         h(1:n3,k)=res%hess(1:n3,izero(k))
+         isqm(  k)=res%freq(izero(k))
+      enddo
+      j=kend
+      do k=1,n3
+         if(abs(res%freq(k)).gt.0.01_wp)then
+            j=j+1
+            if(j.gt.n3) then
+               call env%error('internal error while sorting hessian', source)
+               return
+            end if
+            h(1:n3,j)=res%hess(1:n3,k)
+            isqm(  j)=res%freq(   k)
+         endif
+      enddo
+   end if
    res%hess = h
    res%freq = isqm
    call PREIGF(env%unit,res%freq,res%n3true)

src/printout.f90

@@ -56,7 +56,7 @@ subroutine writecosmofile(np,pa,espe,fname,nat,at,xyz,atom_weight)
 
 end subroutine writecosmofile
 
-subroutine setup_summary(iunit,n,fname,xcontrol,wfx,xrc,exist)
+subroutine setup_summary(iunit,n,fname,xcontrol,wfx,xrc)
    use xtb_mctc_accuracy, only : wp
    use xtb_mctc_global, only : persistentEnv
    use xtb_mctc_systools
@@ -66,12 +66,11 @@ subroutine setup_summary(iunit,n,fname,xcontrol,wfx,xrc,exist)
    implicit none
    type(TWavefunction),intent(in) :: wfx
    integer, intent(in) :: iunit
-   character(len=*),intent(in) :: xrc
+   character(len=*),intent(in),optional :: xrc
    character(len=*),intent(in) :: xcontrol
    character(len=*),intent(in) :: fname
    integer :: i,l,err
    integer,intent(in) :: n
-   logical,intent(in) :: exist
    real(wp) :: dum5
    character(len=:),allocatable :: cdum
    write(iunit,'(a)')
@@ -94,7 +93,7 @@ subroutine setup_summary(iunit,n,fname,xcontrol,wfx,xrc,exist)
       write(iunit,'(10x,a,":",1x,a)') 'xtbhome directory          ',xenv%home
       write(iunit,'(10x,a,":",1x,a)') 'path for xtb               ',xenv%path
       write(iunit,'(10x,a,":",1x,a)') 'xcontrol input file        ',xcontrol
-      if (exist) &
+      if (present(xrc)) &
          write(iunit,'(10x,a,":",1x,a)') 'global configurations file ',xrc
    endif
    ! ----------------------------------------------------------------------

src/prog/main.F90

@@ -387,7 +387,7 @@ subroutine xtbMain(env, argParser)
    mol%uhf = set%nalphabeta
    call initrand
 
-   call setup_summary(env%unit,mol%n,fname,xcontrol,chk%wfn,xrc,exist)
+   call setup_summary(env%unit,mol%n,fname,xcontrol,chk%wfn,xrc)
 
    if(set%fit) acc=0.2 ! higher SCF accuracy during fit
 

src/scf_module.F90

@@ -289,7 +289,7 @@ subroutine scf(env, mol, wfn, basis, pcem, xtbData, solvation, &
    call qsh2qat(basis%ash,wfn%qsh,wfn%q)
 
 !  # atom arrays
-   allocate(qq(mol%n),qlmom(3,mol%n),cm5(mol%n),sqrab(mol%n*(mol%n+1)/2))
+   allocate(qq(mol%n),qlmom(3,mol%n),cm5(mol%n),cm5a(mol%n),sqrab(mol%n*(mol%n+1)/2))
    allocate(dcndr(3,mol%n,mol%n),cn(mol%n),dcndL(3,3,mol%n))
 
 !  initialize the GBSA module (GBSA works with CM5 charges)
@@ -299,7 +299,7 @@ subroutine scf(env, mol, wfn, basis, pcem, xtbData, solvation, &
          return
       end if
       call solvation%update(env, mol%at, mol%xyz)
-      allocate(cm5a(mol%n),dcm5a(3,mol%n,mol%n))
+      allocate(dcm5a(3,mol%n,mol%n))
       gborn=0._wp
       gsasa=0._wp
       ghb=0._wp