Moved from C++ to C. Updated Python files.

.gitignore

@@ -1,3 +1,8 @@
+.vscode
+bin/
+build/
+include/
+lib/
 
 ### C ###
 # Prerequisites
@@ -19,7 +24,6 @@
 *.gch
 *.pch
 
-
 # Libraries
 *.dll
 *.so
@@ -89,6 +93,7 @@ ehthumbs_vista.db
 
 # Folder config file
 Desktop.ini
+.DS_Store
 
 # Recycle Bin used on file shares
 $RECYCLE.BIN/

src/basis/getlibr.py

@@ -1,26 +1,36 @@
-#!/usr/bin/env python3
+#!/usr/bin/python3
 # This script downloads the basis set library data from www.basissetexchange.org
 # into the directory $NWCHEM_TOP/src/basis/libraries.bse
-# to use, set the env. variable NWCHEM_BASIS_LIBRARY=$NWCHEM_TOP/src/basis/libraries.bse/
+# To run, cd $NWCHEM_TOP/src/basis/libraries.bse/ && ../getlibr.py
+# this will update the content of $NWCHEM_TOP/src/basis/libraries.bse
+# To use the updates library, set the env. variable
+NWCHEM_BASIS_LIBRARY=$NWCHEM_TOP/src/basis/libraries.bse/
 # Requires the installation of the python env. from
 # https://github.com/MolSSI-BSE/basis_set_exchange
+# e.g. python3 -m pip install --user basis_set_exchange
 # See https://molssi-bse.github.io/basis_set_exchange/
 #
 # names changed
 # def2-universal-jfit was weigend_coulomb_fitting
 # dgauss-a1-dftjfit   was dgauss_a1_dft_coulomb_fitting
 # dgauss-a2-dftjfit   was dgauss_a2_dft_coulomb_fitting
-#
+
 import basis_set_exchange as bse
 from datetime import datetime
 today = datetime.now().isoformat(timespec='minutes')
 print(today)
 all_bs = bse.get_all_basis_names()
 md = bse.get_metadata()
-for bas_name in all_bs:
-    #get version and list of elements
-    version_bs = md[bas_name]['latest_version']
-    elements_list = md[bas_name]['versions'][version_bs]['elements']
+summary_file = open('summary.txt','w')
+
+def writebs(md, bas_name, summary_file, get_aux=0):
+    md_bas_name = bas_name.lower()
+    md_bas_name = md_bas_name.replace("*","_st_")
+    md_bas_name = md_bas_name.replace("/","_sl_")
+    print(' md_bas_name '+md_bas_name+"\n")
+    print(' bas_name '+bas_name+"\n")
+    version_bs = md[md_bas_name]['latest_version']
+    elements_list = md[md_bas_name]['versions'][version_bs]['elements']
     #open file
     # get rid of asterisks
     file_name = bas_name.replace("*","s")
@@ -33,19 +43,42 @@ for bas_name in all_bs:
     file_name = file_name.replace(" ","_")
     #replace forward slash with underscore
     file_name = file_name.replace("/","_")
+    #lowercase
+    file_name = file_name.lower()
+    if get_aux==1:
+        file_name = file_name + "-autoaux"
     print(' file name is '+file_name+"\n")
     output_file = open(file_name,'w')
     output_file.write('# BSE Version '+bse.version()+'\n')
