diff --git a/Manifest.toml b/Manifest.toml index 3f20f75..c96f130 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -16,7 +16,7 @@ version = "1.2.1" [[deps.ActiveSpaceSolvers]] deps = ["Arpack", "BenchmarkTools", "BlockDavidson", "Documenter", "InCoreIntegrals", "InteractiveUtils", "JLD2", "KrylovKit", "LinearAlgebra", "LinearMaps", "NPZ", "OrderedCollections", "Parameters", "Printf", "Profile", "QCBase", "StaticArrays", "TensorOperations", "TimerOutputs"] -git-tree-sha1 = "f781036635964dfdd47b61e04473a028e2f22904" +git-tree-sha1 = "2e7c4864a25152ea84b2d217e7e8d0339f55b2f9" repo-rev = "main" repo-url = "https://github.com/nmayhall-vt/ActiveSpaceSolvers.jl.git" uuid = "f8e94ea7-44cc-4f7e-92b7-2a8da6f47aca" diff --git a/README.md b/README.md index 66669f4..a440639 100644 --- a/README.md +++ b/README.md @@ -17,15 +17,16 @@ Perform `CMF` (Cluster Mean-Field) calculations. This is simply a variational op 2. Create conda environment to install Julia and will hold the PySCF executable. Install Julia with conda makes sure the correct python version will be found when using PyCall. where `-tauto` let's Julia pick the max number of threads. Use `-t N` to select `N` manually. Removing defaults to 1 thread. ```bash - conda create -n my_env python=3.7 - conda activate my_env + conda create -n my_env + conda activate my_env + conda install python==3.7 + conda config --add channels conda-forge + conda install -c pyscf pyscf + conda install h5py==2.10.0 conda install julia - conda install numpy - pip install pyscf - julia --project=./ -tauto + julia --project=./ -tauto ``` - 3. Build PyCall from Julia REPL ```julia @@ -39,7 +40,7 @@ Perform `CMF` (Cluster Mean-Field) calculations. This is simply a variational op ``` -## Installation with Conda on an Apple M1 chip +## Installation with Conda on an Apple M1 or M2 chip 1. Download ```bash @@ -51,12 +52,14 @@ Perform `CMF` (Cluster Mean-Field) calculations. This is simply a variational op 2. Create conda environment to install Julia and will hold the PySCF executable. Install Julia with conda makes sure the correct python version will be found when using PyCall. where `-tauto` let's Julia pick the max number of threads. Use `-t N` to select `N` manually. Removing defaults to 1 thread. ```bash - CONDA_SUBDIR=osx-64 conda create -n myenv_x86 python=3.7 - conda activate my_env_X86 - conda config --env --set subdir osx-64 + conda create -n env_osx + conda activate env_osx + conda config --env --set subdir osx-64 + conda install python==3.7 + conda config --add channels conda-forge + conda install -c pyscf pyscf + conda install h5py==2.10.0 conda install julia - conda install numpy - pip install pyscf julia --project=./ -tauto ``` diff --git a/examples/rasci_cmf.jl b/examples/rasci_cmf.jl new file mode 100644 index 0000000..e75e9ee --- /dev/null +++ b/examples/rasci_cmf.jl @@ -0,0 +1,84 @@ +using QCBase +using ClusterMeanField +using ActiveSpaceSolvers +using RDM +using InCoreIntegrals +using PyCall +using JLD2 +using NPZ + +molecule = " +C -6.8604981940 -0.4376683883 0.0000000000 +C -7.0417938838 0.9473143671 0.0000000000 +C -5.5765818320 -0.9841959654 0.0000000000 +C -5.9160648559 1.7742824642 0.0000000000 +C -4.4255614400 -0.1700917199 0.0000000000 +C -4.6336867757 1.2242295397 0.0000000000 +C -3.0514687259 -0.7633813258 0.0000000000 +C -2.8520613255 -2.1570405955 0.0000000000 +C -1.5640941394 -2.6887327507 0.0000000000 +C -0.4451033014 -1.8584342989 0.0000000000 +C -0.5927492072 -0.4581981628 0.0000000000 +C -1.9041039939 0.0505338987 0.0000000000 +H -8.0463585128 1.3756739224 0.0000000000 +H -5.4821758879 -2.0696091478 0.0000000000 +H -6.0325833442 2.8604965477 0.0000000000 +H -3.7849879120 1.9076345469 0.0000000000 +H -3.6983932318 -2.8425398416 0.0000000000 +H -1.4298106521 -3.7726953568 0.0000000000 +H 0.5426199813 -2.3168156796 0.0000000000 +H -2.0368540411 1.1280906363 0.0000000000 +H -7.7260669052 -1.1042415323 0.0000000000 +H 0.22700 0.16873 0.00000 +" +atoms = [] +for (li,line) in enumerate(split(rstrip(lstrip(molecule)), "\n")) + l = split(line) + push!(atoms, Atom(li, l[1], parse.(Float64,l[2:4]))) +end + +basis = "cc-pvdz" +# Create FermiCG.Molecule type +pyscf = pyimport("pyscf") +tools = pyimport("pyscf.tools") +fcidump = pyimport("pyscf.tools.fcidump"); + +pymol = pyscf.gto.Mole(atom=molecule, + symmetry = false, spin =0,charge=0, + basis = basis) +pymol.build() + +h0 = npzread("examples/two_benzene/ints_h0.npy") +h1 = npzread("examples/two_benzene/ints_h1.npy") +h2 = npzread("examples/two_benzene/ints_h2.npy") +ints = InCoreInts(h0, h1, h2) + +Cact = npzread("examples/two_benzene/mo_coeffs.npy") + +clusters_in = [(1:6),(7:12)] +n_clusters = 2 +cluster_list = [collect(1:6), collect(7:12)] +clusters = [MOCluster(i,collect(cluster_list[i])) for i = 1:length(cluster_list)] +init_fspace = [ (3,3) for i in 1:n_clusters] +display(clusters) +ansatze = [RASCIAnsatz(6, 3, 3, (1,4,1), max_h=2, max_p=2), RASCIAnsatz(6,3,3,(1,4,1), max_h=2, max_p=2)] #FCI type RASCI calculation +#ansatze = [RASCIAnsatz(6, 3, 3, (1,4,1), max_h=1, max_p=1), RASCIAnsatz(6,3,3,(1,4,1), max_h=1, max_p=1)] #Single excitation RASCI calculation +#ansatze = [RASCIAnsatz(6, 3, 3, (1,4,1), max_h=0, max_p=0), RASCIAnsatz(6,3,3,(1,4,1), max_h=0, max_p=0)] #CAS type RASCI calculation +for i in ansatze + display(i) +end + + + +rdm1 = zeros(size(ints.h1)) +#run cmf_oo +e_cmf, U_cmf, d1 = ClusterMeanField.cmf_oo_diis(ints, clusters, init_fspace, ansatze, RDM1(rdm1, rdm1), maxiter_oo = 100, zero_intra_rots=true, orb_hessian=true, tol_oo=1e-6, tol_ci=1e-8, verbose=0, diis_start=1); +ints_cmf = orbital_rotation(ints,U_cmf) +C_cmf = Cact*U_cmf +pyscf.tools.molden.from_mo(pymol, "examples/two_benzene/C_cmf.molden", C_cmf) + + + + + + diff --git a/examples/two_benzene/generate_integrals.ipynb b/examples/two_benzene/generate_integrals.ipynb new file mode 100644 index 0000000..b5002c7 --- /dev/null +++ b/examples/two_benzene/generate_integrals.ipynb @@ -0,0 +1,941 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "import pyscf\n", + "import pyscf.tools\n", + "import sys\n", + "# git clone https://github.com/nmayhall-vt/OrbitalPartitioning.git\n", + "# give path to OrbitalPartitioning repo\n", + "sys.path.append('../../../OrbitalPartitioning/orbitalpartitioning')\n", + "import orbitalpartitioning\n", + "from orbitalpartitioning import *" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n", + "\n", + "******** ********\n", + "method = RHF\n", + "initial guess = minao\n", + "damping factor = 0\n", + "level_shift factor = 0\n", + "DIIS = \n", + "diis_start_cycle = 1\n", + "diis_space = 8\n", + "SCF conv_tol = 1e-08\n", + "SCF conv_tol_grad = None\n", + "SCF max_cycles = 200\n", + "direct_scf = True\n", + "direct_scf_tol = 1e-13\n", + "chkfile to save SCF result = /var/folders/td/qpnnxwv93pq0t7bbdkh5rvzr0000gn/T/tmpppy4megv\n", + "max_memory 4000 MB (current use 0 MB)\n", + "Set gradient conv threshold to 0.0001\n", + "init E= -463.937039728682\n", + " HOMO = -0.201399173349765 LUMO = 0.00065918833382099\n", + "cycle= 1 E= -460.148732039677 delta_E= 3.79 |g|= 0.563 |ddm|= 4.88\n", + " HOMO = -0.31292892155204 LUMO = 0.0615105158847815\n", + "cycle= 2 E= -460.269257271685 delta_E= -0.121 |g|= 0.16 |ddm|= 1\n", + " HOMO = -0.277271190926687 LUMO = 0.100559181707899\n", + "cycle= 3 E= -460.277736606206 delta_E= -0.00848 |g|= 0.0783 |ddm|= 0.408\n", + " HOMO = -0.289669476214775 LUMO = 0.0900201444921001\n", + "cycle= 4 E= -460.27951161518 delta_E= -0.00178 |g|= 0.0151 |ddm|= 0.0956\n", + " HOMO = -0.291035449002703 LUMO = 0.0890988356785322\n", + "cycle= 5 E= -460.279579645017 delta_E= -6.8e-05 |g|= 0.00361 |ddm|= 0.0194\n", + " HOMO = -0.290627999711558 LUMO = 0.0895776292836337\n", + "cycle= 6 E= -460.279583929374 delta_E= -4.28e-06 |g|= 0.000701 |ddm|= 0.00396\n", + " HOMO = -0.290610925686448 LUMO = 0.0897018966170447\n", + "cycle= 7 E= -460.279584235993 delta_E= -3.07e-07 |g|= 0.00021 |ddm|= 0.00156\n", + " HOMO = -0.290627453534956 LUMO = 0.0897139202990004\n", + "cycle= 8 E= -460.279584259332 delta_E= -2.33e-08 |g|= 8.31e-05 |ddm|= 0.000364\n", + " HOMO = -0.290633736370228 LUMO = 0.0897220087167466\n", + "cycle= 9 E= -460.279584263704 delta_E= -4.37e-09 |g|= 2.29e-05 |ddm|= 0.000154\n", + " HOMO = -0.290636046410558 LUMO = 0.0897215106909744\n", + "Extra cycle E= -460.279584264235 delta_E= -5.31e-10 |g|= 1.41e-05 |ddm|= 4.23e-05\n", + "converged SCF energy = -460.279584264235\n", + " Hartree-Fock Energy: -460.27958426\n", + "**** SCF Summaries ****\n", + "Total Energy = -460.279584264234813\n", + "Nuclear Repulsion Energy = 597.208587579417213\n", + "One-electron Energy = -1805.834403535370257\n", + "Two-electron Energy = 748.346231691718231\n", + "**** MO energy ****\n", + "MO #1 energy= -11.2548525821824 occ= 2\n", + "MO #2 energy= -11.2533807731765 occ= 2\n", + "MO #3 energy= -11.2426857474598 occ= 2\n", + "MO #4 energy= -11.2424198032002 occ= 2\n", + "MO #5 energy= -11.2421424115507 occ= 2\n", + "MO #6 energy= -11.2419920146235 occ= 2\n", + "MO #7 energy= -11.2415095063144 occ= 2\n", + "MO #8 energy= -11.241023644533 occ= 2\n", + "MO #9 energy= -11.2408813861744 occ= 2\n", + "MO #10 energy= -11.2402674567159 occ= 2\n", + "MO #11 energy= -11.2401525093658 occ= 2\n", + "MO #12 energy= -11.2397111488365 occ= 2\n", + "MO #13 energy= -1.1662290305918 occ= 2\n", + "MO #14 energy= -1.13874705111095 occ= 2\n", + "MO #15 energy= -1.05160453809425 occ= 2\n", + "MO #16 energy= -1.02050996310093 occ= 2\n", + "MO #17 energy= -1.00926616309708 occ= 2\n", + "MO #18 energy= -0.981837524791909 occ= 2\n", + "MO #19 energy= -0.851733508930868 occ= 2\n", + "MO #20 energy= -0.846261944215565 occ= 2\n", + "MO #21 energy= -0.804546417751293 occ= 2\n", + "MO #22 energy= -0.790244732794544 occ= 2\n", + "MO #23 energy= -0.720879440973737 occ= 2\n", + "MO #24 energy= -0.67452605097356 occ= 2\n", + "MO #25 energy= -0.662418902206137 occ= 2\n", + "MO #26 energy= -0.638761458792849 occ= 2\n", + "MO #27 energy= -0.615238076726895 occ= 2\n", + "MO #28 energy= -0.596277744386189 occ= 2\n", + "MO #29 energy= -0.594331107759636 occ= 2\n", + "MO #30 energy= -0.585593303656711 occ= 2\n", + "MO #31 energy= -0.572241207624882 occ= 2\n", + "MO #32 energy= -0.523195195047075 occ= 2\n", + "MO #33 energy= -0.513116069668836 occ= 2\n", + "MO #34 energy= -0.507663228907959 occ= 2\n", + "MO #35 energy= -0.486620291107741 occ= 2\n", + "MO #36 energy= -0.48652682347127 occ= 2\n", + "MO #37 energy= -0.472794803697839 occ= 2\n", + "MO #38 energy= -0.376654956388138 occ= 2\n", + "MO #39 energy= -0.339254537013095 occ= 2\n", + "MO #40 energy= -0.332259423685731 occ= 2\n", + "MO #41 energy= -0.290636046410558 occ= 2\n", + "MO #42 energy= 0.0897215106909744 occ= 0\n", + "MO #43 energy= 0.119355813473763 occ= 0\n", + "MO #44 energy= 0.150075001499915 occ= 0\n", + "MO #45 energy= 0.180205581212967 occ= 0\n", + "MO #46 energy= 0.184675594324371 occ= 0\n", + "MO #47 energy= 0.192968062349489 occ= 0\n", + "MO #48 energy= 0.20619177858781 occ= 0\n", + "MO #49 energy= 0.222834247370418 occ= 0\n", + "MO #50 energy= 0.229269882034943 occ= 0\n", + "MO #51 energy= 0.248937272179874 occ= 0\n", + "MO #52 energy= 0.251596095123004 occ= 0\n", + "MO #53 energy= 0.2622706979524 occ= 0\n", + "MO #54 energy= 0.323979302547739 occ= 0\n", + "MO #55 energy= 0.326484551438172 occ= 0\n", + "MO #56 energy= 0.348172291190935 occ= 0\n", + "MO #57 energy= 0.396235628300867 occ= 0\n", + "MO #58 energy= 0.407391046811851 occ= 0\n", + "MO #59 energy= 0.413061412092286 occ= 0\n", + "MO #60 energy= 0.421571135160137 occ= 0\n", + "MO #61 energy= 0.440522194698462 occ= 0\n", + "MO #62 energy= 0.460967067073428 occ= 0\n", + "MO #63 energy= 0.46483577368429 occ= 0\n", + "MO #64 energy= 0.508778229215486 occ= 0\n", + "MO #65 energy= 0.539349386935712 occ= 0\n", + "MO #66 energy= 0.560187015792421 occ= 0\n", + "MO #67 energy= 0.581403546614198 occ= 0\n", + "MO #68 energy= 0.59664729203523 occ= 0\n", + "MO #69 energy= 0.659069131290935 occ= 0\n", + "MO #70 energy= 0.660590258479025 occ= 0\n", + "MO #71 energy= 0.665198289033454 occ= 0\n", + "MO #72 energy= 0.689313305476685 occ= 0\n", + "MO #73 energy= 0.703831371349851 occ= 0\n", + "MO #74 energy= 0.706116042019784 occ= 0\n", + "MO #75 energy= 0.707897590698848 occ= 0\n", + "MO #76 energy= 0.715518293706654 occ= 0\n", + "MO #77 energy= 0.71771357027193 occ= 0\n", + "MO #78 energy= 0.720170933712935 occ= 0\n", + "MO #79 energy= 0.727447295543825 occ= 0\n", + "MO #80 energy= 0.730227216914237 occ= 0\n", + "MO #81 energy= 0.737474791985108 occ= 0\n", + "MO #82 energy= 0.743154512364923 occ= 0\n", + "MO #83 energy= 0.744848837477113 occ= 0\n", + "MO #84 energy= 0.756837902662189 occ= 0\n", + "MO #85 energy= 0.75841943929689 occ= 0\n", + "MO #86 energy= 0.770682243043897 occ= 0\n", + "MO #87 energy= 0.785901706072578 occ= 0\n", + "MO #88 energy= 0.791940952131178 occ= 0\n", + "MO #89 energy= 0.816166705225343 occ= 0\n", + "MO #90 energy= 0.83526223956226 occ= 0\n", + "MO #91 energy= 0.839468823687663 occ= 0\n", + "MO #92 energy= 0.868945887350696 occ= 0\n", + "MO #93 energy= 0.878390638020876 occ= 0\n", + "MO #94 energy= 0.881915218640657 occ= 0\n", + "MO #95 energy= 0.892318992485719 occ= 0\n", + "MO #96 energy= 0.905471726813278 occ= 0\n", + "MO #97 energy= 0.906293922480843 occ= 0\n", + "MO #98 energy= 0.915116535720594 occ= 0\n", + "MO #99 energy= 0.918044242276448 occ= 0\n", + "MO #100 energy= 0.9428311345866 occ= 0\n", + "MO #101 energy= 0.944714162263674 occ= 0\n", + "MO #102 energy= 0.953649302254394 occ= 0\n", + "MO #103 energy= 0.974807974328732 occ= 0\n", + "MO #104 energy= 1.03266123772735 occ= 0\n", + "MO #105 energy= 1.03808372050473 occ= 0\n", + "MO #106 energy= 1.09806604254795 occ= 0\n", + "MO #107 energy= 1.10887049247839 occ= 0\n", + "MO #108 energy= 1.11628821308336 occ= 0\n", + "MO #109 energy= 1.11787671116292 occ= 0\n", + "MO #110 energy= 1.14568395054946 occ= 0\n", + "MO #111 energy= 1.14653333370105 occ= 0\n", + "MO #112 energy= 1.14874815018169 occ= 0\n", + "MO #113 energy= 1.15417283505954 occ= 0\n", + "MO #114 energy= 1.16931205217559 occ= 0\n", + "MO #115 energy= 1.1774123465538 occ= 0\n", + "MO #116 energy= 1.18465433585382 occ= 0\n", + "MO #117 energy= 1.2164573838625 occ= 0\n", + "MO #118 energy= 1.22339874848031 occ= 0\n", + "MO #119 energy= 1.23805066479608 occ= 0\n", + "MO #120 energy= 1.24215827970562 occ= 0\n", + "MO #121 energy= 1.2457782298926 occ= 0\n", + "MO #122 energy= 1.26259652665412 occ= 0\n", + "MO #123 energy= 1.26323284263734 occ= 0\n", + "MO #124 energy= 1.28925895053329 occ= 0\n", + "MO #125 energy= 1.29073829964841 occ= 0\n", + "MO #126 energy= 1.30564756727339 occ= 0\n", + "MO #127 energy= 1.31159345663445 occ= 0\n", + "MO #128 energy= 1.32350625955934 occ= 0\n", + "MO #129 energy= 1.33627870960257 occ= 0\n", + "MO #130 energy= 1.34131998589766 occ= 0\n", + "MO #131 energy= 1.34897612103783 occ= 0\n", + "MO #132 energy= 1.43470591607807 occ= 0\n", + "MO #133 energy= 1.50415533300682 occ= 0\n", + "MO #134 energy= 1.54515152330658 occ= 0\n", + "MO #135 energy= 1.56758895429883 occ= 0\n", + "MO #136 energy= 1.58087642050707 occ= 0\n", + "MO #137 energy= 1.60917143138629 occ= 0\n", + "MO #138 energy= 1.62393703508318 occ= 0\n", + "MO #139 energy= 1.62918981029844 occ= 0\n", + "MO #140 energy= 1.7112570413825 occ= 0\n", + "MO #141 energy= 1.73620875890284 occ= 0\n", + "MO #142 energy= 1.74432400839607 occ= 0\n", + "MO #143 energy= 1.74523577885504 occ= 0\n", + "MO #144 energy= 1.74798363065883 occ= 0\n", + "MO #145 energy= 1.75424221208343 occ= 0\n", + "MO #146 energy= 1.76483854924502 occ= 0\n", + "MO #147 energy= 1.77047342704979 occ= 0\n", + "MO #148 energy= 1.78958484495769 occ= 0\n", + "MO #149 energy= 1.81054369778689 occ= 0\n", + "MO #150 energy= 1.87086783187752 occ= 0\n", + "MO #151 energy= 1.88021937737932 occ= 0\n", + "MO #152 energy= 1.8849118491094 occ= 0\n", + "MO #153 energy= 1.88540789887318 occ= 0\n", + "MO #154 energy= 1.89435270882091 occ= 0\n", + "MO #155 energy= 1.91293615783357 occ= 0\n", + "MO #156 energy= 1.91636477020422 occ= 0\n", + "MO #157 energy= 1.92149986447157 occ= 0\n", + "MO #158 energy= 1.92329694386752 occ= 0\n", + "MO #159 energy= 1.94067525403265 occ= 0\n", + "MO #160 energy= 1.94544356853578 occ= 0\n", + "MO #161 energy= 1.95350294222221 occ= 0\n", + "MO #162 energy= 1.97567624189273 occ= 0\n", + "MO #163 energy= 1.98885631319241 occ= 0\n", + "MO #164 energy= 2.00011875379592 occ= 0\n", + "MO #165 energy= 2.01127267675015 occ= 0\n", + "MO #166 energy= 2.02867619568614 occ= 0\n", + "MO #167 energy= 2.04332043522983 occ= 0\n", + "MO #168 energy= 2.04511544172026 occ= 0\n", + "MO #169 energy= 2.07512036731921 occ= 0\n", + "MO #170 energy= 2.11961525699477 occ= 0\n", + "MO #171 energy= 2.13621150931658 occ= 0\n", + "MO #172 energy= 2.13622282379066 occ= 0\n", + "MO #173 energy= 2.13875335078443 occ= 0\n", + "MO #174 energy= 2.14441938854006 occ= 0\n", + "MO #175 energy= 2.1453324341992 occ= 0\n", + "MO #176 energy= 2.15106359027083 occ= 0\n", + "MO #177 energy= 2.15612802827034 occ= 0\n", + "MO #178 energy= 2.16364451406117 occ= 0\n", + "MO #179 energy= 2.176366929546 occ= 0\n", + "MO #180 energy= 2.19566692706305 occ= 0\n", + "MO #181 energy= 2.22778740461821 occ= 0\n", + "MO #182 energy= 2.24633959689231 occ= 0\n", + "MO #183 energy= 2.27464630660617 occ= 0\n", + "MO #184 energy= 2.28617336713284 occ= 0\n", + "MO #185 energy= 2.29615669046041 occ= 0\n", + "MO #186 energy= 2.29625862443544 occ= 0\n", + "MO #187 energy= 2.29993180649203 occ= 0\n", + "MO #188 energy= 2.32000833664368 occ= 0\n", + "MO #189 energy= 2.35139371377032 occ= 0\n", + "MO #190 energy= 2.37658947674462 occ= 0\n", + "MO #191 energy= 2.40292072582343 occ= 0\n", + "MO #192 energy= 2.44267673036076 occ= 0\n", + "MO #193 energy= 2.48630309189178 occ= 0\n", + "MO #194 energy= 2.55009101482822 occ= 0\n", + "MO #195 energy= 2.57528860569754 occ= 0\n", + "MO #196 energy= 2.62018065240959 occ= 0\n", + "MO #197 energy= 2.63016811304037 occ= 0\n", + "MO #198 energy= 2.70261542643907 occ= 0\n", + "MO #199 energy= 2.7244925665085 occ= 0\n", + "MO #200 energy= 2.74061592444683 occ= 0\n", + "MO #201 energy= 2.77332800647077 occ= 0\n", + "MO #202 energy= 2.78510115043655 occ= 0\n", + "MO #203 energy= 2.79960560797016 occ= 0\n", + "MO #204 energy= 2.81378318576405 occ= 0\n", + "MO #205 energy= 2.83974636791895 occ= 0\n", + "MO #206 energy= 2.92104135473924 occ= 0\n", + "MO #207 energy= 2.98367645289646 occ= 0\n", + "MO #208 energy= 3.01199364577017 occ= 0\n", + "MO #209 energy= 3.03142999676666 occ= 0\n", + "MO #210 energy= 3.05137330689583 occ= 0\n", + "MO #211 energy= 3.0720549239942 occ= 0\n", + "MO #212 energy= 3.09802112817749 occ= 0\n", + "MO #213 energy= 3.24045891567354 occ= 0\n", + "MO #214 energy= 3.24266879464421 occ= 0\n", + "MO #215 energy= 3.28967995268323 occ= 0\n", + "MO #216 energy= 3.39744928359553 occ= 0\n", + "MO #217 energy= 3.84876104575469 occ= 0\n", + "MO #218 energy= 4.01926738279133 occ= 0\n", + " ** Mulliken pop on meta-lowdin orthogonal AOs **\n", + " ** Mulliken pop **\n", + "pop of 0 C 1s 1.99999\n", + "pop of 0 C 2s 0.97227\n", + "pop of 0 C 3s 0.00497\n", + "pop of 0 C 2px 1.05526\n", + "pop of 0 C 2py 1.04219\n", + "pop of 0 C 2pz 0.98923\n", + "pop of 0 C 3px 0.00555\n", + "pop of 0 C 3py 0.00604\n", + "pop of 0 C 3pz 0.00423\n", + "pop of 0 C 3dxy 0.00184\n", + "pop of 0 C 3dyz 0.00103\n", + "pop of 0 C 3dz^2 0.00215\n", + "pop of 0 C 3dxz 0.00098\n", + "pop of 0 C 3dx2-y2 0.00273\n", + "pop of 1 C 1s 1.99999\n", + "pop of 1 C 2s 0.97308\n", + "pop of 1 C 3s 0.00489\n", + "pop of 1 C 2px 1.06672\n", + "pop of 1 C 2py 1.03093\n", + "pop of 1 C 2pz 0.99621\n", + "pop of 1 C 3px 0.00520\n", + "pop of 1 C 3py 0.00644\n", + "pop of 1 C 3pz 0.00476\n", + "pop of 1 C 3dxy 0.00227\n", + "pop of 1 C 3dyz 0.00100\n", + "pop of 1 C 3dz^2 0.00213\n", + "pop of 1 C 3dxz 0.00089\n", + "pop of 1 C 3dx2-y2 0.00231\n", + "pop of 2 C 1s 1.99999\n", + "pop of 2 C 2s 0.97340\n", + "pop of 2 C 3s 0.00492\n", + "pop of 2 C 2px 1.02498\n", + "pop of 2 C 2py 1.06937\n", + "pop of 2 C 2pz 1.00446\n", + "pop of 2 C 3px 0.00593\n", + "pop of 2 C 3py 0.00484\n", + "pop of 2 C 3pz 0.00495\n", + "pop of 2 C 3dxy 0.00258\n", + "pop of 2 C 3dyz 0.00081\n", + "pop of 2 C 3dz^2 0.00216\n", + "pop of 2 C 3dxz 0.00109\n", + "pop of 2 C 3dx2-y2 0.00177\n", + "pop of 3 C 1s 1.99999\n", + "pop of 3 C 2s 0.97226\n", + "pop of 3 C 3s 0.00497\n", + "pop of 3 C 2px 1.02327\n", + "pop of 3 C 2py 1.07414\n", + "pop of 3 C 2pz 0.98916\n", + "pop of 3 C 3px 0.00667\n", + "pop of 3 C 3py 0.00493\n", + "pop of 3 C 3pz 0.00423\n", + "pop of 3 C 3dxy 0.00275\n", + "pop of 3 C 3dyz 0.00082\n", + "pop of 3 C 3dz^2 0.00215\n", + "pop of 3 C 3dxz 0.00119\n", + "pop of 3 C 3dx2-y2 0.00182\n", + "pop of 4 C 1s 1.99999\n", + "pop of 4 C 2s 0.95036\n", + "pop of 4 C 3s 0.00484\n", + "pop of 4 C 2px 1.00998\n", + "pop of 4 C 2py 1.02623\n", + "pop of 4 C 2pz 0.97820\n", + "pop of 4 C 3px 0.00405\n", + "pop of 4 C 3py 0.00456\n", + "pop of 4 C 3pz 0.00317\n", + "pop of 4 C 3dxy 0.00205\n", + "pop of 4 C 3dyz 0.00120\n", + "pop of 4 C 3dz^2 0.00225\n", + "pop of 4 C 3dxz 0.00135\n", + "pop of 4 C 3dx2-y2 0.00205\n", + "pop of 5 C 1s 1.99999\n", + "pop of 5 C 2s 0.97331\n", + "pop of 5 C 3s 0.00492\n", + "pop of 5 C 2px 1.05736\n", + "pop of 5 C 2py 1.03712\n", + "pop of 5 C 2pz 1.00456\n", + "pop of 5 C 3px 0.00504\n", + "pop of 5 C 3py 0.00569\n", + "pop of 5 C 3pz 0.00494\n", + "pop of 5 C 3dxy 0.00178\n", + "pop of 5 C 3dyz 0.00099\n", + "pop of 5 C 3dz^2 0.00216\n", + "pop of 5 C 3dxz 0.00091\n", + "pop of 5 C 3dx2-y2 0.00257\n", + "pop of 6 C 1s 1.99999\n", + "pop of 6 C 2s 0.94825\n", + "pop of 6 C 3s 0.00487\n", + "pop of 6 C 2px 1.01010\n", + "pop of 6 C 2py 1.02518\n", + "pop of 6 C 2pz 0.98503\n", + "pop of 6 C 3px 0.00407\n", + "pop of 6 C 3py 0.00459\n", + "pop of 6 C 3pz 0.00318\n", + "pop of 6 C 3dxy 0.00203\n", + "pop of 6 C 3dyz 0.00119\n", + "pop of 6 C 3dz^2 0.00225\n", + "pop of 6 C 3dxz 0.00130\n", + "pop of 6 C 3dx2-y2 0.00208\n", + "pop of 7 C 1s 1.99999\n", + "pop of 7 C 2s 0.97315\n", + "pop of 7 C 3s 0.00495\n", + "pop of 7 C 2px 1.05887\n", + "pop of 7 C 2py 1.03884\n", + "pop of 7 C 2pz 0.99587\n", + "pop of 7 C 3px 0.00511\n", + "pop of 7 C 3py 0.00575\n", + "pop of 7 C 3pz 0.00479\n", + "pop of 7 C 3dxy 0.00176\n", + "pop of 7 C 3dyz 0.00100\n", + "pop of 7 C 3dz^2 0.00217\n", + "pop of 7 C 3dxz 0.00093\n", + "pop of 7 C 3dx2-y2 0.00261\n", + "pop of 8 C 1s 1.99999\n", + "pop of 8 C 2s 0.97022\n", + "pop of 8 C 3s 0.00501\n", + "pop of 8 C 2px 1.02308\n", + "pop of 8 C 2py 1.07422\n", + "pop of 8 C 2pz 0.99490\n", + "pop of 8 C 3px 0.00671\n", + "pop of 8 C 3py 0.00493\n", + "pop of 8 C 3pz 0.00424\n", + "pop of 8 C 3dxy 0.00274\n", + "pop of 8 C 3dyz 0.00080\n", + "pop of 8 C 3dz^2 0.00215\n", + "pop of 8 C 3dxz 0.00117\n", + "pop of 8 C 3dx2-y2 0.00183\n", + "pop of 9 C 1s 1.99999\n", + "pop of 9 C 2s 0.97476\n", + "pop of 9 C 3s 0.00484\n", + "pop of 9 C 2px 1.06623\n", + "pop of 9 C 2py 1.03134\n", + "pop of 9 C 2pz 0.99105\n", + "pop of 9 C 3px 0.00526\n", + "pop of 9 C 3py 0.00620\n", + "pop of 9 C 3pz 0.00480\n", + "pop of 9 C 3dxy 0.00228\n", + "pop of 9 C 3dyz 0.00103\n", + "pop of 9 C 3dz^2 0.00213\n", + "pop of 9 C 3dxz 0.00085\n", + "pop of 9 C 3dx2-y2 0.00226\n", + "pop of 10 C 1s 1.99999\n", + "pop of 10 C 2s 0.96836\n", + "pop of 10 C 3s 0.00530\n", + "pop of 10 C 2px 1.06446\n", + "pop of 10 C 2py 1.04596\n", + "pop of 10 C 2pz 0.99288\n", + "pop of 10 C 3px 0.00649\n", + "pop of 10 C 3py 0.00649\n", + "pop of 10 C 3pz 0.00462\n", + "pop of 10 C 3dxy 0.00204\n", + "pop of 10 C 3dyz 0.00098\n", + "pop of 10 C 3dz^2 0.00223\n", + "pop of 10 C 3dxz 0.00093\n", + "pop of 10 C 3dx2-y2 0.00268\n", + "pop of 11 C 1s 1.99999\n", + "pop of 11 C 2s 0.97458\n", + "pop of 11 C 3s 0.00488\n", + "pop of 11 C 2px 1.02614\n", + "pop of 11 C 2py 1.06886\n", + "pop of 11 C 2pz 0.99891\n", + "pop of 11 C 3px 0.00575\n", + "pop of 11 C 3py 0.00493\n", + "pop of 11 C 3pz 0.00499\n", + "pop of 11 C 3dxy 0.00246\n", + "pop of 11 C 3dyz 0.00075\n", + "pop of 11 C 3dz^2 0.00215\n", + "pop of 11 C 3dxz 0.00113\n", + "pop of 11 C 3dx2-y2 0.00182\n", + "pop of 12 H 1s 0.89277\n", + "pop of 12 H 2s 0.00968\n", + "pop of 12 H 2px 0.00190\n", + "pop of 12 H 2py 0.00047\n", + "pop of 12 H 2pz 0.00024\n", + "pop of 13 H 1s 0.89818\n", + "pop of 13 H 2s 0.00879\n", + "pop of 13 H 2px 0.00016\n", + "pop of 13 H 2py 0.00219\n", + "pop of 13 H 2pz 0.00025\n", + "pop of 14 H 1s 0.89339\n", + "pop of 14 H 2s 0.00966\n", + "pop of 14 H 2px 0.00018\n", + "pop of 14 H 2py 0.00219\n", + "pop of 14 H 2pz 0.00022\n", + "pop of 15 H 1s 0.89806\n", + "pop of 15 H 2s 0.00878\n", + "pop of 15 H 2px 0.00120\n", + "pop of 15 H 2py 0.00114\n", + "pop of 15 H 2pz 0.00025\n", + "pop of 16 H 1s 0.89680\n", + "pop of 16 H 2s 0.00888\n", + "pop of 16 H 2px 0.00121\n", + "pop of 16 H 2py 0.00114\n", + "pop of 16 H 2pz 0.00024\n", + "pop of 17 H 1s 0.89404\n", + "pop of 17 H 2s 0.00961\n", + "pop of 17 H 2px 0.00020\n", + "pop of 17 H 2py 0.00217\n", + "pop of 17 H 2pz 0.00022\n", + "pop of 18 H 1s 0.89499\n", + "pop of 18 H 2s 0.00973\n", + "pop of 18 H 2px 0.00183\n", + "pop of 18 H 2py 0.00053\n", + "pop of 18 H 2pz 0.00024\n", + "pop of 19 H 1s 0.89993\n", + "pop of 19 H 2s 0.00882\n", + "pop of 19 H 2px 0.00016\n", + "pop of 19 H 2py 0.00216\n", + "pop of 19 H 2pz 0.00025\n", + "pop of 20 H 1s 0.89349\n", + "pop of 20 H 2s 0.00966\n", + "pop of 20 H 2px 0.00145\n", + "pop of 20 H 2py 0.00092\n", + "pop of 20 H 2pz 0.00022\n", + "pop of 21 H 1s 0.87426\n", + "pop of 21 H 2s 0.01239\n", + "pop of 21 H 2px 0.00150\n", + "pop of 21 H 2py 0.00093\n", + "pop of 21 H 2pz 0.00022\n", + " ** Mulliken atomic charges **\n", + "charge of 0C = -0.08845\n", + "charge of 1C = -0.09682\n", + "charge of 2C = -0.10124\n", + "charge of 3C = -0.08835\n", + "charge of 4C = 0.00973\n", + "charge of 5C = -0.10131\n", + "charge of 6C = 0.00590\n", + "charge of 7C = -0.09579\n", + "charge of 8C = -0.09198\n", + "charge of 9C = -0.09302\n", + "charge of 10C = -0.10340\n", + "charge of 11C = -0.09736\n", + "charge of 12H = 0.09494\n", + "charge of 13H = 0.09044\n", + "charge of 14H = 0.09436\n", + "charge of 15H = 0.09057\n", + "charge of 16H = 0.09173\n", + "charge of 17H = 0.09375\n", + "charge of 18H = 0.09269\n", + "charge of 19H = 0.08868\n", + "charge of 20H = 0.09426\n", + "charge of 21H = 0.11070\n", + "Dipole moment(X, Y, Z, Debye): 0.04456, 0.00167, -0.00000\n", + "size of MO coeff (218, 218)\n", + "converged SCF energy = -460.279584264496\n", + "41\n", + "0\n", + "177\n" + ] + } + ], + "source": [ + "molecule = '''\n", + "C -6.8604981940 -0.4376683883 0.0000000000\n", + "C -7.0417938838 0.9473143671 0.0000000000\n", + "C -5.5765818320 -0.9841959654 0.0000000000\n", + "C -5.9160648559 1.7742824642 0.0000000000\n", + "C -4.4255614400 -0.1700917199 0.0000000000\n", + "C -4.6336867757 1.2242295397 0.0000000000\n", + "C -3.0514687259 -0.7633813258 0.0000000000\n", + "C -2.8520613255 -2.1570405955 0.0000000000\n", + "C -1.5640941394 -2.6887327507 0.0000000000\n", + "C -0.4451033014 -1.8584342989 0.0000000000\n", + "C -0.5927492072 -0.4581981628 0.0000000000\n", + "C -1.9041039939 0.0505338987 0.0000000000\n", + "H -8.0463585128 1.3756739224 0.0000000000\n", + "H -5.4821758879 -2.0696091478 0.0000000000\n", + "H -6.0325833442 2.8604965477 0.0000000000\n", + "H -3.7849879120 1.9076345469 0.0000000000\n", + "H -3.6983932318 -2.8425398416 0.0000000000\n", + "H -1.4298106521 -3.7726953568 0.0000000000\n", + "H 0.5426199813 -2.3168156796 0.0000000000\n", + "H -2.0368540411 1.1280906363 0.0000000000\n", + "H -7.7260669052 -1.1042415323 0.0000000000\n", + "H 0.22700 0.16873 0.00000\n", + "'''\n", + "basis = \"cc-pvdz\"\n", + "pymol = pyscf.gto.Mole(\n", + " atom = molecule,\n", + " #symmetry= True,\n", + " spin = 0, # number of unpaired electrons\n", + " charge = 0,\n", + " basis = basis)\n", + "\n", + "\n", + "pymol.build()\n", + "#print(\"symmetry: \",pymol.topgroup)\n", + "mf = pyscf.scf.RHF(pymol)\n", + "mf.verbose = 4\n", + "mf.conv_tol = 1e-8\n", + "#mf.conv_tol_grad = 1e-5\n", + "#mf.init_guess = \"sad\"\n", + "\n", + "mf.run(max_cycle=200)\n", + "\n", + "print(\" Hartree-Fock Energy: %12.8f\" % mf.e_tot)\n", + "mf.analyze()\n", + "# Get data\n", + "F = mf.get_fock()\n", + "C = mf.mo_coeff\n", + "S = mf.get_ovlp()\n", + "print(\"size of MO coeff\", C.shape)\n", + "ndocc = pymol.nelec[1]\n", + "nsing = pymol.nelec[0] - ndocc\n", + "nvirt = pymol.mol.nao - ndocc - nsing\n", + "print(ndocc)\n", + "print(nsing)\n", + "print(nvirt)\n", + "# Just use alpha orbitals\n", + "Cdocc = mf.mo_coeff[:,0:ndocc]\n", + "Csing = mf.mo_coeff[:,ndocc:ndocc+nsing]\n", + "Cvirt = mf.mo_coeff[:,ndocc+nsing:ndocc+nsing+nvirt]\n", + "np.save(\"C_scf\", C)\n", + "\n", + "nbas = Cdocc.shape[0]" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "(0, 'C', '2p', 'z')\n", + "(1, 'C', '2p', 'z')\n", + "(2, 'C', '2p', 'z')\n", + "(3, 'C', '2p', 'z')\n", + "(4, 'C', '2p', 'z')\n", + "(5, 'C', '2p', 'z')\n", + "[[5, 19, 33, 47, 61, 75], [89, 103, 117, 131, 145, 159]]\n", + "6\n", + "6\n" + ] + } + ], + "source": [ + "# Find AO's corresponding to atoms\n", + "full = []\n", + "frag1 = []\n", + "frag2 = []\n", + "for ao_idx,ao in enumerate(mf.mol.ao_labels(fmt=False)):\n", + " #print(ao_idx)\n", + " if ao[0] in (0,1,2,3,4,5):\n", + " if ao[2] in (\"2p\"):\n", + " if ao[3] == \"z\":\n", + " print(ao)\n", + " frag1.append(ao_idx)\n", + " full.append(ao_idx)\n", + " if ao[0] in (6,7,8,9,10,11):\n", + " if ao[2] in (\"2p\"):\n", + " if ao[3] == \"z\":\n", + " frag2.append(ao_idx)\n", + " full.append(ao_idx) \n", + "\n", + "\n", + "frags = [frag1, frag2]\n", + "print(frags)\n", + "print(len(frag1))\n", + "print(len(frag2))" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "(218, 12)" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "# Define projectors\n", + "nbas = Cdocc.shape[0]\n", + "X = scipy.linalg.sqrtm(S)\n", + "X = np.eye(nbas) \n", + "Pfull = X[:,full] # non-orthogonal\n", + "Pf = []\n", + "for f in frags:\n", + " Pf.append(X[:,f])\n", + "display(Pfull.shape)" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + " Partition 218 orbitals into a total of 12 orbitals\n", + " Index Sing. Val. Space \n", + " 0 1.00000000 2*\n", + " 1 1.00000000 2*\n", + " 2 1.00000000 2*\n", + " 3 1.00000000 2*\n", + " 4 1.00000000 2*\n", + " 5 1.00000000 2*\n", + " 6 0.91524412 0*\n", + " 7 0.91023129 0*\n", + " 8 0.90788521 0*\n", + " 9 0.90607946 0*\n", + " 10 0.90555587 0*\n", + " 11 0.90319379 0*\n", + " 12 0.42923301 2\n", + " 13 0.42422702 2\n", + " 14 0.42310756 2\n", + " 15 0.41921885 2\n", + " 16 0.41410022 2\n", + " 17 0.40289975 2\n", + " SVD active space has the following dimensions:\n", + " Orbital Block Environment Active\n", + " 0 35 6\n", + " 1 0 0\n", + " 2 171 6\n" + ] + }, + { + "data": { + "text/plain": [ + "array([[-1.71179287e-16, 4.27395488e-17, 1.62212117e-16, ...,\n", + " -1.46293922e-16, -2.77530984e-17, -7.65855847e-17],\n", + " [-5.73040115e-17, 2.09893551e-16, 1.72926708e-16, ...,\n", + " -3.43271443e-16, -5.32738586e-17, -5.75274797e-17],\n", + " [ 2.63912186e-15, -1.14872390e-15, -3.02457415e-15, ...,\n", + " 9.66902418e-16, 1.14079139e-15, 2.12681624e-15],\n", + " ...,\n", + " [-8.13411488e-17, -1.69089094e-17, 1.71573730e-16, ...,\n", + " -7.08096169e-18, -7.79065050e-17, -4.63592939e-17],\n", + " [ 2.15328339e-16, 1.36589067e-16, 2.10054685e-16, ...,\n", + " -6.16834514e-17, 7.69831224e-17, -6.19860666e-17],\n", + " [ 1.54958226e-03, -6.51683915e-03, 1.15484626e-03, ...,\n", + " -6.53734282e-17, 5.26755623e-17, 4.72958075e-17]])" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + " Should be 1: 0.9999999999999996\n" + ] + } + ], + "source": [ + "(Oact, Sact, Vact), (Cenv, Cerr, _) = svd_subspace_partitioning((Cdocc, Csing, Cvirt), Pfull, S)\n", + "assert(Cerr.shape[1] == 0)\n", + "Cact = np.hstack((Oact, Sact, Vact))\n", + "display(Cact)\n", + "pyscf.tools.molden.from_mo(mf.mol, \"Cact.molden\", Cact)\n", + "print(\" Should be 1: \", np.linalg.det(Cact.T @ S @ Cact))" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n", + " Fragment: [5, 19, 33, 47, 61, 75]\n", + " Partition 12 orbitals into a total of 6 orbitals\n", + " Index Sing. Val. Space \n", + " 0 0.99999643 2*\n", + " 1 0.99994058 2*\n", + " 2 0.98685246 2*\n", + " 3 0.90700805 0*\n", + " 4 0.90607117 0*\n", + " 5 0.90440602 0*\n", + " 6 0.14736696 0\n", + " 7 0.12288367 2\n", + " 8 0.01096706 2\n", + " 9 0.00988094 0\n", + " 10 0.00241800 0\n", + " 11 0.00152824 2\n", + " SVD active space has the following dimensions:\n", + " Orbital Block Environment Active\n", + " 0 3 3\n", + " 1 0 0\n", + " 2 3 3\n", + "\n", + " Fragment: [89, 103, 117, 131, 145, 159]\n", + " Partition 12 orbitals into a total of 6 orbitals\n", + " Index Sing. Val. Space \n", + " 0 0.99999658 2*\n", + " 1 0.99993926 2*\n", + " 2 0.98727403 2*\n", + " 3 0.90689569 0*\n", + " 4 0.90510125 0*\n", + " 5 0.90429286 0*\n", + " 6 0.14495947 0\n", + " 7 0.12518690 2\n", + " 8 0.01084660 2\n", + " 9 0.00999707 0\n", + " 10 0.00237129 0\n", + " 11 0.00156636 2\n", + " SVD active space has the following dimensions:\n", + " Orbital Block Environment Active\n", + " 0 3 3\n", + " 1 0 0\n", + " 2 3 3\n", + "[-0.50056208 -0.33842113 -0.33590706 0.34236119 0.34670818 0.53073629]\n", + "[-0.4998164 -0.33821807 -0.33560099 0.34109111 0.34652785 0.52908879]\n", + "nick: [1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]\n", + " init_fspace = [(3, 3), (3, 3)]\n", + " clusters = [[1, 2, 3, 4, 5, 6], [7, 8, 9, 10, 11, 12]]\n" + ] + } + ], + "source": [ + "# Project active orbitals onto fragments\n", + "init_fspace = []\n", + "clusters = []\n", + "Cfrags = []\n", + "orb_index = 1\n", + "\n", + "\n", + "for fi,f in enumerate(frags):\n", + " print()\n", + " print(\" Fragment: \", f)\n", + " (Of, Sf, Vf), (_, _, _) = svd_subspace_partitioning((Oact, Sact, Vact), Pf[fi], S)\n", + " # (Of, Sf, Vf), (Oact, Sact, Vact) = svd_subspace_partitioning((Oact, Sact, Vact), Pf[fi], S)\n", + "\n", + " Cfrags.append(np.hstack((Of, Sf, Vf)))\n", + " ndocc_f = Of.shape[1]\n", + " init_fspace.append((ndocc_f+Sf.shape[1], ndocc_f))\n", + " nmof = Of.shape[1] + Sf.shape[1] + Vf.shape[1]\n", + " clusters.append(list(range(orb_index, orb_index+nmof)))\n", + " orb_index += nmof\n", + "\n", + "\n", + "# Orthogonalize Fragment orbitals\n", + "Cfrags = sym_ortho(Cfrags, S)\n", + "\n", + "# Pseudo canonicalize fragments\n", + "Cfrags = canonicalize(Cfrags, F)\n", + "\n", + "Cact = np.hstack(Cfrags)\n", + "\n", + "print(\"nick: \", np.linalg.svd(Cact.T @ S @ Cact)[1])\n", + "# Write Molden files for visualization\n", + "pyscf.tools.molden.from_mo(mf.mol, \"Pfull.molden\", Pfull)\n", + "pyscf.tools.molden.from_mo(mf.mol, \"Cact.molden\", Cact)\n", + "pyscf.tools.molden.from_mo(mf.mol, \"Cenv.molden\", Cenv)\n", + "for i in range(len(frags)):\n", + " pyscf.tools.molden.from_mo(mf.mol, \"Cfrag%i.molden\"%i, Cfrags[i])\n", + "\n", + "\n", + "print(\" init_fspace = \", init_fspace)\n", + "print(\" clusters = \", clusters)" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "(218, 35)\n", + "(218, 12)\n", + "(218, 218)\n" + ] + } + ], + "source": [ + "print(Cenv.shape)\n", + "print(Cact.shape)\n", + "d1_embed = 2 * Cenv @ Cenv.T\n", + "\n", + "h0 = pyscf.gto.mole.energy_nuc(mf.mol)\n", + "h = pyscf.scf.hf.get_hcore(mf.mol)\n", + "j, k = pyscf.scf.hf.get_jk(mf.mol, d1_embed, hermi=1)\n", + "\n", + "print(h.shape)\n", + "h0 += np.trace(d1_embed @ ( h + .5*j - .25*k))\n", + "\n", + "h = Cact.T @ h @ Cact;\n", + "j = Cact.T @ j @ Cact;\n", + "k = Cact.T @ k @ Cact;\n", + "nact = h.shape[0]\n", + "\n", + "h2 = pyscf.ao2mo.kernel(pymol, Cact, aosym=\"s4\", compact=False)\n", + "h2.shape = (nact, nact, nact, nact)\n", + "# The use of d1_embed only really makes sense if it has zero electrons in the\n", + "# active space. Let's warn the user if that's not true\n", + "\n", + "S = pymol.intor(\"int1e_ovlp_sph\")\n", + "n_act = np.trace(S @ d1_embed @ S @ Cact @ Cact.T)\n", + "if abs(n_act) > 1e-8 == False:\n", + " print(n_act)\n", + " error(\" I found embedded electrons in the active space?!\")\n", + "\n", + "h1 = h + j - .5*k;\n", + "\n", + "np.save(\"ints_h0\", h0)\n", + "np.save(\"ints_h1\", h1)\n", + "np.save(\"ints_h2\", h2)\n", + "np.save(\"mo_coeffs\", Cact)\n", + "np.save(\"overlap_mat\", S)\n", + "\n", + "Pa = mf.make_rdm1()[0]\n", + "Pb = mf.make_rdm1()[1]\n", + "#np.save(\"Pa\", Cact.T @ S @ Pa @ S @ Cact)\n", + "#np.save(\"Pb\", Cact.T @ S @ Pb @ S @ Cact)" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "mim", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.7.6" + }, + "orig_nbformat": 4 + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/examples/two_benzene/ints_h0.npy b/examples/two_benzene/ints_h0.npy new file mode 100644 index 0000000..8926a65 Binary files /dev/null and b/examples/two_benzene/ints_h0.npy differ diff --git a/examples/two_benzene/ints_h1.npy b/examples/two_benzene/ints_h1.npy new file mode 100644 index 0000000..4e6b36d Binary files /dev/null and b/examples/two_benzene/ints_h1.npy differ diff --git a/examples/two_benzene/ints_h2.npy b/examples/two_benzene/ints_h2.npy new file mode 100644 index 0000000..2cbe9f1 Binary files /dev/null and b/examples/two_benzene/ints_h2.npy differ diff --git a/examples/two_benzene/mo_coeffs.npy b/examples/two_benzene/mo_coeffs.npy new file mode 100644 index 0000000..97e08ee Binary files /dev/null and b/examples/two_benzene/mo_coeffs.npy differ diff --git a/src/incore_cmf.jl b/src/incore_cmf.jl index 645fb58..0f4f4c9 100644 --- a/src/incore_cmf.jl +++ b/src/incore_cmf.jl @@ -123,9 +123,6 @@ function cmf_ci_iteration(ints::InCoreInts{T}, clusters::Vector{MOCluster}, in_r verbose < 2 || @printf(" Slater Det Energy: %12.8f\n", e) else - # - # run PYSCF FCI - #e, d1a,d1b, d2 = pyscf_fci(ints_i,fspace[ci.idx][1],fspace[ci.idx][2], verbose=verbose) if use_pyscf pyscf = pyimport("pyscf") @@ -201,7 +198,7 @@ function cmf_ci_iteration(ints::InCoreInts{T}, clusters::Vector{MOCluster}, in_r end """ - cmf_ci_iteration(ints::InCoreInts{T}, clusters::Vector{MOCluster}, in_rdm1::RDM1{T}, fspace, ansatze::Vector{Vector{A}}; + cmf_ci_iteration(ints::InCoreInts{T}, clusters::Vector{MOCluster}, in_rdm1::RDM1{T}, fspace, ansatze::Vector{<:Ansatz}; use_pyscf = true, verbose = 1, sequential= false, @@ -212,8 +209,8 @@ end Perform single CMF-CI iteration, returning new energy, and density """ -function cmf_ci_iteration(ints::InCoreInts{T}, clusters::Vector{MOCluster}, in_rdm1::RDM1{T}, fspace, ansatze::Vector{Vector{Ansatz}}; - use_pyscf = false, +function cmf_ci_iteration(ints::InCoreInts{T}, clusters::Vector{MOCluster}, in_rdm1::RDM1{T}, fspace, ansatze::Vector{<:Ansatz}; + use_pyscf = true, verbose = 1, sequential= false, spin_avg = true, @@ -228,7 +225,7 @@ function cmf_ci_iteration(ints::InCoreInts{T}, clusters::Vector{MOCluster}, in_r ci = clusters[i] flush(stdout) - ansatz = ansatze[i][1] + ansatz = ansatze[i] verbose < 2 || display(ansatz) ints_i = subset(ints, ci, rdm1) @@ -255,27 +252,38 @@ function cmf_ci_iteration(ints::InCoreInts{T}, clusters::Vector{MOCluster}, in_r db = Matrix(1.0I, no, no) end - if spin_avg - da = .5*(da + db) - db .= da - end - d1 = RDM1(da,db) d2 = RDM2(d1) e = compute_energy(ints_i, d1) verbose < 2 || @printf(" Slater Det Energy: %12.8f\n", e) else - # + + if use_pyscf && typeof(ansatz) == FCIAnsatz + pyscf = pyimport("pyscf") + fci = pyimport("pyscf.fci") + cisolver = pyscf.fci.direct_spin1.FCI() + cisolver.max_cycle = maxiter_ci + cisolver.conv_tol = tol_ci + cisolver.conv_tol_residual = tol_ci + #cisolver.lindep = 1e-8 + nelec = na + nb + e, vfci = cisolver.kernel(ints_i.h1, ints_i.h2, no, (na,nb), ecore=ints_i.h0) + (d1a, d1b), (d2aa, d2ab, d2bb) = cisolver.make_rdm12s(vfci, no, (na,nb)) + + d1 = RDM1(d1a, d1b) + d2 = RDM2(d2aa, d2ab, d2bb) + else - solver = SolverSettings(verbose=1, tol=tol_ci, maxiter=maxiter_ci) - solution = solve(ints_i, ansatz, solver) - d1a, d1b, d2aa, d2bb, d2ab = compute_1rdm_2rdm(solution) + solver = SolverSettings(verbose=1, tol=tol_ci, maxiter=maxiter_ci) + solution = solve(ints_i, ansatz, solver) + d1a, d1b, d2aa, d2bb, d2ab = compute_1rdm_2rdm(solution) - verbose < 2 || display(solution) + verbose < 2 || display(solution) - d1 = RDM1(d1a, d1b) - d2 = RDM2(d2aa, d2ab, d2bb) + d1 = RDM1(d1a, d1b) + d2 = RDM2(d2aa, d2ab, d2bb) + end end if spin_avg @@ -342,6 +350,7 @@ function cmf_ci(ints, clusters, fspace, in_rdm1::RDM1; tol_d1 = 1e-6, tol_ci = 1e-8, verbose = 1, + use_pyscf = true, sequential = false) rdm1 = deepcopy(in_rdm1) energies = [] @@ -363,6 +372,7 @@ function cmf_ci(ints, clusters, fspace, in_rdm1::RDM1; maxiter_ci = maxiter_ci, tol_ci = tol_ci, verbose = verbose, + use_pyscf = use_pyscf, sequential = sequential ) rdm1_curr = assemble_full_rdm(clusters, rdm1_dict) @@ -393,13 +403,13 @@ function cmf_ci(ints, clusters, fspace, in_rdm1::RDM1; end """ - cmf_ci(ints, clusters, fspace, ansatze::Vector{Vector{A}}, in_rdm1::RDM1; + cmf_ci(ints, clusters, fspace, ansatze::Vector{<:Ansatz}, in_rdm1::RDM1; maxiter_ci = 100, maxiter_d1 = 20, tol_d1 = 1e-6, tol_ci = 1e-8, verbose = 1, - sequential = false) where A + sequential = false) Optimize the 1RDM for CMF-CI @@ -420,12 +430,13 @@ Optimize the 1RDM for CMF-CI - `rdm1_dict`: Dictionary of 1RDMs; cluster index --> RDM1 - `rdm2_dict`: Dictionary of 2RDMs; cluster index --> RDM2 """ -function cmf_ci(ints, clusters, fspace, ansatze::Vector{Vector{Ansatz}}, in_rdm1::RDM1; +function cmf_ci(ints, clusters, fspace, ansatze::Vector{<:Ansatz}, in_rdm1::RDM1; maxiter_ci = 100, maxiter_d1 = 20, tol_d1 = 1e-6, tol_ci = 1e-8, verbose = 1, + use_pyscf = true, sequential = false) rdm1 = deepcopy(in_rdm1) energies = [] @@ -447,6 +458,7 @@ function cmf_ci(ints, clusters, fspace, ansatze::Vector{Vector{Ansatz}}, in_rdm1 maxiter_ci = maxiter_ci, tol_ci = tol_ci, verbose = verbose, + use_pyscf = use_pyscf, sequential = sequential ) rdm1_curr = assemble_full_rdm(clusters, rdm1_dict) @@ -508,6 +520,7 @@ function cmf_oo(ints::InCoreInts{T}, clusters::Vector{MOCluster}, fspace, dguess verbose=0, method="bfgs", alpha=nothing, + use_pyscf=true, sequential=false) where T norb = size(ints.h1)[1] @@ -539,6 +552,7 @@ function cmf_oo(ints::InCoreInts{T}, clusters::Vector{MOCluster}, fspace, dguess e, rdm1_dict, _ = cmf_ci(ints_tmp, clusters, fspace, orbital_rotation(d1_curr, U), tol_d1=gconv/10.0, verbose=0, + use_pyscf= use_pyscf, sequential=sequential) d1_tmp = assemble_full_rdm(clusters, rdm1_dict) @@ -619,7 +633,7 @@ function cmf_oo(ints::InCoreInts{T}, clusters::Vector{MOCluster}, fspace, dguess end """ - cmf_oo(ints::InCoreInts{T}, clusters::Vector{MOCluster}, fspace, ansatze::Vector{Vector{Ansatz}}, dguess::RDM1{T}; + cmf_oo(ints::InCoreInts{T}, clusters::Vector{MOCluster}, fspace, ansatze::Vector{<:Ansatz}, dguess::RDM1{T}; max_iter_oo=100, max_iter_ci=100, gconv=1e-6, @@ -642,13 +656,14 @@ Do CMF with orbital optimization - `verbose`: Printing level - `method`: optimization method """ -function cmf_oo(ints::InCoreInts{T}, clusters::Vector{MOCluster}, fspace, ansatze::Vector{Vector{Ansatz}}, dguess::RDM1{T}; +function cmf_oo(ints::InCoreInts{T}, clusters::Vector{MOCluster}, fspace, ansatze::Vector{<:Ansatz}, dguess::RDM1{T}; max_iter_oo=100, max_iter_ci=100, gconv=1e-6, verbose=0, method="bfgs", alpha=nothing, + use_pyscf=true, sequential=false) where T norb = size(ints.h1)[1] @@ -680,6 +695,7 @@ function cmf_oo(ints::InCoreInts{T}, clusters::Vector{MOCluster}, fspace, ansatz e, rdm1_dict, _ = cmf_ci(ints_tmp, clusters, fspace, ansatze, orbital_rotation(d1_curr, U), tol_d1=gconv/10.0, verbose=0, + use_pyscf= use_pyscf, sequential=sequential) d1_tmp = assemble_full_rdm(clusters, rdm1_dict) @@ -885,13 +901,13 @@ function orbital_objective_function(ints, clusters, kappa, fspace, rdm::RDM1; end """ - orbital_objective_function(ints, clusters, kappa, fspace, ansatze::Vector{Vector{Ansatz}}, da, db; + orbital_objective_function(ints, clusters, kappa, fspace, ansatze::Vector{<:Ansatz}, da, db; ci_conv = 1e-9, sequential = false, verbose = 1) Objective function to minimize in OO-CMF """ -function orbital_objective_function(ints, clusters, kappa, fspace, ansatze::Vector{Vector{Ansatz}}, rdm::RDM1{T}; +function orbital_objective_function(ints, clusters, kappa, fspace, ansatze::Vector{<:Ansatz}, rdm::RDM1{T}; ci_conv = 1e-9, sequential = false, verbose = 0) where T @@ -935,13 +951,13 @@ function orbital_gradient_numerical(ints, clusters, kappa, fspace, d::RDM1; end """ - orbital_gradient_numerical(ints, clusters, kappa, fspace, ansatze::Vector{Vector{Ansatz}}, da, db; + orbital_gradient_numerical(ints, clusters, kappa, fspace, ansatze::Vector{<:Ansatz}, da, db; gconv = 1e-8, verbose = 1, stepsize = 1e-6) Compute orbital gradient with finite difference """ -function orbital_gradient_numerical(ints, clusters, kappa, fspace, ansatze::Vector{Vector{Ansatz}}, d::RDM1; +function orbital_gradient_numerical(ints, clusters, kappa, fspace, ansatze::Vector{<:Ansatz}, d::RDM1; ci_conv = 1e-10, verbose = 0, stepsize = 1e-6) @@ -1075,6 +1091,118 @@ function cmf_oo_gd( ints_in::InCoreInts{T}, clusters::Vector{MOCluster}, fspace, #=}}}=# end +""" + cmf_oo_gd( ints_in::InCoreInts{T}, clusters::Vector{MOCluster}, fspace, ansatze::Vector{<:Ansatz}}, dguess::RDM1{T}; + maxiter_oo = 100, + maxiter_ci = 100, + maxiter_d1 = 100, + tol_oo = 1e-6, + tol_d1 = 1e-7, + tol_ci = 1e-8, + verbose = 0, + alpha = .1, + zero_intra_rots = true, + sequential = false) where T + +Do CMF with orbital optimization + +# Arguments + +- `ints::InCoreInts`: integrals for full system +- `clusters::Vector{MOCluster}`: vector of cluster objects +- `fspace::Vector{Vector{Integer}}`: vector of particle number occupations for each cluster specifying the sectors of fock space +- `dguess`: initial guess for 1particle density matrix +- `maxiter_oo`: Max iter for the orbital optimization iterations +- `maxiter_d1`: Max iter for the cmf iteration for the 1RDM +- `maxiter_ci`: Max iter for the CI diagonalization of the cluster states +- `tol_oo`: Convergence threshold for change in orbital gradient +- `tol_ci`: Convergence threshold for the cluster CI problems +- `tol_d1`: Convergence threshold for the CMF 1RDM +- `sequential`: If true use the density matrix of the previous cluster in a cMF iteration to form effective integrals. Improves comvergence, may depend on cluster orderings +- `verbose`: Printing level + +# Returns + +- `e`: Energy +- `U::Matrix`: Orbital rotation matrix from input to output orbitals +- `d1::RDM1`: Optimized 1RDM in the optimized orbital basis +""" +function cmf_oo_gd( ints_in::InCoreInts{T}, clusters::Vector{MOCluster}, fspace, ansatze::Vector{<:Ansatz}, dguess::RDM1{T}; + maxiter_oo = 100, + maxiter_ci = 100, + maxiter_d1 = 100, + tol_oo = 1e-6, + tol_d1 = 1e-7, + tol_ci = 1e-8, + verbose = 0, + alpha = .1, + zero_intra_rots = true, + sequential = false + ) where T + #={{{=# + ints = deepcopy(ints_in) + norb = n_orb(ints) + d1 = deepcopy(dguess) + U = Matrix(1.0I, norb, norb) + e = 0.0 + + function step!(ints, d1, k) + norb = n_orb(ints) + K = unpack_gradient(k, norb) + if zero_intra_rots + # Remove intracluster rotations + for ci in clusters + K[ci.orb_list, ci.orb_list] .= 0 + end + end + + Ui = exp(K) + + tmp = orbital_rotation(ints,Ui) + ints.h1 .= tmp.h1 + ints.h2 .= tmp.h2 + + tmp = orbital_rotation(d1,Ui) + d1.a .= tmp.a + d1.b .= tmp.b + + + e_i, rdm1_dict, rdm2_dict = cmf_ci(ints, clusters, fspace, ansatze, d1, + maxiter_d1 = maxiter_d1, + maxiter_ci = maxiter_ci, + tol_d1 = tol_d1, + tol_ci = tol_ci, + verbose = 0, + sequential = sequential) + + gd1, gd2 = assemble_full_rdm(clusters, rdm1_dict, rdm2_dict) + g = build_orbital_gradient(ints, gd1, gd2) + return e_i, g, Ui, gd1 + end + + # Compute initial gradient + converged = false + step_i = zeros(norb*(norb-1)÷2) + for i in 1:maxiter_oo + ei, gi, Ui, d1 = step!(ints, d1, step_i) + step_i = -alpha*gi + e = ei + U = U*Ui + + converged = norm(gi) < tol_oo + if converged + @printf("*Step: %4i E: %16.15f G: %4.3e\n", i, ei, norm(gi)) + break + else + @printf(" Step: %4i E: %16.15f G: %4.3e\n", i, ei, norm(gi)) + end + end + + return e, U, d1 +#=}}}=# +end + + """ cmf_oo_diis( ints_in::InCoreInts{T}, clusters::Vector{MOCluster}, fspace, dguess::RDM1{T}; maxiter_oo = 100, @@ -1087,7 +1215,7 @@ end max_ss_size = 8, diis_start = 1, alpha = .1, - zero_intra_rots = true + zero_intra_rots = true, sequential = false ) where T @@ -1166,6 +1294,10 @@ function cmf_oo_diis(ints_in::InCoreInts{T}, clusters::Vector{MOCluster}, fspace d1_i = orbital_rotation(d1_i, Ui') d2_i = orbital_rotation(d2_i, Ui') g_i = build_orbital_gradient(ints, d1_i, d2_i) + + if verbose == 1 + display(unpack_gradient(g_i, norb)) + end #g_i = build_orbital_gradient(ints_i, d1_i, d2_i) if zero_intra_rots @@ -1189,11 +1321,11 @@ function cmf_oo_diis(ints_in::InCoreInts{T}, clusters::Vector{MOCluster}, fspace #g_ss = hcat(g_ss, g_i) #k_ss = hcat(k_ss, k_i) nss = size(g_ss,2) - - for i in 1:maxiter_oo - k_i = k_i - alpha*g_i + for i in 1:maxiter_oo + k_i = k_i - alpha*g_i + if nss < max_ss_size nss += 1 g_ss = hcat(g_ss, g_i) @@ -1275,7 +1407,7 @@ function cmf_oo_diis(ints_in::InCoreInts{T}, clusters::Vector{MOCluster}, fspace # # Compute energy and gradient # - e_i, g_i, d1_i = step!(k_i) + e_i, g_i, d1_i = step!(k_i) g_i = reshape(g_i, (norb2,1)) k_i = reshape(k_i, (norb2,1)) d_i = d1_i @@ -1303,7 +1435,7 @@ function cmf_oo_diis(ints_in::InCoreInts{T}, clusters::Vector{MOCluster}, fspace end """ - cmf_oo_diis( ints_in::InCoreInts{T}, clusters::Vector{MOCluster}, fspace, ansatze::Vector{Vector{A}}, dguess::RDM1{T}; + cmf_oo_diis( ints_in::InCoreInts{T}, clusters::Vector{MOCluster}, fspace, ansatze::Vector{<:Ansatz}, dguess::RDM1{T}; maxiter_oo = 100, maxiter_ci = 100, maxiter_d1 = 100, @@ -1314,9 +1446,10 @@ end max_ss_size = 8, diis_start = 1, alpha = .1, - zero_intra_rots = true + zero_intra_rots = true, + orb_hessian = true, sequential = false - ) where A, T + ) where T Do CMF with orbital optimization using DIIS @@ -1345,7 +1478,7 @@ Do CMF with orbital optimization using DIIS - `U::Matrix`: Orbital rotation matrix from input to output orbitals - `d1::RDM1`: Optimized 1RDM in the optimized orbital basis """ -function cmf_oo_diis(ints_in::InCoreInts{T}, clusters::Vector{MOCluster}, fspace, ansatze::Vector{Vector{Ansatz}}, dguess::RDM1{T}; +function cmf_oo_diis(ints_in::InCoreInts{T}, clusters::Vector{MOCluster}, fspace, ansatze::Vector{<:Ansatz}, dguess::RDM1{T}; maxiter_oo = 100, maxiter_ci = 100, maxiter_d1 = 100, @@ -1356,8 +1489,8 @@ function cmf_oo_diis(ints_in::InCoreInts{T}, clusters::Vector{MOCluster}, fspace max_ss_size = 8, diis_start = 1, alpha = .1, - #zero_intra_rots = true, - zero_intra_rots = false, + zero_intra_rots = true, + orb_hessian = true, sequential = false ) where T #={{{=# @@ -1374,6 +1507,7 @@ function cmf_oo_diis(ints_in::InCoreInts{T}, clusters::Vector{MOCluster}, fspace g_ss = zeros(T,norb2,0) converged = false condB = 0.0 + function step!(k) K = unpack_gradient(k, norb) @@ -1397,24 +1531,35 @@ function cmf_oo_diis(ints_in::InCoreInts{T}, clusters::Vector{MOCluster}, fspace g_i = build_orbital_gradient(ints, d1_i, d2_i) #g_i = orbital_gradient_numerical(ints, clusters, k, fspace, ansatze, d1) #g_i = build_orbital_gradient(ints_i, d1_i, d2_i) + if orb_hessian + h = RDM.build_orbital_hessian(ints,d1_i,d2_i) + else + h = nothing + end + if verbose == 1 display(unpack_gradient(g_i, norb)) end - if zero_intra_rots - g_i = unpack_gradient(g_i, norb) - for ci in clusters - g_i[ci.orb_list, ci.orb_list] .= 0 - end - g_i = pack_gradient(g_i, norb) - end + #if zero_intra_rots + # g_i = unpack_gradient(g_i, norb) + # for ci in clusters + # g_i[ci.orb_list, ci.orb_list] .= 0 + # end + # g_i = pack_gradient(g_i, norb) + # + # proj_vec = projection_vector(ansatze, norb) + # h = proj_vec'*h*proj_vec + # g_i = proj_vec'*g_i + #end + e = e_i U = Ui - return e_i, g_i, d1_i + return e_i, g_i, d1_i, h end # First step - e_i, g_i, d1_i = step!(zeros(norb2)) + e_i, g_i, d1_i, h_i = step!(zeros(norb2)) k_i = zeros(norb2) #g_i = reshape(g_i, (norb2,1)) @@ -1423,19 +1568,36 @@ function cmf_oo_diis(ints_in::InCoreInts{T}, clusters::Vector{MOCluster}, fspace #k_ss = hcat(k_ss, k_i) nss = size(g_ss,2) + if zero_intra_rots + proj_vec = projection_vector(ansatze, norb) + end + for i in 1:maxiter_oo - - k_i = k_i - alpha*g_i + #project out invarant orbital rotations + if zero_intra_rots && orb_hessian + tmp_step = (pinv(proj_vec'*h_i*proj_vec))*(proj_vec'*g_i) + step_i = proj_vec*tmp_step + elseif orb_hessian && zero_intra_rots==false + step_i = pinv(h_i)*g_i + elseif orb_hessian==false && zero_intra_rots + tmp_step = proj_vec'*g_i + g_i = proj_vec*tmp_step + step_i = alpha*g_i + else + step_i = alpha*g_i + end + + k_i = k_i - step_i if nss < max_ss_size nss += 1 - g_ss = hcat(g_ss, g_i) + g_ss = hcat(g_ss, step_i) k_ss = hcat(k_ss, k_i) else g_ss[:,1:end-1] .= g_ss[:,2:end] k_ss[:,1:end-1] .= k_ss[:,2:end] - g_ss[:,nss] .= g_i + g_ss[:,nss] .= step_i k_ss[:,nss] .= k_i end @@ -1496,19 +1658,20 @@ function cmf_oo_diis(ints_in::InCoreInts{T}, clusters::Vector{MOCluster}, fspace verbose < 2 || display(g_ss) end - if zero_intra_rots - # Remove intracluster rotations - k_i = unpack_gradient(k_i, norb) - for ci in clusters - k_i[ci.orb_list, ci.orb_list] .= 0 - end - k_i = pack_gradient(k_i, norb) - end + #if zero_intra_rots + # # Remove intracluster rotations + # k_i = unpack_gradient(k_i, norb) + # for ci in clusters + # k_i[ci.orb_list, ci.orb_list] .= 0 + # end + # k_i = pack_gradient(k_i, norb) + #end # # Compute energy and gradient # - e_i, g_i, d1_i = step!(k_i) + #e_i, g_i, d1_i = step!(k_i) + e_i, g_i, d1_i, h_i= step!(k_i) g_i = reshape(g_i, (norb2,1)) k_i = reshape(k_i, (norb2,1)) d_i = d1_i @@ -1516,7 +1679,7 @@ function cmf_oo_diis(ints_in::InCoreInts{T}, clusters::Vector{MOCluster}, fspace if norm(g_i) < tol_oo - @printf("*ooCMF Iter: %4i Total= %16.12f G= %12.2e #SS: %4s\n", i, e_i, norm(g_i), nss) + @printf("*ooCMF Iter: %4i Total= %16.15f G= %12.3e #SS: %4s\n", i, e_i, norm(g_i), nss) break end @@ -1526,7 +1689,7 @@ function cmf_oo_diis(ints_in::InCoreInts{T}, clusters::Vector{MOCluster}, fspace - @printf(" ooCMF Iter: %4i Total= %16.12f G= %12.2e #SS: %4s\n", i, e_i, norm(g_i), nss) + @printf(" ooCMF Iter: %4i Total= %16.15f G= %12.3e #SS: %4s\n", i, e_i, norm(g_i), nss) end U = exp(unpack_gradient(k_i,norb)) @@ -1534,3 +1697,44 @@ function cmf_oo_diis(ints_in::InCoreInts{T}, clusters::Vector{MOCluster}, fspace return e, U, d1 #=}}}=# end + + +function projection_vector(ansatze::Vector{<:Ansatz}, norb; gradient=false) + n_dim = norb*(norb-1)÷2#={{{=# + + + + tmp_mat = Matrix(1I, n_dim, n_dim) + + shift = 0 + invar = Vector{Tuple{Int,Int}}() + + for cluster in ansatze + tmp = ActiveSpaceSolvers.invariant_orbital_rotations(cluster) + for j in 1:length(tmp) + tmp[j] = tmp[j].+shift + end + shift += cluster.no + append!(invar, tmp) + end + + if gradient == true + return invar + end + + fci = ActiveSpaceSolvers.FCIAnsatz(norb, 0, 0) #dummie FCI anstaz to generate all pairs + full_list = ActiveSpaceSolvers.invariant_orbital_rotations(fci) + + keep_list = [] + for (a,b) in enumerate(full_list) + if b in invar + continue + else + push!(keep_list, a) + end + end + + proj_vec = tmp_mat[:,keep_list] + + return proj_vec#=}}}=# +end