diff --git a/src/SCF/fock2.F90 b/src/SCF/fock2.F90
index 6c8f5aab..ea414908 100644
--- a/src/SCF/fock2.F90
+++ b/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
diff --git a/src/interface/mopac_api.F90 b/src/interface/mopac_api.F90
index 2cd25db1..291a59ad 100644
--- a/src/interface/mopac_api.F90
+++ b/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
diff --git a/src/interface/mopac_api_finalize.F90 b/src/interface/mopac_api_finalize.F90
index c5685f4a..6e88aebf 100644
--- a/src/interface/mopac_api_finalize.F90
+++ b/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,27 +279,35 @@ 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
@@ -226,7 +315,7 @@ module subroutine mopac_finalize(properties)
           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
diff --git a/src/interface/mopac_api_initialize.F90 b/src/interface/mopac_api_initialize.F90
index b4fc97a7..34e9f5bf 100644
--- a/src/interface/mopac_api_initialize.F90
+++ b/src/interface/mopac_api_initialize.F90
@@ -29,7 +29,7 @@
     na, nb, nc, & ! internal coordinate connectivity information
     breaks, chains, cell_ijk, nw, nfirst, nlast
   use cosmo_C, only : iseps, & ! flag for use of COSMO model
-    useps, lpka, solv_energy, area, fepsi, ediel
+    noeps, useps, lpka, solv_energy, area, fepsi, ediel, nspa
   use funcon_C, only : fpc, fpc_9 ! fundamental constants used in the MOPAC calculation
   use maps_C, only : latom, lparam, lpara1, latom1, lpara2, latom2, rxn_coord
   use meci_C, only : nmos, lab
@@ -72,14 +72,33 @@
   ! investment of development time. Similarly, the keyword line is set to represent
   ! the calculation that the API is trying to perform, but some of the keyword-dependent
   ! logic has been stripped out of this restricted initialization.
-  module subroutine mopac_initialize(system, status)
+  module subroutine mopac_initialize(system)
     type(mopac_system), intent(in) :: system
-    integer, intent(out) :: status
     double precision, external :: seconds, C_triple_bond_C
     character(100) :: num2str
-    integer :: i, j, nelectron
+    integer :: i, j, nelectron, status
     double precision :: eat
 
+    ! validity checks of data in system
+    if (system%natom <= 0) call mopend("Invalid number of atoms")
+    if (system%natom_move < 0 .or. system%natom_move > system%natom) &
+      call mopend("Invalid number of moveable atoms")
+    if (system%epsilon /= 1.d0 .and. system%nlattice > 0) &
+      call mopend("COSMO solvent not available for periodic systems")
+    if (size(system%atom) < system%natom) call mopend("List of atomic numbers is too small")
+    if (size(system%coord) < 3*system%natom) call mopend("List of atomic coordinates is too small")
+    if (system%nlattice < 0 .or. system%nlattice > 3) call mopend("Invalid number of lattice vectors")
+    if (system%nlattice_move < 0 .or. system%nlattice_move > system%nlattice) &
+      call mopend("Invalid number of moveable lattice vectors")
+    if (system%nlattice > 0) then
+      if(size(system%lattice) < 3*system%nlattice) call mopend("List of lattice vectors is too small")
+    end if
+    if (system%tolerance <= 0.d0) call mopend("Relative numerical tolerance must be a positive number")
+    if (system%max_time <= 0) call mopend("Time limit must be a positive number")
+    if (system%nlattice_move > 0 .and. index(keywrd," FORCETS") /= 0) &
+      call mopend("Lattice vectors cannot move during vibrational calculations")
+    if (moperr) return
+
     use_disk = .false.
     ! load ios, iop, and iod data into tore
     tore = ios + iop + iod
@@ -94,7 +113,7 @@ module subroutine mopac_initialize(system, status)
     ! initialize CODATA fundamental constants
     fpc(:) = fpcref(1,:)
     ! assemble virtual keyword line
-    keywrd = trim(keywrd) // " LET NOSYM"
+    keywrd = trim(keywrd) // " LET NOSYM GEO-OK"
     write(num2str,'(i6)') system%charge
     keywrd = trim(keywrd) // " CHARGE=" // adjustl(num2str)
     if (system%pressure /= 0.d0) then
@@ -112,9 +131,10 @@ module subroutine mopac_initialize(system, status)
     if (MOPAC_OS == "Windows") output_fn = 'NUL'
 #endif
     close(iw)
-    open(unit=iw, file=output_fn, status='UNKNOWN', position='asis', iostat = i)
-    if (i /= 0) then
-        ! TO DO: error handling ???
+    open(unit=iw, file=output_fn, status='UNKNOWN', position='asis', iostat=status)
+    if (status /= 0) then
+      call mopend("Failed to open NULL file in MOPAC_INITIALIZE")
+      return
     end if
     ! initialize job timers
     time0 = seconds(1)
@@ -128,7 +148,13 @@ module subroutine mopac_initialize(system, status)
     mozyme = (index(keywrd," MOZ") /= 0)
     ! initialize common workspaces (Common_arrays_C)
     call setup_mopac_arrays(natoms, 1)
-    if (.not. allocated(lopt)) allocate(lopt(3,maxatoms))
+    if (.not. allocated(lopt)) then
+      allocate(lopt(3,maxatoms), stat=status)
+      if (status /= 0) then
+        call mopend("Failed to allocate memory in MOPAC_INITIALIZE")
+        return
+      end if
+    end if
     ! set semiempirical model
     methods = .false.
     select case (system%model)
@@ -145,7 +171,7 @@ module subroutine mopac_initialize(system, status)
       case (5) ! RM1
         i = 4
       case default
-! TO DO: handle error
+        call mopend("Unknown semiempirical model requested")
     end select
     methods(i) = .true.
     keywrd = trim(keywrd) // methods_keys(i)
@@ -154,22 +180,28 @@ module subroutine mopac_initialize(system, status)
     ! check for use of COSMO model
     if (system%epsilon /= 1.d0) then
       iseps = .true.
+      useps = .true.
+      noeps = .false.
       write(num2str,'(f6.2)') system%epsilon
       keywrd = trim(keywrd) // " EPS=" // adjustl(num2str)
+      fepsi = (system%epsilon-1.d0) / (system%epsilon+0.5d0)
     else
       iseps = .false.
+      useps = .false.
+      noeps = .true.
+      fepsi = 0.d0
     end if
     ! insert geometry information into MOPAC data structures
     id = system%nlattice
     nat(:system%natom) = system%atom(:)
     labels(:system%natom) = system%atom(:)
-    atmass(:natoms) = ams(labels(:))
     geo(:,:system%natom) = reshape(system%coord,[3, system%natom])
     if (id > 0) then
       nat(system%natom+1:system%natom+id) = 107
       labels(system%natom+1:system%natom+id) = 107
       geo(:,system%natom+1:system%natom+id) = reshape(system%lattice,[3, id])
     end if
+    atmass(:natoms) = ams(labels(:))
     ! set optimization flags & xparam
     nvar = 3*system%natom_move + 3*system%nlattice_move
     lopt(:,:system%natom_move) = 1
@@ -229,12 +261,11 @@ module subroutine mopac_initialize(system, status)
     chains = " "
     breaks(1) = -300
     ! initialize variables that need default values (cosmo_C)
-    useps = .false.
     lpka = .false.
     solv_energy = 0.d0
     area = 0.d0
-    fepsi = 0.d0
     ediel = 0.d0
+    nspa = 42
     ! initialize variables that need default values (maps_C)
     latom = 0
     lparam = 0
@@ -274,7 +305,11 @@ module subroutine mopac_initialize(system, status)
     ! initialize system-specific MOPAC workspaces (Common_arrays_C)
     call setup_mopac_arrays(1,2)
     if (allocated(nw)) deallocate(nw)
-    allocate(nw(numat))
+    allocate(nw(numat), stat=status)
+    if (status /= 0) then
+      call mopend("Failed to allocate memory in MOPAC_INITIALIZE")
+      return
+    end if
     j = 1
     do i = 1, numat
       nw(i) = j
@@ -294,8 +329,6 @@ module subroutine mopac_initialize(system, status)
     atheat = atheat + C_triple_bond_C()
     ! setup MOZYME calculations
     if (mozyme) call set_up_MOZYME_arrays()
-    ! successful setup
-    status = 0
   end subroutine mopac_initialize
 
 end submodule mopac_api_initialize
diff --git a/src/interface/mopac_api_operations.F90 b/src/interface/mopac_api_operations.F90
index c95bf22d..5ee48747 100644
--- a/src/interface/mopac_api_operations.F90
+++ b/src/interface/mopac_api_operations.F90
@@ -22,9 +22,8 @@
   interface
 
     ! initialize MOPAC modules before calculations
-    module subroutine mopac_initialize(system, status)
+    module subroutine mopac_initialize(system)
       type(mopac_system), intent(in) :: system
-      integer, intent(out) :: status
     end subroutine mopac_initialize
 
     ! extract data from MOPAC modules after calculations
@@ -33,27 +32,23 @@ module subroutine mopac_finalize(properties)
     end subroutine mopac_finalize
 
     ! save MOPAC density matrices
-    module subroutine mopac_save(state, status)
+    module subroutine mopac_save(state)
       type(mopac_state), intent(out) :: state
-      integer, intent(out) :: status
     end subroutine mopac_save
 
     ! load MOPAC density matrices, or construct initial guesses
-    module subroutine mopac_load(state, status)
+    module subroutine mopac_load(state)
       type(mopac_state), intent(in) :: state
-      integer, intent(out) :: status
     end subroutine mopac_load
 
     ! save MOZYME density matrix
-    module subroutine mozyme_save(state, status)
+    module subroutine mozyme_save(state)
       type(mozyme_state), intent(out) :: state
-      integer, intent(out) :: status
     end subroutine mozyme_save
 
     ! load MOZYME density matrix, or construct initial guess
-    module subroutine mozyme_load(state, status)
+    module subroutine mozyme_load(state)
       type(mozyme_state), intent(in) :: state
-      integer, intent(out) :: status
     end subroutine mozyme_load
 
   end interface
@@ -62,197 +57,118 @@ end subroutine mozyme_load
 
   ! MOPAC electronic ground state calculation
   module subroutine mopac_scf(system, state, properties)
-    !dec$ attributes dllexport :: mopac_scf
-      type(mopac_system), intent(in) :: system
-      type(mopac_state), intent(inout) :: state
-      type(mopac_properties), intent(out) :: properties
-      integer :: status
-
-      keywrd = " 1SCF PULAY GRADIENTS BONDS"
-      call mopac_initialize(system, status)
-      if (status /= 0) then
-        properties%status = status
-        return
-      end if
-      call mopac_load(state, status)
-      if (status /= 0) then
-        properties%status = status
-        return
-      end if
-      ! call computational routine for SCF calculations
-      call compfg (xparam, .true., escf, .true., grad, .true.)
-      ! TO DO: compfg error handling
-      call mopac_save(state, status)
-      if (status /= 0) then
-        properties%status = status
-        return
-      end if
-      call mopac_finalize(properties)
-    end subroutine mopac_scf
-  
-    ! MOPAC geometry relaxation
-    module subroutine mopac_relax(system, state, properties)
-    !dec$ attributes dllexport :: mopac_relax
-      type(mopac_system), intent(in) :: system
-      type(mopac_state), intent(inout) :: state
-      type(mopac_properties), intent(out) :: properties
-      integer :: status
-
-      keywrd = " LBFGS PULAY BONDS"
-      call mopac_initialize(system, status)
-      if (status /= 0) then
-        properties%status = status
-        return
-      end if
-      call mopac_load(state, status)
-      if (status /= 0) then
-        properties%status = status
-        return
-      end if
-      ! call computational routine for LBFGS geometry optimization
-      call lbfgs (xparam, escf)
-      ! TO DO: lbfgs error handling
-      call mopac_save(state, status)
-      if (status /= 0) then
-        properties%status = status
-        return
-      end if
-      call mopac_finalize(properties)
-    end subroutine mopac_relax
-  
-    ! MOPAC vibrational calculation
-    module subroutine mopac_vibe(system, state, properties)
-    !dec$ attributes dllexport :: mopac_vibe
-      type(mopac_system), intent(in) :: system
-      type(mopac_state), intent(inout) :: state
-      type(mopac_properties), intent(out) :: properties
-      integer :: status, i
-
-      keywrd = " FORCETS PULAY BONDS"
-      call mopac_initialize(system, status)
-      if (status /= 0) then
-        properties%status = status
-        return
-      end if
-      call mopac_load(state, status)
-      if (status /= 0) then
-        properties%status = status
-        return
-      end if
-      ! call computational routine to evaluate vibrational matrix
-      call force ()
-      ! TO DO: force error handling
-      ! recompute nvar because it is zero'd after vibrational calculations
+  !dec$ attributes dllexport :: mopac_scf
+    type(mopac_system), intent(in) :: system
+    type(mopac_state), intent(inout) :: state
+    type(mopac_properties), intent(out) :: properties
+
+    keywrd = " 1SCF PULAY GRADIENTS BONDS"
+    call mopac_initialize(system)
+    if (.not. moperr) call mopac_load(state)
+    ! call computational routine for SCF calculations
+    if (.not. moperr) call compfg (xparam, .true., escf, .true., grad, .true.)
+    if (.not. moperr) call mopac_save(state)
+    call mopac_finalize(properties)
+  end subroutine mopac_scf
+
+  ! MOPAC geometry relaxation
+  module subroutine mopac_relax(system, state, properties)
+  !dec$ attributes dllexport :: mopac_relax
+    type(mopac_system), intent(in) :: system
+    type(mopac_state), intent(inout) :: state
+    type(mopac_properties), intent(out) :: properties
+
+    keywrd = " LBFGS PULAY BONDS"
+    call mopac_initialize(system)
+    if (.not. moperr) call mopac_load(state)
+    ! call computational routine for LBFGS geometry optimization
+    if (.not. moperr) call lbfgs (xparam, escf)
+    if (.not. moperr) call mopac_save(state)
+    call mopac_finalize(properties)
+  end subroutine mopac_relax
+
+  ! MOPAC vibrational calculation
+  module subroutine mopac_vibe(system, state, properties)
+  !dec$ attributes dllexport :: mopac_vibe
+    type(mopac_system), intent(in) :: system
+    type(mopac_state), intent(inout) :: state
+    type(mopac_properties), intent(out) :: properties
+    integer :: i
+
+    keywrd = " FORCETS PULAY BONDS"
+    call mopac_initialize(system)
+    if (.not. moperr) call mopac_load(state)
+    ! call computational routine to evaluate vibrational matrix
+    if (.not. moperr) call force()
+    ! recompute nvar because it is zero'd after vibrational calculations
+    if (.not. moperr) then
       nvar = 0
       do i=1, system%natom+system%nlattice
         if (lopt(1,i) == 1) nvar = nvar + 3
       end do