-    output_file.write('# Data downloaded at '+today+'\n')
-    output_file.write('# '+bas_name+' version number '+version_bs+'\n')
-    output_file.write('# Description: '+md[bas_name]['description']+'\n')
-    output_file.write('# Role: '+md[bas_name]['role']+'\n')
-    output_file.write('# '+bse.get_references(bas_name,fmt='txt').replace('\n','\n# '))
-    output_file.write('# \n')
+    output_file.write('# Data downloaded on '+today+'\n')
+
+    if get_aux==0:
+        output_file.write('# '+bas_name+' version number '+version_bs+'\n')
+        output_file.write('# Description: '+md[md_bas_name]['description']+'\n')
+        output_file.write('# Role: '+md[md_bas_name]['role']+'\n')
+        output_file.write('# '+bse.get_references(bas_name,fmt='txt').replace('\n','\n# '))
+        output_file.write('# \n')
+    elif get_aux==1:
+        output_file.write('# '+bas_name+' version number '+version_bs+' AutoAux \n')
+        output_file.write('# Role: JK Fitting \n')
+        output_file.write('# Stoychev GL, Auer AA, Neese F. \n# Automatic Generation of Auxiliary
+Basis Sets.\n# J Chem Theory Comput. 2017 Feb 14;13(2):554-562.\n# doi:
+10.1021/acs.jctc.6b01041.\n')
+        output_file.write('# \n')
+
+    n_elements=0
+    for element in elements_list:
+        n_elements = n_elements + 1
+    if get_aux==1:
+        summary_file.write('Basis set \"'+bas_name+'-autoaux\" (number of atoms
+'+str(n_elements)+')\n')
+    else:
+        summary_file.write('Basis set \"'+bas_name+'\" (number of atoms '+str(n_elements)+')\n')
     for element in elements_list:
         #element='h'
         try:
-            bs_str=bse.get_basis(bas_name, header=False, elements=element, fmt='nwchem', optimize_general=True, uncontract_general=True)
+            bs_str=bse.get_basis(bas_name, header=False, elements=element, fmt='nwchem',
+optimize_general=True, uncontract_general=True, get_aux=get_aux)
         except:
 #            print("failed for"+element)
             pass
@@ -54,11 +87,22 @@ for bas_name in all_bs:
             bs_str=bs_str.replace("END","end")
             bs_str=bs_str.replace("PRINT","")
             element_str=bse.misc.compact_elements([element])
-            bs_str=bs_str.replace("ao basis",element_str+"_"+bas_name)
+            if get_aux==1:
+                bs_str=bs_str.replace("ao basis",element_str+"_"+bas_name+"-autoaux")
+            else:
+                bs_str=bs_str.replace("ao basis",element_str+"_"+bas_name)
             #ECP
             bs_str=bs_str.replace("ECP","ecp \""+element_str+"_"+bas_name+"\"")
             output_file.write(bs_str)
             #
             print(bas_name+" "+element_str)
-print("end")
+    return
 
+for bas_name in all_bs:
+    md_bas_name = bas_name.lower()
+    md_bas_name = md_bas_name.replace("*","_st_")
+    md_bas_name = md_bas_name.replace("/","_sl_")
+    writebs(md, bas_name, summary_file)
+    if md[md_bas_name]['role'] == 'orbital':
+        writebs(md, bas_name, summary_file, get_aux=1)
+print("end")

src/ddscf/comp4_ext.c

@@ -0,0 +1,32 @@
+#include <stdio.h>
+#include "bitops_decls.h"
+#include "bitops_funcs.h"
+
+void comp4_extract(int* m, int i, double s, int nb_per_i) {
+
+    int v;          // Value after compression
+
+#if defined(CRAY)
+  int vv, vvv;
+#endif
+
+    int index, nbits;
+    double fast[] = {0.0, 1.0e-13, 1.0e-12, 1.0e-11, 1.0e-10, 1.0e-9,
+               1.0e-8, 1.0e-7, 1.0e-6, 1.0e-5, 1.0e-4, 1.0e-3, 1.0e-2,
+               1.0e-1, 1.0e0, 1.0e1};
+
+    v = 15;
+    index = (i - 1)/(2*nb_per_i) + 1;
+    nbits = 4*(i - (index-1)*(2*nb_per_i) - 1);
+#if defined(CRAY)
+    vvv =  shiftl(v, nbits);
+    vv = shiftr(iand(m(index), vvv), nbits);
+    v = iand(vv,15);
+#else
+    v = iand(ishft(iand(m(index), ishft(v, nbits)), -nbits),15);
+#endif
+
+    s = fast(v);
+    // printf('%d -> %d %d %d %0.4f');
+
+}

src/ddscf/fock_2e_file.c

@@ -0,0 +1,142 @@
+#include "util.h"
+#include "cscfps.h"
+#include "cfock.h"
+#include <bool.h>
+#include <math.h>
+
+void fock_2e_from_file(int geom, int basis, int nfock, int ablklen,
+                 double jfac[nfock], double kfac[nfock], double tol2e, bool oskel,
+                 double dij[nfock*ablklen], double dik[nfock*ablklen], double dli[nfock*ablklen],
+                 double djk[nfock*ablklen], double dlj[nfock*ablklen], double dlk[nfock*ablklen],
+                 double fij[nfock*ablklen], double fik[nfock*ablklen], double fli[nfock*ablklen],
+                 double fjk[nfock*ablklen], double flj[nfock*ablklen], double flk[nfock*ablklen],
+                 double tmp, int vg_dens[nfock], int vg_fock[nfock]) {
+
+    //$Id$
+
+    /*Accumulate the contribution to the fock matrices from
+      integrals store in the integral file.  Simply read thru
+      the file getting a range of indices, fetch the corresponding
+      density matrix blocks and then read the integrals in that
+      block.
+    */
+
+    double den_tol, denmax, dtol2e;
+    int ilo, jlo, klo, llo;
+    int ihi, jhi, khi, lhi;
+    int ijk_prev[3][2];
+    int blklen;
+
+    bool int2e_get_bf_range, int2e_file_read;
+
+    if (oscfps) pstat_on(ps_fock_io);
+
+    den_tol = fmax(tol2e*0.01, 1e-300); // To avoid a hard zero
+
+    ijk_prev[0][0] = -1;
+    ijk_prev[1][0] = -1;
+    ijk_prev[2][0] = -1;
+    ijk_prev[0][1] = -1;
+    ijk_prev[1][1] = -1;
+    ijk_prev[2][1] = -1;
+
+    blklen = nfock*ablklen;
+    dfill(blklen, 0.0e0, fij, 1);
+    dfill(blklen, 0.0e0, fik, 1);
+    dfill(blklen, 0.0e0, fli, 1);
+    dfill(blklen, 0.0e0, fjk, 1);
+    dfill(blklen, 0.0e0, flj, 1);
+    dfill(blklen, 0.0e0, flk, 1);
+
+    // Loop over blocks of integral labels
+
+    while (int2e_get_bf_range(ilo,ihi,jlo,jhi,klo,khi,llo,lhi)) {
+        // Get matrices for this block of labels
+        fock_init_cmul(ihi-ilo+1,jhi-jlo+1,lhi-llo+1);
+        fock_2e_cache_dens_fock(
+              ilo, jlo, klo, llo,
+              ihi, jhi, khi, lhi,
+              ijk_prev,
+              nfock, vg_dens, vg_fock,
+              jfac, kfac,
+              dij, dik, dli, djk, dlj, dlk,
+              fij, fik, fli, fjk, flj, flk,
+              tmp);
+
+        fock_density_screen(nfock,
+              ilo, jlo, klo, llo,
+              ihi, jhi, khi, lhi,
+              ilo, jlo, klo, llo,
+              ihi, jhi, khi, lhi,
+              dij, dik, dli, djk, dlj, dlk, denmax)
+
+        dtol2e = min(dentolmax, den_tol/max(1e-10,denmax), den_tol/max(1e-10,denmax**2))
+
+        call int2e_file_fock_block(nfock, dtol2e,
+                 dij, dik, dli, djk, dlj, dlk,
+                 fij, fik, fli, fjk, flj, flk)
+
+        // Update F blocks
+
+        call fock_upd_blk(nfock, vg_fock,
+              llo, lhi, ilo, ihi, kfac, fli, tmp)
+        call fock_upd_blk(nfock, vg_fock,
+              llo, lhi, jlo, jhi, kfac, flj, tmp)
+        call fock_upd_blk(nfock, vg_fock,
+              llo, lhi, klo, khi, jfac, flk, tmp)
+    }
+
+    if (ijk_prev[0][0]) != -1) {
+        fock_upd_blk(nfock, vg_fock,
+              ijk_prev[0][0]), ijk_prev[0][1]),
+              ijk_prev(2,1), ijk_prev(2,2),
+              jfac, fij, tmp)
+        fock_upd_blk(nfock, vg_fock,
+              ijk_prev(2,1), ijk_prev(2,2),
+              ijk_prev(3,1), ijk_prev(3,2),
+              kfac, fjk, tmp )
+        fock_upd_blk( nfock, vg_fock,
+              ijk_prev(1,1), ijk_prev(1,2),
+              ijk_prev(3,1), ijk_prev(3,2),
+              kfac, fik, tmp )
+     }
+
+    if (oscfps) pstat_off(ps_fock_io);
+
+}
+
+void fock_2e_rep_from_file(int geom, int basis, int nfock, int nbf,
+           double jfac[nfock], double kfac[nfock], double tol2e, bool oskel,
+           double dens[nfock][nbf*nbf], fock[nfock][nbf*nbf]) {
+
+    double den_tol, denmax;
+    int ilo, jlo, klo, llo;
+    int ihi, jhi, khi, lhi, i, j;
+
+    bool int2e_get_bf_range, int2e_file_read;
+    int idamax;
+
+    if (oscfps) pstat_on(ps_fock_io);
+
+    denmax = 0.0;
+    for (i = 0; i < nfock; i++) {
+        j = idamax(nbf*nbf, dens[i][0], nfock);
+        denmax = max(denmax, abs(dens[i][j]);
+    }
+    // return if DM is null (e.g imaginary part of RTTDFT DM at t=0)
+    if (denmax < 1e-12) return;
+    den_tol = min(dentolmax,tol2e/denmax,tol2e/denmax**2) // Threshold to screen integs only
+
+    if (ga_nodeid() == 0 && util_print('fockfile',print_debug)) {
+        printf("fockfile: tols %d %d %d %d", tol2e, dentolmax, denmax, den_tol);
+    }
+
+    fock_init_cmul(nbf,nbf,nbf) // lookup table for f build
+
+    // Loop over blocks of integral labels
+}
+
+
+
+
+}

