From 7e8da872a143e165aca4ed14374257136cc3b144 Mon Sep 17 00:00:00 2001 From: Adam Parler <35127748+bigperm17@users.noreply.github.com> Date: Sun, 3 Oct 2021 11:49:12 -0700 Subject: [PATCH] Added a few files to begin translating --- src/basis/getlibr.py | 64 ++ src/nwchem.cpp | 403 ++++++++ src/rtdb/rtdb.h | 105 ++ src/tce/oce.py | 2273 +++++++++++++++++++++++++++++++++++++++++ src/tce/splitfiles.py | 32 + src/util/errquit.h | 33 + 6 files changed, 2910 insertions(+) create mode 100644 src/basis/getlibr.py create mode 100644 src/nwchem.cpp create mode 100644 src/rtdb/rtdb.h create mode 100644 src/tce/oce.py create mode 100644 src/tce/splitfiles.py create mode 100644 src/util/errquit.h diff --git a/src/basis/getlibr.py b/src/basis/getlibr.py new file mode 100644 index 0000000..c5aefbf --- /dev/null +++ b/src/basis/getlibr.py @@ -0,0 +1,64 @@ +#!/usr/bin/env 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/ +# Requires the installation of the python env. from +# https://github.com/MolSSI-BSE/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'] + #open file + # get rid of asterisks + file_name = bas_name.replace("*","s") + #get rid of parenthesis + file_name = file_name.replace("(","") + file_name = file_name.replace(")","") + #replace commas with underscore + file_name = file_name.replace(",","_") + #replace whitespace with underscore + file_name = file_name.replace(" ","_") + #replace forward slash with underscore + file_name = file_name.replace("/","_") + 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') + 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) + except: +# print("failed for"+element) + pass + else: + bs_str=bs_str.replace("BASIS","basis") + 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) + #ECP + bs_str=bs_str.replace("ECP","ecp \""+element_str+"_"+bas_name+"\"") + output_file.write(bs_str) + # + print(bas_name+" "+element_str) +print("end") + diff --git a/src/nwchem.cpp b/src/nwchem.cpp new file mode 100644 index 0000000..620af57 --- /dev/null +++ b/src/nwchem.cpp @@ -0,0 +1,403 @@ +#include +#include +#include "errquit.h" +#include "rtdb.h" +#ifdef USE_TCGMSG + #include "tcgmsg.h" +#else + int NODEID; + extern NODEID; +#endif +#include "pstat.h" +#include "util.h" +#include "inp.h" +#include "bgj_common.h" +#include "stdio.h" + int RTDB; + int STACK; + int HEAP; + int GLOBAL; + bool STATUS; + bool OVERIFY, OHARDFAIL; +#ifdef CRAY_T3D + int oldact, fsigctl; +#endif +#ifdef PSCALE + int IO_CODE; +#else + int32_t IO_CODE; +#endif + +// $Id$ + +// ====================================================================================================== +//> \mainpage Northwest Computational Chemistry Package (NWChem) 7.0.2 +//> +//> NWChem is an open-source computational chemistry package distributed under the terms of +//> the Educational Community License (ECL) 2.0 +//> +//> This software and its documentation were developed at the EMSL at Pacific Northwest National Laboratory, +//> a multiprogram national laboratory, operated for the U.S. Department of Energy by Battelle under +//> Contract Number DE-AC05-76RL01830. Support for this work was provided by the Department of Energy +//> Office of Biological and Environmental Research, Office of Basic Energy Science, and the Office of +//> Advanced Scientific Computing. +//> +//> Licensed under the Educational Community License, Version 2.0 (the "License"); you may +//> not use this file except in compliance with the License. You may obtain a copy of the +//> License at https://opensource.org/licenses/ECL-2.0. +//> +//> Unless required by applicable law or agreed to in writing, software distributed under the +//> License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, +//> either express or implied. See the License for the specific language governing +//> permissions and limitations under the License. +//> +//> Further information, including user documentation and forums, may be found at +//> http://www.nwchem-sw.org/. Alternatively, +//> the paper +//> +//> * M. Valiev, E.J. Bylaska, N. Govind, K. Kowalski, T.P. Straatsma, H.J.J. Van Dam, +//> D. Wang, J. Nieplocha, E. Apra, T.L. Windus, W.A. de Jong (2010)
+//> "NWChem: A comprehensive and scalable open-source solution for large scale molecular simulations"
+//> Computer Physics Communications, 181, 1477–1489, DOI: 10.1016/j.cpc.2010.04.018 +//> +//> provides details on the codes capabilities. +//> +//> Copyright (c) 1994-2020 Pacific Northwest National Laboratory, Battelle Memorial Institute +//> +//> Environmental Molecular Sciences Laboratory (EMSL)
+//> Pacific Northwest National Laboratory
+//> Richland, WA 99352 + +// ====================================================================================================== + + +int main(){ + char input_filename[nw_max_path_len], rtdb_name[nw_max_path_len]; + double total_wall, total_cpu; + #ifdef USE_OFFLOAD + int ppnout; + bool offload_enabled; + extern offload_enabled; + int offload_device; + extern offload_device; + #endif + + bool ostartup, ocontinue, orestart; + bool input_parse; + extern input_parse; + #if defined(USE_OPENMP) + int omp_get_max_threads; + extern omp_get_max_threads; + #endif + +// Create parallel processes and initialize IPC layer + + pbeginf(); + +// MXINIT is needed by PeIGS and PFFT to initialize +// the communication fabric they use. + + mxinit(); + +// Initialize timers so they are relative to job start + + total_wall = -util_wallsec(); + total_cpu = -util_cpusec(); + +// Only process 0 opens the input file +// (note that ga_nodeid() does not work yet!) + + if (nodeid() == 0){ + +// Get the name of the input file from the command line + + get_input_filename(input_filename); + + try { + FILE * LuIn = fopen(input_filename,'r'); + } + catch { + errquit('nwchem: failed to open the input file', 0, INPUT_ERR); + } + } + else{ + #if !(defined(KSR) || defined(IBM) || defined(FUJITSU_SOLARIS) ||defined(FUJITSU_VPP) ||defined(FUJITSU_VPP64)) + fclose(LuIn); + #endif + } + +// Look for memory directive in the input ... must eventually +// open the rtdb first so that can get memory directive out of that +// if it is not in the input + +// The user input model has well-defined categories of memory, +// each of which has a specific size. How we use these limits +// depends on the platform we are running on. + + input_mem_size(stack, heap, global, overify, ohardfail); + +// Initialize local memory allocator & global array tools + + ga_initialize_ltd(ma_sizeof(mt_dbl,global,mt_byte)); + // this must happen after GA and before MA + util_setup_gpu_affinity(); + if ( ga_uses_ma() ) { + if (!ma_init(mt_dbl, stack, heap+global)){ + errquit('nwchem.F: ma_init failed (ga_uses_ma=T)',911, MA_ERR); + } + } else{ + if (!ma_init(mt_dbl,stack,heap)) { + errquit('nwchem.F: ma_init failed (ga_uses_ma=F)',911, MA_ERR); + } + } + +// Touch OpenMP here so that any runtime initialization happens up-front. +// This ensures that any printout that the OpenMP runtime generates, +// such as affinity information, appears at the top of the output file. +// Otherwise, it might not appear until e.g. the CCSD module, at which +// point it will pollute the output file in an undesirable way. + +// Do not move this in front of GA/MPI/TCGMSG initialization, since the +// OpenMP runtime may inherit affinity information from MPI that is only +// determined during MPI initialization. + +// Format definition is outside of preprocessor protection to ensure the +// label is not accidentally reused, since that will not be caught by +// testing that does not enable OpenMP. + + g99 format(2x, 'NWChem w/ OpenMP: Maximum threads = ',i4); + #if defined(USE_OPENMP){ + #pragma omp parallel + #pragma omp master + { + if (ga_nodeid() == 0){ + write(luout,g99); + omp_get_max_threads(); + } + } + #endif +// set no. threads for threaded BLAS to 1 + util_blas_set_num_threads(1); + + rtdb_init() + +// More for amusement then efficiency force all MA allocated entities +// to be aligned at the beginning of a 128 byte cache line + +// if (!ma_set_numalign(7)){ +// errquit('nwchem.cpp: ma_set_numalign failed',911, MA_ERR); +// } +// aligned to 64byte record + if (!ma_set_numalign(6)){ + errquit('nwchem.cpp: ma_set_numalign failed', 911, MA_ERR); + } + + +//old:------------------------------------------------------- START --------- +//old:C GA allocations come out of MA space, so lump them together +//old:C and let MA impose the limits on GA sizes instead of actually +//old:C using the global limit. +//old:C +//old: if ( ga_uses_ma() ) then +//old: if (.not. ma_init(mt_dbl, stack, heap+global)) +//old: $ call errquit('nwchem: ma_init failed', -1) +//old: call ga_initialize +//old:C +//old:C GA allocations are separate from MA, so the separate limit +//old:C must be enforced. Note GA only understands bytes. +//old:C +//old: else +//old: if (.not. ma_init(mt_dbl, stack, heap)) +//old: $ call errquit('nwchem: ma_init failed', -1) +//old: call ga_initialize_ltd(ma_sizeof(mt_dbl, global, mt_byte) ) +//old: endif +//old:------------------------------------------------------- END ----------- + +//*** call nxtval_ga_initialize() + +// Trap SIGFPE after GA to override handler + +//*** call ieeetrap() + #if defined(LINUXALPHA) + dec_fpe(); // To avoid underflow problems on Alpha in Texas + #endif + #ifdef CRAY_T3D + // This as a temporary fix for SIGFPE in Texas that does not seem + // to affect the final results + oldact = fsigctl('IGNORE','SIGFPE',0); + #endif + #ifdef LINUX +// uncommenting this line turns on sigfpe trapping under linux +// linux_trapfpe(); + #endif + #ifdef MACX +// uncommenting this line turns on sigfpe trapping under Mac OSX +// macx_trapfpe(); + #endif +// Hard fail is good for development but means that we cannot +// respond to allocation problems. Disable by default. + status = ma_set_auto_verify(overify); + status = ma_set_hard_fail(ohardfail); + status = ma_set_error_print(ohardfail); + +// Initialize pstat + + if (!pstat_init(20,1,' ')){ + errquit('nwchem: pstat_init failed', 0, UNKNOWN_ERR); + } + + input_file_info(input_filename, rtdb_name, ostartup, ocontinue); + +// Now are ready to summarize the environment + + nwchem_banner(input_filename, rtdb_name, ostartup, ocontinue); + +// Actually open the database and store the file prefix + +// Note that only process 0 has the database name ... that is OK. + + if (ostartup){ + if (!rtdb_open(rtdb_name, 'empty', rtdb)){ + errquit('start: rtdb_open empty failed', 0, RTDB_ERR); + } + } else{ + if (!rtdb_open(rtdb_name, 'old', rtdb)){ + errquit('start: rtdb_open old failed', 0, RTDB_ERR); + } + } + + +// initialize nxtask + nxtask_init(rtdb); + +//!! BGJ + bgj_rtdb = rtdb; +//!! BGJ + + + if (ostartup || ocontinue){ + orestart = false; + } else{ + orestart = true; + } + + util_set_rtdb_state(rtdb, ostartup, ocontinue, orestart); + util_file_info_rtdb(rtdb); // Save file info for restart + movecs_ecce_print_on(); + geom_hnd_parallel(true) + perfm_start(); + + #ifdef USE_OFFLOAD + util_getppn(ppnout); + if (ppnout == 0){ + errquit('util_getppn failed',0,UERR); + } + if (ga_nodeid() == 0){ + write(luout,*) ga_nodeid(), ' ppn ', ppnout; + } + if (offload_enabled()){ + if (ga_nodeid() < ppnout){ + write(luout, '(I8,A,I2)') ga_nodeid(), ' offload enabled, GPU: ', + offload_device(); + } + } + ga_sync() + #endif + + if (orestart || ocontinue){ + nw_print_restart_info(rtdb); + } + +// if continue then go right to task stored on rtdb do not further parse +// input. if input is required then user should have used restart + + if (ocontinue){ + task(rtdb); + } + +// Parse input data, shove into database and execute tasks + + g10 if (input_parse(rtdb)){ // while(tasks to do) + util_print_rtdb_load(rtdb, ' '); // High level print + if (util_print('tcgmsg', print_never)){ + setdbg(1); + } else{ + setdbg(0); + } + #ifdef CATAMOUNT + util_allocga(); + #endif + + task(rtdb); + goto g10; // end while + } + +// Close the RTDB + + util_print_rtdb_load(rtdb, ' '); // High level print + if (util_print('rtdbvalues', print_debug)){ + if (!rtdb_print(rtdb,true)){ + errquit('control: rtdb_print failed', 0, RTDB_ERR); + } + } else if (util_print('rtdb', print_high)){ + if (! rtdb_print(rtdb, false)){ + errquit('control: rtdb_print failed', 0, RTDB_ERR); + } + } + + if (!rtdb_close(rtdb, 'keep')){ + errquit('nwchem: rtdb_close failed', rtdb, RTDB_ERR); + } + + if (util_print('rtdb', print_high) || util_print('rtdbvalues', print_high)){ + rtdb_print_usage(); // Called after closing so memory leaks apparent + } + +// Tidy up pstat + + if (!pstat_terminate()){ + errquit('nwchem: pstat_terminate failed', 0, UNKNOWN_ERR); + } + +//** nxtval_ga_terminate() + +// Print memory and other info + + ga_sync(); + if (ga_nodeid == 0){ + if (util_print('ga summary', print_default)){ + ga_summarize(0); + } + if (util_print('ga stats', print_default)){ + ga_print_stats(); + write(LuOut,*); + } + + } + + + + + + + + + + + + + + + + + + + + + + + + + + return 0; +} diff --git a/src/rtdb/rtdb.h b/src/rtdb/rtdb.h new file mode 100644 index 0000000..c50f473 --- /dev/null +++ b/src/rtdb/rtdb.h @@ -0,0 +1,105 @@ +// +// 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; diff --git a/src/tce/oce.py b/src/tce/oce.py new file mode 100644 index 0000000..bba8979 --- /dev/null +++ b/src/tce/oce.py @@ -0,0 +1,2273 @@ +# Operator Contraction Engine v.1.0 +# (c) All rights reserved by Battelle & Pacific Northwest Nat'l Lab (2002) +# $Id$ + +import string + +import copy + +import sys + +def readfromfile(filename): + """Converts the content of a file to a ListOperatorSequences object""" + + result = ListOperatorSequences() + file = open(filename,"r") + alwaystrue = 1 + while (alwaystrue): + line = file.readline() + if (line == ""): + file.close() + return result + else: + line = line[0:len(line)-1] + result.add(stringtooperatorsequence(line)) + +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), + # (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, + # (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} + + sequences = expression[expression.index("{"):] + expression = expression[0:expression.index("{")] + operatorlist = [] + # 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:] + if (sequences[len(sequences)-1] != "}"): + 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: + newsequence = [] + sequence = string.split(sequence) + for index in sequence: + if (index[len(index)-1] == "+"): + dagger = "creation" + index = index[0:len(index)-1] + else: + dagger = "annihilation" + 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)) + 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"))): + newsequence.append(Operator("general",dagger,string.ascii_letters.index(index)+1)) + else: + if (index[0] == "p"): + newsequence.append(Operator("particle",dagger,int(index[1:]))) + elif (index[0] == "h"): + newsequence.append(Operator("hole",dagger,int(index[1:]))) + elif (index[0] == "g"): + newsequence.append(Operator("general",dagger,int(index[1:]))) + else: + print("Syntax error: an operator not recognized") + stop + operatorlist = operatorlist + newsequence + newsequences.append(newsequence) + + breakdown = string.split(expression) + + # get a numerical factor if any + numericalfactor = "" + for element in breakdown: + if (((element[0] >= 'a') and (element[0] <= 'z')) or \ + ((element[0] >= 'A') and (element[0] <= 'Z'))): + break + numericalfactor = string.join([numericalfactor, element]) + if (numericalfactor == ""): + numericalfactor = Factor([1.0],[[]]) + else: + numericalfactor = eval(numericalfactor) + numericalfactor = Factor([numericalfactor],[[]]) + + # get a summation if any + summationindexes = "" + remainder = "" + join = 0 + for element in breakdown: + if ((element[0:3] == "SUM") or (element[0:3] == "sum") or (element[0:3] == "Sum")): + join = 1 + if (join == 1): + summationindexes = string.join([summationindexes,element]) + elif (join == 2): + remainder = string.join([remainder,element]) + if ((join == 1) and (")" in element)): + join = 2 + index = string.find(summationindexes,"sum") + if (index != -1): + summationindexes = summationindexes[0:index] + summationindexes[index+3:] + index = string.find(summationindexes,"SUM") + if (index != -1): + summationindexes = summationindexes[0:index] + summationindexes[index+3:] + index = string.find(summationindexes,"Sum") + if (index != -1): + summationindexes = summationindexes[0:index] + summationindexes[index+3:] + index = string.find(summationindexes,"(") + if (index != -1): + summationindexes = summationindexes[0:index] + summationindexes[index+1:] + index = string.find(summationindexes,")") + if (index != -1): + summationindexes = summationindexes[0:index] + summationindexes[index+1:] + summationindexes = string.split(summationindexes) + summation = Summation([]) + for index in summationindexes: + 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"))): + for indexinthelist in operatorlist: + 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)): + 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)): + summation.indexes.append(indexinthelist) + break + else: + if (index[0] == "p"): + for indexinthelist in operatorlist: + if ((indexinthelist.type == "particle") and (indexinthelist.index == int(index[1:]))): + summation.indexes.append(indexinthelist) + break + elif (index[0] == "h"): + for indexinthelist in operatorlist: + if ((indexinthelist.type == "hole") and (indexinthelist.index == int(index[1:]))): + summation.indexes.append(indexinthelist) + break + elif (index[0] == "g"): + for indexinthelist in operatorlist: + if ((indexinthelist.type == "general") and (indexinthelist.index == int(index[1:]))): + summation.indexes.append(indexinthelist) + break + else: + print("syntax error") + stop + + # get amplitudes + remainder = string.split(remainder,")") + tobeamplitudes = remainder[0:len(remainder)-1] + remainder = remainder[len(remainder)-1] + amplitudes = [] + for tobeamplitude in tobeamplitudes: + newamplitude = Amplitude() + tobeamplitude = string.split(tobeamplitude,"(") + type = tobeamplitude[0] + conjugate = 1 + lastdaggerposition = len(type) + for i in range(len(type)-1,-1,-1): + if (type[i] == "+"): + conjugate = - conjugate + lastdaggerposition = i + if (conjugate == -1): + newamplitude.conjugate = 1 + else: + newamplitude.conjugate = 0 + newamplitude.type = string.strip(type[0:lastdaggerposition]) + index = 0 + for amplitude in amplitudes: + if ((amplitude.type == newamplitude.type) and (amplitude.index > index)): + index = amplitude.index + newamplitude.index = index + 1 + tobeamplitudeindexes = string.split(tobeamplitude[1]) + newamplitude.indexes = [] + for index in tobeamplitudeindexes: + 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"))): + for indexinthelist in operatorlist: + 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)): + 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)): + newamplitude.indexes.append(indexinthelist) + break + else: + if (index[0] == "p"): + for indexinthelist in operatorlist: + 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:]))): + newamplitude.indexes.append(indexinthelist) + break + elif (index[0] == "g"): + for indexinthelist in operatorlist: + if ((indexinthelist.type == "general") and (indexinthelist.index == int(index[1:]))): + newamplitude.indexes.append(indexinthelist) + break + else: + print("syntax error") + stop + amplitudes.append(newamplitude) + + newoperatorsequence = OperatorSequence(numericalfactor,summation,amplitudes,newsequences) + return newoperatorsequence + +def combinepermutations(one,two): + """Connects two permutations of indexes""" + if (len(one) != len(two)): + print("Internal error") + stop + three = [] + for n in range(len(one)/2): + three.append(one[n]) + for n in range(len(one)/2,len(one)): + for m in range(len(two)/2): + if (one[n].isidenticalto(two[m])): + three.append(two[m+len(two)/2]) + return three + +def isidenticalto(one,two): + """Returns true if two permutations of indexes are identical""" + if ((one == []) and (two == [])): + return 1 + if (one == []): + for m in range(len(two)/2): + if (not two[m].isidenticalto(two[m+len(two)/2])): + return 0 + return 1 + if (two == []): + for m in range(len(one)/2): + if (not one[m].isidenticalto(one[m+len(one)/2])): + return 0 + return 1 + if (len(one) != len(two)): + return 0 + 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]))): + found = 1 + if (not found): + return 0 + return 1 + +class Operator: + + def __init__(self,type="unknown",dagger="unknown",index=0): + """Creates a second-quantized hole/particle/general creation/annihilation operator""" + self.type = type + self.dagger = dagger + self.index = index + + def __str__(self): + """Prints the content""" + return self.show() + + def show(self): + """Returns a human-friendly string of the content""" + show = string.join([self.type[0], repr(self.index)], "") + if (self.dagger == "creation"): + show = string.join([show, "+"], "") + return show + + def tex(self): + """Returns a LaTex form of output""" + show = string.join([self.type[0],"_{",repr(self.index),"}"], "") + if (self.dagger == "creation"): + show = string.join([show, "^{\dagger}"], "") + return show + + def duplicate(self): + """Returns a deepcopy of self""" + duplicate = Operator(self.type,self.dagger,self.index) + return duplicate + + 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)): + return 1 + else: + return 0 + + def issimilarto(self,another): + """Checks if two second-quantized operators are similar""" + if ((self.type == another.type) and (self.dagger == another.dagger)): + return 1 + else: + return 0 + + def isin(self,list): + """Returns true if an operator is in the list""" + for index in list: + if (self.isidenticalto(index)): + return 1 + return 0 + + def showwithoutdagger(self): + """Returns a human-friendly string of the content""" + show = string.join([self.type[0], repr(self.index)], "") + return show + + def texwithoutdagger(self): + """Returns a human-friendly string of the content""" + show = string.join([self.type[0],"_{",repr(self.index),"}"], "") + return show + + def isgreaterthan(self,another,operatorsequence): + """Returns true if self should be to the right of another in the canonical order""" + + if ((self.type == 'hole') and (another.type == 'particle')): + return 0 + elif ((self.type == 'hole') and (another.type == 'general')): + return 0 + elif ((self.type == 'particle') and (another.type == 'hole')): + return 1 + elif ((self.type == 'particle') and (another.type == 'general')): + return 0 + elif ((self.type == 'general') and (another.type == 'hole')): + return 1 + elif ((self.type == 'general') and (another.type == 'particle')): + return 1 + + # at this point, self.type = another.type + 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))): + return 0 + 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 + selfconnectivity = [] + anotherconnectivity = [] + for namplitude in range(len(operatorsequence.amplitudes)): + amplitude = operatorsequence.amplitudes[namplitude] + if (amplitude.hastheindex(self)): + selfconnectivity.append(namplitude) + for namplitude in range(len(operatorsequence.amplitudes)): + amplitude = operatorsequence.amplitudes[namplitude] + if (amplitude.hastheindex(another)): + anotherconnectivity.append(namplitude) + selfconnectivity.sort() + anotherconnectivity.sort() + if (selfconnectivity < anotherconnectivity): + return 1 + elif (anotherconnectivity < selfconnectivity): + return 0 + + return 0 + +class Summation: + + def __init__(self,indexes=[]): + """Creates a summation""" + self.indexes = indexes + + def __str__(self): + """Print the amplitude""" + return self.show() + + def show(self): + """Returns a human-friendly string of the content""" + show = "Sum (" + for index in self.indexes: + show = string.join([show, index.showwithoutdagger()]) + show = string.join([show,")"]) + return show + + def tex(self): + """Returns a LaTeX string of the content""" + show = "" + for index in self.indexes: + if (show): + show = string.join([show,","],"") + else: + show = "\\sum_{" + show = string.join([show,index.texwithoutdagger()]) + show = string.join([show,"}"]) + return show + + def duplicate(self): + """Returns a deepcopy of itself""" + duplicate = Summation([]) + for index in self.indexes: + duplicate.indexes.append(index.duplicate()) + return duplicate + + def hasthesameform(self,another): + """Checks if two summations have the same numbers of holes, particles, and generals""" + nself = 0 + for operator in self.indexes: + if (operator.type == "hole"): + nself = nself + 1 + nanother = 0 + for operator in another.indexes: + if (operator.type == "hole"): + nanother = nanother + 1 + if (nself != nanother): + return 0 + nself = 0 + for operator in self.indexes: + if (operator.type == "particle"): + nself = nself + 1 + nanother = 0 + for operator in another.indexes: + if (operator.type == "particle"): + nanother = nanother + 1 + if (nself != nanother): + return 0 + nself = 0 + for operator in self.indexes: + if (operator.type == "general"): + nself = nself + 1 + nanother = 0 + for operator in another.indexes: + if (operator.type == "general"): + nanother = nanother + 1 + if (nself != nanother): + return 0 + return 1 + + def isidenticalto(self,another): + """Returns true if two summations are identical""" + if (len(self.indexes) != len(another.indexes)): + return 0 + else: + for nindex in range(len(self.indexes)): + selfindex = self.indexes[nindex] + anotherindex = another.indexes[nindex] + if (not selfindex.isidenticalto(anotherindex)): + return 0 + return 1 + + def hastheindex(self,another): + """Returns true if the summation has the input index""" + has = 0 + for index in self.indexes: + if (index.isidenticalto(another)): + has = 1 + return has + +class Amplitude: + + def __init__(self,type="unknown",indexes=[],index=0,conjugate=0): + """Creates an integral/amplitude""" + self.type = type + self.indexes = indexes + self.index = index + self.conjugate = conjugate + + def __str__(self): + """Print the amplitude""" + return self.show() + + def show(self): + """Returns a human-friendly string of the content""" + show = self.type + if (self.conjugate): + show = string.join([show, "+"],"") + show = string.join([show, "("]) + for index in self.indexes: + show = string.join([show, index.showwithoutdagger()]) + show = string.join([show,")"]) + return show + + def hastheindex(self,another): + """Returns true if the summation has the input index""" + has = 0 + for index in self.indexes: + if (index.isidenticalto(another)): + has = 1 + return has + + def tex(self): + """Returns a LaTeX string of the content""" + show = self.type + show = string.join([show, "^{"]) + for index in self.indexes[0:len(self.indexes)/2]: + show = string.join([show, index.texwithoutdagger()]) + show = string.join([show,"}"]) + show = string.join([show, "_{"]) + for index in self.indexes[len(self.indexes)/2:len(self.indexes)]: + show = string.join([show, index.texwithoutdagger()]) + show = string.join([show,"}"]) + if (self.conjugate): + show = string.join(["\\left(",show,"\\right)^{\\dagger}"],"") + return show + + def duplicate(self): + """Returns a deepcopy of itself""" + duplicate = Amplitude(self.type,[],self.index,self.conjugate) + for index in self.indexes: + duplicate.indexes.append(index.duplicate()) + return duplicate + + def hasthesameform(self,another): + """Checks if two amplitude sets have the same numbers of holes, particles, and generals""" + if (self.type != another.type): + return 0 + if (len(self.indexes) != len(another.indexes)): + return 0 + if (self.conjugate != another.conjugate): + return 0 + nself = 0 + for operator in self.indexes: + if (operator.type == "hole"): + nself = nself + 1 + nanother = 0 + for operator in another.indexes: + if (operator.type == "hole"): + nanother = nanother + 1 + if (nself != nanother): + return 0 + nself = 0 + for operator in self.indexes: + if (operator.type == "particle"): + nself = nself + 1 + nanother = 0 + for operator in another.indexes: + if (operator.type == "particle"): + nanother = nanother + 1 + if (nself != nanother): + return 0 + nself = 0 + for operator in self.indexes: + if (operator.type == "general"): + nself = nself + 1 + nanother = 0 + for operator in another.indexes: + if (operator.type == "general"): + nanother = nanother + 1 + if (nself != nanother): + return 0 + return 1 + + def isidenticalto(self,another): + """Returns true if two amplitudes are identical""" + if (self.type != another.type): + return 0 + elif (len(self.indexes) != len(another.indexes)): + return 0 + elif (self.conjugate != another.conjugate): + return 0 + else: + for nindex in range(len(self.indexes)): + selfindex = self.indexes[nindex] + anotherindex = another.indexes[nindex] + if (not selfindex.isidenticalto(anotherindex)): + return 0 + return 1 + + def isgreaterthan(self,another,operatorsequence): + """Returns true if self should be to the right of another in the canonical order""" + + # 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)): + 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)): + nanother = nanother + 1 + if (nself > nanother): + return 0 + elif (nself < nanother): + return 1 + # conjugate + if (self.conjugate > another.conjugate): + return 0 + elif (self.conjugate < another.conjugate): + return 1 + # type + if (self.type > another.type): + return 1 + elif (self.type < another.type): + return 0 + # number of indexes + if (len(self.indexes) < len(another.indexes)): + return 1 + elif (len(self.indexes) > len(another.indexes)): + return 0 + # number of external indexes + nself = 0 + iself = 9999999999 + for operator in self.indexes: + if (not operatorsequence.summation.hastheindex(operator)): + nself = nself + 1 + if (iself > operator.index): + iself = operator.index + nanother = 0 + ianother = 9999999999 + for operator in another.indexes: + if (not operatorsequence.summation.hastheindex(operator)): + nanother = nanother + 1 + if (ianother > operator.index): + ianother = operator.index + if (nself < nanother): + return 1 + elif (nself > nanother): + return 0 + # earliest external indexes +# if (nself > 0): +# if (ianother > iself): +# return 0 +# elif (ianother < iself): +# return 1 + # connectivity + selfconnectivity = [] + anotherconnectivity = [] + for operator in self.indexes: + if (operatorsequence.summation.hastheindex(operator)): + for amplitude in operatorsequence.amplitudes: + if (amplitude.hastheindex(operator)): + amplitudesymbol = amplitude.type + repr(len(amplitude.indexes)) + if (amplitude.conjugate): + amplitudesymbol = amplitudesymbol + "+" + selfconnectivity.append(amplitudesymbol) + for operator in another.indexes: + if (operatorsequence.summation.hastheindex(operator)): + for amplitude in operatorsequence.amplitudes: + if (amplitude.hastheindex(operator)): + amplitudesymbol = amplitude.type + repr(len(amplitude.indexes)) + if (amplitude.conjugate): + amplitudesymbol = amplitudesymbol + "+" + anotherconnectivity.append(amplitudesymbol) + selfconnectivity.sort() + anotherconnectivity.sort() + if (selfconnectivity < anotherconnectivity): + return 1 + elif (anotherconnectivity < selfconnectivity): + return 0 + return 0 + + def canonicalize(self,operatorsequence): + """Reorder the indexes in the canonical order""" + + another = self.duplicate() + parity = 1 + done = 0 + while (not done): + done = 1 + # reorder super indexes + for noperatora in range(len(another.indexes)/2): + for noperatorb in range(len(another.indexes)/2): + if (noperatora >= noperatorb): + continue + operatora = another.indexes[noperatora] + operatorb = another.indexes[noperatorb] + if (operatora.isgreaterthan(operatorb,operatorsequence)): + another.indexes[noperatorb] = copy.deepcopy(operatora) + another.indexes[noperatora] = copy.deepcopy(operatorb) + parity = parity * (-1) + done = 0 + done = 0 + while (not done): + done = 1 + # reorder sub indexes + for noperatora in range(len(another.indexes)/2,len(another.indexes)): + for noperatorb in range(len(another.indexes)/2,len(another.indexes)): + if (noperatora >= noperatorb): + continue + operatora = another.indexes[noperatora] + operatorb = another.indexes[noperatorb] + if (operatora.isgreaterthan(operatorb,operatorsequence)): + another.indexes[noperatorb] = copy.deepcopy(operatora) + another.indexes[noperatora] = copy.deepcopy(operatorb) + parity = parity * (-1) + done = 0 + + return [another,parity] + +class Factor: + + def __init__(self,coefficients=[],permutations=[]): + """Creates a numerical and permutation factor of an operator sequence""" + self.coefficients = coefficients + self.permutations = copy.deepcopy(permutations) + + def __str__(self): + """Prints the content""" + return self.show() + + def show(self): + """Returns a human-friendly string of contests""" + show = "[" + for n in range(len(self.coefficients)): + coefficient = self.coefficients[n] + # str() rounds a float after 12 digits, while repr() after 17, + # so the former tends to give a more pleasant expression. +# num = rationaltofractional(coefficient)[0] +# den = rationaltofractional(coefficient)[1] +# if (num >= 0): +# show = string.join([show,"+",repr(num)]) +# elif (num < 0): +# show = string.join([show,"-",repr(-num)]) +# if (den != 1): +# show = string.join([show,"/",repr(den)],"") + if (coefficient >= 0.0): + show = string.join([show,"+",str(coefficient)]) + elif (coefficient < 0.0): + show = string.join([show,"-",str(-coefficient)]) + if (self.permutations[n]): + show = string.join([show,"* P("]) + for noperator in range(len(self.permutations[n])/2): + operator = self.permutations[n][noperator] + show = string.join([show,operator.showwithoutdagger()]) + show = string.join([show,"=>"]) + for noperator in range(len(self.permutations[n])/2,len(self.permutations[n])): + operator = self.permutations[n][noperator] + show = string.join([show,operator.showwithoutdagger()]) + show = string.join([show,")"]) + show = string.join([show,"]"]) + return show + + def tex(self): + """Returns a LaTeX string of contests""" + coefficient = self.coefficients[0] + for n in range(len(self.coefficients)): + if (abs(self.coefficients[n]) != abs(coefficient)): + raise RuntimeError("unrealistic factor") + fraction = abs(int(1.0/coefficient)) + if (1.0/float(fraction) != abs(coefficient)): + print(" !!! WARNING !!! inaccurate arithmatic") + if (fraction == 1): + frac = "" + else: + frac = string.join(["\\frac{1}{",str(fraction),"}"],"") + if (coefficient >= 0.0): + show = string.join(["+",frac]) + elif (coefficient < 0.0): + show = string.join(["-",frac]) + if (len(self.coefficients) > 1): + show = string.join([show,"\\left("],"") + for n in range(len(self.coefficients)): + if (self.coefficients[n]/coefficient > 0.0): + show = string.join([show,"+"],"") + else: + show = string.join([show,"-"],"") + if (self.permutations[n]): + show = string.join([show,"P^{"]) + for nindex in range(len(self.permutations[n])/2,3*len(self.permutations[n])/4): + index = self.permutations[n][nindex] + show = string.join([show,index.texwithoutdagger()]) + for nindex in range(len(self.permutations[n])/4,len(self.permutations[n])/2): + index = self.permutations[n][nindex] + show = string.join([show,index.texwithoutdagger()]) + show = string.join([show,"}_{"]) + for nindex in range(len(self.permutations[n])/4): + index = self.permutations[n][nindex] + show = string.join([show,index.texwithoutdagger()]) + for nindex in range(3*len(self.permutations[n])/4,len(self.permutations[n])): + index = self.permutations[n][nindex] + show = string.join([show,index.texwithoutdagger()]) + show = string.join([show,"}"]) + else: + show = string.join([show,"1"],"") + show = string.join([show,"\\right)"]) + return show + + def multiply(self,factor): + """Multiply a factor to all coefficients""" + for n in range(len(self.coefficients)): + self.coefficients[n] = self.coefficients[n] * factor + + def add(self,another,factor=1.0): + """Add two Factors together""" + for m in range(len(another.coefficients)): + done = 0 + for n in range(len(self.coefficients)): + if (isidenticalto(self.permutations[n],another.permutations[m])): + if ((self.coefficients[n] < 0.0) and (another.coefficients[m] * factor > 0.0)): + print(" ! Warning ! cancellation of terms occurred ") + if ((self.coefficients[n] > 0.0) and (another.coefficients[m] * factor < 0.0)): + print(" ! Warning ! cancellation of terms occurred ") + self.coefficients[n] = self.coefficients[n] + another.coefficients[m] * factor + done = 1 + if (not done): + self.coefficients.append(another.coefficients[m] * factor) + self.permutations.append(another.permutations[m]) + +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""" + self.factor = factor + self.summation = summation + self.amplitudes = amplitudes + self.sequence = sequence + + def __str__(self): + """Prints the sequence of operator contractions""" + return self.show() + + def show(self): + """Returns a human-friendly string of the content""" + show = self.factor.show() + if (self.summation): + if (len(self.summation.indexes) > 0): + show = string.join([show, "*", self.summation.show()]) + for index in self.amplitudes: + show = string.join([show, "*", index.show()]) + if (self.sequence): + show = string.join([show, "* <0|"]) + for sequence in self.sequence: + show = string.join([show, "{"]) + for operator in sequence: + show = string.join([show, operator.show()]) + show = string.join([show, "}"]) + show = string.join([show, "|0>"]) + return show + + def tex(self): + """Returns a LaTeX string of the content""" + show = self.factor.tex() +# if (self.summation): +# if (len(self.summation.indexes) > 0): +# show = string.join([show, self.summation.tex()]) + for index in self.amplitudes: + show = string.join([show, index.tex()]) + if (self.sequence): + show = string.join([show, "\\langle 0 |"]) + for sequence in self.sequence: + show = string.join([show, "\{"]) + for operator in sequence: + show = string.join([show, operator.tex()]) + show = string.join([show, "\}"]) + show = string.join([show, "|0\\rangle"]) + return show + + def duplicate(self): + """Makes a copy of itself""" + duplicate = OperatorSequence() + duplicate.factor = copy.deepcopy(self.factor) + duplicate.summation = copy.deepcopy(self.summation) + duplicate.amplitudes = copy.deepcopy(self.amplitudes) + duplicate.sequence = copy.deepcopy(self.sequence) + return duplicate + + def writetofile(self,filename): + """Writes the output to a given file""" + file = open(filename,"w") + file.write(self.show()) + file.write("\n") + + def removeemptycurly(self): + """Eliminates all empty curly brackets (curly means a sequence of normal ordered operator in {})""" + + hasempty = 0 + for ncurly in range(len(self.sequence)): + curly = self.sequence[ncurly] + if (not curly): + del self.sequence[ncurly] + hasempty = 1 + break + + if (hasempty): + self.removeemptycurly() + else: + return self + + def alreadycontracted(self): + """Checks if an operator sequence object is fully contracted""" + + # first, we delete all empty {} just in case + self.removeemptycurly() + + # already fully contracted? + if (not self.sequence): + return 1 + else: + return 0 + + def isunabletocontract(self): + """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 + nholeannihilation = 0 + nparticlecreation = 0 + nparticleannihilation = 0 + ngeneralcreation = 0 + ngeneralannihilation = 0 + for sequence in self.sequence: + for operator in sequence: + if ((operator.type == "hole") and (operator.dagger == "creation")): + nholecreation = nholecreation + 1 + elif ((operator.type == "hole") and (operator.dagger == "annihilation")): + nholeannihilation = nholeannihilation + 1 + if ((operator.type == "particle") and (operator.dagger == "creation")): + nparticlecreation = nparticlecreation + 1 + elif ((operator.type == "particle") and (operator.dagger == "annihilation")): + nparticleannihilation = nparticleannihilation + 1 + if ((operator.type == "general") and (operator.dagger == "creation")): + ngeneralcreation = ngeneralcreation + 1 + elif ((operator.type == "general") and (operator.dagger == "annihilation")): + ngeneralannihilation = ngeneralannihilation + 1 + + # see if enough operators remain for contractions to survive + uncontractable = 0 + if (nholecreation + ngeneralcreation < nholeannihilation): + uncontractable = 1 + if (nholeannihilation + ngeneralannihilation < nholecreation): + uncontractable = 1 + if (nparticlecreation + ngeneralcreation < nparticleannihilation): + uncontractable = 1 + if (nparticleannihilation + ngeneralannihilation < nparticlecreation): + uncontractable = 1 + return uncontractable + + def performcontraction(self): + """Perform a contraction of the left-most operator""" + + # result will be a list of new operator sequence objects + result = ListOperatorSequences() + + # already fully contracted? + if (self.alreadycontracted()): + newsequence = self.duplicate() + result.add(newsequence) + return result + + # no way to contract? + elif ((len(self.sequence) == 1) or (self.isunabletocontract())): + return result + + # get the left-most operator + leftmost = self.sequence[0][0] + + # loop over other {} + for ncurly in range(len(self.sequence)): + curly = self.sequence[ncurly] + if (ncurly == 0): + continue + for noperator in range(len(self.sequence[ncurly])): + operator = curly[noperator] + + # only allowed contractions are {h+}{h} and {p}{p+} + if (leftmost.dagger == operator.dagger): + continue + elif ((leftmost.type == "hole") and (operator.type == "particle")): + continue + elif ((leftmost.type == "particle") and (operator.type == "hole")): + continue + elif ((leftmost.type == "hole") and (leftmost.dagger == "annihilation")): + continue + elif ((leftmost.type == "particle") and (leftmost.dagger == "creation")): + continue + elif ((operator.type == "hole") and (operator.dagger == "creation")): + continue + elif ((operator.type == "particle") and (operator.dagger == "annihilation")): + continue + + # check if the indexes can be made to match by virtue of summation + exist = "neither" + for index in self.summation.indexes: + if (leftmost.isidenticalto(index)): + exist = "leftmost" + if (exist == "neither"): + for index in self.summation.indexes: + if (operator.isidenticalto(index)): + exist = "operator" + if (exist == "leftmost"): + + # now contraction is possible --- add a new operator sequence object to result + newsequence = self.duplicate() + + # 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)] + + # 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): + break + else: + length = length + len(anothercurly) + parity = (-1)**length + newsequence.factor.multiply(parity) + + # delete the contracted pair from the sequence + del newsequence.sequence[0][0] + del newsequence.sequence[ncurly][noperator] + + # replace any appearance of leftmost by operator + if (operator.type == 'general'): + for namplitude in range(len(newsequence.amplitudes)): + amplitude = newsequence.amplitudes[namplitude] + for nindex in range(len(amplitude.indexes)): + index = amplitude.indexes[nindex] + if (operator.isidenticalto(index)): + newsequence.amplitudes[namplitude].indexes[nindex] = leftmost + for nindex in range(len(newsequence.summation.indexes)): + index = newsequence.summation.indexes[nindex] + if (operator.isidenticalto(index)): + newsequence.summation.indexes[nindex] = leftmost + else: + for namplitude in range(len(newsequence.amplitudes)): + amplitude = newsequence.amplitudes[namplitude] + for nindex in range(len(amplitude.indexes)): + index = amplitude.indexes[nindex] + if (leftmost.isidenticalto(index)): + newsequence.amplitudes[namplitude].indexes[nindex] = operator + for nindex in range(len(newsequence.summation.indexes)): + index = newsequence.summation.indexes[nindex] + if (leftmost.isidenticalto(index)): + newsequence.summation.indexes[nindex] = operator + + # cleanup the empty brackets + newsequence.removeemptycurly() + + # add to the result + result.add(newsequence) + + elif (exist == "operator"): + + # contraction is again possible --- add a new operator sequence object to result + newsequence = self.duplicate() + + # 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)] + + # 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): + break + else: + length = length + len(anothercurly) + parity = (-1)**length + newsequence.factor.multiply(parity) + + # delete the contracted pair from the sequence + del newsequence.sequence[0][0] + del newsequence.sequence[ncurly][noperator] + + # replace any appearance of operator by leftmost + if (leftmost.type == 'general'): + for namplitude in range(len(newsequence.amplitudes)): + amplitude = newsequence.amplitudes[namplitude] + for nindex in range(len(amplitude.indexes)): + index = amplitude.indexes[nindex] + if (leftmost.isidenticalto(index)): + newsequence.amplitudes[namplitude].indexes[nindex] = operator + for nindex in range(len(newsequence.summation.indexes)): + index = newsequence.summation.indexes[nindex] + if (leftmost.isidenticalto(index)): + newsequence.summation.indexes[nindex] = operator + else: + for namplitude in range(len(newsequence.amplitudes)): + amplitude = newsequence.amplitudes[namplitude] + for nindex in range(len(amplitude.indexes)): + index = amplitude.indexes[nindex] + if (operator.isidenticalto(index)): + newsequence.amplitudes[namplitude].indexes[nindex] = leftmost + for nindex in range(len(newsequence.summation.indexes)): + index = newsequence.summation.indexes[nindex] + if (operator.isidenticalto(index)): + newsequence.summation.indexes[nindex] = leftmost + + # cleanup the empty brackets + newsequence.removeemptycurly() + + # add to the result + result.add(newsequence) + + else: + break + + return result + + def performfullcontraction(self): + """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 = ListOperatorSequences() + result.add(self) + + # see if already fully contracted + done = self.alreadycontracted() + + # recursive execution of performcontraction() + iteration = 0 + while (not done): + iteration = iteration + 1 + newresult = ListOperatorSequences() + for halfwaycontracted in result.list: + newaddition = halfwaycontracted.performcontraction() + if (newaddition): + newresult.join(newaddition) + newresult.simplifyone() + numberofterms = len(newresult.list) + print(" ... iteration = %d, number of terms = %d" %(iteration, numberofterms)) + done = 1 + for halfwaycontracted in newresult.list: + if (not halfwaycontracted.alreadycontracted()): + done = 0 + result = newresult.duplicate() + + return result + + def hasthesameform(self,another): + """Checks if two operator sequences have the same form for possible consolidation""" + if (not self.summation.hasthesameform(another.summation)): + return 0 + if (len(self.amplitudes) != len(another.amplitudes)): + return 0 + else: + for namplitude in range(len(self.amplitudes)): + if (not self.amplitudes[namplitude].hasthesameform(another.amplitudes[namplitude])): + return 0 + if (len(self.sequence) != len(another.sequence)): + return 0 + else: + for nsequence in range(len(self.sequence)): + nself = 0 + for operator in self.sequence[nsequence]: + if (operator.type == "hole"): + nself = nself + 1 + nanother = 0 + for operator in another.sequence[nsequence]: + if (operator.type == "hole"): + nanother = nanother + 1 + if (nself != nanother): + return 0 + nself = 0 + for operator in self.sequence[nsequence]: + if (operator.type == "particle"): + nself = nself + 1 + nanother = 0 + for operator in another.sequence[nsequence]: + if (operator.type == "particle"): + nanother = nanother + 1 + if (nself != nanother): + return 0 + nself = 0 + for operator in self.sequence[nsequence]: + if (operator.type == "general"): + nself = nself + 1 + nanother = 0 + for operator in another.sequence[nsequence]: + if (operator.type == "general"): + nanother = nanother + 1 + if (nself != nanother): + return 0 + return 1 + + def isidenticalto(self,another): + """Returns true if two operator sequences are identical except for the factor""" + + if (not self.summation.isidenticalto(another.summation)): + return 0 + if (len(self.amplitudes) != len(another.amplitudes)): + return 0 + if (len(self.sequence) != len(another.sequence)): + return 0 + for namplitude in range(len(self.amplitudes)): + if (not self.amplitudes[namplitude].isidenticalto(another.amplitudes[namplitude])): + return 0 + for nsequence in range(len(self.sequence)): + if (len(self.sequence[nsequence]) != len(another.sequence[nsequence])): + return 0 + else: + for noperator in range(len(self.sequence[nsequence])): + if (not self.sequence[nsequence][noperator].isidenticalto(another.sequence[nsequence][noperator])): + return 0 + return 1 + + def has(self,index): + """Checks if a certain index is included in an operator sequence""" + + # see if the index is in summation indexes + for another in self.summation.indexes: + if (another.isidenticalto(index)): + return 1 + + # see if the index is in amplitude indexes + for amplitude in self.amplitudes: + for another in amplitude.indexes: + if (another.isidenticalto(index)): + return 1 + + # see if the index is in the operator sequences + for sequence in self.sequence: + for another in sequence: + if (another.isidenticalto(index)): + return 1 + + # not included + return 0 + + def relabels(self,another): + """Relabels the operator indexes to help consolidate terms""" + + if (not self.hasthesameform(another)): + return another + + else: + + # find a lone index in summation indexes + for index in self.summation.indexes: + if (not another.has(index)): + 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 + for nyetanother in range(len(another.summation.indexes)): + yetanother = another.summation.indexes[nyetanother] + if (yetanother.isidenticalto(anotherindex)): + another.summation.indexes[nyetanother] = copy.deepcopy(index) + + for namplitude in range(len(another.amplitudes)): + amplitude = another.amplitudes[namplitude] + for nyetanother in range(len(amplitude.indexes)): + yetanother = amplitude.indexes[nyetanother] + if (yetanother.isidenticalto(anotherindex)): + 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) + + return another + + # find a lone index in amplitude indexes + for selfamplitude in self.amplitudes: + for index in selfamplitude.indexes: + if (not another.has(index)): + for anotheramplitude in another.amplitudes: + 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 + for nyetanother in range(len(another.summation.indexes)): + yetanother = another.summation.indexes[nyetanother] + if (yetanother.isidenticalto(anotherindex)): + another.summation.indexes[nyetanother] = copy.deepcopy(index) + + for namplitude in range(len(another.amplitudes)): + amplitude = another.amplitudes[namplitude] + for nyetanother in range(len(amplitude.indexes)): + yetanother = amplitude.indexes[nyetanother] + if (yetanother.isidenticalto(anotherindex)): + 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) + + return another + + # find a lone index in operator sequences + for selfsequence in self.sequence: + for index in selfsequence: + if (not another.has(index)): + for anothersequence in another.sequence: + 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 + for nyetanother in range(len(another.summation.indexes)): + yetanother = another.summation.indexes[nyetanother] + if (yetanother.isidenticalto(anotherindex)): + another.summation.indexes[nyetanother] = copy.deepcopy(index) + + for namplitude in range(len(another.amplitudes)): + amplitude = another.amplitudes[namplitude] + for nyetanother in range(len(amplitude.indexes)): + yetanother = amplitude.indexes[nyetanother] + if (yetanother.isidenticalto(anotherindex)): + 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) + + return another + + return another + + def fullyrelabels(self,another): + """Relabels the operator indexes to help consolidate terms""" + + if (not self.hasthesameform(another)): + return another + else: + done = 0 + while (not done): + another = self.relabels(another) + done = self.hasnomismatch(another) + + return another + + def hasnomismatch(self,another): + """Returns 1 if there is no index in self that does not exist in another""" + + if (not self.hasthesameform(another)): + return 0 + + nomismatch = 1 + for index in self.summation.indexes: + if (not another.has(index)): + nomismatch = 0 + for selfamplitude in self.amplitudes: + for index in selfamplitude.indexes: + if (not another.has(index)): + nomismatch = 0 + for selfsequence in self.sequence: + for index in selfsequence: + if (not another.has(index)): + nomismatch = 0 + + return nomismatch + + def canmerge(self,another): + """Returns 1 if another operator sequence can be merged to itself""" + + # do they have any index mismatch? + if (not self.hasnomismatch(another)): + return 0 + + # do they have the identical operator sequences? + for nsequence in range(len(self.sequence)): + selfsequence = self.sequence[nsequence] + anothersequence = another.sequence[nsequence] + for noperator in range(len(selfsequence)): + selfoperator = selfsequence[noperator] + anotheroperator = anothersequence[noperator] + if (not selfoperator.isidenticalto(anotheroperator)): + return 0 + + # do they have the summation indexes that do not differ by more than just permutation? + for selfindex in self.summation.indexes: + exist = 0 + for anotherindex in another.summation.indexes: + if (anotherindex.isidenticalto(selfindex)): + exist = 1 + if (not exist): + return 0 + + # do they have the amplitude indexes that do not differ by more than just permutation? + for namplitude in range(len(self.amplitudes)): + selfamplitude = self.amplitudes[namplitude].indexes + anotheramplitude = another.amplitudes[namplitude].indexes + for selfindex in selfamplitude: + exist = 0 + for anotherindex in anotheramplitude: + if (anotherindex.isidenticalto(selfindex)): + exist = 1 + if (not exist): + return 0 + + return 1 + + 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 = 1.0 + + # determine the parity for amplitudes + for namplitude in range(len(self.amplitudes)): + selfamplitude = self.amplitudes[namplitude] + anotheramplitude = another.amplitudes[namplitude] + for nselfindexa in range(len(selfamplitude.indexes)): + selfindexa = selfamplitude.indexes[nselfindexa] + for nselfindexb in range(len(selfamplitude.indexes)): + if (nselfindexb <= nselfindexa): + continue + selfindexb = selfamplitude.indexes[nselfindexb] + for nanotherindexa in range(len(anotheramplitude.indexes)): + anotherindexa = anotheramplitude.indexes[nanotherindexa] + if (anotherindexa.isidenticalto(selfindexa)): + for nanotherindexb in range(len(anotheramplitude.indexes)): + anotherindexb = anotheramplitude.indexes[nanotherindexb] + if (anotherindexb.isidenticalto(selfindexb)): + if (nanotherindexb < nanotherindexa): + parity = parity * (-1.0) + + self.factor.add(another.factor, parity) + + return self + + def swapoperators(self,indexa,indexb): + """Swap indexa and indexb everywhere they appear in self""" + + for nindex in range(len(self.summation.indexes)): + index = self.summation.indexes[nindex] + if (index.isidenticalto(indexa)): + self.summation.indexes[nindex] = indexb + elif (index.isidenticalto(indexb)): + self.summation.indexes[nindex] = indexa + + for namplitude in range(len(self.amplitudes)): + amplitude = self.amplitudes[namplitude] + for nindex in range(len(amplitude.indexes)): + index = amplitude.indexes[nindex] + if (index.isidenticalto(indexa)): + self.amplitudes[namplitude].indexes[nindex] = indexb + elif (index.isidenticalto(indexb)): + self.amplitudes[namplitude].indexes[nindex] = indexa + + for nsequence in range(len(self.sequence)): + sequence = self.sequence[nsequence] + for nindex in range(len(sequence)): + index = sequence[nindex] + if (index.isidenticalto(indexa)): + self.sequence[nsequence][nindex] = indexb + elif (index.isidenticalto(indexb)): + self.sequence[nsequence][nindex] = indexa + + return self + + def swapamplitudes(self,namplitudea,namplitudeb): + """Swap two amplitudes in self""" + + swap = copy.deepcopy(self.amplitudes[namplitudea]) + self.amplitudes[namplitudea] = copy.deepcopy(self.amplitudes[namplitudeb]) + self.amplitudes[namplitudeb] = copy.deepcopy(swap) + + return self + + def targetindexpermutation(self): + """Returns a list of all possible permutations and redundancy of target indexes of self""" + + # generate a target tensor + super = [] + sub = [] + for tensor in self.amplitudes: + for nindex in range(len(tensor.indexes)/2): + index = tensor.indexes[nindex] + common = 0 + if (self.summation): + for another in self.summation.indexes: + if (index.isidenticalto(another)): + common = 1 + if (not common): + super.append(tensor.indexes[nindex]) + for nindex in range(len(tensor.indexes)/2,len(tensor.indexes)): + index = tensor.indexes[nindex] + common = 0 + if (self.summation): + for another in self.summation.indexes: + if (index.isidenticalto(another)): + common = 1 + if (not common): + sub.append(tensor.indexes[nindex]) + + # permutation + result = ListOperatorSequences() + result.add(self.duplicate()) + for nsupera in range(len(super)-1): + result = result.targetsuperpermutation(nsupera) + for nsuba in range(len(sub)-1): + result = result.targetsubpermutation(nsuba) + for operatorsequence in result.list: + newsuper = [] + newsub = [] + for tensor in operatorsequence.amplitudes: + for nindex in range(len(tensor.indexes)/2): + index = tensor.indexes[nindex] + common = 0 + if (operatorsequence.summation): + for another in operatorsequence.summation.indexes: + if (index.isidenticalto(another)): + common = 1 + if (not common): + newsuper.append(tensor.indexes[nindex]) + for nindex in range(len(tensor.indexes)/2,len(tensor.indexes)): + index = tensor.indexes[nindex] + common = 0 + if (operatorsequence.summation): + for another in operatorsequence.summation.indexes: + if (index.isidenticalto(another)): + common = 1 + if (not common): + newsub.append(tensor.indexes[nindex]) + for ncoeff in range(len(operatorsequence.factor.coefficients)): + operatorsequence.factor.coefficients[ncoeff] = self.factor.coefficients[ncoeff] + newpermutation = newsuper + newsub + super + sub + if (operatorsequence.factor.permutations[ncoeff] == []): + operatorsequence.factor.permutations[ncoeff] = newpermutation + else: + operatorsequence.factor.permutations[ncoeff] = combinepermutations(newpermutation, operatorsequence.factor.permutations[ncoeff]) + return result + + def canonicalize(self): + """Reorder amplitudes and common indexes in the canonical order""" + # In canonical order, amplitudes are ordered in alphabetical then size-ascending order. + # Then same amplitudes are ordered in ascending order in target index labels. + # Then all common indexes (that are summation indexes) are renamed in the ascending order. + + another = self.duplicate() + + # reorder amplitudes + done = 0 + while (not done): + done = 1 + for namplitudea in range(len(another.amplitudes)): + for namplitudeb in range(len(another.amplitudes)): + if (namplitudea >= namplitudeb): + continue + amplitudea = another.amplitudes[namplitudea] + amplitudeb = another.amplitudes[namplitudeb] + if (amplitudea.isgreaterthan(amplitudeb,another)): + another.swapamplitudes(namplitudea,namplitudeb) + done = 0 + + # reorder indexes + for namplitude in range(len(another.amplitudes)): + amplitude = another.amplitudes[namplitude] + result = amplitude.canonicalize(another) + another.amplitudes[namplitude] = copy.deepcopy(result[0]) + parity = result[1] + another.factor.multiply(parity) + + # relabel summation indexes in the order of appearance + labelsinuse = [] + for amplitude in another.amplitudes: + for operator in amplitude.indexes: + if (not another.summation.hastheindex(operator)): + labelsinuse.append(operator.index) + for sequence in another.sequence: + for operator in sequence: + if (not another.summation.hastheindex(operator)): + labelsinuse.append(operator.index) + oldlabels = [] + newlabels = [] + newlabel = 0 + for amplitude in another.amplitudes: + for operator in amplitude.indexes: + if (another.summation.hastheindex(operator)): + if (operator.index not in oldlabels): + oldlabels.append(operator.index) + newlabel = newlabel + 1 + while (newlabel in labelsinuse): + newlabel = newlabel + 1 + newlabels.append(newlabel) + for operator in another.summation.indexes: + if (operator.index in oldlabels): + operator.index = newlabels[oldlabels.index(operator.index)] + for amplitude in another.amplitudes: + for operator in amplitude.indexes: + if (operator.index in oldlabels): + operator.index = newlabels[oldlabels.index(operator.index)] + for sequence in another.sequence: + for operator in sequence: + if (operator.index in oldlabels): + operator.index = newlabels[oldlabels.index(operator.index)] + + # reorder summation indexes + for nindexa in range(len(another.summation.indexes)): + indexa = another.summation.indexes[nindexa] + for nindexb in range(len(another.summation.indexes)): + indexb = another.summation.indexes[nindexb] + if (nindexa <= nindexb): + continue + if (indexa.index < indexb.index): + swap = another.summation.indexes[nindexa] + another.summation.indexes[nindexa] = copy.deepcopy(another.summation.indexes[nindexb]) + another.summation.indexes[nindexb] = copy.deepcopy(swap) + + return another + + def isacycliccontraction(self): + """Returns 1 if self is a cyclic contraction""" + ncontractions = 0 + for namplitudea in range(len(self.amplitudes)): + for namplitudeb in range(len(self.amplitudes)): + if (namplitudea > namplitudeb): + amplitudea = self.amplitudes[namplitudea] + amplitudeb = self.amplitudes[namplitudeb] + for operator in self.summation.indexes: + if (operator.isin(amplitudea.indexes) and operator.isin(amplitudeb.indexes)): + ncontractions = ncontractions + 1 + break + if (ncontractions > len(self.amplitudes) - 1): + return 1 + else: + return 0 + + def isdisconnected(self,withrespectto=[]): + """Returns 1 if disconnected; if (withrespectto) connectivity among the given amplitude types is tested""" + + if (self.alreadycontracted()): + + # make a connectedness table + connectedness = [0]*len(self.amplitudes) + connectedness[0] = 1 + for iteration in range(len(self.amplitudes)): + for namplitudea in range(len(self.amplitudes)): + if (connectedness[namplitudea] == 1): + amplitudea = self.amplitudes[namplitudea] + if ((withrespectto) and (amplitudea.type not in withrespectto)): + continue + for namplitudeb in range(len(self.amplitudes)): + if (connectedness[namplitudeb] == 0): + amplitudeb = self.amplitudes[namplitudeb] + if ((withrespectto) and (amplitudeb.type not in withrespectto)): + continue + + # see if they have at least one common index + exist = 0 + for indexa in amplitudea.indexes: + for indexb in amplitudeb.indexes: + if (indexa.isidenticalto(indexb)): + exist = 1 + if (exist): + connectedness[namplitudeb] = 1 + + if (0 not in connectedness): + # connected! + return 0 + else: + if (withrespectto): + for namplitudea in range(len(self.amplitudes)): + amplitudea = self.amplitudes[namplitudea] + if (amplitudea.type in withrespectto): + if (connectedness[namplitudea] == 0): + return 1 + return 0 + + else: + return 1 + + else: + return 0 + + def isunlinked(self): + """Returns 1 if unlinked""" + + if (not self.isdisconnected()): + return 0 + + for seed in range(len(self.amplitudes)): + # make a connectedness table + connectedness = [0]*len(self.amplitudes) + connectedness[seed] = 1 + for iteration in range(len(self.amplitudes)): + for namplitudea in range(len(self.amplitudes)): + if (connectedness[namplitudea] == 1): + amplitudea = self.amplitudes[namplitudea] + for namplitudeb in range(len(self.amplitudes)): + if (connectedness[namplitudeb] == 0): + amplitudeb = self.amplitudes[namplitudeb] + + # see if they have at least one common index + exist = 0 + for indexa in amplitudea.indexes: + for indexb in amplitudeb.indexes: + if (indexa.isidenticalto(indexb)): + exist = 1 + if (exist): + connectedness[namplitudeb] = 1 + + if (0 in connectedness): + # disconnected + closed = 1 + for namplitudeb in range(len(self.amplitudes)): + amplitudeb = self.amplitudes[namplitudeb] + if (connectedness[namplitudeb] == 1): + for indexa in amplitudeb.indexes: + if (not self.summation.hastheindex(indexa)): + closed = 0 + if (closed): + # disconnected & closed = unlinked + return 1 + + return 0 + + def iszero(self): + """True if the numerical factor is computationally zero""" + threshold = 1.0e-12 + zero = 1 + for coefficient in self.factor.coefficients: + if (abs(coefficient) > threshold): + zero = 0 + return zero + +class ListOperatorSequences: + + def __init__(self): + """Creates a list of operator sequence objects""" + self.list = [] + + def __str__(self): + """Prints the sequences of operator contractions""" + print("") + for line in self.show(): + print(line) + return "" + + def show(self): + """Returns a human-friendly string of the content""" + show = [] + for operatorsequence in self.list: + if (operatorsequence == "deleted"): + show.append("deleted") + else: + show.append(operatorsequence.show()) + return show + + def tex(self): + """Returns a LaTeX string of the content""" + show = [] + for noperatorsequence in range(len(self.list)): + operatorsequence = self.list[noperatorsequence] + if (operatorsequence != "deleted"): + if (noperatorsequence == 0): + show.append("\\begin{eqnarray}") + show.append(string.join(["&&",operatorsequence.tex(),"\\nonumber\\\\"],"")) + elif (noperatorsequence == len(self.list)-1): + show.append(string.join(["&&",operatorsequence.tex(),"\\nonumber"],"")) + show.append("\\end{eqnarray}") + else: + show.append(string.join(["&&",operatorsequence.tex(),"\\nonumber\\\\"],"")) + return show + + def duplicate(self): + """Makes a copy of itself""" + duplicate = ListOperatorSequences() + for operatorsequence in self.list: + duplicate.list.append(operatorsequence.duplicate()) + return duplicate + + def writetofile(self,filename): + """Writes the output to a given file""" + file = open(filename,"w") + for operatorsequence in self.list: + file.write(operatorsequence.show()) + file.write("\n") + + def appendtofile(self,filename): + """Writes the output to a given file""" + file = open(filename,"a") + for operatorsequence in self.list: + file.write(operatorsequence.show()) + file.write("\n") + + def add(self,newoperatorsequence): + """Adds a new operator sequence member to the list""" + self.list.append(newoperatorsequence) + + def join(self,another): + """Joins two list operator sequences""" + for operatorsequence in another.list: + self.list.append(operatorsequence) + + def performcontraction(self): + """Perform a contraction of the left-most operator""" + + # result will be a list of new operator sequence objects + result = ListOperatorSequences() + + # loop over operator sequences + for operatorsequence in self.list: + + # call performcontraction() + result.join(operatorsequence.performcontraction()) + + return result + + def simplifythreesub(self,verbose=0): + """Simplify the list by consolidating operator sequences using permutation of operators""" + + if (len(self.list) == 1): + return self + + # pick up a pair of operator sequences + for nsequencea in range(len(self.list)): +# if (verbose): +# print('processing ',nsequencea,' / ',range(len(self.list))) + sequencea = self.list[nsequencea] + for nsequenceb in range(len(self.list)): + sequenceb = self.list[nsequenceb] + if (nsequenceb <= nsequencea): + continue + if (sequencea.hasthesameform(sequenceb)): + 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 + del self.list[nsequenceb] + del self.list[nsequencea] + return self + elif (sequencea.hasnomismatch(sequencec)): + permutation = ListOperatorSequences() + permutation.add(sequencec) + # permutation of tensors + for namplitude in range(len(sequencec.amplitudes)): + permutation = permutation.amplitudepermutation(namplitude) + 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 + del self.list[nsequenceb] + del self.list[nsequencea] + return self + # permutation of summation indexes + for noperator in range(len(sequencec.summation.indexes)): + permutation = permutation.operatorpermutation(noperator) + 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 + del self.list[nsequenceb] + del self.list[nsequencea] + return self + + return self + + def simplifyfoursub(self,quick=0): + """Identify the permutation symmetry among the target indexes""" + + if (len(self.list) == 1): + return self + + # pick up a pair of operator sequences + for nsequencea in range(len(self.list)): + sequencea = self.list[nsequencea] + for nsequenceb in range(len(self.list)): + sequenceb = self.list[nsequenceb] + if (nsequenceb <= nsequencea): + continue + if (sequencea.hasthesameform(sequenceb)): + sequencec = sequencea.fullyrelabels(sequenceb) + if (sequencea.hasnomismatch(sequencec)): + permutation = sequencec.targetindexpermutation() + # permutation of target indexes + for sequenced in permutation.list: + if (sequencea.canmerge(sequenced)): + self.add(sequencea.merges(sequenced)) + # Important that the following two statements are executed in this order + del self.list[nsequenceb] + del self.list[nsequencea] + return self + if (quick): + continue + seed = ListOperatorSequences() + seed.add(sequencec) + # permutation of tensors + for namplitude in range(len(sequencec.amplitudes)): + seed = seed.amplitudepermutation(namplitude) + for noperator in range(len(sequencec.summation.indexes)): + seed = seed.operatorpermutation(noperator) + permutation = ListOperatorSequences() + for sequenced in seed.list: + permutation.join(sequenced.targetindexpermutation()) + 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 + del self.list[nsequenceb] + del self.list[nsequencea] + return self + + return self + + def simplify(self,verbose=0): + """Call simplyone through four""" + #self.simplifyone(1) + if (self.containscycliccontractions()): + 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 + # self.simplifyfour(1) + self = copy.deepcopy(self.deletezero()) + return self + + def simplifyone(self,verbose=0): + """Consolidate the identical operator sequences""" + if (len(self.list) == 0): + return self + if (verbose): + print(" ... canonicalizing the expressions") + self = self.canonicalize() + if (len(self.list) == 1): + return self + if (verbose): + print(" ... consolidating terms") + originallength = len(self.list) + # pick up a pair of operator sequences + for nsequencea in range(len(self.list)): + if (verbose): + if ((nsequencea/100)*100 == nsequencea): + print("simplifying:",nsequencea,"/",len(self.list)) + sequencea = self.list[nsequencea] + if (sequencea == "deleted"): + continue + for nsequenceb in range(len(self.list)): + sequenceb = self.list[nsequenceb] + if (nsequenceb <= nsequencea): + continue + if (sequenceb == "deleted"): + continue + if (sequencea.isidenticalto(sequenceb)): + sequencea.factor.add(sequenceb.factor,1) + self.list[nsequencea] = copy.deepcopy(sequencea) + self.list[nsequenceb] = "deleted" + numberofdeleted = self.list.count("deleted") + for dummy in range(numberofdeleted): + self.list.remove("deleted") + if (verbose): + print(" ... %d terms have been consolidated" %(originallength - len(self.list))) + print(" ... number of terms = %d" %(len(self.list))) + return self + + def simplifytwo(self,verbose=0): + """Identify the permutation symmetries of target indexes""" + if (len(self.list) == 0): + return self + self = self.canonicalize() + if (len(self.list) == 1): + return self + print(" ... identifying permutation symmetry among target indexes") + originallength = len(self.list) + # pick up a pair of operator sequences + for nsequencea in range(len(self.list)): + if (verbose): + if ((nsequencea/10)*10 == nsequencea): + print("permutation-simplifying:",nsequencea,"/",len(self.list)) + sequencea = self.list[nsequencea] + if (sequencea == "deleted"): + continue + for nsequenceb in range(len(self.list)): + if (nsequenceb <= nsequencea): + continue + sequenceb = self.list[nsequenceb] + if (sequenceb == "deleted"): + continue + if (sequencea.hasthesameform(sequenceb)): + permutation = sequenceb.targetindexpermutation() + permutation = permutation.canonicalize() + for sequencec in permutation.list: + if (sequencea.isidenticalto(sequencec)): + sequencea.factor.add(sequencec.factor,1) + self.list[nsequencea] = copy.deepcopy(sequencea) + self.list[nsequenceb] = "deleted" + break + numberofdeleted = self.list.count("deleted") + for dummy in range(numberofdeleted): + self.list.remove("deleted") + print(" ... %d terms have been consolidated" %(originallength - len(self.list))) + print(" ... number of terms = %d" %(len(self.list))) + return self + + def simplifythree(self,verbose=0): + """Aggressively consolidate identical terms""" + if (len(self.list) == 0): + return self + elif (len(self.list) == 1): + return self + if (verbose): + print(" ... aggressively consolidating terms") + done = 0 + iteration = 0 + originallength = len(self.list) + while (not done): + iteration = iteration + 1 + if (verbose): + print('iteration ',iteration,' number of terms ',len(self.list)) + beforesimplify = len(self.list) + self = self.simplifythreesub(verbose) + if (len(self.list) < beforesimplify): + done = 0 + else: + done = 1 + print(" ... %d terms have been consolidated" %(originallength - len(self.list))) + print(" ... number of terms = %d" %(len(self.list))) + return self + + def simplifyfour(self,verbose=0): + """Aggressively identify the permutation symmetries of target indexes""" + if (len(self.list) == 0): + return self + elif (len(self.list) == 1): + return self + if (verbose): + print(" ... aggressively identifying permutation symmetry among target indexes") + originallength = len(self.list) + # quick merge to reduce the number of terms + done = 0 + iteration = 0 + while (not done): + iteration = iteration + 1 + beforesimplify = len(self.list) + self = self.simplifyfoursub(1) + if (len(self.list) < beforesimplify): + done = 0 + else: + done = 1 + # more exhaustive merge + done = 0 + iteration = 0 + while (not done): + iteration = iteration + 1 + beforesimplify = len(self.list) + self = self.simplifyfoursub(0) + if (len(self.list) < beforesimplify): + done = 0 + else: + done = 1 + if (originallength - len(self.list) > 0): + print(" ... ***** warning *****") + print(" ... %d terms have been consolidated" %(originallength - len(self.list))) + print(" ... number of terms = %d" %(len(self.list))) + return self + + def performfullcontraction(self): + """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) + self = self.simplifyone() + result = ListOperatorSequences() + + # loop over operator sequences + for operatorsequence in self.list: + + # call performcontraction() + result.join(operatorsequence.performfullcontraction()) + + return result + + def operatorpermutation(self,noperatora=0): + """Return all possible permutation of operators in self""" + + result = ListOperatorSequences() + + for operatorsequence in self.list: + result.add(operatorsequence) + if (operatorsequence.summation.indexes): + operatora = operatorsequence.summation.indexes[noperatora] + for noperatorb in range(len(operatorsequence.summation.indexes)): + if (noperatorb <= noperatora): + continue + operatorb = operatorsequence.summation.indexes[noperatorb] + if (not operatora.issimilarto(operatorb)): + continue + permutation = operatorsequence.duplicate() + permutation.swapoperators(operatora,operatorb) + result.add(permutation) + + return result + + def amplitudepermutation(self,namplitudea=0): + """Return all possible permutation of amplitude in self""" + + result = ListOperatorSequences() + + for operatorsequence in self.list: + result.add(operatorsequence) + amplitudea = operatorsequence.amplitudes[namplitudea] + for namplitudeb in range(len(operatorsequence.amplitudes)): + if (namplitudeb <= namplitudea): + continue + amplitudeb = operatorsequence.amplitudes[namplitudeb] + if (amplitudea.type != amplitudeb.type): + continue + permutation = operatorsequence.duplicate() + permutation.swapamplitudes(namplitudea,namplitudeb) + result.add(permutation) + + return result + + def targetsuperpermutation(self,nsupera=0): + """Return all possible permutation of target super indexes in self""" + + result = ListOperatorSequences() + + for operatorsequence in self.list: + + super = [] + sub = [] + for tensor in operatorsequence.amplitudes: + for nindex in range(len(tensor.indexes)/2): + index = tensor.indexes[nindex] + common = 0 + if (operatorsequence.summation): + for another in operatorsequence.summation.indexes: + if (index.isidenticalto(another)): + common = 1 + if (not common): + super.append(tensor.indexes[nindex]) + for nindex in range(len(tensor.indexes)/2,len(tensor.indexes)): + index = tensor.indexes[nindex] + common = 0 + if (operatorsequence.summation): + for another in operatorsequence.summation.indexes: + if (index.isidenticalto(another)): + common = 1 + if (not common): + sub.append(tensor.indexes[nindex]) + + result.add(operatorsequence) + supera = super[nsupera] + for nsuperb in range(len(super)): + if (nsuperb <= nsupera): + continue + superb = super[nsuperb] + permutation = operatorsequence.duplicate() + permutation.swapoperators(supera,superb) + result.add(permutation) + + return result + + def targetsubpermutation(self,nsuba=0): + """Return all possible permutation of target sub indexes in self""" + + result = ListOperatorSequences() + + for operatorsequence in self.list: + + super = [] + sub = [] + for tensor in operatorsequence.amplitudes: + for nindex in range(len(tensor.indexes)/2): + index = tensor.indexes[nindex] + common = 0 + if (operatorsequence.summation): + for another in operatorsequence.summation.indexes: + if (index.isidenticalto(another)): + common = 1 + if (not common): + super.append(tensor.indexes[nindex]) + for nindex in range(len(tensor.indexes)/2,len(tensor.indexes)): + index = tensor.indexes[nindex] + common = 0 + if (operatorsequence.summation): + for another in operatorsequence.summation.indexes: + if (index.isidenticalto(another)): + common = 1 + if (not common): + sub.append(tensor.indexes[nindex]) + + result.add(operatorsequence) + suba = sub[nsuba] + for nsubb in range(len(sub)): + if (nsubb <= nsuba): + continue + subb = sub[nsubb] + permutation = operatorsequence.duplicate() + permutation.swapoperators(suba,subb) + result.add(permutation) + + return result + + def canonicalize(self): + """Reorder amplitudes and common indexes in the canonical order""" + + for noperatorsequence in range(len(self.list)): + operatorsequence = self.list[noperatorsequence] + self.list[noperatorsequence] = operatorsequence.canonicalize() + return self + + def deletedisconnected(self,withrespectto=[]): + """Deletes disconnected terms""" + + result = ListOperatorSequences() + + originallength = len(self.list) + + # for a fully contracted sequence ... + for noperatorsequence in range(len(self.list)): + operatorsequence = self.list[noperatorsequence] + if (not operatorsequence.isdisconnected(withrespectto)): + result.add(operatorsequence) + + newlength = len(result.list) + + print(" ... %d disconnected terms have been deleted" %(originallength - newlength)) + + return result + + def deleteunlinked(self): + """Deletes unlinked terms""" + + result = ListOperatorSequences() + + originallength = len(self.list) + + # for a fully contracted sequence ... + for noperatorsequence in range(len(self.list)): + operatorsequence = self.list[noperatorsequence] + if (not operatorsequence.isunlinked()): + result.add(operatorsequence) + + newlength = len(result.list) + + print(" ... %d unlinked terms have been deleted" %(originallength - newlength)) + + return result + + def containscycliccontractions(self): + """Returns 1 if self contains a cyclic contraction""" + for operatorsequence in self.list: + if (operatorsequence.isacycliccontraction()): + return 1 + return 0 + + def relabelamplitudes(self,old,new): + """Relabels amplitude""" + for noperatorsequence in range(len(self.list)): + operatorsequence = self.list[noperatorsequence] + for namplitude in range(len(operatorsequence.amplitudes)): + amplitude = operatorsequence.amplitudes[namplitude] + if (amplitude.type == old): + self.list[noperatorsequence].amplitudes[namplitude].type = copy.deepcopy(new) + return self + + def deletezero(self): + """Deletes computationally zero terms""" + + result = ListOperatorSequences() + + originallength = len(self.list) + + # for a fully contracted sequence ... + for noperatorsequence in range(len(self.list)): + operatorsequence = self.list[noperatorsequence] + if (not operatorsequence.iszero()): + result.add(operatorsequence) + + newlength = len(result.list) + + if (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 new file mode 100644 index 0000000..75613a3 --- /dev/null +++ b/src/tce/splitfiles.py @@ -0,0 +1,32 @@ +#!/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") + +nfiles = 0 +for line in source: + if (string.find(line,"SUBROUTINE") != -1): + if (nfiles): + file = open(filename+".F","w") + for newline in filecontent: + file.write(newline) + filename = string.split(line[string.find(line,"SUBROUTINE")+11:],"(")[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") +for newline in filecontent: + file.write(newline) +print("Number of files generated:",nfiles) diff --git a/src/util/errquit.h b/src/util/errquit.h new file mode 100644 index 0000000..9c86292 --- /dev/null +++ b/src/util/errquit.h @@ -0,0 +1,33 @@ +// UERR - Not yet assigned to a category +// UNKNOWN_ERR - Not yet assigned to a category +// MEM_ERROR - 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 +// BASIS_ERR - Error related to basis set +// GEOM_ERR - Error related to geometry +// MA_ERR - local memory error +// GA_ERR - global memory error +// INT_ERR - error related to integrals +// DISK_ERR - error in reading or writing from disk +// CALC_ERR - calculation failed to converge +// FMM_ERR - +// STACK_ERR - error in MA stack +// HEAP_ERR - error in MA heap +const int UERR = 0; +const int UNKNOWN_ERR = 0; +const int MEM_ERR = 10; +const int STACK_ERR = 11; +const int HEAP_ERR = 12; +const int RTDB_ERR = 20; +const int INPUT_ERR = 30; +const int CAPMIS_ERR = 40; +const int BASIS_ERR = 50; +const int GEOM_ERR = 60; +const int GA_ERR = 70; +const int MA_ERR = 80; +const int INT_ERR = 90; +const int DISK_ERR = 100; +const int CALC_ERR = 110; +const int FMM_ERR = 120; +// $Id$