-      ! call computational routine for standard properties calculations
-      call compfg (xparam, .true., escf, .true., grad, .true.)
-      ! TO DO: compfg error handling
-      call mopac_save(state, status)
-      if (status /= 0) then
-        properties%status = status
-        return
-      end if
-      call mopac_finalize(properties)
-    end subroutine mopac_vibe
-  
-    ! MOZYME electronic ground state calculation
-    module subroutine mozyme_scf(system, state, properties)
-    !dec$ attributes dllexport :: mozyme_scf
-      type(mopac_system), intent(in) :: system
-      type(mozyme_state), intent(inout) :: state
-      type(mopac_properties), intent(out) :: properties
-      integer :: status
-
-      keywrd = " MOZYME 1SCF GRADIENTS ALLBONDS"
-      call mopac_initialize(system, status)
-      if (status /= 0) then
-        properties%status = status
-        return
-      end if
-      call mozyme_load(state, status)
-      if (status /= 0) then
-        properties%status = status
-        return
-      end if
-      ! call computational routine for SCF calculations
-      call compfg (xparam, .true., escf, .true., grad, .true.)
-      ! TO DO: compfg error handling
-      call mozyme_save(state, status)
-      if (status /= 0) then
-        properties%status = status
-        return
-      end if
-      call mopac_finalize(properties)
-    end subroutine mozyme_scf
-  
-    ! MOZYME geometry relaxation
-    module subroutine mozyme_relax(system, state, properties)
-    !dec$ attributes dllexport :: mozyme_relax
-      type(mopac_system), intent(in) :: system
-      type(mozyme_state), intent(inout) :: state
-      type(mopac_properties), intent(out) :: properties
-      integer :: status
-
-      keywrd = " MOZYME LBFGS ALLBONDS"
-      call mopac_initialize(system, status)
-      if (status /= 0) then
-        properties%status = status
-        return
-      end if
-      call mozyme_load(state, status)
-      if (status /= 0) then
-        properties%status = status
-        return
-      end if
-      ! call computational routine for LBFGS geometry optimization
-      call lbfgs (xparam, escf)
-      ! TO DO: lbfgs error handling
-      call mozyme_save(state, status)
-      if (status /= 0) then
-        properties%status = status
-        return
-      end if
-      call mopac_finalize(properties)
-    end subroutine mozyme_relax
-
-    ! MOZYME vibrational calculation
-    module subroutine mozyme_vibe(system, state, properties)
-    !dec$ attributes dllexport :: mozyme_vibe
-        type(mopac_system), intent(in) :: system
-        type(mozyme_state), intent(inout) :: state
-        type(mopac_properties), intent(out) :: properties
-        integer :: status, i
-  
-        keywrd = " MOZYME FORCETS PULAY BONDS"
-        call mopac_initialize(system, status)
-        if (status /= 0) then
-          properties%status = status
-          return
-        end if
-        call mozyme_load(state, status)
-        if (status /= 0) then
-          properties%status = status
-          return
-        end if
-        ! call computational routine to evaluate vibrational matrix
-        call force ()
-        ! TO DO: force error handling
-        ! recompute nvar because it is zero'd after vibrational calculations
-        nvar = 0
-        do i=1, system%natom+system%nlattice
-          if (lopt(1,i) == 1) nvar = nvar + 3
-        end do
-        ! call computational routine for standard properties calculations
-        call compfg (xparam, .true., escf, .true., grad, .true.)
-        call mozyme_save(state, status)
-        if (status /= 0) then
-          properties%status = status
-          return
-        end if
-        call mopac_finalize(properties)
-      end subroutine mozyme_vibe
+    end if
+    ! call computational routine for standard properties calculations
+    if (.not. moperr) call compfg (xparam, .true., escf, .true., grad, .true.)
+    if (.not. moperr) call mopac_save(state)
+    call mopac_finalize(properties)
+  end subroutine mopac_vibe
+
+  ! MOZYME electronic ground state calculation
+  module subroutine mozyme_scf(system, state, properties)
+  !dec$ attributes dllexport :: mozyme_scf
+    type(mopac_system), intent(in) :: system
+    type(mozyme_state), intent(inout) :: state
+    type(mopac_properties), intent(out) :: properties
+
+    keywrd = " MOZYME 1SCF GRADIENTS ALLBONDS"
+    call mopac_initialize(system)
+    if (.not. moperr) call mozyme_load(state)
+    ! call computational routine for SCF calculations
+    if (.not. moperr) call compfg (xparam, .true., escf, .true., grad, .true.)
+    if (.not. moperr) call mozyme_save(state)
+    call mopac_finalize(properties)
+  end subroutine mozyme_scf
+
+  ! MOZYME geometry relaxation
+  module subroutine mozyme_relax(system, state, properties)
+  !dec$ attributes dllexport :: mozyme_relax
+    type(mopac_system), intent(in) :: system
+    type(mozyme_state), intent(inout) :: state
+    type(mopac_properties), intent(out) :: properties
+
+    keywrd = " MOZYME LBFGS ALLBONDS"
+    call mopac_initialize(system)
+    if (.not. moperr) call mozyme_load(state)
+    ! call computational routine for LBFGS geometry optimization
+    if (.not. moperr) call lbfgs (xparam, escf)
+    if (.not. moperr) call mozyme_save(state)
+    call mopac_finalize(properties)
+  end subroutine mozyme_relax
+
+  ! MOZYME vibrational calculation
+  module subroutine mozyme_vibe(system, state, properties)
+  !dec$ attributes dllexport :: mozyme_vibe
+    type(mopac_system), intent(in) :: system
+    type(mozyme_state), intent(inout) :: state
+    type(mopac_properties), intent(out) :: properties
+    integer :: i
+
+    keywrd = " MOZYME FORCETS PULAY BONDS"
+    call mopac_initialize(system)
+    if (.not. moperr) call mozyme_load(state)
+    ! call computational routine to evaluate vibrational matrix
+    if (.not. moperr) call force()
+    ! recompute nvar because it is zero'd after vibrational calculations
+    if (.not. moperr) then
+      nvar = 0
+      do i=1, system%natom+system%nlattice
+        if (lopt(1,i) == 1) nvar = nvar + 3
+      end do
+    end if
+    ! call computational routine for standard properties calculations
+    if (.not. moperr) call compfg (xparam, .true., escf, .true., grad, .true.)
+    if (.not. moperr) call mozyme_save(state)
+    call mopac_finalize(properties)
+  end subroutine mozyme_vibe
 
 end submodule mopac_api_operations
diff --git a/src/interface/mopac_api_saveload.F90 b/src/interface/mopac_api_saveload.F90
index e2af4dec..cb87131f 100644
--- a/src/interface/mopac_api_saveload.F90
+++ b/src/interface/mopac_api_saveload.F90
@@ -18,58 +18,69 @@
   use Common_arrays_C, only: pa, pb, nbonds, ibonds
   use molkst_C, only: keywrd, uhf, mpack, numat
   use MOZYME_C, only: iorbs, noccupied, ncf, nvirtual, nce, icocc_dim, &
-    icocc, icvir_dim, icvir, cocc_dim, cocc, cvir_dim, cvir
+    icocc, icvir_dim, icvir, cocc_dim, cocc, cvir_dim, cvir, nnce, nncf, ncocc, ncvir
   implicit none
 
 contains
 
   ! save MOPAC density matrices
-  module subroutine mopac_save(state, status)
+  module subroutine mopac_save(state)
     type(mopac_state), intent(out) :: state
-    integer, intent(out) :: status
+    integer :: status
 
     state%save_state = .true.
     state%mpack = mpack
 
     if (allocated(state%pa)) deallocate(state%pa)
     allocate(state%pa(mpack), stat=status)
-    if (status /= 0) return
+    if (status /= 0) then
+      call mopend("Failed to allocate memory in MOPAC_SAVE")
+      return
+    end if
     state%pa = pa
     if (uhf) then
       if (allocated(state%pb)) deallocate(state%pb)
       allocate(state%pb(mpack), stat=status)
-      if (status /= 0) return
+      if (status /= 0) then
+        call mopend("Failed to allocate memory in MOPAC_SAVE")
+        return
+      end if
       state%pb = pb
     end if
   end subroutine mopac_save
 
   ! load MOPAC density matrices, or construct initial guesses
-  module subroutine mopac_load(state, status)
+  module subroutine mopac_load(state)
     type(mopac_state), intent(in) :: state
-    integer, intent(out) :: status
+    integer :: status
 
-    status = 0
     if(state%save_state) then
       ! TO DO: compatibility tests
       keywrd = trim(keywrd) // " OLDENS"
       mpack = state%mpack
       if (allocated(pa)) deallocate(pa)
       allocate(pa(mpack), stat=status)
-      if (status /= 0) return
+      if (status /= 0) then
+        call mopend("Failed to allocate memory in MOPAC_LOAD")
+        return
+      end if  
       pa = state%pa
       if(uhf) then
         if (allocated(pb)) deallocate(pb)
         allocate(pb(mpack), stat=status)
-        if (status /= 0) return
+        if (status /= 0) then
+          call mopend("Failed to allocate memory in MOPAC_LOAD")
+          return
+        end if
         pb = state%pb
       end if
     end if
   end subroutine mopac_load
 
   ! save MOZYME density matrix
-  module subroutine mozyme_save(state, status)
+  module subroutine mozyme_save(state)
     type(mozyme_state), intent(out) :: state
-    integer, intent(out) :: status
+    integer :: status
 
     state%save_state = .true.
     state%numat = numat
@@ -89,23 +100,50 @@ module subroutine mozyme_save(state, status)
     if (allocated(state%cocc)) deallocate(state%cocc)
     if (allocated(state%cvir)) deallocate(state%cvir)
     allocate(state%nbonds(numat), stat=status)
-    if (status /= 0) return
+    if (status /= 0) then
+      call mopend("Failed to allocate memory in MOZYME_SAVE")
+      return
+    end if
     allocate(state%ibonds(9,numat), stat=status)
-    if (status /= 0) return
+    if (status /= 0) then
+      call mopend("Failed to allocate memory in MOZYME_SAVE")
+      return
+    end if
     allocate(state%iorbs(numat), stat=status)
-    if (status /= 0) return
+    if (status /= 0) then
+      call mopend("Failed to allocate memory in MOZYME_SAVE")
+      return
+    end if
     allocate(state%ncf(noccupied), stat=status)
-    if (status /= 0) return
+    if (status /= 0) then
+      call mopend("Failed to allocate memory in MOZYME_SAVE")
+      return
+    end if
     allocate(state%nce(nvirtual), stat=status)
-    if (status /= 0) return
+    if (status /= 0) then
+      call mopend("Failed to allocate memory in MOZYME_SAVE")
+      return
+    end if
     allocate(state%icocc(icocc_dim), stat=status)
-    if (status /= 0) return
+    if (status /= 0) then
+      call mopend("Failed to allocate memory in MOZYME_SAVE")
+      return
+    end if
     allocate(state%icvir(icvir_dim), stat=status)
-    if (status /= 0) return
+    if (status /= 0) then
+      call mopend("Failed to allocate memory in MOZYME_SAVE")
+      return
+    end if
     allocate(state%cocc(cocc_dim), stat=status)
-    if (status /= 0) return
+    if (status /= 0) then
+      call mopend("Failed to allocate memory in MOZYME_SAVE")
+      return
+    end if
     allocate(state%cvir(cvir_dim), stat=status)
-    if (status /= 0) return
+    if (status /= 0) then
+      call mopend("Failed to allocate memory in MOZYME_SAVE")
+      return
+    end if
     state%nbonds = nbonds
     state%ibonds = ibonds
     state%iorbs = iorbs
@@ -118,11 +156,10 @@ module subroutine mozyme_save(state, status)
   end subroutine mozyme_save
 
   ! load MOZYME density matrix, or construct initial guess
-  module subroutine mozyme_load(state, status)
+  module subroutine mozyme_load(state)
     type(mozyme_state), intent(in) :: state
-    integer, intent(out) :: status
+    integer :: status, i, j, k, l
 
-    status = 0
     if(state%save_state) then
       ! TO DO: compatibility tests
       keywrd = trim(keywrd) // " OLDENS"
@@ -142,24 +179,75 @@ module subroutine mozyme_load(state, status)
       if (allocated(icvir)) deallocate(icvir)
       if (allocated(cocc)) deallocate(cocc)
       if (allocated(cvir)) deallocate(cvir)
+      if (allocated(nncf)) deallocate(nncf)
+      if (allocated(nnce)) deallocate(nnce)
+      if (allocated(ncocc)) deallocate(ncocc)
+      if (allocated(ncvir)) deallocate(ncvir)
       allocate(nbonds(numat), stat=status)
-      if (status /= 0) return
+      if (status /= 0) then
+        call mopend("Failed to allocate memory in MOZYME_LOAD")
+        return
+      end if
       allocate(ibonds(9,numat), stat=status)
-      if (status /= 0) return
+      if (status /= 0) then
+        call mopend("Failed to allocate memory in MOZYME_LOAD")
+        return
+      end if
       allocate(iorbs(numat), stat=status)
-      if (status /= 0) return
+      if (status /= 0) then
+        call mopend("Failed to allocate memory in MOZYME_LOAD")
+        return
+      end if
       allocate(ncf(noccupied), stat=status)
-      if (status /= 0) return
+      if (status /= 0) then
+        call mopend("Failed to allocate memory in MOZYME_LOAD")
+        return
+      end if
       allocate(nce(nvirtual), stat=status)
-      if (status /= 0) return
+      if (status /= 0) then
+        call mopend("Failed to allocate memory in MOZYME_LOAD")
+        return
+      end if
       allocate(icocc(icocc_dim), stat=status)
-      if (status /= 0) return
+      if (status /= 0) then
+        call mopend("Failed to allocate memory in MOZYME_LOAD")
+        return
+      end if
       allocate(icvir(icvir_dim), stat=status)
-      if (status /= 0) return
+      if (status /= 0) then
+        call mopend("Failed to allocate memory in MOZYME_LOAD")
+        return
+      end if
       allocate(cocc(cocc_dim), stat=status)
-      if (status /= 0) return
+      if (status /= 0) then
+        call mopend("Failed to allocate memory in MOZYME_LOAD")
+        return
+      end if
       allocate(cvir(cvir_dim), stat=status)
