Merge pull request #222 from openmopac/api-bugfix

API bug fix

README.md

@@ -95,7 +95,8 @@ There is a [new documentation website](https://openmopac.github.io) under develo
 ## Interfaces
 
 While MOPAC is primarily a self-contained command-line program whose behavior is specified by an input file, it also has other modes of
-operation, some of which only require the MOPAC shared library and not the executable.
+operation, some of which only require the MOPAC shared library and not the executable. Note that API calls to the MOPAC library are not
+thread safe. Each thread must load its own instance of the MOPAC library, such as by running independent calling programs.
 
 ### MDI Engine
 

src/MOZYME/iter_for_MOZYME.F90

@@ -15,7 +15,7 @@
 ! along with this program.  If not, see <https://www.gnu.org/licenses/>.
 
 subroutine iter_for_MOZYME (ee)
-    use molkst_C, only: norbs, step_num, numcal, nscf, escf, &
+    use molkst_C, only: norbs, step_num, numcal, numcal0, nscf, escf, &
        & numat,  enuclr, atheat, emin, keywrd, moperr, line, use_disk
 !
     use chanel_C, only: iw, iend, end_fn
@@ -519,7 +519,7 @@ subroutine iter_for_MOZYME (ee)
           goto 80
         end if
       end if
-      if (bigscf .or. numcal /= 1) then
+      if (bigscf .or. numcal /= 1+numcal0) then
           call diagg (f, nocc1, nvir1,  idiagg,  partp, indi)
         idiagg = idiagg + 1
       else
@@ -610,7 +610,7 @@ subroutine iter_for_MOZYME (ee)
         backspace (iw)
       end if
       call isitsc (escf, selcon, emin, iemin, iemax, okscf, niter, itrmax)
-      if ( .not. bigscf .and. numcal == 1) then
+      if ( .not. bigscf .and. numcal == 1+numcal0) then
         exit
       else if (okscf .and. niter > 1 .and. (emin /= 0.d0 .or. niter > 3)) then
         exit

src/MOZYME/tidy.F90

@@ -16,7 +16,7 @@
 
 subroutine tidy (nmos_loc, nc, ic, n01, c, n02, nnc_loc, ncmo, ln, mn, mode)
     use MOZYME_C, only: iorbs, jopt, thresh, numred
-    use molkst_C, only: numat, step_num, norbs, moperr, keywrd, numcal, use_disk
+    use molkst_C, only: numat, step_num, step_num0, norbs, moperr, keywrd, numcal, use_disk
     use chanel_C, only: iw
     implicit none
     integer, intent (in) :: mode,  n02, nmos_loc
@@ -115,7 +115,7 @@ subroutine tidy (nmos_loc, nc, ic, n01, c, n02, nnc_loc, ncmo, ln, mn, mode)
       !
       !   Do the LMOs of the SCF need to be put at the start of the storage?
       !
-      if (isnew /= 0 .or. step_num <= 1 .or. step_num == imode(mode)) exit
+      if (isnew /= 0 .or. step_num <= 1+step_num0 .or. step_num == imode(mode)) exit
       isnew = 2
       if (numred >= numat-1) then
         isnew = 1

src/input/getdat.F90

@@ -65,16 +65,6 @@ subroutine getdat(input, output)
             call getarg (run, jobnam)
 #endif
             natoms = 1
-            do i = len_trim(jobnam), 1, -1   !  Remove any unprintable characters from the end of the file-name
-              if (ichar(jobnam(i:i)) > 39 .and. ichar(jobnam(i:i)) < 126 .or. jobnam(i:i) =="'") exit
-            end do
-            jobnam(i + 1:) = " "
-!
-!  Replace backslash with forward-slash
-!
-            do i = 1, len_trim(jobnam)
-              if (jobnam(i:i) == backslash) jobnam(i:i) = "/"
-            end do
           else if (i == 0) then
             if (is_PARAM) then
               write(line,'(2a)')" PARAM is the parameter optimization program for use with MOPAC"
@@ -99,6 +89,17 @@ subroutine getdat(input, output)
           natoms = 1
         end if
       end if
+!  Remove any unprintable characters from the end of the file-name
+      do i = len_trim(jobnam), 1, -1
+        if (ichar(jobnam(i:i)) > 39 .and. ichar(jobnam(i:i)) < 126 .or. jobnam(i:i) =="'") exit
+      end do
+      jobnam(i + 1:) = " "
+!
+!  Replace backslash with forward-slash
+!
+      do i = 1, len_trim(jobnam)
+        if (jobnam(i:i) == backslash) jobnam(i:i) = "/"
+      end do
       if (natoms == 0) return
 !
 ! Check for the data set in the order: <file>.mop, <file>.dat, <file>
@@ -294,7 +295,7 @@ subroutine getdat(input, output)
       i = index(keywrd, "GEO-DAT")
       if (i /= 0) keywrd(i:i+6) = "GEO_DAT"
       i = index(keywrd, "GEO-REF")
-      if (i /= 0) keywrd(i:i+6) = "GEO_REF"
+      if (i /= 0) keywrd(i:i+6) = "GEO_REF"
       double_plus = (index(keywrd, " ++ ") /= 0)
       if (index(keywrd, " GEO_DAT") + index(keywrd, " SETUP")/= 0) then
         nlines = nlines + 3

src/input/getgeo.F90

@@ -22,7 +22,7 @@ subroutine getgeo(iread, labels, geo, xyz, lopt, na, nb, nc, int)
       use parameters_C, only : ams
 !
       use molkst_C, only : natoms, keywrd, numat, maxtxt, line, moperr, &
-        numcal, id, units, Angstroms, arc_hof_1, arc_hof_2, keywrd_txt, pdb_label
+        numcal, numcal0, id, units, Angstroms, arc_hof_1, arc_hof_2, keywrd_txt, pdb_label
 !
       use chanel_C, only : iw, ir, input_fn, end_fn, iend
 !
@@ -162,11 +162,11 @@ subroutine getgeo(iread, labels, geo, xyz, lopt, na, nb, nc, int)
       if (line == '$end') go to 20
       if (line(1:1) == '*') go to 20
       if (line == ' ') then
-        if(natoms == 0 .and. numcal == 1) then
+        if(natoms == 0 .and. numcal == 1+numcal0) then
 !
 !  Check:  Is this an ARC file?
 !
-          numcal = 2
+          numcal = 2+numcal0
           rewind (iread)
           sum = 0.d0
           do i = 1, 10000
@@ -572,7 +572,7 @@ subroutine getgeo(iread, labels, geo, xyz, lopt, na, nb, nc, int)
           write(iw,'(/10x,a,i5)')"Faulty atom:", natoms
           write(iw,'(/10x,a)')"Faulty line: """//trim(line)//""""
           call mopend("Unless MINI is used, optimization flags must be 1, 0, or -1")
-          numcal = 2
+          numcal = 2+numcal0
           if ((lopt(1,natoms) > 10 .or. lopt(2,natoms) > 10 .or. lopt(3,natoms) > 10) .and. natoms > 1) &
             write(iw,'(/10x,a)')" If the geometry is in Gaussian format, add keyword ""AIGIN"" and re-run"
           return
@@ -626,7 +626,7 @@ subroutine getgeo(iread, labels, geo, xyz, lopt, na, nb, nc, int)
 !***********************************************************************
   120 continue
       if (natoms == 0) then
-        if (numcal == 1) call mopend (' Error detected while reading geometry')
+        if (numcal == 1+numcal0) call mopend (' Error detected while reading geometry')
         return
       end if
       if ( .not. Angstroms) then

src/input/gettxt.F90

@@ -16,7 +16,7 @@
 
       subroutine gettxt
       use chanel_C, only: ir, iw, isetup, input_fn
-      use molkst_C, only: keywrd, keywrd_quoted, koment, title, refkey, gui, numcal, line, &
+      use molkst_C, only: keywrd, keywrd_quoted, koment, title, refkey, gui, numcal, numcal0, line, &
         moperr, allkey, backslash
       implicit none
 !-----------------------------------------------
@@ -107,7 +107,7 @@ subroutine gettxt
         if (.not. exists) then
           if (setup_present .and. .not. zero_scf) then
             write (line, '(A)') "SETUP FILE """//trim(filen)//""" MISSING."
-            numcal = 2
+            numcal = 2+numcal0
             if (.not. gui )write(0,'(//30x,a)')' SETUP FILE "'//trim(filen)//'" MISSING'
             call mopend (trim(line))
             return
@@ -401,7 +401,7 @@ subroutine gettxt
       go to 60
 50    continue
       if (zero_scf) go to 60
-      numcal = 2
+      numcal = 2+numcal0
       call mopend ('SETUP FILE "'//trim(filen)//'" MISSING')
       write(iw,'(a)') " (Setup file name: '"//trim(filen)//"')"
       return
@@ -411,7 +411,7 @@ subroutine gettxt
 100   continue
       call split_keywords(oldkey)
 
-      if (numcal > 1) then
+      if (numcal > 1+numcal0) then
         if (index(keywrd,"OLDGEO") /= 0) return ! User forgot to add extra lines for title and comment
         if (aux) keywrd = " AUX"
         line = "JOB ENDED NORMALLY"

src/input/readmo.F90

@@ -27,7 +27,7 @@ subroutine readmo
 !
       USE symmetry_C, ONLY: idepfn, locdep, depmul, locpar
 !
-      use molkst_C, only : ndep, numat, numcal, natoms, nvar, keywrd, dh, &
+      use molkst_C, only : ndep, numat, numcal, numcal0, natoms, nvar, keywrd, dh, &
       & verson, is_PARAM, line, nl_atoms, l_feather, backslash, &
       & moperr, maxatoms, koment, title, method_pm6, refkey, l_feather_1, &
       isok, method_pm6_dh2, caltyp, keywrd_quoted, &
@@ -464,13 +464,13 @@ subroutine readmo
             intern = .false.
           else
             call getgeo (ir, labels, geo, coord, lopt, na, nb, nc, intern)
-            if (numcal == 1 .and. natoms == 0) then
+            if (numcal == 1+numcal0 .and. natoms == 0) then
               i = index(keywrd, "GEO_DAT")
               if (i /= 0) then
                 write(line,'(2a)')" GEO_DAT file """//trim(line_1)//""" exists, but does not contain any atoms."
                 write(0,'(//10x,a,//)')trim(line)
                 call mopend(trim(line))
-              else if (.not. gui .and. numcal < 2) then
+              else if (.not. gui .and. numcal < 2+numcal0) then
                 write(line,'(2a)')" Data set '"//trim(job_fn)//" exists, but does not contain any atoms."
                 write(0,'(//10x,a,//)')trim(line)
                 call mopend(trim(line))
@@ -498,7 +498,7 @@ subroutine readmo
               coorda(:,:numat) = geo(:,:numat)
               numat_old = numat
             else if (natoms /= -3) then
-              if (moperr .and. numcal == 1) return
+              if (moperr .and. numcal == 1+numcal0) return
               if (maxtxt > txtmax) txtmax = maxtxt
               txtatm1(:natoms) = txtatm(:natoms)
               if (index(keywrd, " RESID") /= 0) txtatm1(:numat)(22:22) = " "
@@ -623,7 +623,7 @@ subroutine readmo
             end do
           end if
           if (natoms < 0 ) then
-            if (numcal == 1) rewind ir
+            if (numcal == 1+numcal0) rewind ir
             if (.not.isok) then
               write (iw, '(A)') &
                 ' Use AIGIN to allow more geometries to be used'
@@ -634,7 +634,7 @@ subroutine readmo
               stop
             end if
             isok = .FALSE.
-            if (numcal > 2) then
+            if (numcal > 2+numcal0) then
               naigin = naigin + 1
               write (iw, '(2/,2A)') '   GAUSSIAN INPUT REQUIRES', &
                 ' STAND-ALONE JOB'
@@ -647,15 +647,15 @@ subroutine readmo
             go to 10
           end if
         end if
-        if (natoms == 0 .and. numcal == 1) then
+        if (natoms == 0 .and. numcal == 1+numcal0) then
           call mopend ('NO ATOMS IN SYSTEM')
           return
         end if
       else
 !
 !   Use the old geometry, if one exists
 !
-        if (numcal == 1) then
+        if (numcal == 1+numcal0) then
           write(line,'(a)')" Keyword OLDGEO cannot be used in the first calculation - there is no old geometry"
           write(iw,'(//10x,a)')trim(line)
           call to_screen(trim(line))
@@ -689,7 +689,7 @@ subroutine readmo
           call mopend(trim(line))
           return
         else
-          if (numcal == 1 .and. numat > 50) write(0,'(10x,a)')idate//"  Job: '"//trim(jobnam)//"' started successfully"
+          if (numcal == 1+numcal0 .and. numat > 50) write(0,'(10x,a)')idate//"  Job: '"//trim(jobnam)//"' started successfully"
         end if
       end if
       maxci = 10000

src/molkst_C.F90

@@ -110,6 +110,7 @@ module molkst_C
                  !  Each stage is limited to the same electronic structure.  Most calculations
                  !  will only have one stage, e.g. geometry optimization or force constants.
                  !
+  &  job_no0, numcal0, step_num0, & ! Needed for repeated API calls to run_mopac
   &  mpack,    & !  Number of elements in a lower-half-triangle = (norbs*(norbs+1))/2
   &  n2elec,   & !  Number of two-electron integrals
   &  nscf,     & !  Number of SCF calculations done

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, dummy, errtxt
+        use molkst_C, only : line, job_no, job_no0, natoms, dummy, errtxt
         implicit none
         integer, intent (in) :: ntxt
         character, intent (in) :: txt*(*)
@@ -83,7 +83,7 @@ subroutine summary(txt, ntxt)
            end if
         end if
         if (ntxt == 1) then
-          if (natoms == 0 .and. job_no == 1) then
+          if (natoms == 0 .and. job_no == job_no0) then
             write(iw,'(/10x, a)')"Job failed to run because no atoms were detected in the system"
             write(iw,'(10x, a, /)')"The start of the data-set is as follows:"
             rewind (ir)

src/output/to_screen.F90

@@ -63,7 +63,7 @@ subroutine current_version (text)
   method_am1, method_mndo, method_pm3, method_rm1, method_mndod, method_pm6, &
   method_pm7, nvar, koment, keywrd, zpe, id, density, natoms, formula, press, voigt, &
   uhf, nalpha, nbeta,  gnorm, mozyme, mol_weight, ilim, &
-  line, nscf, time0, sz, ss2, no_pKa, title, jobnam, job_no, fract
+  line, nscf, time0, sz, ss2, no_pKa, title, jobnam, job_no, job_no0, fract
 !
   use MOZYME_C, only : ncf, ncocc, noccupied, icocc_dim, cocc_dim, nvirtual, icvir_dim, &
   nncf, iorbs, cocc, icocc, ncvir, nnce, nce, icvir, cvir, tyres, size_mres, &
@@ -1340,7 +1340,7 @@ subroutine current_version (text)
 !
 !  Don't print processor-independent CPU times for quick jobs - that would waste too much time.
 !
-    if (time0 > 1.d0 .and. job_no < 4) then
+    if (time0 > 1.d0 .and. job_no < 4+job_no0) then
 !
 !  Deliberately run a time-consuming calculation to work out CPU speed
 !  The "j" index is set to use up 1.0 seconds on the development computer.

src/output/writmo.F90

@@ -18,7 +18,7 @@ subroutine  writmo
       use cosmo_C, only : iseps, area, fepsi, cosvol, ediel, solv_energy
 !
       use molkst_C, only : numat, nclose, nopen, fract, nalpha, nelecs, nbeta, &
-      & norbs, nvar, gnorm, iflepo, enuclr,elect, ndep, nscf, numcal, escf, &
+      & norbs, nvar, gnorm, iflepo, enuclr,elect, ndep, nscf, numcal, numcal0, escf, &
       & keywrd, os, verson, time0, moperr, last, iscf, id, pressure, mol_weight, &
       jobnam, line, mers, uhf, method_indo, &
       density, formula, mozyme, mpack, stress, &
@@ -971,7 +971,7 @@ subroutine  writmo
           call pdbout(31)
           close (31)
         end if
-        if (numcal == 2) then
+        if (numcal == 2+numcal0) then
           if (index(keywrd, "OLDGEO") /= 0) then
 !
 ! Write a warning that OLDGEO has been used, so user is aware that multiple ARC files are present