Merge pull request #199 from openmopac/xyz-bug-fix

XYZ & packaging bug fixes

CMakeLists.txt

@@ -286,7 +286,7 @@ if(CMAKE_Fortran_COMPILER_ID STREQUAL "Intel")
     if(INTEL_OMP_DLL STREQUAL "INTEL_OMP_DLL-NOTFOUND")
       message("WARNING: Cannot locate Intel's OpenMP dll for packaging")
     endif()
-    install(PROGRAMS ${INTEL_OMP_DLL} DESTINATION ${CMAKE_INSTALL_BINDIR} COMPONENT qtifw EXCLUDE_FROM_ALL)
+    install(PROGRAMS ${INTEL_OMP_DLL} DESTINATION ${CMAKE_INSTALL_BINDIR} COMPONENT redist EXCLUDE_FROM_ALL)
   else()
     find_library(INTEL_OMP_LIB iomp5 HINT CMAKE_Fortran_COMPILER PATH_SUFFIXES ../../compiler/lib ../../compiler/lib/intel64_lin)
     if(INTEL_OMP_LIB STREQUAL "INTEL_OMP_LIB-NOTFOUND")

src/Common_arrays_C.F90

@@ -192,9 +192,6 @@ module Common_arrays_C
                    !  Units         dimensionless
                    !  Min inclusive -1.0
                    !  Max inclusive  1.0
-  &  workmat1,   & !  unpacked matrix workspace #1
-  &  workmat2,   & !  unpacked matrix workspace #2
-  &  workmat3,   & !  unpacked matrix workspace #3
   &  dxyz2,      & !  dxyz, but 2-D
   &  ptot2,      & !
      xmat,       &

src/SCF/iter.F90

@@ -19,7 +19,7 @@ subroutine iter(ee, fulscf, rand)
        &  c, nat, nfirst, nlast, eigb, pdiag, f, w, wk, fb
       use iter_C, only : pold, pold2, pbold, pbold2, &
       & pold3, pbold3, vec_ai, vec_bi, fock_ai, fock_bi, p_ai, &
-      p_bi, h_ai, h_bi, vecl_ai, vecl_bi
+      p_bi, h_ai, h_bi, vecl_ai, vecl_bi, pulay_work1, pulay_work2, pulay_work3
       USE parameters_C, only :
       USE funcon_C, only : fpc_9
       USE maps_C, ONLY: latom
@@ -812,6 +812,14 @@ subroutine iter(ee, fulscf, rand)
 !                                                                      *
 !***********************************************************************
            if (okpuly .and. makea .and. iredy>1) then