-      if (status /= 0) return
+      if (status /= 0) then
+        call mopend("Failed to allocate memory in MOZYME_LOAD")
+        return
+      end if
+      allocate(nncf(noccupied), stat=status)
+      if (status /= 0) then
+        call mopend("Failed to allocate memory in MOZYME_LOAD")
+        return
+      end if
+      allocate(nnce(nvirtual), stat=status)
+      if (status /= 0) then
+        call mopend("Failed to allocate memory in MOZYME_LOAD")
+        return
+      end if
+      allocate(ncocc(noccupied), stat=status)
+      if (status /= 0) then
+        call mopend("Failed to allocate memory in MOZYME_LOAD")
+        return
+      end if
+      allocate(ncvir(nvirtual), stat=status)
+      if (status /= 0) then
+        call mopend("Failed to allocate memory in MOZYME_LOAD")
+        return
+      end if
       nbonds = state%nbonds
       ibonds = state%ibonds
       iorbs = state%iorbs
@@ -169,6 +257,35 @@ module subroutine mozyme_load(state, status)
       icvir = state%icvir
       cocc = state%cocc
       cvir = state%cvir
+      ! reconstruct nncf, nnce, ncocc, & ncvir
+      j = 0
+      do i = 1, noccupied
+        nncf(i) = j
+        j = j + ncf(i)
+      end do
+      j = 0
+      do i = 1, nvirtual
+        nnce(i) = j
+        j = j + nce(i)
+      end do
+      k = 0
+      do i = 1, noccupied
+        ncocc(i) = k
+        l = 0
+        do j = nncf(i) + 1, nncf(i) + ncf(i)
+          l = l + iorbs(icocc(j))
+        end do
+        k = k + l
+      end do
+      k = 0
+      do i = 1, nvirtual
+        ncvir(i) = k
+        l = 0
+        do j = nnce(i) + 1, nnce(i) + nce(i)
+          l = l + iorbs(icvir(j))
+        end do
+        k = k + l
+      end do
     else
       call geochk()
     end if
diff --git a/src/mopend.F90 b/src/mopend.F90
index eae64815..f21c0515 100644
--- a/src/mopend.F90
+++ b/src/mopend.F90
@@ -42,7 +42,7 @@ subroutine summary(txt, ntxt)
 !
 !
         use chanel_C, only : iw, ir
-        use molkst_C, only : line, job_no, natoms
+        use molkst_C, only : line, job_no, natoms, dummy, errtxt
         implicit none
         integer, intent (in) :: ntxt
         character, intent (in) :: txt*(*)
@@ -54,6 +54,23 @@ subroutine summary(txt, ntxt)
         integer :: nmessages = 0, i, j, max_txt
         logical :: first = .true., opend
         save
+! vvv MOPAC API hijacking of the error handler vvv
+        ! report number of errors in dummy
+        if (ntxt == 0) then
+          dummy = nmessages
+          return
+        ! record error message in errtxt
+        else if (ntxt < 0) then
+          ! reset error handler
+          if (ntxt < -nmessages) then
+            first = .true.
+            nmessages = 0
+            return
+          end if
+          errtxt = trim(messages(-ntxt))
+          return
+        end if
+! ^^^ MOPAC API hijacking of the error handler ^^^
         if (first) then
            messages(1)(:18) = "JOB ENDED NORMALLY"
            first = .false.
diff --git a/tests/mopac_api_test.F90 b/tests/mopac_api_test.F90
index 92f14f10..2238cb18 100644
--- a/tests/mopac_api_test.F90
+++ b/tests/mopac_api_test.F90
@@ -26,6 +26,10 @@ program mopac_api_test
     call test_mozyme_scf1(nfail)
     call test_mozyme_relax1(nfail)
     call test_mozyme_vibe1(nfail)
+    call test_mopac_restart1(nfail)
+    call test_mozyme_restart1(nfail)
+    call test_cosmo1(nfail)
+    call test_crystal1(nfail)
     call exit(nfail)
 end program mopac_api_test
 
@@ -37,7 +41,7 @@ subroutine test_mopac_scf1(nfail)
     type(mopac_state) :: test_restore
     type(mopac_properties) :: test_target
     type(mopac_properties) :: test_out
-    character(20) :: test_name
+    character(50) :: test_name
     integer :: i
     test_name = 'H2O SCF'
 
@@ -101,7 +105,7 @@ subroutine test_mopac_scf1(nfail)
     test_target%bond_order(6) = 0.895d0
     test_target%bond_order(7) = 1.791d0
     test_target%calc_vibe = .false.
-    test_target%status = 0
+    test_target%nerror = 0
 
     call mopac_scf(test_in, test_restore, test_out)
     call test_output(test_name, test_in, test_target, test_out, nfail)
@@ -115,7 +119,7 @@ subroutine test_mopac_relax1(nfail)
     type(mopac_state) :: test_restore
     type(mopac_properties) :: test_target
     type(mopac_properties) :: test_out
-    character(20) :: test_name
+    character(50) :: test_name
     integer :: i
     test_name = 'H2O relax'
 
@@ -187,7 +191,7 @@ subroutine test_mopac_relax1(nfail)
     test_target%bond_order(6) = 0.894d0
     test_target%bond_order(7) = 1.788d0
     test_target%calc_vibe = .false.
-    test_target%status = 0
+    test_target%nerror = 0
 
     call mopac_relax(test_in, test_restore, test_out)
     call test_output(test_name, test_in, test_target, test_out, nfail)
@@ -201,7 +205,7 @@ subroutine test_mopac_vibe1(nfail)
     type(mopac_state) :: test_restore
     type(mopac_properties) :: test_target
     type(mopac_properties) :: test_out
-    character(20) :: test_name
+    character(50) :: test_name
     test_name = 'H2O vibe'
 
     ! 1 - geometry relaxation of H2O
@@ -264,7 +268,7 @@ subroutine test_mopac_vibe1(nfail)
     test_target%bond_order(6) = 0.894d0
     test_target%bond_order(7) = 1.788d0
     test_target%calc_vibe = .true.
-    test_target%status = 0
+    test_target%nerror = 0
     allocate(test_target%freq(9))
     test_target%freq(1) = -18.0d0
     test_target%freq(2) = -14.6d0
@@ -370,7 +374,7 @@ subroutine test_mozyme_scf1(nfail)
     type(mozyme_state) :: test_restore
     type(mopac_properties) :: test_target
     type(mopac_properties) :: test_out
-    character(20) :: test_name
+    character(50) :: test_name
     integer :: i
     test_name = 'H2O SCF MOZYME'
 
@@ -434,7 +438,7 @@ subroutine test_mozyme_scf1(nfail)
     test_target%bond_order(6) = 0.896d0
     test_target%bond_order(7) = 1.793d0
     test_target%calc_vibe = .false.
-    test_target%status = 0
+    test_target%nerror = 0
 
     call mozyme_scf(test_in, test_restore, test_out)
     call test_output(test_name, test_in, test_target, test_out, nfail)
@@ -448,7 +452,7 @@ subroutine test_mozyme_relax1(nfail)
     type(mozyme_state) :: test_restore
     type(mopac_properties) :: test_target
     type(mopac_properties) :: test_out
-    character(20) :: test_name
+    character(50) :: test_name
     integer :: i
     test_name = 'H2O relax MOZYME'
 
@@ -520,7 +524,7 @@ subroutine test_mozyme_relax1(nfail)
     test_target%bond_order(6) = 0.894d0
     test_target%bond_order(7) = 1.788d0
     test_target%calc_vibe = .false.
-    test_target%status = 0
+    test_target%nerror = 0
 
     call mozyme_relax(test_in, test_restore, test_out)
     call test_output(test_name, test_in, test_target, test_out, nfail)
@@ -535,7 +539,7 @@ subroutine test_mozyme_vibe1(nfail)
     type(mozyme_state) :: test_restore
     type(mopac_properties) :: test_target
     type(mopac_properties) :: test_out
-    character(20) :: test_name
+    character(50) :: test_name
     test_name = 'H2O vibe MOZYME'
 
     ! 1 - geometry relaxation of H2O
@@ -599,7 +603,7 @@ subroutine test_mozyme_vibe1(nfail)
     test_target%bond_order(6) = 0.894d0
     test_target%bond_order(7) = 1.788d0
     test_target%calc_vibe = .true.
-    test_target%status = 0
+    test_target%nerror = 0
     allocate(test_target%freq(9))
     test_target%freq(1) = -18.0d0
     test_target%freq(2) = -14.6d0
@@ -697,10 +701,1505 @@ subroutine test_mozyme_vibe1(nfail)
     call test_output(test_name, test_in, test_target, test_out, nfail)
 end subroutine test_mozyme_vibe1
 
