Merge pull request #211 from openmopac/fix-forcets

FORCETS bug fixes

src/SCF/fock2.F90

@@ -30,8 +30,8 @@ subroutine fock2(f, ptot, p, w, wj, wk, numat, nfirst, nlast, mode)
 !   D u m m y   A r g u m e n t s
 !-----------------------------------------------
       integer  :: numat, mode
-      integer  :: nfirst(numat_ref)
-      integer  :: nlast(numat_ref)
+      integer  :: nfirst(abs(numat))
+      integer  :: nlast(abs(numat))
       double precision  :: f(mpack)
       double precision , intent(in) :: ptot(mpack)
       double precision  :: p(mpack)
@@ -87,6 +87,8 @@ subroutine fock2(f, ptot, p, w, wj, wk, numat, nfirst, nlast, mode)
 !
 !   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
           ifact(i) = (i*(i - 1))/2
           i1fact(i) = ifact(i) + i

src/forces/force.F90

@@ -136,7 +136,7 @@ subroutine force()
       end if
       if (id > 0) then
         natoms = numat + id
-        nvar = nvar - 3*id ! No need to check for TS - solids can't have transition states.
+        if(.not. ts) nvar = nvar - 3*id
         labels(numat + 1:numat + id) = 107
       end if
       na = 0
@@ -329,7 +329,7 @@ subroutine force()
       end do
       if (debug) then
         write (iw, '(2/10X,'' FULL FORCE MATRIX, INVOKED BY "DFORCE"'')')
-        if (index(keywrd, " NOREOR") == 0) then
+        if (index(keywrd, " NOREOR") == 0 .and. .not. ts) then
           write(iw,'(/10x,a)')" Caution: NOREOR is NOT present, therefore system will be oriented"
           write(iw,'(10x,a)')" so that the moments of inertia are along the Cartesian axes."
         end if

src/output/matou1.F90

@@ -125,6 +125,7 @@ subroutine matou1(a, b, ncx, nr, ndim, iflag)
       if (numat == 0) go to 50
       if (nlast(numat) /= nr) go to 50
       do i = 1, numat
+        if (iflag == 5 .and. 3*i > nr) exit
         jlo = nfirst(i)
         jhi = nlast(i)
         l = nat(i)

src/properties/mullik.F90

@@ -251,7 +251,7 @@ subroutine mullik()
       nlower = (norbs*(norbs + 1))/2
       call density_for_GPU (vecs, fract, nclose, nopen, 2.d0, nlower, norbs, 2, pb, 3)
 !
-      pb = pb*store
+      pb(:ifact(norbs+1)) = pb(:ifact(norbs+1))*store(:ifact(norbs+1))
       summ = 0.D0
       do i = 1, norbs
         sum = 0