# ____________________________________________________________________________ # # :: R u n E D D B :: 11-Feb-2024 # # This R-script implements the Electron Density of Delocalized Bonds method, # originally introduced in 2014 (DOI:10.1039/c4cp02932a) and revised in 2016 # (10.1016/J.COMPTC.2016.02.003), 2017 (10.1016/J.COMPTC.2016.12.003), 2017 # (10.1039/C7CP06114E), extended in 2019 (10.1002/OPEN.201900014), reviewed # in 2021 (10.1016/B978-0-12-822723-7.00008-X). Project supported by the EU # Framework Programme for Research and Innovation Horizon 2020 under the # Marie Sklodowska-Curie Grant Agreement No 797335 ("MulArEffect") and the # National Agency for Academic Exchange (Poland) under the Bekker programme # Grant Agreement No PPN/BEK/2019/1/00219. # # Website: www.aromaticity.eu Contact: dariusz.szczepanik@uj.edu.pl # ____________________________________________________________________________ # # LIST OF IMPORTANT KEYWORDS AND OPTIONS # ____________________________________________________________________________ # # To check for updates run: RunEDDB -u # To manually select input files run: RunEDDB -i file.fchk file.49 # # To visualize the EDDB(r) function set it TRUE ( -o [filename] ) SAVE_EDDB_TO_GAUSSIAN_FCHK = TRUE # To visualize alpha and beta components of the EDDB(r) function set it TRUE ( -oo [filename] ) SAVE_SPIN_EDDBS_SEPARATELY = FALSE # To dissect EDDB into different components (sigma,pi,delta,etc.) set it TRUE ( -d [list of NOBDs] ) ORBITAL_DISSECTION_OF_EDDB = FALSE # To diagonalize atomic blocks of the EDDB density matrix set it TRUE ( -a ... set TRUE ) ATOMICORBIT_DECOMP_OF_EDDB = FALSE # To skip the EDDB calculations and perform the NPA analysis set it TRUE ( -c ... set TRUE ) SKIP_EDDB__AND_PERFORM_NPA = FALSE # To reduce the full basis into the natural minimum basis (NMB) set it TRUE ( -m ... set TRUE ) MINIMAL_BASIS_CALCULATIONS = FALSE # To perform the EDDB analysis only for formally bonded atoms set it TRUE ( -b ... set TRUE ) CONNECTIVITY_OUT_FROM_FCHK = FALSE # To remove the UHF spin contamination and perform the EDDB analysis for # paired (closed subshell) and unpaired electrons (open subshell) set it TRUE ( -s ... set TRUE ) REMOVE_UHF_SPIN_CONTAMINAT = FALSE # To enable quiet mode set it TRUE ( -q ... set TRUE ) ENABLE_QUIET_TERMINAL_MODE = FALSE # ____________________________________________________________________________ # # TYPES OF THE EDDB FUNCTION # ____________________________________________________________________________ # # -g ............................................. select the EDDB_G function # -h ............................................. select the EDDB_H function # -f [molecuar fragment, e.g. 1:4,7,8] ........... select the EDDB_F function # -e [molecuar fragment, e.g. 1:4,7,8] ........... select the EDDB_E function # -p [delocalization pathway, e.g. 1-2-3-4-5-1] .. select the EDDB_P function # # EDDB_G(r) - electrons delocalized through the system of all chemical bonds # in a molecule including X-H bonds ([G]lobal delocalization). # # EDDB_H(r) - electrons delocalized through the system of all chemical bonds # in a molecule excluding X-H bonds ([H]eavy atoms only). To # study global pi-delocalization it is recommended to use this # function rather than EDDB_G(r) because hydrogen-bond orbitals # tend to conjugate with the adjacent sigma-bonds noticeably inc- # reasing delocalization in sigma system; both funtions give almost # quantitatively the same predictions of the pi-delocalization # if the manual dissection of pi-orbitals is applied (option -d). # # EDDB_F(r) - electrons delocalized through the system of chemical bonds in a # EDDB_E(r) particular molecular [F]ragment. The difference between EDDB_F(r) # and EDDB_E(r) is that in the latter the density of electrons # delocalized in selected fragment is simply cut down from the # EDDB_G(r) function, and thus it includes the [E]xternal (non- # local) resonance effects. # # EDDB_P(r) - electrons delocalized through selected [P]athway. # ____________________________________________________________________________ # # WARNING: MODIFICATION TO THE CODE BELOW IS AT YOUR OWN RISK! # ____________________________________________________________________________ # # ############################################################################ # ##########################Kk0N##############Xk0N############################ # ##########################Xxd0W####WWW#####WOoxX############################ # ############################N0KXXXXXNWWWWWWKON############################## # ###########################WNXKOxxkO00KXXXNNNW############################## # ###########N0XWW#########WXKXXX0o:clllooodkKNWW############NXN############## # ###########XxdKW#######WNKO0XXK0xxxdoooddccd0NWWW##########OcdX############# # ############WN0KNNNNWWNNK0OOOOO0XW#WWWW#WOccoOXNWWWWWW####XO0XW############# # ############WNX00OO0KKXX0kxddO0NW########W0l;:lkOKXXNWWWNNKXW############### # ##########WNXNNKkoloddxxolclOXN###########WOc:,,;ldxxkOOO0XNWW############## # #########N0OKNNKOddxdddoolodkXW############Xxoocclllccol:cxKNWW############# # ########WOoxO000O0N#WWWWWXX0k0WW#######WWWWKxk0KNNXXXXX0l;cd0NWW############ # ##WNNWWN0oldkkO0KN#########XOkk0NWW###WNXOkOXW#########W0l:coOXNW########### # ##WOxOK0dcoOXXXXW##########0xxdx00000O0KKkx0NW###########KocoOKXX000xkXW#### # ###WNN##NxlkXNWWW##########XkO0kdxO00Okdod0NNW##########WKkOKXNNXKXKkxKW#### # ########WKod0NNWWWW######WNKOOO0KN#####NOxxOXW#########W0xxOXNNNW########### # #########WOoxKNNXNW###WNKkxOKNW###########NkkNWWW####WNOolx0XXXW############ # ##########XxdKXX0k0KOOOKd;l0NNW###########NxoK0kOXX0O0Ol:lkKXXW############# # ##########KkOXNX0xOKXX0xl;c0NNW###########NxdKOddk00kxo:;lONWWW############# # ########WXkx0KXKOOXW##WXOkxkKNW###########NOOKK0KXW##Xkl;cxKNWW############# # ###WW###XdoxO0OOKW########W0x0WWWW####WNNX00XW#########Xl;lkXWW############# # ##XOOKK0dldkO0KXW#########WOc:okKXXXXXXKOxdON###########Xo;lkKNWW##WNW###### # ##W00XNXxokKNNNWW#########Nkdoc:cx00OOdooco0NW##########W0llkKXNK0KkokN##### # ########KdxKNNNWWW#########0xO0O0XW#WWKxddxOXW#########WKxx0KNNNNWWX00N##### # ########WOox0NNNNWWWWWWNNNXOONWW#######WNXXkONWW#WWW#WN0ddx0XNNW############ # #########W0oox0KKKXNNNNNX0OKWW############N0kk0XXNNNNNKkoldOKKXW############ # ###########XklclodkO0KXNNXXNNW############0xkxkO0KXXXNXKxooxOXW############# # ############WkoxxxdddodkOKNNWWW#########WKOO0KKOkxxxkkkxdoxKW############### # ##########NKKKN###WNXK0xodOXNWWWWWWWWWWWKkkOKK0OkkkOO0KKKOON################ # ##########XOONW########W0olk0XNNNNNNNNNKxdxO00KNW########WOxXW############## # #########################XxclxOO0KXNNNNXkddxOXW###########Kx0W############## # ##########################WOc:ccoxOO0KK0kdoON############################### # #########################WN0OKKOkxxkOO0000OKW############################### # ########################WOx0N######W######NOONW############################# # #########################X0XW#############W0kXW############################# # ############################################################################ # ____________________________________________________________________________ # # Setting up program parameters and variables # ____________________________________________________________________________ # Program_Version = "11-Feb-2024" if (!exists("SAVE_EDDB_TO_GAUSSIAN_FCHK")) SAVE_EDDB_TO_GAUSSIAN_FCHK = FALSE if (!exists("SAVE_SPIN_EDDBS_SEPARATELY")) SAVE_SPIN_EDDBS_SEPARATELY = FALSE if (!exists("ORBITAL_DISSECTION_OF_EDDB")) ORBITAL_DISSECTION_OF_EDDB = FALSE if (!exists("SKIP_EDDB__AND_PERFORM_NPA")) SKIP_EDDB__AND_PERFORM_NPA = FALSE if (!exists("ATOMICORBIT_DECOMP_OF_EDDB")) ATOMICORBIT_DECOMP_OF_EDDB = FALSE if (!exists("REMOVE_CORE_AND_EMPTY_NAOS")) REMOVE_CORE_AND_EMPTY_NAOS = FALSE if (!exists("MINIMAL_BASIS_CALCULATIONS")) MINIMAL_BASIS_CALCULATIONS = FALSE if (!exists("ENABLE_QUIET_TERMINAL_MODE")) ENABLE_QUIET_TERMINAL_MODE = FALSE if (!exists("REMOVE_UHF_SPIN_CONTAMINAT")) REMOVE_UHF_SPIN_CONTAMINAT = FALSE if (!exists("CONNECTIVITY_OUT_FROM_FCHK")) CONNECTIVITY_OUT_FROM_FCHK = FALSE if (!exists("IGNORE_DSYEVR_LAPACK_ERROR")) IGNORE_DSYEVR_LAPACK_ERROR = TRUE if (!exists("DEFAULT_EDDB_FUNCTION_TYPE")) DEFAULT_EDDB_FUNCTION_TYPE = NA if (!exists("NATURAL_OCCUPATION_THRSHLD")) NATURAL_OCCUPATION_THRSHLD = 0.0000050 if (!exists("NOBD_OCCUPATIONS_SCAL_FACT")) NOBD_OCCUPATIONS_SCAL_FACT = 0.0367493 if ( ENABLE_QUIET_TERMINAL_MODE ) INTERACTIVE_MODE = FALSE else INTERACTIVE_MODE = TRUE DSYEVR_Err_Counter = 0 if (!exists("setMKLthreads")) setMKLthreads = function(mklarg) { invisible() } setMKLthreads(1) STD_IN = file("stdin") PERIODIC_TABLE = c( "H ","He","Li","Be","B ","C ","N ","O ","F ","Ne","Na","Mg","Al","Si", "P ","S ","Cl","Ar","K ","Ca","Sc","Ti","V ","Cr","Mn","Fe","Co","Ni", "Cu","Zn","Ga","Ge","As","Se","Br","Kr","Rb","Sr","Y ","Zr","Nb","Mo", "Tc","Ru","Rh","Pd","Ag","Cd","In","Sn","Sb","Te","I ","Xe","Cs","Ba", "La","Ce","Pr","Nd","Pm","Sm","Eu","Gd","Tb","Dy","Ho","Er","Tm","Yb", "Lu","Hf","Ta","W ","Re","Os","Ir","Pt","Au","Hg","Tl","Pb","Bi","Po", "At","Rn","Fr","Ra","Ac","Th","Pa","U ","Np","Pu","Am","Cm","Bk","Cf", "Es","Fm","Md","No","Lr","Rf","Db","Sg","Bh","Hs","Mt","Ds","Rg","Cn", "Nh","Fl","Mc","Lv","Ts","Og") cat( paste0(" ", paste0(rep("_",74),collapse=""), "\n\n", collapse="") ) cat( paste0(" ", paste0(rep(" ",26),collapse=""),":: R u n E D D B ::", paste0(rep(" ",18),collapse=""),Program_Version,"\n\n",collapse="") ) cat( " This R-script implements the Electron Density of Delocalized Bonds method,\n" ) cat( " originally introduced in 2014 (DOI:10.1039/c4cp02932a) and revised in 2016\n" ) cat( " (10.1016/J.COMPTC.2016.02.003), 2017 (10.1016/J.COMPTC.2016.12.003), 2017\n" ) cat( " (10.1039/C7CP06114E), extended in 2019 (10.1002/OPEN.201900014), reviewed\n" ) cat( " in 2021 (10.1016/B978-0-12-822723-7.00008-X). Project supported by the EU\n" ) cat( " Marie Sklodowska-Curie Grant Agreement No 797335 (\"MulArEffect\") and the\n") cat( " National Agency for Academic Exchange (Poland) under the Bekker programme\n") cat( " Grant Agreement No PPN/BEK/2019/1/00219.\n\n") cat( " Website: www.aromaticity.eu Contact: dariusz.szczepanik@uj.edu.pl" ) cat( paste0("\n ", paste0(rep("_",74),collapse=""), "\n\n", collapse="" ) ) # ____________________________________________________________________________ # # Parsing command line options # ____________________________________________________________________________ # RscriptARGS=commandArgs(trailingOnly = TRUE) ARGS_options = c("-i","-g","-gg","-h","-hh","-f","-ff","-e","-ee","-p","-a","-c","-d","-m","-o","-oo","-u","-q","-b") ARGS_options_idxs = suppressWarnings( which(RscriptARGS %in% ARGS_options) ) if (length(ARGS_options_idxs) > 0) { if (!is.na( match("-i", tolower(RscriptARGS) ) ) ) { if ( !( (match("-i", tolower(RscriptARGS))+1) %in% ARGS_options_idxs ) ) { if (file.exists( basename(RscriptARGS[ (match("-i", tolower(RscriptARGS))+1) ]) )) {ListOfFchkFiles = basename(RscriptARGS[ (match("-i", tolower(RscriptARGS))+1) ]) } else { cat(" o File ", basename(RscriptARGS[ (match("-i", tolower(RscriptARGS))+1) ]) ," not found!\n\n Press to terminate the RunEDDB script.\n") ; if(INTERACTIVE_MODE) nullek=readLines(STD_IN,n=1) else cat("\n") ; close(STD_IN) ; quit(save="no") } } if ( !( (match("-i", tolower(RscriptARGS))+2) %in% ARGS_options_idxs ) ) { if ( file.exists( basename(RscriptARGS[ (match("-i", tolower(RscriptARGS))+2) ]) ) ) { ListOfNBOFiles = basename(RscriptARGS[ (match("-i", tolower(RscriptARGS))+2) ]) } else { cat(" o File ", basename(RscriptARGS[ (match("-i", tolower(RscriptARGS))+2) ]) ," not found!\n\n Press to terminate the RunEDDB script.\n") ; if(INTERACTIVE_MODE) nullek=readLines(STD_IN,n=1) else cat("\n") ; close(STD_IN) ; quit(save="no") } } } if (!is.na( match("-o", tolower(RscriptARGS) ) ) ) { SAVE_EDDB_TO_GAUSSIAN_FCHK = TRUE if ( !( (match("-o", tolower(RscriptARGS))+1) %in% ARGS_options_idxs ) ) { if (!is.na( RscriptARGS[ (match("-o", tolower(RscriptARGS))+1) ] )) SAVE_EDDB_TO_FCHK_FILENAME = RscriptARGS[ (match("-o", tolower(RscriptARGS))+1) ] } } if (!is.na( match("-oo", tolower(RscriptARGS) ) ) ) { SAVE_EDDB_TO_GAUSSIAN_FCHK = TRUE SAVE_SPIN_EDDBS_SEPARATELY = TRUE if ( !( (match("-oo", tolower(RscriptARGS))+1) %in% ARGS_options_idxs ) ) { if (!is.na( RscriptARGS[ (match("-oo", tolower(RscriptARGS))+1) ] )) SAVE_EDDB_TO_FCHK_FILENAME = RscriptARGS[ (match("-oo", tolower(RscriptARGS))+1) ] } } if (!is.na( match("-d", tolower(RscriptARGS) ) ) ) { ORBITAL_DISSECTION_OF_EDDB = TRUE if ( !( (match("-d", tolower(RscriptARGS))+1) %in% ARGS_options_idxs ) ) { if (!is.na( RscriptARGS[ (match("-d", tolower(RscriptARGS))+1) ] )) SELECTED_NOBD_List = RscriptARGS[ (match("-d", tolower(RscriptARGS))+1) ] } if ( !( (match("-d", tolower(RscriptARGS))+2) %in% ARGS_options_idxs ) ) { if (!is.na( RscriptARGS[ (match("-d", tolower(RscriptARGS))+2) ] )) SELECTED_NOBD_Beta_List = RscriptARGS[ (match("-d", tolower(RscriptARGS))+2) ] } } if (!is.na( match("-c", tolower(RscriptARGS) ) ) ) { SKIP_EDDB__AND_PERFORM_NPA = TRUE } if (!is.na( match("-a", tolower(RscriptARGS) ) ) ) { ATOMICORBIT_DECOMP_OF_EDDB = TRUE } if (!is.na( match("-m", tolower(RscriptARGS) ) ) ) { MINIMAL_BASIS_CALCULATIONS = TRUE } if (!is.na( match("-s", tolower(RscriptARGS) ) ) ) { REMOVE_UHF_SPIN_CONTAMINAT = TRUE } if (!is.na( match("-q", tolower(RscriptARGS) ) ) ) { INTERACTIVE_MODE = FALSE } if (!INTERACTIVE_MODE) { txtProgressBar = function(min=0, max=68, initial=0, char="#", width=68, style=3) { invisible() } setTxtProgressBar = function(progressBar,pbinc) { invisible() } } if (!is.na( match("-u", tolower(RscriptARGS) ) ) ) { cat( " o Checking for updates... ") versionFile = suppressWarnings(tempfile()) suppressWarnings(download.file("http://www.aromaticity.eu/RunEDDB.version", versionFile, quiet=TRUE)) currentVersion = suppressWarnings ( scan( text=readChar(versionFile,file.info(versionFile)$size,useBytes=T), what="", sep="", quiet=TRUE ) ) if (!grepl( tolower(currentVersion) , tolower(Program_Version) ,fixed=FALSE) ) { cat(paste0("\n ", paste0(rep("-=",37),collapse=""), collapse="")) cat("\n =-= There is a new version of RunEDDB available on www.aromaticity.eu -=-") cat(paste0("\n ", paste0(rep("-=",37),collapse=""), "", collapse="")) Sys.sleep(5) } else cat(paste0("DONE! This is the latest version of RunEDDB.", collapse="" ) ) cat(paste0("\n o Normal termination of RunEDDB at ",format(Sys.time(), "%a %b %d %X %Y"),".\n", collapse="" ) ) if (INTERACTIVE_MODE) { cat("\n Press to quit the script.\n", collapse="") ; nullek=readLines(STD_IN,n=1) } else { cat("\n\n" ) } close(STD_IN) quit(save="no") } if (!is.na( match("-p", tolower(RscriptARGS) ) ) ) { DEFAULT_EDDB_FUNCTION_TYPE = 5 if ( !( (match("-p", tolower(RscriptARGS))+1) %in% ARGS_options_idxs ) ) { if (!is.na( RscriptARGS[ (match("-p", tolower(RscriptARGS))+1) ] )) SELECTED_Pathway_readline = RscriptARGS[ (match("-p", tolower(RscriptARGS))+1) ] } } if (!is.na( match("-f", tolower(RscriptARGS) ) ) ) { DEFAULT_EDDB_FUNCTION_TYPE = 3 if ( !( (match("-f", tolower(RscriptARGS))+1) %in% ARGS_options_idxs ) ) { if (!is.na( RscriptARGS[ (match("-f", tolower(RscriptARGS))+1) ] )) SELECTED_Fragment_readline = RscriptARGS[ (match("-f", tolower(RscriptARGS))+1) ] } } if (!is.na( match("-ff", tolower(RscriptARGS) ) ) ) { DEFAULT_EDDB_FUNCTION_TYPE = 30 if ( !( (match("-ff", tolower(RscriptARGS))+1) %in% ARGS_options_idxs ) ) { if (!is.na( RscriptARGS[ (match("-ff", tolower(RscriptARGS))+1) ] )) SELECTED_Fragment_readline = RscriptARGS[ (match("-ff", tolower(RscriptARGS))+1) ] } } if (!is.na( match("-e", tolower(RscriptARGS) ) ) ) { DEFAULT_EDDB_FUNCTION_TYPE = 4 if ( !( (match("-e", tolower(RscriptARGS))+1) %in% ARGS_options_idxs ) ) { if (!is.na( RscriptARGS[ (match("-e", tolower(RscriptARGS))+1) ] )) SELECTED_Fragment_readline = RscriptARGS[ (match("-e", tolower(RscriptARGS))+1) ] } } if (!is.na( match("-ee", tolower(RscriptARGS) ) ) ) { DEFAULT_EDDB_FUNCTION_TYPE = 40 if ( !( (match("-ee", tolower(RscriptARGS))+1) %in% ARGS_options_idxs ) ) { if (!is.na( RscriptARGS[ (match("-ee", tolower(RscriptARGS))+1) ] )) SELECTED_Fragment_readline = RscriptARGS[ (match("-ee", tolower(RscriptARGS))+1) ] } } if (!is.na( match("-h", tolower(RscriptARGS) ) ) ) { DEFAULT_EDDB_FUNCTION_TYPE = 2 } if (!is.na( match("-hh", tolower(RscriptARGS) ) ) ) { DEFAULT_EDDB_FUNCTION_TYPE = 20 } if (!is.na( match("-g", tolower(RscriptARGS) ) ) ) { DEFAULT_EDDB_FUNCTION_TYPE = 1 } if (!is.na( match("-gg", tolower(RscriptARGS) ) ) ) { DEFAULT_EDDB_FUNCTION_TYPE = 10 } if (!is.na( match("-b", tolower(RscriptARGS) ) ) ) { CONNECTIVITY_OUT_FROM_FCHK = TRUE } } # ____________________________________________________________________________ # # Reading *.fchk and *.49 input files # ____________________________________________________________________________ # if ( !exists("ListOfFchkFiles") ) ListOfFchkFiles = list.files(pattern = "\\.fchk$", ignore.case=TRUE) if ( length(ListOfFchkFiles) <1 ) { cat(" o Unable to find a formatted checkpoint file (.fchk) in the directory!\n\n Press to terminate the script.\n") if(INTERACTIVE_MODE) nullek=readLines(STD_IN,n=1) else cat("\n") ; close(STD_IN) ; quit(save="no") } else if ( length(ListOfFchkFiles) >1 ) { cat(" o There is more than one formatted checkpoint file in the directory!\n") if (!INTERACTIVE_MODE) { cat("\n Interactive mode disabled! Script terminated.\n\n") ; close(STD_IN) ; quit(save="no") } for (i in 1:length(ListOfFchkFiles)) { cat(sprintf( "%5d) %s\n",i,ListOfFchkFiles[i])) } SELECTED_FileNumber = NA while (is.na(SELECTED_FileNumber) || SELECTED_FileNumber<1 || SELECTED_FileNumber>length(ListOfFchkFiles) ) { cat(sprintf(paste0(" o Select file (1-%", as.integer(log10( length(ListOfFchkFiles) ))+1, "d) and press :\n > ",collapse=""), length(ListOfFchkFiles) ) ) SELECTED_FileNumber = suppressWarnings( as.integer( readLines(STD_IN,n=1) ) ) } ListOfFchkFiles = ListOfFchkFiles[SELECTED_FileNumber] } FCHK_file_name = ListOfFchkFiles[1] # if ( !exists("ListOfNBOFiles") ) ListOfNBOFiles = list.files(pattern = "\\.49$", ignore.case=TRUE) if ( length(ListOfNBOFiles) <1 ) { cat(" o Unable to find the DMNAO.49 file in the directory!\n\n") cat(" Press to terminate the script.\n") if(INTERACTIVE_MODE) nullek=readLines(STD_IN,n=1) else cat("\n") close(STD_IN) quit(save="no") } else if ( length(ListOfNBOFiles) >1 ) { cat(" o There is more than one NBO file with the DMNAO matrix in the directory!\n") if (!INTERACTIVE_MODE) { cat("\n Interactive mode disabled! Script terminated.\n\n") ; close(STD_IN) ; quit(save="no") } for (i in 1:length(ListOfNBOFiles)) { cat(sprintf( "%5d) %s\n",i,ListOfNBOFiles[i])) } SELECTED_FileNumber = NA while (is.na(SELECTED_FileNumber) || SELECTED_FileNumber<1 || SELECTED_FileNumber>length(ListOfNBOFiles) ) { cat(sprintf(paste0(" o Select file (1-%", as.integer(log10( length(ListOfNBOFiles) ))+1, "d) and press \n > ",collapse=""), length(ListOfNBOFiles) ) ) SELECTED_FileNumber = suppressWarnings( as.integer( readLines(STD_IN,n=1) ) ) } ListOfNBOFiles = ListOfNBOFiles[SELECTED_FileNumber] } NBO_file_name = ListOfNBOFiles[1] # ____________________________________________________________________________ # # Parsing formatted checkpoint file .fchk # ____________________________________________________________________________ # cat(paste0(" o Parsing ",FCHK_file_name,"",collapse="")) FCHK_file = file(FCHK_file_name,"rb") FCHK_content = suppressWarnings(scan(text=readChar(FCHK_file,file.info(FCHK_file_name)$size,useBytes=T),what="",sep="\n",quiet=TRUE)) suppressWarnings(close(FCHK_file)) FCHK_content_len = length(FCHK_content) # # Number of atoms # Number_of_atoms = NA fchktxttyp = "Number of atoms I" fchktxttyp = scan(text=fchktxttyp,what="",sep="",quiet=TRUE) ; fchktxttyp_len = length(fchktxttyp) fchktxttyp_found = FALSE ; idx = 1 while( (!fchktxttyp_found) && idx < FCHK_content_len ) { fchkline = scan(text=FCHK_content[idx],what="",sep="",quiet=TRUE) ; fchkline_len = length(fchkline) if ( suppressWarnings( grepl( tolower( paste0( fchktxttyp, collapse="" ) ), tolower( paste0( fchkline[1:fchktxttyp_len], collapse="") ), fixed=FALSE ) ) ) { fchktxttyp_found = TRUE Number_of_atoms = suppressWarnings( as.integer( fchkline[fchkline_len] ) ) } if ( (!fchktxttyp_found) ) { if ( fchkline_len > 3 && suppressWarnings( grepl( "in=", tolower( paste0( fchkline[(fchkline_len-2):fchkline_len], collapse="") ), fixed=FALSE ) ) ) idx = idx + max( 1, suppressWarnings( as.integer( fchkline[fchkline_len] ) %/% 6 ) ) else if ( fchkline_len > 3 && suppressWarnings( grepl( "rn=", tolower( paste0( fchkline[(fchkline_len-2):fchkline_len], collapse="") ), fixed=FALSE ) ) ) idx = idx + max( 1, suppressWarnings( as.integer( fchkline[fchkline_len] ) %/% 5 ) ) else idx = idx + 1 } } if ( is.na( Number_of_atoms ) || ( idx >= FCHK_content_len ) ) { cat(" FAILED!\n Number of atoms not found.\n\n Press to terminate the RunEDDB script.\n") if(INTERACTIVE_MODE) nullek=readLines(STD_IN,n=1) else cat("\n") ; close(STD_IN) ; quit(save="no") } if ( Number_of_atoms < 3 ) { cat(" FAILED!\n Number of atoms must be greater than 2.\n\n Press to terminate the RunEDDB script.\n") if(INTERACTIVE_MODE) nullek=readLines(STD_IN,n=1) else cat("\n") ; close(STD_IN) ; quit(save="no") } # # Number of alpha electrons # Number_of_alpha_electrons = NA fchktxttyp = "Number of alpha electrons I" fchktxttyp = scan(text=fchktxttyp,what="",sep="",quiet=TRUE) ; fchktxttyp_len = length(fchktxttyp) fchktxttyp_found = FALSE ; idx = 1 while( (!fchktxttyp_found) && idx < FCHK_content_len ) { fchkline = scan(text=FCHK_content[idx],what="",sep="",quiet=TRUE) ; fchkline_len = length(fchkline) if ( suppressWarnings( grepl( tolower( paste0( fchktxttyp, collapse="" ) ), tolower( paste0( fchkline[1:fchktxttyp_len], collapse="") ), fixed=FALSE ) ) ) { fchktxttyp_found = TRUE Number_of_alpha_electrons = suppressWarnings( as.integer( fchkline[fchkline_len] ) ) } if ( (!fchktxttyp_found) ) { if ( fchkline_len > 3 && suppressWarnings( grepl( "in=", tolower( paste0( fchkline[(fchkline_len-2):fchkline_len], collapse="") ), fixed=FALSE ) ) ) idx = idx + max( 1, suppressWarnings( as.integer( fchkline[fchkline_len] ) %/% 6 ) ) else if ( fchkline_len > 3 && suppressWarnings( grepl( "rn=", tolower( paste0( fchkline[(fchkline_len-2):fchkline_len], collapse="") ), fixed=FALSE ) ) ) idx = idx + max( 1, suppressWarnings( as.integer( fchkline[fchkline_len] ) %/% 5 ) ) else idx = idx + 1 } } if ( is.na( Number_of_alpha_electrons ) || ( idx >= FCHK_content_len ) ) { cat(" FAILED!\n Number of alpha electrons not found.\n\n Press to terminate the RunEDDB script.\n") if(INTERACTIVE_MODE) nullek=readLines(STD_IN,n=1) else cat("\n") ; close(STD_IN) ; quit(save="no") } # # Number of beta electrons # Number_of_beta_electrons = NA fchktxttyp = "Number of beta electrons I" fchktxttyp = scan(text=fchktxttyp,what="",sep="",quiet=TRUE) ; fchktxttyp_len = length(fchktxttyp) fchktxttyp_found = FALSE ; idx = 1 while( (!fchktxttyp_found) && idx < FCHK_content_len ) { fchkline = scan(text=FCHK_content[idx],what="",sep="",quiet=TRUE) ; fchkline_len = length(fchkline) if ( suppressWarnings( grepl( tolower( paste0( fchktxttyp, collapse="" ) ), tolower( paste0( fchkline[1:fchktxttyp_len], collapse="") ), fixed=FALSE ) ) ) { fchktxttyp_found = TRUE Number_of_beta_electrons = suppressWarnings( as.integer( fchkline[fchkline_len] ) ) } if ( (!fchktxttyp_found) ) { if ( fchkline_len > 3 && suppressWarnings( grepl( "in=", tolower( paste0( fchkline[(fchkline_len-2):fchkline_len], collapse="") ), fixed=FALSE ) ) ) idx = idx + max( 1, suppressWarnings( as.integer( fchkline[fchkline_len] ) %/% 6 ) ) else if ( fchkline_len > 3 && suppressWarnings( grepl( "rn=", tolower( paste0( fchkline[(fchkline_len-2):fchkline_len], collapse="") ), fixed=FALSE ) ) ) idx = idx + max( 1, suppressWarnings( as.integer( fchkline[fchkline_len] ) %/% 5 ) ) else idx = idx + 1 } } if ( is.na( Number_of_beta_electrons ) || ( idx >= FCHK_content_len ) ) { cat(" FAILED!\n Number of beta electrons not found.\n\n Press to terminate the RunEDDB script.\n") if(INTERACTIVE_MODE) nullek=readLines(STD_IN,n=1) else cat("\n") ; close(STD_IN) ; quit(save="no") } # Number_of_electrons = Number_of_alpha_electrons + Number_of_beta_electrons # # Number of basis functions # Number_of_basis_functions = NA fchktxttyp = "Number of basis functions I" fchktxttyp = scan(text=fchktxttyp,what="",sep="",quiet=TRUE) ; fchktxttyp_len = length(fchktxttyp) fchktxttyp_found = FALSE ; idx = 1 while( (!fchktxttyp_found) && idx < FCHK_content_len ) { fchkline = scan(text=FCHK_content[idx],what="",sep="",quiet=TRUE) ; fchkline_len = length(fchkline) if ( suppressWarnings( grepl( tolower( paste0( fchktxttyp, collapse="" ) ), tolower( paste0( fchkline[1:fchktxttyp_len], collapse="") ), fixed=FALSE ) ) ) { fchktxttyp_found = TRUE Number_of_basis_functions = suppressWarnings( as.integer( fchkline[fchkline_len] ) ) } if ( (!fchktxttyp_found) ) { if ( fchkline_len > 3 && suppressWarnings( grepl( "in=", tolower( paste0( fchkline[(fchkline_len-2):fchkline_len], collapse="") ), fixed=FALSE ) ) ) idx = idx + max( 1, suppressWarnings( as.integer( fchkline[fchkline_len] ) %/% 6 ) ) else if ( fchkline_len > 3 && suppressWarnings( grepl( "rn=", tolower( paste0( fchkline[(fchkline_len-2):fchkline_len], collapse="") ), fixed=FALSE ) ) ) idx = idx + max( 1, suppressWarnings( as.integer( fchkline[fchkline_len] ) %/% 5 ) ) else idx = idx + 1 } } if ( is.na( Number_of_basis_functions ) || ( idx >= FCHK_content_len ) ) { cat(" FAILED!\n Number of basis functions not found.\n\n Press to terminate the RunEDDB script.\n") if(INTERACTIVE_MODE) nullek=readLines(STD_IN,n=1) else cat("\n") ; close(STD_IN) ; quit(save="no") } # # Number of independent functions # Number_of_independent_functions = NA fchktxttyp = "Number of independent functions I" fchktxttyp = scan(text=fchktxttyp,what="",sep="",quiet=TRUE) ; fchktxttyp_len = length(fchktxttyp) fchktxttyp_found = FALSE ; idx = 1 while( (!fchktxttyp_found) && idx < FCHK_content_len ) { fchkline = scan(text=FCHK_content[idx],what="",sep="",quiet=TRUE) ; fchkline_len = length(fchkline) if ( suppressWarnings( grepl( tolower( paste0( fchktxttyp, collapse="" ) ), tolower( paste0( fchkline[1:fchktxttyp_len], collapse="") ), fixed=FALSE ) ) ) { fchktxttyp_found = TRUE Number_of_independent_functions = suppressWarnings( as.integer( fchkline[fchkline_len] ) ) } if ( (!fchktxttyp_found) ) { if ( fchkline_len > 3 && suppressWarnings( grepl( "in=", tolower( paste0( fchkline[(fchkline_len-2):fchkline_len], collapse="") ), fixed=FALSE ) ) ) idx = idx + max( 1, suppressWarnings( as.integer( fchkline[fchkline_len] ) %/% 6 ) ) else if ( fchkline_len > 3 && suppressWarnings( grepl( "rn=", tolower( paste0( fchkline[(fchkline_len-2):fchkline_len], collapse="") ), fixed=FALSE ) ) ) idx = idx + max( 1, suppressWarnings( as.integer( fchkline[fchkline_len] ) %/% 5 ) ) else idx = idx + 1 } } if ( is.na( Number_of_independent_functions ) || ( idx >= FCHK_content_len ) ) { cat(" FAILED!\n Number of independent functions not found.\n\n Press to terminate the RunEDDB script.\n") if(INTERACTIVE_MODE) nullek=readLines(STD_IN,n=1) else cat("\n") ; close(STD_IN) ; quit(save="no") } # # Array of atomic numbers # Array_of_atomic_numbers = NA fchktxttyp = "Atomic numbers I N=" fchktxttyp = scan(text=fchktxttyp,what="",sep="",quiet=TRUE) ; fchktxttyp_len = length(fchktxttyp) fchktxttyp_found = FALSE ; idx = 1 while( (!fchktxttyp_found) && idx < FCHK_content_len ) { fchkline = scan(text=FCHK_content[idx],what="",sep="",quiet=TRUE) ; fchkline_len = length(fchkline) if ( suppressWarnings( grepl( tolower( paste0( fchktxttyp, collapse="" ) ), tolower( paste0( fchkline[1:fchktxttyp_len], collapse="") ), fixed=FALSE ) ) ) { fchktxttyp_found = TRUE Array_of_atomic_numbers = suppressWarnings( as.integer( scan( text=FCHK_content[( idx + 1 ):( idx + ceiling(Number_of_atoms/6 ) )] , what="" , sep="" , quiet=TRUE ) ) ) } if ( (!fchktxttyp_found) ) { if ( fchkline_len > 3 && suppressWarnings( grepl( "in=", tolower( paste0( fchkline[(fchkline_len-2):fchkline_len], collapse="") ), fixed=FALSE ) ) ) idx = idx + max( 1, suppressWarnings( as.integer( fchkline[fchkline_len] ) %/% 6 ) ) else if ( fchkline_len > 3 && suppressWarnings( grepl( "rn=", tolower( paste0( fchkline[(fchkline_len-2):fchkline_len], collapse="") ), fixed=FALSE ) ) ) idx = idx + max( 1, suppressWarnings( as.integer( fchkline[fchkline_len] ) %/% 5 ) ) else idx = idx + 1 } } if ( is.na( Array_of_atomic_numbers[1] ) || ( idx >= FCHK_content_len ) ) { cat(" FAILED!\n Array of atomic numbers not found.\n\n Press to terminate the RunEDDB script.\n") if(INTERACTIVE_MODE) nullek=readLines(STD_IN,n=1) else cat("\n") ; close(STD_IN) ; quit(save="no") } # # Matrix of formal connectivity # if (CONNECTIVITY_OUT_FROM_FCHK) { Matrix_of_formal_connectivity = NA fchktxttyp = "MxBond I" fchktxttyp = scan(text=fchktxttyp,what="",sep="",quiet=TRUE) ; fchktxttyp_len = length(fchktxttyp) fchktxttyp_found = FALSE ; idx = 1 while( (!fchktxttyp_found) && idx < FCHK_content_len ) { fchkline = scan(text=FCHK_content[idx],what="",sep="",quiet=TRUE) ; fchkline_len = length(fchkline) if ( suppressWarnings( grepl( tolower( paste0( fchktxttyp, collapse="" ) ), tolower( paste0( fchkline[1:fchktxttyp_len], collapse="") ), fixed=FALSE ) ) ) { fchktxttyp_found = TRUE MxBond = as.integer(scan( text=FCHK_content[idx] , what="" , sep="" , quiet=TRUE )[3]) idx = idx + 2 Matrix_of_formal_connectivity = matrix ( 0 , Number_of_atoms, Number_of_atoms ) diag(Matrix_of_formal_connectivity) = as.integer(scan( text=FCHK_content[idx:(idx + ceiling(Number_of_atoms/6) - 1)] , what="" , sep="" , quiet=TRUE )) idx = idx + ceiling(Number_of_atoms/6) + 1 IBond = as.integer(scan( text=FCHK_content[idx:(idx + ceiling(Number_of_atoms*MxBond/6) - 1)] , what="" , sep="" , quiet=TRUE )) for (i in 1:length(IBond)) if (IBond[i]>0) { Matrix_of_formal_connectivity[ceiling(i/MxBond),IBond[i]] = 1 Matrix_of_formal_connectivity[IBond[i],ceiling(i/MxBond)] = 1 } } idx = idx + 1 } if ( is.na( Matrix_of_formal_connectivity ) || ( idx >= FCHK_content_len ) ) { cat(" FAILED!\n Atom-atom connectivity matrix not found.\n\n Press to terminate the RunEDDB script.\n") if(INTERACTIVE_MODE) nullek=readLines(STD_IN,n=1) else cat("\n") ; close(STD_IN) ; quit(save="no") } Matrix_of_formal_connectivity = Matrix_of_formal_connectivity > 0 } # ____________________________________________________________________________ # # Parsing NBO file .49 # ____________________________________________________________________________ # NBODensType = 1 cat(paste0(" and ",NBO_file_name,"...",collapse="")) NBO_file = file(NBO_file_name,"rb") NBO_content = suppressWarnings(scan(text=readChar(NBO_file,file.info(NBO_file_name)$size,useBytes=T),skip=0,nlines=4,what="",quiet=TRUE)) NBO_content_len = length(NBO_content) # # NAOs in the AO basis # if (suppressWarnings(is.na(sum(as.numeric(NBO_content[(NBO_content_len-4):NBO_content_len]))))) { cat(" FAILED!\n The AONAO matrix not found.\n\n Press to terminate the RunEDDB script.\n") ; if(INTERACTIVE_MODE) nullek=readLines(STD_IN,n=1) else cat("\n") ; close(STD_IN) ; quit(save="no") } if (suppressWarnings(!is.na(sum(as.numeric(NBO_content[1:(NBO_content_len-5)]))))) { NumberOfElementsToSkip = 0 } else { NumberOfElementsToSkip = NBO_content_len - 5 if (!grepl("NAOsintheAObasis:",paste0(NBO_content[NBO_content_len-10],NBO_content[NBO_content_len-9],NBO_content[NBO_content_len-8],NBO_content[NBO_content_len-7],NBO_content[NBO_content_len-6],collapse=""),fixed=FALSE)) { cat(" FAILED!\n The AONAO matrix not found.\n\n Press to terminate the RunEDDB script.\n") ; if(INTERACTIVE_MODE) nullek=readLines(STD_IN,n=1) else cat("\n") ; close(STD_IN) ; quit(save="no") } } close(NBO_file) ; NBO_file = file(NBO_file_name,"rb") NBO_content = suppressWarnings(scan(text=readChar(NBO_file,file.info(NBO_file_name)$size,useBytes=T),what="",quiet=TRUE)) close(NBO_file) NBO_content_len = length(NBO_content) FirstElementIndex = NumberOfElementsToSkip + 1 idx = FirstElementIndex while (grepl("\\.",NBO_content[idx],fixed=FALSE)) idx = idx + Number_of_basis_functions if (suppressWarnings(is.na(NBO_content[idx]))) { cat(" FAILED!\n The AONAO matrix not found.\n\n Press to terminate the RunEDDB script.\n") ; if(INTERACTIVE_MODE) nullek=readLines(STD_IN,n=1) else cat("\n") ; close(STD_IN) ; quit(save="no") } while (!grepl("\\.",NBO_content[idx],fixed=FALSE)) idx = idx - 1 if (((idx - FirstElementIndex + 1) %% Number_of_basis_functions) > 0 ) { cat("FAILED!\n Cannot read the AONAO matrix (columns to close).\n\n Press to terminate the RunEDDB script.\n") ; if(INTERACTIVE_MODE) nullek=readLines(STD_IN,n=1) else cat("\n") ; close(STD_IN) ; quit(save="no") } Number_of_NAOs = (idx - FirstElementIndex + 1) %/% Number_of_basis_functions AOtoNAO_matrix = matrix(suppressWarnings(as.numeric(NBO_content[FirstElementIndex:idx])),Number_of_basis_functions,Number_of_NAOs) idx = idx + 1 NAOtoAtomArray = matrix(0,Number_of_atoms,2) for (i in idx:(idx + Number_of_NAOs -1)) NAOtoAtomArray[as.integer(NBO_content[i]),2] = i - idx + 1 for (i in (idx + Number_of_NAOs -1):idx) NAOtoAtomArray[as.integer(NBO_content[i]),1] = i - idx + 1 if (suppressWarnings(is.na(sum(NAOtoAtomArray)))) { cat(" FAILED!\n The AONAO matrix not found.\n\n Press to terminate the RunEDDB script.\n") ; if(INTERACTIVE_MODE) nullek=readLines(STD_IN,n=1) else cat("\n") ; close(STD_IN) ; quit(save="no") } idx = idx + 2*Number_of_NAOs isNAOaNMB = suppressWarnings(as.logical(as.integer(NBO_content[idx:(idx + Number_of_NAOs - 1)]))) idx = idx + Number_of_NAOs # # NAO density matrix # while (!(grepl("\\.",NBO_content[idx],fixed=FALSE) && grepl("\\.",NBO_content[idx+1],fixed=FALSE) && grepl("\\.",NBO_content[idx+2],fixed=FALSE) && grepl("\\.",NBO_content[idx+3],fixed=FALSE) && grepl("\\.",NBO_content[idx+4],fixed=FALSE))) idx = idx + 1 idx = idx + ((Number_of_NAOs*Number_of_NAOs)-Number_of_NAOs)/2 + Number_of_NAOs while (suppressWarnings(idx to terminate the RunEDDB script.\n") ; if(INTERACTIVE_MODE) nullek=readLines(STD_IN,n=1) else cat("\n") ; close(STD_IN) ; quit(save="no") } idx = idx + 4 isSpinUnrestricted = FALSE if ( grepl( "alphaspin" , tolower( paste0( NBO_content[idx:(idx+1)], collapse="" ) ), fixed=FALSE ) ) { isSpinUnrestricted = TRUE idx = idx + 2 } Matrix_of_Alpha_Density = matrix(NA,Number_of_NAOs,Number_of_NAOs) Matrix_of_Alpha_Density[upper.tri(Matrix_of_Alpha_Density,diag=TRUE)] = suppressWarnings(as.numeric(NBO_content[idx:(idx+((Number_of_NAOs*Number_of_NAOs)-Number_of_NAOs)/2+Number_of_NAOs-1)])) Matrix_of_Alpha_Density = as.matrix(Matrix::forceSymmetric(Matrix_of_Alpha_Density,uplo="U")) if (is.na(sum(diag(Matrix_of_Alpha_Density)))) { cat(" FAILED!\n The DMNAO matrix not found.\n\n Press to terminate the RunEDDB script.\n") ; if(INTERACTIVE_MODE) nullek=readLines(STD_IN,n=1) else cat("\n") ; close(STD_IN) ; quit(save="no") } if (isSpinUnrestricted) { idx = idx + ((Number_of_NAOs*Number_of_NAOs)-Number_of_NAOs)/2 + Number_of_NAOs + 2 Matrix_of_Beta_Density = matrix(NA,Number_of_NAOs,Number_of_NAOs) Matrix_of_Beta_Density[upper.tri(Matrix_of_Beta_Density,diag=TRUE)] = suppressWarnings(as.numeric(NBO_content[idx:(idx+((Number_of_NAOs*Number_of_NAOs)-Number_of_NAOs)/2+Number_of_NAOs-1)])) Matrix_of_Beta_Density = as.matrix(Matrix::forceSymmetric(Matrix_of_Beta_Density,uplo="U")) if (is.na(sum(diag(Matrix_of_Beta_Density)))) { cat(" FAILED!\n The DMNAO matrix not found.\n\n Press to terminate the RunEDDB script.\n") ; if(INTERACTIVE_MODE) nullek=readLines(STD_IN,n=1) else cat("\n") ; close(STD_IN) ; quit(save="no") } } cat(" DONE!\n") # ____________________________________________________________________________ # # Dealing with post-HF (correlated) wavefunctions # ____________________________________________________________________________ # isPostHF = FALSE if (!isSpinUnrestricted) { Natural_Orbitals = eigen(Matrix_of_Alpha_Density,symmetric=TRUE) Natural_Occupations = Natural_Orbitals$values Natural_Orbitals = Natural_Orbitals$vectors QDmat = matrix(1,1,2) for (i in 2:Number_of_NAOs) { if ((Natural_Occupations[i-1]-Natural_Occupations[i])>=NATURAL_OCCUPATION_THRSHLD) { QDmat[nrow(QDmat),2] = i-1 QDmat = rbind( QDmat, c(i,Number_of_NAOs) ) } } QDval = rep(0,nrow(QDmat)) for (i in 1:nrow(QDmat)) { QDval[i] = mean(Natural_Occupations[QDmat[i,1]:QDmat[i,2]]) } QDval = QDval[1:sum(as.integer(QDval>=NATURAL_OCCUPATION_THRSHLD))] QDnos = length(QDval) if (QDnos<2) QDmat = matrix(c(QDmat[1,1],QDmat[1,2]),1,2) else QDmat[1:QDnos,] if (QDnos>1) isPostHF = TRUE } else { Natural_Orbitals_Alpha = eigen(Matrix_of_Alpha_Density,symmetric=TRUE) Natural_Occupations_Alpha = Natural_Orbitals_Alpha$values Natural_Orbitals_Alpha = Natural_Orbitals_Alpha$vectors Natural_Orbitals_Beta = eigen(Matrix_of_Beta_Density,symmetric=TRUE) Natural_Occupations_Beta = Natural_Orbitals_Beta$values Natural_Orbitals_Beta = Natural_Orbitals_Beta$vectors QDmat_Alpha = matrix(1,1,2) QDmat_Beta = matrix(1,1,2) for (i in 2:Number_of_NAOs) { if ((Natural_Occupations_Alpha[i-1]-Natural_Occupations_Alpha[i])>=NATURAL_OCCUPATION_THRSHLD) { QDmat_Alpha[nrow(QDmat_Alpha),2] = i-1 QDmat_Alpha = rbind( QDmat_Alpha, c(i,Number_of_NAOs) ) } if ((Natural_Occupations_Beta[i-1]-Natural_Occupations_Beta[i])>=NATURAL_OCCUPATION_THRSHLD) { QDmat_Beta[nrow(QDmat_Beta),2] = i-1 QDmat_Beta = rbind( QDmat_Beta, c(i,Number_of_NAOs) ) } } QDval_Alpha = rep(0,nrow(QDmat_Alpha)) QDval_Beta = rep(0,nrow(QDmat_Beta)) for (i in 1:nrow(QDmat_Alpha)) { QDval_Alpha[i] = mean(Natural_Occupations_Alpha[QDmat_Alpha[i,1]:QDmat_Alpha[i,2]]) } for (i in 1:nrow(QDmat_Beta)) { QDval_Beta[i] = mean(Natural_Occupations_Beta[QDmat_Beta[i,1]:QDmat_Beta[i,2]]) } QDval_Alpha = QDval_Alpha[1:sum(as.integer( QDval_Alpha >= NATURAL_OCCUPATION_THRSHLD ))] QDval_Beta = QDval_Beta[1:sum(as.integer( QDval_Beta >= NATURAL_OCCUPATION_THRSHLD ))] QDnos_Alpha = length(QDval_Alpha) QDnos_Beta = length(QDval_Beta) QDnos = max(QDnos_Alpha , QDnos_Beta) if (QDnos_Alpha<2) QDmat_Alpha = matrix(c(QDmat_Alpha[1,1],QDmat_Alpha[1,2]),1,2) else QDmat_Alpha[1:QDnos_Alpha,] if (QDnos_Beta<2) QDmat_Beta = matrix(c(QDmat_Beta[1,1],QDmat_Beta[1,2]),1,2) else QDmat_Beta[1:QDnos_Beta,] if (QDnos>1) isPostHF = TRUE } if (!isPostHF) { if (isSpinUnrestricted) { suppressWarnings( rm( Natural_Occupations, Natural_Orbitals, QDmat, QDval ) ) } else { suppressWarnings( rm ( Natural_Occupations_Alpha, Natural_Occupations_Beta, Natural_Orbitals_Alpha, Natural_Orbitals_Beta, QDmat_Alpha, QDmat_Beta, QDval_Alpha, QDval_Beta ) ) } } else { cat(" o 1eRDM is of post-HF type (fractional NO occupation numbers).\n") } # ____________________________________________________________________________ # # Removing spin contamination from the spinless 1eRDM # ____________________________________________________________________________ # if (isSpinUnrestricted && REMOVE_UHF_SPIN_CONTAMINAT && ( isPostHF || ( NBODensType != 1 ) ) ) { REMOVE_UHF_SPIN_CONTAMINAT = FALSE cat(" o Removing spin contamination aborted (HF/DFT-type DMNAO only)!\n") } if (isSpinUnrestricted && REMOVE_UHF_SPIN_CONTAMINAT) { cat(" o Removing spin contamination from 1eRDM...") RDM_eigenvectors = eigen(Matrix_of_Alpha_Density + Matrix_of_Beta_Density, symmetric=TRUE) RDM_eigenvalues = RDM_eigenvectors$values if ( length( RDM_eigenvalues[ abs(RDM_eigenvalues) >= NATURAL_OCCUPATION_THRSHLD ] ) != ( Number_of_alpha_electrons - Number_of_beta_electrons ) ) { RDM_eigenvectors = eigen(Matrix_of_Alpha_Density - Matrix_of_Beta_Density,symmetric=TRUE) RDM_eigenvalues = diag(RDM_eigenvectors$values, Number_of_NAOs) NUnPEl = rep(TRUE, Number_of_NAOs) for (i in 1:Number_of_NAOs) if ( ( 1-abs(RDM_eigenvalues[i,i]) ) < NATURAL_OCCUPATION_THRSHLD ) { NUnPEl[i] = FALSE ; RDM_eigenvalues[i,i] = 0 } RDM_eigenvectors = eigen( ( RDM_eigenvectors$vectors %*% RDM_eigenvalues %*% t(RDM_eigenvectors$vectors) ) + Matrix_of_Alpha_Density + Matrix_of_Beta_Density, symmetric=TRUE ) RDM_eigenvalues = RDM_eigenvectors$values } RDM_eigenvectors = RDM_eigenvectors$vectors Number_of_alpha_electrons = sum ( as.integer( round(RDM_eigenvalues,digits=0) > 0 ) ) Number_of_beta_electrons = sum ( as.integer( round(RDM_eigenvalues,digits=0) == 2 ) ) if ( (Number_of_alpha_electrons + Number_of_beta_electrons) != Number_of_electrons ) { cat(" FAILED!\n\n Press to terminate the RunEDDB script.\n") ; if(INTERACTIVE_MODE) nullek=readLines(STD_IN,n=1) else cat("\n") ; close(STD_IN) ; quit(save="no") } Matrix_of_Alpha_Density = RDM_eigenvectors[,1:Number_of_alpha_electrons] %*% t( RDM_eigenvectors[,1:Number_of_alpha_electrons] ) Matrix_of_Beta_Density = RDM_eigenvectors[,1:Number_of_beta_electrons] %*% t( RDM_eigenvectors[,1:Number_of_beta_electrons] ) cat(" DONE!\n") suppressWarnings( rm( RDM_eigenvectors, RDM_eigenvalues, NUnPEl ) ) } if (!isSpinUnrestricted) DMtrace = sum(diag(Matrix_of_Alpha_Density)) else DMtrace = sum(diag(Matrix_of_Alpha_Density + Matrix_of_Beta_Density)) ftrc = abs( 1 - DMtrace / (Number_of_electrons) ) if ( (ftrc >= NATURAL_OCCUPATION_THRSHLD) && ( NBODensType == 1 ) ) { cat(sprintf(paste0(" The density matrix contains %", as.integer(log10( DMtrace ))+6, ".4fe but expected %" , as.integer(log10( (Number_of_electrons) ))+6 , ".4fe.",collapse=""), DMtrace , Number_of_electrons ) ) cat("\n\n Press to terminate the RunEDDB script.\n") if(INTERACTIVE_MODE) nullek=readLines(STD_IN,n=1) else cat("\n") ; close(STD_IN) ; quit(save="no") } Number_of_electrons_in_DM = DMtrace # ____________________________________________________________________________ # # Reduction of the number of NAO (NAO->NMB) # ____________________________________________________________________________ # DMtraceNMB = 0 if (isTRUE(MINIMAL_BASIS_CALCULATIONS)) { cat(sprintf(paste0(" o Switching from NAO (%", as.integer(log10( DMtrace ))+6 ,".4fe) to NMB (", collapse=""), DMtrace) ) VectorOfAtomLabels = rep(0,Number_of_NAOs) for (i in 1:Number_of_atoms) VectorOfAtomLabels[NAOtoAtomArray[i,1]:NAOtoAtomArray[i,2]] = i VectorOfAtomLabels = VectorOfAtomLabels[which(isNAOaNMB)] Number_of_NAOs = length(VectorOfAtomLabels) NAOtoAtomArray = matrix(0,Number_of_atoms,2) for (i in 1:Number_of_NAOs) NAOtoAtomArray[VectorOfAtomLabels[i],2] = i for (i in Number_of_NAOs:1) NAOtoAtomArray[VectorOfAtomLabels[i],1] = i Matrix_of_Alpha_Density = Matrix_of_Alpha_Density[,which(isNAOaNMB)] Matrix_of_Alpha_Density = Matrix_of_Alpha_Density[which(isNAOaNMB),] AOtoNAO_matrix = AOtoNAO_matrix[, which(isNAOaNMB)] DMtraceNMB = sum(diag( Matrix_of_Alpha_Density )) if (isSpinUnrestricted) { Matrix_of_Beta_Density = Matrix_of_Beta_Density[,which(isNAOaNMB)] Matrix_of_Beta_Density = Matrix_of_Beta_Density[which(isNAOaNMB),] DMtraceNMB = DMtraceNMB + sum(diag( Matrix_of_Beta_Density )) } if (isPostHF) { if(!isSpinUnrestricted) { Natural_Orbitals = Natural_Orbitals[which(isNAOaNMB),] } else { Natural_Orbitals_Alpha = Natural_Orbitals_Alpha[which(isNAOaNMB),] Natural_Orbitals_Beta = Natural_Orbitals_Beta[which(isNAOaNMB),] } } cat(sprintf(paste0("%", as.integer(log10( DMtraceNMB ))+6, ".4fe, %" , as.integer(log10( DMtraceNMB/DMtrace ))+3, ".1f%%) basis... DONE!\n", collapse=""), DMtraceNMB, 100*DMtraceNMB/DMtrace ) ) EDDB_basis = "NMB" Number_of_electrons_in_DM = DMtraceNMB } else EDDB_basis = "NAO" # DMtraceNVB = 0 if (isTRUE(REMOVE_CORE_AND_EMPTY_NAOS)) # Must be reimplemented <===================================== { isNMBaNVB = rep( TRUE, Number_of_NAOs ) for (i in 1:Number_of_NAOs) { OrbOcc = Matrix_of_Alpha_Density[i,i] if (isSpinUnrestricted) OrbOcc = OrbOcc + Matrix_of_Beta_Density[i,i] if ( ( ( 2 - OrbOcc ) < 10*NATURAL_OCCUPATION_THRSHLD ) || ( OrbOcc < 10*NATURAL_OCCUPATION_THRSHLD ) ) isNMBaNVB[i] = FALSE } NumCorNMB = sum( as.integer( !isNMBaNVB ) ) VectorOfAtomLabels = rep(0,Number_of_NAOs) for (i in 1:Number_of_atoms) VectorOfAtomLabels[NAOtoAtomArray[i,1]:NAOtoAtomArray[i,2]] = i VectorOfAtomLabels = VectorOfAtomLabels[which(isNMBaNVB)] Number_of_NAOs = length(VectorOfAtomLabels) NAOtoAtomArray = matrix(0,Number_of_atoms,2) for (i in 1:Number_of_NAOs) NAOtoAtomArray[VectorOfAtomLabels[i],2] = i for (i in Number_of_NAOs:1) NAOtoAtomArray[VectorOfAtomLabels[i],1] = i Matrix_of_Alpha_Density = Matrix_of_Alpha_Density[,which(isNMBaNVB)] Matrix_of_Alpha_Density = Matrix_of_Alpha_Density[which(isNMBaNVB),] AOtoNAO_matrix = AOtoNAO_matrix[, which(isNMBaNVB)] DMtraceNVB = sum(diag( Matrix_of_Alpha_Density )) if (isSpinUnrestricted) { Matrix_of_Beta_Density = Matrix_of_Beta_Density[,which(isNMBaNVB)] Matrix_of_Beta_Density = Matrix_of_Beta_Density[which(isNMBaNVB),] DMtraceNVB = DMtraceNVB + sum(diag( Matrix_of_Beta_Density )) } if (isPostHF) { if(!isSpinUnrestricted) { Natural_Orbitals = Natural_Orbitals[which(isNMBaNVB),] } else { Natural_Orbitals_Alpha = Natural_Orbitals_Alpha[which(isNMBaNVB),] Natural_Orbitals_Beta = Natural_Orbitals_Beta[which(isNMBaNVB),] } } if (NumCorNMB>0) { cat(sprintf(paste0(" o Removing %",as.integer(log10( NumCorNMB ))+1,"d core/empty orbitals (%", as.integer(log10( Number_of_electrons_in_DM-DMtraceNVB ))+6, ".4fe) from the ",EDDB_basis," basis... DONE!\n", collapse=""), NumCorNMB, Number_of_electrons_in_DM-DMtraceNVB ) ) } Number_of_electrons_in_DM = DMtraceNVB } # suppressWarnings( rm( FCHK_file, FirstElementIndex, i, idx, ListOfNBOFiles, NBO_content_len, NumberOfElementsToSkip, DMtrace, ftrc, ListOfFchkFiles, NBO_content, NBO_file, SELECTED_FileNumber, DMtraceNMB, isNAOaNMB, VectorOfAtomLabels, DMtraceNVB ) ) if (!SAVE_EDDB_TO_GAUSSIAN_FCHK) { suppressWarnings( rm( FCHK_content, FCHK_content_len ) ) } # ____________________________________________________________________________ # # Definition of the EDDB function # ____________________________________________________________________________ # Perform_BOP_analysis <- function( DM, BM, lbAT, Natoms, BOtrsh , pbinc) { nNAOs = nrow( DM ) DDM = matrix(0, nNAOs, nNAOs) aDDM = matrix(0, nNAOs, nNAOs) for (ATOM_X in 1:Natoms) { if (diag(BM)[ATOM_X] > 1) { isBOclc = rep(FALSE, Natoms) numbo = rep(0, Natoms) for (q in 1:Natoms) { if ((BM[ATOM_X,q]==1)&&(q!=ATOM_X)) { numbo[q] = min(lbAT[q,2] - lbAT[q,1] + 1,lbAT[ATOM_X,2] - lbAT[ATOM_X,1] + 1) isBOclc[q]=TRUE } } BOXv = matrix(0, nNAOs, sum(numbo)) BOXn = matrix(0, sum(numbo), sum(numbo)) BOXx = matrix(0, sum(numbo), sum(numbo)) aBOXv = matrix(0, nNAOs, sum(numbo)) for (q in 1:Natoms) { if (isBOclc[q]) { nA = lbAT[q,2] - lbAT[q,1] + 1 nX = lbAT[ATOM_X,2] - lbAT[ATOM_X,1] + 1 nAX=nA+nX As = 1 Ak = nA Xs = Ak+1 Xk = Ak + nX AX = matrix(0,nAX,nAX) AX[As:Ak,Xs:Xk] = round(DM[lbAT[q,1]:lbAT[q,2], lbAT[ATOM_X,1]:lbAT[ATOM_X,2]], digits=7) AX[Xs:Xk,As:Ak] = t(AX[As:Ak,Xs:Xk]) eigenAX = eigen(AX,symmetric=TRUE) AXv = as.matrix(eigenAX$vectors[, 1:numbo[q]]) AXn = as.matrix(eigenAX$values[1:numbo[q]]^2) aAXv = as.matrix((as.matrix(eigenAX$vectors[, (nAX-numbo[q]+1):nAX]))[, numbo[q]:1]) if (q==1) { diag(BOXn)[1:numbo[1]] = AXn BOXv[lbAT[q,1]:lbAT[q,2], 1:numbo[1]] = AXv[As:Ak,] BOXv[lbAT[ATOM_X,1]:lbAT[ATOM_X,2], 1:numbo[1]] = AXv[Xs:Xk,] aBOXv[lbAT[q,1]:lbAT[q,2], 1:numbo[1]] = aAXv[As:Ak,] aBOXv[lbAT[ATOM_X,1]:lbAT[ATOM_X,2], 1:numbo[1]] = aAXv[Xs:Xk,] } else { diag(BOXn)[(sum(numbo[1:(q-1)])+1):(sum(numbo[1:q]))] = AXn BOXv[lbAT[q,1]:lbAT[q,2], (sum(numbo[1:(q-1)])+1):(sum(numbo[1:q]))] = AXv[As:Ak,] BOXv[lbAT[ATOM_X,1]:lbAT[ATOM_X,2], (sum(numbo[1:(q-1)])+1):(sum(numbo[1:q]))] = AXv[Xs:Xk,] aBOXv[lbAT[q,1]:lbAT[q,2], (sum(numbo[1:(q-1)])+1):(sum(numbo[1:q]))] = aAXv[As:Ak,] aBOXv[lbAT[ATOM_X,1]:lbAT[ATOM_X,2], (sum(numbo[1:(q-1)])+1):(sum(numbo[1:q]))] = aAXv[Xs:Xk,] } if ( sum(AXn) < NATURAL_OCCUPATION_THRSHLD ) { BM[q,ATOM_X] = 0 ; BM[ATOM_X,q] = 0 } } } for (qm in 1:(Natoms-1)) { for (qn in (qm+1):Natoms) { if ((BM[qm,ATOM_X]==1)&&(BM[qn,ATOM_X]==1)) { ATOM_A = qm ATOM_B = qn nA = lbAT[ATOM_A,2] - lbAT[ATOM_A,1] + 1 nB = lbAT[ATOM_B,2] - lbAT[ATOM_B,1] + 1 nX = lbAT[ATOM_X,2] - lbAT[ATOM_X,1] + 1 nAXB = nA + nB + nX As = 1 Ak = nA Xs = Ak+1 Xk = Ak + nX Bs = Xk + 1 Bk = Xk + nB AXB = matrix(0,nAXB,nAXB) AXB[As:Ak,Xs:Xk] = round(DM[lbAT[ATOM_A, 1]:lbAT[ATOM_A, 2],lbAT[ATOM_X, 1]:lbAT[ATOM_X, 2]], digits=7) AXB[Xs:Xk,As:Ak] = t(AXB[As:Ak,Xs:Xk]) AXB[Xs:Xk,Bs:Bk] = round(DM[lbAT[ATOM_X, 1]:lbAT[ATOM_X, 2],lbAT[ATOM_B, 1]:lbAT[ATOM_B, 2]], digits=7) AXB[Bs:Bk,Xs:Xk] = t(AXB[Xs:Xk,Bs:Bk]) if ( IGNORE_DSYEVR_LAPACK_ERROR ) eigenAXB = try( eigen(AXB, symmetric=TRUE), silent=TRUE ) else eigenAXB = eigen(AXB, symmetric=TRUE) # Possible DSYEVR errors for unknown reason... REVISION NEEDED !!!! AXB_vec = try( eigenAXB$vectors[,1:nX], silent=TRUE ) AXB_val = try( diag(eigenAXB$values[1:nX]^2,nX,nX), silent=TRUE) if ( !is.numeric(AXB_val) ) DSYEVR_Err_Counter <<- DSYEVR_Err_Counter + 1 else { nn = diag(AXB_val) AX_val = matrix(0,numbo[ATOM_A],numbo[ATOM_A]) if (ATOM_A==1) diag(AX_val) = diag(BOXn)[1:numbo[1]] else diag(AX_val) = diag(BOXn)[(sum(numbo[1:(ATOM_A-1)])+1):(sum(numbo[1:ATOM_A]))] if (ATOM_A==1) AX_vec = BOXv[c(lbAT[ATOM_A,1]:lbAT[ATOM_A,2], lbAT[ATOM_X,1]:lbAT[ATOM_X,2], lbAT[ATOM_B,1]:lbAT[ATOM_B,2]),1:numbo[1]] else AX_vec = BOXv[c(lbAT[ATOM_A,1]:lbAT[ATOM_A,2], lbAT[ATOM_X,1]:lbAT[ATOM_X,2], lbAT[ATOM_B,1]:lbAT[ATOM_B,2]), (sum(numbo[1:(ATOM_A-1)])+1):(sum(numbo[1:ATOM_A]))] BX_val = matrix(0,numbo[ATOM_B], numbo[ATOM_B]) if (ATOM_B==1) diag(BX_val) = diag(BOXn)[1:numbo[1]] else diag(BX_val) = diag(BOXn)[(sum(numbo[1:(ATOM_B-1)])+1):(sum(numbo[1:ATOM_B]))] if (ATOM_B==1) BX_vec = BOXv[c(lbAT[ATOM_A,1]:lbAT[ATOM_A,2], lbAT[ATOM_X,1]:lbAT[ATOM_X,2], lbAT[ATOM_B,1]:lbAT[ATOM_B,2]), 1:numbo[1]] else BX_vec = BOXv[c(lbAT[ATOM_A,1]:lbAT[ATOM_A,2], lbAT[ATOM_X,1]:lbAT[ATOM_X,2], lbAT[ATOM_B,1]:lbAT[ATOM_B,2]), (sum(numbo[1:(ATOM_B-1)])+1):(sum(numbo[1:ATOM_B]))] for (d in 1:nX) { if(AXB_val[d,d]>=BOtrsh) { AXB_val[d,d] = AXB_val[d,d]/BOtrsh } else { AXB_val[d,d] = 1 } } for (d in 1:numbo[ATOM_A]) { if(AX_val[d,d]>=BOtrsh) { AX_val[d,d] = AX_val[d,d]/BOtrsh } else { AX_val[d,d] = 1 } } for (d in 1:numbo[ATOM_B]) { if(BX_val[d,d]>=BOtrsh) { BX_val[d,d] = BX_val[d,d]/BOtrsh } else { BX_val[d,d] = 1 } } CC = matrix(0,numbo[ATOM_A]+numbo[ATOM_B], nX) AXBvAXBv = AXB_vec%*%AXB_val CC[1:numbo[ATOM_A],] = t(AX_vec%*%AX_val)%*%AXBvAXBv CC[(numbo[ATOM_A]+1):(numbo[ATOM_A]+numbo[ATOM_B]),] = t(BX_vec%*%BX_val)%*%AXBvAXBv XCCX = t(CC)%*%CC eigenXCCX = eigen(XCCX, symmetric=TRUE) XCCX_vec = eigenXCCX$vectors XCCX_val = diag(eigenXCCX$values,nX,nX) for (i in 1:nX) { if (XCCX_val[i,i]1) { for (i in 2:nX) { Thu[i] = Tht[i] - Tht[i-1] } } deloc = Thu * nn for (ii in 1:length(deloc)) { if (deloc[ii]<0) { smtp = 0 smii = 1 for (jj in 1:length(deloc)) { if ((jj!=ii)&&(Thu[jj]>0)) { tmpsm = abs( (t(CCo[1:numbo[ATOM_A], ii])%*%CCo[1:numbo[ATOM_A], jj])*(t(CCo[(numbo[ATOM_A] + 1):(numbo[ATOM_A] + numbo[ATOM_B]), ii])%*%CCo[(numbo[ATOM_A] + 1):(numbo[ATOM_A]+numbo[ATOM_B]), jj])) if (tmpsm>smtp) { smtp = tmpsm smii = jj } } } deloc[smii] = deloc[smii] + deloc[ii] deloc[ii] = 0 } } for (oi in (sum(numbo[1:ATOM_A]) - numbo[ATOM_A]+1):sum(numbo[1:ATOM_A])) { smboxx = sum(deloc*CCo[(oi - (sum(numbo[1:ATOM_A]) - numbo[ATOM_A])), ]^2) if ((BOXx[oi,oi]BOXn[oi,oi])) { BOXx[oi,oi] = BOXn[oi,oi] } } } for (oi in (sum(numbo[1:ATOM_B]) - numbo[ATOM_B] + 1):sum(numbo[1:ATOM_B])) { smboxx = sum(deloc*CCo[(oi - (sum(numbo[1:ATOM_B]) - numbo[ATOM_B])+numbo[ATOM_A]), ]^2) if ((BOXx[oi,oi]BOXn[oi,oi])) { BOXx[oi,oi] = BOXn[oi,oi] } } } } } } } DDM = DDM + BOXv%*%BOXx%*%t(BOXv) aDDM = aDDM + aBOXv%*%BOXx%*%t(aBOXv) setTxtProgressBar(progressBar,pbinc) pbinc = pbinc + 1 } } return( DDM + aDDM ) } # ____________________________________________________________________________ # # Saving the results into a new .fchk file # ____________________________________________________________________________ # Save_EDDB_to_FCHK <- function( FCHKin, FCHKout, EDDB, SDDB, NOBDn, nscal, NOBDv, AONAO, nAOs, nMOs, isoo ) { if ( isSpinUnrestricted && isoo) { Mtmp = EDDB EDDB = ( EDDB + SDDB ) / 2 # Now EDDB is EDDB_Alpha SDDB = ( Mtmp - SDDB ) / 2 # Now SDDB is EDDB_Beta Mtmp = eigen(EDDB,symmetric=TRUE) NOBDn = Mtmp$values ; NOBDv = Mtmp$vectors ; Mtmp = 0 # Now NOBD refers to EDDB_Alpha } NOBDt = rep(0 , max( nMOs , length( NOBDn ) ) ) NOBDt[1:length( NOBDn )] = NOBDn * nscal NOBDn = NOBDt[1:nMOs] NOBDv_AO = AONAO %*% NOBDv if (ncol( NOBDv_AO ) > nMOs ) NOBDv_AO = NOBDv_AO[, 1:nMOs ] if (ncol( NOBDv_AO ) < nMOs ) NOBDv_AO = cbind( NOBDv_AO, matrix(0, nAOs, ( nMOs - ncol( NOBDv_AO ) ) ) ) NOBDv_AO = as.vector(NOBDv_AO) # FCHKin_len = length( FCHKin ) fchktxttyp = "Alpha Orbital Energies R" fchktxttyp = scan(text=fchktxttyp,what="",sep="",quiet=TRUE) ; fchktxttyp_len = length(fchktxttyp) fchktxttyp_found = FALSE ; idx = 1 while( (!fchktxttyp_found) && idx < FCHKin_len ) { fchkline = scan(text=FCHKin[idx],what="",sep="",quiet=TRUE) ; fchkline_len = length(fchkline) if ( suppressWarnings( grepl( tolower( paste0( fchktxttyp, collapse="" ) ), tolower( paste0( fchkline[1:fchktxttyp_len], collapse="") ), fixed=FALSE ) ) ) { fchktxttyp_found = TRUE idx = idx + 1 } if ( (!fchktxttyp_found) ) { if ( fchkline_len > 3 && suppressWarnings( grepl( "in=", tolower( paste0( fchkline[(fchkline_len-2):fchkline_len], collapse="") ), fixed=FALSE ) ) ) idx = idx + max( 1, suppressWarnings( as.integer( fchkline[fchkline_len] ) %/% 6 ) ) else if ( fchkline_len > 3 && suppressWarnings( grepl( "rn=", tolower( paste0( fchkline[(fchkline_len-2):fchkline_len], collapse="") ), fixed=FALSE ) ) ) idx = idx + max( 1, suppressWarnings( as.integer( fchkline[fchkline_len] ) %/% 5 ) ) else idx = idx + 1 } } if ( idx >= FCHKin_len ) { cat(" FAILED!\n Alpha orbital energies not found.\n\n Press to terminate the RunEDDB script.\n") if(INTERACTIVE_MODE) nullek=readLines(STD_IN,n=1) else cat("\n") ; close(STD_IN) ; quit(save="no") } nFullLines = nMOs%/%5 ; isAddLine = nMOs%%5 > 0 for ( i in 1:nFullLines ) { FCHKin[idx] = paste0( sprintf( "%16.8E",round( NOBDn[ ((i-1)*5+1):(5*i) ],digits=15)), collapse="" ) idx = idx + 1 } if ( isTRUE( isAddLine ) ) FCHKin[idx] = paste0( sprintf( "%16.8E",round( NOBDn[ ( nMOs - nMOs%%5 + 1 ):nMOs ],digits=15)), collapse="" ) # fchktxttyp = "Alpha MO coefficients R" fchktxttyp = scan(text=fchktxttyp,what="",sep="",quiet=TRUE) ; fchktxttyp_len = length(fchktxttyp) fchktxttyp_found = FALSE ; idx = 1 while( (!fchktxttyp_found) && idx < FCHKin_len ) { fchkline = scan(text=FCHKin[idx],what="",sep="",quiet=TRUE) ; fchkline_len = length(fchkline) if ( suppressWarnings( grepl( tolower( paste0( fchktxttyp, collapse="" ) ), tolower( paste0( fchkline[1:fchktxttyp_len], collapse="") ), fixed=FALSE ) ) ) { fchktxttyp_found = TRUE idx = idx + 1 } if ( (!fchktxttyp_found) ) { if ( fchkline_len > 3 && suppressWarnings( grepl( "in=", tolower( paste0( fchkline[(fchkline_len-2):fchkline_len], collapse="") ), fixed=FALSE ) ) ) idx = idx + max( 1, suppressWarnings( as.integer( fchkline[fchkline_len] ) %/% 6 ) ) else if ( fchkline_len > 3 && suppressWarnings( grepl( "rn=", tolower( paste0( fchkline[(fchkline_len-2):fchkline_len], collapse="") ), fixed=FALSE ) ) ) idx = idx + max( 1, suppressWarnings( as.integer( fchkline[fchkline_len] ) %/% 5 ) ) else idx = idx + 1 } } if ( idx >= FCHKin_len ) { cat(" FAILED!\n Alpha MO coefficients not found.\n\n Press to terminate the RunEDDB script.\n") if(INTERACTIVE_MODE) nullek=readLines(STD_IN,n=1) else cat("\n") ; close(STD_IN) ; quit(save="no") } nFullLines = ( nMOs*nAOs )%/%5 ; isAddLine = ( nMOs*nAOs )%%5 > 0 for ( i in 1:nFullLines ) { FCHKin[idx] = paste0( sprintf( "%16.8E",round( NOBDv_AO[ ((i-1)*5+1):(5*i) ],digits=15)), collapse="" ) idx = idx + 1 } if ( isTRUE( isAddLine ) ) FCHKin[idx] = paste0( sprintf( "%16.8E",round( NOBDv_AO[ ( ( nMOs*nAOs ) - ( nMOs*nAOs )%%5 + 1 ):( nMOs*nAOs ) ],digits=15)), collapse="" ) # if ( isSpinUnrestricted && isoo) { Mtmp = eigen(SDDB,symmetric=TRUE) NOBDn = Mtmp$values ; NOBDv = Mtmp$vectors ; Mtmp = 0 # Now NOBD refers to EDDB_Beta NOBDt = rep(0 , max( nMOs , length( NOBDn ) ) ) NOBDt[1:length( NOBDn )] = NOBDn * nscal NOBDn = NOBDt[1:nMOs] NOBDv_AO = AONAO %*% NOBDv if (ncol( NOBDv_AO ) > nMOs ) NOBDv_AO = NOBDv_AO[, 1:nMOs ] if (ncol( NOBDv_AO ) < nMOs ) NOBDv_AO = cbind( NOBDv_AO, matrix(0, nAOs, ( nMOs - ncol( NOBDv_AO ) ) ) ) NOBDv_AO = as.vector(NOBDv_AO) # fchktxttyp = "Beta Orbital Energies R" fchktxttyp = scan(text=fchktxttyp,what="",sep="",quiet=TRUE) ; fchktxttyp_len = length(fchktxttyp) fchktxttyp_found = FALSE ; idx = 1 while( (!fchktxttyp_found) && idx < FCHKin_len ) { fchkline = scan(text=FCHKin[idx],what="",sep="",quiet=TRUE) ; fchkline_len = length(fchkline) if ( suppressWarnings( grepl( tolower( paste0( fchktxttyp, collapse="" ) ), tolower( paste0( fchkline[1:fchktxttyp_len], collapse="") ), fixed=FALSE ) ) ) { fchktxttyp_found = TRUE idx = idx + 1 } if ( (!fchktxttyp_found) ) { if ( fchkline_len > 3 && suppressWarnings( grepl( "in=", tolower( paste0( fchkline[(fchkline_len-2):fchkline_len], collapse="") ), fixed=FALSE ) ) ) idx = idx + max( 1, suppressWarnings( as.integer( fchkline[fchkline_len] ) %/% 6 ) ) else if ( fchkline_len > 3 && suppressWarnings( grepl( "rn=", tolower( paste0( fchkline[(fchkline_len-2):fchkline_len], collapse="") ), fixed=FALSE ) ) ) idx = idx + max( 1, suppressWarnings( as.integer( fchkline[fchkline_len] ) %/% 5 ) ) else idx = idx + 1 } } if ( idx >= FCHKin_len ) { cat(" FAILED!\n Beta orbital energies not found.\n\n Press to terminate the RunEDDB script.\n") if(INTERACTIVE_MODE) nullek=readLines(STD_IN,n=1) else cat("\n") ; close(STD_IN) ; quit(save="no") } nFullLines = nMOs%/%5 ; isAddLine = nMOs%%5 > 0 for ( i in 1:nFullLines ) { FCHKin[idx] = paste0( sprintf( "%16.8E",round( NOBDn[ ((i-1)*5+1):(5*i) ],digits=15)), collapse="" ) idx = idx + 1 } if ( isTRUE( isAddLine ) ) FCHKin[idx] = paste0( sprintf( "%16.8E",round( NOBDn[ ( nMOs - nMOs%%5 + 1 ):nMOs ],digits=15)), collapse="" ) # fchktxttyp = "Beta MO coefficients R" fchktxttyp = scan(text=fchktxttyp,what="",sep="",quiet=TRUE) ; fchktxttyp_len = length(fchktxttyp) fchktxttyp_found = FALSE ; idx = 1 while( (!fchktxttyp_found) && idx < FCHKin_len ) { fchkline = scan(text=FCHKin[idx],what="",sep="",quiet=TRUE) ; fchkline_len = length(fchkline) if ( suppressWarnings( grepl( tolower( paste0( fchktxttyp, collapse="" ) ), tolower( paste0( fchkline[1:fchktxttyp_len], collapse="") ), fixed=FALSE ) ) ) { fchktxttyp_found = TRUE idx = idx + 1 } if ( (!fchktxttyp_found) ) { if ( fchkline_len > 3 && suppressWarnings( grepl( "in=", tolower( paste0( fchkline[(fchkline_len-2):fchkline_len], collapse="") ), fixed=FALSE ) ) ) idx = idx + max( 1, suppressWarnings( as.integer( fchkline[fchkline_len] ) %/% 6 ) ) else if ( fchkline_len > 3 && suppressWarnings( grepl( "rn=", tolower( paste0( fchkline[(fchkline_len-2):fchkline_len], collapse="") ), fixed=FALSE ) ) ) idx = idx + max( 1, suppressWarnings( as.integer( fchkline[fchkline_len] ) %/% 5 ) ) else idx = idx + 1 } } if ( idx >= FCHKin_len ) { cat(" FAILED!\n Beta MO coefficients not found.\n\n Press to terminate the RunEDDB script.\n") if(INTERACTIVE_MODE) nullek=readLines(STD_IN,n=1) else cat("\n") ; close(STD_IN) ; quit(save="no") } nFullLines = ( nMOs*nAOs )%/%5 ; isAddLine = ( nMOs*nAOs )%%5 > 0 for ( i in 1:nFullLines ) { FCHKin[idx] = paste0( sprintf( "%16.8E",round( NOBDv_AO[ ((i-1)*5+1):(5*i) ],digits=15)), collapse="" ) idx = idx + 1 } if ( isTRUE( isAddLine ) ) FCHKin[idx] = paste0( sprintf( "%16.8E",round( NOBDv_AO[ ( ( nMOs*nAOs ) - ( nMOs*nAOs )%%5 + 1 ):( nMOs*nAOs ) ],digits=15)), collapse="" ) } # fchktxttyp = "Total SCF Density R" fchktxttyp = scan(text=fchktxttyp,what="",sep="",quiet=TRUE) ; fchktxttyp_len = length(fchktxttyp) fchktxttyp_found = FALSE ; idx = 1 while( (!fchktxttyp_found) && idx < FCHKin_len ) { fchkline = scan(text=FCHKin[idx],what="",sep="",quiet=TRUE) ; fchkline_len = length(fchkline) if ( suppressWarnings( grepl( tolower( paste0( fchktxttyp, collapse="" ) ), tolower( paste0( fchkline[1:fchktxttyp_len], collapse="") ), fixed=FALSE ) ) ) { fchktxttyp_found = TRUE idx = idx + 1 } if ( (!fchktxttyp_found) ) { if ( fchkline_len > 3 && suppressWarnings( grepl( "in=", tolower( paste0( fchkline[(fchkline_len-2):fchkline_len], collapse="") ), fixed=FALSE ) ) ) idx = idx + max( 1, suppressWarnings( as.integer( fchkline[fchkline_len] ) %/% 6 ) ) else if ( fchkline_len > 3 && suppressWarnings( grepl( "rn=", tolower( paste0( fchkline[(fchkline_len-2):fchkline_len], collapse="") ), fixed=FALSE ) ) ) idx = idx + max( 1, suppressWarnings( as.integer( fchkline[fchkline_len] ) %/% 5 ) ) else idx = idx + 1 } } if ( idx >= FCHKin_len ) { cat(" FAILED!\n Total SCF Density not found.\n\n Press to terminate the RunEDDB script.\n") if(INTERACTIVE_MODE) nullek=readLines(STD_IN,n=1) else cat("\n") ; close(STD_IN) ; quit(save="no") } EDDB_AO = AONAO %*% EDDB %*% t( AONAO ) EDDB_AO = as.vector( EDDB_AO[upper.tri( EDDB_AO, diag=TRUE )] ); EDDB_AO_len = length( EDDB_AO ) nFullLines = EDDB_AO_len%/%5 ; isAddLine = EDDB_AO_len%%5 > 0 for ( i in 1:nFullLines ) { FCHKin[idx] = paste0( sprintf( "%16.8E",round( EDDB_AO[ ((i-1)*5+1):(5*i) ],digits=15)), collapse="" ) idx = idx + 1 } if ( isTRUE( isAddLine ) ) FCHKin[idx] = paste0( sprintf( "%16.8E",round( EDDB_AO[ ( EDDB_AO_len - EDDB_AO_len%%5 + 1 ):EDDB_AO_len ],digits=15)), collapse="" ) # if (isSpinUnrestricted) { fchktxttyp = "Spin SCF Density R" fchktxttyp = scan(text=fchktxttyp,what="",sep="",quiet=TRUE) ; fchktxttyp_len = length(fchktxttyp) fchktxttyp_found = FALSE ; idx = 1 while( (!fchktxttyp_found) && idx < FCHKin_len ) { fchkline = scan(text=FCHKin[idx],what="",sep="",quiet=TRUE) ; fchkline_len = length(fchkline) if ( suppressWarnings( grepl( tolower( paste0( fchktxttyp, collapse="" ) ), tolower( paste0( fchkline[1:fchktxttyp_len], collapse="") ), fixed=FALSE ) ) ) { fchktxttyp_found = TRUE idx = idx + 1 } if ( (!fchktxttyp_found) ) { if ( fchkline_len > 3 && suppressWarnings( grepl( "in=", tolower( paste0( fchkline[(fchkline_len-2):fchkline_len], collapse="") ), fixed=FALSE ) ) ) idx = idx + max( 1, suppressWarnings( as.integer( fchkline[fchkline_len] ) %/% 6 ) ) else if ( fchkline_len > 3 && suppressWarnings( grepl( "rn=", tolower( paste0( fchkline[(fchkline_len-2):fchkline_len], collapse="") ), fixed=FALSE ) ) ) idx = idx + max( 1, suppressWarnings( as.integer( fchkline[fchkline_len] ) %/% 5 ) ) else idx = idx + 1 } } if ( idx >= FCHKin_len ) { cat(" FAILED!\n Spin SCF Density not found.\n\n Press to terminate the RunEDDB script.\n") if(INTERACTIVE_MODE) nullek=readLines(STD_IN,n=1) else cat("\n") ; close(STD_IN) ; quit(save="no") } EDDB_AO = AONAO %*% SDDB %*% t( AONAO ) EDDB_AO = as.vector( EDDB_AO[upper.tri( EDDB_AO, diag=TRUE )] ); EDDB_AO_len = length( EDDB_AO ) nFullLines = EDDB_AO_len%/%5 ; isAddLine = EDDB_AO_len%%5 > 0 for ( i in 1:nFullLines ) { FCHKin[idx] = paste0( sprintf( "%16.8E",round( EDDB_AO[ ((i-1)*5+1):(5*i) ],digits=15)), collapse="" ) idx = idx + 1 } if ( isTRUE( isAddLine ) ) FCHKin[idx] = paste0( sprintf( "%16.8E",round( EDDB_AO[ ( EDDB_AO_len - EDDB_AO_len%%5 + 1 ):EDDB_AO_len ],digits=15)), collapse="" ) } cat(FCHKin, file=FCHKout, sep="\n", fill=FALSE, labels=NULL, append=FALSE ) invisible() } # ____________________________________________________________________________ # cat("\n __________________ELECTRON_DENSITY_OF_DELOCALIZED_BONDS___________________ ") cat("\n | |") cat("\n | 1) for the entire molecular system (including hydrogens)......EDDB_G(r) |") cat("\n | 2) for the entire molecular system (excluding hydrogens)......EDDB_H(r) |") cat("\n | 3) for the molecular fragment (excluding non-local effects)...EDDB_F(r) |") cat("\n | 4) for the molecular fragment (including non-local effects)...EDDB_E(r) |") cat("\n | 5) for the particular delocalization pathway..................EDDB_P(r) |") cat("\n | 0) skip the EDDB calculations and perform the standard NPA analysis. |") cat("\n |__________________________________________________________________________|\n\n") # if ((!(is.na(DEFAULT_EDDB_FUNCTION_TYPE)))&&(DEFAULT_EDDB_FUNCTION_TYPE %in% c( 0,1:5,10,20,30,40 ))) { cat( paste0(" o Select type of the EDDB function and press :\n > ",DEFAULT_EDDB_FUNCTION_TYPE,"\n",collapse="") ) } else DEFAULT_EDDB_FUNCTION_TYPE = NA while (is.na(DEFAULT_EDDB_FUNCTION_TYPE) || !(DEFAULT_EDDB_FUNCTION_TYPE %in% c( 0,1:5,10,20,30,40 ) ) ) { cat(" o Select type of the EDDB function and press :\n > ") if (!INTERACTIVE_MODE) { cat("\n Interactive mode disabled! Script terminated.\n\n") ; close(STD_IN) ; quit(save="no") } DEFAULT_EDDB_FUNCTION_TYPE = suppressWarnings( as.integer( scan( text=readLines(STD_IN,n=1) , what="" , sep="" , quiet=TRUE )[1] ) ) } if (DEFAULT_EDDB_FUNCTION_TYPE==0) SKIP_EDDB__AND_PERFORM_NPA = TRUE if (DEFAULT_EDDB_FUNCTION_TYPE==1) EDDB_type = "EDDB_G" if (DEFAULT_EDDB_FUNCTION_TYPE==10) EDDB_type = "EDDB_G" if (DEFAULT_EDDB_FUNCTION_TYPE==2) EDDB_type = "EDDB_H" if (DEFAULT_EDDB_FUNCTION_TYPE==20) EDDB_type = "EDDB_H" if (DEFAULT_EDDB_FUNCTION_TYPE==3) EDDB_type = "EDDB_F" if (DEFAULT_EDDB_FUNCTION_TYPE==30) EDDB_type = "EDDB_F" if (DEFAULT_EDDB_FUNCTION_TYPE==4) EDDB_type = "EDDB_E" if (DEFAULT_EDDB_FUNCTION_TYPE==40) EDDB_type = "EDDB_E" if (DEFAULT_EDDB_FUNCTION_TYPE==5) EDDB_type = "EDDB_P" if (SKIP_EDDB__AND_PERFORM_NPA) EDDB_type = paste0("ED_",EDDB_basis,collapse="") if (DEFAULT_EDDB_FUNCTION_TYPE %in% c( 1,2,3,4 ) ) { Matrix_of_atomic_connectivity = matrix(0,Number_of_atoms,Number_of_atoms) for (i in 2:Number_of_atoms) { for (j in 1:(i-1)) { if (!isSpinUnrestricted) { Bond_order = sum( Matrix_of_Alpha_Density[ NAOtoAtomArray[i,1]:NAOtoAtomArray[i,2] , NAOtoAtomArray[j,1]:NAOtoAtomArray[j,2] ]^2 ) } else { Bond_order = 2* ( sum( Matrix_of_Alpha_Density[ NAOtoAtomArray[i,1]:NAOtoAtomArray[i,2] , NAOtoAtomArray[j,1]:NAOtoAtomArray[j,2] ]^2 ) + sum( Matrix_of_Beta_Density[ NAOtoAtomArray[i,1]:NAOtoAtomArray[i,2] , NAOtoAtomArray[j,1]:NAOtoAtomArray[j,2] ]^2 ) ) } if ( Bond_order >= 10*NATURAL_OCCUPATION_THRSHLD ) { Matrix_of_atomic_connectivity[i,j] = 1 Matrix_of_atomic_connectivity[j,i] = 1 } } } for (i in 1:Number_of_atoms) Matrix_of_atomic_connectivity[i,i] = sum( Matrix_of_atomic_connectivity[,i] ) } if (DEFAULT_EDDB_FUNCTION_TYPE %in% c( 0,10,20,30,40 ) ) { Matrix_of_atomic_connectivity = matrix(1,Number_of_atoms,Number_of_atoms) diag(Matrix_of_atomic_connectivity) = Number_of_atoms - 1 } if (DEFAULT_EDDB_FUNCTION_TYPE %in% c( 2, 20 ) ) { for (i in 1:Number_of_atoms) if (Array_of_atomic_numbers[i]==1) { Matrix_of_atomic_connectivity[i,] = 0 Matrix_of_atomic_connectivity[,i] = 0 } for (i in 1:Number_of_atoms) Matrix_of_atomic_connectivity[i,i] = sum( Matrix_of_atomic_connectivity[,i] ) - 1 } if (DEFAULT_EDDB_FUNCTION_TYPE %in% c( 3,30,4,40 ) ) { cat(" o Define molecular fragment ( e.g. 1:6,9,10 ) and press ") ContinueLoop = TRUE LoopCounter = 0 while( ContinueLoop ) { ContinueLoop = FALSE LoopCounter = LoopCounter + 1 cat("\n > ") if ((LoopCounter==1)&&(exists("SELECTED_Fragment_readline"))) cat(paste0(SELECTED_Fragment_readline,"\n",collapse="")) if (((LoopCounter==1)&&(!exists("SELECTED_Fragment_readline"))&&(!INTERACTIVE_MODE))||((LoopCounter>1)&&(!INTERACTIVE_MODE))) { cat("\n Interactive mode disabled! Script terminated.\n\n") ; close(STD_IN) ; quit(save="no") } if ((LoopCounter>1)||(!exists("SELECTED_Fragment_readline"))) SELECTED_Fragment_readline = suppressWarnings( readLines(STD_IN,n=1) ) SELECTED_Fragment_splited = strsplit( SELECTED_Fragment_readline, ',', fixed=TRUE)[[1]] if (length(SELECTED_Fragment_splited)<1) { ContinueLoop = TRUE } else { SELECTED_Fragment = rep( FALSE , Number_of_atoms ) for (i in 1:length(SELECTED_Fragment_splited)) { if (!is.na(suppressWarnings( as.integer( SELECTED_Fragment_splited[i] )))) { if (as.integer( SELECTED_Fragment_splited[i] )<=Number_of_atoms) SELECTED_Fragment[as.integer( SELECTED_Fragment_splited[i] )] = TRUE else ContinueLoop = TRUE } else { SplitedItems = strsplit( SELECTED_Fragment_splited[i] , ':', fixed=TRUE)[[1]] if (length(SplitedItems)!=2) { ContinueLoop = TRUE } else { if ( (is.na(suppressWarnings( as.integer( SplitedItems[1] )))) || (is.na(suppressWarnings( as.integer( SplitedItems[1] )))) ) { ContinueLoop = TRUE } else { if ( (as.integer( SplitedItems[1] ) <= Number_of_atoms) && (as.integer( SplitedItems[2] ) <= Number_of_atoms) ) SELECTED_Fragment[as.integer( SplitedItems[1] ):as.integer( SplitedItems[2] )] = TRUE else ContinueLoop = TRUE } } } } if (sum(as.integer(SELECTED_Fragment)) < 3 ) ContinueLoop = TRUE } } SELECTED_Fragment = (1:Number_of_atoms)[SELECTED_Fragment] if (DEFAULT_EDDB_FUNCTION_TYPE %in% c( 3,30 ) ) { for (i in 1:Number_of_atoms) { if (!(i %in% SELECTED_Fragment)) { Matrix_of_atomic_connectivity[,i] = 0 Matrix_of_atomic_connectivity[i,] = 0 } } } diag(Matrix_of_atomic_connectivity) = 0 for (i in 1:Number_of_atoms) Matrix_of_atomic_connectivity[i,i] = sum( Matrix_of_atomic_connectivity[,i] ) } if (DEFAULT_EDDB_FUNCTION_TYPE %in% c( 5 ) ) { cat(" o Define cyclic delocalization pathway ( e.g. 1-2-3-4-5-6-1 ) and press ") ContinueLoop = TRUE LoopCounter = 0 while( ContinueLoop ) { ContinueLoop = FALSE LoopCounter = LoopCounter + 1 SingleCircuit = FALSE cat("\n > ") if ((LoopCounter==1)&&(exists("SELECTED_Pathway_readline"))) cat(paste0(SELECTED_Pathway_readline,"\n",collapse="")) if (((LoopCounter==1)&&(!exists("SELECTED_Pathway_readline"))&&(!INTERACTIVE_MODE))||((LoopCounter>1)&&(!INTERACTIVE_MODE))) { cat("\n Interactive mode disabled! Script terminated.\n\n") ; close(STD_IN) ; quit(save="no") } if ((LoopCounter>1)||(!exists("SELECTED_Pathway_readline"))) SELECTED_Pathway_readline = suppressWarnings( readLines(STD_IN,n=1) ) SELECTED_Pathway_splited = strsplit( SELECTED_Pathway_readline, ',', fixed=TRUE)[[1]] if (length(SELECTED_Pathway_splited)<1) { ContinueLoop = TRUE } else { Matrix_of_atomic_connectivity = matrix( 0, Number_of_atoms, Number_of_atoms ) for (i in 1:length(SELECTED_Pathway_splited)) { SplitedItems = strsplit( SELECTED_Pathway_splited[i] , '-', fixed=TRUE)[[1]] if ( (length(SplitedItems)<2) || (is.na( sum( suppressWarnings( as.integer( SplitedItems ) ) ) ) ) ) { ContinueLoop = TRUE } else { for (j in 2:length(SplitedItems)) { Matrix_of_atomic_connectivity[ as.integer( SplitedItems[j] ) , as.integer( SplitedItems[j-1] ) ] = 1 Matrix_of_atomic_connectivity[ as.integer( SplitedItems[j-1] ) , as.integer( SplitedItems[j] ) ] = 1 } if ( as.integer( SplitedItems[1] ) == as.integer( SplitedItems[length( SplitedItems )] ) ) SingleCircuit = TRUE } } if (length(SELECTED_Pathway_splited)==1) SingleCircuit = TRUE else SingleCircuit = FALSE diag(Matrix_of_atomic_connectivity) = 0 for (i in 1:Number_of_atoms) Matrix_of_atomic_connectivity[i,i] = sum( Matrix_of_atomic_connectivity[,i] ) if ( sum( as.integer(diag(Matrix_of_atomic_connectivity)>1) ) < 1 ) ContinueLoop = TRUE } } } # # Using formal connectivity array taken from the formatted checkpoint file # if (CONNECTIVITY_OUT_FROM_FCHK && (DEFAULT_EDDB_FUNCTION_TYPE!=5) ) { for (i in 1:Number_of_atoms) for (j in 1:Number_of_atoms) if (!Matrix_of_formal_connectivity[i,j]) Matrix_of_atomic_connectivity[i,j] = 0 diag(Matrix_of_atomic_connectivity) = 0 for (i in 1:Number_of_atoms) Matrix_of_atomic_connectivity[i,i] = sum( Matrix_of_atomic_connectivity[,i] ) cat(paste0(" o Using formal connectivity array from ",FCHK_file_name,".\n", collapse="")) } # ____________________________________________________________________________ # if (SKIP_EDDB__AND_PERFORM_NPA) { cat(paste0(" o Calculation of the ", EDDB_type , "(r) function (1eRDM in ",EDDB_basis," basis)...", collapse="")) if (!isSpinUnrestricted) { if (isPostHF) suppressWarnings( rm( Natural_Occupations, Natural_Orbitals, QDmat, QDval ) ) EDDB_Total = Matrix_of_Alpha_Density } else { if (isPostHF) suppressWarnings( rm ( Natural_Occupations_Alpha, Natural_Occupations_Beta, Natural_Orbitals_Alpha, Natural_Orbitals_Beta, QDmat_Alpha, QDmat_Beta, QDval_Alpha, QDval_Beta ) ) EDDB_Alpha = Matrix_of_Alpha_Density EDDB_Beta = Matrix_of_Beta_Density EDDB_Total = EDDB_Alpha + EDDB_Beta } } else { if ( sum( as.integer(diag(Matrix_of_atomic_connectivity)>1) ) < 1 ) { cat(paste0(" o For the selected criteria ", EDDB_type , " predicts no delocalized bonds.\n\n Press to terminate the RunEDDB script.\n", collapse="")) if(INTERACTIVE_MODE) nullek=readLines(STD_IN,n=1) else cat("\n") ; close(STD_IN) ; quit(save="no") } cat(paste0(" o Calculation of the ", EDDB_type , "(r) function in the " , EDDB_basis , " basis...", collapse="")) if (INTERACTIVE_MODE) cat("\n") Matrix_of_Alpha_Density_Cur = 0 EDDB_Total = 0 EDDB_Total_Cur = 0 if (isSpinUnrestricted) { Matrix_of_Beta_Density_Cur = 0 EDDB_Alpha = 0 EDDB_Beta = 0 EDDB_Alpha_Cur = 0 EDDB_Beta_Cur = 0 } progressBar = txtProgressBar(min=0, max=( 1 + as.integer(isSpinUnrestricted) )*QDnos*length( (1:Number_of_atoms)[ diag(Matrix_of_atomic_connectivity)>1 ] ), initial=0, char=".", width=68, style=3) for (iQD in 1:QDnos ) { pbinc = 1 + (iQD-1)*length( (1:Number_of_atoms)[ diag(Matrix_of_atomic_connectivity)>1 ] ) if (!isSpinUnrestricted) { if (!isPostHF) Matrix_of_Alpha_Density_Cur = Matrix_of_Alpha_Density else Matrix_of_Alpha_Density_Cur = Matrix_of_Alpha_Density_Cur + 2*Natural_Orbitals[,QDmat[iQD,1]:QDmat[iQD,2]] %*% t( Natural_Orbitals[,QDmat[iQD,1]:QDmat[iQD,2]] ) EDDB_Total_Pre = EDDB_Total_Cur EDDB_Total_Cur = Matrix_of_Alpha_Density_Cur %*% Perform_BOP_analysis ( 0.5*Matrix_of_Alpha_Density_Cur , Matrix_of_atomic_connectivity, NAOtoAtomArray, Number_of_atoms, NATURAL_OCCUPATION_THRSHLD, pbinc ) %*% Matrix_of_Alpha_Density_Cur EDDB_Total = EDDB_Total + 0.5*QDval[iQD]*(EDDB_Total_Cur - EDDB_Total_Pre) } else { if (!isPostHF) Matrix_of_Alpha_Density_Cur = Matrix_of_Alpha_Density else Matrix_of_Alpha_Density_Cur = Matrix_of_Alpha_Density_Cur + Natural_Orbitals_Alpha[,QDmat_Alpha[iQD,1]:QDmat_Alpha[iQD,2]] %*% t( Natural_Orbitals_Alpha[,QDmat_Alpha[iQD,1]:QDmat_Alpha[iQD,2]] ) if (!isPostHF) Matrix_of_Beta_Density_Cur = Matrix_of_Beta_Density else Matrix_of_Beta_Density_Cur = Matrix_of_Beta_Density_Cur + Natural_Orbitals_Beta[,QDmat_Beta[iQD,1]:QDmat_Beta[iQD,2]] %*% t( Natural_Orbitals_Beta[,QDmat_Beta[iQD,1]:QDmat_Beta[iQD,2]] ) EDDB_Alpha_Pre = EDDB_Alpha_Cur EDDB_Alpha_Cur = 2 * Matrix_of_Alpha_Density_Cur %*% Perform_BOP_analysis ( Matrix_of_Alpha_Density_Cur , Matrix_of_atomic_connectivity, NAOtoAtomArray, Number_of_atoms, NATURAL_OCCUPATION_THRSHLD, pbinc ) %*% Matrix_of_Alpha_Density_Cur if (iQD <= QDnos_Alpha) EDDB_Alpha = EDDB_Alpha + QDval_Alpha[iQD]*(EDDB_Alpha_Cur - EDDB_Alpha_Pre) pbinc = 1 + iQD*length( (1:Number_of_atoms)[ diag(Matrix_of_atomic_connectivity)>1 ] ) EDDB_Beta_Pre = EDDB_Beta_Cur EDDB_Beta_Cur = 2 * Matrix_of_Beta_Density_Cur %*% Perform_BOP_analysis ( Matrix_of_Beta_Density_Cur , Matrix_of_atomic_connectivity, NAOtoAtomArray, Number_of_atoms, NATURAL_OCCUPATION_THRSHLD, pbinc ) %*% Matrix_of_Beta_Density_Cur if (iQD <= QDnos_Beta) EDDB_Beta = EDDB_Beta + QDval_Beta[iQD]*(EDDB_Beta_Cur - EDDB_Beta_Pre) } } if (isSpinUnrestricted) EDDB_Total = EDDB_Alpha + EDDB_Beta suppressWarnings( rm( Matrix_of_Alpha_Density_Cur, EDDB_Total_Cur, EDDB_Total_Pre ) ) if (isSpinUnrestricted) suppressWarnings( rm( Matrix_of_Beta_Density_Cur, EDDB_Alpha_Cur, EDDB_Alpha_Pre, EDDB_Beta_Cur, EDDB_Beta_Pre ) ) if (isPostHF) { if (!isSpinUnrestricted) suppressWarnings( rm( Natural_Occupations, Natural_Orbitals, QDmat, QDval ) ) else suppressWarnings( rm ( Natural_Occupations_Alpha, Natural_Occupations_Beta, Natural_Orbitals_Alpha, Natural_Orbitals_Beta, QDmat_Alpha, QDmat_Beta, QDval_Alpha, QDval_Beta ) ) } } if (!isSpinUnrestricted) { for (i in 1:Number_of_atoms) { if ( diag(Matrix_of_atomic_connectivity)[i] < 1 ) { EDDB_Total[ , NAOtoAtomArray[i,1]:NAOtoAtomArray[i,2] ] = 0 EDDB_Total[ NAOtoAtomArray[i,1]:NAOtoAtomArray[i,2] , ] = 0 } } if (DEFAULT_EDDB_FUNCTION_TYPE %in% c( 4,40 ) ) { for (i in 1:Number_of_atoms) { if (!(i %in% SELECTED_Fragment)) { EDDB_Total[ , NAOtoAtomArray[i,1]:NAOtoAtomArray[i,2] ] = 0 EDDB_Total[ NAOtoAtomArray[i,1]:NAOtoAtomArray[i,2] , ] = 0 } } } sumEDDB_Total = sum( diag( EDDB_Total ) ) } else { for (i in 1:Number_of_atoms) { if ( diag(Matrix_of_atomic_connectivity)[i] < 1 ) { EDDB_Alpha[ , NAOtoAtomArray[i,1]:NAOtoAtomArray[i,2] ] = 0 EDDB_Alpha[ NAOtoAtomArray[i,1]:NAOtoAtomArray[i,2] , ] = 0 EDDB_Beta[ , NAOtoAtomArray[i,1]:NAOtoAtomArray[i,2] ] = 0 EDDB_Beta[ NAOtoAtomArray[i,1]:NAOtoAtomArray[i,2] , ] = 0 } } if (DEFAULT_EDDB_FUNCTION_TYPE %in% c( 4,40 ) ) { for (i in 1:Number_of_atoms) { if (!(i %in% SELECTED_Fragment)) { EDDB_Alpha[ , NAOtoAtomArray[i,1]:NAOtoAtomArray[i,2] ] = 0 EDDB_Alpha[ NAOtoAtomArray[i,1]:NAOtoAtomArray[i,2] , ] = 0 EDDB_Beta[ , NAOtoAtomArray[i,1]:NAOtoAtomArray[i,2] ] = 0 EDDB_Beta[ NAOtoAtomArray[i,1]:NAOtoAtomArray[i,2] , ] = 0 } } } EDDB_Total = EDDB_Alpha + EDDB_Beta sumEDDB_Alpha = sum( diag( EDDB_Alpha ) ) sumEDDB_Beta = sum( diag( EDDB_Beta ) ) sumEDDB_Total = sum( diag( EDDB_Total ) ) } # if ( DEFAULT_EDDB_FUNCTION_TYPE %in% c( 5 ) ) { frstat = ((1:Number_of_atoms)[ diag(Matrix_of_atomic_connectivity)>0 ])[1] MinEDDBPerAtom = sum ( diag( EDDB_Total[ NAOtoAtomArray[frstat,1]:NAOtoAtomArray[frstat,2], NAOtoAtomArray[frstat,1]:NAOtoAtomArray[frstat,2] ] ) ) for (i in ((1:Number_of_atoms)[ diag(Matrix_of_atomic_connectivity)>0 ]) ) { ithat = sum ( diag( EDDB_Total[ NAOtoAtomArray[i,1]:NAOtoAtomArray[i,2], NAOtoAtomArray[i,1]:NAOtoAtomArray[i,2] ] ) ) MinEDDBPerAtom = min( MinEDDBPerAtom, ithat ) } } if ( (!INTERACTIVE_MODE) || SKIP_EDDB__AND_PERFORM_NPA ) cat(" DONE!") if (exists("DSYEVR_Err_Counter") && ( DSYEVR_Err_Counter > 0 ) ) { cat(sprintf(paste0("\n WARNING: The DSYEVR subroutine returned %",as.integer(log10( DSYEVR_Err_Counter ))+1,"d errors during calculations!", collapse=""), DSYEVR_Err_Counter ) ) } # ____________________________________________________________________________ # if (!exists("NOBD_val")) { NOBD_val = eigen( EDDB_Total, symmetric=TRUE, only.values=TRUE )$values } cat(paste0("\n o List of non-zero eigenvalues of the ",EDDB_type,"(r) function:\n", collapse="") ) if (!SKIP_EDDB__AND_PERFORM_NPA) { cat("\n NOBD ") } else { cat("\n _NO_ ") } cat("___1__ ___2__ ___3__ ___4__ ___5__ ___6__ ___7__ ___8__ ___9__ __10__\n\n 1 ") NOBDocc_to_display = sum ( as.integer( NOBD_val >= 10*NATURAL_OCCUPATION_THRSHLD ) ) for (i in 1: NOBDocc_to_display) { cat(sprintf("%6.4f ", NOBD_val[i])) if ((i%%10==0)&&(i0) cat("\n") # if (isSpinUnrestricted && ( SAVE_SPIN_EDDBS_SEPARATELY || NBODensType == 2 ) ) { NOBD_Alpha_val = eigen(EDDB_Alpha,symmetric=TRUE,only.values=TRUE)$values NOBD_Beta_val = eigen(EDDB_Beta, symmetric=TRUE,only.values=TRUE)$values NOBDocc_Alpha_to_display = sum ( as.integer( NOBD_Alpha_val >= 10*NATURAL_OCCUPATION_THRSHLD ) ) NOBDocc_Beta_to_display = sum ( as.integer( NOBD_Beta_val >= 10*NATURAL_OCCUPATION_THRSHLD ) ) cat("\n Alpha ___1__ ___2__ ___3__ ___4__ ___5__ ___6__ ___7__ ___8__ ___9__ __10__\n\n 1 ") for (i in 1: NOBDocc_Alpha_to_display) { cat(sprintf("%6.4f ", NOBD_Alpha_val[i])) if ((i%%10==0)&&(i0) cat("\n") cat("\n Beta ___1__ ___2__ ___3__ ___4__ ___5__ ___6__ ___7__ ___8__ ___9__ __10__\n\n 1 ") for (i in 1: NOBDocc_Beta_to_display) { cat(sprintf("%6.4f ", NOBD_Beta_val[i])) if ((i%%10==0)&&(i0) cat("\n") } # cat( paste0(" ", paste0(rep("_",74),collapse=""),"\n", collapse="")) # ____________________________________________________________________________ # if (isSpinUnrestricted && SAVE_SPIN_EDDBS_SEPARATELY) { # ABLNOBD-dissection module need to be rewritten to deal with option -oo ATOMICORBIT_DECOMP_OF_EDDB = FALSE } # if (SAVE_EDDB_TO_GAUSSIAN_FCHK&&(!ORBITAL_DISSECTION_OF_EDDB)&&(!ATOMICORBIT_DECOMP_OF_EDDB)) { if (!exists("SAVE_EDDB_TO_FCHK_FILENAME")) SAVE_EDDB_TO_FCHK_FILENAME = paste0(EDDB_type,".fchk",collapse="") if (!SKIP_EDDB__AND_PERFORM_NPA) cat( paste0("\n o Saving ", EDDB_type , "(r) and its eigenvectors (NOBDs) to ", SAVE_EDDB_TO_FCHK_FILENAME, "...", collapse="") ) else cat( paste0("\n o Saving ED(r) and its eigenvectors (NOs) to ", SAVE_EDDB_TO_FCHK_FILENAME, "...", collapse="") ) NOBD_vec = eigen( EDDB_Total, symmetric=TRUE ) NOBD_val = NOBD_vec$values NOBD_vec = NOBD_vec$vectors FCHK_output = file(SAVE_EDDB_TO_FCHK_FILENAME, open="wt") if (!isSpinUnrestricted) Save_EDDB_to_FCHK( FCHK_content, FCHK_output, EDDB_Total, NA, NOBD_val, NOBD_OCCUPATIONS_SCAL_FACT, NOBD_vec, AOtoNAO_matrix, Number_of_basis_functions, Number_of_independent_functions, SAVE_SPIN_EDDBS_SEPARATELY ) else Save_EDDB_to_FCHK( FCHK_content, FCHK_output, EDDB_Total, EDDB_Alpha-EDDB_Beta, NOBD_val, NOBD_OCCUPATIONS_SCAL_FACT, NOBD_vec, AOtoNAO_matrix, Number_of_basis_functions, Number_of_independent_functions, SAVE_SPIN_EDDBS_SEPARATELY ) suppressWarnings(close(FCHK_output)) cat(" DONE!") } # ____________________________________________________________________________ # OrbitalDissectionAborted = FALSE if (ORBITAL_DISSECTION_OF_EDDB) { if (!SKIP_EDDB__AND_PERFORM_NPA) cat("\n o Enter the list of spinless NOBDs to dissect and/or press ") else cat("\n o Enter the list of spinless NOs to dissect and/or press ") ContinueLoop = TRUE LoopCounter = 0 while( ContinueLoop ) { ContinueLoop = FALSE LoopCounter = LoopCounter + 1 cat("\n > ") if ((LoopCounter==1)&&(exists("SELECTED_NOBD_List"))) cat(paste0(SELECTED_NOBD_List,"\n",collapse="")) if (((LoopCounter==1)&&(!exists("SELECTED_NOBD_List"))&&(!INTERACTIVE_MODE))||((LoopCounter>1)&&(!INTERACTIVE_MODE))) { cat("\n Interactive mode disabled! Script terminated.\n\n") ; close(STD_IN) ; quit(save="no") } if ((LoopCounter>1)||(!exists("SELECTED_NOBD_List"))) SELECTED_NOBD_List = suppressWarnings( readLines(STD_IN,n=1) ) SELECTED_NOBD_List_splited = strsplit( SELECTED_NOBD_List, ',', fixed=TRUE)[[1]] if (length(SELECTED_NOBD_List_splited)<1) { ContinueLoop = FALSE OrbitalDissectionAborted = TRUE cat(" Orbital-dissection aborted!") } else { SELECTED_NOBD_List = rep( FALSE , Number_of_NAOs ) for (i in 1:length(SELECTED_NOBD_List_splited)) { if (!is.na(suppressWarnings( as.integer( SELECTED_NOBD_List_splited[i] )))) { if (as.integer( SELECTED_NOBD_List_splited[i] )<=Number_of_NAOs) SELECTED_NOBD_List[as.integer( SELECTED_NOBD_List_splited[i] )] = TRUE else ContinueLoop = TRUE } else { SplitedItemsd = strsplit( SELECTED_NOBD_List_splited[i] , ':', fixed=TRUE)[[1]] if (length(SplitedItemsd)!=2) { ContinueLoop = TRUE } else { if ( (is.na(suppressWarnings( as.integer( SplitedItemsd[1] )))) || (is.na(suppressWarnings( as.integer( SplitedItemsd[1] )))) ) { ContinueLoop = TRUE } else { if ( (as.integer( SplitedItemsd[1] ) <= Number_of_NAOs) && (as.integer( SplitedItemsd[2] ) <= Number_of_NAOs) ) SELECTED_NOBD_List[as.integer( SplitedItemsd[1] ):as.integer( SplitedItemsd[2] )] = TRUE else ContinueLoop = TRUE } } } } } } } if (ORBITAL_DISSECTION_OF_EDDB && (!OrbitalDissectionAborted) ) { cat(paste0(" o Calculation of the dissected ", EDDB_type, "(r) function...", collapse="") ) if (!exists("NOBD_vec")) { NOBD_vec = eigen( EDDB_Total, symmetric=TRUE ) NOBD_val = NOBD_vec$values NOBD_vec = NOBD_vec$vectors } NOBD_val = NOBD_val * as.integer( SELECTED_NOBD_List ) EDDB_Total = NOBD_vec %*% diag(NOBD_val,Number_of_NAOs) %*% t( NOBD_vec ) sumEDDB_Total = sum( diag( EDDB_Total ) ) Matrix_of_Alpha_Density = t(NOBD_vec) %*% Matrix_of_Alpha_Density %*% NOBD_vec for (i in 1:Number_of_NAOs) { if (!SELECTED_NOBD_List[i]) { Matrix_of_Alpha_Density[,i] = 0 ; Matrix_of_Alpha_Density[i,] = 0 } } Matrix_of_Alpha_Density = NOBD_vec %*% Matrix_of_Alpha_Density %*% t(NOBD_vec) if (isSpinUnrestricted) { EDDB_Alpha = t( NOBD_vec ) %*% EDDB_Alpha %*% NOBD_vec EDDB_Beta = t( NOBD_vec ) %*% EDDB_Beta %*% NOBD_vec for (i in 1:Number_of_NAOs) { if (!SELECTED_NOBD_List[i]) { EDDB_Alpha[,i] = 0 ; EDDB_Alpha[i,] = 0 EDDB_Beta[,i] = 0 ; EDDB_Beta[i,] = 0 } } EDDB_Alpha = NOBD_vec %*% EDDB_Alpha %*% t(NOBD_vec) EDDB_Beta = NOBD_vec %*% EDDB_Beta %*% t(NOBD_vec) sumEDDB_Alpha = sum ( diag( EDDB_Alpha ) ) sumEDDB_Beta = sum ( diag( EDDB_Beta ) ) Matrix_of_Beta_Density = t(NOBD_vec) %*% Matrix_of_Beta_Density %*% NOBD_vec for (i in 1:Number_of_NAOs) { if (!SELECTED_NOBD_List[i]) { Matrix_of_Beta_Density[,i] = 0 ; Matrix_of_Beta_Density[i,] = 0 } } Matrix_of_Beta_Density = NOBD_vec %*% Matrix_of_Beta_Density %*% t(NOBD_vec) } cat(" DONE!") } if (ORBITAL_DISSECTION_OF_EDDB && SAVE_EDDB_TO_GAUSSIAN_FCHK && (!ATOMICORBIT_DECOMP_OF_EDDB)) { if (!exists("NOBD_vec")) { NOBD_vec = eigen( EDDB_Total, symmetric=TRUE ) NOBD_val = NOBD_vec$values NOBD_vec = NOBD_vec$vectors } if (!exists("SAVE_EDDB_TO_FCHK_FILENAME")) SAVE_EDDB_TO_FCHK_FILENAME = paste0(EDDB_type, ".fchk", collapse="") if (OrbitalDissectionAborted) cat( paste0("\n o Saving ", EDDB_type , "(r) function to ", SAVE_EDDB_TO_FCHK_FILENAME, "...", collapse="") ) else cat( paste0("\n o Saving dissected ", EDDB_type , "(r) function to ", SAVE_EDDB_TO_FCHK_FILENAME, "...", collapse="") ) FCHK_output = file(SAVE_EDDB_TO_FCHK_FILENAME, open="wt") if (!isSpinUnrestricted) Save_EDDB_to_FCHK( FCHK_content, FCHK_output, EDDB_Total, NA, NOBD_val, NOBD_OCCUPATIONS_SCAL_FACT, NOBD_vec, AOtoNAO_matrix, Number_of_basis_functions, Number_of_independent_functions, SAVE_SPIN_EDDBS_SEPARATELY ) else Save_EDDB_to_FCHK( FCHK_content, FCHK_output, EDDB_Total, EDDB_Alpha-EDDB_Beta, NOBD_val, NOBD_OCCUPATIONS_SCAL_FACT, NOBD_vec, AOtoNAO_matrix, Number_of_basis_functions, Number_of_independent_functions, SAVE_SPIN_EDDBS_SEPARATELY ) suppressWarnings(close(FCHK_output)) # ____________________________________________________________________________ # # Save the dissected full electron density into new EDED.fchk # ____________________________________________________________________________ # if (!is.na( match("-od", tolower(RscriptARGS) ) ) ) { if (isSpinUnrestricted) EDED_Total = Matrix_of_Alpha_Density + Matrix_of_Beta_Density else EDED_Total = Matrix_of_Alpha_Density if (!exists("EDED_vec")) { EDED_vec = eigen( EDED_Total, symmetric=TRUE ) EDED_val = EDED_vec$values EDED_vec = EDED_vec$vectors } FCHK_output2 = file("EDED.fchk", open="wt") if (!isSpinUnrestricted) Save_EDDB_to_FCHK( FCHK_content, FCHK_output, EDED_Total, NA, EDED_val, NOBD_OCCUPATIONS_SCAL_FACT, EDED_vec, AOtoNAO_matrix, Number_of_basis_functions, Number_of_independent_functions, SAVE_SPIN_EDDBS_SEPARATELY ) else Save_EDDB_to_FCHK( FCHK_content, FCHK_output, EDED_Total, Matrix_of_Alpha_Density-Matrix_of_Beta_Density, EDED_val, NOBD_OCCUPATIONS_SCAL_FACT, EDED_vec, AOtoNAO_matrix, Number_of_basis_functions, Number_of_independent_functions, SAVE_SPIN_EDDBS_SEPARATELY ) suppressWarnings(close(FCHK_output2)) } # ____________________________________________________________________________ # cat(" DONE!") } if ( DEFAULT_EDDB_FUNCTION_TYPE %in% c( 5 ) ) { frstat = ((1:Number_of_atoms)[ diag(Matrix_of_atomic_connectivity)>0 ])[1] MinEDDBPerAtom = sum ( diag( EDDB_Total[ NAOtoAtomArray[frstat,1]:NAOtoAtomArray[frstat,2], NAOtoAtomArray[frstat,1]:NAOtoAtomArray[frstat,2] ] ) ) for (i in ((1:Number_of_atoms)[ diag(Matrix_of_atomic_connectivity)>0 ]) ) { ithat = sum ( diag( EDDB_Total[ NAOtoAtomArray[i,1]:NAOtoAtomArray[i,2], NAOtoAtomArray[i,1]:NAOtoAtomArray[i,2] ] ) ) MinEDDBPerAtom = min( MinEDDBPerAtom, ithat ) } } # ____________________________________________________________________________ # if (ATOMICORBIT_DECOMP_OF_EDDB) { NOBD_val = rep(0, Number_of_NAOs) NOBD_vec = matrix(0, Number_of_NAOs, Number_of_NAOs) if (DEFAULT_EDDB_FUNCTION_TYPE %in% c( 0,1,10 )) AtomsToPrint = 1:Number_of_atoms if (DEFAULT_EDDB_FUNCTION_TYPE %in% c( 2,20 )) { AtomsToPrint = 1:Number_of_atoms for (i in 1:Number_of_atoms) if (Array_of_atomic_numbers[i]==1) AtomsToPrint[i] = 0 AtomsToPrint = AtomsToPrint[ AtomsToPrint > 0 ] } if (DEFAULT_EDDB_FUNCTION_TYPE %in% c( 3,30,4,40 )) AtomsToPrint = SELECTED_Fragment if (DEFAULT_EDDB_FUNCTION_TYPE==5) AtomsToPrint = (1:Number_of_atoms)[ diag(Matrix_of_atomic_connectivity)>0 ] NumAtomsInc = length(AtomsToPrint) cat(paste0("\n o List of non-zero eigenvalues of the atomic-block localized ",EDDB_type,"(r):",collapse="")) if (!SKIP_EDDB__AND_PERFORM_NPA) cat("\n\n _Atom_ ABLNOBDsIdx ") else cat("\n\n _Atom_ _NOs_index_ ") cat("___1__ ___2__ ___3__ ___4__ ___5__ ___6__ ___7__ ___8__\n\n") for (i in 1:Number_of_atoms) if (i %in% AtomsToPrint) { NOBD_vec_i = eigen(EDDB_Total[ NAOtoAtomArray[i,1]:NAOtoAtomArray[i,2], NAOtoAtomArray[i,1]:NAOtoAtomArray[i,2] ]) NOBD_val_i = NOBD_vec_i$values NOBD_vec_i = NOBD_vec_i$vectors cat(sprintf( " %-4d%s %-5d-%-5d", i, PERIODIC_TABLE[ Array_of_atomic_numbers[i] ], NAOtoAtomArray[i,1], NAOtoAtomArray[i,2] ) ) for ( j in 1:min(8, length( NOBD_val_i )) ) if (NOBD_val_i[j] >= NATURAL_OCCUPATION_THRSHLD) cat(sprintf( "%7.4f", NOBD_val_i[j] ) ) cat("\n") NOBD_val[ NAOtoAtomArray[i,1]:NAOtoAtomArray[i,2] ] = NOBD_val_i NOBD_vec[ NAOtoAtomArray[i,1]:NAOtoAtomArray[i,2], NAOtoAtomArray[i,1]:NAOtoAtomArray[i,2]] = NOBD_vec_i } cat( paste0(" ", paste0(rep("_",74),collapse=""),"\n", collapse="")) } if (ATOMICORBIT_DECOMP_OF_EDDB && SAVE_EDDB_TO_GAUSSIAN_FCHK) { SAVE_SPIN_EDDBS_SEPARATELY = FALSE # This will be changed after reimplementation of this module if (!exists("SAVE_EDDB_TO_FCHK_FILENAME")) SAVE_EDDB_TO_FCHK_FILENAME = paste0(EDDB_type, ".fchk", collapse="") cat( paste0("\n o Saving ", EDDB_type , "(r) and its ABL-NOBDs to ", SAVE_EDDB_TO_FCHK_FILENAME, "...", collapse="") ) FCHK_output = file(SAVE_EDDB_TO_FCHK_FILENAME, open="wt") if (!isSpinUnrestricted) Save_EDDB_to_FCHK( FCHK_content, FCHK_output, EDDB_Total, NA, NOBD_val, NOBD_OCCUPATIONS_SCAL_FACT, NOBD_vec, AOtoNAO_matrix, Number_of_basis_functions, Number_of_independent_functions, SAVE_SPIN_EDDBS_SEPARATELY ) else Save_EDDB_to_FCHK( FCHK_content, FCHK_output, EDDB_Total, EDDB_Alpha-EDDB_Beta, NOBD_val, NOBD_OCCUPATIONS_SCAL_FACT, NOBD_vec, AOtoNAO_matrix, Number_of_basis_functions, Number_of_independent_functions, SAVE_SPIN_EDDBS_SEPARATELY ) suppressWarnings(close(FCHK_output)) cat(" DONE!") } # ____________________________________________________________________________ # cat(paste0("\n o Results of the NPA and ", EDDB_type," population analyses in the ",EDDB_basis," basis:\n", collapse="") ) cat( paste0(" ", paste0(rep("_",74),collapse="") , "\n", collapse="") ) if (DEFAULT_EDDB_FUNCTION_TYPE %in% c( 0,1,10 )) AtomsToPrint = 1:Number_of_atoms if (DEFAULT_EDDB_FUNCTION_TYPE %in% c( 2,20 )) { AtomsToPrint = 1:Number_of_atoms for (i in 1:Number_of_atoms) if (Array_of_atomic_numbers[i]==1) AtomsToPrint[i] = 0 AtomsToPrint = AtomsToPrint[ AtomsToPrint > 0 ] } if (DEFAULT_EDDB_FUNCTION_TYPE %in% c( 3,30,4,40 )) AtomsToPrint = SELECTED_Fragment if (DEFAULT_EDDB_FUNCTION_TYPE==5) AtomsToPrint = (1:Number_of_atoms)[ diag(Matrix_of_atomic_connectivity)>0 ] NumAtomsInc = length(AtomsToPrint) for (i in 1:Number_of_atoms ) if (!(i %in% AtomsToPrint)) { Matrix_of_Alpha_Density[NAOtoAtomArray[i,1]:NAOtoAtomArray[i,2],] = 0 Matrix_of_Alpha_Density[,NAOtoAtomArray[i,1]:NAOtoAtomArray[i,2]] = 0 } if (isSpinUnrestricted) { for (i in 1:Number_of_atoms ) if (!(i %in% AtomsToPrint)) { Matrix_of_Beta_Density[NAOtoAtomArray[i,1]:NAOtoAtomArray[i,2],] = 0 Matrix_of_Beta_Density[,NAOtoAtomArray[i,1]:NAOtoAtomArray[i,2]] = 0 } } sumED_Total = sum( diag( Matrix_of_Alpha_Density ) ) if (isSpinUnrestricted) { sumED_Alpha = sumED_Total sumED_Beta = sum( diag( Matrix_of_Beta_Density ) ) sumED_Total = sumED_Alpha + sumED_Beta } if (isSpinUnrestricted) { if (REMOVE_UHF_SPIN_CONTAMINAT) { cat(" | ::Total:: | ::Paired Electrons:: | ::Unpaired Electrons:: |\n") } else { cat(" | ::Total:: | ::Alpha Spin:: | ::Beta Spin:: |\n") } } else { cat(" | | | |\n") } cat( paste0(" | Atom NPA ",EDDB_type," | Atom NPA ",EDDB_type," | Atom NPA ",EDDB_type," |\n",collapse="" ) ) cat( paste0(" | '''''' '''''' '''''' | '''''' '''''' '''''' | '''''' '''''' '''''' |\n",collapse="" ) ) if (isSpinUnrestricted) { if (REMOVE_UHF_SPIN_CONTAMINAT) { for (i in 1:Number_of_atoms ) if (i %in% AtomsToPrint) { cat(sprintf(" | %-4d%s%8.4f%8.4f | %-4d%s%8.4f%8.4f | %-4d%s%8.4f%8.4f |\n", i , PERIODIC_TABLE[ Array_of_atomic_numbers[i] ] , sum( diag(Matrix_of_Alpha_Density)[ NAOtoAtomArray[i,1]:NAOtoAtomArray[i,2] ] ) + sum( diag(Matrix_of_Beta_Density)[ NAOtoAtomArray[i,1]:NAOtoAtomArray[i,2] ] ), sum( diag(EDDB_Total)[ NAOtoAtomArray[i,1]:NAOtoAtomArray[i,2] ] ), i , PERIODIC_TABLE[ Array_of_atomic_numbers[i] ], 2*sum( diag(Matrix_of_Beta_Density)[ NAOtoAtomArray[i,1]:NAOtoAtomArray[i,2] ] ), 2*sum( diag(EDDB_Beta)[ NAOtoAtomArray[i,1]:NAOtoAtomArray[i,2] ] ), i , PERIODIC_TABLE[ Array_of_atomic_numbers[i] ], sum( diag(Matrix_of_Alpha_Density)[ NAOtoAtomArray[i,1]:NAOtoAtomArray[i,2] ] ) - sum( diag(Matrix_of_Beta_Density)[ NAOtoAtomArray[i,1]:NAOtoAtomArray[i,2] ] ), sum( diag(EDDB_Alpha)[ NAOtoAtomArray[i,1]:NAOtoAtomArray[i,2] ] ) - sum( diag(EDDB_Beta)[ NAOtoAtomArray[i,1]:NAOtoAtomArray[i,2] ] ) ) ) } } else { for (i in 1:Number_of_atoms ) if (i %in% AtomsToPrint) { cat(sprintf(" | %-4d%s%8.4f%8.4f | %-4d%s%8.4f%8.4f | %-4d%s%8.4f%8.4f |\n", i , PERIODIC_TABLE[ Array_of_atomic_numbers[i] ] , sum( diag(Matrix_of_Alpha_Density)[ NAOtoAtomArray[i,1]:NAOtoAtomArray[i,2] ] ) + sum( diag(Matrix_of_Beta_Density)[ NAOtoAtomArray[i,1]:NAOtoAtomArray[i,2] ] ), sum( diag(EDDB_Total)[ NAOtoAtomArray[i,1]:NAOtoAtomArray[i,2] ] ), i , PERIODIC_TABLE[ Array_of_atomic_numbers[i] ], sum( diag(Matrix_of_Alpha_Density)[ NAOtoAtomArray[i,1]:NAOtoAtomArray[i,2] ] ), sum( diag(EDDB_Alpha)[ NAOtoAtomArray[i,1]:NAOtoAtomArray[i,2] ] ), i , PERIODIC_TABLE[ Array_of_atomic_numbers[i] ], sum( diag(Matrix_of_Beta_Density)[ NAOtoAtomArray[i,1]:NAOtoAtomArray[i,2] ] ), sum( diag(EDDB_Beta)[ NAOtoAtomArray[i,1]:NAOtoAtomArray[i,2] ] ) ) ) } } cat(" |") } else { cat(" |") ix = 1 for (i in 1:Number_of_atoms ) if (i %in% AtomsToPrint) { cat(sprintf(" %-4d%s%8.4f%8.4f |", i , PERIODIC_TABLE[ Array_of_atomic_numbers[i] ] , sum( diag(Matrix_of_Alpha_Density)[ NAOtoAtomArray[i,1]:NAOtoAtomArray[i,2] ] ), sum( diag(EDDB_Total)[ NAOtoAtomArray[i,1]:NAOtoAtomArray[i,2] ] ) ) ) if (ix%%3==0) cat("\n |") ix = ix + 1 } if ( (ix-1)%%3 > 0 ) { for (i in 1:(3-(ix-1)%%3) ) cat(" |") cat("\n |") } } cat( paste0(paste0(rep("_",24),collapse="") , "|", paste0(rep("_",24),collapse="") , "|", paste0(rep("_",24),collapse="") , "|\n", collapse="") ) cat( paste0(" |",paste0(rep(" ",74),collapse=""), "|\n", collapse="") ) cat(sprintf(paste0(" | Total population of electrons (from NPA):%10.4fe (%7.4fe / atom ) |\n", collapse=""), sumED_Total , sumED_Total / NumAtomsInc ) ) if (isSpinUnrestricted && !REMOVE_UHF_SPIN_CONTAMINAT) { cat(sprintf(" | + Alpha spin:%10.4fe (%7.4fe / atom ) |\n", sumED_Alpha , sumED_Alpha / NumAtomsInc ) ) cat(sprintf(" | + Beta spin:%10.4fe (%7.4fe / atom ) |\n", sumED_Beta , sumED_Beta / NumAtomsInc ) ) } if (isSpinUnrestricted && REMOVE_UHF_SPIN_CONTAMINAT) { cat(sprintf(" | + Paired electrons:%10.4fe (%7.4fe / atom ) |\n", 2*sumED_Beta , 2*sumED_Beta / NumAtomsInc ) ) cat(sprintf(" | + Unpaired electrons:%10.4fe (%7.4fe / atom ) |\n", sumED_Alpha-sumED_Beta , (sumED_Alpha-sumED_Beta) / NumAtomsInc ) ) } if (!SKIP_EDDB__AND_PERFORM_NPA) { cat(sprintf(" | Total population of delocalized electrons:%10.4fe (%7.4fe / atom ) |\n", sumEDDB_Total , sumEDDB_Total / NumAtomsInc ) ) if (isSpinUnrestricted && !REMOVE_UHF_SPIN_CONTAMINAT) { cat(sprintf(" | + Alpha spin:%10.4fe (%7.4fe / atom ) |\n", sumEDDB_Alpha , sumEDDB_Alpha / NumAtomsInc ) ) cat(sprintf(" | + Beta spin:%10.4fe (%7.4fe / atom ) |\n", sumEDDB_Beta , sumEDDB_Beta / NumAtomsInc ) ) } if (isSpinUnrestricted && REMOVE_UHF_SPIN_CONTAMINAT) { cat(sprintf(" | + Contribution from paired electrons:%10.4fe (%7.4fe / atom ) |\n", 2*sumEDDB_Beta , 2*sumEDDB_Beta / NumAtomsInc ) ) cat(sprintf(" | + Contribution from unpaired electrons:%10.4fe (%7.4fe / atom ) |\n", sumEDDB_Alpha-sumEDDB_Beta , (sumEDDB_Alpha-sumEDDB_Beta) / NumAtomsInc ) ) } if ( DEFAULT_EDDB_FUNCTION_TYPE %in% c( 5 ) ) { cat(sprintf(" | o Cyclically delocalized electrons:%10.4fe (%7.4fe / atom ) |\n", MinEDDBPerAtom*NumAtomsInc , MinEDDBPerAtom ) ) } } cat( paste0(" |", paste0(rep("_",74),collapse="") , "|\n", collapse="") ) # cat(paste0("\n o Normal termination of RunEDDB at ",format(Sys.time(), "%a %b %d %X %Y"),".\n", collapse="" ) ) if (INTERACTIVE_MODE) { cat("\n Press to quit the script.\n", collapse="") ; nullek=readLines(STD_IN,n=1) } else { cat("\n\n" ) } close(STD_IN) quit(save="no") # ____________________________________________________________________________ # # ###########################NNW################WNW########################### # ###########################WNXXNWW######WWWWWXXNW########################### # #############################KxkKNNNWWWN0kOKNXN############################# # #############################NxcldkO000Ol;;lkXW############################# # #############################WOc:lkXWWWW0l;;oK############################## # ###########################WWWKkOKN######NKkxKWW############################ # ########################WNXXXNNNW#########N0kOXWW###########WW############## # ##############NXNNNNWWWNKOkO0XWWWWW#####WXOxxk0XNWWWWNNNWWWNXN############## # ##############W0xOKXWWWXOdoodxOKNNNNNNNNNKOxdddxkOOOOkO0XNWWW############### # ##############XxoxO0KK0ko:,,,;:ldxkOKXNNWWXOdc:::;;;;:ldk0NWWW############## # #############W0oodxxkxdoolcc::;,,;coxO0KKKKkocccll;'##';lxOXWW############## # #############NOxxxxxolox0XNNNNKOl,,;coodxkkOKXNNWNKxc'##,ckKWW############## # ###########WX0O0XXK0kxOXW#######Nkc,,;;:lx0N#######WXx:##:d0NW############## # WW########N0xkOKNNNXXNW##########WX0kxxk0XW##########Nx,';lkKNW##########WWW # WNXXNNWWWNkoloxO0XNNXNW##############################Wx,';coxOKXNNNXXNWWNXXW # #W0k0XWWN0o;,;coxOXNNNW#############################WXxccllooddddoodxOKNWWW# # #Nkdk0KK0kkkkkdlld0XNNWWW#########################NK0OO0KK0kxkOkc'#';lxKNW## # #W0ddkkkxxOXW#WXdlxOXWWWWWW#####################W0oldk0XWWWNNW#WKo,##;okKN## # ##Nkk0K00KXW###WkloxOXNNXXKXNW###############WNNKo,,cokKNNNW#####0l;';dKNW## # ##XOKNWNNWW####NxccodkOOOOKNWW################WWXx;',:lxOKNW####NOo::o0NW### # #NOx0XWWNNNNWNNKkdollloxOKW######################N0d:;;:oOXWWWWNkc;;cdOXWW## # #Nkodk0KXXKOOO0XNNKOxxOXW##########################WXkc';oOKNWWN0koccoxOKNW# # WN0xdxxkOOdodx0XWWNXKNW##############################Wx,':okKNWWWXOdoodkOKNW # NNNWWWNWWWKdloxOXNNXKNW###########WWNWWWW############Wk,#,cdk0XXNWWWWNNWWNNN # ###########NkooxOXNNNNW#########WX0O0XXNNWW#########WXx;,:odkO0XW########### # ############W0ddkXWWNNNWW######NOddk0XNNWWWWW#####WXOdllox0KKKNW############ # #############WOox0NWWWNNNWNXXX0xoloxkKXXWWWWWNNNNX0occoxOKNWWW############## # ##############0oox0XWWWWNXK0000Oxooodxkk0KNNNNNNXXOxoodxOKNNNW############## # ##############Nxcldk0KKXK0O0KXNX0koloooooodkO0XNNNXOxoodxkO0XW############## # ##############N0kdddxxxxdddxk0KK000KXXXXKklcldk0KXKOxdoddkO0XW############## # #############WNWWWWNNNNX0kdoodxOKN#######WKd:codkOKXXNNNNWWWNW############## # #########################WNK0OxkX#########WKddk0XNW######################### # ##############################Kx0NW#####NKKkkX############################## # #############################NOkKNNXKXXKxok0XW############################## # #############################Xxx0XNXKKK0xokKXNW############################# # ############################WKOk0KNW####NK00KKN############################# # ###########################NXNW##############WNNW########################### # ##########################WNN#################WNW###########################