+subroutine test_mopac_restart1(nfail)
+    use mopac_api
+    implicit none
+    integer, intent(inout) :: nfail
+    type(mopac_system) :: test_in
+    type(mopac_state) :: test_restore
+    type(mopac_properties) :: test_target
+    type(mopac_properties) :: test_out
+    character(50) :: test_name
+    integer :: i
+    test_name = 'H2O SCF restart'
+
+    ! SCF calculation of H2O
+    test_in%natom = 3
+    test_in%natom_move = 3
+    allocate(test_in%atom(3))
+    test_in%atom(1) = 1
+    test_in%atom(2) = 1
+    test_in%atom(3) = 8
+    allocate(test_in%coord(3*3))
+    test_in%coord(1) = 0.76d0
+    test_in%coord(2) = 0.59d0
+    test_in%coord(3) = 0.0d0
+    test_in%coord(4) = -0.76d0
+    test_in%coord(5) = 0.59d0
+    test_in%coord(6) = 0.0d0
+    test_in%coord(7) = 0.0d0
+    test_in%coord(8) = 0.0d0
+    test_in%coord(9) = 0.0d0
+    test_target%heat = -57.76975d0
+    allocate(test_target%coord_update(3*3))
+    test_target%coord_update = test_in%coord
+    allocate(test_target%coord_deriv(3*3))
+    test_target%coord_deriv(1) = 2.307865d0
+    test_target%coord_deriv(2) = 2.742432d0
+    test_target%coord_deriv(3) = 0.0d0
+    test_target%coord_deriv(4) = -2.307865d0
+    test_target%coord_deriv(5) = 2.711610d0
+    test_target%coord_deriv(6) = 0.0d0
+    test_target%coord_deriv(7) = 0.0d0
+    test_target%coord_deriv(8) = -5.454042d0
+    test_target%coord_deriv(9) = 0.0d0
+    allocate(test_target%charge(3))
+    test_target%charge(1) = 0.322260d0
+    test_target%charge(2) = 0.322260d0
+    test_target%charge(3) = -0.644520d0
+    test_target%dipole(1) = 0.0d0
+    test_target%dipole(2) = 2.147d0
+    test_target%dipole(3) = 0.0d0
+    test_target%stress(:) = 0.0d0
+    allocate(test_target%bond_index(4))
+    test_target%bond_index(1) = 1
+    test_target%bond_index(2) = 3
+    test_target%bond_index(3) = 5
+    test_target%bond_index(4) = 8
+    allocate(test_target%bond_atom(7))
+    test_target%bond_atom(1) = 1
+    test_target%bond_atom(2) = 3
+    test_target%bond_atom(3) = 2
+    test_target%bond_atom(4) = 3
+    test_target%bond_atom(5) = 1
+    test_target%bond_atom(6) = 2
+    test_target%bond_atom(7) = 3
+    allocate(test_target%bond_order(7))
+    test_target%bond_order(1) = 0.896d0
+    test_target%bond_order(2) = 0.895d0
+    test_target%bond_order(3) = 0.896d0
+    test_target%bond_order(4) = 0.895d0
+    test_target%bond_order(5) = 0.895d0
+    test_target%bond_order(6) = 0.895d0
+    test_target%bond_order(7) = 1.791d0
+    test_target%calc_vibe = .false.
+    test_target%nerror = 0
+
+    call mopac_scf(test_in, test_restore, test_out)
+    call mopac_scf(test_in, test_restore, test_out)
+    call test_output(test_name, test_in, test_target, test_out, nfail)
+end subroutine test_mopac_restart1
+
+subroutine test_mozyme_restart1(nfail)
+    use mopac_api
+    implicit none
+    integer, intent(inout) :: nfail
+    type(mopac_system) :: test_in
+    type(mozyme_state) :: test_restore
+    type(mopac_properties) :: test_target
+    type(mopac_properties) :: test_out
+    character(50) :: test_name
+    integer :: i
+    test_name = 'H2O SCF MOZYME restart'
+
+    ! SCF calculation of H2O
+    test_in%natom = 3
+    test_in%natom_move = 3
+    ! API-based save/load is captured at slightly different point in SCF cycle,
+    ! needs tigher tolerance for consistent results
+    test_in%tolerance = 0.1d0
+    allocate(test_in%atom(3))
+    test_in%atom(1) = 1
+    test_in%atom(2) = 1
+    test_in%atom(3) = 8
+    allocate(test_in%coord(3*3))
+    test_in%coord(1) = 0.76d0
+    test_in%coord(2) = 0.59d0
+    test_in%coord(3) = 0.0d0
+    test_in%coord(4) = -0.76d0
+    test_in%coord(5) = 0.59d0
+    test_in%coord(6) = 0.0d0
+    test_in%coord(7) = 0.0d0
+    test_in%coord(8) = 0.0d0
+    test_in%coord(9) = 0.0d0
+    test_target%heat = -57.76970d0
+    allocate(test_target%coord_update(3*3))
+    test_target%coord_update = test_in%coord
+    allocate(test_target%coord_deriv(3*3))
+    test_target%coord_deriv(1) = 2.296d0
+    test_target%coord_deriv(2) = 2.662d0
+    test_target%coord_deriv(3) = 0.001d0
+    test_target%coord_deriv(4) = -2.295d0
+    test_target%coord_deriv(5) = 2.631d0
+    test_target%coord_deriv(6) = 0.001d0
+    test_target%coord_deriv(7) = -0.001d0
+    test_target%coord_deriv(8) = -5.293d0
+    test_target%coord_deriv(9) = -0.001d0
+    allocate(test_target%charge(3))
+    test_target%charge(1) = 0.321963d0
+    test_target%charge(2) = 0.321963d0
+    test_target%charge(3) = -0.643926d0
+    test_target%dipole(1) = 0.0d0
+    test_target%dipole(2) = 2.146d0
+    test_target%dipole(3) = 0.0d0
+    test_target%stress(:) = 0.0d0
+    allocate(test_target%bond_index(4))
+    test_target%bond_index(1) = 1
+    test_target%bond_index(2) = 3
+    test_target%bond_index(3) = 5
+    test_target%bond_index(4) = 8
+    allocate(test_target%bond_atom(7))
+    test_target%bond_atom(1) = 1
+    test_target%bond_atom(2) = 3
+    test_target%bond_atom(3) = 2
+    test_target%bond_atom(4) = 3
+    test_target%bond_atom(5) = 1
+    test_target%bond_atom(6) = 2
+    test_target%bond_atom(7) = 3
+    allocate(test_target%bond_order(7))
+    test_target%bond_order(1) = 0.896d0
+    test_target%bond_order(2) = 0.895d0
+    test_target%bond_order(3) = 0.896d0
+    test_target%bond_order(4) = 0.895d0
+    test_target%bond_order(5) = 0.895d0
+    test_target%bond_order(6) = 0.895d0
+    test_target%bond_order(7) = 1.791d0
+    test_target%calc_vibe = .false.
+    test_target%nerror = 0
+
+    call mozyme_scf(test_in, test_restore, test_out)
+    call mozyme_scf(test_in, test_restore, test_out)
+    call test_output(test_name, test_in, test_target, test_out, nfail)
+end subroutine test_mozyme_restart1
+
+subroutine test_cosmo1(nfail)
+    use mopac_api
+    implicit none
+    integer, intent(inout) :: nfail
+    type(mopac_system) :: test_in
+    type(mopac_state) :: test_restore
+    type(mopac_properties) :: test_target
+    type(mopac_properties) :: test_out
+    character(50) :: test_name
+    integer :: i
+    test_name = 'H2O SCF COSMO'
+
+    ! SCF calculation of H2O
+    test_in%natom = 3
+    test_in%natom_move = 3
+    test_in%epsilon = 78.4d0
+    allocate(test_in%atom(3))
+    test_in%atom(1) = 1
+    test_in%atom(2) = 1
+    test_in%atom(3) = 8
+    allocate(test_in%coord(3*3))
+    test_in%coord(1) = 0.76d0
+    test_in%coord(2) = 0.59d0
+    test_in%coord(3) = 0.0d0
+    test_in%coord(4) = -0.76d0
+    test_in%coord(5) = 0.59d0
+    test_in%coord(6) = 0.0d0
+    test_in%coord(7) = 0.0d0
+    test_in%coord(8) = 0.0d0
+    test_in%coord(9) = 0.0d0
+    test_target%heat = -65.21458d0
+    allocate(test_target%coord_update(3*3))
+    test_target%coord_update = test_in%coord
+    allocate(test_target%coord_deriv(3*3))
+    test_target%coord_deriv(1) = -0.773696d0
+    test_target%coord_deriv(2) = -1.302274d0
+    test_target%coord_deriv(3) = 0.045812d0
+    test_target%coord_deriv(4) = 0.734278d0
+    test_target%coord_deriv(5) = -1.261587d0
+    test_target%coord_deriv(6) = -0.029041d0
+    test_target%coord_deriv(7) = 0.039417d0
+    test_target%coord_deriv(8) = 2.563860d0
+    test_target%coord_deriv(9) = -0.016772d0
+    allocate(test_target%charge(3))
+    test_target%charge(1) = 0.374234d0
+    test_target%charge(2) = 0.374006d0
+    test_target%charge(3) = -0.748240d0
+    test_target%dipole(1) = 0.0d0
+    test_target%dipole(2) = 2.43d0
+    test_target%dipole(3) = 0.0d0
+    test_target%stress(:) = 0.0d0
+    allocate(test_target%bond_index(4))
+    test_target%bond_index(1) = 1
+    test_target%bond_index(2) = 3
+    test_target%bond_index(3) = 5
+    test_target%bond_index(4) = 8
+    allocate(test_target%bond_atom(7))
+    test_target%bond_atom(1) = 1
+    test_target%bond_atom(2) = 3
+    test_target%bond_atom(3) = 2
+    test_target%bond_atom(4) = 3
+    test_target%bond_atom(5) = 1
+    test_target%bond_atom(6) = 2
+    test_target%bond_atom(7) = 3
+    allocate(test_target%bond_order(7))
+    test_target%bond_order(1) = 0.860d0
+    test_target%bond_order(2) = 0.859d0
+    test_target%bond_order(3) = 0.860d0
+    test_target%bond_order(4) = 0.859d0
+    test_target%bond_order(5) = 0.859d0
+    test_target%bond_order(6) = 0.859d0
+    test_target%bond_order(7) = 1.718d0
+    test_target%calc_vibe = .false.
+    test_target%nerror = 0
+
+    call mopac_scf(test_in, test_restore, test_out)
+    call test_output(test_name, test_in, test_target, test_out, nfail)
+end subroutine test_cosmo1
+
+
+subroutine test_crystal1(nfail)
+    use mopac_api
+    implicit none
+    integer, intent(inout) :: nfail
+    type(mopac_system) :: test_in
+    type(mopac_state) :: test_restore
+    type(mopac_properties) :: test_target
+    type(mopac_properties) :: test_out
+    character(50) :: test_name
+    integer :: i
+    test_name = 'LiF SCF'
+
+    ! SCF calculation of LiF
+    test_in%natom = 54
+    test_in%natom_move = 54
+    test_in%nlattice = 3
+    test_in%nlattice_move = 3
+    allocate(test_in%atom(54))
+    test_in%atom(1) = 3
+    test_in%atom(2) = 9
+    test_in%atom(3) = 3
+    test_in%atom(4) = 9
+    test_in%atom(5) = 3
+    test_in%atom(6) = 9
+    test_in%atom(7) = 3
+    test_in%atom(8) = 9
+    test_in%atom(9) = 3
+    test_in%atom(10) = 9
+    test_in%atom(11) = 3
+    test_in%atom(12) = 9
+    test_in%atom(13) = 3
+    test_in%atom(14) = 9
+    test_in%atom(15) = 3
+    test_in%atom(16) = 9
+    test_in%atom(17) = 3
+    test_in%atom(18) = 9
+    test_in%atom(19) = 3
+    test_in%atom(20) = 9
+    test_in%atom(21) = 3
+    test_in%atom(22) = 9
+    test_in%atom(23) = 3
+    test_in%atom(24) = 9
+    test_in%atom(25) = 3
+    test_in%atom(26) = 9
+    test_in%atom(27) = 3
+    test_in%atom(28) = 9
+    test_in%atom(29) = 3
+    test_in%atom(30) = 9
+    test_in%atom(31) = 3
+    test_in%atom(32) = 9
+    test_in%atom(33) = 3
+    test_in%atom(34) = 9
+    test_in%atom(35) = 3
+    test_in%atom(36) = 9
+    test_in%atom(37) = 3
+    test_in%atom(38) = 9
+    test_in%atom(39) = 3
+    test_in%atom(40) = 9
+    test_in%atom(41) = 3
+    test_in%atom(42) = 9
+    test_in%atom(43) = 3
+    test_in%atom(44) = 9
+    test_in%atom(45) = 3
+    test_in%atom(46) = 9
+    test_in%atom(47) = 3
+    test_in%atom(48) = 9
+    test_in%atom(49) = 3
+    test_in%atom(50) = 9
+    test_in%atom(51) = 3
+    test_in%atom(52) = 9
+    test_in%atom(53) = 3
+    test_in%atom(54) = 9
+    allocate(test_in%coord(3*54))
+    test_in%coord(1) = 0.0d0
+    test_in%coord(2) = 0.0d0
+    test_in%coord(3) = 0.0d0
+    test_in%coord(4) = 1.795d0
+    test_in%coord(5) = 1.795d0
+    test_in%coord(6) = 1.795d0
+    test_in%coord(7) = 0.0d0
+    test_in%coord(8) = 1.795d0
+    test_in%coord(9) = 1.795d0
+    test_in%coord(10) = 1.795d0
+    test_in%coord(11) = 3.59d0
+    test_in%coord(12) = 3.59d0
+    test_in%coord(13) = 0.0d0
+    test_in%coord(14) = 3.59d0
+    test_in%coord(15) = 3.59d0
+    test_in%coord(16) = 1.795d0
+    test_in%coord(17) = 5.385d0
+    test_in%coord(18) = 5.385d0
+    test_in%coord(19) = 1.795d0
+    test_in%coord(20) = 0.0d0
+    test_in%coord(21) = 1.795d0
+    test_in%coord(22) = 3.59d0
+    test_in%coord(23) = 1.795d0
+    test_in%coord(24) = 3.59d0
+    test_in%coord(25) = 1.795d0
+    test_in%coord(26) = 1.795d0
+    test_in%coord(27) = 3.59d0
+    test_in%coord(28) = 3.59d0
+    test_in%coord(29) = 3.59d0
+    test_in%coord(30) = 5.385d0
+    test_in%coord(31) = 1.795d0
+    test_in%coord(32) = 3.59d0
+    test_in%coord(33) = 5.385d0
+    test_in%coord(34) = 3.59d0
+    test_in%coord(35) = 5.385d0
+    test_in%coord(36) = 7.18d0
+    test_in%coord(37) = 3.59d0
+    test_in%coord(38) = 0.0d0
+    test_in%coord(39) = 3.59d0
+    test_in%coord(40) = 5.385d0
+    test_in%coord(41) = 1.795d0
+    test_in%coord(42) = 5.385d0
+    test_in%coord(43) = 3.59d0
+    test_in%coord(44) = 1.795d0
+    test_in%coord(45) = 5.385d0
+    test_in%coord(46) = 5.385d0
+    test_in%coord(47) = 3.59d0
+    test_in%coord(48) = 7.18d0
+    test_in%coord(49) = 3.59d0
+    test_in%coord(50) = 3.59d0
+    test_in%coord(51) = 7.18d0
+    test_in%coord(52) = 5.385d0
+    test_in%coord(53) = 5.385d0
+    test_in%coord(54) = 8.975d0
+    test_in%coord(55) = 1.795d0
+    test_in%coord(56) = 1.795d0
+    test_in%coord(57) = 0.0d0
+    test_in%coord(58) = 3.59d0
+    test_in%coord(59) = 3.59d0
+    test_in%coord(60) = 1.795d0
+    test_in%coord(61) = 1.795d0
+    test_in%coord(62) = 3.59d0
+    test_in%coord(63) = 1.795d0
+    test_in%coord(64) = 3.59d0
+    test_in%coord(65) = 5.385d0
+    test_in%coord(66) = 3.59d0
+    test_in%coord(67) = 1.795d0
+    test_in%coord(68) = 5.385d0
+    test_in%coord(69) = 3.59d0
+    test_in%coord(70) = 3.59d0
+    test_in%coord(71) = 7.18d0
+    test_in%coord(72) = 5.385d0
+    test_in%coord(73) = 3.59d0
+    test_in%coord(74) = 1.795d0
+    test_in%coord(75) = 1.795d0
+    test_in%coord(76) = 5.385d0
+    test_in%coord(77) = 3.59d0
+    test_in%coord(78) = 3.59d0
+    test_in%coord(79) = 3.59d0
+    test_in%coord(80) = 3.59d0
+    test_in%coord(81) = 3.59d0
+    test_in%coord(82) = 5.385d0
+    test_in%coord(83) = 5.385d0
+    test_in%coord(84) = 5.385d0
+    test_in%coord(85) = 3.59d0
+    test_in%coord(86) = 5.385d0
+    test_in%coord(87) = 5.385d0
+    test_in%coord(88) = 5.385d0
+    test_in%coord(89) = 7.18d0
+    test_in%coord(90) = 7.18d0
+    test_in%coord(91) = 5.385d0
+    test_in%coord(92) = 1.795d0
+    test_in%coord(93) = 3.59d0
+    test_in%coord(94) = 7.18d0
+    test_in%coord(95) = 3.59d0
+    test_in%coord(96) = 5.385d0
+    test_in%coord(97) = 5.385d0
+    test_in%coord(98) = 3.59d0
+    test_in%coord(99) = 5.385d0
+    test_in%coord(100) = 7.18d0
+    test_in%coord(101) = 5.385d0
+    test_in%coord(102) = 7.18d0
+    test_in%coord(103) = 5.385d0
+    test_in%coord(104) = 5.385d0
+    test_in%coord(105) = 7.18d0
+    test_in%coord(106) = 7.18d0
+    test_in%coord(107) = 7.18d0
+    test_in%coord(108) = 8.975d0
+    test_in%coord(109) = 3.59d0
+    test_in%coord(110) = 3.59d0
+    test_in%coord(111) = 0.0d0
+    test_in%coord(112) = 5.385d0
+    test_in%coord(113) = 5.385d0
+    test_in%coord(114) = 1.795d0
+    test_in%coord(115) = 3.59d0
+    test_in%coord(116) = 5.385d0
+    test_in%coord(117) = 1.795d0
+    test_in%coord(118) = 5.385d0
+    test_in%coord(119) = 7.18d0
+    test_in%coord(120) = 3.59d0
+    test_in%coord(121) = 3.59d0
+    test_in%coord(122) = 7.18d0
+    test_in%coord(123) = 3.59d0
+    test_in%coord(124) = 5.385d0
+    test_in%coord(125) = 8.975d0
+    test_in%coord(126) = 5.385d0
+    test_in%coord(127) = 5.385d0
+    test_in%coord(128) = 3.59d0
+    test_in%coord(129) = 1.795d0
+    test_in%coord(130) = 7.18d0
+    test_in%coord(131) = 5.385d0
+    test_in%coord(132) = 3.59d0
+    test_in%coord(133) = 5.385d0
+    test_in%coord(134) = 5.385d0
+    test_in%coord(135) = 3.59d0
+    test_in%coord(136) = 7.18d0
+    test_in%coord(137) = 7.18d0
+    test_in%coord(138) = 5.385d0
+    test_in%coord(139) = 5.385d0
+    test_in%coord(140) = 7.18d0
+    test_in%coord(141) = 5.385d0
+    test_in%coord(142) = 7.18d0
+    test_in%coord(143) = 8.975d0
+    test_in%coord(144) = 7.18d0
+    test_in%coord(145) = 7.18d0
+    test_in%coord(146) = 3.59d0
+    test_in%coord(147) = 3.59d0
+    test_in%coord(148) = 8.975d0
+    test_in%coord(149) = 5.385d0
+    test_in%coord(150) = 5.385d0
+    test_in%coord(151) = 7.18d0
+    test_in%coord(152) = 5.385d0
+    test_in%coord(153) = 5.385d0
+    test_in%coord(154) = 8.975d0
+    test_in%coord(155) = 7.18d0
+    test_in%coord(156) = 7.18d0
+    test_in%coord(157) = 7.18d0
+    test_in%coord(158) = 7.18d0
+    test_in%coord(159) = 7.18d0
+    test_in%coord(160) = 8.975d0
+    test_in%coord(161) = 8.975d0
+    test_in%coord(162) = 8.975d0
+    allocate(test_in%lattice(3*3))
+    test_in%lattice(1) = 5.385d0
+    test_in%lattice(2) = 5.385d0
+    test_in%lattice(3) = 0.0d0
+    test_in%lattice(4) = 5.385d0
+    test_in%lattice(5) = 0.0d0
+    test_in%lattice(6) = 5.385d0
+    test_in%lattice(7) = 0.0d0
+    test_in%lattice(8) = 5.385d0
+    test_in%lattice(9) = 5.385d0
+    test_target%heat = -4588.63517d0
+    allocate(test_target%coord_update(3*54))
+    test_target%coord_update = test_in%coord
+    allocate(test_target%coord_deriv(3*54))
+    test_target%coord_deriv(1) = -0.003072d0
+    test_target%coord_deriv(2) = -0.001544d0
+    test_target%coord_deriv(3) = -0.000184d0
+    test_target%coord_deriv(4) = -0.002616d0
+    test_target%coord_deriv(5) = -0.000781d0
+    test_target%coord_deriv(6) = -0.000527d0
+    test_target%coord_deriv(7) = -0.003046d0
+    test_target%coord_deriv(8) = -0.003246d0
+    test_target%coord_deriv(9) = -0.001886d0
+    test_target%coord_deriv(10) = -0.003656d0
+    test_target%coord_deriv(11) = -0.002632d0
+    test_target%coord_deriv(12) = -0.002377d0
+    test_target%coord_deriv(13) = -0.004086d0
+    test_target%coord_deriv(14) = -0.004968d0
+    test_target%coord_deriv(15) = -0.003608d0
+    test_target%coord_deriv(16) = -0.003656d0
+    test_target%coord_deriv(17) = -0.004478d0
+    test_target%coord_deriv(18) = -0.004224d0
+    test_target%coord_deriv(19) = -0.003531d0
+    test_target%coord_deriv(20) = 0.002452d0
+    test_target%coord_deriv(21) = -0.001774d0
+    test_target%coord_deriv(22) = -0.000814d0
+    test_target%coord_deriv(23) = -0.000424d0
+    test_target%coord_deriv(24) = -0.001538d0
+    test_target%coord_deriv(25) = -0.003531d0
+    test_target%coord_deriv(26) = 0.000725d0
+    test_target%coord_deriv(27) = -0.003496d0
+    test_target%coord_deriv(28) = -0.001854d0
+    test_target%coord_deriv(29) = -0.002273d0
+    test_target%coord_deriv(30) = -0.003385d0
+    test_target%coord_deriv(31) = -0.004570d0
+    test_target%coord_deriv(32) = -0.000998d0
+    test_target%coord_deriv(33) = -0.005222d0
+    test_target%coord_deriv(34) = -0.001854d0
+    test_target%coord_deriv(35) = -0.004120d0
+    test_target%coord_deriv(36) = -0.005232d0
+    test_target%coord_deriv(37) = -0.000966d0
+    test_target%coord_deriv(38) = 0.000375d0
+    test_target%coord_deriv(39) = -0.001107d0
+    test_target%coord_deriv(40) = -0.000700d0
+    test_target%coord_deriv(41) = 0.003466d0
+    test_target%coord_deriv(42) = -0.002285d0
+    test_target%coord_deriv(43) = -0.000966d0
+    test_target%coord_deriv(44) = -0.001352d0
+    test_target%coord_deriv(45) = -0.002832d0
+    test_target%coord_deriv(46) = -0.001741d0
+    test_target%coord_deriv(47) = 0.001616d0
+    test_target%coord_deriv(48) = -0.004131d0
+    test_target%coord_deriv(49) = -0.002006d0
+    test_target%coord_deriv(50) = -0.003074d0
+    test_target%coord_deriv(51) = -0.004558d0
+    test_target%coord_deriv(52) = -0.001740d0
+    test_target%coord_deriv(53) = -0.000230d0
+    test_target%coord_deriv(54) = -0.005979d0
+    test_target%coord_deriv(55) = -0.003615d0
+    test_target%coord_deriv(56) = -0.003741d0
+    test_target%coord_deriv(57) = 0.009705d0
+    test_target%coord_deriv(58) = 0.000734d0
+    test_target%coord_deriv(59) = 0.001680d0
+    test_target%coord_deriv(60) = -0.002221d0
+    test_target%coord_deriv(61) = -0.003615d0
+    test_target%coord_deriv(62) = -0.005463d0
+    test_target%coord_deriv(63) = 0.007978d0
+    test_target%coord_deriv(64) = -0.000306d0
+    test_target%coord_deriv(65) = -0.000167d0
+    test_target%coord_deriv(66) = -0.004070d0
+    test_target%coord_deriv(67) = -0.004655d0
+    test_target%coord_deriv(68) = -0.007189d0
+    test_target%coord_deriv(69) = 0.006255d0
+    test_target%coord_deriv(70) = -0.000306d0
+    test_target%coord_deriv(71) = -0.002013d0
+    test_target%coord_deriv(72) = -0.005917d0
+    test_target%coord_deriv(73) = -0.004098d0
+    test_target%coord_deriv(74) = 0.000230d0
+    test_target%coord_deriv(75) = 0.008089d0
+    test_target%coord_deriv(76) = 0.002535d0
+    test_target%coord_deriv(77) = 0.002039d0
+    test_target%coord_deriv(78) = -0.003233d0
+    test_target%coord_deriv(79) = -0.004098d0
+    test_target%coord_deriv(80) = -0.001493d0
+    test_target%coord_deriv(81) = 0.006365d0
+    test_target%coord_deriv(82) = 0.001495d0
+    test_target%coord_deriv(83) = 0.000193d0
+    test_target%coord_deriv(84) = -0.005079d0
+    test_target%coord_deriv(85) = -0.005138d0
+    test_target%coord_deriv(86) = -0.003218d0
+    test_target%coord_deriv(87) = 0.004640d0
+    test_target%coord_deriv(88) = 0.001495d0
+    test_target%coord_deriv(89) = -0.001654d0
+    test_target%coord_deriv(90) = -0.006926d0
+    test_target%coord_deriv(91) = -0.001534d0
+    test_target%coord_deriv(92) = -0.001845d0
+    test_target%coord_deriv(93) = 0.008754d0
+    test_target%coord_deriv(94) = 0.002650d0
+    test_target%coord_deriv(95) = 0.005927d0
+    test_target%coord_deriv(96) = -0.003979d0
+    test_target%coord_deriv(97) = -0.001534d0
+    test_target%coord_deriv(98) = -0.003569d0
+    test_target%coord_deriv(99) = 0.007028d0
+    test_target%coord_deriv(100) = 0.001610d0
+    test_target%coord_deriv(101) = 0.004081d0
+    test_target%coord_deriv(102) = -0.005826d0
+    test_target%coord_deriv(103) = -0.002574d0
+    test_target%coord_deriv(104) = -0.005294d0
+    test_target%coord_deriv(105) = 0.005303d0
+    test_target%coord_deriv(106) = 0.001610d0
+    test_target%coord_deriv(107) = 0.002233d0
+    test_target%coord_deriv(108) = -0.007674d0
+    test_target%coord_deriv(109) = 0.003433d0
+    test_target%coord_deriv(110) = 0.003270d0
+    test_target%coord_deriv(111) = 0.002281d0
+    test_target%coord_deriv(112) = 0.002916d0
+    test_target%coord_deriv(113) = 0.001338d0
+    test_target%coord_deriv(114) = 0.006969d0
+    test_target%coord_deriv(115) = 0.003433d0
+    test_target%coord_deriv(116) = 0.001544d0
+    test_target%coord_deriv(117) = 0.000555d0
+    test_target%coord_deriv(118) = 0.001876d0
+    test_target%coord_deriv(119) = -0.000509d0
+    test_target%coord_deriv(120) = 0.005122d0
+    test_target%coord_deriv(121) = 0.002393d0
+    test_target%coord_deriv(122) = -0.000181d0
+    test_target%coord_deriv(123) = -0.001169d0
+    test_target%coord_deriv(124) = 0.001876d0
+    test_target%coord_deriv(125) = -0.002357d0
+    test_target%coord_deriv(126) = 0.003275d0
+    test_target%coord_deriv(127) = 0.002947d0
+    test_target%coord_deriv(128) = 0.007238d0
+    test_target%coord_deriv(129) = 0.000669d0
+    test_target%coord_deriv(130) = 0.004719d0
+    test_target%coord_deriv(131) = 0.001697d0
+    test_target%coord_deriv(132) = 0.005959d0
+    test_target%coord_deriv(133) = 0.002947d0
+    test_target%coord_deriv(134) = 0.005512d0
+    test_target%coord_deriv(135) = -0.001055d0
+    test_target%coord_deriv(136) = 0.003679d0
+    test_target%coord_deriv(137) = -0.000149d0
+    test_target%coord_deriv(138) = 0.004112d0
+    test_target%coord_deriv(139) = 0.001907d0
+    test_target%coord_deriv(140) = 0.003787d0
+    test_target%coord_deriv(141) = -0.002780d0
+    test_target%coord_deriv(142) = 0.003679d0
+    test_target%coord_deriv(143) = -0.001997d0
+    test_target%coord_deriv(144) = 0.002265d0
+    test_target%coord_deriv(145) = 0.005513d0
+    test_target%coord_deriv(146) = 0.005162d0
+    test_target%coord_deriv(147) = 0.001333d0
+    test_target%coord_deriv(148) = 0.004834d0
+    test_target%coord_deriv(149) = 0.005587d0
+    test_target%coord_deriv(150) = 0.005213d0
+    test_target%coord_deriv(151) = 0.005513d0
+    test_target%coord_deriv(152) = 0.003437d0
+    test_target%coord_deriv(153) = -0.000392d0
+    test_target%coord_deriv(154) = 0.003794d0
+    test_target%coord_deriv(155) = 0.003740d0
+    test_target%coord_deriv(156) = 0.003366d0
+    test_target%coord_deriv(157) = 0.004473d0
+    test_target%coord_deriv(158) = 0.001712d0
+    test_target%coord_deriv(159) = -0.002117d0
+    test_target%coord_deriv(160) = 0.003794d0
+    test_target%coord_deriv(161) = 0.001892d0
+    test_target%coord_deriv(162) = 0.001519d0
+    allocate(test_target%lattice_update(3*3))
+    test_target%lattice_update = test_in%lattice
+    allocate(test_target%lattice_deriv(3*3))
+    test_target%lattice_deriv(1) = -2.329702d0
+    test_target%lattice_deriv(2) = -2.324697d0
+    test_target%lattice_deriv(3) = 2.296600d0
+    test_target%lattice_deriv(4) = -2.314736d0
+    test_target%lattice_deriv(5) = 2.301304d0
+    test_target%lattice_deriv(6) = -2.299118d0
+    test_target%lattice_deriv(7) = 2.302052d0
+    test_target%lattice_deriv(8) = -2.296322d0
+    test_target%lattice_deriv(9) = -2.281867d0
+    allocate(test_target%charge(54))
+    test_target%charge(1) = 0.855660d0
+    test_target%charge(2) = -0.855642d0
+    test_target%charge(3) = 0.855660d0
+    test_target%charge(4) = -0.855642d0
+    test_target%charge(5) = 0.855660d0
+    test_target%charge(6) = -0.855642d0
+    test_target%charge(7) = 0.855654d0
+    test_target%charge(8) = -0.855644d0
+    test_target%charge(9) = 0.855654d0
+    test_target%charge(10) = -0.855644d0
+    test_target%charge(11) = 0.855654d0
+    test_target%charge(12) = -0.855644d0
+    test_target%charge(13) = 0.855654d0
+    test_target%charge(14) = -0.855642d0
+    test_target%charge(15) = 0.855654d0
+    test_target%charge(16) = -0.855642d0
+    test_target%charge(17) = 0.855654d0
+    test_target%charge(18) = -0.855642d0
+    test_target%charge(19) = 0.855642d0
+    test_target%charge(20) = -0.855647d0
+    test_target%charge(21) = 0.855642d0
+    test_target%charge(22) = -0.855647d0
+    test_target%charge(23) = 0.855642d0
+    test_target%charge(24) = -0.855647d0
+    test_target%charge(25) = 0.855637d0
+    test_target%charge(26) = -0.855649d0
+    test_target%charge(27) = 0.855637d0
+    test_target%charge(28) = -0.855649d0
+    test_target%charge(29) = 0.855637d0
+    test_target%charge(30) = -0.855649d0
+    test_target%charge(31) = 0.855636d0
+    test_target%charge(32) = -0.855647d0
+    test_target%charge(33) = 0.855636d0
+    test_target%charge(34) = -0.855647d0
+    test_target%charge(35) = 0.855636d0
+    test_target%charge(36) = -0.855647d0
+    test_target%charge(37) = 0.855642d0
+    test_target%charge(38) = -0.855641d0
+    test_target%charge(39) = 0.855642d0
+    test_target%charge(40) = -0.855641d0
+    test_target%charge(41) = 0.855642d0
+    test_target%charge(42) = -0.855641d0
+    test_target%charge(43) = 0.855636d0
+    test_target%charge(44) = -0.855643d0
+    test_target%charge(45) = 0.855636d0
+    test_target%charge(46) = -0.855643d0
+    test_target%charge(47) = 0.855636d0
+    test_target%charge(48) = -0.855643d0
+    test_target%charge(49) = 0.855636d0
+    test_target%charge(50) = -0.855641d0
+    test_target%charge(51) = 0.855636d0
+    test_target%charge(52) = -0.855641d0
+    test_target%charge(53) = 0.855636d0
+    test_target%charge(54) = -0.855641d0
+    test_target%dipole(:) = 0.0d0
+    test_target%stress(1) = -0.550d0
+    test_target%stress(2) = -0.553d0
+    test_target%stress(3) = -0.554d0
+    test_target%stress(4) = 0.001d0
+    test_target%stress(5) = 0.000d0
+    test_target%stress(6) = 0.002d0
+    allocate(test_target%bond_index(55))
+    do i = 1, 55
+        test_target%bond_index(i) = 1 + 7*(i-1)
+    end do
+    allocate(test_target%bond_atom(54*7))
+    test_target%bond_atom(1) = 1
+    test_target%bond_atom(2) = 6
+    test_target%bond_atom(3) = 14
+    test_target%bond_atom(4) = 18
+    test_target%bond_atom(5) = 38
+    test_target%bond_atom(6) = 42
+    test_target%bond_atom(7) = 50
+    test_target%bond_atom(8) = 2
+    test_target%bond_atom(9) = 3
+    test_target%bond_atom(10) = 7
+    test_target%bond_atom(11) = 9
+    test_target%bond_atom(12) = 19
+    test_target%bond_atom(13) = 21
+    test_target%bond_atom(14) = 25
+    test_target%bond_atom(15) = 2
+    test_target%bond_atom(16) = 3
+    test_target%bond_atom(17) = 14
+    test_target%bond_atom(18) = 16
+    test_target%bond_atom(19) = 38
+    test_target%bond_atom(20) = 40
+    test_target%bond_atom(21) = 52
+    test_target%bond_atom(22) = 4
+    test_target%bond_atom(23) = 5
+    test_target%bond_atom(24) = 9
+    test_target%bond_atom(25) = 11
+    test_target%bond_atom(26) = 21
+    test_target%bond_atom(27) = 23
+    test_target%bond_atom(28) = 27
+    test_target%bond_atom(29) = 4
+    test_target%bond_atom(30) = 5
+    test_target%bond_atom(31) = 16
+    test_target%bond_atom(32) = 18
+    test_target%bond_atom(33) = 40
+    test_target%bond_atom(34) = 42
+    test_target%bond_atom(35) = 54
+    test_target%bond_atom(36) = 1
+    test_target%bond_atom(37) = 6
+    test_target%bond_atom(38) = 7
+    test_target%bond_atom(39) = 11
+    test_target%bond_atom(40) = 19
+    test_target%bond_atom(41) = 23
+    test_target%bond_atom(42) = 29
+    test_target%bond_atom(43) = 2
+    test_target%bond_atom(44) = 6
+    test_target%bond_atom(45) = 7
+    test_target%bond_atom(46) = 12
+    test_target%bond_atom(47) = 38
+    test_target%bond_atom(48) = 44
+    test_target%bond_atom(49) = 48
+    test_target%bond_atom(50) = 8
+    test_target%bond_atom(51) = 9
+    test_target%bond_atom(52) = 13
+    test_target%bond_atom(53) = 15
+    test_target%bond_atom(54) = 25
+    test_target%bond_atom(55) = 27
+    test_target%bond_atom(56) = 31
+    test_target%bond_atom(57) = 2
+    test_target%bond_atom(58) = 4
+    test_target%bond_atom(59) = 8
+    test_target%bond_atom(60) = 9
+    test_target%bond_atom(61) = 40
+    test_target%bond_atom(62) = 44
+    test_target%bond_atom(63) = 46
+    test_target%bond_atom(64) = 10
+    test_target%bond_atom(65) = 11
+    test_target%bond_atom(66) = 15
+    test_target%bond_atom(67) = 17
+    test_target%bond_atom(68) = 27
+    test_target%bond_atom(69) = 29
+    test_target%bond_atom(70) = 33
+    test_target%bond_atom(71) = 4
+    test_target%bond_atom(72) = 6
+    test_target%bond_atom(73) = 10
+    test_target%bond_atom(74) = 11
+    test_target%bond_atom(75) = 42
+    test_target%bond_atom(76) = 46
+    test_target%bond_atom(77) = 48
+    test_target%bond_atom(78) = 7
+    test_target%bond_atom(79) = 12
+    test_target%bond_atom(80) = 13
+    test_target%bond_atom(81) = 17
+    test_target%bond_atom(82) = 25
+    test_target%bond_atom(83) = 29
+    test_target%bond_atom(84) = 35
+    test_target%bond_atom(85) = 8
+    test_target%bond_atom(86) = 12
+    test_target%bond_atom(87) = 13
+    test_target%bond_atom(88) = 18
+    test_target%bond_atom(89) = 44
+    test_target%bond_atom(90) = 50
+    test_target%bond_atom(91) = 54
+    test_target%bond_atom(92) = 1
+    test_target%bond_atom(93) = 3
+    test_target%bond_atom(94) = 14
+    test_target%bond_atom(95) = 15
+    test_target%bond_atom(96) = 19
+    test_target%bond_atom(97) = 31
+    test_target%bond_atom(98) = 33
+    test_target%bond_atom(99) = 8
+    test_target%bond_atom(100) = 10
+    test_target%bond_atom(101) = 14
+    test_target%bond_atom(102) = 15
+    test_target%bond_atom(103) = 46
+    test_target%bond_atom(104) = 50
+    test_target%bond_atom(105) = 52
+    test_target%bond_atom(106) = 3
+    test_target%bond_atom(107) = 5
+    test_target%bond_atom(108) = 16
+    test_target%bond_atom(109) = 17
+    test_target%bond_atom(110) = 21
+    test_target%bond_atom(111) = 33
+    test_target%bond_atom(112) = 35
+    test_target%bond_atom(113) = 10
+    test_target%bond_atom(114) = 12
+    test_target%bond_atom(115) = 16
+    test_target%bond_atom(116) = 17
+    test_target%bond_atom(117) = 48
+    test_target%bond_atom(118) = 52
+    test_target%bond_atom(119) = 54
+    test_target%bond_atom(120) = 1
+    test_target%bond_atom(121) = 5
+    test_target%bond_atom(122) = 13
+    test_target%bond_atom(123) = 18
+    test_target%bond_atom(124) = 23
+    test_target%bond_atom(125) = 31
+    test_target%bond_atom(126) = 35
+    test_target%bond_atom(127) = 2
+    test_target%bond_atom(128) = 6
+    test_target%bond_atom(129) = 14
+    test_target%bond_atom(130) = 19
+    test_target%bond_atom(131) = 24
+    test_target%bond_atom(132) = 32
+    test_target%bond_atom(133) = 36
+    test_target%bond_atom(134) = 20
+    test_target%bond_atom(135) = 21
+    test_target%bond_atom(136) = 25
+    test_target%bond_atom(137) = 27
+    test_target%bond_atom(138) = 37
+    test_target%bond_atom(139) = 39
+    test_target%bond_atom(140) = 43
+    test_target%bond_atom(141) = 2
+    test_target%bond_atom(142) = 4
+    test_target%bond_atom(143) = 16
+    test_target%bond_atom(144) = 20
+    test_target%bond_atom(145) = 21
+    test_target%bond_atom(146) = 32
+    test_target%bond_atom(147) = 34
+    test_target%bond_atom(148) = 22
+    test_target%bond_atom(149) = 23
+    test_target%bond_atom(150) = 27
+    test_target%bond_atom(151) = 29
+    test_target%bond_atom(152) = 39
+    test_target%bond_atom(153) = 41
+    test_target%bond_atom(154) = 45
+    test_target%bond_atom(155) = 4
+    test_target%bond_atom(156) = 6
+    test_target%bond_atom(157) = 18
+    test_target%bond_atom(158) = 22
+    test_target%bond_atom(159) = 23
+    test_target%bond_atom(160) = 34
+    test_target%bond_atom(161) = 36
+    test_target%bond_atom(162) = 19
+    test_target%bond_atom(163) = 24
+    test_target%bond_atom(164) = 25
+    test_target%bond_atom(165) = 29
+    test_target%bond_atom(166) = 37
+    test_target%bond_atom(167) = 41
+    test_target%bond_atom(168) = 47
+    test_target%bond_atom(169) = 2
+    test_target%bond_atom(170) = 8
+    test_target%bond_atom(171) = 12
+    test_target%bond_atom(172) = 20
+    test_target%bond_atom(173) = 24
+    test_target%bond_atom(174) = 25
+    test_target%bond_atom(175) = 30
+    test_target%bond_atom(176) = 26
+    test_target%bond_atom(177) = 27
+    test_target%bond_atom(178) = 31
+    test_target%bond_atom(179) = 33
+    test_target%bond_atom(180) = 43
+    test_target%bond_atom(181) = 45
+    test_target%bond_atom(182) = 49
+    test_target%bond_atom(183) = 4
+    test_target%bond_atom(184) = 8
+    test_target%bond_atom(185) = 10
+    test_target%bond_atom(186) = 20
+    test_target%bond_atom(187) = 22
+    test_target%bond_atom(188) = 26
+    test_target%bond_atom(189) = 27
+    test_target%bond_atom(190) = 28
+    test_target%bond_atom(191) = 29
+    test_target%bond_atom(192) = 33
+    test_target%bond_atom(193) = 35
+    test_target%bond_atom(194) = 45
+    test_target%bond_atom(195) = 47
+    test_target%bond_atom(196) = 51
+    test_target%bond_atom(197) = 6
+    test_target%bond_atom(198) = 10
+    test_target%bond_atom(199) = 12
+    test_target%bond_atom(200) = 22
+    test_target%bond_atom(201) = 24
+    test_target%bond_atom(202) = 28
+    test_target%bond_atom(203) = 29
+    test_target%bond_atom(204) = 25
+    test_target%bond_atom(205) = 30
+    test_target%bond_atom(206) = 31
+    test_target%bond_atom(207) = 35
+    test_target%bond_atom(208) = 43
+    test_target%bond_atom(209) = 47
+    test_target%bond_atom(210) = 53
+    test_target%bond_atom(211) = 8
+    test_target%bond_atom(212) = 14
+    test_target%bond_atom(213) = 18
+    test_target%bond_atom(214) = 26
+    test_target%bond_atom(215) = 30
+    test_target%bond_atom(216) = 31
+    test_target%bond_atom(217) = 36
+    test_target%bond_atom(218) = 19
+    test_target%bond_atom(219) = 21
+    test_target%bond_atom(220) = 32
+    test_target%bond_atom(221) = 33
+    test_target%bond_atom(222) = 37
+    test_target%bond_atom(223) = 49
+    test_target%bond_atom(224) = 51
+    test_target%bond_atom(225) = 10
+    test_target%bond_atom(226) = 14
+    test_target%bond_atom(227) = 16
+    test_target%bond_atom(228) = 26
+    test_target%bond_atom(229) = 28
+    test_target%bond_atom(230) = 32
+    test_target%bond_atom(231) = 33
+    test_target%bond_atom(232) = 21
+    test_target%bond_atom(233) = 23
+    test_target%bond_atom(234) = 34
+    test_target%bond_atom(235) = 35
+    test_target%bond_atom(236) = 39
+    test_target%bond_atom(237) = 51
+    test_target%bond_atom(238) = 53
+    test_target%bond_atom(239) = 12
+    test_target%bond_atom(240) = 16
+    test_target%bond_atom(241) = 18
+    test_target%bond_atom(242) = 28
+    test_target%bond_atom(243) = 30
+    test_target%bond_atom(244) = 34
+    test_target%bond_atom(245) = 35
+    test_target%bond_atom(246) = 19
+    test_target%bond_atom(247) = 23
+    test_target%bond_atom(248) = 31
+    test_target%bond_atom(249) = 36
+    test_target%bond_atom(250) = 41
+    test_target%bond_atom(251) = 49
+    test_target%bond_atom(252) = 53
+    test_target%bond_atom(253) = 20
+    test_target%bond_atom(254) = 24
+    test_target%bond_atom(255) = 32
+    test_target%bond_atom(256) = 37
+    test_target%bond_atom(257) = 42
+    test_target%bond_atom(258) = 50
+    test_target%bond_atom(259) = 54
+    test_target%bond_atom(260) = 1
+    test_target%bond_atom(261) = 3
+    test_target%bond_atom(262) = 7
+    test_target%bond_atom(263) = 38
+    test_target%bond_atom(264) = 39
+    test_target%bond_atom(265) = 43
+    test_target%bond_atom(266) = 45
+    test_target%bond_atom(267) = 20
+    test_target%bond_atom(268) = 22
+    test_target%bond_atom(269) = 34
+    test_target%bond_atom(270) = 38
+    test_target%bond_atom(271) = 39
+    test_target%bond_atom(272) = 50
+    test_target%bond_atom(273) = 52
+    test_target%bond_atom(274) = 3
+    test_target%bond_atom(275) = 5
+    test_target%bond_atom(276) = 9
+    test_target%bond_atom(277) = 40
+    test_target%bond_atom(278) = 41
+    test_target%bond_atom(279) = 45
+    test_target%bond_atom(280) = 47
+    test_target%bond_atom(281) = 22
+    test_target%bond_atom(282) = 24
+    test_target%bond_atom(283) = 36
+    test_target%bond_atom(284) = 40
+    test_target%bond_atom(285) = 41
+    test_target%bond_atom(286) = 52
+    test_target%bond_atom(287) = 54
+    test_target%bond_atom(288) = 1
+    test_target%bond_atom(289) = 5
+    test_target%bond_atom(290) = 11
+    test_target%bond_atom(291) = 37
+    test_target%bond_atom(292) = 42
+    test_target%bond_atom(293) = 43
+    test_target%bond_atom(294) = 47
+    test_target%bond_atom(295) = 20
+    test_target%bond_atom(296) = 26
+    test_target%bond_atom(297) = 30
+    test_target%bond_atom(298) = 38
+    test_target%bond_atom(299) = 42
+    test_target%bond_atom(300) = 43
+    test_target%bond_atom(301) = 48
+    test_target%bond_atom(302) = 7
+    test_target%bond_atom(303) = 9
+    test_target%bond_atom(304) = 13
+    test_target%bond_atom(305) = 44
+    test_target%bond_atom(306) = 45
+    test_target%bond_atom(307) = 49
+    test_target%bond_atom(308) = 51
+    test_target%bond_atom(309) = 22
+    test_target%bond_atom(310) = 26
+    test_target%bond_atom(311) = 28
+    test_target%bond_atom(312) = 38
+    test_target%bond_atom(313) = 40
+    test_target%bond_atom(314) = 44
+    test_target%bond_atom(315) = 45
+    test_target%bond_atom(316) = 9
+    test_target%bond_atom(317) = 11
+    test_target%bond_atom(318) = 15
+    test_target%bond_atom(319) = 46
+    test_target%bond_atom(320) = 47
+    test_target%bond_atom(321) = 51
+    test_target%bond_atom(322) = 53
+    test_target%bond_atom(323) = 24
+    test_target%bond_atom(324) = 28
+    test_target%bond_atom(325) = 30
+    test_target%bond_atom(326) = 40
+    test_target%bond_atom(327) = 42
+    test_target%bond_atom(328) = 46
+    test_target%bond_atom(329) = 47
+    test_target%bond_atom(330) = 7
+    test_target%bond_atom(331) = 11
+    test_target%bond_atom(332) = 17
+    test_target%bond_atom(333) = 43
+    test_target%bond_atom(334) = 48
+    test_target%bond_atom(335) = 49
+    test_target%bond_atom(336) = 53
+    test_target%bond_atom(337) = 26
+    test_target%bond_atom(338) = 32
+    test_target%bond_atom(339) = 36
+    test_target%bond_atom(340) = 44
+    test_target%bond_atom(341) = 48
+    test_target%bond_atom(342) = 49
+    test_target%bond_atom(343) = 54
+    test_target%bond_atom(344) = 1
+    test_target%bond_atom(345) = 13
+    test_target%bond_atom(346) = 15
+    test_target%bond_atom(347) = 37
+    test_target%bond_atom(348) = 39
+    test_target%bond_atom(349) = 50
+    test_target%bond_atom(350) = 51
+    test_target%bond_atom(351) = 28
+    test_target%bond_atom(352) = 32
+    test_target%bond_atom(353) = 34
+    test_target%bond_atom(354) = 44
+    test_target%bond_atom(355) = 46
+    test_target%bond_atom(356) = 50
+    test_target%bond_atom(357) = 51
+    test_target%bond_atom(358) = 3
+    test_target%bond_atom(359) = 15
+    test_target%bond_atom(360) = 17
+    test_target%bond_atom(361) = 39
+    test_target%bond_atom(362) = 41
+    test_target%bond_atom(363) = 52
+    test_target%bond_atom(364) = 53
+    test_target%bond_atom(365) = 30
+    test_target%bond_atom(366) = 34
+    test_target%bond_atom(367) = 36
+    test_target%bond_atom(368) = 46
+    test_target%bond_atom(369) = 48
+    test_target%bond_atom(370) = 52
+    test_target%bond_atom(371) = 53
+    test_target%bond_atom(372) = 5
+    test_target%bond_atom(373) = 13
+    test_target%bond_atom(374) = 17
+    test_target%bond_atom(375) = 37
+    test_target%bond_atom(376) = 41
+    test_target%bond_atom(377) = 49
+    test_target%bond_atom(378) = 54
+    allocate(test_target%bond_order(54*7))
+    test_target%bond_order(1) = 0.268d0
+    test_target%bond_order(2) = 0.037d0
+    test_target%bond_order(3) = 0.037d0
+    test_target%bond_order(4) = 0.037d0
+    test_target%bond_order(5) = 0.037d0
+    test_target%bond_order(6) = 0.037d0
+    test_target%bond_order(7) = 0.037d0
+    test_target%bond_order(8) = 0.283d0
+    test_target%bond_order(9) = 0.037d0
+    test_target%bond_order(10) = 0.037d0
+    test_target%bond_order(11) = 0.037d0
+    test_target%bond_order(12) = 0.037d0
+    test_target%bond_order(13) = 0.037d0
+    test_target%bond_order(14) = 0.037d0
+    test_target%bond_order(15) = 0.037d0
+    test_target%bond_order(16) = 0.268d0
+    test_target%bond_order(17) = 0.037d0
+    test_target%bond_order(18) = 0.037d0
+    test_target%bond_order(19) = 0.037d0
+    test_target%bond_order(20) = 0.037d0
+    test_target%bond_order(21) = 0.037d0
+    test_target%bond_order(22) = 0.283d0
+    test_target%bond_order(23) = 0.037d0
+    test_target%bond_order(24) = 0.037d0
+    test_target%bond_order(25) = 0.037d0
+    test_target%bond_order(26) = 0.037d0
+    test_target%bond_order(27) = 0.037d0
+    test_target%bond_order(28) = 0.037d0
+    test_target%bond_order(29) = 0.037d0
+    test_target%bond_order(30) = 0.268d0
+    test_target%bond_order(31) = 0.037d0
+    test_target%bond_order(32) = 0.037d0
+    test_target%bond_order(33) = 0.037d0
+    test_target%bond_order(34) = 0.037d0
+    test_target%bond_order(35) = 0.037d0
+    test_target%bond_order(36) = 0.037d0
+    test_target%bond_order(37) = 0.283d0
+    test_target%bond_order(38) = 0.037d0
+    test_target%bond_order(39) = 0.037d0
+    test_target%bond_order(40) = 0.037d0
+    test_target%bond_order(41) = 0.037d0
+    test_target%bond_order(42) = 0.037d0
+    test_target%bond_order(43) = 0.037d0
+    test_target%bond_order(44) = 0.037d0
+    test_target%bond_order(45) = 0.268d0
+    test_target%bond_order(46) = 0.037d0
+    test_target%bond_order(47) = 0.037d0
+    test_target%bond_order(48) = 0.037d0
+    test_target%bond_order(49) = 0.037d0
+    test_target%bond_order(50) = 0.283d0
+    test_target%bond_order(51) = 0.037d0
+    test_target%bond_order(52) = 0.037d0
+    test_target%bond_order(53) = 0.037d0
+    test_target%bond_order(54) = 0.037d0
+    test_target%bond_order(55) = 0.037d0
+    test_target%bond_order(56) = 0.037d0
+    test_target%bond_order(57) = 0.037d0
+    test_target%bond_order(58) = 0.037d0
+    test_target%bond_order(59) = 0.037d0
+    test_target%bond_order(60) = 0.268d0
+    test_target%bond_order(61) = 0.037d0
+    test_target%bond_order(62) = 0.037d0
+    test_target%bond_order(63) = 0.037d0
+    test_target%bond_order(64) = 0.283d0
+    test_target%bond_order(65) = 0.037d0
+    test_target%bond_order(66) = 0.037d0
+    test_target%bond_order(67) = 0.037d0
+    test_target%bond_order(68) = 0.037d0
+    test_target%bond_order(69) = 0.037d0
+    test_target%bond_order(70) = 0.037d0
+    test_target%bond_order(71) = 0.037d0
+    test_target%bond_order(72) = 0.037d0
+    test_target%bond_order(73) = 0.037d0
+    test_target%bond_order(74) = 0.268d0
+    test_target%bond_order(75) = 0.037d0
+    test_target%bond_order(76) = 0.037d0
+    test_target%bond_order(77) = 0.037d0
+    test_target%bond_order(78) = 0.037d0
+    test_target%bond_order(79) = 0.283d0
+    test_target%bond_order(80) = 0.037d0
+    test_target%bond_order(81) = 0.037d0
+    test_target%bond_order(82) = 0.037d0
+    test_target%bond_order(83) = 0.037d0
+    test_target%bond_order(84) = 0.037d0
+    test_target%bond_order(85) = 0.037d0
+    test_target%bond_order(86) = 0.037d0
+    test_target%bond_order(87) = 0.268d0
+    test_target%bond_order(88) = 0.037d0
+    test_target%bond_order(89) = 0.037d0
+    test_target%bond_order(90) = 0.037d0
+    test_target%bond_order(91) = 0.037d0
+    test_target%bond_order(92) = 0.037d0
+    test_target%bond_order(93) = 0.037d0
+    test_target%bond_order(94) = 0.283d0
+    test_target%bond_order(95) = 0.037d0
+    test_target%bond_order(96) = 0.037d0
+    test_target%bond_order(97) = 0.037d0
+    test_target%bond_order(98) = 0.037d0
+    test_target%bond_order(99) = 0.037d0
+    test_target%bond_order(100) = 0.037d0
+    test_target%bond_order(101) = 0.037d0
+    test_target%bond_order(102) = 0.268d0
+    test_target%bond_order(103) = 0.037d0
+    test_target%bond_order(104) = 0.037d0
+    test_target%bond_order(105) = 0.037d0
+    test_target%bond_order(106) = 0.037d0
+    test_target%bond_order(107) = 0.037d0
+    test_target%bond_order(108) = 0.283d0
+    test_target%bond_order(109) = 0.037d0
+    test_target%bond_order(110) = 0.037d0
+    test_target%bond_order(111) = 0.037d0
+    test_target%bond_order(112) = 0.037d0
+    test_target%bond_order(113) = 0.037d0
+    test_target%bond_order(114) = 0.037d0
+    test_target%bond_order(115) = 0.037d0
+    test_target%bond_order(116) = 0.268d0
+    test_target%bond_order(117) = 0.037d0
+    test_target%bond_order(118) = 0.037d0
+    test_target%bond_order(119) = 0.037d0
+    test_target%bond_order(120) = 0.037d0
+    test_target%bond_order(121) = 0.037d0
+    test_target%bond_order(122) = 0.037d0
+    test_target%bond_order(123) = 0.283d0
+    test_target%bond_order(124) = 0.037d0
+    test_target%bond_order(125) = 0.037d0
+    test_target%bond_order(126) = 0.037d0
+    test_target%bond_order(127) = 0.037d0
+    test_target%bond_order(128) = 0.037d0
+    test_target%bond_order(129) = 0.037d0
+    test_target%bond_order(130) = 0.268d0
+    test_target%bond_order(131) = 0.037d0
+    test_target%bond_order(132) = 0.037d0
+    test_target%bond_order(133) = 0.037d0
+    test_target%bond_order(134) = 0.283d0
+    test_target%bond_order(135) = 0.037d0
+    test_target%bond_order(136) = 0.037d0
+    test_target%bond_order(137) = 0.037d0
+    test_target%bond_order(138) = 0.037d0
+    test_target%bond_order(139) = 0.037d0
+    test_target%bond_order(140) = 0.037d0
+    test_target%bond_order(141) = 0.037d0
+    test_target%bond_order(142) = 0.037d0
+    test_target%bond_order(143) = 0.037d0
+    test_target%bond_order(144) = 0.037d0
+    test_target%bond_order(145) = 0.268d0
+    test_target%bond_order(146) = 0.037d0
+    test_target%bond_order(147) = 0.037d0
+    test_target%bond_order(148) = 0.283d0
+    test_target%bond_order(149) = 0.037d0
+    test_target%bond_order(150) = 0.037d0
+    test_target%bond_order(151) = 0.037d0
+    test_target%bond_order(152) = 0.037d0
+    test_target%bond_order(153) = 0.037d0
+    test_target%bond_order(154) = 0.037d0
+    test_target%bond_order(155) = 0.037d0
+    test_target%bond_order(156) = 0.037d0
+    test_target%bond_order(157) = 0.037d0
+    test_target%bond_order(158) = 0.037d0
+    test_target%bond_order(159) = 0.268d0
+    test_target%bond_order(160) = 0.037d0
+    test_target%bond_order(161) = 0.037d0
+    test_target%bond_order(162) = 0.037d0
+    test_target%bond_order(163) = 0.283d0
+    test_target%bond_order(164) = 0.037d0
+    test_target%bond_order(165) = 0.037d0
+    test_target%bond_order(166) = 0.037d0
+    test_target%bond_order(167) = 0.037d0
+    test_target%bond_order(168) = 0.037d0
+    test_target%bond_order(169) = 0.037d0
+    test_target%bond_order(170) = 0.037d0
+    test_target%bond_order(171) = 0.037d0
+    test_target%bond_order(172) = 0.037d0
+    test_target%bond_order(173) = 0.037d0
+    test_target%bond_order(174) = 0.268d0
+    test_target%bond_order(175) = 0.037d0
+    test_target%bond_order(176) = 0.283d0
+    test_target%bond_order(177) = 0.037d0
+    test_target%bond_order(178) = 0.037d0
+    test_target%bond_order(179) = 0.037d0
+    test_target%bond_order(180) = 0.037d0
+    test_target%bond_order(181) = 0.037d0
+    test_target%bond_order(182) = 0.037d0
+    test_target%bond_order(183) = 0.037d0
+    test_target%bond_order(184) = 0.037d0
+    test_target%bond_order(185) = 0.037d0
+    test_target%bond_order(186) = 0.037d0
+    test_target%bond_order(187) = 0.037d0
+    test_target%bond_order(188) = 0.037d0
+    test_target%bond_order(189) = 0.268d0
+    test_target%bond_order(190) = 0.283d0
+    test_target%bond_order(191) = 0.037d0
+    test_target%bond_order(192) = 0.037d0
+    test_target%bond_order(193) = 0.037d0
+    test_target%bond_order(194) = 0.037d0
+    test_target%bond_order(195) = 0.037d0
+    test_target%bond_order(196) = 0.037d0
+    test_target%bond_order(197) = 0.037d0
+    test_target%bond_order(198) = 0.037d0
+    test_target%bond_order(199) = 0.037d0
+    test_target%bond_order(200) = 0.037d0
+    test_target%bond_order(201) = 0.037d0
+    test_target%bond_order(202) = 0.037d0
+    test_target%bond_order(203) = 0.268d0
+    test_target%bond_order(204) = 0.037d0
+    test_target%bond_order(205) = 0.283d0
+    test_target%bond_order(206) = 0.037d0
+    test_target%bond_order(207) = 0.037d0
+    test_target%bond_order(208) = 0.037d0
+    test_target%bond_order(209) = 0.037d0
+    test_target%bond_order(210) = 0.037d0
+    test_target%bond_order(211) = 0.037d0
+    test_target%bond_order(212) = 0.037d0
+    test_target%bond_order(213) = 0.037d0
+    test_target%bond_order(214) = 0.037d0
+    test_target%bond_order(215) = 0.037d0
+    test_target%bond_order(216) = 0.268d0
+    test_target%bond_order(217) = 0.037d0
+    test_target%bond_order(218) = 0.037d0
+    test_target%bond_order(219) = 0.037d0
+    test_target%bond_order(220) = 0.283d0
+    test_target%bond_order(221) = 0.037d0
+    test_target%bond_order(222) = 0.037d0
+    test_target%bond_order(223) = 0.037d0
+    test_target%bond_order(224) = 0.037d0
+    test_target%bond_order(225) = 0.037d0
+    test_target%bond_order(226) = 0.037d0
+    test_target%bond_order(227) = 0.037d0
+    test_target%bond_order(228) = 0.037d0
+    test_target%bond_order(229) = 0.037d0
+    test_target%bond_order(230) = 0.037d0
+    test_target%bond_order(231) = 0.268d0
+    test_target%bond_order(232) = 0.037d0
+    test_target%bond_order(233) = 0.037d0
+    test_target%bond_order(234) = 0.283d0
+    test_target%bond_order(235) = 0.037d0
+    test_target%bond_order(236) = 0.037d0
+    test_target%bond_order(237) = 0.037d0
+    test_target%bond_order(238) = 0.037d0
+    test_target%bond_order(239) = 0.037d0
+    test_target%bond_order(240) = 0.037d0
+    test_target%bond_order(241) = 0.037d0
+    test_target%bond_order(242) = 0.037d0
+    test_target%bond_order(243) = 0.037d0
+    test_target%bond_order(244) = 0.037d0
+    test_target%bond_order(245) = 0.268d0
+    test_target%bond_order(246) = 0.037d0
+    test_target%bond_order(247) = 0.037d0
+    test_target%bond_order(248) = 0.037d0
+    test_target%bond_order(249) = 0.283d0
+    test_target%bond_order(250) = 0.037d0
+    test_target%bond_order(251) = 0.037d0
+    test_target%bond_order(252) = 0.037d0
+    test_target%bond_order(253) = 0.037d0
+    test_target%bond_order(254) = 0.037d0
+    test_target%bond_order(255) = 0.037d0
+    test_target%bond_order(256) = 0.268d0
+    test_target%bond_order(257) = 0.037d0
+    test_target%bond_order(258) = 0.037d0
+    test_target%bond_order(259) = 0.037d0
+    test_target%bond_order(260) = 0.037d0
+    test_target%bond_order(261) = 0.037d0
+    test_target%bond_order(262) = 0.037d0
+    test_target%bond_order(263) = 0.283d0
+    test_target%bond_order(264) = 0.037d0
+    test_target%bond_order(265) = 0.037d0
+    test_target%bond_order(266) = 0.037d0
+    test_target%bond_order(267) = 0.037d0
+    test_target%bond_order(268) = 0.037d0
+    test_target%bond_order(269) = 0.037d0
+    test_target%bond_order(270) = 0.037d0
+    test_target%bond_order(271) = 0.268d0
+    test_target%bond_order(272) = 0.037d0
+    test_target%bond_order(273) = 0.037d0
+    test_target%bond_order(274) = 0.037d0
+    test_target%bond_order(275) = 0.037d0
+    test_target%bond_order(276) = 0.037d0
+    test_target%bond_order(277) = 0.283d0
+    test_target%bond_order(278) = 0.037d0
+    test_target%bond_order(279) = 0.037d0
+    test_target%bond_order(280) = 0.037d0
+    test_target%bond_order(281) = 0.037d0
+    test_target%bond_order(282) = 0.037d0
+    test_target%bond_order(283) = 0.037d0
+    test_target%bond_order(284) = 0.037d0
+    test_target%bond_order(285) = 0.268d0
+    test_target%bond_order(286) = 0.037d0
+    test_target%bond_order(287) = 0.037d0
+    test_target%bond_order(288) = 0.037d0
+    test_target%bond_order(289) = 0.037d0
+    test_target%bond_order(290) = 0.037d0
+    test_target%bond_order(291) = 0.037d0
+    test_target%bond_order(292) = 0.283d0
+    test_target%bond_order(293) = 0.037d0
+    test_target%bond_order(294) = 0.037d0
+    test_target%bond_order(295) = 0.037d0
+    test_target%bond_order(296) = 0.037d0
+    test_target%bond_order(297) = 0.037d0
+    test_target%bond_order(298) = 0.037d0
+    test_target%bond_order(299) = 0.037d0
+    test_target%bond_order(300) = 0.268d0
+    test_target%bond_order(301) = 0.037d0
+    test_target%bond_order(302) = 0.037d0
+    test_target%bond_order(303) = 0.037d0
+    test_target%bond_order(304) = 0.037d0
+    test_target%bond_order(305) = 0.283d0
+    test_target%bond_order(306) = 0.037d0
+    test_target%bond_order(307) = 0.037d0
+    test_target%bond_order(308) = 0.037d0
+    test_target%bond_order(309) = 0.037d0
+    test_target%bond_order(310) = 0.037d0
+    test_target%bond_order(311) = 0.037d0
+    test_target%bond_order(312) = 0.037d0
+    test_target%bond_order(313) = 0.037d0
+    test_target%bond_order(314) = 0.037d0
+    test_target%bond_order(315) = 0.268d0
+    test_target%bond_order(316) = 0.037d0
+    test_target%bond_order(317) = 0.037d0
+    test_target%bond_order(318) = 0.037d0
+    test_target%bond_order(319) = 0.283d0
+    test_target%bond_order(320) = 0.037d0
+    test_target%bond_order(321) = 0.037d0
+    test_target%bond_order(322) = 0.037d0
+    test_target%bond_order(323) = 0.037d0
+    test_target%bond_order(324) = 0.037d0
+    test_target%bond_order(325) = 0.037d0
+    test_target%bond_order(326) = 0.037d0
+    test_target%bond_order(327) = 0.037d0
+    test_target%bond_order(328) = 0.037d0
+    test_target%bond_order(329) = 0.268d0
+    test_target%bond_order(330) = 0.037d0
+    test_target%bond_order(331) = 0.037d0
+    test_target%bond_order(332) = 0.037d0
+    test_target%bond_order(333) = 0.037d0
+    test_target%bond_order(334) = 0.283d0
+    test_target%bond_order(335) = 0.037d0
+    test_target%bond_order(336) = 0.037d0
+    test_target%bond_order(337) = 0.037d0
+    test_target%bond_order(338) = 0.037d0
+    test_target%bond_order(339) = 0.037d0
+    test_target%bond_order(340) = 0.037d0
+    test_target%bond_order(341) = 0.037d0
+    test_target%bond_order(342) = 0.268d0
+    test_target%bond_order(343) = 0.037d0
+    test_target%bond_order(344) = 0.037d0
+    test_target%bond_order(345) = 0.037d0
+    test_target%bond_order(346) = 0.037d0
+    test_target%bond_order(347) = 0.037d0
+    test_target%bond_order(348) = 0.037d0
+    test_target%bond_order(349) = 0.283d0
+    test_target%bond_order(350) = 0.037d0
+    test_target%bond_order(351) = 0.037d0
+    test_target%bond_order(352) = 0.037d0
+    test_target%bond_order(353) = 0.037d0
+    test_target%bond_order(354) = 0.037d0
+    test_target%bond_order(355) = 0.037d0
+    test_target%bond_order(356) = 0.037d0
+    test_target%bond_order(357) = 0.268d0
+    test_target%bond_order(358) = 0.037d0
+    test_target%bond_order(359) = 0.037d0
+    test_target%bond_order(360) = 0.037d0
+    test_target%bond_order(361) = 0.037d0
+    test_target%bond_order(362) = 0.037d0
+    test_target%bond_order(363) = 0.283d0
+    test_target%bond_order(364) = 0.037d0
+    test_target%bond_order(365) = 0.037d0
+    test_target%bond_order(366) = 0.037d0
+    test_target%bond_order(367) = 0.037d0
+    test_target%bond_order(368) = 0.037d0
+    test_target%bond_order(369) = 0.037d0
+    test_target%bond_order(370) = 0.037d0
+    test_target%bond_order(371) = 0.268d0
+    test_target%bond_order(372) = 0.037d0
+    test_target%bond_order(373) = 0.037d0
+    test_target%bond_order(374) = 0.037d0
+    test_target%bond_order(375) = 0.037d0
+    test_target%bond_order(376) = 0.037d0
+    test_target%bond_order(377) = 0.037d0
+    test_target%bond_order(378) = 0.283d0
+    test_target%calc_vibe = .false.
+    test_target%nerror = 0
+
+    call mopac_scf(test_in, test_restore, test_out)
+    call test_output(test_name, test_in, test_target, test_out, nfail)
+end subroutine test_crystal1
+
 subroutine test_output(name, input, target, output, nfail)
     use mopac_api
     implicit none