src/nwchem.c

@@ -265,7 +265,7 @@ int main(){
             errquit('start: rtdb_open old failed', 0, RTDB_ERR);
         }
     }
-    
+
 
 //    initialize nxtask
     nxtask_init(rtdb);
@@ -344,7 +344,7 @@ int main(){
             errquit('control: rtdb_print failed', 0, RTDB_ERR);
         }
     }
-    
+
     if (!rtdb_close(rtdb, 'keep')){
         errquit('nwchem: rtdb_close failed', rtdb, RTDB_ERR);
     }
@@ -372,7 +372,7 @@ int main(){
             ga_print_stats();
             write(LuOut,*);
         }
-        
+
     }
 
 

src/rtdb/rtdb.h

@@ -1,105 +1,183 @@
-//
-//    Header file for intial FORTRAN interface to RTDB
-//    (see the C header file rtdb.h for more detail)
-//
-//    All functions return .TRUE. on success, .FALSE. on failure
-//
-//    All functions are also mirrored by routines rtdb_* -> rtdb_par_*
-//    in which process 0 performs the operation and all other processes
-//    are broadcast the result of a read and discard writes.
-//
-//    rtdb_max_key  ... an integer parameter that defines the maximum
-//                      length of a character string key
-//
-//    rtdb_max_file ... an integer parameter that defines the maximum
-//                      length of a file name
-//
-//
-//    logical function rtdb_parallel(mode)
-//    logical mode              [input]
-//
-//
-//    logical function rtdb_open(filename, mode, handle)
-//    character *(*) filename   [input]
-//    character *(*) mode       [input]
-//    integer handle            [output]
-//
-//    logical function rtdb_clone(handle, suffix)
-//    integer handle            [input]
-//    character*(*) suffix        [input]
-//
-//    logical function rtdb_close(handle, mode)
-//    integer handle            [input]
-//    character*(*) mode        [input]
-//
-//    logical function rtdb_put(handle, name, ma_type, nelem, array)
-//    integer handle            [input]
-//    character *(*) name       [input]
-//    integer ma_type           [input]
-//    integer nelem             [input]
-//    <ma_type>array(nelem)     [input]
-//
-//    logical function rtdb_get_info(handle, name, ma_type, nelem, date)
-//    integer handle            [input]
-//    character *(*) name       [input]
-//    integer ma_type           [output]
-//    integer nelem             [output]
-//    character*26 date         [output]
-//
-//    logical function rtdb_get(handle, name, ma_type, nelem, array)
-//    integer handle            [input]
-//    character *(*) name       [input]
-//    integer ma_type           [input]
-//    integer nelem             [input]
-//    <ma_type>array(nelem)     [output]
-//
-//    logical function rtdb_ma_get(handle, name, ma_type, nelem, ma_handle)
-//    integer handle            [input]
-//    character *(*) name       [input]
-//    integer ma_type           [output]
-//    integer nelem             [output]
-//    integer ma_handle         [output]
-//
-//    logical function rtdb_cput(handle, name, nelem, buf)
-//    integer handle            [input]
-//    character *(*) name       [input]
-//    character *(*) buf        [input]
-//
-//    logical function rtdb_cget(handle, name, nelem, buf)
-//    integer handle            [input]
-//    character *(*) name       [input]
-//    character *(*) buf        [output]
-//
-//    logical function rtdb_print(handle, print_values)
-//    integer handle            [input]
-//    logical print_values      [input]
-//
-//    logical function rtdb_first(handle, name)
-//    integer handle            [input]
-//    character *(*) name       [output]
-//
-//    logical function rtdb_next(handle, name)
-//    integer handle            [input]
-//    character *(*) name       [output]
-//
-//    logical function rtdb_delete(handle, name)
-//    integer handle            [input]
-//    character *(*) name       [input]
-//
-bool rtdb_open, rtdb_close, rtdb_put, rtdb_get, rtdb_ma_get,
-     rtdb_cput, rtdb_cget, rtdb_print, rtdb_get_info,
-     rtdb_first, rtdb_next, rtdb_delete, rtdb_parallel,
-     rtdb_clone,rtdb_getfname,rtdb_cget_size;
-//$Id$
-extern rtdb_open, rtdb_close, rtdb_put, rtdb_get, rtdb_ma_get,
-     rtdb_cput, rtdb_cget, rtdb_print, rtdb_get_info,
-     rtdb_first, rtdb_next, rtdb_delete, rtdb_parallel,
-     rtdb_clone,rtdb_getfname,rtdb_cget_size;
-//
-//    Check these values against rtdb_f2c.c
-//
-const int rtdb_max_key=255;
-const int rtdb_max_file=255;
-//
-const bool rtdb_seq_mode = false;
-const bool rtdb_par_mode = true;
+#ifndef _RTDB_H
+#define _RTDB_H
+
+/*
+  All routines return TRUE (1) on success, FALSE (0) on failure.
+
+  int rtdb_parallel(const int mode)
+
+    Set the parallel access mode of all databases to mode and
+    return the previous setting
+
+
+  int rtdb_open(const char *filename, const char *mode, int *handle)
+
+    Filename = path to file associated with the data base
+    mode     = 'new'     Open only if it does not exist already
+               'old',    Open only if it does exist already
+               'unknown' Create new or open existing (preserving contents)
+               'empty'   Create new or open existing (deleting contents)
+               'scratch' Create new or open existing (deleting contents)
+                         and automatically delete upon closing.  Also, items
+                         cached in memory are not written to disk.
+
+    handle   = returns handle by which all future references to the
+               data base are made
+
+
+
+  int rtdb_clone(const int handle, const char *suffix)
+
+    Copy the data base file
+
+    handle   = handle to RTDB
+    suffix
+
+
+  int rtdb_close(const int handle, const char *mode)
+
+    Close the data base
+
+    handle   = handle to RTDB
+    mode     = 'keep'    Preserve the data base file to enable restart
+               'delete'  Delete the data base file freeing all resources
+
+               mode is overridden by opening the data base with
+               mode='scratch' in which instance it is always deleted
+               upon closing
+
+
+  int rtdb_get_info(const int handle, const char *name, int *ma_type,
+                    int *nelem, char date[26])
+
+    Get info about an entry from the data base
+
+    handle   = handle to RTDB
+    name     = entry name (null terminated character string)
+    ma_type  = returns MA type of the entry
+    nelem    = returns no. of elements of the given type
+    date     = returns date of insertion (null terminated character string)
+
+
+  int rtdb_put(const int handle, const char *name, const int ma_type,
+               const int nelem, const void *array)
+
+    Insert an entry into the data base replacing previous entry
+
+    handle   = handle to RTDB
+    name     = entry name (null terminated character string)
+    ma_type  = MA type of the entry
+    nelem    = no. of elements of the given type
+    array    = data to be inserted
+
+
+  int rtdb_get(const int handle, const char *name, const int ma_type,
+               const int nelem, void *array)
+
+    Get an entry from the data base
+
+    handle   = handle to RTDB
+    name     = entry name (null terminated character string)
+    ma_type  = MA type of the entry which must match entry type
+    nelem    = size of array in units of ma_type
+    array    = user provided buffer that returns data
+
+
+  int rtdb_ma_get(const int handle, const char *name, int *ma_type,
+                  int *nelem, int *ma_handle)
+
+    Get an entry from the data base returning an MA handle
+
+    handle   = handle to RTDB
+    name     = entry name (null terminated character string)
+    ma_type  = returns MA type of the entry
+    nelem    = returns no. of elements of type ma_type in data
+    ma_handle= returns MA handle to data
+
+
+  int rtdb_first(const int handle, const int namelen, char *name)
+
+    Return the name of the first (user inserted) entry in the data base.
+    The order is effectively random.
+
+    handle  = handle to RTDB
+    namelen = size of user provided buffer name
+    name    = name of entry is returned in this buffer
+
+
+  int rtdb_next(const int handle, const int namelen, char *name)
+
+    Return the name of the next (user inserted) entry in the data base.
+    The order is effectively random.
+
+    handle  = handle to RTDB
+    namelen = size of user provided buffer name
+    name    = name of entry is returned in this buffer
+
+
+  int rtdb_print(const int handle, const int print_values)
+
+    Print the contents of the data base to stdout
+
+    handle  = handle to RTDB
+    print_values = boolean flag ... if true values as well as
+                   keys are printed out.
+
+
+  int rtdb_delete(const int handle, const char *name)
+
+    Delete the entry from the database.
+    Return
+          1 if key was present and successfully deleted
+
+         0 if key was not present, or if an error occured
+
+    handle  = handle to RTDB
+    name    = name of entry to delete
+
+*/
+
+int rtdb_open(const char *, const char *, int *);
+int rtdb_clone(const int, const char *);
+int rtdb_getfname(const int,  char [36]);
+int rtdb_close(const int, const char *);
+int rtdb_put(const int, const char *, const int, const int,
+		    const void *);
+int rtdb_get(const int, const char *, const int, const int,
+		    bool);
+int rtdb_get_info(const int, const char *, int *, int *, char [26]);
+int rtdb_ma_get(const int, const char *, int *, int *, int *);
+int rtdb_first(const int, const int, char *);
+int rtdb_next(const int, const int, char *);
+int rtdb_print(const int, const int);
+int rtdb_delete(const int, const char *);
+int rtdb_parallel(const int);
+
+/*
+  Following are 'sequential' versions of the above
+  for internal use only
+*/
+
+int rtdb_seq_open(const char *, const char *, int *);
+int rtdb_seq_copy(const int, const char *);
+int rtdb_seq_getfname(const int,  char [36]);
+int rtdb_seq_close(const int, const char *);
+int rtdb_seq_put(const int, const char *, const int, const int,
+		    const void *);
+int rtdb_seq_get(const int, const char *, const int, const int,
+		    void *);
+int rtdb_seq_get_info(const int, const char *, int *, int *, char [26]);
+int rtdb_seq_ma_get(const int, const char *, int *, int *, int *);
+int rtdb_seq_first(const int, const int, char *);
+int rtdb_seq_next(const int, const int, char *);
+int rtdb_seq_print(const int, const int);
+int rtdb_seq_delete(const int, const char *);
+
+#define RTDB_SEQ_MODE 0 //* Sequential mode
+#define RTDB_PAR_MODE 1 //* Parallel mode
+
+#if (defined(CRAY) || defined(WIN32)) &&!defined(__crayx1) &&!defined(__MINGW32__)
+#include "rtdb.cray.h"
+#endif
+
+#endif

