Merge pull request #213 from openmopac/api-error-handling

API error handler

src/SCF/fock2.F90

@@ -82,14 +82,12 @@ subroutine fock2(f, ptot, p, w, wj, wk, numat, nfirst, nlast, mode)
         if (allocated(ptot2))  deallocate(ptot2)
         if (allocated(ifact))  deallocate(ifact)
         if (allocated(i1fact)) deallocate(i1fact)
-        allocate(ptot2(numat,81), ifact(3 + norbs), i1fact(3 + norbs))
+        allocate(ptot2(max(2,numat),81), ifact(max(18,norbs)), i1fact(max(18,norbs)))
         icalcn = numcal
 !
 !   SET UP ARRAY OF LOWER HALF TRIANGLE INDICES (PASCAL'S TRIANGLE)
 !
-        ifact = 0  ! the bookkeeping used in dhc can read past ifact(norbs), which
-        i1fact = 0 ! explains the extended size of these arrays and their need for zero-padding
-        do i = 1, norbs
+        do i = 1, max(18,norbs)
           ifact(i) = (i*(i - 1))/2
           i1fact(i) = ifact(i) + i
         end do

src/interface/mopac_api.F90

@@ -42,8 +42,6 @@ module mopac_api
     integer, dimension (:), allocatable :: atom
     ! (x,y,z) coordinates of each atom (Angstroms) [3*natom]
     double precision, dimension (:), allocatable :: coord
-    ! flag to determine if each atom is allowed to move [natom]
-    logical, dimension (:), allocatable :: move_atom
     ! number of lattice vectors / translation vectors / periodic dimensions
     integer :: nlattice = 0
     ! number of lattice vectors that are allowed to move (first nlattice_move vectors in array)
@@ -77,7 +75,7 @@ module mopac_api
     ! (x,y,z) displacement vectors of normal modes [3*natom_move,3*natom_move]
     double precision, dimension (:,:), allocatable :: disp
     ! bond-order matrix in compressed sparse column (CSC) matrix format
-    ! with insignificant bond orders (<0.001) truncated
+    ! with insignificant bond orders (<0.01) truncated
     ! diagonal matrix entries are atomic valencies
     ! > first index of each atom in CSC bond-order matrix [natom+1]
     integer, dimension (:), allocatable :: bond_index
@@ -89,11 +87,12 @@ module mopac_api
     double precision, dimension (:), allocatable :: lattice_update
     ! (x,y,z) heat gradients for each moveable lattice vector (kcal/mol/Angstrom) [3*nlattice_move]
     double precision, dimension (:), allocatable :: lattice_deriv
-    ! stress tensor (Gigapascals) in Voigt form (xx, yy, zz, yz, xz, xy), if nlattice_move == 3
+    ! stress tensor (Gigapascals) in Voigt form (xx, yy, zz, yz, xz, xy), if available
     double precision, dimension (6) :: stress
-    ! status of MOPAC job
-    integer :: status
-    ! TO DO: compile list of status values & their meaning
+    ! number of MOPAC error messages (negative value indicates that allocation of error_msg failed)
+    integer :: nerror
+    ! text of MOPAC error messages [nerror,120]
+    character*120, dimension(:), allocatable :: error_msg
   end type
 
   ! data that describes the electronic state using standard molecular orbitals
@@ -141,11 +140,11 @@ module mopac_api
     ! > size of array cocc
     integer :: cocc_dim
     ! > atomic orbital coefficients of the occupied LMOs [cocc_dim]
-    integer, dimension (:), allocatable :: cocc
+    double precision, dimension (:), allocatable :: cocc
     ! > size of array cvir
     integer :: cvir_dim
     ! > atomic orbital coefficients of the virtual LMOs [cvir_dim]
-    integer, dimension (:), allocatable :: cvir
+    double precision, dimension (:), allocatable :: cvir
   end type
 
   interface

src/interface/mopac_api_finalize.F90

@@ -25,12 +25,15 @@
     loc, & ! indices of atoms and coordinates marked for optimization
     bondab ! bond order matrix in packed triangular format
   use molkst_C, only : escf, & ! heat of formation
+    moperr, & ! error status
     numat, & ! number of real atoms
     nvar, & ! number of coordinates to be optimized
     id, & ! number of lattice vectors
     keywrd, & ! keyword string to adjust MOPAC behavior
     voigt, & ! Voigt stress tensor
     mozyme, & ! logical flag for MOZYME calculations
+    dummy, & ! dummy integer
+    errtxt, & ! most recent error message
     use_disk ! logical flag to enable disk access
   use MOZYME_C, only : iorbs ! number of atomic orbitals for each atom
   use parameters_C, only : tore ! number of valence electrons per element
@@ -44,13 +47,11 @@
   ! save properties and clean up after a MOPAC/MOZYME calculation
   module subroutine mopac_finalize(properties)
     type(mopac_properties), intent(out) :: properties
-    integer, external :: ijbo
-    double precision, external :: dipole, dipole_for_MOZYME
-    integer :: i, j, k, kk, kl, ku, io, jo, natom_move, nlattice_move
-    double precision :: valenc, sum, dumy(3)
+    integer :: status, i
 
     ! close dummy output file to free up /dev/null
     close(iw)
+
     ! deallocate any prior arrays
     if (allocated(properties%coord_update)) deallocate(properties%coord_update)
     if (allocated(properties%coord_deriv)) deallocate(properties%coord_deriv)
@@ -62,19 +63,63 @@ module subroutine mopac_finalize(properties)
     if (allocated(properties%bond_order)) deallocate(properties%bond_order)
     if (allocated(properties%lattice_update)) deallocate(properties%lattice_update)
     if (allocated(properties%lattice_deriv)) deallocate(properties%lattice_deriv)
+
+    ! record properties
+    if (.not. moperr) call mopac_record(properties)
+
+    ! collect error messages
+    if (moperr) then
+      call summary("",0)
+      properties%nerror = dummy
+      allocate(properties%error_msg(properties%nerror), stat=status)
+      if (status /= 0) then
+        properties%nerror = -properties%nerror
+      else
+        do i=1, properties%nerror
+          call summary("",-i)
+          properties%error_msg(i) = trim(errtxt)
+        end do
+      end if
+      call summary("",-abs(properties%nerror)-1)
+    else
+      properties%nerror = 0
+    end if
+
+    ! deallocate memory
+    call setup_mopac_arrays(0,0)
+    if (mozyme) call delete_MOZYME_arrays()
+    ! turn use_disk back on
+    use_disk = .true.
+  end subroutine mopac_finalize
+
+  subroutine mopac_record(properties)
+    type(mopac_properties), intent(out) :: properties
+    integer, external :: ijbo
+    double precision, external :: dipole, dipole_for_MOZYME
+    integer :: status, i, j, k, kk, kl, ku, io, jo, natom_move, nlattice_move
+    double precision :: valenc, sum, dumy(3)
+
     ! trigger charge & dipole calculation
     call chrge (p, q)
     q(:numat) = tore(nat(:numat)) - q(:numat)
-    if (mozyme) then
-      sum = dipole_for_MOZYME(dumy, 2)
-      properties%dipole = dumy
+    if (id == 0) then
+      if (mozyme) then
+        sum = dipole_for_MOZYME(dumy, 2)
+        properties%dipole = dumy
+      else
+        sum = dipole(p, xparam, dumy, 1)
+        properties%dipole = dip(:3,3)
+      end if
     else
-      sum = dipole(p, xparam, dumy, 1)
-      properties%dipole = dip(:3,3)
+      properties%dipole = 0.d0
     end if
     ! save basic properties
     properties%heat = escf
-    allocate(properties%charge(numat))
+    allocate(properties%charge(numat), stat=status)
+    if (status /= 0) then
+      call mopend("Failed to allocate memory in MOPAC_FINALIZE")
+      return
+    end if
     properties%charge = q(:numat)
     properties%stress = voigt
     natom_move = nvar/3
@@ -86,28 +131,56 @@ module subroutine mopac_finalize(properties)
       end if
     end do
     ! save properties of moveable coordinates
-    allocate(properties%coord_update(3*natom_move))
+    allocate(properties%coord_update(3*natom_move), stat=status)
+    if (status /= 0) then
+      call mopend("Failed to allocate memory in MOPAC_FINALIZE")
+      return
+    end if
     properties%coord_update = xparam(:3*natom_move)
-    allocate(properties%coord_deriv(3*natom_move))
+    allocate(properties%coord_deriv(3*natom_move), stat=status)
+    if (status /= 0) then
+      call mopend("Failed to allocate memory in MOPAC_FINALIZE")
+      return
+    end if
     properties%coord_deriv = grad(:3*natom_move)
     if (nlattice_move > 0) then
-      allocate(properties%lattice_update(3*nlattice_move))
+      allocate(properties%lattice_update(3*nlattice_move), stat=status)
+      if (status /= 0) then
+        call mopend("Failed to allocate memory in MOPAC_FINALIZE")
+        return
+      end if
       properties%lattice_update = xparam(3*natom_move+1:)
-      allocate(properties%lattice_deriv(3*nlattice_move))
+      allocate(properties%lattice_deriv(3*nlattice_move), stat=status)
+      if (status /= 0) then
+        call mopend("Failed to allocate memory in MOPAC_FINALIZE")
+        return
+      end if
       properties%lattice_deriv = grad(3*natom_move+1:)
     end if
     ! save vibrational properties if available
-    if (index(keywrd, " FORCE") /= 0) then
+    if (index(keywrd, " FORCETS") /= 0) then
       properties%calc_vibe = .true.
-      allocate(properties%freq(nvar))
-      allocate(properties%disp(nvar,nvar))
+      allocate(properties%freq(nvar), stat=status)
+      if (status /= 0) then
+        call mopend("Failed to allocate memory in MOPAC_FINALIZE")
+        return
+      end if
+      allocate(properties%disp(nvar,nvar), stat=status)
+      if (status /= 0) then
+        call mopend("Failed to allocate memory in MOPAC_FINALIZE")
+        return
+      end if
       properties%freq = freq
       properties%disp = reshape(cnorml,[nvar, nvar])
     else
       properties%calc_vibe = .false.
     end if
     ! prune bond order matrix
-    allocate(properties%bond_index(numat+1))
+    allocate(properties%bond_index(numat+1), stat=status)
+    if (status /= 0) then
+      call mopend("Failed to allocate memory in MOPAC_FINALIZE")
+      return
+    end if
     if (mozyme) then
       ! 1st pass to populate bond_index
       properties%bond_index(1) = 1
@@ -129,7 +202,7 @@ module subroutine mopac_finalize(properties)
             valenc = valenc + 2.d0 * p(kk)
           end do
         end if
-        if (valenc > 0.001d0) then
+        if (valenc > 0.01d0) then
           properties%bond_index(i+1) = properties%bond_index(i+1) + 1
         end if
         do j = 1, numat
@@ -141,15 +214,23 @@ module subroutine mopac_finalize(properties)
             do k = kl, ku
               sum = sum + p(k) ** 2
             end do
-            if (sum > 0.001d0) then
+            if (sum > 0.01d0) then
               properties%bond_index(i+1) = properties%bond_index(i+1) + 1
             end if
           end if
         end do
       end do
       ! 2nd pass to populate bond_atom and bond_order
-      allocate(properties%bond_atom(properties%bond_index(numat+1)))
-      allocate(properties%bond_order(properties%bond_index(numat+1)))
+      allocate(properties%bond_atom(properties%bond_index(numat+1)), stat=status)
+      if (status /= 0) then
+        call mopend("Failed to allocate memory in MOPAC_FINALIZE")
+        return
+      end if
+      allocate(properties%bond_order(properties%bond_index(numat+1)), stat=status)
+      if (status /= 0) then
+        call mopend("Failed to allocate memory in MOPAC_FINALIZE")
+        return
+      end if
       do i = 1, numat
         io = iorbs(i)
         valenc = 0.d0
@@ -177,12 +258,12 @@ module subroutine mopac_finalize(properties)
             do k = kl, ku
               sum = sum + p(k) ** 2
             end do
-            if (sum > 0.001d0) then
+            if (sum > 0.01d0) then
               properties%bond_atom(kk) = j
               properties%bond_order(kk) = sum
               kk = kk + 1
             end if
-          else if (valenc > 0.001d0) then
+          else if (valenc > 0.01d0) then
             properties%bond_atom(kk) = j
             properties%bond_order(kk) = valenc
             kk = kk + 1
@@ -198,35 +279,43 @@ module subroutine mopac_finalize(properties)
         ku = i*(i-1)/2 + 1
         kl = (i+1)*(i+2)/2 - 1
         do j = 1, i
-          if (bondab(ku) > 0.001d0) then
+          if (bondab(ku) > 0.01d0) then
             properties%bond_index(i+1) = properties%bond_index(i+1) + 1
           end if
           ku = ku + 1
         end do
         do j = i+1, numat
-          if (bondab(kl) > 0.001d0) then
+          if (bondab(kl) > 0.01d0) then
             properties%bond_index(i+1) = properties%bond_index(i+1) + 1
           end if
           kl = kl + j
         end do
       end do
       ! 2nd pass to populate bond_atom and bond_order
-      allocate(properties%bond_atom(properties%bond_index(numat+1)))
-      allocate(properties%bond_order(properties%bond_index(numat+1)))
+      allocate(properties%bond_atom(properties%bond_index(numat+1)), stat=status)
+      if (status /= 0) then
+        call mopend("Failed to allocate memory in MOPAC_FINALIZE")
+        return
+      end if
+      allocate(properties%bond_order(properties%bond_index(numat+1)), stat=status)
+      if (status /= 0) then
+        call mopend("Failed to allocate memory in MOPAC_FINALIZE")
+        return
+      end if
       do i = 1, numat
         ku = i*(i-1)/2 + 1
         kl = (i+1)*(i+2)/2 - 1
         kk = properties%bond_index(i)
         do j = 1, i
-          if (bondab(ku) > 0.001d0) then
+          if (bondab(ku) > 0.01d0) then
             properties%bond_atom(kk) = j
             properties%bond_order(kk) = bondab(ku)
             kk = kk + 1
           end if
           ku = ku + 1
         end do
         do j = i+1, numat
-          if (bondab(kl) > 0.001d0) then
+          if (bondab(kl) > 0.01d0) then
             properties%bond_atom(kk) = j
             properties%bond_order(kk) = bondab(kl)
             kk = kk + 1
@@ -235,13 +324,6 @@ module subroutine mopac_finalize(properties)
         end do
       end do
     end if
-    ! deallocate memory
-    call setup_mopac_arrays(0,0)
-    if (mozyme) call delete_MOZYME_arrays()
-    ! turn use_disk back on
-    use_disk = .true.
-    ! mark the job as successful
-    properties%status = 0
-  end subroutine mopac_finalize
+  end subroutine mopac_record
 
 end submodule mopac_api_finalize