-    character(20), intent(in) :: name
+    character(50), intent(in) :: name
     type(mopac_system), intent(in) :: input
     type(mopac_properties), intent(in) :: target, output
     integer, intent(inout) :: nfail
@@ -716,13 +2215,13 @@ subroutine test_output(name, input, target, output, nfail)
     ! compare heat
     if(abs(target%heat - output%heat) > heat_tol) then
         nfail = nfail + 1
-        write(*,*) "heat mismatch in test '", name, "':", target%heat, "vs", output%heat
+        write(*,*) "heat mismatch in test '", trim(name), "':", target%heat, "vs", output%heat
     end if
     ! compare dipole
     do i = 1, 3
         if(abs(target%dipole(i) - output%dipole(i)) > charge_tol) then
             nfail = nfail + 1
-            write(*,*) "dipole(", i, ") mismatch in test '", name, "':", target%dipole(i), &
+            write(*,*) "dipole(", i, ") mismatch in test '", trim(name), "':", target%dipole(i), &
             "vs", output%dipole(i)
         end if
     end do
@@ -730,7 +2229,7 @@ subroutine test_output(name, input, target, output, nfail)
     do i = 1, input%natom
         if(abs(target%charge(i) - output%charge(i)) > charge_tol) then
             nfail = nfail + 1