src/tce/oce.py

@@ -26,14 +26,19 @@ def readfromfile(filename):
 def stringtooperatorsequence(expression):
    """Converts a string to an operatorsequence object"""
    # Syntax of the string is rather loosely defined as:
-   # (1) Numerical factor (with no permutation allowed) (optional), summation (optional), amplitudes (optional), normal ordered sequence,
-   # (2) Numerical factor can be an arithmatic expression such as (1.0/4.0),
+   # (1) Numerical factor (with no permutation allowed) (optional), summation (optional),
+amplitudes (optional), normal ordered sequence,
+   # (2) Numerical factor can be an arithmetic expression such as (1.0/4.0),
    # (3) Summation starts with either "SUM" or "sum" followed by a parenthesis of indexes,
-   # (4) Indexes can be either in one-letter notation (a-h, A-H for virtuals, i-o, I-O for occupieds, p-z, P-Z for either, case matters) 
-   #     or in OCE notation (p1,p2 for virtuals, h3,h4 for occupieds, g5,g6 for either, no overlap in numbering)
-   # (5) Amplitudes start with "t" or any name followed by a dagger ("+") indicating complex conjugate (optional) and a parenthesis of indexes,
+   # (4) Indexes can be either in one-letter notation (a-h, A-H for virtuals, i-o, I-O for
+occupieds, p-z, P-Z for either, case matters)
+   #     or in OCE notation (p1,p2 for virtuals, h3,h4 for occupieds, g5,g6 for either, no
+overlap in numbering)
+   # (5) Amplitudes start with "t" or any name followed by a dagger ("+") indicating complex
+conjugate (optional) and a parenthesis of indexes,
    # (6) Normal ordered operator sequence must exist even when it is empty "{}".
