Fix travil

.travis.yml

@@ -27,7 +27,7 @@ python:
 script: 
   - ./configure  --production ./config/gfortran.cfg
   - source ./quantum_package.rc
-  - qp_module.py install Full_CI Hartree_Fock
+  - qp_install_module.py install Full_CI Hartree_Fock CAS_SD MRCC_CASSD
   - ninja
   - cd ocaml ; make ; cd -
   - cd testing_no_regression ; ./unit_test.py

config/gfortran.cfg

@@ -10,7 +10,7 @@
 #
 #
 [COMMON]
-FC           : gfortran -g -ffree-line-length-none -I .
+FC           : gfortran -g -ffree-line-length-none -I . -static-libgcc 
 LAPACK_LIB   : -llapack -lblas
 IRPF90       : irpf90
 IRPF90_FLAGS : --ninja --align=32
@@ -22,7 +22,7 @@ IRPF90_FLAGS : --ninja --align=32
 # 0 : Deactivate
 #
 [OPTION]
-MODE    : OPT        ; [ OPT | PROFILE | DEBUG ] : Chooses the section below
+MODE    : DEBUG        ; [ OPT | PROFILE | DEBUG ] : Chooses the section below
 CACHE   : 1          ; Enable cache_compile.py
 OPENMP  : 1          ; Append OpenMP flags
 

configure

@@ -33,8 +33,8 @@ import sys
 from os.path import join
 
 if not any(i in ["--production", "--development"] for i in sys.argv):
-    print __doc__
-    sys.exit()
+    sys.argv += ["--development"]
+
 if len(sys.argv) != 3:
     print __doc__
     sys.exit()
@@ -46,6 +46,7 @@ if len(sys.argv) != 3:
 
 QP_ROOT = os.getcwd()
 QP_ROOT_BIN = join(QP_ROOT, "bin")
+QP_ROOT_LIB = join(QP_ROOT, "lib")
 QP_ROOT_INSTALL = join(QP_ROOT, "install")
 
 os.environ["PATH"] = os.environ["PATH"] + ":" + QP_ROOT_BIN
@@ -62,8 +63,11 @@ d_dependency = {
     "resultsFile": ["python"],
     "emsl": ["python"],
     "gcc": [],
+    "g++": [],
+    "zeromq" : [ "g++" ],
+    "f77zmq" : [ "zeromq", "python" ],
     "python": [],
-    "ninja": ["gcc", "python"],
+    "ninja": ["g++", "python"],
     "make": [],
     "p_graphviz": ["python"]
 }
@@ -92,12 +96,12 @@ curl = Info(
 zlib = Info(
     url='http://zlib.net/zlib-1.2.8.tar.gz',
     description=' zlib',
-    default_path=join(QP_ROOT_INSTALL, "zlib"))
+    default_path=join(QP_ROOT_LIB, "libz.a"))
 