+            if (.not. Allocated (pulay_work1)) then
+              allocate (pulay_work1(norbs, norbs), pulay_work2(norbs, norbs), &
+              & pulay_work3(norbs, norbs), stat=i)
+              if (i /= 0) then
+                call memory_error("Pulay converger in Iter")
+                return
+              end if
+            end if
 #ifdef GPU
               if (lgpu) then
                  call pulay_for_gpu (f, pa, norbs, pold, pold2, pold3, &
@@ -925,6 +933,14 @@ subroutine iter(ee, fulscf, rand)
 !                                                                      *
 !***********************************************************************
             if (okpuly .and. makeb .and. iredy>1) then
+              if (.not. Allocated (pulay_work1)) then
+                allocate (pulay_work1(norbs, norbs), pulay_work2(norbs, norbs), &
+                & pulay_work3(norbs, norbs), stat=i)
+                if (i /= 0) then
+                  call memory_error("Pulay converger in Iter")
+                  return
+                end if
+              end if
 #ifdef GPU
               if (lgpu) then
                  call pulay_for_gpu (fb, pb, norbs, pbold, pbold2, pbold3, &
@@ -1095,7 +1111,7 @@ subroutine delete_iter_arrays
 !  become corrupt if INTERP is called.
 !
       use iter_C, only : vec_ai, vec_bi, fock_ai, fock_bi, p_ai, &
-      p_bi, h_ai, h_bi, vecl_ai, vecl_bi
+      p_bi, h_ai, h_bi, vecl_ai, vecl_bi, pulay_work1, pulay_work2, pulay_work3
       implicit none
         if (allocated(vec_ai)) deallocate(vec_ai)
         if (allocated(vec_bi)) deallocate(vec_bi)
@@ -1107,6 +1123,9 @@ subroutine delete_iter_arrays
         if (allocated(h_bi)) deallocate(h_bi)
         if (allocated(vecl_ai)) deallocate(vecl_ai)
         if (allocated(vecl_bi)) deallocate(vecl_bi)
+        if (allocated(pulay_work1)) deallocate(pulay_work1)
+        if (allocated(pulay_work2)) deallocate(pulay_work2)
+        if (allocated(pulay_work3)) deallocate(pulay_work3)
       end subroutine delete_iter_arrays
 
 

src/SCF/iter_C.F90

@@ -22,6 +22,6 @@ module iter_C
 !
       double precision, dimension(:), allocatable :: h_ai, vecl_ai, h_bi, vecl_bi
       double precision, dimension(:,:), allocatable :: vec_ai, fock_ai, &
-       & p_ai, vec_bi, fock_bi, p_bi
+       & p_ai, vec_bi, fock_bi, p_bi, pulay_work1, pulay_work2, pulay_work3
 
       end module iter_C

src/SCF/pulay.F90

@@ -17,7 +17,7 @@
       subroutine pulay(f, p, n, fppf, fock, emat, lfock, nfock, msize, start, pl)
       use chanel_C, only : iw
       use molkst_C, only : numcal, keywrd, mpack, npulay
-      use Common_arrays_C, only : workmat1, workmat2, workmat3
+      use iter_C, only : pulay_work1, pulay_work2, pulay_work3
       implicit none
       integer  :: n
       integer , intent(inout) :: lfock
@@ -95,10 +95,10 @@ subroutine pulay(f, p, n, fppf, fock, emat, lfock, nfock, msize, start, pl)
 !
 !      call mamult (p, f, fppf(lbase+1), n, 0.D0)
 !      call mamult (f, p, fppf(lbase+1), n, -1.D0)
-      call unpack_matrix(p, workmat1, n)
-      call unpack_matrix(f, workmat2, n)
-      call sym_commute(workmat1, workmat2, workmat3, n)
-      call pack_matrix(workmat3, fppf(lbase+1), n)
+      call unpack_matrix(p, pulay_work1, n)
+      call unpack_matrix(f, pulay_work2, n)
+      call sym_commute(pulay_work1, pulay_work2, pulay_work3, n)
+      call pack_matrix(pulay_work3, fppf(lbase+1), n)
 !
 !   FPPF NOW CONTAINS THE RESULT OF FP - PF.
 !

src/input/getgeo.F90

@@ -206,18 +206,18 @@ subroutine getgeo(iread, labels, geo, xyz, lopt, na, nb, nc, int)
         end if
         if (ii == 0) goto 120
         ii = 0
-      end if
-!
-!   If two quotation marks are side-by-side, force a space in between them
-!
-      do
-        i = index(line, '""')
-        if (i /= 0) then
-          string = trim(line)
-          line = string(:i)//" "//trim(string(i + 1:))
-        else
-          exit
-        end if
+      end if
+!
+!   If two quotation marks are side-by-side, force a space in between them
+!
+      do
+        i = index(line, '""')
+        if (i /= 0) then
+          string = trim(line)
+          line = string(:i)//" "//trim(string(i + 1:))
+        else
+          exit
+        end if
       end do
       ltl = len_trim(line)
       icomma = ichar(comma)
@@ -793,6 +793,7 @@ subroutine getgeo(iread, labels, geo, xyz, lopt, na, nb, nc, int)
             k = k + lopt(j,i)
           end do
         end do
+        j = na(3)
 !
 !  Get rid of dummy atoms
 !
@@ -811,7 +812,7 @@ subroutine getgeo(iread, labels, geo, xyz, lopt, na, nb, nc, int)
 !   If everything is marked for optimization then unconditionally mark the first
 !   three atoms for optimization
 !
-        if (k >= 3*numat - 6) lopt(:,:min(3, numat)) = 1
+        if (k >= 3*numat - 6 .and. j /= 0) lopt(:,:min(3, numat)) = 1
         natoms = numat
         if (saddle .or. (index(keywrd, " LOCATE-TS") /= 0)) then
           if (index(keywrd, " LOCATE-TS") == 0) lopt(:,:numat) = 1  ! In a saddle calculation, all parameters must be optimizable.

src/input/readmo.F90

@@ -338,6 +338,13 @@ subroutine readmo
           lopt(3,3) = 0
         end if
         if (index(keywrd, " XYZ") /= 0) then
+          k = 0
+          do i = 1, natoms
+            do j = 1, 3
+              k = k + lopt(j,i)
+            end do
+          end do
+          j = na(3)
           numat = 0
           do i = 1, natoms
             if (labels(i) /= 99) then
@@ -353,7 +360,7 @@ subroutine readmo
 !   If everything is marked for optimization then unconditionally mark the first
 !   three atoms for optimization
 !
-          if (k >= 3*numat - 6) lopt(:,:min(3, numat)) = 1
+          if (k >= 3*numat - 6 .and. j /= 0) lopt(:,:min(3, numat)) = 1
           natoms = numat
         end if
 
@@ -453,7 +460,7 @@ subroutine readmo
           if (index(line, " PDB ") + index(line, " PDB(")/= 0) then
             call getpdb(geo)
             coorda(:,:numat) = geo(:,:numat)
-            numat_old = numat
+            numat_old = numat
             intern = .false.
           else
             call getgeo (ir, labels, geo, coord, lopt, na, nb, nc, intern)
@@ -771,7 +778,7 @@ subroutine readmo
           dh = "H   "
         end if
         keywrd = trim(keywrd_txt)
-      end if
+      end if
       l_feather = (method_PM7 .or. method_pm6_org .or. index(keywrd, " MOZ") /= 0 .and. (index(keywrd, " PM6") /= 0))
       l_feather_1 = (index(keywrd, " MACRO") /= 0)
       write (iw, &

src/setup_mopac_arrays.F90

@@ -28,7 +28,7 @@ subroutine setup_mopac_arrays(n, mode)
      & h, w, p, pa, pb, f, c, eigs, dxyz, grad, eigb, txtatm, fb, cb, &
      & wk, errfn, aicorr, nbonds, ibonds, na_store, q, nw, lopt, &
      & hesinv, gnext1, gmin1, ifact, i1fact, ptot2, geoa, l_atom, coorda, &
-       txtatm1, workmat1, workmat2, workmat3
+       txtatm1
 !
   USE maps_C, only : react
 !
@@ -143,17 +143,12 @@ subroutine setup_mopac_arrays(n, mode)
           if (allocated(eigs))     deallocate(eigs)
           if (allocated(q))        deallocate(q)
           if (allocated(eigb))     deallocate(eigb)
-          if (allocated(workmat1)) deallocate(workmat1)
-          if (allocated(workmat2)) deallocate(workmat2)
-          if (allocated(workmat3)) deallocate(workmat3)
           allocate(h(mpack), p(mpack), pa(mpack), pb(mpack), stat=i)
           j = j + i
           allocate(pold(npulay*mpack),  pold2(npulay*mpack), f(mpack), stat=i)
           j = j + i
           allocate(c(norbs, norbs), eigs(norbs + 1), q(numat), stat=i)
           j = j + i
-          allocate(workmat1(norbs, norbs), workmat2(norbs, norbs), workmat3(norbs, norbs), stat=i)
-          j = j + i
           allocate(eigb(norbs + 1), pold3(max(mpack, 400)), stat=i)
           j = j + i
           allocate(errfn(3*natoms*l123), aicorr(nvar), stat=i)
@@ -468,9 +463,6 @@ subroutine setup_mopac_arrays(n, mode)
     if (allocated(cc0))        deallocate(cc0, stat = i)
     if (allocated(wk0))        deallocate(wk0, stat = i)
     if (allocated(q0))         deallocate(q0, stat = i)
-    if (allocated(workmat1))   deallocate(workmat1, stat = i)
-    if (allocated(workmat2))   deallocate(workmat2, stat = i)
-    if (allocated(workmat3))   deallocate(workmat3, stat = i)
     end if
   end subroutine setup_mopac_arrays
   subroutine memory_error(txt)

tests/compare_output.py

@@ -31,12 +31,12 @@
 EIGVEC_THRESHOLD = 1e-3
 
 # regular expression pattern for a time stamp or other signifier of timing output, "CLOCK" or "TIME" or "SECONDS", & system-dependent versioning
-skip_criteria = re.compile('([A-Z][a-z][a-z] [A-Z][a-z][a-z] [ 0-9][0-9] [0-9][0-9]:[0-9][0-9]:[0-9][0-9] [0-9][0-9][0-9][0-9])'
-                           '|(CLOCK)|(TIME)|(SECONDS)|(Version)|(THE VIBRATIONAL FREQUENCY)|(ITERATION)|(SCF CALCULATIONS)|(Stewart)'
-                           '|(remaining)|(\*  THREADS)|(\*  ISOTOPE)|(\*  DENOUT)|(\*  OLDENS)|(\*  SETUP)|(ITER.)|(\*\* )|(web-site)|(MOPAC)|(GRADIENT NORM)')
+skip_criteria = re.compile(r'([A-Z][a-z][a-z] [A-Z][a-z][a-z] [ 0-9][0-9] [0-9][0-9]:[0-9][0-9]:[0-9][0-9] [0-9][0-9][0-9][0-9])'
+                           r'|(CLOCK)|(TIME)|(SECONDS)|(Version)|(THE VIBRATIONAL FREQUENCY)|(ITERATION)|(SCF CALCULATIONS)|(Stewart)'
+                           r'|(remaining)|(\*  THREADS)|(\*  ISOTOPE)|(\*  DENOUT)|(\*  OLDENS)|(\*  SETUP)|(ITER.)|(\*\* )|(web-site)|(MOPAC)|(GRADIENT NORM)')
 
 # regular expression pattern for an eigenvector block
-eigen_criteria = re.compile('(Root No.)|(ROOT NO.)')
+eigen_criteria = re.compile(r'(Root No.)|(ROOT NO.)')
 
 def is_float(string):
     '''check if a string contains a float'''
@@ -213,6 +213,14 @@ def compare_mopac_out_file(out_line, out_list, ref_line, ref_list, heat_error_th
     '''Compares the output to the given reference'''
 
     if len(ref_list) != len(out_list):
+        dump = open("ref.dump", "w")
+        for item in ref_list:
+            dump.write(f"{item}\n")
+        dump.close()
+        dump = open("out.dump", "w")
+        for item in out_list:
+            dump.write(f"{item}\n")
+        dump.close()
         assert len(ref_list) == len(out_list), f'ERROR: output file size mismatch, {len(ref_list)} vs. {len(out_list)}'
         #print(f'WARNING: output file size mismatch, {len(ref_list)} vs. {len(out_list)}')