-            write(*,*) "charge(", i, ") mismatch in test '", name, "':", target%charge(i), &
+            write(*,*) "charge(", i, ") mismatch in test '", trim(name), "':", target%charge(i), &
             "vs", output%charge(i)
         end if
     end do
@@ -738,7 +2237,7 @@ subroutine test_output(name, input, target, output, nfail)
     do i = 1, 3*input%natom_move
         if(abs(target%coord_update(i) - output%coord_update(i)) > coord_tol) then
             nfail = nfail + 1
-            write(*,*) "coord_update(", i, ") mismatch in test '", name, "':", &
+            write(*,*) "coord_update(", i, ") mismatch in test '", trim(name), "':", &
             target%coord_update(i), "vs", output%coord_update(i)
         end if
     end do
@@ -746,20 +2245,20 @@ subroutine test_output(name, input, target, output, nfail)
     do i = 1, 3*input%natom_move
         if(abs(target%coord_deriv(i) - output%coord_deriv(i)) > deriv_tol) then
             nfail = nfail + 1
-            write(*,*) "coord_deriv(", i, ") mismatch in test '", name, "':", &
+            write(*,*) "coord_deriv(", i, ") mismatch in test '", trim(name), "':", &
             target%coord_deriv(i), "vs", output%coord_deriv(i)
         end if
     end do
     ! compare vibrational properties
     if(target%calc_vibe .neqv. output%calc_vibe) then
         nfail = nfail + 1