-   # (7) An example is: (1.0/16.0) Sum (p q r s c d k l) v(p q r s) t(c d k l) {i+ j+ b a}{p+ q+ s r}{c+ d+ l k}
+   # (7) An example is: (1.0/16.0) Sum (p q r s c d k l) v(p q r s) t(c d k l) {i+ j+ b a}{p+ q+
+s r}{c+ d+ l k}
 
    sequences = expression[expression.index("{"):]
    expression = expression[0:expression.index("{")]
@@ -41,9 +46,11 @@ def stringtooperatorsequence(expression):
    # first we decipher normal ordered operator sequence and define operators with/without daggers
    newsequences = []
    while (string.find(sequences,"{") != -1):
-      sequences = sequences[0:string.find(sequences,"{")] + sequences[string.find(sequences,"{")+1:]
+      sequences = sequences[0:string.find(sequences,"{")] +
+sequences[string.find(sequences,"{")+1:]
    if (sequences[len(sequences)-1] != "}"):
-      raise RuntimeError("Syntax error: the string must end with a normal ordered operator sequence")
+      raise RuntimeError("Syntax error: the string must end with a normal ordered operator
+sequence")
    sequences = sequences[0:len(sequences)-1]
    sequences = string.split(sequences,"}")
    for sequence in sequences:
@@ -58,7 +65,8 @@ def stringtooperatorsequence(expression):
          if (len(index) == 1):
             # one letter notation (vir: a-h & A-H; occ: i-o & I-O; gen: p-z & P-Z)
             if (((index >= "a") and (index <= "h")) or ((index >= "A") and (index <= "H"))):
-               newsequence.append(Operator("particle",dagger,string.ascii_letters.index(index)+1))
+
+newsequence.append(Operator("particle",dagger,string.ascii_letters.index(index)+1))
             elif (((index >= "i") and (index <= "o")) or ((index >= "I") and (index <= "O"))):
                newsequence.append(Operator("hole",dagger,string.ascii_letters.index(index)+1))
             elif (((index >= "p") and (index <= "z")) or ((index >= "P") and (index <= "Z"))):
@@ -71,8 +79,7 @@ def stringtooperatorsequence(expression):
             elif (index[0] == "g"):
                newsequence.append(Operator("general",dagger,int(index[1:])))
             else:
-               print("Syntax error: an operator not recognized")
-               stop
+               raise SyntaxError(" an operator not recognized")
       operatorlist = operatorlist + newsequence
       newsequences.append(newsequence)
 
@@ -126,23 +133,27 @@ def stringtooperatorsequence(expression):
          # one letter notation (vir: a-h & A-H; occ: i-o & I-O; gen: p-z & P-Z)
          if (((index >= "a") and (index <= "h")) or ((index >= "A") and (index <= "H"))):
             for indexinthelist in operatorlist:
-               if ((indexinthelist.type == "particle") and (indexinthelist.index == string.ascii_letters.index(index)+1)):
+               if ((indexinthelist.type == "particle") and (indexinthelist.index ==
+string.ascii_letters.index(index)+1)):
                   summation.indexes.append(indexinthelist)
                   break
          elif (((index >= "i") and (index <= "o")) or ((index >= "I") and (index <= "O"))):
             for indexinthelist in operatorlist:
-               if ((indexinthelist.type == "hole") and (indexinthelist.index == string.ascii_letters.index(index)+1)):
+               if ((indexinthelist.type == "hole") and (indexinthelist.index ==
+string.ascii_letters.index(index)+1)):
                   summation.indexes.append(indexinthelist)
                   break
          elif (((index >= "p") and (index <= "z")) or ((index >= "P") and (index <= "Z"))):
             for indexinthelist in operatorlist:
-               if ((indexinthelist.type == "general") and (indexinthelist.index == string.ascii_letters.index(index)+1)):
+               if ((indexinthelist.type == "general") and (indexinthelist.index ==
+string.ascii_letters.index(index)+1)):
                   summation.indexes.append(indexinthelist)
                   break
       else:
          if (index[0] == "p"):
             for indexinthelist in operatorlist:
-               if ((indexinthelist.type == "particle") and (indexinthelist.index == int(index[1:]))):
+               if ((indexinthelist.type == "particle") and (indexinthelist.index ==
+int(index[1:]))):
                   summation.indexes.append(indexinthelist)
                   break
          elif (index[0] == "h"):
@@ -152,12 +163,12 @@ def stringtooperatorsequence(expression):
                   break
          elif (index[0] == "g"):
             for indexinthelist in operatorlist:
-               if ((indexinthelist.type == "general") and (indexinthelist.index == int(index[1:]))):
+               if ((indexinthelist.type == "general") and (indexinthelist.index ==
+int(index[1:]))):
                   summation.indexes.append(indexinthelist)
                   break
          else:
-            print("syntax error")
-            stop
+            raise SyntaxError(" ")
 
    # get amplitudes
    remainder = string.split(remainder,")")
@@ -191,38 +202,43 @@ def stringtooperatorsequence(expression):
             # one letter notation (vir: a-h & A-H; occ: i-o & I-O; gen: p-z & P-Z)
             if (((index >= "a") and (index <= "h")) or ((index >= "A") and (index <= "H"))):
                for indexinthelist in operatorlist:
-                  if ((indexinthelist.type == "particle") and (indexinthelist.index == string.ascii_letters.index(index)+1)):
+                  if ((indexinthelist.type == "particle") and (indexinthelist.index ==
+string.ascii_letters.index(index)+1)):
                      newamplitude.indexes.append(indexinthelist)
                      break
             elif (((index >= "i") and (index <= "o")) or ((index >= "I") and (index <= "O"))):
                for indexinthelist in operatorlist:
-                  if ((indexinthelist.type == "hole") and (indexinthelist.index == string.ascii_letters.index(index)+1)):
+                  if ((indexinthelist.type == "hole") and (indexinthelist.index ==
+string.ascii_letters.index(index)+1)):
                      newamplitude.indexes.append(indexinthelist)
                      break
             elif (((index >= "p") and (index <= "z")) or ((index >= "P") and (index <= "Z"))):
                for indexinthelist in operatorlist:
-                  if ((indexinthelist.type == "general") and (indexinthelist.index == string.ascii_letters.index(index)+1)):
+                  if ((indexinthelist.type == "general") and (indexinthelist.index ==
+string.ascii_letters.index(index)+1)):
                      newamplitude.indexes.append(indexinthelist)
                      break
          else:
             if (index[0] == "p"):
                for indexinthelist in operatorlist:
-                  if ((indexinthelist.type == "particle") and (indexinthelist.index == int(index[1:]))):
+                  if ((indexinthelist.type == "particle") and (indexinthelist.index ==
+int(index[1:]))):
                      newamplitude.indexes.append(indexinthelist)
                      break
             elif (index[0] == "h"):
                for indexinthelist in operatorlist:
-                  if ((indexinthelist.type == "hole") and (indexinthelist.index == int(index[1:]))):
+                  if ((indexinthelist.type == "hole") and (indexinthelist.index ==
+int(index[1:]))):
                      newamplitude.indexes.append(indexinthelist)
                      break
             elif (index[0] == "g"):
                for indexinthelist in operatorlist:
-                  if ((indexinthelist.type == "general") and (indexinthelist.index == int(index[1:]))):
+                  if ((indexinthelist.type == "general") and (indexinthelist.index ==
+int(index[1:]))):
                      newamplitude.indexes.append(indexinthelist)
                      break
             else:
-               print("syntax error")
-               stop
+               raise SyntaxError(" ")
       amplitudes.append(newamplitude)
 
    newoperatorsequence = OperatorSequence(numericalfactor,summation,amplitudes,newsequences)
@@ -231,8 +247,7 @@ def stringtooperatorsequence(expression):
 def combinepermutations(one,two):
    """Connects two permutations of indexes"""
    if (len(one) != len(two)):
-      print("Internal error")
-      stop
+      raise SyntaxError(" ")
    three = []
    for n in range(len(one)/2):
       three.append(one[n])
@@ -261,7 +276,8 @@ def isidenticalto(one,two):
    for n in range(len(one)/2):
       found = 0
       for m in range(len(two)/2):
-         if ((one[n].isidenticalto(two[m])) and (one[n+len(one)/2].isidenticalto(two[m+len(two)/2]))):
+         if ((one[n].isidenticalto(two[m])) and
+(one[n+len(one)/2].isidenticalto(two[m+len(two)/2]))):
             found = 1
       if (not found):
          return 0
@@ -300,7 +316,8 @@ class Operator:
 
    def isidenticalto(self,another):
       """Checks if two second-quantized operators are identical"""
-      if ((self.type == another.type) and (self.dagger == another.dagger) and (self.index == another.index)):
+      if ((self.type == another.type) and (self.dagger == another.dagger) and (self.index ==
+another.index)):
          return 1
       else:
          return 0
@@ -346,14 +363,17 @@ class Operator:
          return 1
 
       # at this point, self.type = another.type
-      if ((not operatorsequence.summation.hastheindex(self)) and (not operatorsequence.summation.hastheindex(another))):
+      if ((not operatorsequence.summation.hastheindex(self)) and (not
+operatorsequence.summation.hastheindex(another))):
          if (self.index > another.index):
             return 1
          else:
             return 0
-      elif (operatorsequence.summation.hastheindex(self) and (not operatorsequence.summation.hastheindex(another))):
+      elif (operatorsequence.summation.hastheindex(self) and (not
+operatorsequence.summation.hastheindex(another))):
          return 0
-      elif ((not operatorsequence.summation.hastheindex(self)) and operatorsequence.summation.hastheindex(another)):
+      elif ((not operatorsequence.summation.hastheindex(self)) and
+operatorsequence.summation.hastheindex(another)):
          return 1
       else:
          # at this point, self.type = another.type and both are summed over
@@ -583,11 +603,13 @@ class Amplitude:
       # count the number of like amplitudes in operatorsequence
       nself = 0
       for amplitude in operatorsequence.amplitudes:
-         if ((amplitude.type == self.type) and (len(amplitude.indexes) == len(self.indexes)) and (amplitude.conjugate == self.conjugate)):
+         if ((amplitude.type == self.type) and (len(amplitude.indexes) == len(self.indexes)) and
+(amplitude.conjugate == self.conjugate)):
             nself = nself + 1
       nanother = 0
       for amplitude in operatorsequence.amplitudes:
-         if ((amplitude.type == another.type) and (len(amplitude.indexes) == len(another.indexes)) and (amplitude.conjugate == another.conjugate)):
+         if ((amplitude.type == another.type) and (len(amplitude.indexes) ==
+len(another.indexes)) and (amplitude.conjugate == another.conjugate)):
             nanother = nanother + 1
       if (nself > nanother):
          return 0
@@ -749,7 +771,7 @@ class Factor:
             raise RuntimeError("unrealistic factor")
       fraction = abs(int(1.0/coefficient))
       if (1.0/float(fraction) != abs(coefficient)):
-         print(" !!! WARNING !!! inaccurate arithmatic")
+         print(" !!! WARNING !!! inaccurate arithmetic")
       if (fraction == 1):
          frac = ""
       else:
@@ -810,7 +832,8 @@ class Factor:
 class OperatorSequence:
 
    def __init__(self,factor=[],summation=[],amplitudes=[],sequence=[]):
-      """Creates a sequence of normal ordered second-quantized operators with some numerical factor, amplitudes, and summation"""
+      """Creates a sequence of normal ordered second-quantized operators with some numerical
+factor, amplitudes, and summation"""
       self.factor = factor
       self.summation = summation
       self.amplitudes = amplitudes
@@ -872,7 +895,8 @@ class OperatorSequence:
       file.write("\n")
 
    def removeemptycurly(self):
-      """Eliminates all empty curly brackets (curly means a sequence of normal ordered operator in {})"""
+      """Eliminates all empty curly brackets (curly means a sequence of normal ordered operator
+in {})"""
 
       hasempty = 0
       for ncurly in range(len(self.sequence)):
@@ -900,7 +924,8 @@ class OperatorSequence:
          return 0
 
    def isunabletocontract(self):
-      """Counts the number of operators and determine if it is possible to give nonzero contraction at the end"""
+      """Counts the number of operators and determine if it is possible to give nonzero
+contraction at the end"""
 
       # count the number of hole/particle/general creation/annihilation operators
       nholecreation = 0
@@ -996,9 +1021,11 @@ class OperatorSequence:
                # delete leftmost from the summation indexes
                for index in newsequence.summation.indexes:
                   if (leftmost.isidenticalto(index)):
-                     del newsequence.summation.indexes[newsequence.summation.indexes.index(index)]
+                     del
+newsequence.summation.indexes[newsequence.summation.indexes.index(index)]
 
-               # count the number of operators between leftmost and the current operator and determine the parity
+               # count the number of operators between leftmost and the current operator and
+determine the parity
                length = len(self.sequence[0][1:]) + curly.index(operator)
                for anothercurly in self.sequence[1:]:
                   if (anothercurly == curly):
@@ -1050,9 +1077,11 @@ class OperatorSequence:
                # delete operator from the summation indexes
                for index in newsequence.summation.indexes:
                   if (operator.isidenticalto(index)):
-                     del newsequence.summation.indexes[newsequence.summation.indexes.index(index)]
+                     del
+newsequence.summation.indexes[newsequence.summation.indexes.index(index)]
 
-               # count the number of operators between leftmost and the current operator and determine the parity
+               # count the number of operators between leftmost and the current operator and
+determine the parity
                length = len(self.sequence[0][1:]) + curly.index(operator)
                for anothercurly in self.sequence[1:]:
                   if (anothercurly == curly):
@@ -1102,12 +1131,14 @@ class OperatorSequence:
       return result
 
    def performfullcontraction(self):
-      """Performs full contraction of a given operator sequence and returns a list of tensor contractions"""
+      """Performs full contraction of a given operator sequence and returns a list of tensor
+contractions"""
 
       print(self.show())
       print(" ... commencing full operator contraction")
 
-      # result will be a list of tensor contractions (operator sequence objects with empty operator sequence)
+      # result will be a list of tensor contractions (operator sequence objects with empty
+operator sequence)
       result = ListOperatorSequences()
       result.add(self)
 
