From 02608cf25a502f6777f33583bf00a4752d0319f4 Mon Sep 17 00:00:00 2001 From: Adam Parler Date: Wed, 29 Nov 2023 23:29:08 -0800 Subject: [PATCH 1/2] Added some scf files --- src/ddscf/comp4_ext.c | 32 +++++++++ src/ddscf/fock_2e_file.c | 142 +++++++++++++++++++++++++++++++++++++++ 2 files changed, 174 insertions(+) create mode 100644 src/ddscf/comp4_ext.c create mode 100644 src/ddscf/fock_2e_file.c diff --git a/src/ddscf/comp4_ext.c b/src/ddscf/comp4_ext.c new file mode 100644 index 0000000..b97284e --- /dev/null +++ b/src/ddscf/comp4_ext.c @@ -0,0 +1,32 @@ +#include +#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'); + +} diff --git a/src/ddscf/fock_2e_file.c b/src/ddscf/fock_2e_file.c new file mode 100644 index 0000000..fa0021c --- /dev/null +++ b/src/ddscf/fock_2e_file.c @@ -0,0 +1,142 @@ +#include "util.h" +#include "cscfps.h" +#include "cfock.h" +#include +#include + +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 +} + + + + +} From f0ec69874e0f3dbe58de80f94cfb141c2957cbb6 Mon Sep 17 00:00:00 2001 From: Adam Parler Date: Wed, 6 Dec 2023 02:37:36 -0800 Subject: [PATCH 2/2] Moved from C++ to C. Updated Python files. --- .gitignore | 7 +- src/basis/getlibr.py | 76 +++++++-- src/{nwchem.cpp => nwchem.c} | 6 +- src/rtdb/rtdb.h | 288 ++++++++++++++++++++++------------- src/tce/oce.py | 187 ++++++++++++++--------- src/tce/splitfiles.py | 17 +-- src/util/errquit.h | 6 +- 7 files changed, 384 insertions(+), 203 deletions(-) rename src/{nwchem.cpp => nwchem.c} (99%) diff --git a/.gitignore b/.gitignore index 71319bc..0062803 100644 --- a/.gitignore +++ b/.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/ diff --git a/src/basis/getlibr.py b/src/basis/getlibr.py index c5aefbf..215cf88 100644 --- a/src/basis/getlibr.py +++ b/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") diff --git a/src/nwchem.cpp b/src/nwchem.c similarity index 99% rename from src/nwchem.cpp rename to src/nwchem.c index 620af57..e98437a 100644 --- a/src/nwchem.cpp +++ b/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,*); } - + } diff --git a/src/rtdb/rtdb.h b/src/rtdb/rtdb.h index c50f473..069b115 100644 --- a/src/rtdb/rtdb.h +++ b/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] -// 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] -// 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 diff --git a/src/tce/oce.py b/src/tce/oce.py index bba8979..36153e3 100644 --- a/src/tce/oce.py +++ b/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 diff --git a/src/tce/splitfiles.py b/src/tce/splitfiles.py index 75613a3..66ebdb8 100644 --- a/src/tce/splitfiles.py +++ b/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) diff --git a/src/util/errquit.h b/src/util/errquit.h index 9c86292..53cc291 100644 --- a/src/util/errquit.h +++ b/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