-        write(*,*) "calc_vibe mistmatch in test '", name, "':", target%calc_vibe, &
+        write(*,*) "calc_vibe mistmatch in test '", trim(name), "':", target%calc_vibe, &
         "vs", output%calc_vibe
     else if(target%calc_vibe .eqv. .true.) then
         do i = 1, 3*input%natom_move
             if(abs(target%freq(i) - output%freq(i)) > freq_tol) then
                 nfail = nfail + 1
-                write(*,*) "freq(", i, ") mismatch in test '", name, "':", target%freq(i), &
+                write(*,*) "freq(", i, ") mismatch in test '", trim(name), "':", target%freq(i), &
                 "vs", output%freq(i)
             end if
         end do
@@ -778,7 +2277,7 @@ subroutine test_output(name, input, target, output, nfail)
             do j=1, 3*input%natom_move
                 if(abs(fmat1(i,j) - fmat2(i,j)) > freq_tol*freq_tol) then
                     nfail = nfail + 1
-                    write(*,*) "fmat(", i, ",", j, ") mismatch in test '", name, "':", fmat1(i,j), &
+                    write(*,*) "fmat(", i, ",", j, ") mismatch in test '", trim(name), "':", fmat1(i,j), &
                     "vs", fmat2(i,j)
                 end if
             end do