@@ -1197,7 +1228,8 @@ class OperatorSequence:
             return 0
          else:
             for noperator in range(len(self.sequence[nsequence])):
-               if (not self.sequence[nsequence][noperator].isidenticalto(another.sequence[nsequence][noperator])):
+               if (not
+self.sequence[nsequence][noperator].isidenticalto(another.sequence[nsequence][noperator])):
                   return 0
       return 1
 
@@ -1238,7 +1270,8 @@ class OperatorSequence:
                for anotherindex in another.summation.indexes:
                   if ((not self.has(anotherindex)) and (index.issimilarto(anotherindex))):
 
-                     # at this point, we know that we should relabel anotherindex by index everywhere in another
+                     # at this point, we know that we should relabel anotherindex by index
+everywhere in another
                      for nyetanother in range(len(another.summation.indexes)):
                         yetanother = another.summation.indexes[nyetanother]
                         if (yetanother.isidenticalto(anotherindex)):
@@ -1249,7 +1282,8 @@ class OperatorSequence:
                         for nyetanother in range(len(amplitude.indexes)):
                            yetanother = amplitude.indexes[nyetanother]
                            if (yetanother.isidenticalto(anotherindex)):
-                              another.amplitudes[namplitude].indexes[nyetanother] = copy.deepcopy(index)
+                              another.amplitudes[namplitude].indexes[nyetanother] =
+copy.deepcopy(index)
 
                      for nsequence in range(len(another.sequence)):
                         sequence = another.sequence[nsequence]