-path = Info(
+patch = Info(
     url='ftp://ftp.gnu.org/gnu/patch/patch-2.7.5.tar.gz',
-    description=' path',
-    default_path=join(QP_ROOT, "lib", "libz.a"))
+    description=' patch',
+    default_path=join(QP_ROOT_BIN, "patch"))
 
 irpf90 = Info(
     url='{head}/LCPQ/irpf90/{tail}'.format(**path_github),
@@ -116,12 +120,11 @@ resultsFile = Info(
 
 ninja = Info(
     url='{head}/martine/ninja/{tail}'.format(**path_github),
-    description=' nina',
+    description=' ninja',
     default_path=join(QP_ROOT_BIN, "ninja"))
 
 emsl = Info(
-    url='{head}/LCPQ/EMSL_Basis_Set_Exchange_Local/{tail}'.format(**
-                                                                  path_github),
+    url='{head}/LCPQ/EMSL_Basis_Set_Exchange_Local/{tail}'.format(**path_github),
     description=' EMSL basis set library',
     default_path=join(QP_ROOT_INSTALL, "emsl"))
 
@@ -130,6 +133,16 @@ ezfio = Info(
     description=' EZFIO',
     default_path=join(QP_ROOT_INSTALL, "EZFIO"))
 
+zeromq = Info(
+    url='http://download.zeromq.org/zeromq-4.1.3.tar.gz',
+    description=' ZeroMQ',
+    default_path=join(QP_ROOT_LIB, "libzmq.a"))
+
+f77zmq = Info(
+    url='{head}/zeromq/f77_zmq/{tail}'.format(**path_github),
+    description=' F77-ZeroMQ',
+    default_path=join(QP_ROOT_LIB, "libf77zmq.a"))
+
 p_graphviz = Info(
     url='https://github.com/xflr6/graphviz/archive/master.tar.gz',
     description=' Python library for graphviz',
@@ -137,8 +150,9 @@ p_graphviz = Info(
 
 d_info = dict()
 
-for m in ["ocaml", "m4", "curl", "zlib", "path", "irpf90", "docopt",
-          "resultsFile", "ninja", "emsl", "ezfio", "p_graphviz"]:
+for m in ["ocaml", "m4", "curl", "zlib", "patch", "irpf90", "docopt",
+          "resultsFile", "ninja", "emsl", "ezfio", "p_graphviz",
+          "zeromq", "f77zmq" ]:
     exec ("d_info['{0}']={0}".format(m))
 
 
@@ -189,8 +203,7 @@ def check_output(*popenargs, **kwargs):
 
 def checking(d_dependency):
     """
-    For each key in d_dependency check if it
-    is avalabie or not
+    For each key in d_dependency check if it is avalabie 
     """
 
     def check_python():
@@ -260,7 +273,7 @@ def checking(d_dependency):
     l_installed = dict()
     l_needed = []
 
-    # Check all the other
+    # Check all the others
     length = max(map(len, d_dependency))
 
     for i in d_dependency.keys():
@@ -275,7 +288,7 @@ def checking(d_dependency):
             l_needed.append(i)
     print ""
 
-    # Expend the need_stuff for all the genealogy
+    # Expand the needed stuff for all the genealogy
     l_install_descendant = get_list_descendant(d_dependency, l_installed,
                                                l_needed)
 
@@ -328,7 +341,7 @@ _|_ | | _>  |_ (_| | | (_|  |_ | (_) | |
     d_print = {
         "install_ninja": "Install ninja...",
         "build": "Creating build.ninja...",
-        "install": "Installing the dependencies with Ninja..."
+        "install": "Installing the dependencies using Ninja..."
     }
 
     length = max(map(len, d_print.values()))
@@ -372,7 +385,7 @@ _|_ | | _>  |_ (_| | | (_|  |_ | (_) | |
         descr = d_info[need].description
         default_path = d_info[need].default_path
 
-        # Build to dowload
+        # Build to download
         l_build += ["build {0}: download".format(archive_path),
                     "   url = {0}".format(url), "   descr = {0}".format(descr),
                     ""]
@@ -404,7 +417,16 @@ _|_ | | _>  |_ (_| | | (_|  |_ | (_) | |
         path_ninja = find_path("ninja", l_installed)
         subprocess.check_call("cd install ;{0}".format(path_ninja), shell=True)
     except:
-        raise
+        prefix = os.path.join('install', '_build')
+        for filename in os.listdir(prefix):
+            if filename.endswith(".log"):
+                with open( os.path.join(prefix,filename) ,'r') as f:
+                    print "\n\n"
+                    print "=-=-=-=-=-=- %s =-=-=-=-=-=-" %(filename)
+                    print f.read()
+                    print "=-=-=-=-=-=-=-=-=-=-=-=-=-=-\n\n"
+        print "Error in installation of dependencies"
+        sys.exit(1)
     else:
         print r"""
  _________
@@ -437,11 +459,12 @@ def create_ninja_and_rc(l_installed):
     print str_info("qp_root"),
     python_path = [join(QP_ROOT, "scripts"), join(QP_ROOT, "install")]
 
-    l_python = [join(QP_ROOT, "scripts")]
+    l_python = [join("${QP_ROOT}", "scripts")]
     for dir_ in python_path:
         for folder in os.listdir(dir_):
             path = join(dir_, folder)
             if os.path.isdir(path):
+                path = path.replace(QP_ROOT,"${QP_ROOT}")
                 l_python.append(path)
 
     path_ezfio = find_path('ezfio', l_installed, var_for_qp_root=True)
@@ -450,9 +473,9 @@ def create_ninja_and_rc(l_installed):
 
     l_rc = [
         'export QP_ROOT={0}'.format(QP_ROOT),
-        'export QP_EZFIO={0}'.format(path_ezfio),
-        'export IRPF90={0}'.format(path_irpf90),
-        'export NINJA={0}'.format(path_ninja),
+        'export QP_EZFIO={0}'.format(path_ezfio.replace(QP_ROOT,"${QP_ROOT}")),
+        'export IRPF90={0}'.format(path_irpf90.replace(QP_ROOT,"${QP_ROOT}")),
+        'export NINJA={0}'.format(path_ninja.replace(QP_ROOT,"${QP_ROOT}")),
         'export QP_PYTHON={0}'.format(":".join(l_python)), "",
         'export PYTHONPATH="${QP_EZFIO}/Python":"${QP_PYTHON}":"${PYTHONPATH}"',
         'export PATH="${QP_PYTHON}":"${QP_ROOT}"/bin:"${QP_ROOT}"/ocaml:"${PATH}"',
@@ -526,3 +549,4 @@ if __name__ == '__main__':
     create_ninja_and_rc(l_installed)
 
     recommendation()
+

install/scripts/install_f77zmq.sh

@@ -0,0 +1,23 @@
+#!/bin/bash -x
+
+TARGET=f77zmq
+
+function _install()
+{
+  cd ..
+  QP_ROOT=$PWD
+  cd -
+  export C_INCLUDE_PATH="${C_INCLUDE_PATH}":"${QP_ROOT}"/lib
+  set -e
+  set -u
+  export ZMQ_H="${QP_ROOT}"/lib/zmq.h
+  cd "${BUILD}"
+  make -j 8                        || exit 1
+  mv libf77zmq.a "${QP_ROOT}"/lib  || exit 1
+  mv libf77zmq.so "${QP_ROOT}"/lib || exit 1
+  cp f77_zmq.h "${QP_ROOT}"/src/ZMQ/
+  cd -
+  return 0
+}
+
+source scripts/build.sh

install/scripts/install_zeromq.sh

@@ -0,0 +1,27 @@
+#!/bin/bash -x
+
+TARGET=zeromq
+
+function _install()
+{
+  cd ..
+  QP_ROOT=$PWD
+  cd -
+  export C_INCLUDE_PATH="${C_INCLUDE_PATH}":./
+  set -e
+  set -u
+  ORIG=$(pwd)
+  cd "${BUILD}"
+  ./configure --without-libsodium || exit 1
+  make -j 8 || exit 1
+  rm -f -- "${QP_ROOT}"/lib/libzmq.a "${QP_ROOT}"/lib/libzmq.so "${QP_ROOT}"/lib/libzmq.so.5
+  cp .libs/libzmq.a  "${QP_ROOT}"/lib 
+  cp .libs/libzmq.so "${QP_ROOT}"/lib/libzmq.so.5
+  cp include/{zmq.h,zmq_utils.h} "${QP_ROOT}"/lib
+  cd "${QP_ROOT}"/lib
+  ln -s libzmq.so.5 libzmq.so
+  cd ${ORIG}
+  return 0
+}
+
+source scripts/build.sh

ocaml/Atom.ml

@@ -8,8 +8,8 @@ type t =
   coord   : Point3d.t ;
 } with sexp
 
-(** Read xyz coordinates of the atom with unit u *)
-let of_string u s =
+(** Read xyz coordinates of the atom *)
+let of_string ~units s =
   let buffer = s
   |> String.split ~on:' '
   |> List.filter ~f:(fun x -> x <> "")
@@ -18,21 +18,21 @@ let of_string u s =
   | [ name; charge; x; y; z ] ->
     { element = Element.of_string name ;
       charge  = Charge.of_string  charge ;
-      coord   = Point3d.of_string u (String.concat [x; y; z] ~sep:" ")
+      coord   = Point3d.of_string ~units (String.concat [x; y; z] ~sep:" ")
     }
   | [ name; x; y; z ] ->
     let e = Element.of_string name in
     { element = e ;
       charge  = Element.to_charge e;
-      coord   = Point3d.of_string u (String.concat [x; y; z] ~sep:" ")
+      coord   = Point3d.of_string ~units (String.concat [x; y; z] ~sep:" ")
     }
   | _ -> raise (AtomError s)
 ;;
   
-let to_string u a =
+let to_string ~units a =
   [ Element.to_string a.element ;
     Charge.to_string  a.charge ;
-    Point3d.to_string u a.coord ]
+    Point3d.to_string ~units a.coord ]
   |> String.concat ~sep:"   "
 ;;
 

ocaml/Atom.mli

@@ -5,5 +5,5 @@ type t = { element : Element.t; charge : Charge.t; coord : Point3d.t; }
 val t_of_sexp : Sexplib.Sexp.t -> t
 val sexp_of_t : t -> Sexplib.Sexp.t
 
-val of_string : Units.units -> string -> t
-val to_string : Units.units -> t -> string
+val of_string : units:Units.units -> string -> t
+val to_string : units:Units.units -> t -> string

ocaml/Element.ml

@@ -1,4 +1,5 @@
-open Core.Std;;
+open Core.Std
+open Qptypes
 
 exception ElementError of string
 
@@ -8,7 +9,7 @@ type t =
 |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
-with sexp;;
+with sexp
 
 let of_string x = 
   match (String.capitalize (String.lowercase x)) with
@@ -50,7 +51,7 @@ let of_string x =
 |  "Br"  |  "Bromine"     ->  Br
 |  "Kr"  |  "Krypton"     ->  Kr
 | x -> raise (ElementError ("Element "^x^" unknown"))
-;;
+
 
 let to_string = function
 | X   -> "X"
@@ -90,7 +91,7 @@ let to_string = function
 | Se  -> "Se"
 | Br  -> "Br"
 | Kr  -> "Kr"
-;;
+
 
 let to_long_string = function
 | X   -> "Dummy"
@@ -130,7 +131,7 @@ let to_long_string = function
 | Se  -> "Selenium"
 | Br  -> "Bromine"
 | Kr  -> "Krypton"
-;;
+
 
 let to_charge c = 
   let result = match c with
@@ -172,7 +173,7 @@ let to_charge c =
   | Br  -> 35
   | Kr  -> 36
   in Charge.of_int result
-;;
+
 
 let of_charge c = match (Charge.to_int c) with
 |  0   ->  X   
@@ -213,5 +214,134 @@ let of_charge c = match (Charge.to_int c) with
 |  35  ->  Br  
 |  36  ->  Kr  
 | x -> raise (ElementError ("Element of charge "^(string_of_int x)^" unknown"))
-;;
+
+
+let covalent_radius x = 
+  let result = function
+  | X   -> 0.
+  | H   -> 0.37
+  | He  -> 0.70
+  | Li  -> 1.23
+  | Be  -> 0.89
+  | B   -> 0.90
+  | C   -> 0.85
+  | N   -> 0.74
+  | O   -> 0.74
+  | F   -> 0.72
+  | Ne  -> 0.70
+  | Na  -> 1.00
+  | Mg  -> 1.36
+  | Al  -> 1.25
+  | Si  -> 1.17
+  | P   -> 1.10
+  | S   -> 1.10
+  | Cl  -> 0.99
+  | Ar  -> 0.70
+  | K   -> 2.03
+  | Ca  -> 1.74
+  | Sc  -> 1.44
+  | Ti  -> 1.32
+  | V   -> 1.22
+  | Cr  -> 0.00
+  | Mn  -> 1.16
+  | Fe  -> 0.00
+  | Co  -> 1.15
+  | Ni  -> 1.17
+  | Cu  -> 1.25
+  | Zn  -> 1.25
+  | Ga  -> 1.20
+  | Ge  -> 1.21
+  | As  -> 1.16
+  | Se  -> 0.70
+  | Br  -> 1.24
+  | Kr  -> 1.91
+  in
+  Units.angstrom_to_bohr *. (result x)
+  |> Positive_float.of_float
+
+let vdw_radius x = 
+  let result = function
+  | X    -> 0.
+  | H    -> 1.20
+  | He  ->  1.70
+  | Li  ->  1.70
+  | Be  ->  1.70
+  | B   ->  1.70
+  | C   ->  1.70
+  | N   ->  1.55
+  | O   ->  1.52
+  | F   ->  1.47
+  | Ne  ->  1.70
+  | Na  ->  1.70
+  | Mg  ->  1.70
+  | Al  ->  1.94
+  | Si  ->  2.10
+  | P   ->  1.80
+  | S   ->  1.80
+  | Cl  ->  1.75
+  | Ar  ->  1.70
+  | K   ->  1.70
+  | Ca  ->  1.70
+  | Sc  ->  1.70
+  | Ti  ->  1.70
+  | V   ->  1.98
+  | Cr  ->  1.94
+  | Mn  ->  1.93
+  | Fe  ->  1.93
+  | Co  ->  1.92
+  | Ni  ->  1.70
+  | Cu  ->  1.70
+  | Zn  ->  1.70
+  | Ga  ->  2.02
+  | Ge  ->  1.70
+  | As  ->  1.96
+  | Se  ->  1.70
+  | Br  ->  2.10
+  | Kr  ->  1.70
+  in
+  Units.angstrom_to_bohr *. (result x)
+  |> Positive_float.of_float
+
+let mass x =
+  let result = function
+  | X   ->    0.
+  | H   ->    1.0079
+  | He  ->    4.00260
+  | Li  ->    6.941
+  | Be  ->    9.01218
+  | B   ->   10.81
+  | C   ->   12.011
+  | N   ->   14.0067
+  | O   ->   15.9994
+  | F   ->   18.998403
+  | Ne  ->   20.179
+  | Na  ->   22.98977
+  | Mg  ->   24.305
+  | Al  ->   26.98154
+  | Si  ->   28.0855
+  | P   ->   30.97376
+  | S   ->   32.06
+  | Cl  ->   35.453
+  | Ar  ->   39.948
+  | K   ->   39.0983
+  | Ca  ->   40.08
+  | Sc  ->   44.9559
+  | Ti  ->   47.90
+  | V   ->   50.9415
+  | Cr  ->   51.996
+  | Mn  ->   54.9380
+  | Fe  ->   55.9332
+  | Co  ->   58.9332
+  | Ni  ->   58.70
+  | Cu  ->   63.546
+  | Zn  ->   65.38
+  | Ga  ->   69.72
+  | Ge  ->   72.59
+  | As  ->   74.9216
+  | Se  ->   78.96
+  | Br  ->   79.904
+  | Kr  ->   83.80
+  in
+  result x
+  |> Positive_float.of_float
 

ocaml/Element.mli

@@ -13,6 +13,8 @@ val of_string : string -> t
 val to_string : t -> string
 val to_long_string : t -> string
 
-(** get the positive charge *)
+(** Properties *)
 val to_charge : t -> Charge.t
 val of_charge : Charge.t -> t
+val covalent_radius : t -> Qptypes.Positive_float.t
+val vdw_radius : t -> Qptypes.Positive_float.t

ocaml/Input_determinants_by_hand.ml

@@ -7,8 +7,6 @@ module Determinants_by_hand : sig
     { n_int                  : N_int_number.t;
       bit_kind               : Bit_kind.t;
       n_det                  : Det_number.t;
-      n_states               : States_number.t;
-      n_states_diag          : States_number.t;
       expected_s2            : Positive_float.t;
       psi_coef               : Det_coef.t array;
       psi_det                : Determinant.t array;
@@ -23,8 +21,6 @@ end = struct
     { n_int                  : N_int_number.t;
       bit_kind               : Bit_kind.t;
       n_det                  : Det_number.t;
-      n_states               : States_number.t;
-      n_states_diag          : States_number.t;
       expected_s2            : Positive_float.t;
       psi_coef               : Det_coef.t array;
       psi_det                : Determinant.t array;
@@ -146,11 +142,12 @@ end = struct
     |> Array.map ~f:Det_coef.of_float
   ;;
 
-  let write_psi_coef ~n_det ~n_states c =
+  let write_psi_coef ~n_det c =
     let n_det = Det_number.to_int n_det
     and c = Array.to_list c
             |> List.map ~f:Det_coef.to_float
-    and n_states = States_number.to_int n_states
+    and n_states = 
+      read_n_states () |> States_number.to_int
     in
     Ezfio.ezfio_array_of_list ~rank:2 ~dim:[| n_det ; n_states |] ~data:c
     |> Ezfio.set_determinants_psi_coef 
@@ -214,8 +211,6 @@ end = struct
       { n_int                  = read_n_int ()                ;
         bit_kind               = read_bit_kind ()             ;
         n_det                  = read_n_det ()                ;
-        n_states               = read_n_states ()             ;
-        n_states_diag          = read_n_states_diag ()        ;
         expected_s2            = read_expected_s2 ()          ;
         psi_coef               = read_psi_coef ()             ;
         psi_det                = read_psi_det ()              ;
@@ -227,8 +222,6 @@ end = struct
   let write { n_int                ;
               bit_kind             ;
               n_det                ;
-              n_states             ;
-              n_states_diag        ;
               expected_s2          ;
               psi_coef             ;
               psi_det              ;
@@ -236,10 +229,8 @@ end = struct
      write_n_int n_int ;
      write_bit_kind bit_kind;
      write_n_det n_det;
-     write_n_states n_states;
-     write_n_states_diag ~n_states:n_states n_states_diag;
      write_expected_s2 expected_s2;
-     write_psi_coef ~n_det:n_det psi_coef ~n_states:n_states;
+     write_psi_coef ~n_det:n_det psi_coef ;
      write_psi_det ~n_int:n_int ~n_det:n_det psi_det;
   ;;
 
@@ -249,7 +240,7 @@ end = struct
     let mo_tot_num = MO_number.of_int mo_tot_num ~max:mo_tot_num in
     let det_text = 
       let nstates =
-        States_number.to_int b.n_states
+        read_n_states () |> States_number.to_int
       and ndet =
         Det_number.to_int b.n_det
       in
@@ -284,12 +275,6 @@ If true, input the expected value of S^2 ::
 
   expected_s2 = %s
 
-Number of requested states, and number of states used for the
-Davidson diagonalization ::
-
-  n_states      = %s
-  n_states_diag = %s
-
 Number of determinants ::
 
   n_det = %s
@@ -299,8 +284,6 @@ Determinants ::
 %s
 "
      (b.expected_s2   |> Positive_float.to_string)
-     (b.n_states      |> States_number.to_string)
-     (b.n_states_diag |> States_number.to_string)
      (b.n_det         |> Det_number.to_string)
      det_text
      |> Rst_string.of_string
@@ -313,8 +296,6 @@ Determinants ::
 n_int                  = %s
 bit_kind               = %s
 n_det                  = %s
-n_states               = %s
-n_states_diag          = %s
 expected_s2            = %s
 psi_coef               = %s
 psi_det                = %s
@@ -322,8 +303,6 @@ psi_det                = %s
      (b.n_int         |> N_int_number.to_string)
      (b.bit_kind      |> Bit_kind.to_string)
      (b.n_det         |> Det_number.to_string)
-     (b.n_states      |> States_number.to_string)
-     (b.n_states_diag |> States_number.to_string)
      (b.expected_s2   |> Positive_float.to_string)
      (b.psi_coef  |> Array.to_list |> List.map ~f:Det_coef.to_string
       |> String.concat ~sep:", ")

ocaml/Input_nuclei.ml

@@ -147,7 +147,7 @@ nucl_coord       = %s
     (b.nucl_charge |> Array.to_list |> List.map
       ~f:(Charge.to_string) |> String.concat ~sep:", " )
     (b.nucl_coord  |> Array.to_list |> List.map
-      ~f:(Point3d.to_string Units.Bohr) |> String.concat ~sep:"\n" )
+      ~f:(Point3d.to_string ~units:Units.Bohr) |> String.concat ~sep:"\n" )
   ;;
 
 
@@ -161,7 +161,7 @@ nucl_coord       = %s
          Printf.sprintf "  %-3s  %d   %s"
           (b.nucl_label.(i)  |> Element.to_string)
           (b.nucl_charge.(i) |> Charge.to_int )
-          (b.nucl_coord.(i)  |> Point3d.to_string Units.Angstrom) )
+          (b.nucl_coord.(i)  |> Point3d.to_string ~units:Units.Angstrom) )
       ) |> String.concat ~sep:"\n"
      in
      Printf.sprintf "

ocaml/Molecule.ml

@@ -10,27 +10,36 @@ type t = {
 } with sexp
 
 let get_charge { nuclei  ; elec_alpha ; elec_beta } =
-  let result = (Elec_alpha_number.to_int elec_alpha) +
-     (Elec_beta_number.to_int elec_beta) in
+  let result =
+     (Elec_alpha_number.to_int elec_alpha) +
+     (Elec_beta_number.to_int elec_beta)
+  in
   let rec nucl_charge = function
   | a::rest -> (Charge.to_float a.Atom.charge) +. nucl_charge rest
   | [] -> 0.
   in
   Charge.of_float (nucl_charge nuclei  -. (Float.of_int result))
-;;
+
 
 let get_multiplicity m = 
-  let elec_alpha = m.elec_alpha in
+  let elec_alpha =
+     m.elec_alpha
+  in
   Multiplicity.of_alpha_beta elec_alpha m.elec_beta
-;;
+
 
 let get_nucl_num m =
-  let nmax = (List.length m.nuclei) in
+  let nmax =
+    List.length m.nuclei 
+  in
   Nucl_number.of_int nmax ~max:nmax
-;;
+
 
 let name m = 
-  let cm = Charge.to_int (get_charge m) in
+  let cm =
+    get_charge m
+    |> Charge.to_int
+  in
   let c = 
      match cm with
      | 0 -> ""
@@ -39,8 +48,12 @@ let name m =
      | i when i>1 -> Printf.sprintf " (%d+)" i
      | i -> Printf.sprintf " (%d-)" (-i)
   in
-  let mult = Multiplicity.to_string (get_multiplicity m) in
-  let { nuclei  ; elec_alpha ; elec_beta } = m in
+  let mult =
+    get_multiplicity m
+    |> Multiplicity.to_string
+  in
+  let { nuclei  ; elec_alpha ; elec_beta } = m
+  in
   let rec build_list accu = function
   | a::rest ->
       begin
@@ -53,7 +66,9 @@ let name m =
   in
   let rec build_name accu = function
   | (a, n)::rest ->
-    let a = Element.to_string a in
+    let a =
+      Element.to_string a
+    in
     begin
       match n with 
       | 1 -> build_name (a::accu) rest
@@ -64,19 +79,25 @@ let name m =
     end
   | [] -> accu
   in
-  let result = build_list [] nuclei |> build_name [c ; ", " ; mult]
+  let result =
+     build_list [] nuclei |> build_name [c ; ", " ; mult]
   in
   String.concat (result)
-;;
+
 
 let to_string m =
-  let { nuclei  ; elec_alpha ; elec_beta } = m in
-  let n = List.length nuclei in
-  let title = name m in
-  [ Int.to_string n ; title ] @ (List.map ~f:(fun x -> Atom.to_string
-  Units.Angstrom x) nuclei)
+  let { nuclei  ; elec_alpha ; elec_beta } = m
+  in
+  let n =
+    List.length nuclei
+  in
+  let title =
+     name m
+  in
+  [ Int.to_string n ; title ] @ 
+   (List.map ~f:(fun x -> Atom.to_string Units.Angstrom x) nuclei)
   |> String.concat ~sep:"\n"
-;;
+
 
 let of_xyz_string
     ?(charge=(Charge.of_int 0)) ?(multiplicity=(Multiplicity.of_int 1))
@@ -94,7 +115,9 @@ let of_xyz_string
       ) + 1 - (Charge.to_int charge)
       |> Elec_number.of_int 
   in
-  let (na,nb) = Multiplicity.to_alpha_beta ne multiplicity in
+  let (na,nb) =
+     Multiplicity.to_alpha_beta ne multiplicity
+  in
   let result = 
   { nuclei = l ;
     elec_alpha = na ;
@@ -109,7 +132,7 @@ let of_xyz_string
      raise (MultiplicityError msg);
   else () ;
   result
-;;
+
 
 
 let of_xyz_file
@@ -121,8 +144,33 @@ let of_xyz_file
   let (_,buffer) = String.lsplit2_exn buffer ~on:'\n' in
   of_xyz_string ~charge:charge ~multiplicity:multiplicity 
     ~units:units buffer
-;;
 
-include To_md5;;
+
+
+let distance_matrix molecule =
+  let coord =
+    molecule.nuclei
+    |> List.map ~f:(fun x -> x.Atom.coord)
+    |> Array.of_list
+  in
+  let n = 
+    Array.length coord
+  in
+  let result =
+    Array.make_matrix ~dimx:n  ~dimy:n 0.
+  in
+  for i = 0 to (n-1)
+  do
+    for j = 0 to (n-1)
+    do
+      result.(i).(j) <- Point3d.distance coord.(i) coord.(j)
+    done;
+  done;
+  result
+  
+  
+  
+
+include To_md5
 let to_md5 = to_md5 sexp_of_t
-;;
+

ocaml/Molecule.mli

@@ -34,5 +34,9 @@ val of_xyz_string :
   ?multiplicity:Multiplicity.t ->
   ?units:Units.units -> string -> t
 
+(** Creates the distance matrix between all the atoms *)
+val distance_matrix :
+  t -> (float array) array
+
 (** Computes the MD5 hash *)
 val to_md5 : t -> Qptypes.MD5.t

ocaml/Point3d.ml

@@ -7,9 +7,16 @@ type t = {
   z : float ;
 } with sexp
 
+let of_tuple ~units (x,y,z) =
+  let f = match units with 
+    | Units.Bohr     -> 1.
+    | Units.Angstrom -> Units.angstrom_to_bohr
+  in
+  { x = x *. f ; y = y *. f ; z = z *. f }
+
 (** Read x y z coordinates in string s with units u *)
-let of_string u s =
-  let f = match u with 
+let of_string ~units s =
+  let f = match units with 
     | Units.Bohr     -> 1.
     | Units.Angstrom -> Units.angstrom_to_bohr 
   in
@@ -22,7 +29,6 @@ let of_string u s =
   { x = l.(0) *. f ;
     y = l.(1) *. f ;
     z = l.(2) *. f }
-;;
            
 
 let distance2 p1 p2 =
@@ -30,17 +36,18 @@ let distance2 p1 p2 =
   and { x=x2 ; y=y2 ; z=z2 } = p2 in
   (x2-.x1)*.(x2-.x1) +. (y2-.y1)*.(y2-.y1) +. (z2-.z1)*.(z2-.z1) 
   |> Positive_float.of_float
-;;
 
-let distance p1 p2 = sqrt (Positive_float.to_float (distance2 p1 p2))
-;;
 
-let to_string u p =
-  let f = match u with 
+let distance p1 p2 =
+  sqrt (Positive_float.to_float (distance2 p1 p2))
+
+
+let to_string ~units p =
+  let f = match units with 
     | Units.Bohr     -> 1.
     | Units.Angstrom -> Units.bohr_to_angstrom
   in
   let { x=x ; y=y ; z=z } = p in
   Printf.sprintf "%16.8f %16.8f %16.8f" (x*.f) (y*.f) (z*.f)
-;;
+
   

ocaml/Point3d.mli

@@ -4,11 +4,14 @@ type t =
   z : float;
 } with sexp
 
+(** Create from a tuple of floats *)
+val of_tuple  : units:Units.units -> float*float*float -> t
+
 (** Create from an xyz string *)
-val of_string : Units.units -> string -> t
+val of_string : units:Units.units -> string -> t
 
 (** Convert to a string for printing *)
-val to_string : Units.units -> t -> string
+val to_string : units:Units.units -> t -> string
 
 (** Computes the squared distance between 2 points *)
 val distance2 : t -> t -> Qptypes.Positive_float.t

ocaml/Qputils.ml

@@ -12,7 +12,6 @@ let rec transpose = function
 ;;
 *)
 
-
 let input_to_sexp s =
     let result = 
       String.split_lines s 

ocaml/qp_create_ezfio_from_xyz.ml

@@ -11,64 +11,139 @@ let spec =
      ~doc:"string Name of basis set."
   +> flag "c" (optional_with_default 0 int)
      ~doc:"int Total charge of the molecule. Default is 0."
+  +> flag "d" (optional_with_default 0. float)
+     ~doc:"float Add dummy atoms. x * (covalent radii of the atoms)"
   +> flag "m" (optional_with_default 1 int)
      ~doc:"int Spin multiplicity (2S+1) of the molecule. Default is 1."
   +> flag "p" no_arg
      ~doc:" Using pseudopotentials"
-  +> anon ("xyz_file" %: string)
-;;
+  +> anon ("xyz_file" %: file )
 
-let run ?o b c m p xyz_file =
+
+let dummy_centers ~threshold ~molecule ~nuclei =
+  let d =
+    Molecule.distance_matrix molecule
+  in
+  let n =
+    Array.length d
+  in
+  let nuclei =
+    Array.of_list nuclei
+  in
+  let rec aux accu = function
+  | (-1,_) -> accu
+  | (i,-1) -> aux accu (i-1,i-1)
+  | (i,j) when (i>j) ->
+    let new_accu = 
+      let x,y =
+          Element.covalent_radius (nuclei.(i)).Atom.element |> Positive_float.to_float, 
+          Element.covalent_radius (nuclei.(j)).Atom.element |> Positive_float.to_float
+      in 
+      let r =
+        ( x +. y ) *. threshold
+      in
+      if d.(i).(j) < r then
+        (i,x,j,y,d.(i).(j)) :: accu
+      else
+        accu
+    in aux new_accu (i,j-1)
+  | (i,j) when (i=j) -> aux accu (i,j-1)
+  | _ -> assert false
+  in
+  aux [] (n-1,n-1)
+  |> List.map ~f:(fun (i,x,j,y,r) ->
+    let f = 
+      x /. (x +. y)
+    in
+    let u = 
+      Point3d.of_tuple ~units:Units.Bohr
+      ( nuclei.(i).Atom.coord.Point3d.x +. 
+          (nuclei.(j).Atom.coord.Point3d.x -. nuclei.(i).Atom.coord.Point3d.x) *. f, 
+        nuclei.(i).Atom.coord.Point3d.y +. 
+          (nuclei.(j).Atom.coord.Point3d.y -. nuclei.(i).Atom.coord.Point3d.y) *. f, 
+        nuclei.(i).Atom.coord.Point3d.z +. 
+          (nuclei.(j).Atom.coord.Point3d.z -. nuclei.(i).Atom.coord.Point3d.z) *. f) 
+    in
+    Atom.{ element = Element.X ; charge = Charge.of_int 0 ; coord = u }
+    )
+     
+    
+let list_basis () =
+  let basis_list = 
+    Qpackage.root ^ "/install/emsl/EMSL_api.py list_basis" 
+    |> Unix.open_process_in 
+    |> In_channel.input_lines
+    |> List.map ~f:(fun x -> 
+       match String.split x ~on:'\'' with
+       | [] -> ""
+       | a :: [] 
+       | _ :: a :: _ -> String.strip a
+      )
+  in 
+  List.sort basis_list ~cmp:String.ascending
+  |> String.concat ~sep:"\t" 
+
+
+let run ?o b c d m p xyz_file =
 
   (* Read molecule *)
   let molecule =
     (Molecule.of_xyz_file xyz_file ~charge:(Charge.of_int c)
       ~multiplicity:(Multiplicity.of_int m) )
   in
-  let nuclei = molecule.Molecule.nuclei in
+  let dummy =
+    dummy_centers ~threshold:d ~molecule ~nuclei:molecule.Molecule.nuclei
+  in
+(*
+  List.iter dummy ~f:(fun x ->
+    Printf.printf "%s\n" (Atom.to_string ~units:Units.Angstrom x)
+  );
+*)
+  let nuclei =
+    molecule.Molecule.nuclei @ dummy
+  in
+
+
+  let basis_table =
+     Hashtbl.Poly.create ()
+  in
 
-  let basis_table = Hashtbl.Poly.create () in
   (* Open basis set channels *)
   let basis_channel element =
-    let key = Element.to_string element in
+    let key =
+      Element.to_string element
+    in
     match Hashtbl.find basis_table key with
     | Some in_channel -> 
         in_channel
     | None -> 
-        begin
-         Printf.printf "%s is not defined in basis %s.\nEnter alternate basis : %!"
+         let msg = 
+           Printf.sprintf "%s is not defined in basis %s.%!"
            (Element.to_long_string element) b ;
-          let bas =
-             match In_channel.input_line stdin with
-             | Some line -> String.strip line |> String.lowercase
-             | None -> failwith "Aborted"
          in
-          let new_channel = In_channel.create 
-            (Qpackage.root ^ "/data/basis/" ^ bas)
-          in
-          Hashtbl.add_exn basis_table ~key:key ~data:new_channel;
-          new_channel
-        end
+         failwith msg
   in
   
   let temp_filename =
     Filename.temp_file "qp_create_" ".basis"
   in
-  let rec build_basis = function
-  | [] -> ()
-  | elem_and_basis_name :: rest -> 
-    begin
-      match (String.lsplit2 ~on:':' elem_and_basis_name) with
-      | None -> (* Principal basis *)
-        let basis = elem_and_basis_name in
+  let () = 
+    Sys.remove temp_filename
+  in
+
+  let fetch_channel basis =
     let command =
       if (p) then  
         Qpackage.root ^ "/scripts/get_basis.sh \"" ^ temp_filename 
-              ^ "\" \"" ^ basis ^"\" pseudo"
+          ^ "." ^ basis ^ "\" \"" ^ basis ^"\" pseudo"
       else
         Qpackage.root ^ "/scripts/get_basis.sh \"" ^ temp_filename 
-              ^ "\" \"" ^ basis ^"\""
+          ^ "." ^ basis ^ "\" \"" ^ basis ^"\""
     in
+    match Sys.is_file basis with
+    | `Yes -> 
+          In_channel.create basis
+    | _ -> 
       begin
         let filename = 
           Unix.open_process_in command
@@ -79,8 +154,26 @@ let run ?o b c m p xyz_file =
           In_channel.create filename
         in
         Unix.unlink filename;
+        new_channel
+      end
+  in
+
+  let rec build_basis = function
+  | [] -> ()
+  | elem_and_basis_name :: rest -> 
+    begin
+      match (String.lsplit2 ~on:':' elem_and_basis_name) with
+      | None -> (* Principal basis *)
+        begin
+          let basis = 
+            elem_and_basis_name 
+          in
+          let new_channel =
+            fetch_channel basis
+          in
           List.iter nuclei ~f:(fun elem->
-            let key = Element.to_string elem.Atom.element
+            let key = 
+              Element.to_string elem.Atom.element
             in
             match Hashtbl.add basis_table ~key:key ~data:new_channel with
             | `Ok -> ()
@@ -89,25 +182,18 @@ let run ?o b c m p xyz_file =
         end
       | Some (key, basis) -> (*Aux basis *)
         begin
-          let elem  = Element.of_string key
-          and basis = String.lowercase basis
+          let elem  = 
+            Element.of_string key
+          and basis =
+            String.lowercase basis
           in
-          let key = Element.to_string elem
-          in
-          let command =
-              Qpackage.root ^ "/scripts/get_basis.sh \"" ^ temp_filename ^
-                "\" \"" ^ basis ^ "\" " ^ key
-          in
-          begin
-            let filename = 
-              Unix.open_process_in command
-              |> In_channel.input_all
-              |> String.strip
+          let key = 
+             Element.to_string elem
           in
           let new_channel =
-              In_channel.create filename
+            fetch_channel basis
           in
-            Unix.unlink filename;
+          begin
             match Hashtbl.add basis_table ~key:key ~data:new_channel with
             | `Ok -> ()
             | `Duplicate -> failwith ("Duplicate definition of basis for "^(Element.to_long_string elem))
@@ -164,28 +250,30 @@ let run ?o b c m p xyz_file =
   (* Write Basis set *)
   let basis =
 
-    let nmax = Nucl_number.get_max () in
+    let nmax =
+       Nucl_number.get_max ()
+    in
     let rec do_work (accu:(Atom.t*Nucl_number.t) list) (n:int) = function
     | [] -> accu
     | e::tail ->
-      let new_accu = (e,(Nucl_number.of_int ~max:nmax n))::accu in
+      let new_accu =
+        (e,(Nucl_number.of_int ~max:nmax n))::accu
+      in
       do_work new_accu (n+1) tail
     in
     let result = do_work [] 1  nuclei
     |> List.rev
     |> List.map ~f:(fun (x,i) ->
        try 
-         Basis.read_element (basis_channel x.Atom.element) i x.Atom.element
+         let e =
+           match x.Atom.element with
+           | Element.X -> Element.H
+           | e -> e
+         in
+         Basis.read_element (basis_channel x.Atom.element) i e
        with
-       | End_of_file -> 
-         begin
-           let alt_channel = basis_channel x.Atom.element in
-           try
-             Basis.read_element alt_channel i x.Atom.element 
-           with
-           End_of_file -> failwith
-             ("Element "^(Element.to_string x.Atom.element)^" not found")
-         end
+       | End_of_file -> failwith
+             ("Element "^(Element.to_string x.Atom.element)^" not found in basis set.")
        ) 
     |> List.concat
     in
@@ -264,28 +352,37 @@ let run ?o b c m p xyz_file =
   | None -> failwith "Error in basis"
   | Some x -> Input.Ao_basis.write x
  
-;;
+
 
 let command = 
     Command.basic 
     ~summary: "Quantum Package command"
     ~readme:(fun () -> "
-Creates an EZFIO directory from a standard xyz file.
-The basis set is defined as a single string if all the
-atoms are taken from the same basis set, otherwise specific
-elements can be defined as follows:
+
+=== Available basis sets ===
+
+" ^ (list_basis ()) ^ "
+
+============================
+
+Creates an EZFIO directory from a standard xyz file.  The basis set is defined
+as a single string if all the atoms are taken from the same basis set,
+otherwise specific elements can be defined as follows:
 
  -b \"cc-pcvdz | H:cc-pvdz | C:6-31g\"
 
+If a file with the same name as the basis set exists, this file will be read.
+Otherwise, the basis set is obtained from the database.
+
 " )
     spec
-    (fun o b c m p xyz_file () ->
-       run ?o b c m p xyz_file )
-;;
+    (fun o b c d m p xyz_file () ->
+       run ?o b c d m p xyz_file )
+
 
 let () =
     Command.run command
-;;
+
 
 
 

ocaml/qp_edit.ml

@@ -0,0 +1,308 @@
+open Qputils;;
+open Qptypes;;
+open Core.Std;;
+
+(** Interactive editing of the input.
+
+WARNING 
+This file is autogenerad by
+`${QP_ROOT}/script/ezfio_interface/ei_handler.py`
+*)
+
+
+(** Keywords used to define input sections *)
+type keyword = 
+| Ao_basis
+| Determinants_by_hand
+| Electrons
+| Mo_basis
+| Nuclei
+| Determinants
+| Hartree_fock
+| Integrals_bielec
+| Perturbation
+| Properties
+| Pseudo
+;;
+
+
+let keyword_to_string = function
+| Ao_basis              -> "AO basis"
+| Determinants_by_hand  -> "Determinants_by_hand"
+| Electrons             -> "Electrons"
+| Mo_basis              -> "MO basis"
+| Nuclei                -> "Molecule"
+| Determinants -> "Determinants"
+| Hartree_fock -> "Hartree_fock"
+| Integrals_bielec -> "Integrals_bielec"
+| Perturbation -> "Perturbation"
+| Properties -> "Properties"
+| Pseudo -> "Pseudo"
+;;
+
+
+
+(** Create the header of the temporary file *)
+let file_header filename = 
+  Printf.sprintf  "
+==================================================================
+                       Quantum Package
+==================================================================
+
+Editing file `%s`
+
+" filename
+;;
+  
+
+(** Creates the header of a section *)
+let make_header kw =
+  let s = keyword_to_string kw in
+  let l = String.length s in
+  "\n\n"^s^"\n"^(String.init l ~f:(fun _ -> '='))^"\n\n"
+;;
+
+
+(** Returns the rst string of section [s] *)
+let get s = 
+  let header = (make_header s) in
+  let f (read,to_rst) = 
+    match read () with
+    | Some text -> header ^ (Rst_string.to_string (to_rst text))
+    | None      -> ""
+  in
+  let rst = 
+    try
+      begin
+         let open Input in
+         match s with
+         | Mo_basis ->
+           f Mo_basis.(read, to_rst)
+         | Electrons ->
+           f Electrons.(read, to_rst)
+         | Nuclei ->
+           f Nuclei.(read, to_rst)
+         | Ao_basis ->
+           f Ao_basis.(read, to_rst)
+         | Determinants_by_hand ->
+           f Determinants_by_hand.(read, to_rst)
+         | Determinants ->
+           f Determinants.(read, to_rst)
+         | Hartree_fock ->
+           f Hartree_fock.(read, to_rst)
+         | Integrals_bielec ->
+           f Integrals_bielec.(read, to_rst)
+         | Perturbation ->
+           f Perturbation.(read, to_rst)
+         | Properties ->
+           f Properties.(read, to_rst)
+         | Pseudo ->
+           f Pseudo.(read, to_rst)
+      end
+    with
+    | Sys_error msg -> (Printf.eprintf "Info: %s\n%!" msg ; "")
+  in 
+  rst
+;;
+
+
+(** Applies the changes from the string [str] corresponding to section [s] *)
+let set str s = 
+  let header = (make_header s) in
+  match String.substr_index ~pos:0 ~pattern:header str with
+  | None -> ()
+  | Some idx -> 
+    begin
+      let index_begin = idx + (String.length header) in
+      let index_end   = 
+        match ( String.substr_index ~pos:(index_begin+(String.length header)+1)
+          ~pattern:"==" str) with
+          | Some i -> i
+          | None -> String.length str
+      in
+      let l = index_end - index_begin in
+      let str = String.sub ~pos:index_begin ~len:l str
+      |> Rst_string.of_string
+      in
+      let write (of_rst,w) s =
+        try
+        match of_rst str with
+        | Some data -> w data
+        | None -> ()
+        with
+        | _ -> (Printf.eprintf "Info: Read error in %s\n%!"
+               (keyword_to_string s); ignore (of_rst str) )
+      in
+      let open Input in
+        match s with
+        | Determinants -> write Determinants.(of_rst, write) s
+        | Hartree_fock -> write Hartree_fock.(of_rst, write) s
+        | Integrals_bielec -> write Integrals_bielec.(of_rst, write) s
+        | Perturbation -> write Perturbation.(of_rst, write) s
+        | Properties -> write Properties.(of_rst, write) s
+        | Pseudo -> write Pseudo.(of_rst, write) s
+        | Electrons        -> write Electrons.(of_rst, write) s
+        | Determinants_by_hand     -> write Determinants_by_hand.(of_rst, write) s
+        | Nuclei           -> write Nuclei.(of_rst, write) s
+        | Ao_basis         -> () (* TODO *)
+        | Mo_basis         -> () (* TODO *)
+    end 
+;;
+
+
+(** Creates the temporary file for interactive editing *)
+let create_temp_file ezfio_filename fields =
+  let temp_filename  = Filename.temp_file "qp_edit_" ".rst" in
+  begin
+      Out_channel.with_file temp_filename ~f:(fun out_channel ->
+        (file_header ezfio_filename) :: (List.map ~f:get fields) 
+        |> String.concat ~sep:"\n" 
+        |> Out_channel.output_string out_channel 
+      )
+  end
+  ; temp_filename
+;;
+
+
+
+let run check_only ezfio_filename =
+
+  (* Open EZFIO *)
+  if (not (Sys.file_exists_exn ezfio_filename)) then
+    failwith (ezfio_filename^" does not exists");
+
+  Ezfio.set_file ezfio_filename;
+
+  (*
+  let output = (file_header ezfio_filename) :: (
+    List.map ~f:get [
+      Ao_basis ; 
+      Mo_basis ; 
+    ])
+   in
+  String.concat output
+  |> print_string
+  *)
+  
+  let tasks = [
+      Nuclei ;
+      Ao_basis;
+      Electrons ;
+      Determinants ; 
+      Hartree_fock ; 
+      Integrals_bielec ; 
+      Perturbation ; 
+      Properties ; 
+      Pseudo ; 
+      Mo_basis;
+      Determinants_by_hand ;
+  ]
+  in
+
+  (* Create the temp file *)
+  let temp_filename = create_temp_file ezfio_filename tasks in
+
+  (* Open the temp file with external editor *)
+  let editor = 
+    match Sys.getenv "EDITOR" with
+    | Some editor -> editor
+    | None -> "vi"
+  in
+
+  match check_only with
+  | true  -> ()
+  | false -> 
+    Printf.sprintf "%s %s ; tput sgr0 2> /dev/null" editor temp_filename 
+    |> Sys.command_exn 
+  ;
+
+  (* Re-read the temp file *)
+  let temp_string  =
+    In_channel.with_file temp_filename ~f:(fun in_channel ->
+      In_channel.input_all in_channel) 
+  in
+  List.iter ~f:(fun x -> set temp_string x) tasks;
+
+  (* Remove temp_file *)
+  Sys.remove temp_filename;
+;;
+
+
+(** Create a backup file in case of an exception *)
+let create_backup ezfio_filename =
+  Printf.sprintf "
+    rm -f %s/backup.tgz ;
+    tar -zcf .backup.tgz %s && mv .backup.tgz %s/backup.tgz
+  "
+    ezfio_filename ezfio_filename ezfio_filename
+  |> Sys.command_exn
+;;
+
+
+(** Restore the backup file when an exception occuprs *)
+let restore_backup ezfio_filename =
+  Printf.sprintf "tar -zxf %s/backup.tgz"
+    ezfio_filename 
+  |> Sys.command_exn
+;;
+
+
+let spec =
+  let open Command.Spec in
+  empty 
+  +> flag "-c" no_arg
+     ~doc:"Checks the input data"
+(*
+  +> flag "o" (optional string)
+     ~doc:"Prints output data"
+*)
+  +> anon ("ezfio_file" %: string)
+;;
+
+let command = 
+    Command.basic 
+    ~summary: "Quantum Package command"
+    ~readme:(fun () ->
+      "
+Edit input data
+      ")
+    spec
+(*    (fun i o ezfio_file () -> *)
+    (*fun ezfio_file () -> 
+       try 
+           run ezfio_file
+       with
+       | _ msg -> print_string ("\n\nError\n\n"^msg^"\n\n")
+    *)
+    (fun c ezfio_file () ->
+       try
+          run c ezfio_file ;
+          (* create_backup ezfio_file; *)
+       with
+       | Failure exc 
+       | Invalid_argument exc as e -> 
+         begin
+           Printf.eprintf "=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-\n\n";
+           Printf.eprintf "%s\n\n" exc;
+           Printf.eprintf "=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-\n\n";
+           (* restore_backup ezfio_file; *) 
+           raise e
+         end
+       | Assert_failure (file, line, ch) as e ->
+         begin
+           Printf.eprintf "=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-\n\n";
+           Printf.eprintf "Assert error in file $QP_ROOT/ocaml/%s, line %d, character %d\n\n" file line ch;
+           Printf.eprintf "=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-\n\n";
+           (* restore_backup ezfio_file; *)
+           raise e
+         end
+    )
+;;
+
+let () =
+  Command.run command;
+  exit 0
+;;
+
+
+

ocaml/qp_run.ml

@@ -42,7 +42,7 @@ let run exe ezfio_file =
 let spec = 
   let open Command.Spec in
   empty
-  +> anon ("exectuable" %: string)
+  +> anon ("executable" %: string)
   +> anon ("ezfio_file" %: string)
 ;;
 

ocaml/qptypes_generator.ml

@@ -83,6 +83,12 @@ let input_data = "
   assert (x >= 0.) ;
   assert (x <= 1.) ;
 
+* Energy : float
+  assert (x <=0.) ;
+
+* S2 : float
+  assert (x >=0.) ;
+
 * PT2_energy : float
   assert (x >=0.) ;
 

ocaml/test_molecule.ml

@@ -39,6 +39,15 @@ H        0.54386314      0.00000000     -0.92559535
   let m = Molecule.of_xyz_file "c2h6.xyz"  in
   print_string (Molecule.to_string m);
 
+  print_string "\nDistance matrix\n";
+  print_string   "---------------\n";
+  let d = 
+    Molecule.distance_matrix m
+  in
+  Array.iter d ~f:(fun x ->
+    Array.iter x ~f:(fun y -> Printf.printf "%12.8f " y);
+    print_newline ();
+  )  
 ;;
 
 test_molecule ();;

plugins/CASSCF/EZFIO.cfg

@@ -0,0 +1,10 @@
+[energy]
+type: double precision
+doc: "Calculated CAS-SCF energy"
+interface: ezfio
+
+[energy_pt2]
+type: double precision
+doc: "Calculated selected CAS-SCF energy with PT2 correction"
+interface: ezfio
+

plugins/CASSCF/H_apply.irp.f

@@ -0,0 +1,39 @@
+use bitmasks
+BEGIN_SHELL [ /usr/bin/env python ]
+from generate_h_apply import *
+
+s = H_apply("CAS_SD")
+print s
+
+s = H_apply("CAS_SD_selected_no_skip")
+s.set_selection_pt2("epstein_nesbet_2x2")
+s.unset_skip()
+print s
+
+s = H_apply("CAS_SD_selected")
+s.set_selection_pt2("epstein_nesbet_2x2")
+print s
+
+s = H_apply("CAS_SD_PT2")
+s.set_perturbation("epstein_nesbet_2x2")
+print s
+
+
+s = H_apply("CAS_S",do_double_exc=False)
+print s
+
+s = H_apply("CAS_S_selected_no_skip",do_double_exc=False)
+s.set_selection_pt2("epstein_nesbet_2x2")
+s.unset_skip()
+print s
+
+s = H_apply("CAS_S_selected",do_double_exc=False)
+s.set_selection_pt2("epstein_nesbet_2x2")
+print s
+
+s = H_apply("CAS_S_PT2",do_double_exc=False)
+s.set_perturbation("epstein_nesbet_2x2")
+print s
+
+END_SHELL
+

plugins/CASSCF/NEEDED_CHILDREN_MODULES

@@ -0,0 +1 @@
+Generators_CAS Perturbation Selectors_full

plugins/CASSCF/README.rst

@@ -0,0 +1,20 @@
+======
+CASSCF
+======
+
+This module is not a "real" CAS-SCF. It is an orbital optimization step done by :
+
+1) Doing the CAS+SD
+2) Taking one-electron density matrix
+3) Cancelling all active-active rotations
+4) Finding the order which matches with the input MOs
+
+
+Needed Modules
+==============
+.. Do not edit this section It was auto-generated
+.. by the `update_README.py` script.
+Documentation
+=============
+.. Do not edit this section It was auto-generated
+.. by the `update_README.py` script.

plugins/CASSCF/casscf.irp.f

@@ -0,0 +1,220 @@
+program casscf
+  implicit none
+  BEGIN_DOC
+! Optimize MOs and CI coefficients of the CAS
+  END_DOC
+  
+  double precision, allocatable  :: pt2(:), norm_pert(:), H_pert_diag(:)
+  integer(bit_kind), allocatable  :: generators_bitmask_save(:,:,:,:) 
+
+  integer                        :: degree, N_generators_bitmask_save, N_det_ci
+  double precision               :: E_old, E_CI
+  double precision               :: selection_criterion_save, selection_criterion_min_save
+  
+  integer                        :: N_det_old
+  integer                        :: i, j, k, l
+  integer                        :: i_bit, j_bit, i_int, j_int
+  integer(bit_kind), allocatable :: bit_tmp(:), cas_bm(:)
+  character*(64) :: label
+
+  allocate( pt2(N_states), norm_pert(N_states),H_pert_diag(N_states) )
+  allocate( generators_bitmask_save(N_int,2,6,N_generators_bitmask) )
+  allocate( bit_tmp(N_int), cas_bm(N_int) )
+
+  PROVIDE N_det_cas
+  N_det_old = 0
+  pt2 = 1.d0
+  E_CI = 1.d0
+  E_old = 0.d0
+  diag_algorithm = "Lapack"
+  selection_criterion_save = selection_criterion
+  selection_criterion_min_save = selection_criterion_min
+
+
+  cas_bm = 0_bit_kind
+  do i=1,N_cas_bitmask
+    do j=1,N_int
+      cas_bm(j) = ior(cas_bm(j), cas_bitmask(j,1,i))
+      cas_bm(j) = ior(cas_bm(j), cas_bitmask(j,2,i))
+    enddo
+  enddo
+
+  ! Save CAS-SD bitmask
+  generators_bitmask_save = generators_bitmask
+  N_generators_bitmask_save = N_generators_bitmask
+
+  ! Set the CAS bitmask 
+  do i=1,6
+    generators_bitmask(:,:,i,:) = cas_bitmask
+  enddo
+  N_generators_bitmask = N_cas_bitmask
+  SOFT_TOUCH generators_bitmask N_generators_bitmask
+
+
+  ! If the number of dets already in the file is larger than the requested
+  ! number of determinants, truncate the wf
+  if (N_det > N_det_max) then
+    call diagonalize_CI
+    call save_wavefunction
+    psi_det = psi_det_sorted
+    psi_coef = psi_coef_sorted
+    N_det = N_det_max
+    soft_touch N_det psi_det psi_coef
+    call diagonalize_CI
+    call save_wavefunction
+    print *,  'N_det    = ', N_det
+    print *,  'N_states = ', N_states
+    print *,  'PT2      = ', pt2
+    print *,  'E        = ', CI_energy
+    print *,  'E+PT2    = ', CI_energy+pt2
+    print *,  '-----'
+  endif
+
+  ! Start MCSCF iteration
+
+  ! CAS-CI
+  ! ------
+
+  E_old = E_CI
+
+  ! Reset the selection criterion
+  selection_criterion     = selection_criterion_save
+  selection_criterion_min = selection_criterion_min_save
+  SOFT_TOUCH selection_criterion_min selection_criterion selection_criterion_factor
+
+  ! Set the CAS bitmask 
+  do i=1,6
+    generators_bitmask(:,:,i,:) = cas_bitmask
+  enddo
+  N_generators_bitmask = N_cas_bitmask
+  SOFT_TOUCH generators_bitmask N_generators_bitmask
+
+  do while (N_det < N_det_max.and.maxval(abs(pt2(1:N_states))) > pt2_max)
+    N_det_old = N_det
+    call H_apply_CAS_SD_selected_no_skip(pt2, norm_pert, H_pert_diag,  N_states)
+
+    PROVIDE  psi_coef
+    PROVIDE  psi_det
+    PROVIDE  psi_det_sorted
+
+    if (N_det > N_det_max) then
+      psi_det = psi_det_sorted
+      psi_coef = psi_coef_sorted
+      N_det = N_det_max
+      soft_touch N_det psi_det psi_coef
+    endif
+    call diagonalize_CI
+    call save_wavefunction
+    print *,  '======'
+    print *,  'CAS-CI'
+    print *,  '======'
+    print *,  ''
+    print *,  'N_det      = ', N_det
+    print *,  'N_states   = ', N_states
+    print *,  'PT2        = ', pt2
+    print *,  'E(CAS)     = ', CI_energy
+    print *,  'E(CAS)+PT2 = ', CI_energy+pt2
+    print *,  '-----'
+    print *,  ''
+    E_CI = sum(CI_energy(1:N_states)+pt2(1:N_states))/dble(N_states)
+
+    call ezfio_set_casscf_energy(CI_energy(1))
+    if (abort_all) then
+      exit
+    endif
+    if (N_det == N_det_old) then
+      exit
+    endif
+
+  enddo
+
+  ! Super-CI
+  ! --------
+
+  selection_criterion_min = 1.d-12
+  selection_criterion     = 1.d-12
+
+  ! Set the CAS bitmask
+  generators_bitmask = generators_bitmask_save
+  N_generators_bitmask = N_generators_bitmask_save
+  SOFT_TOUCH generators_bitmask N_generators_bitmask selection_criterion selection_criterion_min selection_criterion_factor
+
+  N_det_ci = N_det
+
+  call H_apply_CAS_SD_selected(pt2, norm_pert, H_pert_diag,  N_states)
+  
+
+  do i=1,mo_tot_num
+   i_int = ishft(i-1,-bit_kind_shift)+1
+   i_bit = j-ishft(i_int-1,bit_kind_shift)-1
+   bit_tmp(:) = 0_bit_kind
+   bit_tmp(i_int) = ibset(0_bit_kind,i_bit)
+   if (iand(bit_tmp(i_int), cas_bm(i_int)) == 0_bit_kind) then
+     ! Not a CAS MO
+     cycle
+   endif
+
+   do j=1,mo_tot_num
+     if (j == i) then
+       cycle
+     endif
+     j_int = ishft(j-1,-bit_kind_shift)+1
+     j_bit = j-ishft(j_int-1,bit_kind_shift)-1
+     bit_tmp(:) = 0_bit_kind
+     bit_tmp(j_int) = ibset(0_bit_kind,j_bit)
+     if (iand(bit_tmp(j_int), cas_bm(j_int)) == 0_bit_kind) then
+       ! Not a CAS MO
+       cycle
+     endif
+     ! Now, both i and j are MOs of the CAS. De-couple them in the DM
+     one_body_dm_mo(i,j) = 0.d0
+   enddo
+
+  enddo
+
+  SOFT_TOUCH one_body_dm_mo
+
+  double precision               :: mx, ov
+  double precision, allocatable  :: mo_coef_old(:,:)
+  integer, allocatable           :: iorder(:)
+  logical, allocatable           :: selected(:)
+  allocate( mo_coef_old(size(mo_coef,1), size(mo_coef,2)), iorder(mo_tot_num), selected(mo_tot_num) )
+  mo_coef_old = mo_coef
+  label = "Canonical"
+  call mo_as_eigvectors_of_mo_matrix(one_body_dm_mo,size(one_body_dm_mo,1),size(one_body_dm_mo,2),label,-1)
+  selected = .False.
+  do j=1,mo_tot_num
+   mx = -1.d0
+   iorder(j) = j
+   do i=1,mo_tot_num
+    if (selected(i)) then
+      cycle
+    endif
+    ov = 0.d0
+    do l=1,ao_num
+     do k=1,ao_num
+      ov = ov + mo_coef_old(k,j) * ao_overlap(k,l) * mo_coef(l,i)
+     enddo
+    enddo
+    ov= dabs(ov)
+    if (ov > mx) then
+      mx = ov
+      iorder(j) = i
+    endif
+   enddo
+   selected( iorder(j) ) = .True.
+  enddo
+  mo_coef_old = mo_coef
+  do i=1,mo_tot_num
+    mo_coef(:,i) = mo_coef_old(:,iorder(i))
+  enddo
+
+  call save_mos
+
+  call write_double(6,E_CI,"Energy(CAS)")
+
+  deallocate( mo_coef_old )
+  deallocate( pt2, norm_pert,H_pert_diag )
+  deallocate( generators_bitmask_save )
+  deallocate( bit_tmp, cas_bm, iorder )
+end

plugins/CAS_SD/H_apply.irp.f

@@ -5,7 +5,6 @@ from generate_h_apply import *
 s = H_apply("CAS_SD")
 print s
 
-
 s = H_apply("CAS_SD_selected_no_skip")
 s.set_selection_pt2("epstein_nesbet_2x2")
 s.unset_skip()
@@ -19,5 +18,22 @@ s = H_apply("CAS_SD_PT2")
 s.set_perturbation("epstein_nesbet_2x2")
 print s
 
+
+s = H_apply("CAS_S",do_double_exc=False)
+print s
+
+s = H_apply("CAS_S_selected_no_skip",do_double_exc=False)
+s.set_selection_pt2("epstein_nesbet_2x2")
+s.unset_skip()
+print s
+
+s = H_apply("CAS_S_selected",do_double_exc=False)
+s.set_selection_pt2("epstein_nesbet_2x2")
+print s
+
+s = H_apply("CAS_S_PT2",do_double_exc=False)
+s.set_perturbation("epstein_nesbet_2x2")
+print s
+
 END_SHELL
 

plugins/CAS_SD/README.rst

@@ -118,69 +118,101 @@ Documentation
   Undocumented
 
 
-`h_apply_cas_sd <http://github.com/LCPQ/quantum_package/tree/master/plugins/CAS_SD/H_apply.irp.f_shell_22#L414>`_
+h_apply_cas_sd
   Calls H_apply on the HF determinant and selects all connected single and double
   excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
 
 
-`h_apply_cas_sd_diexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/CAS_SD/H_apply.irp.f_shell_22#L1>`_
+h_apply_cas_sd_diexc
+  Undocumented
+
+
+h_apply_cas_sd_diexcorg
   Generate all double excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_cas_sd_monoexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/CAS_SD/H_apply.irp.f_shell_22#L269>`_
+h_apply_cas_sd_diexcp
+  Undocumented
+
+
+h_apply_cas_sd_monoexc
   Generate all single excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_cas_sd_pt2 <http://github.com/LCPQ/quantum_package/tree/master/plugins/CAS_SD/H_apply.irp.f_shell_22#L2610>`_
+h_apply_cas_sd_pt2
   Calls H_apply on the HF determinant and selects all connected single and double
   excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
 
 
-`h_apply_cas_sd_pt2_diexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/CAS_SD/H_apply.irp.f_shell_22#L2118>`_
+h_apply_cas_sd_pt2_diexc
+  Undocumented
+
+
+h_apply_cas_sd_pt2_diexcorg
   Generate all double excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_cas_sd_pt2_monoexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/CAS_SD/H_apply.irp.f_shell_22#L2427>`_
+h_apply_cas_sd_pt2_diexcp
+  Undocumented
+
+
+h_apply_cas_sd_pt2_monoexc
   Generate all single excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_cas_sd_selected <http://github.com/LCPQ/quantum_package/tree/master/plugins/CAS_SD/H_apply.irp.f_shell_22#L1872>`_
+h_apply_cas_sd_selected
   Calls H_apply on the HF determinant and selects all connected single and double
   excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
 
 
-`h_apply_cas_sd_selected_diexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/CAS_SD/H_apply.irp.f_shell_22#L1346>`_
+h_apply_cas_sd_selected_diexc
+  Undocumented
+
+
+h_apply_cas_sd_selected_diexcorg
   Generate all double excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_cas_sd_selected_monoexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/CAS_SD/H_apply.irp.f_shell_22#L1675>`_
+h_apply_cas_sd_selected_diexcp
+  Undocumented
+
+
+h_apply_cas_sd_selected_monoexc
   Generate all single excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_cas_sd_selected_no_skip <http://github.com/LCPQ/quantum_package/tree/master/plugins/CAS_SD/H_apply.irp.f_shell_22#L1128>`_
+h_apply_cas_sd_selected_no_skip
   Calls H_apply on the HF determinant and selects all connected single and double
   excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
 
 
-`h_apply_cas_sd_selected_no_skip_diexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/CAS_SD/H_apply.irp.f_shell_22#L602>`_
+h_apply_cas_sd_selected_no_skip_diexc
+  Undocumented
+
+
+h_apply_cas_sd_selected_no_skip_diexcorg
   Generate all double excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_cas_sd_selected_no_skip_monoexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/CAS_SD/H_apply.irp.f_shell_22#L931>`_
+h_apply_cas_sd_selected_no_skip_diexcp
+  Undocumented
+
+
+h_apply_cas_sd_selected_no_skip_monoexc
   Generate all single excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.

plugins/CAS_SD/cas_s.irp.f

@@ -0,0 +1,95 @@
+program full_ci
+  implicit none
+  integer                        :: i,k
+  integer                        :: N_det_old
+
+  
+  double precision, allocatable  :: pt2(:), norm_pert(:), H_pert_diag(:)
+  integer                        :: N_st, degree
+  N_st = N_states
+  allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st))
+  character*(64)                 :: perturbation
+  PROVIDE N_det_cas
+  
+  N_det_old = 0
+  pt2 = 1.d0
+  diag_algorithm = "Lapack"
+  if (N_det > N_det_max) then
+    call diagonalize_CI
+    call save_wavefunction
+    psi_det = psi_det_sorted
+    psi_coef = psi_coef_sorted
+    N_det = N_det_max
+    soft_touch N_det psi_det psi_coef
+    call diagonalize_CI
+    call save_wavefunction
+    print *,  'N_det    = ', N_det
+    print *,  'N_states = ', N_states
+    print *,  'PT2      = ', pt2
+    print *,  'E        = ', CI_energy
+    print *,  'E+PT2    = ', CI_energy+pt2
+    print *,  '-----'
+  endif
+
+  do while (N_det < N_det_max.and.maxval(abs(pt2(1:N_st))) > pt2_max)
+    N_det_old = N_det
+    call H_apply_CAS_S(pt2, norm_pert, H_pert_diag,  N_st)
+
+    PROVIDE  psi_coef
+    PROVIDE  psi_det
+    PROVIDE  psi_det_sorted
+
+    if (N_det > N_det_max) then
+       psi_det = psi_det_sorted
+       psi_coef = psi_coef_sorted
+       N_det = N_det_max
+       soft_touch N_det psi_det psi_coef
+    endif
+    call diagonalize_CI
+    call save_wavefunction
+    print *,  'N_det    = ', N_det
+    print *,  'N_states = ', N_states
+    print *,  'PT2      = ', pt2
+    print *,  'E        = ', CI_energy
+    print *,  'E+PT2    = ', CI_energy+pt2
+    print *,  '-----'
+    call ezfio_set_cas_sd_energy(CI_energy(1))
+    if (abort_all) then
+      exit
+    endif
+    if (N_det == N_det_old) then
+      exit
+    endif
+  enddo
+  call diagonalize_CI
+
+   if(do_pt2_end)then
+    print*,'Last iteration only to compute the PT2'
+    threshold_selectors = 1.d0
+    threshold_generators = 0.999d0
+    call H_apply_CAS_S_PT2(pt2, norm_pert, H_pert_diag,  N_st)
+
+    print *,  'Final step'
+    print *,  'N_det    = ', N_det
+    print *,  'N_states = ', N_states
+    print *,  'PT2      = ', pt2
+    print *,  'E        = ', CI_energy
+    print *,  'E+PT2    = ', CI_energy+pt2
+    print *,  '-----'
+    call ezfio_set_cas_sd_energy_pt2(CI_energy(1)+pt2(1))
+   endif
+
+
+  integer :: exc_max, degree_min
+  exc_max = 0
+  print *,  'CAS determinants : ', N_det_cas
+  do i=1,min(N_det_cas,10)
+    do k=i,N_det_cas
+      call get_excitation_degree(psi_cas(1,1,k),psi_cas(1,1,i),degree,N_int)
+      exc_max = max(exc_max,degree)
+    enddo
+    call debug_det(psi_cas(1,1,i),N_int)
+    print *,  ''
+  enddo
+  print *,  'Max excitation degree in the CAS :', exc_max
+end

plugins/CAS_SD/cas_s_selected.irp.f

@@ -0,0 +1,89 @@
+program full_ci
+  implicit none
+  integer                        :: i,k
+
+  
+  double precision, allocatable  :: pt2(:), norm_pert(:), H_pert_diag(:)
+  integer                        :: N_st, degree
+  N_st = N_states
+  allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st))
+  character*(64)                 :: perturbation
+  PROVIDE N_det_cas
+  
+  pt2 = 1.d0
+  diag_algorithm = "Lapack"
+  if (N_det > N_det_max) then
+    call diagonalize_CI
+    call save_wavefunction
+    psi_det = psi_det_sorted
+    psi_coef = psi_coef_sorted
+    N_det = N_det_max
+    soft_touch N_det psi_det psi_coef
+    call diagonalize_CI
+    call save_wavefunction
+    print *,  'N_det    = ', N_det
+    print *,  'N_states = ', N_states
+    print *,  'PT2      = ', pt2
+    print *,  'E        = ', CI_energy
+    print *,  'E+PT2    = ', CI_energy+pt2
+    print *,  '-----'
+  endif
+
+  do while (N_det < N_det_max.and.maxval(abs(pt2(1:N_st))) > pt2_max)
+    call H_apply_CAS_S_selected(pt2, norm_pert, H_pert_diag,  N_st)
+
+    PROVIDE  psi_coef
+    PROVIDE  psi_det
+    PROVIDE  psi_det_sorted
+
+    if (N_det > N_det_max) then
+       psi_det = psi_det_sorted
+       psi_coef = psi_coef_sorted
+       N_det = N_det_max
+       soft_touch N_det psi_det psi_coef
+    endif
+    call diagonalize_CI
+    call save_wavefunction
+    print *,  'N_det    = ', N_det
+    print *,  'N_states = ', N_states
+    print *,  'PT2      = ', pt2
+    print *,  'E        = ', CI_energy
+    print *,  'E+PT2    = ', CI_energy+pt2
+    print *,  '-----'
+    call ezfio_set_cas_sd_energy(CI_energy(1))
+    if (abort_all) then
+      exit
+    endif
+  enddo
+  call diagonalize_CI
+
+   if(do_pt2_end)then
+    print*,'Last iteration only to compute the PT2'
+    threshold_selectors = 1.d0
+    threshold_generators = 0.999d0
+    call H_apply_CAS_S_PT2(pt2, norm_pert, H_pert_diag,  N_st)
+
+    print *,  'Final step'
+    print *,  'N_det    = ', N_det
+    print *,  'N_states = ', N_states
+    print *,  'PT2      = ', pt2
+    print *,  'E        = ', CI_energy
+    print *,  'E+PT2    = ', CI_energy+pt2
+    print *,  '-----'
+    call ezfio_set_cas_sd_energy_pt2(CI_energy(1)+pt2(1))
+   endif
+
+
+  integer :: exc_max, degree_min
+  exc_max = 0
+  print *,  'CAS determinants : ', N_det_cas
+  do i=1,min(N_det_cas,10)
+    do k=i,N_det_cas
+      call get_excitation_degree(psi_cas(1,1,k),psi_cas(1,1,i),degree,N_int)
+      exc_max = max(exc_max,degree)
+    enddo
+    call debug_det(psi_cas(1,1,i),N_int)
+    print *,  ''
+  enddo
+  print *,  'Max excitation degree in the CAS :', exc_max
+end

plugins/CID/README.rst

@@ -203,3 +203,48 @@ Documentation
   particles.
   Assume N_int is already provided.
 
+Needed Modules
+==============
+.. Do not edit this section It was auto-generated
+.. by the `update_README.py` script.
+
+
+.. image:: tree_dependency.png
+
+* `Selectors_full <http://github.com/LCPQ/quantum_package/tree/master/plugins/Selectors_full>`_
+* `SingleRefMethod <http://github.com/LCPQ/quantum_package/tree/master/plugins/SingleRefMethod>`_
+
+Documentation
+=============
+.. Do not edit this section It was auto-generated
+.. by the `update_README.py` script.
+
+
+`cid <http://github.com/LCPQ/quantum_package/tree/master/plugins/CID/cid_lapack.irp.f#L1>`_
+  Undocumented
+
+
+h_apply_cid
+  Calls H_apply on the HF determinant and selects all connected single and double
+  excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
+
+
+h_apply_cid_diexc
+  Undocumented
+
+
+h_apply_cid_diexcorg
+  Generate all double excitations of key_in using the bit masks of holes and
+  particles.
+  Assume N_int is already provided.
+
+
+h_apply_cid_diexcp
+  Undocumented
+
+
+h_apply_cid_monoexc
+  Generate all single excitations of key_in using the bit masks of holes and
+  particles.
+  Assume N_int is already provided.
+

plugins/CISD/.gitignore

@@ -20,6 +20,7 @@ Pseudo
 Selectors_full
 SingleRefMethod
 Utils
+cisd
 cisd_lapack
 ezfio_interface.irp.f
 irpf90.make

plugins/CISD/EZFIO.cfg

@@ -0,0 +1,10 @@
+[energy]
+type: double precision
+doc: Variational CISD energy
+interface: ezfio
+
+[energy_pt2]
+type: double precision
+doc: Estimated CISD energy (including PT2)
+interface: ezfio
+

plugins/CISD/README.rst

@@ -59,22 +59,26 @@ Documentation
 .. by the `update_README.py` script.
 
 
-`cisd <http://github.com/LCPQ/quantum_package/tree/master/plugins/CISD/cisd_lapack.irp.f#L1>`_
-  Undocumented
-
-
-`h_apply_cisd <http://github.com/LCPQ/quantum_package/tree/master/plugins/CISD/H_apply.irp.f_shell_8#L414>`_
+h_apply_cisd
   Calls H_apply on the HF determinant and selects all connected single and double
   excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
 
 
-`h_apply_cisd_diexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/CISD/H_apply.irp.f_shell_8#L1>`_
+h_apply_cisd_diexc
+  Undocumented
+
+
+h_apply_cisd_diexcorg
   Generate all double excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_cisd_monoexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/CISD/H_apply.irp.f_shell_8#L269>`_
+h_apply_cisd_diexcp
+  Undocumented
+
+
+h_apply_cisd_monoexc
   Generate all single excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.

plugins/CISD_SC2_selected/README.rst

@@ -72,3 +72,122 @@ Needed Modules
 
 * `CISD_selected <http://github.com/LCPQ/quantum_package/tree/master/src/CISD_selected>`_
 
+Needed Modules
+==============
+.. Do not edit this section It was auto-generated
+.. by the `update_README.py` script.
+
+
+.. image:: tree_dependency.png
+
+* `CISD_selected <http://github.com/LCPQ/quantum_package/tree/master/plugins/CISD_selected>`_
+
+Documentation
+=============
+.. Do not edit this section It was auto-generated
+.. by the `update_README.py` script.
+
+
+`cisd_sc2_selected <http://github.com/LCPQ/quantum_package/tree/master/plugins/CISD_SC2_selected/cisd_sc2_selection.irp.f#L1>`_
+  Undocumented
+
+
+h_apply_cisd
+  Calls H_apply on the HF determinant and selects all connected single and double
+  excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
+
+
+h_apply_cisd_diexc
+  Undocumented
+
+
+h_apply_cisd_diexcorg
+  Generate all double excitations of key_in using the bit masks of holes and
+  particles.
+  Assume N_int is already provided.
+
+
+h_apply_cisd_diexcp
+  Undocumented
+
+
+h_apply_cisd_monoexc
+  Generate all single excitations of key_in using the bit masks of holes and
+  particles.
+  Assume N_int is already provided.
+
+
+h_apply_pt2
+  Calls H_apply on the HF determinant and selects all connected single and double
+  excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
+
+
+h_apply_pt2_diexc
+  Undocumented
+
+
+h_apply_pt2_diexcorg
+  Generate all double excitations of key_in using the bit masks of holes and
+  particles.
+  Assume N_int is already provided.
+
+
+h_apply_pt2_diexcp
+  Undocumented
+
+
+h_apply_pt2_en_sc2
+  Calls H_apply on the HF determinant and selects all connected single and double
+  excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
+
+
+h_apply_pt2_en_sc2_diexc
+  Undocumented
+
+
+h_apply_pt2_en_sc2_diexcorg
+  Generate all double excitations of key_in using the bit masks of holes and
+  particles.
+  Assume N_int is already provided.
+
+
+h_apply_pt2_en_sc2_diexcp
+  Undocumented
+
+
+h_apply_pt2_en_sc2_monoexc
+  Generate all single excitations of key_in using the bit masks of holes and
+  particles.
+  Assume N_int is already provided.
+
+
+h_apply_pt2_monoexc
+  Generate all single excitations of key_in using the bit masks of holes and
+  particles.
+  Assume N_int is already provided.
+
+
+h_apply_sc2_selected
+  Calls H_apply on the HF determinant and selects all connected single and double
+  excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
+
+
+h_apply_sc2_selected_diexc
+  Undocumented
+
+
+h_apply_sc2_selected_diexcorg
+  Generate all double excitations of key_in using the bit masks of holes and
+  particles.
+  Assume N_int is already provided.
+
+
+h_apply_sc2_selected_diexcp
+  Undocumented
+
+
+h_apply_sc2_selected_monoexc
+  Generate all single excitations of key_in using the bit masks of holes and
+  particles.
+  Assume N_int is already provided.
+

plugins/CISD_selected/README.rst

@@ -196,22 +196,26 @@ Documentation
 .. by the `update_README.py` script.
 
 
-`cisd <http://github.com/LCPQ/quantum_package/tree/master/plugins/CISD_selected/cisd_selection.irp.f#L1>`_
-  Undocumented
-
-
-`h_apply_cisd <http://github.com/LCPQ/quantum_package/tree/master/plugins/CISD_selected/H_apply.irp.f_shell_8#L414>`_
+h_apply_cisd
   Calls H_apply on the HF determinant and selects all connected single and double
   excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
 
 
-`h_apply_cisd_diexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/CISD_selected/H_apply.irp.f_shell_8#L1>`_
+h_apply_cisd_diexc
+  Undocumented
+
+
+h_apply_cisd_diexcorg
   Generate all double excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_cisd_monoexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/CISD_selected/H_apply.irp.f_shell_8#L269>`_
+h_apply_cisd_diexcp
+  Undocumented
+
+
+h_apply_cisd_monoexc
   Generate all single excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
@@ -221,154 +225,226 @@ Documentation
   Undocumented
 
 
-`h_apply_cisd_selection_delta_rho_one_point <http://github.com/LCPQ/quantum_package/tree/master/plugins/CISD_selected/H_apply.irp.f_shell_10#L5931>`_
+h_apply_cisd_selection_delta_rho_one_point
   Calls H_apply on the HF determinant and selects all connected single and double
   excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
 
 
-`h_apply_cisd_selection_delta_rho_one_point_diexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/CISD_selected/H_apply.irp.f_shell_10#L5405>`_
+h_apply_cisd_selection_delta_rho_one_point_diexc
+  Undocumented
+
+
+h_apply_cisd_selection_delta_rho_one_point_diexcorg
   Generate all double excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_cisd_selection_delta_rho_one_point_monoexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/CISD_selected/H_apply.irp.f_shell_10#L5734>`_
+h_apply_cisd_selection_delta_rho_one_point_diexcp
+  Undocumented
+
+
+h_apply_cisd_selection_delta_rho_one_point_monoexc
   Generate all single excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_cisd_selection_dipole_moment_z <http://github.com/LCPQ/quantum_package/tree/master/plugins/CISD_selected/H_apply.irp.f_shell_10#L5159>`_
+h_apply_cisd_selection_dipole_moment_z
   Calls H_apply on the HF determinant and selects all connected single and double
   excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
 
 
-`h_apply_cisd_selection_dipole_moment_z_diexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/CISD_selected/H_apply.irp.f_shell_10#L4633>`_
+h_apply_cisd_selection_dipole_moment_z_diexc
+  Undocumented
+
+
+h_apply_cisd_selection_dipole_moment_z_diexcorg
   Generate all double excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_cisd_selection_dipole_moment_z_monoexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/CISD_selected/H_apply.irp.f_shell_10#L4962>`_
+h_apply_cisd_selection_dipole_moment_z_diexcp
+  Undocumented
+
+
+h_apply_cisd_selection_dipole_moment_z_monoexc
   Generate all single excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_cisd_selection_epstein_nesbet <http://github.com/LCPQ/quantum_package/tree/master/plugins/CISD_selected/H_apply.irp.f_shell_10#L3615>`_
+h_apply_cisd_selection_epstein_nesbet
   Calls H_apply on the HF determinant and selects all connected single and double
   excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
 
 
-`h_apply_cisd_selection_epstein_nesbet_2x2 <http://github.com/LCPQ/quantum_package/tree/master/plugins/CISD_selected/H_apply.irp.f_shell_10#L4387>`_
+h_apply_cisd_selection_epstein_nesbet_2x2
   Calls H_apply on the HF determinant and selects all connected single and double
   excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
 
 
-`h_apply_cisd_selection_epstein_nesbet_2x2_diexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/CISD_selected/H_apply.irp.f_shell_10#L3861>`_
+h_apply_cisd_selection_epstein_nesbet_2x2_diexc
+  Undocumented
+
+
+h_apply_cisd_selection_epstein_nesbet_2x2_diexcorg
   Generate all double excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_cisd_selection_epstein_nesbet_2x2_monoexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/CISD_selected/H_apply.irp.f_shell_10#L4190>`_
+h_apply_cisd_selection_epstein_nesbet_2x2_diexcp
+  Undocumented
+
+
+h_apply_cisd_selection_epstein_nesbet_2x2_monoexc
   Generate all single excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_cisd_selection_epstein_nesbet_diexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/CISD_selected/H_apply.irp.f_shell_10#L3089>`_
+h_apply_cisd_selection_epstein_nesbet_diexc
+  Undocumented
+
+
+h_apply_cisd_selection_epstein_nesbet_diexcorg
   Generate all double excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_cisd_selection_epstein_nesbet_monoexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/CISD_selected/H_apply.irp.f_shell_10#L3418>`_
+h_apply_cisd_selection_epstein_nesbet_diexcp
+  Undocumented
+
+
+h_apply_cisd_selection_epstein_nesbet_monoexc
   Generate all single excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_cisd_selection_epstein_nesbet_sc2 <http://github.com/LCPQ/quantum_package/tree/master/plugins/CISD_selected/H_apply.irp.f_shell_10#L2843>`_
+h_apply_cisd_selection_epstein_nesbet_sc2
   Calls H_apply on the HF determinant and selects all connected single and double
   excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
 
 
-`h_apply_cisd_selection_epstein_nesbet_sc2_diexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/CISD_selected/H_apply.irp.f_shell_10#L2317>`_
+h_apply_cisd_selection_epstein_nesbet_sc2_diexc
+  Undocumented
+
+
+h_apply_cisd_selection_epstein_nesbet_sc2_diexcorg
   Generate all double excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_cisd_selection_epstein_nesbet_sc2_monoexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/CISD_selected/H_apply.irp.f_shell_10#L2646>`_
+h_apply_cisd_selection_epstein_nesbet_sc2_diexcp
+  Undocumented
+
+
+h_apply_cisd_selection_epstein_nesbet_sc2_monoexc
   Generate all single excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_cisd_selection_epstein_nesbet_sc2_no_projected <http://github.com/LCPQ/quantum_package/tree/master/plugins/CISD_selected/H_apply.irp.f_shell_10#L2071>`_
+h_apply_cisd_selection_epstein_nesbet_sc2_no_projected
   Calls H_apply on the HF determinant and selects all connected single and double
   excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
 
 
-`h_apply_cisd_selection_epstein_nesbet_sc2_no_projected_diexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/CISD_selected/H_apply.irp.f_shell_10#L1545>`_
+h_apply_cisd_selection_epstein_nesbet_sc2_no_projected_diexc
+  Undocumented
+
+
+h_apply_cisd_selection_epstein_nesbet_sc2_no_projected_diexcorg
   Generate all double excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_cisd_selection_epstein_nesbet_sc2_no_projected_monoexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/CISD_selected/H_apply.irp.f_shell_10#L1874>`_
+h_apply_cisd_selection_epstein_nesbet_sc2_no_projected_diexcp
+  Undocumented
+
+
+h_apply_cisd_selection_epstein_nesbet_sc2_no_projected_monoexc
   Generate all single excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_cisd_selection_epstein_nesbet_sc2_projected <http://github.com/LCPQ/quantum_package/tree/master/plugins/CISD_selected/H_apply.irp.f_shell_10#L1299>`_
+h_apply_cisd_selection_epstein_nesbet_sc2_projected
   Calls H_apply on the HF determinant and selects all connected single and double
   excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
 
 
-`h_apply_cisd_selection_epstein_nesbet_sc2_projected_diexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/CISD_selected/H_apply.irp.f_shell_10#L773>`_
+h_apply_cisd_selection_epstein_nesbet_sc2_projected_diexc
+  Undocumented
+
+
+h_apply_cisd_selection_epstein_nesbet_sc2_projected_diexcorg
   Generate all double excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_cisd_selection_epstein_nesbet_sc2_projected_monoexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/CISD_selected/H_apply.irp.f_shell_10#L1102>`_
+h_apply_cisd_selection_epstein_nesbet_sc2_projected_diexcp
+  Undocumented
+
+
+h_apply_cisd_selection_epstein_nesbet_sc2_projected_monoexc
   Generate all single excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_cisd_selection_h_core <http://github.com/LCPQ/quantum_package/tree/master/plugins/CISD_selected/H_apply.irp.f_shell_10#L527>`_
+h_apply_cisd_selection_h_core
   Calls H_apply on the HF determinant and selects all connected single and double
   excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
 
 
-`h_apply_cisd_selection_h_core_diexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/CISD_selected/H_apply.irp.f_shell_10#L1>`_
+h_apply_cisd_selection_h_core_diexc
+  Undocumented
+
+
+h_apply_cisd_selection_h_core_diexcorg
   Generate all double excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_cisd_selection_h_core_monoexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/CISD_selected/H_apply.irp.f_shell_10#L330>`_
+h_apply_cisd_selection_h_core_diexcp
+  Undocumented
+
+
+h_apply_cisd_selection_h_core_monoexc
   Generate all single excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_cisd_selection_moller_plesset <http://github.com/LCPQ/quantum_package/tree/master/plugins/CISD_selected/H_apply.irp.f_shell_10#L6703>`_
+h_apply_cisd_selection_moller_plesset
   Calls H_apply on the HF determinant and selects all connected single and double
   excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
 
 
-`h_apply_cisd_selection_moller_plesset_diexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/CISD_selected/H_apply.irp.f_shell_10#L6177>`_
+h_apply_cisd_selection_moller_plesset_diexc
+  Undocumented
+
+
+h_apply_cisd_selection_moller_plesset_diexcorg
   Generate all double excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_cisd_selection_moller_plesset_monoexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/CISD_selected/H_apply.irp.f_shell_10#L6506>`_
+h_apply_cisd_selection_moller_plesset_diexcp
+  Undocumented
+
+
+h_apply_cisd_selection_moller_plesset_monoexc
   Generate all single excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.

plugins/CISD_selected/cisd_selection.irp.f

@@ -41,8 +41,8 @@ program cisd
   N_det = min(N_det,N_det_max)
   touch N_det psi_det psi_coef
   call diagonalize_CI
-  deallocate(pt2,norm_pert,H_pert_diag)
   call save_wavefunction
   call ezfio_set_cisd_selected_energy(CI_energy)
   call ezfio_set_cisd_selected_energy_pt2(CI_energy+pt2)
+  deallocate(pt2,norm_pert,H_pert_diag)
 end

plugins/DDCI_selected/README.rst

@@ -57,3 +57,74 @@ Needed Modules
 * `Selectors_full <http://github.com/LCPQ/quantum_package/tree/master/src/Selectors_full>`_
 * `Generators_CAS <http://github.com/LCPQ/quantum_package/tree/master/src/Generators_CAS>`_
 
+Needed Modules
+==============
+.. Do not edit this section It was auto-generated
+.. by the `update_README.py` script.
+
+
+.. image:: tree_dependency.png
+
+* `Perturbation <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation>`_
+* `Selectors_full <http://github.com/LCPQ/quantum_package/tree/master/plugins/Selectors_full>`_
+* `Generators_CAS <http://github.com/LCPQ/quantum_package/tree/master/plugins/Generators_CAS>`_
+
+Documentation
+=============
+.. Do not edit this section It was auto-generated
+.. by the `update_README.py` script.
+
+
+`ddci <http://github.com/LCPQ/quantum_package/tree/master/plugins/DDCI_selected/ddci.irp.f#L1>`_
+  Undocumented
+
+
+h_apply_ddci_pt2
+  Calls H_apply on the HF determinant and selects all connected single and double
+  excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
+
+
+h_apply_ddci_pt2_diexc
+  Undocumented
+
+
+h_apply_ddci_pt2_diexcorg
+  Generate all double excitations of key_in using the bit masks of holes and
+  particles.
+  Assume N_int is already provided.
+
+
+h_apply_ddci_pt2_diexcp
+  Undocumented
+
+
+h_apply_ddci_pt2_monoexc
+  Generate all single excitations of key_in using the bit masks of holes and
+  particles.
+  Assume N_int is already provided.
+
+
+h_apply_ddci_selection
+  Calls H_apply on the HF determinant and selects all connected single and double
+  excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
+
+
+h_apply_ddci_selection_diexc
+  Undocumented
+
+
+h_apply_ddci_selection_diexcorg
+  Generate all double excitations of key_in using the bit masks of holes and
+  particles.
+  Assume N_int is already provided.
+
+
+h_apply_ddci_selection_diexcp
+  Undocumented
+
+
+h_apply_ddci_selection_monoexc
+  Generate all single excitations of key_in using the bit masks of holes and
+  particles.
+  Assume N_int is already provided.
+

plugins/Full_CI/.gitignore

@@ -0,0 +1,32 @@
+# Automatically created by $QP_ROOT/scripts/module/module_handler.py 
+.ninja_deps
+.ninja_log
+AO_Basis
+Bitmask
+Determinants
+Electrons
+Ezfio_files
+Generators_full
+Hartree_Fock
+IRPF90_man
+IRPF90_temp
+Integrals_Bielec
+Integrals_Monoelec
+MOGuess
+MO_Basis
+Makefile
+Makefile.depend
+Nuclei
+Perturbation
+Properties
+Pseudo
+Selectors_full
+Utils
+ezfio_interface.irp.f
+full_ci
+full_ci_no_skip
+irpf90.make
+irpf90_entities
+tags
+target_pt2
+var_pt2_ratio

plugins/Full_CI/README.rst

@@ -27,137 +27,201 @@ Documentation
   Undocumented
 
 
-`h_apply_fci <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L527>`_
+h_apply_fci
   Calls H_apply on the HF determinant and selects all connected single and double
   excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
 
 
-`h_apply_fci_diexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L1>`_
+h_apply_fci_diexc
+  Undocumented
+
+
+h_apply_fci_diexcorg
   Generate all double excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_fci_mono <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L2744>`_
+h_apply_fci_diexcp
+  Undocumented
+
+
+h_apply_fci_mono
   Calls H_apply on the HF determinant and selects all connected single and double
   excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
 
 
-`h_apply_fci_mono_diexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L2216>`_
+h_apply_fci_mono_diexc
+  Undocumented
+
+
+h_apply_fci_mono_diexcorg
   Generate all double excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_fci_mono_monoexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L2545>`_
+h_apply_fci_mono_diexcp
+  Undocumented
+
+
+h_apply_fci_mono_monoexc
   Generate all single excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_fci_monoexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L330>`_
+h_apply_fci_monoexc
   Generate all single excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_fci_no_skip <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L1998>`_
+h_apply_fci_no_skip
   Calls H_apply on the HF determinant and selects all connected single and double
   excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
 
 
-`h_apply_fci_no_skip_diexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L1472>`_
+h_apply_fci_no_skip_diexc
+  Undocumented
+
+
+h_apply_fci_no_skip_diexcorg
   Generate all double excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_fci_no_skip_monoexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L1801>`_
+h_apply_fci_no_skip_diexcp
+  Undocumented
+
+
+h_apply_fci_no_skip_monoexc
   Generate all single excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_fci_pt2 <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L1265>`_
+h_apply_fci_pt2
   Calls H_apply on the HF determinant and selects all connected single and double
   excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
 
 
-`h_apply_fci_pt2_diexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L773>`_
+h_apply_fci_pt2_diexc
+  Undocumented
+
+
+h_apply_fci_pt2_diexcorg
   Generate all double excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_fci_pt2_monoexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L1082>`_
+h_apply_fci_pt2_diexcp
+  Undocumented
+
+
+h_apply_fci_pt2_monoexc
   Generate all single excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_pt2_mono_delta_rho <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L4258>`_
+h_apply_pt2_mono_delta_rho
   Calls H_apply on the HF determinant and selects all connected single and double
   excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
 
 
-`h_apply_pt2_mono_delta_rho_diexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L3764>`_
+h_apply_pt2_mono_delta_rho_diexc
+  Undocumented
+
+
+h_apply_pt2_mono_delta_rho_diexcorg
   Generate all double excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_pt2_mono_delta_rho_monoexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L4073>`_
+h_apply_pt2_mono_delta_rho_diexcp
+  Undocumented
+
+
+h_apply_pt2_mono_delta_rho_monoexc
   Generate all single excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_pt2_mono_di_delta_rho <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L5729>`_
+h_apply_pt2_mono_di_delta_rho
   Calls H_apply on the HF determinant and selects all connected single and double
   excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
 
 
-`h_apply_pt2_mono_di_delta_rho_diexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L5237>`_
+h_apply_pt2_mono_di_delta_rho_diexc
+  Undocumented
+
+
+h_apply_pt2_mono_di_delta_rho_diexcorg
   Generate all double excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_pt2_mono_di_delta_rho_monoexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L5546>`_
+h_apply_pt2_mono_di_delta_rho_diexcp
+  Undocumented
+
+
+h_apply_pt2_mono_di_delta_rho_monoexc
   Generate all single excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_select_mono_delta_rho <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L3518>`_
+h_apply_select_mono_delta_rho
   Calls H_apply on the HF determinant and selects all connected single and double
   excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
 
 
-`h_apply_select_mono_delta_rho_diexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L2990>`_
+h_apply_select_mono_delta_rho_diexc
+  Undocumented
+
+
+h_apply_select_mono_delta_rho_diexcorg
   Generate all double excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_select_mono_delta_rho_monoexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L3319>`_
+h_apply_select_mono_delta_rho_diexcp
+  Undocumented
+
+
+h_apply_select_mono_delta_rho_monoexc
   Generate all single excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_select_mono_di_delta_rho <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L4991>`_
+h_apply_select_mono_di_delta_rho
   Calls H_apply on the HF determinant and selects all connected single and double
   excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
 
 
-`h_apply_select_mono_di_delta_rho_diexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L4465>`_
+h_apply_select_mono_di_delta_rho_diexc
+  Undocumented
+
+
+h_apply_select_mono_di_delta_rho_diexcorg
   Generate all double excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_select_mono_di_delta_rho_monoexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L4794>`_
+h_apply_select_mono_di_delta_rho_diexcp
+  Undocumented
+
+
+h_apply_select_mono_di_delta_rho_monoexc
   Generate all single excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.

plugins/Full_CI/target_pt2.irp.f

@@ -27,6 +27,8 @@ program var_pt2_ratio_run
     call diagonalize_CI
     ratio = (CI_energy(1) - HF_energy) / (CI_energy(1)+pt2(1) - HF_energy)
     if (N_det > 20000) then
+      N_det = 20000
+      TOUCH N_det
       exit
     endif
   enddo

plugins/Hartree_Fock/EZFIO.cfg

@@ -10,6 +10,12 @@ doc: Maximum number of SCF iterations
 interface: ezfio,provider,ocaml
 default: 200
 
+[level_shift]
+type: Positive_float
+doc: Energy shift on the virtual MOs to improve SCF convergence
+interface: ezfio,provider,ocaml
+default: 0.5
+
 [mo_guess_type]
 type: MO_guess
 doc: Initial MO guess. Can be [ Huckel | HCore ]

plugins/Hartree_Fock/Fock_matrix.irp.f

@@ -73,6 +73,11 @@
      enddo
      
    endif
+   ! Introduce level shift here
+   do i = elec_alpha_num+1, mo_tot_num
+     Fock_matrix_mo(i,i) += level_shift
+   enddo
+
    do i = 1, mo_tot_num
      Fock_matrix_diag_mo(i) = Fock_matrix_mo(i,i) 
    enddo
@@ -108,9 +113,10 @@ END_PROVIDER
  END_DOC
  
  integer                        :: i,j,k,l,k1,r,s
+ integer                        :: i0,j0,k0,l0
  integer*8                      :: p,q
- double precision               :: integral
- double precision               :: ao_bielec_integral
+ double precision               :: integral, c0, c1, c2
+ double precision               :: ao_bielec_integral, local_threshold
  double precision, allocatable  :: ao_bi_elec_integral_alpha_tmp(:,:)
  double precision, allocatable  :: ao_bi_elec_integral_beta_tmp(:,:)
  !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: ao_bi_elec_integral_beta_tmp
@@ -121,8 +127,9 @@ END_PROVIDER
  if (do_direct_integrals) then
 
    !$OMP PARALLEL DEFAULT(NONE)                                      &
-       !$OMP PRIVATE(i,j,l,k1,k,integral,ii,jj,kk,ll,i8,keys,values,p,q,r,s, &
-       !$OMP ao_bi_elec_integral_alpha_tmp,ao_bi_elec_integral_beta_tmp)&
+       !$OMP PRIVATE(i,j,l,k1,k,integral,ii,jj,kk,ll,i8,keys,values,p,q,r,s,i0,j0,k0,l0, &
+       !$OMP ao_bi_elec_integral_alpha_tmp,ao_bi_elec_integral_beta_tmp, c0, c1, c2, &
+       !$OMP local_threshold)&
        !$OMP SHARED(ao_num,ao_num_align,HF_density_matrix_ao_alpha,HF_density_matrix_ao_beta,&
        !$OMP ao_integrals_map,ao_integrals_threshold, ao_bielec_integral_schwartz, &
        !$OMP ao_overlap_abs, ao_bi_elec_integral_alpha, ao_bi_elec_integral_beta)
@@ -152,14 +159,16 @@ END_PROVIDER
               < ao_integrals_threshold) then
              cycle
            endif
-           if (ao_bielec_integral_schwartz(k,l)*ao_bielec_integral_schwartz(i,j)  &
-              < ao_integrals_threshold) then
-             cycle
-           endif
-           values(1) = ao_bielec_integral(k,l,i,j)
-           if (abs(values(1)) < ao_integrals_threshold) then
+           local_threshold = ao_bielec_integral_schwartz(k,l)*ao_bielec_integral_schwartz(i,j)
+           if (local_threshold < ao_integrals_threshold) then
              cycle
            endif
+           i0 = i
+           j0 = j
+           k0 = k
+           l0 = l
+           values(1) = 0.d0
+           local_threshold = ao_integrals_threshold/local_threshold
            do k2=1,8
              if (kk(k2)==0) then
                cycle
@@ -168,12 +177,21 @@ END_PROVIDER
              j = jj(k2)
              k = kk(k2)
              l = ll(k2)
-             integral = (HF_density_matrix_ao_alpha(k,l)+HF_density_matrix_ao_beta(k,l)) * values(1)
+             c0 = HF_density_matrix_ao_alpha(k,l)+HF_density_matrix_ao_beta(k,l)
+             c1 = HF_density_matrix_ao_alpha(k,i)
+             c2 = HF_density_matrix_ao_beta(k,i)
+             if ( dabs(c0)+dabs(c1)+dabs(c2) < local_threshold) then
+               cycle
+             endif
+             if (values(1) == 0.d0) then
+               values(1) = ao_bielec_integral(k0,l0,i0,j0)
+             endif
+             integral = c0 * values(1)
              ao_bi_elec_integral_alpha_tmp(i,j) += integral
              ao_bi_elec_integral_beta_tmp (i,j) += integral
              integral = values(1)
-             ao_bi_elec_integral_alpha_tmp(l,j) -= HF_density_matrix_ao_alpha(k,i) * integral
-             ao_bi_elec_integral_beta_tmp (l,j) -= HF_density_matrix_ao_beta (k,i) * integral
+             ao_bi_elec_integral_alpha_tmp(l,j) -= c1 * integral
+             ao_bi_elec_integral_beta_tmp (l,j) -= c2 * integral
            enddo
    enddo
    !$OMP END DO NOWAIT
@@ -315,7 +333,9 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_ao, (ao_num_align, ao_num) ]
  ! Fock matrix in AO basis set
  END_DOC
  
- if (elec_alpha_num == elec_beta_num) then
+ if ( (elec_alpha_num == elec_beta_num).and. &
+      (level_shift == 0.) ) &
+ then
    integer                        :: i,j
    do j=1,ao_num
      !DIR$ VECTOR ALIGNED
@@ -325,29 +345,46 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_ao, (ao_num_align, ao_num) ]
    enddo
  else
    double precision, allocatable  :: T(:,:), M(:,:)
+   integer                        :: ierr
    ! F_ao = S C F_mo C^t S
-   allocate (T(ao_num_align,ao_num),M(ao_num_align,ao_num))
-   call dgemm('N','N', ao_num,ao_num,ao_num, 1.d0,                   &
+   allocate (T(ao_num_align,ao_num),M(ao_num_align,ao_num),stat=ierr)
+   if (ierr /=0 ) then
+      print *,  irp_here, ' : allocation failed'
+   endif
+
+!  ao_overlap (ao_num,ao_num) . mo_coef (ao_num,mo_tot_num)
+!  -> M(ao_num,mo_tot_num)
+   call dgemm('N','N', ao_num,mo_tot_num,ao_num, 1.d0,               &
        ao_overlap, size(ao_overlap,1),                               &
        mo_coef, size(mo_coef,1),                                     &
        0.d0,                                                         &
        M, size(M,1))
+
+!  M(ao_num,mo_tot_num) . Fock_matrix_mo (mo_tot_num,mo_tot_num)
+!  -> T(ao_num,mo_tot_num)
    call dgemm('N','N', ao_num,mo_tot_num,mo_tot_num, 1.d0,           &
        M, size(M,1),                                                 &
        Fock_matrix_mo, size(Fock_matrix_mo,1),                       &
        0.d0,                                                         &
        T, size(T,1))
-   call dgemm('N','T', mo_tot_num,ao_num,mo_tot_num, 1.d0,           &
+
+!  T(ao_num,mo_tot_num) . mo_coef^T (mo_tot_num,ao_num)
+!  -> M(ao_num,ao_num)
+   call dgemm('N','T', ao_num,ao_num,mo_tot_num, 1.d0,           &
        T, size(T,1),                                                 &
        mo_coef, size(mo_coef,1),                                     &
        0.d0,                                                         &
        M, size(M,1))
+
+!  M(ao_num,ao_num) . ao_overlap (ao_num,ao_num)
+!  -> Fock_matrix_ao(ao_num,ao_num)
    call dgemm('N','N', ao_num,ao_num,ao_num, 1.d0,                   &
        M, size(M,1),                                                 &
        ao_overlap, size(ao_overlap,1),                               &
        0.d0,                                                         &
        Fock_matrix_ao, size(Fock_matrix_ao,1))
 
+
    deallocate(T)
  endif
 END_PROVIDER
@@ -360,23 +397,39 @@ subroutine Fock_mo_to_ao(FMO,LDFMO,FAO,LDFAO)
   double precision, intent(out)  :: FAO(LDFAO,*)
   
   double precision, allocatable  :: T(:,:), M(:,:)
+  integer                        :: ierr
   ! F_ao = S C F_mo C^t S
-  allocate (T(ao_num_align,ao_num),M(ao_num_align,ao_num))
-  call dgemm('N','N', ao_num,ao_num,ao_num, 1.d0,                    &
+  allocate (T(ao_num_align,ao_num),M(ao_num_align,ao_num),stat=ierr)
+  if (ierr /=0 ) then
+     print *,  irp_here, ' : allocation failed'
+  endif
+
+!  ao_overlap (ao_num,ao_num) . mo_coef (ao_num,mo_tot_num)
+!  -> M(ao_num,mo_tot_num)
+  call dgemm('N','N', ao_num,mo_tot_num,ao_num, 1.d0,                    &
       ao_overlap, size(ao_overlap,1),                                &
       mo_coef, size(mo_coef,1),                                      &
       0.d0,                                                          &
       M, size(M,1))
+
+!  M(ao_num,mo_tot_num) . FMO (mo_tot_num,mo_tot_num)
+!  -> T(ao_num,mo_tot_num)
   call dgemm('N','N', ao_num,mo_tot_num,mo_tot_num, 1.d0,            &
       M, size(M,1),                                                  &
       FMO, size(FMO,1),                                              &
       0.d0,                                                          &
       T, size(T,1))
-  call dgemm('N','T', mo_tot_num,ao_num,mo_tot_num, 1.d0,            &
+
+!  T(ao_num,mo_tot_num) . mo_coef^T (mo_tot_num,ao_num)
+!  -> M(ao_num,ao_num)
+  call dgemm('N','T', ao_num,ao_num,mo_tot_num, 1.d0,            &
       T, size(T,1),                                                  &
       mo_coef, size(mo_coef,1),                                      &
       0.d0,                                                          &
       M, size(M,1))
+
+!  M(ao_num,ao_num) . ao_overlap (ao_num,ao_num)
+!  -> Fock_matrix_ao(ao_num,ao_num)
   call dgemm('N','N', ao_num,ao_num,ao_num, 1.d0,                    &
       M, size(M,1),                                                  &
       ao_overlap, size(ao_overlap,1),                                &

plugins/Hartree_Fock/README.rst

@@ -32,11 +32,11 @@ Documentation
 .. by the `update_README.py` script.
 
 
-`ao_bi_elec_integral_alpha <http://github.com/LCPQ/quantum_package/tree/master/plugins/Hartree_Fock/Fock_matrix.irp.f#L102>`_
+`ao_bi_elec_integral_alpha <http://github.com/LCPQ/quantum_package/tree/master/plugins/Hartree_Fock/Fock_matrix.irp.f#L107>`_
   Alpha Fock matrix in AO basis set
 
 
-`ao_bi_elec_integral_beta <http://github.com/LCPQ/quantum_package/tree/master/plugins/Hartree_Fock/Fock_matrix.irp.f#L103>`_
+`ao_bi_elec_integral_beta <http://github.com/LCPQ/quantum_package/tree/master/plugins/Hartree_Fock/Fock_matrix.irp.f#L108>`_
   Alpha Fock matrix in AO basis set
 
 
@@ -62,23 +62,23 @@ Documentation
   Diagonal Fock matrix in the MO basis
 
 
-`fock_matrix_alpha_ao <http://github.com/LCPQ/quantum_package/tree/master/plugins/Hartree_Fock/Fock_matrix.irp.f#L83>`_
+`fock_matrix_alpha_ao <http://github.com/LCPQ/quantum_package/tree/master/plugins/Hartree_Fock/Fock_matrix.irp.f#L88>`_
   Alpha Fock matrix in AO basis set
 
 
-`fock_matrix_alpha_mo <http://github.com/LCPQ/quantum_package/tree/master/plugins/Hartree_Fock/Fock_matrix.irp.f#L254>`_
+`fock_matrix_alpha_mo <http://github.com/LCPQ/quantum_package/tree/master/plugins/Hartree_Fock/Fock_matrix.irp.f#L272>`_
   Fock matrix on the MO basis
 
 
-`fock_matrix_ao <http://github.com/LCPQ/quantum_package/tree/master/plugins/Hartree_Fock/Fock_matrix.irp.f#L312>`_
+`fock_matrix_ao <http://github.com/LCPQ/quantum_package/tree/master/plugins/Hartree_Fock/Fock_matrix.irp.f#L330>`_
   Fock matrix in AO basis set
 
 
-`fock_matrix_beta_ao <http://github.com/LCPQ/quantum_package/tree/master/plugins/Hartree_Fock/Fock_matrix.irp.f#L84>`_
+`fock_matrix_beta_ao <http://github.com/LCPQ/quantum_package/tree/master/plugins/Hartree_Fock/Fock_matrix.irp.f#L89>`_
   Alpha Fock matrix in AO basis set
 
 
-`fock_matrix_beta_mo <http://github.com/LCPQ/quantum_package/tree/master/plugins/Hartree_Fock/Fock_matrix.irp.f#L274>`_
+`fock_matrix_beta_mo <http://github.com/LCPQ/quantum_package/tree/master/plugins/Hartree_Fock/Fock_matrix.irp.f#L292>`_
   Fock matrix on the MO basis
 
 
@@ -114,7 +114,7 @@ Documentation
   .br
 
 
-`fock_mo_to_ao <http://github.com/LCPQ/quantum_package/tree/master/plugins/Hartree_Fock/Fock_matrix.irp.f#L355>`_
+`fock_mo_to_ao <http://github.com/LCPQ/quantum_package/tree/master/plugins/Hartree_Fock/Fock_matrix.irp.f#L375>`_
   Undocumented
 
 
@@ -134,7 +134,7 @@ Documentation
   S^-1 Beta density matrix in the AO basis x S^-1
 
 
-`hf_energy <http://github.com/LCPQ/quantum_package/tree/master/plugins/Hartree_Fock/Fock_matrix.irp.f#L293>`_
+`hf_energy <http://github.com/LCPQ/quantum_package/tree/master/plugins/Hartree_Fock/Fock_matrix.irp.f#L311>`_
   Hartree-Fock energy
 
 
@@ -142,7 +142,11 @@ Documentation
   Build the MOs using the extended Huckel model
 
 
-`mo_guess_type <http://github.com/LCPQ/quantum_package/tree/master/plugins/Hartree_Fock/ezfio_interface.irp.f#L28>`_
+`level_shift <http://github.com/LCPQ/quantum_package/tree/master/plugins/Hartree_Fock/ezfio_interface.irp.f#L28>`_
+  Energy shift on the virtual MOs to improve SCF convergence
+
+
+`mo_guess_type <http://github.com/LCPQ/quantum_package/tree/master/plugins/Hartree_Fock/ezfio_interface.irp.f#L50>`_
   Initial MO guess. Can be [ Huckel | HCore ]
 
 
@@ -161,6 +165,6 @@ Documentation
   optional: mo_basis.mo_coef
 
 
-`thresh_scf <http://github.com/LCPQ/quantum_package/tree/master/plugins/Hartree_Fock/ezfio_interface.irp.f#L46>`_
+`thresh_scf <http://github.com/LCPQ/quantum_package/tree/master/plugins/Hartree_Fock/ezfio_interface.irp.f#L68>`_
   Threshold on the convergence of the Hartree Fock energy
 

plugins/Hartree_Fock/SCF.irp.f

@@ -42,16 +42,13 @@ subroutine run
   BEGIN_DOC
 ! Run SCF calculation
   END_DOC
-  double precision               :: SCF_energy_before,SCF_energy_after,diag_H_mat_elem,get_mo_bielec_integral
+  double precision               :: SCF_energy_before,SCF_energy_after,diag_H_mat_elem
   double precision               :: E0
   integer                        :: i_it, i, j, k
    
   E0 = HF_energy 
 
-  thresh_SCF = 1.d-10
-  call damping_SCF
   mo_label = "Canonical"
-  TOUCH mo_label mo_coef
-  call save_mos
+  call damping_SCF
   
 end

plugins/Hartree_Fock/damping_SCF.irp.f

@@ -86,7 +86,7 @@ subroutine damping_SCF
       if ((E_half > E).and.(E_new < E)) then
         lambda = 1.d0
         exit
-      else if ((E_half > E).and.(lambda > 5.d-2)) then
+      else if ((E_half > E).and.(lambda > 5.d-4)) then
         lambda = 0.5d0 * lambda
         E_new = E_half
       else
@@ -119,7 +119,7 @@ subroutine damping_SCF
   write(output_hartree_fock,'(A4,X,A16, X, A16, X, A16, X, A4 )'), '====','================','================','================', '===='
   write(output_hartree_fock,*)
   
-  call mo_as_eigvectors_of_mo_matrix(Fock_matrix_mo,size(Fock_matrix_mo,1),size(Fock_matrix_mo,2),mo_label)
+  call mo_as_eigvectors_of_mo_matrix(Fock_matrix_mo,size(Fock_matrix_mo,1),size(Fock_matrix_mo,2),mo_label,1)
 
   call write_double(output_hartree_fock, E_min, 'Hartree-Fock energy')
   call ezfio_set_hartree_fock_energy(E_min)

plugins/Hartree_Fock/diagonalize_fock.irp.f

@@ -10,6 +10,10 @@
    integer, allocatable           :: iwork(:)
    double precision, allocatable  :: work(:), F(:,:), S(:,:)
    
+
+   if (mo_tot_num == ao_num) then
+       ! Solve H.C = E.S.C in AO basis set
+
        allocate(F(ao_num_align,ao_num), S(ao_num_align,ao_num) )
        do j=1,ao_num
          do i=1,ao_num
@@ -29,10 +33,6 @@
 
        call dsygvd(1,'v','u',ao_num,F,size(F,1),S,size(S,1),&
          diagonal_Fock_matrix_mo, work, lwork, iwork, liwork, info)
-!    call dsygv(1, 'v', 'u',ao_num,F,size(F,1),S,size(S,1),&
-!     diagonal_Fock_matrix_mo, work, lwork, info)
-
-
 
        if (info /= 0) then
          print *,  irp_here//' failed : ', info
@@ -42,15 +42,10 @@
        liwork = iwork(1)
        deallocate(work,iwork)
        allocate(work(lwork), iwork(liwork) )
-!   deallocate(work)
-!   allocate(work(lwork))
 
        call dsygvd(1,'v','u',ao_num,F,size(F,1),S,size(S,1),&
          diagonal_Fock_matrix_mo, work, lwork, iwork, liwork, info)
 
-!    call dsygv(1, 'v', 'u',ao_num,F,size(F,1),S,size(S,1),&
-!     diagonal_Fock_matrix_mo, work, lwork, info)
-
        if (info /= 0) then
          print *,  irp_here//' failed : ', info
          stop 1
@@ -62,6 +57,56 @@
        enddo
 
        deallocate(work, iwork, F, S)
+
+  else 
+
+       ! Solve H.C = E.C in MO basis set
+       
+       allocate( F(mo_tot_num_align,mo_tot_num) )
+       do j=1,mo_tot_num
+         do i=1,mo_tot_num
+           F(i,j) = Fock_matrix_mo(i,j)
+         enddo
+       enddo
+       
+       n = mo_tot_num
+       lwork = 1+6*n + 2*n*n
+       liwork = 3 + 5*n
+       
+       allocate(work(lwork), iwork(liwork) )
+       
+       lwork = -1
+       liwork = -1
+       
+       call dsyevd( 'V', 'U', mo_tot_num, F,                         &
+           size(F,1), diagonal_Fock_matrix_mo,                       &
+           work, lwork, iwork, liwork, info)
+       
+       if (info /= 0) then
+         print *,  irp_here//' failed : ', info
+         stop 1
+       endif
+       lwork = int(work(1))
+       liwork = iwork(1)
+       deallocate(work,iwork)
+       allocate(work(lwork), iwork(liwork) )
+       
+       call dsyevd( 'V', 'U', mo_tot_num, F,                         &
+           size(F,1), diagonal_Fock_matrix_mo,                       &
+           work, lwork, iwork, liwork, info)
+       
+       if (info /= 0) then
+         print *,  irp_here//' failed : ', info
+         stop 1
+       endif
+
+       call dgemm('N','N',ao_num,mo_tot_num,mo_tot_num, 1.d0, &
+         mo_coef, size(mo_coef,1), F, size(F,1), &
+         0.d0, eigenvectors_Fock_matrix_mo, size(eigenvectors_Fock_matrix_mo,1))
+       deallocate(work, iwork, F)
+
+  endif
+
 END_PROVIDER
  
 BEGIN_PROVIDER [double precision, diagonal_Fock_matrix_mo_sum, (mo_tot_num)]
@@ -73,19 +118,19 @@ BEGIN_PROVIDER [double precision, diagonal_Fock_matrix_mo_sum, (mo_tot_num)]
  END_DOC
  integer :: i,j
  double precision :: accu
+ do j = 1,elec_alpha_num
+  accu = 0.d0
   do i = 1, elec_alpha_num
-  accu = 0.d0
-  do j = 1, elec_alpha_num
    accu += 2.d0 * mo_bielec_integral_jj_from_ao(i,j) - mo_bielec_integral_jj_exchange_from_ao(i,j)
   enddo
-  diagonal_Fock_matrix_mo_sum(i) = accu + mo_mono_elec_integral(i,i)
+  diagonal_Fock_matrix_mo_sum(j) = accu + mo_mono_elec_integral(j,j)
  enddo
- do i = elec_alpha_num+1,mo_tot_num
+ do j = elec_alpha_num+1,mo_tot_num
   accu = 0.d0
-  do j = 1, elec_alpha_num
+  do i = 1, elec_alpha_num
    accu += 2.d0 * mo_bielec_integral_jj_from_ao(i,j) - mo_bielec_integral_jj_exchange_from_ao(i,j)
   enddo
-  diagonal_Fock_matrix_mo_sum(i) = accu + mo_mono_elec_integral(i,i)
+  diagonal_Fock_matrix_mo_sum(j) = accu + mo_mono_elec_integral(j,j)
  enddo
 
 END_PROVIDER

plugins/Hartree_Fock/huckel.irp.f

@@ -4,7 +4,7 @@ subroutine huckel_guess
 ! Build the MOs using the extended Huckel model
   END_DOC
   integer                        :: i,j
-  double precision               :: tmp_matrix(ao_num_align,ao_num),accu
+  double precision               :: accu
   double precision               :: c
   character*(64)                 :: label
 
@@ -13,20 +13,18 @@ subroutine huckel_guess
   label = "Guess"
   call mo_as_eigvectors_of_mo_matrix(mo_mono_elec_integral,          &
                                      size(mo_mono_elec_integral,1),  &
-                                     size(mo_mono_elec_integral,2),label)
+                                     size(mo_mono_elec_integral,2),label,1)
   TOUCH mo_coef
 
   c = 0.5d0 * 1.75d0
 
   do j=1,ao_num
+    !DIR$ VECTOR ALIGNED
     do i=1,ao_num
-      if (i.ne.j) then
-        Fock_matrix_ao(i,j) = c*ao_overlap(i,j)*(ao_mono_elec_integral(i,i) + &
-                                                 ao_mono_elec_integral(j,j))
-      else
-        Fock_matrix_ao(i,j) = Fock_matrix_alpha_ao(i,j)
-      endif
+      Fock_matrix_ao(i,j) = c*ao_overlap(i,j)*(ao_mono_elec_integral_diag(i) + &
+                                                 ao_mono_elec_integral_diag(j))
     enddo
+    Fock_matrix_ao(j,j) = Fock_matrix_alpha_ao(j,j)
   enddo
   TOUCH Fock_matrix_ao
   mo_coef = eigenvectors_fock_matrix_mo

plugins/MP2/EZFIO.cfg

@@ -0,0 +1,5 @@
+[energy]
+type: double precision
+doc: MP2 energy
+interface: ezfio
+

plugins/MP2/H_apply.irp.f

@@ -6,5 +6,9 @@ from perturbation import perturbations
 s = H_apply("mp2")
 s.set_perturbation("Moller_plesset")
 print s
+
+s = H_apply("mp2_selection")
+s.set_selection_pt2("Moller_plesset")
+print s
 END_SHELL
 

plugins/MP2/README.rst

@@ -40,3 +40,49 @@ Needed Modules
 * `Selectors_full <http://github.com/LCPQ/quantum_package/tree/master/src/Selectors_full>`_
 * `SingleRefMethod <http://github.com/LCPQ/quantum_package/tree/master/src/SingleRefMethod>`_
 
+Needed Modules
+==============
+.. Do not edit this section It was auto-generated
+.. by the `update_README.py` script.
+
+
+.. image:: tree_dependency.png
+
+* `Perturbation <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation>`_
+* `Selectors_full <http://github.com/LCPQ/quantum_package/tree/master/plugins/Selectors_full>`_
+* `SingleRefMethod <http://github.com/LCPQ/quantum_package/tree/master/plugins/SingleRefMethod>`_
+
+Documentation
+=============
+.. Do not edit this section It was auto-generated
+.. by the `update_README.py` script.
+
+
+h_apply_mp2
+  Calls H_apply on the HF determinant and selects all connected single and double
+  excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
+
+
+h_apply_mp2_diexc
+  Undocumented
+
+
+h_apply_mp2_diexcorg
+  Generate all double excitations of key_in using the bit masks of holes and
+  particles.
+  Assume N_int is already provided.
+
+
+h_apply_mp2_diexcp
+  Undocumented
+
+
+h_apply_mp2_monoexc
+  Generate all single excitations of key_in using the bit masks of holes and
+  particles.
+  Assume N_int is already provided.
+
+
+`mp2 <http://github.com/LCPQ/quantum_package/tree/master/plugins/MP2/mp2.irp.f#L1>`_
+  Undocumented
+

plugins/MP2/mp2_wf.irp.f

@@ -0,0 +1,31 @@
+program mp2_wf
+  implicit none
+  BEGIN_DOC
+! Save the MP2 wave function
+  END_DOC
+  integer                        :: i,k
+
+  
+  double precision, allocatable  :: pt2(:), norm_pert(:), H_pert_diag(:)
+  integer                        :: N_st, iter
+  N_st = N_states
+  allocate (pt2(N_st), norm_pert(N_st), H_pert_diag(N_st))
+  
+  pt2 = 1.d0
+  selection_criterion = 1.e-12
+  selection_criterion_min = 1.e-12
+  TOUCH selection_criterion_min selection_criterion selection_criterion_factor
+  call H_apply_mp2_selection(pt2, norm_pert, H_pert_diag,  N_st)
+  psi_det = psi_det_sorted
+  psi_coef = psi_coef_sorted
+  touch N_det psi_det psi_coef
+  print*,'N_det = ',N_det
+  print*,'-----'
+  print *,  'PT2                          = ', pt2(1)
+  print *,  'E                            = ', HF_energy
+  print *,  'E_before     +PT2            = ', HF_energy+pt2(1)
+  N_det = min(N_det,N_det_max)
+  call save_wavefunction
+  call ezfio_set_mp2_energy(HF_energy+pt2(1))
+  deallocate(pt2,norm_pert,H_pert_diag)
+end

plugins/MRCC_Utils/README.rst

@@ -1,171 +1,3 @@
-===========
-MRCC Module
-===========
-
-Needed Modules
-==============
-
-.. Do not edit this section. It was auto-generated from the
-.. by the `update_README.py` script.
-
-.. image:: tree_dependency.png
-
-* `Perturbation <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation>`_
-* `Selectors_full <http://github.com/LCPQ/quantum_package/tree/master/src/Selectors_full>`_
-* `Generators_full <http://github.com/LCPQ/quantum_package/tree/master/src/Generators_full>`_
-* `Psiref_Utils <http://github.com/LCPQ/quantum_package/tree/master/src/Psiref_Utils>`_
-
-Documentation
-=============
-
-.. Do not edit this section. It was auto-generated from the
-.. by the `update_README.py` script.
-
-`ci_eigenvectors_dressed <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/mrcc_utils.irp.f#L166>`_
-  Eigenvectors/values of the CI matrix
-
-
-`ci_eigenvectors_s2_dressed <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/mrcc_utils.irp.f#L167>`_
-  Eigenvectors/values of the CI matrix
-
-
-`ci_electronic_energy_dressed <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/mrcc_utils.irp.f#L165>`_
-  Eigenvectors/values of the CI matrix
-
-
-`ci_energy_dressed <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/mrcc_utils.irp.f#L232>`_
-  N_states lowest eigenvalues of the dressed CI matrix
-
-
-`davidson_diag_hjj_mrcc <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/davidson.irp.f#L56>`_
-  Davidson diagonalization with specific diagonal elements of the H matrix
-  .br
-  H_jj : specific diagonal H matrix elements to diagonalize de Davidson
-  .br
-  dets_in : bitmasks corresponding to determinants
-  .br
-  u_in : guess coefficients on the various states. Overwritten
-  on exit
-  .br
-  dim_in : leftmost dimension of u_in
-  .br
-  sze : Number of determinants
-  .br
-  N_st : Number of eigenstates
-  .br
-  iunit : Unit for the I/O
-  .br
-  Initial guess vectors are not necessarily orthonormal
-
-
-`davidson_diag_mrcc <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/davidson.irp.f#L1>`_
-  Davidson diagonalization.
-  .br
-  dets_in : bitmasks corresponding to determinants
-  .br
-  u_in : guess coefficients on the various states. Overwritten
-  on exit
-  .br
-  dim_in : leftmost dimension of u_in
-  .br
-  sze : Number of determinants
-  .br
-  N_st : Number of eigenstates
-  .br
-  iunit : Unit number for the I/O
-  .br
-  Initial guess vectors are not necessarily orthonormal
-
-
-`delta_ii <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/mrcc_utils.irp.f#L104>`_
-  Dressing matrix in N_det basis
-
-
-`delta_ij <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/mrcc_utils.irp.f#L103>`_
-  Dressing matrix in N_det basis
-
-
-`diagonalize_ci_dressed <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/mrcc_utils.irp.f#L247>`_
-  Replace the coefficients of the CI states by the coefficients of the
-  eigenstates of the CI matrix
-
-
-`find_triples_and_quadruples <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/mrcc_dress.irp.f#L206>`_
-  Undocumented
-
-
-`h_apply_mrcc <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/H_apply.irp.f_shell_27#L422>`_
-  Calls H_apply on the HF determinant and selects all connected single and double
-  excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
-
-
-`h_apply_mrcc_diexc <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/H_apply.irp.f_shell_27#L1>`_
-  Generate all double excitations of key_in using the bit masks of holes and
-  particles.
-  Assume N_int is already provided.
-
-
-`h_apply_mrcc_monoexc <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/H_apply.irp.f_shell_27#L273>`_
-  Generate all single excitations of key_in using the bit masks of holes and
-  particles.
-  Assume N_int is already provided.
-
-
-`h_matrix_dressed <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/mrcc_utils.irp.f#L140>`_
-  Dressed H with Delta_ij
-
-
-`h_u_0_mrcc <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/davidson.irp.f#L360>`_
-  Computes v_0 = H|u_0>
-  .br
-  n : number of determinants
-  .br
-  H_jj : array of <j|H|j>
-
-
-`lambda_mrcc <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/mrcc_utils.irp.f#L5>`_
-  cm/<Psi_0|H|D_m> or perturbative 1/Delta_E(m)
-
-
-`lambda_mrcc_tmp <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/mrcc_utils.irp.f#L81>`_
-  Undocumented
-
-
-`lambda_pert <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/mrcc_utils.irp.f#L6>`_
-  cm/<Psi_0|H|D_m> or perturbative 1/Delta_E(m)
-
-
-`mrcc_dress <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/mrcc_dress.irp.f#L15>`_
-  Undocumented
-
-
-`mrcc_dress_simple <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/mrcc_dress.irp.f#L160>`_
-  Undocumented
-
-
-`mrcc_iterations <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/mrcc_general.irp.f#L7>`_
-  Undocumented
-
-
-`oscillations <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/mrcc_utils.irp.f#L86>`_
-  Undocumented
-
-
-`pert_determinants <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/mrcc_utils.irp.f#L1>`_
-  Undocumented
-
-
-`psi_ref_lock <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/mrcc_dress.irp.f#L3>`_
-  Locks on ref determinants to fill delta_ij
-
-
-`run_mrcc <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/mrcc_general.irp.f#L1>`_
-  Undocumented
-
-
-`set_generators_bitmasks_as_holes_and_particles <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/mrcc_general.irp.f#L69>`_
-  Undocumented
-
 Needed Modules
 ==============
 .. Do not edit this section It was auto-generated
@@ -258,10 +90,6 @@ Documentation
   N_states lowest eigenvalues of the dressed CI matrix
 
 
-`create_minilist <http://github.com/LCPQ/quantum_package/tree/master/plugins/MRCC_Utils/mrcc_dress.irp.f#L17>`_
-  Undocumented
-
-
 `davidson_diag_hjj_mrcc <http://github.com/LCPQ/quantum_package/tree/master/plugins/MRCC_Utils/davidson.irp.f#L59>`_
   Davidson diagonalization with specific diagonal elements of the H matrix
   .br
@@ -374,7 +202,7 @@ Documentation
   Find A.C = B
 
 
-`find_triples_and_quadruples <http://github.com/LCPQ/quantum_package/tree/master/plugins/MRCC_Utils/mrcc_dress.irp.f#L271>`_
+`find_triples_and_quadruples <http://github.com/LCPQ/quantum_package/tree/master/plugins/MRCC_Utils/mrcc_dress.irp.f#L245>`_
   Undocumented
 
 
@@ -400,6 +228,22 @@ Documentation
   Undocumented
 
 
+`gen_det_idx <http://github.com/LCPQ/quantum_package/tree/master/plugins/MRCC_Utils/mrcc_dress.irp.f#L237>`_
+  Undocumented
+
+
+`gen_det_shortcut <http://github.com/LCPQ/quantum_package/tree/master/plugins/MRCC_Utils/mrcc_dress.irp.f#L235>`_
+  Undocumented
+
+
+`gen_det_sorted <http://github.com/LCPQ/quantum_package/tree/master/plugins/MRCC_Utils/mrcc_dress.irp.f#L234>`_
+  Undocumented
+
+
+`gen_det_version <http://github.com/LCPQ/quantum_package/tree/master/plugins/MRCC_Utils/mrcc_dress.irp.f#L236>`_
+  Undocumented
+
+
 `get_pseudo_inverse <http://github.com/LCPQ/quantum_package/tree/master/plugins/MRCC_Utils/LinearAlgebra.irp.f#L95>`_
   Find C = A^-1
 
@@ -435,26 +279,26 @@ Documentation
   Undocumented
 
 
-`h_apply_mrcc <http://github.com/LCPQ/quantum_package/tree/master/plugins/MRCC_Utils/H_apply.irp.f_shell_27#L553>`_
+h_apply_mrcc
   Calls H_apply on the HF determinant and selects all connected single and double
   excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
 
 
-`h_apply_mrcc_diexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/MRCC_Utils/H_apply.irp.f_shell_27#L3>`_
+h_apply_mrcc_diexc
   Undocumented
 
 
-`h_apply_mrcc_diexcorg <http://github.com/LCPQ/quantum_package/tree/master/plugins/MRCC_Utils/H_apply.irp.f_shell_27#L127>`_
+h_apply_mrcc_diexcorg
   Generate all double excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_mrcc_diexcp <http://github.com/LCPQ/quantum_package/tree/master/plugins/MRCC_Utils/H_apply.irp.f_shell_27#L100>`_
+h_apply_mrcc_diexcp
   Undocumented
 
 
-`h_apply_mrcc_monoexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/MRCC_Utils/H_apply.irp.f_shell_27#L399>`_
+h_apply_mrcc_monoexc
   Generate all single excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
@@ -464,15 +308,7 @@ Documentation
   Dressed H with Delta_ij
 
 
-`h_u_0_mrcc <http://github.com/LCPQ/quantum_package/tree/master/plugins/MRCC_Utils/davidson.irp.f#L371>`_
-  Computes v_0 = H|u_0>
-  .br
-  n : number of determinants
-  .br
-  H_jj : array of <j|H|j>
-
-
-`h_u_0_mrcc_org <http://github.com/LCPQ/quantum_package/tree/master/plugins/MRCC_Utils/davidson.irp.f#L475>`_
+`h_u_0_mrcc <http://github.com/LCPQ/quantum_package/tree/master/plugins/MRCC_Utils/davidson.irp.f#L367>`_
   Computes v_0 = H|u_0>
   .br
   n : number of determinants
@@ -781,10 +617,6 @@ Documentation
   Undocumented
 
 
-`mrcc_dress_simple <http://github.com/LCPQ/quantum_package/tree/master/plugins/MRCC_Utils/mrcc_dress.irp.f#L225>`_
-  Undocumented
-
-
 `mrcc_iterations <http://github.com/LCPQ/quantum_package/tree/master/plugins/MRCC_Utils/mrcc_general.irp.f#L7>`_
   Undocumented
 

plugins/MRCC_Utils/davidson.irp.f

@@ -105,10 +105,7 @@ subroutine davidson_diag_hjj_mrcc(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nin
   double precision               :: to_print(2,N_st)
   double precision               :: cpu, wall
   
-  integer(bit_kind)              :: dets_in_sorted(Nint,2,sze)
-  integer                        :: idx(sze), shortcut(0:sze+1),sh,ii,tmp
-  
-  PROVIDE det_connections
+  !PROVIDE det_connections
 
   call write_time(iunit)
   call wall_time(wall)
@@ -154,8 +151,6 @@ subroutine davidson_diag_hjj_mrcc(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nin
   ! ==============
   
   
-  dets_in_sorted(:,:,:) = dets_in(:,:,:)
-  call sort_dets_ab(dets_in_sorted, idx, shortcut, sze, Nint)
   
   k_pairs=0
   do l=1,N_st
@@ -215,8 +210,7 @@ subroutine davidson_diag_hjj_mrcc(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nin
       ! ----------------------
       
       do k=1,N_st
-        call H_u_0_mrcc(W(1,k,iter),U(1,k,iter),H_jj,sze,dets_in_sorted,shortcut,idx,Nint,istate)
-        !call H_u_0_mrcc_org(W(1,k,iter),U(1,k,iter),H_jj,sze,dets_in,Nint,istate)
+        call H_u_0_mrcc(W(1,k,iter),U(1,k,iter),H_jj,sze,dets_in,Nint,istate)
       enddo
       
       ! Compute h_kl = <u_k | W_l> = <u_k| H |u_l>
@@ -368,7 +362,9 @@ subroutine davidson_diag_hjj_mrcc(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nin
   abort_here = abort_all
 end
 
-subroutine H_u_0_mrcc(v_0,u_0,H_jj,n,keys_tmp,shortcut,sort_idx,Nint,istate)
+
+
+subroutine H_u_0_mrcc(v_0,u_0,H_jj,n,keys_tmp,Nint,istate)
   use bitmasks
   implicit none
   BEGIN_DOC
@@ -389,8 +385,11 @@ subroutine H_u_0_mrcc(v_0,u_0,H_jj,n,keys_tmp,shortcut,sort_idx,Nint,istate)
   integer                        :: i,j,k,l, jj,ii
   integer                        :: i0, j0
   
-  integer,intent(in)             :: shortcut(0:n+1), sort_idx(n)
-  integer                        :: tmp, warp(2,0:n+1), sh, ni
+  integer                        :: shortcut(0:n+1), sort_idx(n)
+  integer(bit_kind)              :: sorted(Nint,n), version(Nint,n)
+  
+  
+  integer                        :: sh, sh2, ni, exa, ext, org_i, org_j, endi, pass
 !   
   
   ASSERT (Nint > 0)
@@ -401,8 +400,10 @@ subroutine H_u_0_mrcc(v_0,u_0,H_jj,n,keys_tmp,shortcut,sort_idx,Nint,istate)
   
  
    !$OMP PARALLEL DEFAULT(NONE)                                       &
-      !$OMP PRIVATE(i,hij,j,k,idx,jj,vt,ii,warp,tmp,sh)                             &
-      !$OMP SHARED(n_det_ref,n_det_non_ref,idx_ref,idx_non_ref,n,H_jj,u_0,keys_tmp,Nint,v_0,istate,delta_ij,shortcut,sort_idx)
+      !$OMP PRIVATE(i,hij,j,k,idx,jj,vt,ii,sh, sh2, ni, exa, ext, org_i, org_j, endi, pass)                             &
+      !$OMP SHARED(n_det_ref,n_det_non_ref,idx_ref,idx_non_ref,n,H_jj,u_0,keys_tmp,Nint,v_0,istate,delta_ij,sorted,shortcut,sort_idx,version)
+      
+      
       
  !$OMP DO SCHEDULE(static)
   do i=1,n  
@@ -412,43 +413,73 @@ subroutine H_u_0_mrcc(v_0,u_0,H_jj,n,keys_tmp,shortcut,sort_idx,Nint,istate)
 
   allocate(idx(0:n), vt(n))
   Vt = 0.d0
+  
+  
+  !$OMP SINGLE
+  call sort_dets_ab_v(keys_tmp, sorted, sort_idx, shortcut, version, n, Nint)
+  !$OMP END SINGLE
+  
+
   !$OMP DO SCHEDULE(dynamic)
   do sh=1,shortcut(0)
-    warp(1,0) = 0
-    do ii=1,sh!shortcut(0)
-      tmp = 0
+  do sh2=1,sh
+    exa = 0
     do ni=1,Nint
-        tmp = popcnt(xor(keys_tmp(ni,1, shortcut(ii)), keys_tmp(ni,1,shortcut(sh))))
+      exa += popcnt(xor(version(ni,sh), version(ni,sh2)))
     end do
-      if(tmp <= 4) then
-        tmp = warp(1,0) + 1
-        warp(1,0) = tmp
-        warp(1,tmp) = shortcut(ii)
-        warp(2,tmp) = shortcut(ii+1)-1
+    if(exa > 2) then
+      cycle
+    end if
+    
+    do i=shortcut(sh),shortcut(sh+1)-1
+      if(sh==sh2) then
+        endi = i-1
+      else
+        endi = shortcut(sh2+1)-1
+      end if
+      
+      do j=shortcut(sh2),endi
+        ext = exa
+        do ni=1,Nint
+          ext += popcnt(xor(sorted(ni,i), sorted(ni,j)))
+        end do
+        if(ext <= 4) then
+          org_i = sort_idx(i)
+          org_j = sort_idx(j)
+          
+          call i_H_j(keys_tmp(1,1,org_j),keys_tmp(1,1,org_i),Nint,hij)
+          vt (org_i) = vt (org_i) + hij*u_0(org_j)
+          vt (org_j) = vt (org_j) + hij*u_0(org_i)
         end if
       end do
-    
-    do ii=shortcut(sh),shortcut(sh+1)-1
-      idx(0) = ii
-      
-      
-      !call filter_connected_davidson_mwen(keys_tmp,shortcut,keys_tmp(1,1,ii),Nint,ii-1,idx)
-      call filter_connected_davidson_warp(keys_tmp,warp,keys_tmp(1,1,ii),Nint,ii-1,idx)
-      i = sort_idx(ii)
-      
-      do jj=1,idx(0)
-        j = sort_idx(idx(jj))
-        !j = idx(jj)
-        if ( (dabs(u_0(j)) > 1.d-7).or.((dabs(u_0(i)) > 1.d-7)) ) then
-          call i_H_j(keys_tmp(1,1,idx(jj)),keys_tmp(1,1,ii),Nint,hij)
-          vt (i) = vt (i) + hij*u_0(j)
-          vt (j) = vt (j) + hij*u_0(i)
-        endif
     end do
   end do
   enddo
   !$OMP END DO
   
+  !$OMP SINGLE
+  call sort_dets_ba_v(keys_tmp, sorted, sort_idx, shortcut, version, n, Nint)
+  !$OMP END SINGLE
+ 
+  !$OMP DO SCHEDULE(dynamic)
+  do sh=1,shortcut(0)
+    do i=shortcut(sh),shortcut(sh+1)-1
+      do j=shortcut(sh),i-1
+        ext = 0
+        do ni=1,Nint
+          ext += popcnt(xor(sorted(ni,i), sorted(ni,j)))
+        end do
+        if(ext == 4) then
+          org_i = sort_idx(i)
+          org_j = sort_idx(j)
+          call i_H_j(keys_tmp(1,1,org_j),keys_tmp(1,1,org_i),Nint,hij)
+          vt (org_i) = vt (org_i) + hij*u_0(org_j)
+          vt (org_j) = vt (org_j) + hij*u_0(org_i)
+        end if
+      end do
+    end do
+  enddo
+  !$OMP END DO
   
   
   !$OMP DO SCHEDULE(guided)
@@ -470,80 +501,3 @@ subroutine H_u_0_mrcc(v_0,u_0,H_jj,n,keys_tmp,shortcut,sort_idx,Nint,istate)
   !$OMP END PARALLEL
 end
 
-
-
-subroutine H_u_0_mrcc_org(v_0,u_0,H_jj,n,keys_tmp,Nint,istate)
-  use bitmasks
-  implicit none
-  BEGIN_DOC
-  ! Computes v_0 = H|u_0>
-  !   
-  ! n : number of determinants
-  !     
-  ! H_jj : array of <j|H|j>
-  END_DOC
-  integer, intent(in)            :: n,Nint,istate
-  double precision, intent(out)  :: v_0(n)
-  double precision, intent(in)   :: u_0(n)
-  double precision, intent(in)   :: H_jj(n)
-  integer(bit_kind),intent(in)   :: keys_tmp(Nint,2,n)
-  integer, allocatable           :: idx(:)
-  double precision               :: hij
-  double precision, allocatable  :: vt(:)
-  integer                        :: i,j,k,l, jj,ii
-  integer                        :: i0, j0
-  
-
-  
-
-  
-  ASSERT (Nint > 0)
-  ASSERT (Nint == N_int)
-  ASSERT (n>0)
-  PROVIDE ref_bitmask_energy delta_ij 
-  integer, parameter             :: block_size = 157
-  !$OMP PARALLEL DEFAULT(NONE)                                       &
-      !$OMP PRIVATE(i,hij,j,k,idx,jj,ii,vt)                             &
-      !$OMP SHARED(n_det_ref,n_det_non_ref,idx_ref,idx_non_ref,n,H_jj,u_0,keys_tmp,Nint,v_0,istate,delta_ij)
-  !$OMP DO SCHEDULE(static)
-  do i=1,n  
-    v_0(i) = H_jj(i) * u_0(i)
-  enddo 
-  !$OMP END DO
-  allocate(idx(0:n), vt(n))
-  Vt = 0.d0
-  !$OMP DO SCHEDULE(guided)
-  do i=1,n
-    idx(0) = i
-    call filter_connected_davidson(keys_tmp,keys_tmp(1,1,i),Nint,i-1,idx)
-    do jj=1,idx(0)
-      j = idx(jj)
-!     if ( (dabs(u_0(j)) > 1.d-7).or.((dabs(u_0(i)) > 1.d-7)) ) then
-        call i_H_j(keys_tmp(1,1,j),keys_tmp(1,1,i),Nint,hij)
-        hij = hij 
-        vt (i) = vt (i) + hij*u_0(j)
-        vt (j) = vt (j) + hij*u_0(i)
-!     endif
-    enddo
-  enddo
-  !$OMP END DO
-
-  !$OMP DO SCHEDULE(guided)
-  do ii=1,n_det_ref
-    i = idx_ref(ii)
-    do jj = 1, n_det_non_ref
-        j = idx_non_ref(jj)
-        vt (i) = vt (i) + delta_ij(ii,jj,istate)*u_0(j)
-        vt (j) = vt (j) + delta_ij(ii,jj,istate)*u_0(i)
-    enddo
-  enddo
-  !$OMP END DO
-  !$OMP CRITICAL
-  do i=1,n
-    v_0(i) = v_0(i) + vt(i)
-  enddo
-  !$OMP END CRITICAL
-  deallocate(idx,vt)
-  !$OMP END PARALLEL
-end
-

plugins/MRCC_Utils/mrcc_dress.irp.f

@@ -14,54 +14,6 @@ BEGIN_PROVIDER [ integer(omp_lock_kind), psi_ref_lock, (psi_det_size) ]
 END_PROVIDER
 
 
-subroutine create_minilist(key_mask, fullList, miniList, idx_miniList, N_fullList, N_miniList, Nint)
-  use bitmasks
-  implicit none
-  
-  integer(bit_kind), intent(in)            :: fullList(Nint, 2, N_fullList)
-  integer, intent(in)                      :: N_fullList
-  integer(bit_kind),intent(out)            :: miniList(Nint, 2, N_fullList)
-  integer,intent(out)                      :: idx_miniList(N_fullList), N_miniList
-  integer, intent(in)                      :: Nint
-  integer(bit_kind)                        :: key_mask(Nint, 2)
-  integer                                  :: ni, i, n_a, n_b, e_a, e_b
-  
-  
-  n_a = 0
-  n_b = 0
-  do ni=1,nint
-    n_a = n_a + popcnt(key_mask(ni,1))
-    n_b = n_b + popcnt(key_mask(ni,2))
-  end do
-  
-  if(n_a == 0) then
-    N_miniList = N_fullList
-    miniList(:,:,:) = fullList(:,:,:)
-    do i=1,N_fullList
-      idx_miniList(i) = i
-    end do
-    return
-  end if
-  
-  N_miniList = 0
-  
-  do i=1,N_fullList
-    e_a = n_a
-    e_b = n_b
-    do ni=1,nint
-      e_a -= popcnt(iand(fullList(ni, 1, i), key_mask(ni, 1)))
-      e_b -= popcnt(iand(fullList(ni, 2, i), key_mask(ni, 2)))
-    end do
-    
-    if(e_a + e_b <= 2) then
-      N_miniList = N_miniList + 1
-      miniList(:,:,N_miniList) = fullList(:,:,i)
-      idx_miniList(N_miniList) = i
-    end if
-  end do
-end subroutine
-
-
 subroutine mrcc_dress(delta_ij_, delta_ii_, Ndet_ref, Ndet_non_ref,i_generator,n_selected,det_buffer,Nint,iproc,key_mask)
  use bitmasks
  implicit none
@@ -75,11 +27,10 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Ndet_ref, Ndet_non_ref,i_generator,n
   integer                        :: i,j,k,l
   integer                        :: degree_alpha(psi_det_size)
   integer                        :: idx_alpha(0:psi_det_size)
-  logical                        :: good
+  logical                        :: good, fullMatch
 
   integer(bit_kind)              :: tq(Nint,2,n_selected)
   integer                        :: N_tq, c_ref ,degree
-  integer                        :: connected_to_ref
 
   double precision               :: hIk, hla, hIl, dIk(N_states), dka(N_states), dIa(N_states)
   double precision, allocatable  :: dIa_hla(:,:)
@@ -91,11 +42,24 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Ndet_ref, Ndet_non_ref,i_generator,n
   integer                        :: iint, ipos
   integer                        :: i_state, k_sd, l_sd, i_I, i_alpha
   
-  integer(bit_kind)              :: miniList(Nint, 2, N_det_non_ref), key_mask(Nint, 2)
-  integer                        :: idx_miniList(N_det_non_ref), N_miniList
+  integer(bit_kind),allocatable  :: miniList(:,:,:)
+  integer(bit_kind),intent(in)   :: key_mask(Nint, 2)
+  integer,allocatable            :: idx_miniList(:)
+  integer                        :: N_miniList, ni, leng
   
   
-  call find_triples_and_quadruples(i_generator,n_selected,det_buffer,Nint,tq,N_tq)
+  leng = max(N_det_generators, N_det_non_ref)
+  allocate(miniList(Nint, 2, leng), idx_miniList(leng))
+  
+  !create_minilist_find_previous(key_mask, fullList, miniList, N_fullList, N_miniList, fullMatch, Nint)
+  call create_minilist_find_previous(key_mask, psi_det_generators, miniList, i_generator-1, N_miniList, fullMatch, Nint)
+  
+  if(fullMatch) then
+    return
+  end if
+  
+  
+  call find_triples_and_quadruples(i_generator,n_selected,det_buffer,Nint,tq,N_tq,miniList,N_minilist)
 
   allocate (dIa_hla(N_states,Ndet_non_ref))
 
@@ -214,61 +178,24 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Ndet_ref, Ndet_non_ref,i_generator,n
     enddo
   enddo
   deallocate (dIa_hla)
+  deallocate(miniList, idx_miniList)
 end
 
 
 
+ BEGIN_PROVIDER [ integer(bit_kind), gen_det_sorted,  (N_int,2,N_det_generators,2) ]
+&BEGIN_PROVIDER [ integer, gen_det_shortcut, (0:N_det_generators,2) ]
+&BEGIN_PROVIDER [ integer, gen_det_version, (N_int, N_det_generators,2) ]
+&BEGIN_PROVIDER [ integer, gen_det_idx, (N_det_generators,2) ]
+  gen_det_sorted(:,:,:,1) = psi_det_generators(:,:,:N_det_generators)
+  gen_det_sorted(:,:,:,2) = psi_det_generators(:,:,:N_det_generators)
+  call sort_dets_ab_v(gen_det_sorted(:,:,:,1), gen_det_idx(:,1), gen_det_shortcut(0:,1), gen_det_version(:,:,1), N_det_generators, N_int)
+  call sort_dets_ba_v(gen_det_sorted(:,:,:,2), gen_det_idx(:,2), gen_det_shortcut(0:,2), gen_det_version(:,:,2), N_det_generators, N_int)
+END_PROVIDER
 
 
+subroutine find_triples_and_quadruples(i_generator,n_selected,det_buffer,Nint,tq,N_tq,miniList,N_miniList)
 
-
-subroutine mrcc_dress_simple(delta_ij_non_ref_,Ndet_non_ref,i_generator,n_selected,det_buffer,Nint,iproc)
- use bitmasks
- implicit none
-
-  integer, intent(in)            :: i_generator,n_selected, Nint, iproc
-  integer, intent(in) :: Ndet_non_ref
-  double precision, intent(inout) :: delta_ij_non_ref_(Ndet_non_ref,Ndet_non_ref,*)
-
-  integer(bit_kind), intent(in)  :: det_buffer(Nint,2,n_selected)
-  integer                        :: i,j,k,m
-  integer                        :: new_size
-  integer                        :: degree(psi_det_size)
-  integer                        :: idx(0:psi_det_size)
-  logical                        :: good
-
-  integer(bit_kind)              :: tq(Nint,2,n_selected)
-  integer                        :: N_tq, c_ref
-  integer                        :: connected_to_ref
-
-  call find_triples_and_quadruples(i_generator,n_selected,det_buffer,Nint,tq,N_tq)
-
-  ! Compute <k|H|a><a|H|j> / (E0 - Haa)
-  double precision :: hka, haa
-  double precision :: haj
-  double precision :: f(N_states)
-
-  do i=1,N_tq
-    call get_excitation_degree_vector(psi_non_ref,tq(1,1,i),degree,Nint,Ndet_non_ref,idx)
-    call i_h_j(tq(1,1,i),tq(1,1,i),Nint,haa)
-    do m=1,N_states
-      f(m) = 1.d0/(ci_electronic_energy(m)-haa)
-    enddo
-    do k=1,idx(0)
-      call i_h_j(tq(1,1,i),psi_non_ref(1,1,idx(k)),Nint,hka)
-      do j=k,idx(0)
-        call i_h_j(tq(1,1,i),psi_non_ref(1,1,idx(j)),Nint,haj)
-        do m=1,N_states
-          delta_ij_non_ref_(idx(k), idx(j),m) += haj*hka* f(m)
-          delta_ij_non_ref_(idx(j), idx(k),m) += haj*hka* f(m)
-        enddo
-      enddo 
-    enddo
-  enddo
-end
-
-
-subroutine find_triples_and_quadruples(i_generator,n_selected,det_buffer,Nint,tq,N_tq)
  use bitmasks
  implicit none
 
@@ -283,16 +210,22 @@ subroutine find_triples_and_quadruples(i_generator,n_selected,det_buffer,Nint,tq
 
   integer(bit_kind), intent(out) :: tq(Nint,2,n_selected)
   integer, intent(out)           :: N_tq
-  integer                        :: c_ref
-  integer                        :: connected_to_ref
+  
+  
+  integer                        :: nt,ni
+  logical, external              :: is_connected_to
+  
+  
+  integer(bit_kind),intent(in)  :: miniList(Nint,2,N_det_generators)
+  integer,intent(in)            :: N_miniList
+
   
   
   N_tq = 0
-  do i=1,N_selected
-    c_ref = connected_to_ref(det_buffer(1,1,i),psi_det_generators,Nint, &
-       i_generator,N_det_generators)
   
-    if (c_ref /= 0) then
+  
+  i_loop : do i=1,N_selected
+    if(is_connected_to(det_buffer(1,1,i), miniList, Nint, N_miniList)) then
       cycle
     end if
 
@@ -316,8 +249,7 @@ subroutine find_triples_and_quadruples(i_generator,n_selected,det_buffer,Nint,tq
         enddo
       endif
     endif
-  enddo
-
+  enddo i_loop
 end
 
 

plugins/Molden/README.rst

@@ -38,3 +38,39 @@ Needed Modules
 * `MO_Basis <http://github.com/LCPQ/quantum_package/tree/master/src/MO_Basis>`_
 * `Utils <http://github.com/LCPQ/quantum_package/tree/master/src/Utils>`_
 
+Needed Modules
+==============
+.. Do not edit this section It was auto-generated
+.. by the `update_README.py` script.
+
+
+.. image:: tree_dependency.png
+
+* `MO_Basis <http://github.com/LCPQ/quantum_package/tree/master/src/MO_Basis>`_
+* `Utils <http://github.com/LCPQ/quantum_package/tree/master/src/Utils>`_
+
+Documentation
+=============
+.. Do not edit this section It was auto-generated
+.. by the `update_README.py` script.
+
+
+`print_mos <http://github.com/LCPQ/quantum_package/tree/master/plugins/Molden/print_mo.irp.f#L1>`_
+  Undocumented
+
+
+`write_ao_basis <http://github.com/LCPQ/quantum_package/tree/master/plugins/Molden/print_mo.irp.f#L63>`_
+  Undocumented
+
+
+`write_geometry <http://github.com/LCPQ/quantum_package/tree/master/plugins/Molden/print_mo.irp.f#L45>`_
+  Undocumented
+
+
+`write_intro_gamess <http://github.com/LCPQ/quantum_package/tree/master/plugins/Molden/print_mo.irp.f#L26>`_
+  Undocumented
+
+
+`write_mo_basis <http://github.com/LCPQ/quantum_package/tree/master/plugins/Molden/print_mo.irp.f#L112>`_
+  Undocumented
+

plugins/Perturbation/Moller_plesset.irp.f

@@ -1,42 +0,0 @@
-subroutine pt2_moller_plesset(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,n_st)
-  use bitmasks
-  implicit none
-  integer, intent(in)            :: Nint,ndet,n_st
-  integer(bit_kind), intent(in)  :: det_pert(Nint,2)
-  double precision , intent(out) :: c_pert(n_st),e_2_pert(n_st),H_pert_diag(N_st)
-  double precision               :: i_H_psi_array(N_st)
-  
-  BEGIN_DOC
-  ! compute the standard Moller-Plesset perturbative first order coefficient and second order energetic contribution
-  !
-  ! for the various n_st states.
-  !
-  ! c_pert(i) = <psi(i)|H|det_pert>/(difference of orbital energies) 
-  !
-  ! e_2_pert(i) = <psi(i)|H|det_pert>^2/(difference of orbital energies) 
-  !
-  END_DOC
-  
-  integer                        :: i,j
-  double precision               :: diag_H_mat_elem
-  integer                        :: exc(0:2,2,2)
-  integer                        :: degree
-  double precision               :: phase,delta_e,h
-  integer                        :: h1,h2,p1,p2,s1,s2
-  ASSERT (Nint == N_int)
-  ASSERT (Nint > 0)
-  call get_excitation(ref_bitmask,det_pert,exc,degree,phase,Nint)
-  call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
-  delta_e = Fock_matrix_diag_mo(h1) + Fock_matrix_diag_mo(h2) - &
-           (Fock_matrix_diag_mo(p1) + Fock_matrix_diag_mo(p2))
-  delta_e = 1.d0/delta_e
-
-  call i_H_psi(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det,psi_selectors_size,n_st,i_H_psi_array)
-  h = diag_H_mat_elem(det_pert,Nint)
-  do i =1,n_st
-    H_pert_diag(i) = h
-    c_pert(i) = i_H_psi_array(i) *delta_e
-    e_2_pert(i) = c_pert(i) * i_H_psi_array(i)
-  enddo
-  
-end

plugins/Perturbation/README.rst

@@ -107,92 +107,92 @@ Documentation
   Undocumented
 
 
-`perturb_buffer_by_mono_delta_rho_one_point <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation/perturbation.irp.f_shell_13#L161>`_
+perturb_buffer_by_mono_delta_rho_one_point
   Applly pertubration ``delta_rho_one_point`` to the buffer of determinants generated in the H_apply
   routine.
 
 
-`perturb_buffer_by_mono_dipole_moment_z <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation/perturbation.irp.f_shell_13#L266>`_
+perturb_buffer_by_mono_dipole_moment_z
   Applly pertubration ``dipole_moment_z`` to the buffer of determinants generated in the H_apply
   routine.
 
 
-`perturb_buffer_by_mono_epstein_nesbet <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation/perturbation.irp.f_shell_13#L371>`_
+perturb_buffer_by_mono_epstein_nesbet
   Applly pertubration ``epstein_nesbet`` to the buffer of determinants generated in the H_apply
   routine.
 
 
-`perturb_buffer_by_mono_epstein_nesbet_2x2 <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation/perturbation.irp.f_shell_13#L476>`_
+perturb_buffer_by_mono_epstein_nesbet_2x2
   Applly pertubration ``epstein_nesbet_2x2`` to the buffer of determinants generated in the H_apply
   routine.
 
 
-`perturb_buffer_by_mono_epstein_nesbet_sc2 <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation/perturbation.irp.f_shell_13#L791>`_
+perturb_buffer_by_mono_epstein_nesbet_sc2
   Applly pertubration ``epstein_nesbet_sc2`` to the buffer of determinants generated in the H_apply
   routine.
 
 
-`perturb_buffer_by_mono_epstein_nesbet_sc2_no_projected <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation/perturbation.irp.f_shell_13#L686>`_
+perturb_buffer_by_mono_epstein_nesbet_sc2_no_projected
   Applly pertubration ``epstein_nesbet_sc2_no_projected`` to the buffer of determinants generated in the H_apply
   routine.
 
 
-`perturb_buffer_by_mono_epstein_nesbet_sc2_projected <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation/perturbation.irp.f_shell_13#L581>`_
+perturb_buffer_by_mono_epstein_nesbet_sc2_projected
   Applly pertubration ``epstein_nesbet_sc2_projected`` to the buffer of determinants generated in the H_apply
   routine.
 
 
-`perturb_buffer_by_mono_h_core <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation/perturbation.irp.f_shell_13#L896>`_
+perturb_buffer_by_mono_h_core
   Applly pertubration ``h_core`` to the buffer of determinants generated in the H_apply
   routine.
 
 
-`perturb_buffer_by_mono_moller_plesset <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation/perturbation.irp.f_shell_13#L56>`_
+perturb_buffer_by_mono_moller_plesset
   Applly pertubration ``moller_plesset`` to the buffer of determinants generated in the H_apply
   routine.
 
 
-`perturb_buffer_delta_rho_one_point <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation/perturbation.irp.f_shell_13#L110>`_
+perturb_buffer_delta_rho_one_point
   Applly pertubration ``delta_rho_one_point`` to the buffer of determinants generated in the H_apply
   routine.
 
 
-`perturb_buffer_dipole_moment_z <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation/perturbation.irp.f_shell_13#L215>`_
+perturb_buffer_dipole_moment_z
   Applly pertubration ``dipole_moment_z`` to the buffer of determinants generated in the H_apply
   routine.
 
 
-`perturb_buffer_epstein_nesbet <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation/perturbation.irp.f_shell_13#L320>`_
+perturb_buffer_epstein_nesbet
   Applly pertubration ``epstein_nesbet`` to the buffer of determinants generated in the H_apply
   routine.
 
 
-`perturb_buffer_epstein_nesbet_2x2 <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation/perturbation.irp.f_shell_13#L425>`_
+perturb_buffer_epstein_nesbet_2x2
   Applly pertubration ``epstein_nesbet_2x2`` to the buffer of determinants generated in the H_apply
   routine.
 
 
-`perturb_buffer_epstein_nesbet_sc2 <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation/perturbation.irp.f_shell_13#L740>`_
+perturb_buffer_epstein_nesbet_sc2
   Applly pertubration ``epstein_nesbet_sc2`` to the buffer of determinants generated in the H_apply
   routine.
 
 
-`perturb_buffer_epstein_nesbet_sc2_no_projected <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation/perturbation.irp.f_shell_13#L635>`_
+perturb_buffer_epstein_nesbet_sc2_no_projected
   Applly pertubration ``epstein_nesbet_sc2_no_projected`` to the buffer of determinants generated in the H_apply
   routine.
 
 
-`perturb_buffer_epstein_nesbet_sc2_projected <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation/perturbation.irp.f_shell_13#L530>`_
+perturb_buffer_epstein_nesbet_sc2_projected
   Applly pertubration ``epstein_nesbet_sc2_projected`` to the buffer of determinants generated in the H_apply
   routine.
 
 
-`perturb_buffer_h_core <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation/perturbation.irp.f_shell_13#L845>`_
+perturb_buffer_h_core
   Applly pertubration ``h_core`` to the buffer of determinants generated in the H_apply
   routine.
 
 
-`perturb_buffer_moller_plesset <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation/perturbation.irp.f_shell_13#L5>`_
+perturb_buffer_moller_plesset
   Applly pertubration ``moller_plesset`` to the buffer of determinants generated in the H_apply
   routine.
 
@@ -239,7 +239,7 @@ Documentation
   .br
 
 
-`pt2_epstein_nesbet <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation/epstein_nesbet.irp.f#L1>`_
+`pt2_epstein_nesbet <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation/pt2_equations.irp.f_template_363#L3>`_
   compute the standard Epstein-Nesbet perturbative first order coefficient and second order energetic contribution
   .br
   for the various N_st states.
@@ -250,7 +250,7 @@ Documentation
   .br
 
 
-`pt2_epstein_nesbet_2x2 <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation/epstein_nesbet.irp.f#L45>`_
+`pt2_epstein_nesbet_2x2 <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation/pt2_equations.irp.f_template_363#L59>`_
   compute the Epstein-Nesbet 2x2 diagonalization coefficient and energetic contribution
   .br
   for the various N_st states.
@@ -261,7 +261,7 @@ Documentation
   .br
 
 
-`pt2_epstein_nesbet_sc2 <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation/pert_sc2.irp.f#L186>`_
+`pt2_epstein_nesbet_sc2 <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation/pt2_equations.irp.f_template_363#L353>`_
   compute the standard Epstein-Nesbet perturbative first order coefficient and second order energetic contribution
   .br
   for the various N_st states, but with the CISD_SC2 energies and coefficients
@@ -272,7 +272,7 @@ Documentation
   .br
 
 
-`pt2_epstein_nesbet_sc2_no_projected <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation/pert_sc2.irp.f#L87>`_
+`pt2_epstein_nesbet_sc2_no_projected <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation/pt2_equations.irp.f_template_363#L275>`_
   compute the Epstein-Nesbet perturbative first order coefficient and second order energetic contribution
   .br
   for the various N_st states,
@@ -296,7 +296,7 @@ Documentation
   H_pert_diag = <HF|H|det_pert> c_pert
 
 
-`pt2_epstein_nesbet_sc2_projected <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation/pert_sc2.irp.f#L2>`_
+`pt2_epstein_nesbet_sc2_projected <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation/pt2_equations.irp.f_template_363#L181>`_
   compute the Epstein-Nesbet perturbative first order coefficient and second order energetic contribution
   .br
   for the various N_st states,
@@ -336,7 +336,7 @@ Documentation
   than pt2_max in absolute value
 
 
-`pt2_moller_plesset <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation/Moller_plesset.irp.f#L1>`_
+`pt2_moller_plesset <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation/pt2_equations.irp.f_template_363#L120>`_
   compute the standard Moller-Plesset perturbative first order coefficient and second order energetic contribution
   .br
   for the various n_st states.
@@ -352,7 +352,7 @@ Documentation
   provided.
 
 
-`repeat_all_e_corr <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation/pert_sc2.irp.f#L156>`_
+`repeat_all_e_corr <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation/pert_sc2.irp.f#L1>`_
   Undocumented
 
 

plugins/Perturbation/delta_rho_perturbation.irp.f

@@ -1,4 +1,4 @@
-subroutine pt2_delta_rho_one_point(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,n_st)
+subroutine pt2_delta_rho_one_point(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,n_st,minilist,idx_minilist,N_minilist)
   use bitmasks
   implicit none
   integer, intent(in)            :: Nint,ndet,n_st
@@ -7,6 +7,10 @@ subroutine pt2_delta_rho_one_point(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,nde
   double precision               :: i_O1_psi_array(N_st)
   double precision               :: i_H_psi_array(N_st)
   
+  integer, intent(in)            :: N_minilist
+  integer, intent(in)            :: idx_minilist(0:N_det_selectors)
+  integer(bit_kind), intent(in)  :: minilist(Nint,2,N_det_selectors)
+  
   BEGIN_DOC
   ! compute the perturbatibe contribution to the Integrated Spin density at z = z_one point of one determinant
   !
@@ -46,7 +50,8 @@ subroutine pt2_delta_rho_one_point(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,nde
 ! endif
   call i_O1_psi_alpha_beta(mo_integrated_delta_rho_one_point,det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_O1_psi_array)
 
-  call i_H_psi(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array)
+  !call i_H_psi(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array)
+  call i_H_psi_minilist(det_pert,minilist,idx_minilist,N_minilist,psi_selectors_coef,Nint,N_minilist,psi_selectors_size,N_st,i_H_psi_array)
 
   h   = diag_H_mat_elem(det_pert,Nint)
   oii = diag_O1_mat_elem_alpha_beta(mo_integrated_delta_rho_one_point,det_pert,N_int)

plugins/Perturbation/dipole_moment.irp.f

@@ -1,4 +1,4 @@
-subroutine pt2_dipole_moment_z(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,n_st)
+subroutine pt2_dipole_moment_z(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,n_st,minilist,idx_minilist,N_minilist)
   use bitmasks
   implicit none
   integer, intent(in)            :: Nint,ndet,n_st
@@ -7,6 +7,10 @@ subroutine pt2_dipole_moment_z(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,n_
   double precision               :: i_O1_psi_array(N_st)
   double precision               :: i_H_psi_array(N_st)
   
+  integer, intent(in)            :: N_minilist
+  integer, intent(in)            :: idx_minilist(0:N_det_selectors)
+  integer(bit_kind), intent(in)  :: minilist(Nint,2,N_det_selectors)
+  
   BEGIN_DOC
   ! compute the perturbatibe contribution to the dipole moment of one determinant
   !
@@ -46,7 +50,9 @@ subroutine pt2_dipole_moment_z(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,n_
 ! endif
 
   call i_O1_psi(mo_dipole_z,det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_O1_psi_array)
-  call i_H_psi(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array)
+  !call i_H_psi(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array)
+  call i_H_psi_minilist(det_pert,minilist,idx_minilist,N_minilist,psi_selectors_coef,Nint,N_minilist,psi_selectors_size,N_st,i_H_psi_array)
+  
   h = diag_H_mat_elem(det_pert,Nint)
   oii = diag_O1_mat_elem(mo_dipole_z,det_pert,N_int)
 

plugins/Perturbation/epstein_nesbet.irp.f

@@ -1,93 +0,0 @@
-subroutine pt2_epstein_nesbet(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,N_st)
-  use bitmasks
-  implicit none
-  integer, intent(in)            :: Nint,ndet,N_st
-  integer(bit_kind), intent(in)  :: det_pert(Nint,2)
-  double precision , intent(out) :: c_pert(N_st),e_2_pert(N_st),H_pert_diag(N_st)
-  double precision               :: i_H_psi_array(N_st)
-  
-  BEGIN_DOC
-  ! compute the standard Epstein-Nesbet perturbative first order coefficient and second order energetic contribution
-  !
-  ! for the various N_st states.
-  !
-  ! c_pert(i) = <psi(i)|H|det_pert>/( E(i) - <det_pert|H|det_pert> )
-  !
-  ! e_2_pert(i) = <psi(i)|H|det_pert>^2/( E(i) - <det_pert|H|det_pert> )
-  !
-  END_DOC
-  
-  integer                        :: i,j
-  double precision               :: diag_H_mat_elem, h
-  PROVIDE  selection_criterion
-
-  ASSERT (Nint == N_int)
-  ASSERT (Nint > 0)
-  call i_H_psi(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array)
-  h = diag_H_mat_elem(det_pert,Nint)
-  do i =1,N_st
-    if(CI_electronic_energy(i)>h.and.CI_electronic_energy(i).ne.0.d0)then
-      c_pert(i) = -1.d0
-      e_2_pert(i) = selection_criterion*selection_criterion_factor*2.d0
-    else if  (dabs(CI_electronic_energy(i) - h) > 1.d-6) then
-        c_pert(i) = i_H_psi_array(i) / (CI_electronic_energy(i) - h)
-        H_pert_diag(i) = h*c_pert(i)*c_pert(i)
-        e_2_pert(i) = c_pert(i) * i_H_psi_array(i)
-    else
-      c_pert(i) = -1.d0
-      e_2_pert(i) = -dabs(i_H_psi_array(i))
-      H_pert_diag(i) = h
-    endif
-  enddo
-  
-end
-
-subroutine pt2_epstein_nesbet_2x2(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,N_st)
-  use bitmasks
-  implicit none
-  integer, intent(in)            :: Nint,ndet,N_st
-  integer(bit_kind), intent(in)  :: det_pert(Nint,2)
-  double precision , intent(out) :: c_pert(N_st),e_2_pert(N_st),H_pert_diag(N_st)
-  double precision               :: i_H_psi_array(N_st)
-  
-  BEGIN_DOC
-  ! compute the Epstein-Nesbet 2x2 diagonalization coefficient and energetic contribution
-  !
-  ! for the various N_st states.
-  !
-  ! e_2_pert(i) = 0.5 * (( <det_pert|H|det_pert> -  E(i) )  - sqrt( ( <det_pert|H|det_pert> -  E(i)) ^2 + 4 <psi(i)|H|det_pert>^2  )
-  !
-  ! c_pert(i) = e_2_pert(i)/ <psi(i)|H|det_pert>
-  !
-  END_DOC
-  
-  integer                        :: i,j
-  double precision               :: diag_H_mat_elem,delta_e, h
-  ASSERT (Nint == N_int)
-  ASSERT (Nint > 0)
-  PROVIDE CI_electronic_energy
-
-  call i_H_psi(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array)
-  h = diag_H_mat_elem(det_pert,Nint)
-  do i =1,N_st
-    if (i_H_psi_array(i) /= 0.d0) then
-      delta_e = h - CI_electronic_energy(i)
-      if (delta_e > 0.d0) then
-        e_2_pert(i) = 0.5d0 * (delta_e - dsqrt(delta_e * delta_e + 4.d0 * i_H_psi_array(i) * i_H_psi_array(i)))
-      else
-        e_2_pert(i) = 0.5d0 * (delta_e + dsqrt(delta_e * delta_e + 4.d0 * i_H_psi_array(i) * i_H_psi_array(i)))
-      endif
-      if (dabs(i_H_psi_array(i)) > 1.d-6) then
-        c_pert(i) = e_2_pert(i)/i_H_psi_array(i)
-      else
-        c_pert(i) = 0.d0
-      endif
-      H_pert_diag(i) = h*c_pert(i)*c_pert(i)
-    else
-      e_2_pert(i) = 0.d0
-      c_pert(i) = 0.d0
-      H_pert_diag(i) = 0.d0
-    endif
-  enddo
-
-end

plugins/Perturbation/pert_sc2.irp.f

@@ -1,158 +1,3 @@
-
-subroutine pt2_epstein_nesbet_SC2_projected(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,N_st)
-  use bitmasks
-  implicit none
-  integer, intent(in)            :: Nint,ndet,N_st
-  integer(bit_kind), intent(in)  :: det_pert(Nint,2)
-  double precision , intent(out) :: c_pert(N_st),e_2_pert(N_st),H_pert_diag(N_st)
-  double precision               :: i_H_psi_array(N_st)
-  integer                        :: idx_repeat(0:ndet)
-  
-  BEGIN_DOC
-  ! compute the Epstein-Nesbet perturbative first order coefficient and second order energetic contribution
-  !
-  ! for the various N_st states, 
-  ! 
-  ! but  with the correction in the denominator 
-  ! 
-  ! comming from the interaction of that determinant with all the others determinants 
-  ! 
-  ! that can be repeated by repeating all the double excitations
-  !
-  ! : you repeat all the correlation energy already taken into account in CI_electronic_energy(1)
-  ! 
-  ! that could be repeated to this determinant.
-  !
-  ! In addition, for the perturbative energetic contribution you have the standard second order
-  !
-  ! e_2_pert = <psi_i|H|det_pert>^2/(Delta_E)
-  !
-  ! and also the purely projected contribution 
-  !
-  ! H_pert_diag = <HF|H|det_pert> c_pert
-  END_DOC
-  
-  integer                        :: i,j,degree,l
-  double precision               :: diag_H_mat_elem,accu_e_corr,hij,h0j,h,delta_E
-  double precision               :: repeat_all_e_corr,accu_e_corr_tmp,e_2_pert_fonda
-
-  ASSERT (Nint == N_int)
-  ASSERT (Nint > 0)
-
-  call i_H_psi_SC2(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array,idx_repeat)
-  accu_e_corr = 0.d0
-  !$IVDEP
-  do i = 1, idx_repeat(0)
-   accu_e_corr = accu_e_corr + E_corr_per_selectors(idx_repeat(i))
-  enddo
-  h =  diag_H_mat_elem(det_pert,Nint) + accu_e_corr
-  delta_E = 1.d0/(CI_SC2_electronic_energy(1) - h)
-
-
-  c_pert(1) = i_H_psi_array(1) /(CI_SC2_electronic_energy(1) - h)
-  e_2_pert(1) = i_H_psi_array(1) * c_pert(1)
-
-  do i =2,N_st
-    H_pert_diag(i) = h
-    if  (dabs(CI_SC2_electronic_energy(i) - h) > 1.d-6) then
-      c_pert(i) = i_H_psi_array(i) / (-dabs(CI_SC2_electronic_energy(i) - h))
-      e_2_pert(i) = (c_pert(i) * i_H_psi_array(i))
-    else
-      c_pert(i) = i_H_psi_array(i)
-      e_2_pert(i) = -dabs(i_H_psi_array(i))
-    endif
-  enddo
-
-  degree = popcnt(xor( ref_bitmask(1,1), det_pert(1,1))) +                      &
-           popcnt(xor( ref_bitmask(1,2), det_pert(1,2)))
-  !DEC$ NOUNROLL
-  do l=2,Nint
-    degree = degree+ popcnt(xor( ref_bitmask(l,1), det_pert(l,1))) +            &
-                     popcnt(xor( ref_bitmask(l,2), det_pert(l,2)))
-  enddo
-  if(degree==4)then
-  ! <psi|delta_H|psi>
-   e_2_pert_fonda = e_2_pert(1) 
-   H_pert_diag(1) = e_2_pert(1) * c_pert(1) * c_pert(1)
-   do i = 1, N_st
-    do j = 1, idx_repeat(0)
-     e_2_pert(i) += e_2_pert_fonda * psi_selectors_coef(idx_repeat(j),i) * psi_selectors_coef(idx_repeat(j),i)
-    enddo
-   enddo
-  endif
-  
-end
-
-
-subroutine pt2_epstein_nesbet_SC2_no_projected(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,N_st)
-  use bitmasks
-  implicit none
-  integer, intent(in)            :: Nint,ndet,N_st
-  integer(bit_kind), intent(in)  :: det_pert(Nint,2)
-  double precision , intent(out) :: c_pert(N_st),e_2_pert(N_st),H_pert_diag(N_st)
-  double precision               :: i_H_psi_array(N_st)
-  integer                        :: idx_repeat(0:ndet)
-  
-  BEGIN_DOC
-  ! compute the Epstein-Nesbet perturbative first order coefficient and second order energetic contribution
-  !
-  ! for the various N_st states, 
-  ! 
-  ! but  with the correction in the denominator 
-  ! 
-  ! comming from the interaction of that determinant with all the others determinants 
-  ! 
-  ! that can be repeated by repeating all the double excitations
-  !
-  ! : you repeat all the correlation energy already taken into account in CI_electronic_energy(1)
-  ! 
-  ! that could be repeated to this determinant.
-  !
-  ! In addition, for the perturbative energetic contribution you have the standard second order
-  !
-  ! e_2_pert = <psi_i|H|det_pert>^2/(Delta_E)
-  !
-  ! and also the purely projected contribution 
-  !
-  ! H_pert_diag = <HF|H|det_pert> c_pert
-  END_DOC
-  
-  integer                        :: i,j,degree,l
-  double precision               :: diag_H_mat_elem,accu_e_corr,hij,h0j,h,delta_E
-  double precision               :: repeat_all_e_corr,accu_e_corr_tmp,e_2_pert_fonda
-
-  ASSERT (Nint == N_int)
-  ASSERT (Nint > 0)
-
-  call i_H_psi_SC2(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array,idx_repeat)
-  accu_e_corr = 0.d0
-  !$IVDEP
-  do i = 1, idx_repeat(0)
-   accu_e_corr = accu_e_corr + E_corr_per_selectors(idx_repeat(i))
-  enddo
-  h =  diag_H_mat_elem(det_pert,Nint) + accu_e_corr
-  delta_E = 1.d0/(CI_SC2_electronic_energy(1) - h)
-
-
-  c_pert(1) = i_H_psi_array(1) /(CI_SC2_electronic_energy(1) - h)
-  e_2_pert(1) = i_H_psi_array(1) * c_pert(1)
-
-  do i =2,N_st
-    H_pert_diag(i) = h
-    if  (dabs(CI_SC2_electronic_energy(i) - h) > 1.d-6) then
-      c_pert(i) = i_H_psi_array(i) / (-dabs(CI_SC2_electronic_energy(i) - h))
-      e_2_pert(i) = (c_pert(i) * i_H_psi_array(i))
-    else
-      c_pert(i) = i_H_psi_array(i)
-      e_2_pert(i) = -dabs(i_H_psi_array(i))
-    endif
-  enddo
-end
-
-
-
-
-
 double precision function repeat_all_e_corr(key_in)
  implicit none
  integer(bit_kind), intent(in)  :: key_in(N_int,2)
@@ -182,47 +27,3 @@ double precision function repeat_all_e_corr(key_in)
 
 end
 
-
-subroutine pt2_epstein_nesbet_sc2(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,N_st)
-  use bitmasks
-  implicit none
-  integer, intent(in)            :: Nint,ndet,N_st
-  integer(bit_kind), intent(in)  :: det_pert(Nint,2)
-  double precision , intent(out) :: c_pert(N_st),e_2_pert(N_st),H_pert_diag(N_st)
-  double precision               :: i_H_psi_array(N_st)
-  
-  BEGIN_DOC
-  ! compute the standard Epstein-Nesbet perturbative first order coefficient and second order energetic contribution
-  !
-  ! for the various N_st states, but with the CISD_SC2 energies and coefficients
-  !
-  ! c_pert(i) = <psi(i)|H|det_pert>/( E(i) - <det_pert|H|det_pert> )
-  !
-  ! e_2_pert(i) = <psi(i)|H|det_pert>^2/( E(i) - <det_pert|H|det_pert> )
-  !
-  END_DOC
-  
-  integer                        :: i,j
-  double precision               :: diag_H_mat_elem, h
-  PROVIDE  selection_criterion
-
-  ASSERT (Nint == N_int)
-  ASSERT (Nint > 0)
-  call i_H_psi(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array)
-  h = diag_H_mat_elem(det_pert,Nint)
-  do i =1,N_st
-    if(CI_SC2_electronic_energy(i)>h.and.CI_SC2_electronic_energy(i).ne.0.d0)then
-      c_pert(i) = -1.d0
-      e_2_pert(i) = selection_criterion*selection_criterion_factor*2.d0
-    else if  (dabs(CI_SC2_electronic_energy(i) - h) > 1.d-6) then
-        c_pert(i) = i_H_psi_array(i) / (CI_SC2_electronic_energy(i) - h)
-        H_pert_diag(i) = h*c_pert(i)*c_pert(i)
-        e_2_pert(i) = c_pert(i) * i_H_psi_array(i)
-    else
-      c_pert(i) = -1.d0
-      e_2_pert(i) = -dabs(i_H_psi_array(i))
-      H_pert_diag(i) = h
-    endif
-  enddo
-  
-end

plugins/Perturbation/pert_single.irp.f

@@ -1,4 +1,4 @@
-subroutine pt2_h_core(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,N_st)
+subroutine pt2_h_core(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,N_st,minilist,idx_minilist,N_minilist)
   use bitmasks
   implicit none
   integer, intent(in)            :: Nint,ndet,N_st
@@ -6,6 +6,10 @@ subroutine pt2_h_core(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,N_st)
   double precision , intent(out) :: c_pert(N_st),e_2_pert(N_st),H_pert_diag(N_st)
   double precision               :: i_H_psi_array(N_st)
   
+  integer, intent(in)            :: N_minilist
+  integer, intent(in)            :: idx_minilist(0:N_det_selectors)
+  integer(bit_kind), intent(in)  :: minilist(Nint,2,N_det_selectors)
+  
   BEGIN_DOC
   ! compute the standard Epstein-Nesbet perturbative first order coefficient and second order energetic contribution
   !

plugins/Perturbation/perturbation.template.f

@@ -2,7 +2,7 @@ BEGIN_SHELL [ /usr/bin/env python ]
 import perturbation
 END_SHELL
 
-subroutine perturb_buffer_$PERT(i_generator,buffer,buffer_size,e_2_pert_buffer,coef_pert_buffer,sum_e_2_pert,sum_norm_pert,sum_H_pert_diag,N_st,Nint)
+subroutine perturb_buffer_$PERT(i_generator,buffer,buffer_size,e_2_pert_buffer,coef_pert_buffer,sum_e_2_pert,sum_norm_pert,sum_H_pert_diag,N_st,Nint,key_mask,fock_diag_tmp)
   implicit none
   BEGIN_DOC
   !  Applly pertubration ``$PERT`` to the buffer of determinants generated in the H_apply
@@ -11,23 +11,47 @@ subroutine perturb_buffer_$PERT(i_generator,buffer,buffer_size,e_2_pert_buffer,c
   
   integer, intent(in)            :: Nint, N_st, buffer_size, i_generator
   integer(bit_kind), intent(in)  :: buffer(Nint,2,buffer_size)
+  integer(bit_kind),intent(in)    :: key_mask(Nint,2)
+  double precision, intent(in)    :: fock_diag_tmp(2,0:mo_tot_num)
   double precision, intent(inout) :: sum_norm_pert(N_st),sum_e_2_pert(N_st)
   double precision, intent(inout) :: coef_pert_buffer(N_st,buffer_size),e_2_pert_buffer(N_st,buffer_size),sum_H_pert_diag(N_st)
   double precision               :: c_pert(N_st), e_2_pert(N_st),  H_pert_diag(N_st)
-  integer                        :: i,k, c_ref
+  integer                        :: i,k, c_ref, ni, ex
   integer, external              :: connected_to_ref
   logical, external              :: is_in_wavefunction
   
+  integer(bit_kind), allocatable :: minilist(:,:,:)
+  integer, allocatable           :: idx_minilist(:)
+  integer                        :: N_minilist
+  
+  integer(bit_kind), allocatable :: minilist_gen(:,:,:)
+  integer :: N_minilist_gen
+  logical :: fullMatch
+  logical, external :: is_connected_to
+
+  allocate( minilist(Nint,2,N_det_selectors),                        &
+      minilist_gen(Nint,2,N_det_generators),                         &
+      idx_minilist(N_det_selectors) )
+  
+
   ASSERT (Nint > 0)
   ASSERT (Nint == N_int)
   ASSERT (buffer_size >= 0)
   ASSERT (minval(sum_norm_pert) >= 0.d0)
   ASSERT (N_st > 0)
+  
+  call create_minilist(key_mask, psi_selectors, miniList, idx_miniList, N_det_selectors, N_minilist, Nint)
+  call create_minilist_find_previous(key_mask, psi_det_generators, miniList_gen, i_generator-1, N_minilist_gen, fullMatch, Nint)
+  
+  if(fullMatch) then
+    deallocate( minilist, minilist_gen, idx_minilist )
+    return
+  end if
+  
+  
   do i=1,buffer_size
     
-    c_ref = connected_to_ref(buffer(1,1,i),psi_det_generators,Nint,i_generator,N_det_generators)
-
-    if (c_ref /= 0) then
+    if(is_connected_to(buffer(1,1,i), miniList_gen, Nint, N_minilist_gen)) then
       cycle
     end if
     
@@ -35,25 +59,25 @@ subroutine perturb_buffer_$PERT(i_generator,buffer,buffer_size,e_2_pert_buffer,c
       cycle
     endif
     
-    integer :: degree
-    call get_excitation_degree(HF_bitmask,buffer(1,1,i),degree,N_int)
-    call pt2_$PERT(buffer(1,1,i),         &
-        c_pert,e_2_pert,H_pert_diag,Nint,N_det_selectors,n_st)
+    call pt2_$PERT(psi_det_generators(1,1,i_generator),buffer(1,1,i), fock_diag_tmp,        &
+         c_pert,e_2_pert,H_pert_diag,Nint,N_minilist,n_st,minilist,idx_minilist,N_minilist) 
 
     do k = 1,N_st
       e_2_pert_buffer(k,i)   = e_2_pert(k)
       coef_pert_buffer(k,i)  = c_pert(k)
-      sum_norm_pert(k) += c_pert(k) * c_pert(k)
-      sum_e_2_pert(k) += e_2_pert(k)
-      sum_H_pert_diag(k) +=  H_pert_diag(k)
+      sum_norm_pert(k)       = sum_norm_pert(k)   + c_pert(k) * c_pert(k)
+      sum_e_2_pert(k)        = sum_e_2_pert(k)    + e_2_pert(k)
+      sum_H_pert_diag(k)     = sum_H_pert_diag(k) + H_pert_diag(k)
     enddo
     
   enddo 
+  deallocate( minilist, minilist_gen, idx_minilist )
   
 end
 
 
-subroutine perturb_buffer_by_mono_$PERT(i_generator,buffer,buffer_size,e_2_pert_buffer,coef_pert_buffer,sum_e_2_pert,sum_norm_pert,sum_H_pert_diag,N_st,Nint)
+
+subroutine perturb_buffer_by_mono_$PERT(i_generator,buffer,buffer_size,e_2_pert_buffer,coef_pert_buffer,sum_e_2_pert,sum_norm_pert,sum_H_pert_diag,N_st,Nint,key_mask,fock_diag_tmp)
   implicit none
   BEGIN_DOC
   !  Applly pertubration ``$PERT`` to the buffer of determinants generated in the H_apply
@@ -62,18 +86,44 @@ subroutine perturb_buffer_by_mono_$PERT(i_generator,buffer,buffer_size,e_2_pert_
   
   integer, intent(in)            :: Nint, N_st, buffer_size, i_generator
   integer(bit_kind), intent(in)  :: buffer(Nint,2,buffer_size)
+  integer(bit_kind),intent(in)    :: key_mask(Nint,2)
+  double precision, intent(in)    :: fock_diag_tmp(2,0:mo_tot_num)
   double precision, intent(inout) :: sum_norm_pert(N_st),sum_e_2_pert(N_st)
   double precision, intent(inout) :: coef_pert_buffer(N_st,buffer_size),e_2_pert_buffer(N_st,buffer_size),sum_H_pert_diag(N_st)
   double precision               :: c_pert(N_st), e_2_pert(N_st),  H_pert_diag(N_st)
-  integer                        :: i,k, c_ref
+  integer                        :: i,k, c_ref, ni, ex
   integer, external              :: connected_to_ref_by_mono
   logical, external              :: is_in_wavefunction
   
+  integer(bit_kind), allocatable :: minilist(:,:,:)
+  integer, allocatable           :: idx_minilist(:)
+  integer                        :: N_minilist
+  
+  integer(bit_kind), allocatable :: minilist_gen(:,:,:)
+  integer :: N_minilist_gen
+  logical :: fullMatch
+  logical, external :: is_connected_to
+
+  allocate( minilist(Nint,2,N_det_selectors),                        &
+      minilist_gen(Nint,2,N_det_generators),                         &
+      idx_minilist(N_det_selectors) )
+  
+
   ASSERT (Nint > 0)
   ASSERT (Nint == N_int)
   ASSERT (buffer_size >= 0)
   ASSERT (minval(sum_norm_pert) >= 0.d0)
   ASSERT (N_st > 0)
+  
+  call create_minilist(key_mask, psi_selectors, miniList, idx_miniList, N_det_selectors, N_minilist, Nint)
+  call create_minilist_find_previous(key_mask, psi_det_generators, miniList_gen, i_generator-1, N_minilist_gen, fullMatch, Nint)
+  
+  if(fullMatch) then
+    deallocate( minilist, minilist_gen, idx_minilist )
+    return
+  end if
+  
+  
   do i=1,buffer_size
     
     c_ref = connected_to_ref_by_mono(buffer(1,1,i),psi_det_generators,Nint,i_generator,N_det)
@@ -86,19 +136,19 @@ subroutine perturb_buffer_by_mono_$PERT(i_generator,buffer,buffer_size,e_2_pert_
       cycle
     endif
     
-   integer :: degree
-    call pt2_$PERT(buffer(1,1,i),         &
-        c_pert,e_2_pert,H_pert_diag,Nint,N_det_selectors,n_st)
+    call pt2_$PERT(psi_det_generators(1,1,i_generator),buffer(1,1,i), fock_diag_tmp,        &
+         c_pert,e_2_pert,H_pert_diag,Nint,N_minilist,n_st,minilist,idx_minilist,N_minilist) 
 
     do k = 1,N_st
       e_2_pert_buffer(k,i)   = e_2_pert(k)
       coef_pert_buffer(k,i)  = c_pert(k)
-      sum_norm_pert(k) += c_pert(k) * c_pert(k)
-      sum_e_2_pert(k) += e_2_pert(k)
-      sum_H_pert_diag(k) +=  H_pert_diag(k)
+      sum_norm_pert(k)       = sum_norm_pert(k)   + c_pert(k) * c_pert(k)
+      sum_e_2_pert(k)        = sum_e_2_pert(k)    + e_2_pert(k)
+      sum_H_pert_diag(k)     = sum_H_pert_diag(k) + H_pert_diag(k)
     enddo
     
   enddo 
+  deallocate( minilist, minilist_gen, idx_minilist )
   
 end
 

plugins/Perturbation/pt2_equations.irp.f

@@ -0,0 +1,367 @@
+BEGIN_TEMPLATE
+
+subroutine pt2_epstein_nesbet ($arguments)
+  use bitmasks
+  implicit none
+  $declarations
+  
+  BEGIN_DOC
+  ! compute the standard Epstein-Nesbet perturbative first order coefficient and second order energetic contribution
+  !
+  ! for the various N_st states.
+  !
+  ! c_pert(i) = <psi(i)|H|det_pert>/( E(i) - <det_pert|H|det_pert> )
+  !
+  ! e_2_pert(i) = <psi(i)|H|det_pert>^2/( E(i) - <det_pert|H|det_pert> )
+  !
+  END_DOC
+  
+  integer                        :: i,j
+  double precision               :: diag_H_mat_elem_fock, h
+  double precision               :: i_H_psi_array(N_st)
+  PROVIDE  selection_criterion
+
+  ASSERT (Nint == N_int)
+  ASSERT (Nint > 0)
+  !call i_H_psi(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array)
+  call i_H_psi_minilist(det_pert,minilist,idx_minilist,N_minilist,psi_selectors_coef,Nint,N_minilist,psi_selectors_size,N_st,i_H_psi_array)
+  
+  
+  h = diag_H_mat_elem_fock(det_ref,det_pert,fock_diag_tmp,Nint)
+  do i =1,N_st
+    if(CI_electronic_energy(i)>h.and.CI_electronic_energy(i).ne.0.d0)then
+      c_pert(i) = -1.d0
+      e_2_pert(i) = selection_criterion*selection_criterion_factor*2.d0
+    else if  (dabs(CI_electronic_energy(i) - h) > 1.d-6) then
+        c_pert(i) = i_H_psi_array(i) / (CI_electronic_energy(i) - h)
+        H_pert_diag(i) = h*c_pert(i)*c_pert(i)
+        e_2_pert(i) = c_pert(i) * i_H_psi_array(i)
+    else
+      c_pert(i) = -1.d0
+      e_2_pert(i) = -dabs(i_H_psi_array(i))
+      H_pert_diag(i) = h
+    endif
+  enddo
+  
+end
+
+subroutine pt2_epstein_nesbet_2x2 ($arguments)
+  use bitmasks
+  implicit none
+  $declarations
+  
+  BEGIN_DOC
+  ! compute the Epstein-Nesbet 2x2 diagonalization coefficient and energetic contribution
+  !
+  ! for the various N_st states.
+  !
+  ! e_2_pert(i) = 0.5 * (( <det_pert|H|det_pert> -  E(i) )  - sqrt( ( <det_pert|H|det_pert> -  E(i)) ^2 + 4 <psi(i)|H|det_pert>^2  )
+  !
+  ! c_pert(i) = e_2_pert(i)/ <psi(i)|H|det_pert>
+  !
+  END_DOC
+  
+  integer                        :: i,j
+  double precision               :: diag_H_mat_elem_fock,delta_e, h
+  double precision               :: i_H_psi_array(N_st)
+  ASSERT (Nint == N_int)
+  ASSERT (Nint > 0)
+  PROVIDE CI_electronic_energy
+
+  !call i_H_psi(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array)
+  call i_H_psi_minilist(det_pert,minilist,idx_minilist,N_minilist,psi_selectors_coef,Nint,N_minilist,psi_selectors_size,N_st,i_H_psi_array)
+  
+  h = diag_H_mat_elem_fock(det_ref,det_pert,fock_diag_tmp,Nint)
+  do i =1,N_st
+    if (i_H_psi_array(i) /= 0.d0) then
+      delta_e = h - CI_electronic_energy(i)
+      if (delta_e > 0.d0) then
+        e_2_pert(i) = 0.5d0 * (delta_e - dsqrt(delta_e * delta_e + 4.d0 * i_H_psi_array(i) * i_H_psi_array(i)))
+      else
+        e_2_pert(i) = 0.5d0 * (delta_e + dsqrt(delta_e * delta_e + 4.d0 * i_H_psi_array(i) * i_H_psi_array(i)))
+      endif
+      if (dabs(i_H_psi_array(i)) > 1.d-6) then
+        c_pert(i) = e_2_pert(i)/i_H_psi_array(i)
+      else
+        c_pert(i) = 0.d0
+      endif
+      H_pert_diag(i) = h*c_pert(i)*c_pert(i)
+    else
+      e_2_pert(i) = 0.d0
+      c_pert(i) = 0.d0
+      H_pert_diag(i) = 0.d0
+    endif
+  enddo
+
+end
+
+subroutine pt2_moller_plesset ($arguments)
+  use bitmasks
+  implicit none
+  $declarations
+  
+  BEGIN_DOC
+  ! compute the standard Moller-Plesset perturbative first order coefficient and second order energetic contribution
+  !
+  ! for the various n_st states.
+  !
+  ! c_pert(i) = <psi(i)|H|det_pert>/(difference of orbital energies) 
+  !
+  ! e_2_pert(i) = <psi(i)|H|det_pert>^2/(difference of orbital energies) 
+  !
+  END_DOC
+  
+  integer                        :: i,j
+  double precision               :: diag_H_mat_elem_fock
+  integer                        :: exc(0:2,2,2)
+  integer                        :: degree
+  double precision               :: phase,delta_e,h
+  double precision               :: i_H_psi_array(N_st)
+  integer                        :: h1,h2,p1,p2,s1,s2
+  ASSERT (Nint == N_int)
+  ASSERT (Nint > 0)
+  call get_excitation(ref_bitmask,det_pert,exc,degree,phase,Nint)
+  if (degree == 2) then
+    call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
+    delta_e = Fock_matrix_diag_mo(h1) + Fock_matrix_diag_mo(h2) - &
+            (Fock_matrix_diag_mo(p1) + Fock_matrix_diag_mo(p2))
+    delta_e = 1.d0/delta_e
+  else if (degree == 1) then
+    call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
+    delta_e = Fock_matrix_diag_mo(h1) - Fock_matrix_diag_mo(p1) 
+    delta_e = 1.d0/delta_e
+  else
+    delta_e = 0.d0
+  endif
+
+  call i_H_psi(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det,psi_selectors_size,n_st,i_H_psi_array)
+  h = diag_H_mat_elem_fock(det_ref,det_pert,fock_diag_tmp,Nint)
+  do i =1,n_st
+    H_pert_diag(i) = h
+    c_pert(i) = i_H_psi_array(i) *delta_e
+    e_2_pert(i) = c_pert(i) * i_H_psi_array(i)
+  enddo
+  
+end
+
+
+subroutine pt2_epstein_nesbet_SC2_projected ($arguments)
+  use bitmasks
+  implicit none
+  $declarations
+  BEGIN_DOC
+  ! compute the Epstein-Nesbet perturbative first order coefficient and second order energetic contribution
+  !
+  ! for the various N_st states, 
+  ! 
+  ! but  with the correction in the denominator 
+  ! 
+  ! comming from the interaction of that determinant with all the others determinants 
+  ! 
+  ! that can be repeated by repeating all the double excitations
+  !
+  ! : you repeat all the correlation energy already taken into account in CI_electronic_energy(1)
+  ! 
+  ! that could be repeated to this determinant.
+  !
+  ! In addition, for the perturbative energetic contribution you have the standard second order
+  !
+  ! e_2_pert = <psi_i|H|det_pert>^2/(Delta_E)
+  !
+  ! and also the purely projected contribution 
+  !
+  ! H_pert_diag = <HF|H|det_pert> c_pert
+  END_DOC
+  
+  double precision               :: i_H_psi_array(N_st)
+  integer                        :: idx_repeat(0:ndet)
+  integer                        :: i,j,degree,l
+  double precision               :: diag_H_mat_elem_fock,accu_e_corr,hij,h0j,h,delta_E
+  double precision               :: repeat_all_e_corr,accu_e_corr_tmp,e_2_pert_fonda
+
+  ASSERT (Nint == N_int)
+  ASSERT (Nint > 0)
+
+  call i_H_psi_SC2(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array,idx_repeat)
+  accu_e_corr = 0.d0
+  !$IVDEP
+  do i = 1, idx_repeat(0)
+   accu_e_corr = accu_e_corr + E_corr_per_selectors(idx_repeat(i))
+  enddo
+  h = diag_H_mat_elem_fock(det_ref,det_pert,fock_diag_tmp,Nint)
+  h = h + accu_e_corr
+  delta_E = 1.d0/(CI_SC2_electronic_energy(1) - h)
+
+
+  c_pert(1) = i_H_psi_array(1) /(CI_SC2_electronic_energy(1) - h)
+  e_2_pert(1) = i_H_psi_array(1) * c_pert(1)
+
+  do i =2,N_st
+    H_pert_diag(i) = h
+    if  (dabs(CI_SC2_electronic_energy(i) - h) > 1.d-6) then
+      c_pert(i) = i_H_psi_array(i) / (-dabs(CI_SC2_electronic_energy(i) - h))
+      e_2_pert(i) = (c_pert(i) * i_H_psi_array(i))
+    else
+      c_pert(i) = i_H_psi_array(i)
+      e_2_pert(i) = -dabs(i_H_psi_array(i))
+    endif
+  enddo
+
+  degree = popcnt(xor( ref_bitmask(1,1), det_pert(1,1))) +                      &
+           popcnt(xor( ref_bitmask(1,2), det_pert(1,2)))
+  !DEC$ NOUNROLL
+  do l=2,Nint
+    degree = degree+ popcnt(xor( ref_bitmask(l,1), det_pert(l,1))) +            &
+                     popcnt(xor( ref_bitmask(l,2), det_pert(l,2)))
+  enddo
+  if(degree==4)then
+  ! <psi|delta_H|psi>
+   e_2_pert_fonda = e_2_pert(1) 
+   H_pert_diag(1) = e_2_pert(1) * c_pert(1) * c_pert(1)
+   do i = 1, N_st
+    do j = 1, idx_repeat(0)
+     e_2_pert(i) += e_2_pert_fonda * psi_selectors_coef(idx_repeat(j),i) * psi_selectors_coef(idx_repeat(j),i)
+    enddo
+   enddo
+  endif
+  
+end
+
+
+subroutine pt2_epstein_nesbet_SC2_no_projected ($arguments)
+  use bitmasks
+  implicit none
+  $declarations
+  BEGIN_DOC
+  ! compute the Epstein-Nesbet perturbative first order coefficient and second order energetic contribution
+  !
+  ! for the various N_st states, 
+  ! 
+  ! but  with the correction in the denominator 
+  ! 
+  ! comming from the interaction of that determinant with all the others determinants 
+  ! 
+  ! that can be repeated by repeating all the double excitations
+  !
+  ! : you repeat all the correlation energy already taken into account in CI_electronic_energy(1)
+  ! 
+  ! that could be repeated to this determinant.
+  !
+  ! In addition, for the perturbative energetic contribution you have the standard second order
+  !
+  ! e_2_pert = <psi_i|H|det_pert>^2/(Delta_E)
+  !
+  ! and also the purely projected contribution 
+  !
+  ! H_pert_diag = <HF|H|det_pert> c_pert
+  END_DOC
+  
+  double precision               :: i_H_psi_array(N_st)
+  integer                        :: idx_repeat(0:ndet)
+  integer                        :: i,j,degree,l
+  double precision               :: diag_H_mat_elem_fock,accu_e_corr,hij,h0j,h,delta_E
+  double precision               :: repeat_all_e_corr,accu_e_corr_tmp,e_2_pert_fonda
+
+  ASSERT (Nint == N_int)
+  ASSERT (Nint > 0)
+
+  call i_H_psi_SC2(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array,idx_repeat)
+  accu_e_corr = 0.d0
+  !$IVDEP
+  do i = 1, idx_repeat(0)
+   accu_e_corr = accu_e_corr + E_corr_per_selectors(idx_repeat(i))
+  enddo
+  h = diag_H_mat_elem_fock(det_ref,det_pert,fock_diag_tmp,Nint)
+  h = h + accu_e_corr
+  delta_E = 1.d0/(CI_SC2_electronic_energy(1) - h)
+
+
+  c_pert(1) = i_H_psi_array(1) /(CI_SC2_electronic_energy(1) - h)
+  e_2_pert(1) = i_H_psi_array(1) * c_pert(1)
+
+  do i =2,N_st
+    H_pert_diag(i) = h
+    if  (dabs(CI_SC2_electronic_energy(i) - h) > 1.d-6) then
+      c_pert(i) = i_H_psi_array(i) / (-dabs(CI_SC2_electronic_energy(i) - h))
+      e_2_pert(i) = (c_pert(i) * i_H_psi_array(i))
+    else
+      c_pert(i) = i_H_psi_array(i)
+      e_2_pert(i) = -dabs(i_H_psi_array(i))
+    endif
+  enddo
+end
+
+
+
+
+
+subroutine pt2_epstein_nesbet_sc2 ($arguments)
+  use bitmasks
+  implicit none
+  $declarations
+  BEGIN_DOC
+  ! compute the standard Epstein-Nesbet perturbative first order coefficient and second order energetic contribution
+  !
+  ! for the various N_st states, but with the CISD_SC2 energies and coefficients
+  !
+  ! c_pert(i) = <psi(i)|H|det_pert>/( E(i) - <det_pert|H|det_pert> )
+  !
+  ! e_2_pert(i) = <psi(i)|H|det_pert>^2/( E(i) - <det_pert|H|det_pert> )
+  !
+  END_DOC
+  
+  integer                        :: i,j
+  double precision               :: i_H_psi_array(N_st)
+  double precision               :: diag_H_mat_elem_fock, h
+  PROVIDE  selection_criterion
+
+  ASSERT (Nint == N_int)
+  ASSERT (Nint > 0)
+  !call i_H_psi(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array)
+  call i_H_psi_minilist(det_pert,minilist,idx_minilist,N_minilist,psi_selectors_coef,Nint,N_minilist,psi_selectors_size,N_st,i_H_psi_array)
+
+  
+  h = diag_H_mat_elem_fock(det_ref,det_pert,fock_diag_tmp,Nint)
+  do i =1,N_st
+    if(CI_SC2_electronic_energy(i)>h.and.CI_SC2_electronic_energy(i).ne.0.d0)then
+      c_pert(i) = -1.d0
+      e_2_pert(i) = selection_criterion*selection_criterion_factor*2.d0
+    else if  (dabs(CI_SC2_electronic_energy(i) - h) > 1.d-6) then
+        c_pert(i) = i_H_psi_array(i) / (CI_SC2_electronic_energy(i) - h)
+        H_pert_diag(i) = h*c_pert(i)*c_pert(i)
+        e_2_pert(i) = c_pert(i) * i_H_psi_array(i)
+    else
+      c_pert(i) = -1.d0
+      e_2_pert(i) = -dabs(i_H_psi_array(i))
+      H_pert_diag(i) = h
+    endif
+  enddo
+  
+end
+
+
+
+SUBST [ arguments, declarations ]
+
+det_ref,det_pert,fock_diag_tmp,c_pert,e_2_pert,H_pert_diag,Nint,ndet,N_st,minilist,idx_minilist,N_minilist ;
+
+    integer, intent(in)             :: Nint
+    integer, intent(in)             :: ndet
+    integer, intent(in)             :: N_st
+    integer, intent(in)             :: N_minilist
+    integer(bit_kind), intent(in)   :: det_ref (Nint,2)
+    integer(bit_kind), intent(in)   :: det_pert(Nint,2)
+    double precision , intent(in)   :: fock_diag_tmp(2,mo_tot_num+1)
+    double precision , intent(out)  :: c_pert(N_st)
+    double precision , intent(out)  :: e_2_pert(N_st)
+    double precision, intent(out)   :: H_pert_diag(N_st)
+    integer, intent(in)             :: idx_minilist(0:N_det_selectors)
+    integer(bit_kind), intent(in)   :: minilist(Nint,2,N_det_selectors)
+;;
+
+
+END_TEMPLATE
+
+! Note : If the arguments are changed here, they should also be changed accordingly in
+! the perturbation.template.f file.
+

plugins/QmcChem/.gitignore

@@ -21,3 +21,4 @@ irpf90.make
 irpf90_entities
 save_for_qmcchem
 tags
+target_pt2_qmc

plugins/QmcChem/README.rst

@@ -66,6 +66,14 @@ Documentation
   title="f(|r-r_A|) = \int Y_{lm}^{C} (|r-r_C|, \Omega_C) \chi_i^{A} (r-r_A) d\Omega_C" />
 
 
+`compute_energy <http://github.com/LCPQ/quantum_package/tree/master/plugins/QmcChem/target_pt2_qmc.irp.f#L80>`_
+  Compute an energy when a threshold is applied
+
+
+`e_curve <http://github.com/LCPQ/quantum_package/tree/master/plugins/QmcChem/target_pt2_qmc.irp.f#L1>`_
+  Undocumented
+
+
 `mo_pseudo_grid <http://github.com/LCPQ/quantum_package/tree/master/plugins/QmcChem/pot_ao_pseudo_ints.irp.f#L56>`_
   Grid points for f(|r-r_A|) = \int Y_{lm}^{C} (|r-r_C|, \Omega_C) \phi_i^{A} (r-r_A) d\Omega_C
   .br

plugins/QmcChem/save_for_qmcchem.irp.f

@@ -1,6 +1,7 @@
 program save_for_qmc
  read_wf = .True.
  TOUCH read_wf
+ print *,  "N_det = ", N_det
  call write_spindeterminants 
  if (do_pseudo) then
    call write_pseudopotential

plugins/QmcChem/target_pt2_qmc.irp.f

@@ -0,0 +1,121 @@
+program e_curve
+  use bitmasks
+  implicit none
+  integer :: i,j,k, nab, m, l, n_up, n_dn, n
+  double precision :: norm, E, hij, num, ci, cj
+  integer, allocatable :: iorder(:)
+  double precision , allocatable :: norm_sort(:), psi_bilinear_matrix_values_save(:)
+  nab = n_det_alpha_unique+n_det_beta_unique
+
+  allocate ( norm_sort(0:nab), iorder(0:nab), psi_bilinear_matrix_values_save(N_det) )
+
+
+  norm_sort(0) = 0.d0
+  iorder(0) = 0
+  do i=1,n_det_alpha_unique
+   norm_sort(i) = det_alpha_norm(i)
+   iorder(i) = i
+  enddo
+  
+  do i=1,n_det_beta_unique
+   norm_sort(i+n_det_alpha_unique) = det_beta_norm(i)
+   iorder(i+n_det_alpha_unique) = -i
+  enddo
+  
+  call dsort(norm_sort(1),iorder(1),nab)
+
+  if (.not.read_wf) then
+    stop 'Please set read_wf to true'
+  endif
+
+  psi_bilinear_matrix_values_save = psi_bilinear_matrix_values(:,1)
+  print *,  '=========================================================='
+  print '(A8,2X,A8,2X,A12,2X,A10,2X,A12)',  'Thresh.', 'Ndet', 'Cost', 'Norm', 'E'
+  print *,  '=========================================================='
+  integer(bit_kind), allocatable :: det_i(:,:), det_j(:,:)
+
+  double precision :: thresh, E_min, E_max, E_prev
+  thresh = 0.d0
+  call compute_energy(psi_bilinear_matrix_values_save,E_max,m,norm)
+  call i_h_j(psi_det_sorted(1,1,1), psi_det_sorted(1,1,1), N_int, E_min)
+  print *,  E_min, E_max
+
+  n_up = nab
+  n_dn = 0
+  do while (n_up > n_dn)
+    n = n_dn + (n_up-n_dn)/2
+    psi_bilinear_matrix_values_save = psi_bilinear_matrix_values(:,1)
+    do j=1,n
+      i = iorder(j)
+      if (i<0) then
+        do k=1,n_det
+          if (psi_bilinear_matrix_columns(k) == -i) then
+            psi_bilinear_matrix_values_save(k) = 0.d0
+          endif
+        enddo
+      else
+        do k=1,n_det
+          if (psi_bilinear_matrix_rows(k) == i) then
+            psi_bilinear_matrix_values_save(k) = 0.d0
+          endif
+        enddo
+      endif
+    enddo
+    call compute_energy(psi_bilinear_matrix_values_save,E,m,norm)
+    print '(E9.1,2X,I8,2X,F10.2,2X,F10.6,2X,F12.6)',  norm_sort(n), m, &
+      dble( elec_alpha_num**3 + elec_alpha_num**2 * m ) / &
+      dble( elec_alpha_num**3 + elec_alpha_num**2 * n ), norm, E
+    if (E < target_energy) then
+       n_dn = n+1
+    else
+       n_up = n
+    endif
+  enddo
+  print *,  '=========================================================='
+  print *,  norm_sort(n), target_energy
+
+  deallocate (iorder, norm_sort, psi_bilinear_matrix_values_save)
+end
+
+subroutine compute_energy(psi_bilinear_matrix_values_save, E, m, norm)
+  implicit none
+  BEGIN_DOC
+  ! Compute an energy when a threshold is applied
+  END_DOC
+  double precision, intent(in) :: psi_bilinear_matrix_values_save(n_det)
+  integer(bit_kind), allocatable :: det_i(:,:), det_j(:,:)
+  integer :: i,j, k, l, m
+  double precision :: num, norm, ci, cj, hij, E
+
+  
+    num = 0.d0
+    norm = 0.d0
+    m = 0
+    !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(k,l,det_i,det_j,ci,cj,hij) REDUCTION(+:norm,m,num)
+    allocate( det_i(N_int,2), det_j(N_int,2))
+    !$OMP DO schedule(guided)
+    do k=1,n_det
+      if (psi_bilinear_matrix_values_save(k) == 0.d0) then
+        cycle
+      endif
+      ci = psi_bilinear_matrix_values_save(k)
+      det_i(:,1) = psi_det_alpha_unique(:,psi_bilinear_matrix_rows(k))
+      det_i(:,2) = psi_det_beta_unique(:,psi_bilinear_matrix_columns(k))
+      do l=1,n_det
+        if (psi_bilinear_matrix_values_save(l) == 0.d0) then
+          cycle
+        endif
+        cj = psi_bilinear_matrix_values_save(l)
+        det_j(:,1) = psi_det_alpha_unique(:,psi_bilinear_matrix_rows(l))
+        det_j(:,2) = psi_det_beta_unique(:,psi_bilinear_matrix_columns(l))
+        call i_h_j(det_i, det_j, N_int, hij)
+        num = num + ci*cj*hij
+      enddo
+     norm = norm + ci*ci
+     m = m+1
+    enddo
+    !$OMP END DO
+    deallocate (det_i,det_j)
+    !$OMP END PARALLEL
+    E = num / norm + nuclear_repulsion
+end

plugins/loc_cele/README.rst

@@ -20,3 +20,25 @@ Needed Modules
 
 * `MO_Basis <http://github.com/LCPQ/quantum_package/tree/master/src/MO_Basis>`_
 
+Needed Modules
+==============
+.. Do not edit this section It was auto-generated
+.. by the `update_README.py` script.
+
+
+.. image:: tree_dependency.png
+
+* `MO_Basis <http://github.com/LCPQ/quantum_package/tree/master/src/MO_Basis>`_
+
+Documentation
+=============
+.. Do not edit this section It was auto-generated
+.. by the `update_README.py` script.
+
+
+`loc_rasorb <http://github.com/LCPQ/quantum_package/tree/master/plugins/loc_cele/loc_cele.irp.f#L1>`_
+  This program performs a localization of the active orbitals
+  of a CASSCF wavefunction, reading the orbitals from a RASORB
+  file of molcas.
+  id1=max is the number of MO in a given symmetry.
+

plugins/loc_cele/loc_cele.irp.f

@@ -280,7 +280,7 @@
 
       enddo !big loop over symmetry
 
-      10 format (4E18.12)
+      10 format (4E20.12)
 
 
 !  Now we copyt the newcmo into the mo_coef

scripts/compilation/qp_create_ninja.py

@@ -35,7 +35,9 @@ except ImportError:
 
 from qp_path import QP_ROOT, QP_SRC, QP_EZFIO
 
-EZFIO_LIB = join(QP_ROOT, "lib", "libezfio.a")
+LIB = "" # join(QP_ROOT, "lib", "rdtsc.o") 
+EZFIO_LIB = join(QP_ROOT, "lib", "libezfio_irp.a") 
+ZMQ_LIB   = join(QP_ROOT, "lib", "libf77zmq.a") + " "  + join(QP_ROOT, "lib", "libzmq.a") + " -lstdc++ -lrt"
 ROOT_BUILD_NINJA = join(QP_ROOT, "config", "build.ninja")
 
 header = r"""#
@@ -94,7 +96,7 @@ def ninja_create_env_variable(pwd_config_file):
         l_string.append(str_)
 
     lib_lapack = get_compilation_option(pwd_config_file, "LAPACK_LIB")
-    l_string.append("{0} = {1} {2}".format("LIB", lib_lapack, EZFIO_LIB))
+    l_string.append("LIB = {0} {1} {2} {3}".format(LIB, lib_lapack, EZFIO_LIB, ZMQ_LIB))
 
     l_string.append("")
 
@@ -260,7 +262,7 @@ def ninja_ezfio_rule():
     l_flag = ["export {0}='${0}'".format(flag)
               for flag in ["FC", "FCFLAGS", "IRPF90"]]
 
-    install_lib_ezfio = join(QP_ROOT, 'install', 'EZFIO', "lib", "libezfio.a")
+    install_lib_ezfio = join(QP_ROOT, 'install', 'EZFIO', "lib", "libezfio_irp.a")
     l_cmd = ["cd {0}".format(QP_EZFIO)] + l_flag
     l_cmd += ["rm -f make.config ; ninja && ln -sf {0} {1}".format(install_lib_ezfio, EZFIO_LIB)]
 
@@ -707,7 +709,7 @@ def ninja_dot_tree_rule():
     l_string = [
         "rule build_dot_tree", "   command = {0}".format(" ; ".join(l_cmd)),
         "   generator = 1",
-        "   description = Generate Png representtion of the Tree Dependencies of $module_rel",
+        "   description = Generating Png representation of the Tree Dependencies of $module_rel",
         ""
     ]
 

scripts/generate_h_apply.py

@@ -60,7 +60,7 @@ class H_apply(object):
     s["omp_master"]       = "!$OMP MASTER"
     s["omp_end_master"]   = "!$OMP END MASTER"
     s["omp_barrier"]      = "!$OMP BARRIER"
-    s["omp_do"]           = "!$OMP DO SCHEDULE (static)"
+    s["omp_do"]           = "!$OMP DO SCHEDULE (static,1)"
     s["omp_enddo"]        = "!$OMP ENDDO NOWAIT"
 
     d = { True : '.True.', False : '.False.'}
@@ -99,7 +99,7 @@ class H_apply(object):
   deallocate(H_jj,iorder)
     """
 
-    s["size_max"] = str(1024*128) 
+    s["size_max"] = "256"
     s["copy_buffer"] = """call copy_H_apply_buffer_to_wf
   if (s2_eig) then
     call make_s2_eigenfunction
@@ -131,10 +131,10 @@ class H_apply(object):
   def filter_vvvv_excitation(self):
     self["filter_vvvv_excitation"] = """
      key_union_hole_part = 0_bit_kind
-     call set_bite_to_integer(i_a,key_union_hole_part,N_int)
-     call set_bite_to_integer(j_a,key_union_hole_part,N_int)
-     call set_bite_to_integer(i_b,key_union_hole_part,N_int)
-     call set_bite_to_integer(j_b,key_union_hole_part,N_int)
+     call set_bit_to_integer(i_a,key_union_hole_part,N_int)
+     call set_bit_to_integer(j_a,key_union_hole_part,N_int)
+     call set_bit_to_integer(i_b,key_union_hole_part,N_int)
+     call set_bit_to_integer(j_b,key_union_hole_part,N_int)
      do jtest_vvvv = 1, N_int
       if(iand(key_union_hole_part(jtest_vvvv),virt_bitmask(jtest_vvvv,1).ne.key_union_hole_part(jtest_vvvv)))then
        b_cycle = .False.
@@ -157,7 +157,6 @@ class H_apply(object):
 
   def set_filter_2h_2p(self):
     self["filter2h2p"] = """
-!    ! DIR$ FORCEINLINE
      if (is_a_two_holes_two_particles(key)) cycle
     """
 
@@ -201,11 +200,11 @@ class H_apply(object):
       """
       self.data["size_max"] = "256" 
       self.data["initialization"] = """
-      PROVIDE CI_electronic_energy psi_selectors_coef psi_selectors E_corr_per_selectors psi_det_sorted_bit
+      PROVIDE psi_selectors_coef psi_selectors E_corr_per_selectors psi_det_sorted_bit
       """
       self.data["keys_work"] = """
       call perturb_buffer_%s(i_generator,keys_out,key_idx,e_2_pert_buffer,coef_pert_buffer,sum_e_2_pert, &
-       sum_norm_pert,sum_H_pert_diag,N_st,N_int)
+       sum_norm_pert,sum_H_pert_diag,N_st,N_int,key_mask,fock_diag_tmp)
       """%(pert,)
       self.data["finalization"] = """
       """
@@ -219,7 +218,7 @@ class H_apply(object):
   double precision, intent(inout):: norm_pert(N_st) 
   double precision, intent(inout):: H_pert_diag(N_st)
   double precision               :: delta_pt2(N_st), norm_psi(N_st), pt2_old(N_st)
-  PROVIDE CI_electronic_energy N_det_generators 
+  PROVIDE N_det_generators 
   do k=1,N_st
     pt2(k) = 0.d0
     norm_pert(k) = 0.d0
@@ -266,7 +265,7 @@ class H_apply(object):
       double precision, intent(inout) :: select_max_out"""
 
       self.data["params_post"] += ", select_max(min(i_generator,size(select_max,1)))"
-      self.data["size_max"] = str(1024*128) 
+      self.data["size_max"] = "256"
       self.data["copy_buffer"] = """
       call copy_H_apply_buffer_to_wf
       if (s2_eig) then

scripts/get_basis.sh

@@ -8,7 +8,8 @@
 # Prints in stdout the name of a temporary file containing the basis set.
 #
 
-
+#DEBUG:
+#echo $0 $@ 1>&2
 
 if [[ -z ${QP_ROOT} ]]
 then

scripts/module/create_gitignore.sh

@@ -13,7 +13,6 @@ then
 fi
 source ${QP_ROOT}/scripts/qp_include.sh
 
-
 function do_gitingore()
 { 
   cat << EOF > .gitignore

scripts/module/module_handler.py

@@ -88,7 +88,7 @@ def get_l_module_descendant(d_child, l_module):
             except KeyError:
                 print >> sys.stderr, "Error: "
                 print >> sys.stderr, "`{0}` is not a submodule".format(module)
-                print >> sys.stderr, "Check the typo (orthograph, case, '/', etc.) "
+                print >> sys.stderr, "Check the typo (spelling, case, '/', etc.) "
                 sys.exit(1)
 
     return list(set(l))
@@ -180,6 +180,11 @@ class ModuleHandler():
     def create_png(self, l_module):
         """Create the png of the dependency tree for a l_module"""
 
+        # Don't update if we are not in the main repository
+        from is_master_repository import is_master_repository
+        if not is_master_repository:
+            return
+
         basename = "tree_dependency"
         path = '{0}.png'.format(basename)
 
@@ -289,6 +294,12 @@ if __name__ == '__main__':
                         pass
 
             if arguments["create_git_ignore"]:
+
+                # Don't update if we are not in the main repository
+                from is_master_repository import is_master_repository
+                if not is_master_repository:
+                    sys.exit()
+
                 path = os.path.join(module_abs, ".gitignore")
 
                 with open(path, "w+") as f:

scripts/module/qp_module.py

@@ -57,10 +57,48 @@ def save_new_module(path, l_child):
         f.write(D_KEY["needed_module"])
         f.write(D_KEY["documentation"])
 
+    with open(os.path.join(path, "%s.main.irp.f"%(module_name) ), "w") as f:
+        f.write("program {0}".format(module_name) )
+        f.write("""  implicit none
+  BEGIN_DOC
+! TODO
+  END_DOC
+  print *, '  _/                              '
+  print *, ' -:\_?,     _Jm####La             '
+  print *, 'J"(:" >  _]#AZ#Z#UUZ##,           '
+  print *, '_,::./   %(|i%12XmX1*1XL      _?, '
+  print *, '  \..\ _\(vmWQwodY+ia%lnL  _",/ ( '
+  print *, '   .:< ]J=mQD?WXn<uQWmmvd, -.-:=!\'
+  print *, '   "{Z jC]QW|=3Zv)Bi3BmXv3  =   _7'
+  print *, '    ]h[Z6)WQ;)jZs]C;|$BZv+, : ./  '
+  print *, '    -#sJX%$Wmm#ev]hinW#Xi:` c ;   '
+  print *, '     #X#X23###1}vI$WWmX1>|,)nr"   '
+  print *, '     4XZ#Xov1v}=)vnXAX1nnv;1n"    '
+  print *, '     ]XX#ZXoovvvivnnnlvvo2*i7     '
+  print *, '     "23Z#1S2oo2XXSnnnoSo2>v"     '
+  print *, '      miX#L -~`""!!1}oSoe|i7      '
+  print *, '      4cn#m,        v221=|v[      '
+  print *, '      ]hI3Zma,;..__wXSe=+vo       '
+  print *, '      ]Zov*XSUXXZXZXSe||vo2       '
+  print *, '      ]Z#><iiii|i||||==vn2(       '
+  print *, '      ]Z#i<ii||+|=||=:{no2[       '
+  print *, '      ]ZUsiiiiivi|=||=vo22[       '
+  print *, '      ]XZvlliiIi|i=|+|vooo        '
+  print *, '      =v1llli||||=|||||lii(       '
+  print *, '      ]iillii||||||||=>=|<        '
+  print *, '      -ziiiii||||||+||==+>        '
+  print *, '       -%|+++||=|=+|=|==/         '
+  print *, '        -a>====+|====-:-          '
+  print *, '          "~,- --   /-            '
+  print *, '            -.     )>             '
+  print *, '           .~      +-             '
+  print *, '           . .... : .             '
+  print *, '            -------~              '
+  print *, ''
+end
+""")
 
-if __name__ == '__main__':
-    arguments = docopt(__doc__)
-
+def main(arguments):
     if arguments["list"]:
 
         if arguments["--installed"]:
@@ -107,12 +145,14 @@ if __name__ == '__main__':
         save_new_module(path, l_child_reduce)
 
         print "    [ OK ]"
-        print "Your module is created in the `plugins` directory."
-        print "You need to create some `.irp.f` to be able to install it."
 #        print "` {0} install {1} `".format(os.path.basename(__file__), name)
         print ""
+        arguments["create"]=False
+        arguments["install"]=True
+        main(arguments)
 
     elif arguments["download"]:
+        print "Not yet implemented"
         pass
 #        d_local = get_dict_child([QP_SRC])
 #        d_remote = get_dict_child(arguments["<path_folder>"])
@@ -205,3 +245,8 @@ if __name__ == '__main__':
             except OSError:
                 print "%s is a core module which can't be removed" % module
 
+
+if __name__ == '__main__':
+    arguments = docopt(__doc__)
+    main(arguments)
+

scripts/module/qp_update_readme.py

@@ -44,7 +44,7 @@ def get_url(path_module_rel):
     elif is_module(path_module_rel):
         url = "http://github.com/LCPQ/quantum_package/tree/master/src"
     else:
-        print "{0} Is not a valide module nor plugin".format(path_module_rel)
+        print "{0} Is not a valid module nor plugin".format(path_module_rel)
         sys.exit(1)
 
     return os.path.join(url, path_module_rel)
@@ -155,10 +155,14 @@ def update_documentation(d_readmen, root_module):
         l_doc = []
 
         for irp in d_info[path]:
-
             url = os.path.join(get_url(os.path.basename(path)), irp.file)
             doc = extract_doc(root_module, irp.provider)
 
+            if ".irp.f_shell_" in irp.file:
+                l_doc += ["{0}".format(irp.provider),
+                          doc,
+                          ""]
+            else:
                 l_doc += ["`{0} <{1}#L{2}>`_".format(irp.provider, url, irp.line),
                           doc,
                           ""]
@@ -168,7 +172,15 @@ def update_documentation(d_readmen, root_module):
 
         d_readme[path]["documentation"] = "\n".join(l_doc_section)
 
+
 if __name__ == '__main__':
+
+    # Update documentation only if the remote repository is
+    # the main repository
+    from is_master_repository import is_master_repository
+    if not is_master_repository:
+        sys.exit(0)
+
     arguments = docopt(__doc__)
 
     if arguments["--root_module"]:
@@ -188,8 +200,8 @@ if __name__ == '__main__':
         fetch_splitted_data(d_readme, l_module_readme)
     except IOError:
         print l_module_readme, "is not a module and/or",
-        print "have not a `README.rst` file inside"
-        print "Abort..."
+        print "has not a `README.rst` file inside"
+        print "Aborting..."
         sys.exit(1)
 
     update_needed(d_readme)

scripts/pseudo/elts_num_ele.py

@@ -1,4 +1,5 @@
-name_to_elec = {"H": 1,
+name_to_elec = {"X": 0,
+                "H": 1,
                 "He": 2,
                 "Li": 3,
                 "Be": 4,

scripts/pseudo/put_pseudo_in_ezfio.py

@@ -58,6 +58,8 @@ def get_pseudo_str(l_atom):
     str_ = ""
 
     for a in l_atom:
+
+        if a is not 'X':
             l_cmd_atom = ["--atom", a]
 
             l_cmd_head = [EMSL_path, "get_basis_data",
@@ -69,6 +71,22 @@ def get_pseudo_str(l_atom):
             stdout, _ = process.communicate()
             str_ += stdout.strip() + "\n"
 
+        else: # Dummy atoms
+            str_ += """Element Symbol: X
+Number of replaced protons: 0
+Number of projectors: 0
+
+Pseudopotential data:
+
+Local component: 
+Coeff.		r^n	Exp.
+0.0  -1   0.
+0.0   1   0.
+0.0   0   0.
+
+Non-local component: 
+Coeff.		r^n	Exp.		Proj.
+"""
     return str_
 
 

scripts/utility/is_master_repository.py

@@ -0,0 +1,14 @@
+#!/usr/bin/env python
+
+import subprocess
+pipe = subprocess.Popen("git config --get remote.origin.url", \
+              shell=True, stdout=subprocess.PIPE)
+result = pipe.stdout.read()
+is_master_repository = "LCPQ/quantum_package" in result
+
+if __name__ == "__main__":
+    import sys
+    if is_master_repository:
+        sys.exit(0)
+    else:
+        sys.exit(-1)

src/AO_Basis/ao_overlap.irp.f

@@ -59,6 +59,7 @@
    enddo
   enddo
   !$OMP END PARALLEL DO
+
 END_PROVIDER
 
 

src/Bitmask/bitmasks.irp.f

@@ -262,13 +262,7 @@ END_PROVIDER
  logical                        :: exists
  integer                        :: j,i
  integer :: i_hole,i_part,i_gen
- PROVIDE ezfio_filename
-!do j = 1, N_int
-! inact_bitmask(j,1) = xor(generators_bitmask(j,1,1,1),cas_bitmask(j,1,1))
-! inact_bitmask(j,2) = xor(generators_bitmask(j,2,1,1),cas_bitmask(j,2,1))
-! virt_bitmask(j,1) = xor(generators_bitmask(j,1,2,1),cas_bitmask(j,1,1))
-! virt_bitmask(j,2) = xor(generators_bitmask(j,2,2,1),cas_bitmask(j,2,1))
-!enddo
+
  n_inact_orb = 0
  n_virt_orb = 0
  if(N_generators_bitmask == 1)then

src/Determinants/EZFIO.cfg

@@ -40,6 +40,12 @@ doc: Force the wave function to be an eigenfunction of S^2
 interface: ezfio,provider,ocaml
 default: False
 
+[threshold_davidson]
+type: Threshold
+doc: Thresholds of Davidson's algorithm
+interface: ezfio,provider,ocaml
+default: 1.e-8
+
 [threshold_generators]
 type: Threshold
 doc: Thresholds on generators (fraction of the norm) 
@@ -53,9 +59,10 @@ interface: ezfio,provider,ocaml
 default:        0.999
 
 [n_states_diag]
-type: integer
+type: States_number
 doc: n_states_diag
-interface: ezfio,provider
+default: 1
+interface: ezfio,provider,ocaml
 
 [n_int]
 interface: ezfio
@@ -89,24 +96,25 @@ doc: psi_det
 type: integer*8
 size: (determinants.n_int*determinants.bit_kind/8,2,determinants.n_det)
 
-[det_num]
-interface: ezfio,provider
-doc: det_num
-type: integer
-
 [det_occ]
 interface: ezfio,provider
 doc: det_occ
 type: integer          
-size:  (electrons.elec_alpha_num,determinants.det_num,2)
+size:  (electrons.elec_alpha_num,determinants.n_det,2)
 
 [det_coef]
 interface: ezfio,provider
 doc: det_coef
 type: double precision 
-size:  (determinants.det_num)
+size:  (determinants.n_det)
 
 [expected_s2]
 interface: ezfio,provider
-doc: expcted_s2
+doc: Expected value of S^2
 type: double precision
+
+[target_energy]
+interface: ezfio,provider,ocaml
+doc: Energy that should be obtained when truncating the wave function (optional)
+type: Energy
+default: 0.

src/Determinants/Fock_diag.irp.f

@@ -0,0 +1,84 @@
+subroutine build_fock_tmp(fock_diag_tmp,det_ref,Nint)
+  use bitmasks
+  implicit none
+  BEGIN_DOC
+! Build the diagonal of the Fock matrix corresponding to a generator
+! determinant. F_00 is <i|H|i> = E0.
+  END_DOC
+  integer, intent(in)            :: Nint
+  integer(bit_kind), intent(in)  :: det_ref(Nint,2)
+  double precision, intent(out)  :: fock_diag_tmp(2,mo_tot_num+1)
+  
+  integer                        :: occ(Nint*bit_kind_size,2)
+  integer                        :: ne(2), i, j, ii, jj
+  double precision               :: E0
+  
+  ! Compute Fock matrix diagonal elements
+  call bitstring_to_list_ab(det_ref,occ,Ne,Nint)
+
+  fock_diag_tmp = 0.d0
+  E0 = 0.d0
+
+  ! Occupied MOs
+  do ii=1,elec_alpha_num
+    i = occ(ii,1)
+    fock_diag_tmp(1,i) = fock_diag_tmp(1,i) + mo_mono_elec_integral(i,i)
+    E0 = E0 + mo_mono_elec_integral(i,i)
+    do jj=1,elec_alpha_num
+      j = occ(jj,1)
+      if (i==j) cycle
+      fock_diag_tmp(1,i) = fock_diag_tmp(1,i) + mo_bielec_integral_jj_anti(i,j)
+      E0 = E0 + 0.5d0*mo_bielec_integral_jj_anti(i,j)
+    enddo
+    do jj=1,elec_beta_num
+      j = occ(jj,2)
+      fock_diag_tmp(1,i) = fock_diag_tmp(1,i) + mo_bielec_integral_jj(i,j)
+      E0 = E0 + mo_bielec_integral_jj(i,j)
+    enddo
+  enddo
+  do ii=1,elec_beta_num
+    i = occ(ii,2)
+    fock_diag_tmp(2,i) = fock_diag_tmp(2,i) + mo_mono_elec_integral(i,i)
+    E0 = E0 + mo_mono_elec_integral(i,i)
+    do jj=1,elec_beta_num
+      j = occ(jj,2)
+      if (i==j) cycle
+      fock_diag_tmp(2,i) = fock_diag_tmp(2,i) + mo_bielec_integral_jj_anti(i,j)
+      E0 = E0 + 0.5d0*mo_bielec_integral_jj_anti(i,j)
+    enddo
+    do jj=1,elec_alpha_num
+      j = occ(jj,1)
+      fock_diag_tmp(2,i) = fock_diag_tmp(2,i) + mo_bielec_integral_jj(i,j)
+    enddo
+  enddo
+
+  ! Virtual MOs
+  do i=1,mo_tot_num
+    if (fock_diag_tmp(1,i) /= 0.d0) cycle
+    fock_diag_tmp(1,i) = fock_diag_tmp(1,i) + mo_mono_elec_integral(i,i)
+    do jj=1,elec_alpha_num
+      j = occ(jj,1)
+      fock_diag_tmp(1,i) = fock_diag_tmp(1,i) + mo_bielec_integral_jj_anti(i,j)
+    enddo
+    do jj=1,elec_beta_num
+      j = occ(jj,2)
+      fock_diag_tmp(1,i) = fock_diag_tmp(1,i) + mo_bielec_integral_jj(i,j)
+    enddo
+  enddo
+  do i=1,mo_tot_num
+    if (fock_diag_tmp(2,i) /= 0.d0) cycle
+    fock_diag_tmp(2,i) = fock_diag_tmp(2,i) + mo_mono_elec_integral(i,i)
+    do jj=1,elec_beta_num
+      j = occ(jj,2)
+      fock_diag_tmp(2,i) = fock_diag_tmp(2,i) + mo_bielec_integral_jj_anti(i,j)
+    enddo
+    do jj=1,elec_alpha_num
+      j = occ(jj,1)
+      fock_diag_tmp(2,i) = fock_diag_tmp(2,i) + mo_bielec_integral_jj(i,j)
+    enddo
+  enddo
+
+  fock_diag_tmp(1,mo_tot_num+1) = E0
+  fock_diag_tmp(2,mo_tot_num+1) = E0
+
+end