@@ -789,24 +2288,24 @@ subroutine test_output(name, input, target, output, nfail)
     do i = 1, input%natom+1
         if(target%bond_index(i) /= output%bond_index(i)) then
             nfail = nfail + 1
-            write(*,*) "bond_index(", i, ") mismatch in test '", name, "':", target%bond_index(i), "vs", output%bond_index(i)
+            write(*,*) "bond_index(", i, ") mismatch in test '", trim(name), "':", target%bond_index(i), "vs", output%bond_index(i)
         end if
     end do
     do i = 1, target%bond_index(input%natom+1)-1
         if(target%bond_atom(i) /= output%bond_atom(i)) then
             nfail = nfail + 1
-            write(*,*) "bond_atom(", i, ") mismatch in test '", name, "':", target%bond_atom(i), "vs", output%bond_atom(i)
+            write(*,*) "bond_atom(", i, ") mismatch in test '", trim(name), "':", target%bond_atom(i), "vs", output%bond_atom(i)
         end if
         if(abs(target%bond_order(i) - output%bond_order(i)) > charge_tol) then
             nfail = nfail + 1
-            write(*,*) "bond_order(", i, ") mismatch in test '", name, "':", target%bond_order(i), "vs", output%bond_order(i)
+            write(*,*) "bond_order(", i, ") mismatch in test '", trim(name), "':", target%bond_order(i), "vs", output%bond_order(i)
         end if
     end do
     ! compare updated coordinates
     do i = 1, 3*input%nlattice_move
         if(abs(target%lattice_update(i) - output%lattice_update(i)) > coord_tol) then
             nfail = nfail + 1
-            write(*,*) "lattice_update(", i, ") mismatch in test '", name, "':", &
+            write(*,*) "lattice_update(", i, ") mismatch in test '", trim(name), "':", &
             target%lattice_update(i), "vs", output%lattice_update(i)
         end if
     end do
@@ -814,7 +2313,7 @@ subroutine test_output(name, input, target, output, nfail)
     do i = 1, 3*input%nlattice_move
         if(abs(target%lattice_deriv(i) - output%lattice_deriv(i)) > deriv_tol) then
             nfail = nfail + 1
-            write(*,*) "lattice_deriv(", i, ") mismatch in test '", name, "':", &
+            write(*,*) "lattice_deriv(", i, ") mismatch in test '", trim(name), "':", &
             target%lattice_deriv(i), "vs", output%lattice_deriv(i)
         end if
     end do
@@ -822,13 +2321,13 @@ subroutine test_output(name, input, target, output, nfail)
     do i = 1, 6
         if(abs(target%stress(i) - output%stress(i)) > stress_tol) then
             nfail = nfail + 1
-            write(*,*) "stress(", i, ") mismatch in test '", name, "':", target%stress(i), &
+            write(*,*) "stress(", i, ") mismatch in test '", trim(name), "':", target%stress(i), &
             "vs", output%stress(i)
         end if
     end do
-    ! compare status
-    if(target%status /= output%status) then
+    ! compare error status
+    if(target%nerror /= output%nerror) then
         nfail = nfail + 1
-        write(*,*) "status mismatch in test '", name, "':", target%status, "vs", output%status
+        write(*,*) "nerror mismatch in test '", trim(name), "':", target%nerror, "vs", output%nerror
     end if
 end subroutine test_output