@@ -1268,7 +1302,8 @@ class OperatorSequence:
                      for anotherindex in anotheramplitude.indexes:
                         if ((not self.has(anotherindex)) and (index.issimilarto(anotherindex))):
 
-                           # at this point, we know that we should relabel anotherindex by index everywhere in another
+                           # at this point, we know that we should relabel anotherindex by index
+everywhere in another
                            for nyetanother in range(len(another.summation.indexes)):
                               yetanother = another.summation.indexes[nyetanother]
                               if (yetanother.isidenticalto(anotherindex)):
@@ -1279,14 +1314,16 @@ class OperatorSequence:
                               for nyetanother in range(len(amplitude.indexes)):
                                  yetanother = amplitude.indexes[nyetanother]
                                  if (yetanother.isidenticalto(anotherindex)):
-                                    another.amplitudes[namplitude].indexes[nyetanother] = copy.deepcopy(index)
+                                    another.amplitudes[namplitude].indexes[nyetanother] =
+copy.deepcopy(index)
 
                            for nsequence in range(len(another.sequence)):
                               sequence = another.sequence[nsequence]
                               for nyetanother in range(len(sequence)):
                                  yetanother = sequence[nyetanother]
                                  if (yetanother.isidenticalto(anotherindex)):
-                                    another.sequence[nsequence][nyetanother] = copy.deepcopy(index)
+                                    another.sequence[nsequence][nyetanother] =
+copy.deepcopy(index)
 
                            return another
 
@@ -1298,7 +1335,8 @@ class OperatorSequence:
                      for anotherindex in anothersequence:
                         if ((not self.has(anotherindex)) and (index.issimilarto(anotherindex))):
 
-                           # at this point, we know that we should relabel anotherindex by index everywhere in another
+                           # at this point, we know that we should relabel anotherindex by index
+everywhere in another
                            for nyetanother in range(len(another.summation.indexes)):
                               yetanother = another.summation.indexes[nyetanother]
                               if (yetanother.isidenticalto(anotherindex)):
@@ -1309,14 +1347,16 @@ class OperatorSequence:
                               for nyetanother in range(len(amplitude.indexes)):
                                  yetanother = amplitude.indexes[nyetanother]
                                  if (yetanother.isidenticalto(anotherindex)):
-                                    another.amplitudes[namplitude].indexes[nyetanother] = copy.deepcopy(index)
+                                    another.amplitudes[namplitude].indexes[nyetanother] =
+copy.deepcopy(index)
 
                            for nsequence in range(len(another.sequence)):
                               sequence = another.sequence[nsequence]
                               for nyetanother in range(len(sequence)):
                                  yetanother = sequence[nyetanother]
                                  if (yetanother.isidenticalto(anotherindex)):
-                                    another.sequence[nsequence][nyetanother] = copy.deepcopy(index)
+                                    another.sequence[nsequence][nyetanother] =
+copy.deepcopy(index)
 
                            return another
 
@@ -1399,7 +1439,8 @@ class OperatorSequence:
    def merges(self,another):
       """Merges another operator sequence to itself when possible"""
 
-      # parity of a permutation can be computed as the product of parities of all pairwise permutations
+      # parity of a permutation can be computed as the product of parities of all pairwise
+permutations
       parity = 1.0
 
       # determine the parity for amplitudes
@@ -1525,7 +1566,8 @@ class OperatorSequence:
             if (operatorsequence.factor.permutations[ncoeff] == []):
                operatorsequence.factor.permutations[ncoeff] = newpermutation
             else:
-               operatorsequence.factor.permutations[ncoeff] = combinepermutations(newpermutation, operatorsequence.factor.permutations[ncoeff])
+               operatorsequence.factor.permutations[ncoeff] = combinepermutations(newpermutation,
+operatorsequence.factor.permutations[ncoeff])
       return result
 
    def canonicalize(self):
@@ -1601,7 +1643,8 @@ class OperatorSequence:
                continue
             if (indexa.index < indexb.index):
                swap = another.summation.indexes[nindexa]
-               another.summation.indexes[nindexa] = copy.deepcopy(another.summation.indexes[nindexb])
+               another.summation.indexes[nindexa] =
+copy.deepcopy(another.summation.indexes[nindexb])
                another.summation.indexes[nindexb] = copy.deepcopy(swap)
 
       return another
@@ -1624,7 +1667,8 @@ class OperatorSequence:
          return 0
 
    def isdisconnected(self,withrespectto=[]):
-      """Returns 1 if disconnected; if (withrespectto) connectivity among the given amplitude types is tested"""
+      """Returns 1 if disconnected; if (withrespectto) connectivity among the given amplitude
+types is tested"""
 
       if (self.alreadycontracted()):
 
@@ -1732,7 +1776,7 @@ class ListOperatorSequences:
       print("")
       for line in self.show():
          print(line)
-      return ""
+      return""
 
    def show(self):
       """Returns a human-friendly string of the content"""
@@ -1813,7 +1857,7 @@ class ListOperatorSequences:
       # pick up a pair of operator sequences
       for nsequencea in range(len(self.list)):
 #        if (verbose):
-#           print('processing ',nsequencea,' / ',range(len(self.list)))
+#           print 'processing ',nsequencea,' / ',range(len(self.list))
          sequencea = self.list[nsequencea]
          for nsequenceb in range(len(self.list)):
             sequenceb = self.list[nsequenceb]
@@ -1823,7 +1867,8 @@ class ListOperatorSequences:
                sequencec = sequencea.fullyrelabels(sequenceb)
                if (sequencea.canmerge(sequencec)):
                   self.add(sequencea.merges(sequencec))
-                  # It is extremely important that the following two statements are executed in this order
+                  # It is extremely important that the following two statements are executed in
+this order
                   del self.list[nsequenceb]
                   del self.list[nsequencea]
                   return self
@@ -1836,7 +1881,8 @@ class ListOperatorSequences:
                   for sequenced in permutation.list:
                      if (sequencea.canmerge(sequenced)):
                         self.add(sequencea.merges(sequenced))
-                        # It is extremely important that the following two statements are executed in this order
+                        # It is extremely important that the following two statements are
+executed in this order
                         del self.list[nsequenceb]
                         del self.list[nsequencea]
                         return self
@@ -1846,7 +1892,8 @@ class ListOperatorSequences:
                   for sequenced in permutation.list:
                      if (sequencea.canmerge(sequenced)):
                         self.add(sequencea.merges(sequenced))
-                        # It is extremely important that the following two statements are executed in this order
+                        # It is extremely important that the following two statements are
+executed in this order
                         del self.list[nsequenceb]
                         del self.list[nsequencea]
                         return self
@@ -1893,7 +1940,8 @@ class ListOperatorSequences:
                   for sequenced in permutation.list:
                      if (sequencea.canmerge(sequenced)):
                         self.add(sequencea.merges(sequenced))
-                        # It is extremely important that the following two statements are executed in this order
+                        # It is extremely important that the following two statements are
+executed in this order
                         del self.list[nsequenceb]
                         del self.list[nsequencea]
                         return self
@@ -1907,7 +1955,7 @@ class ListOperatorSequences:
          print(" ! Warning! a cyclic contraction is found")
 #        self.simplifythree(verbose)
       self.simplifytwo(verbose)
-      # the followings do not seem to affect the result, yet it costs enormous memory & time
+      # the following do not seem to affect the result, yet it costs enormous memory & time
       # self.simplifyfour(1)
       self = copy.deepcopy(self.deletezero())
       return self
@@ -2052,9 +2100,11 @@ class ListOperatorSequences:
       return self
 
    def performfullcontraction(self):
-      """Performs full contraction of a list of operator sequences and returns a list of tensor contractions"""
+      """Performs full contraction of a list of operator sequences and returns a list of tensor
+contractions"""
 
-      # result will be a list of tensor contractions (operator sequence objects with empty operator sequence)
+      # result will be a list of tensor contractions (operator sequence objects with empty
+operator sequence)
       self = self.simplifyone()
       result = ListOperatorSequences()
 
@@ -2268,6 +2318,7 @@ class ListOperatorSequences:
       newlength = len(result.list)
 
       if (originallength != newlength):
-         print(" !!! WARNING !!! %d computationally zero terms have been deleted" %(originallength - newlength))
+         print(" !!! WARNING !!! %d computationally zero terms have been deleted"
+%(originallength - newlength))
 
       return result

src/tce/splitfiles.py

@@ -1,32 +1,31 @@
-#!/usr/bin/env python3
 # Usage: python splitfiles.py < inputfile.F
 # (c) All rights reserved by Battelle & Pacific Northwest Nat'l Lab (2002)
 # $Id$
 
-import string
-import copy
 import sys
 
 source = sys.stdin.readlines()
 
 if (not source):
-   print("Usage: python splitfiles.py < inputfile.F")
+    print("Usage: python splitfiles.py < inputfile.F")
 
 nfiles = 0
+filename = " "
+filecontent = " "
 for line in source:
-   if (string.find(line,"SUBROUTINE") != -1):
+   if (line.find("SUBROUTINE") != -1):
       if (nfiles):
-         file = open(filename+".F","w")
+         file = open(filename+".F", "w")
          for newline in filecontent:
             file.write(newline)
-      filename = string.split(line[string.find(line,"SUBROUTINE")+11:],"(")[0]
+      filename = line[line.find("SUBROUTINE")+11:].split("(", 99999)[0]
       print(filename+".o\\")
       nfiles = nfiles + 1
       filecontent = [line]
    else:
       filecontent.append(line)
 # don't forget to dump the last subroutine
-file = open(filename+".F","w")
+file = open(filename+".F", "w")
 for newline in filecontent:
    file.write(newline)
-print("Number of files generated:",nfiles)
+print("Number of files generated:", nfiles)

src/util/errquit.h

@@ -1,6 +1,8 @@
+#ifndef _ERRQUIT_H
+#define _ERRQUIT_H
 // UERR - Not yet assigned to a category
 // UNKNOWN_ERR - Not yet assigned to a category
-// MEM_ERROR - Generic Memory error
+// MEM_ERR - Generic Memory error
 // RTDB_ERR - Error in the Runtime Database
 // INPUT_ERR - Error resulting from inproper user input
 // CAPMIS_ERR - Features that have not been implemented yet
@@ -31,3 +33,5 @@ const int DISK_ERR    = 100;
 const int CALC_ERR    = 110;
 const int FMM_ERR     = 120;
 // $Id$
+
+#endif