Merge pull request #6 from LCPQ/master

merge with LCPQ

configure

@@ -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
@@ -63,6 +64,8 @@ d_dependency = {
     "emsl": ["python"],
     "gcc": [],
     "g++": [],
+    "zeromq" : [ "g++" ],
+    "f77zmq" : [ "zeromq", "python" ],
     "python": [],
     "ninja": ["g++", "python"],
     "make": [],
@@ -93,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',
+    description=' patch',
-    default_path=join(QP_ROOT, "lib", "libz.a"))
+    default_path=join(QP_ROOT_BIN, "patch"))
 
 irpf90 = Info(
     url='{head}/LCPQ/irpf90/{tail}'.format(**path_github),
@@ -117,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(**
+    url='{head}/LCPQ/EMSL_Basis_Set_Exchange_Local/{tail}'.format(**path_github),
-                                                                  path_github),
     description=' EMSL basis set library',
     default_path=join(QP_ROOT_INSTALL, "emsl"))
 
@@ -131,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',
@@ -138,8 +150,9 @@ p_graphviz = Info(
 
 d_info = dict()
 
-for m in ["ocaml", "m4", "curl", "zlib", "path", "irpf90", "docopt",
+for m in ["ocaml", "m4", "curl", "zlib", "patch", "irpf90", "docopt",
-          "resultsFile", "ninja", "emsl", "ezfio", "p_graphviz"]:
+          "resultsFile", "ninja", "emsl", "ezfio", "p_graphviz",
+          "zeromq", "f77zmq" ]:
     exec ("d_info['{0}']={0}".format(m))
 
 
@@ -190,8 +203,7 @@ def check_output(*popenargs, **kwargs):
 
 def checking(d_dependency):
     """
-    For each key in d_dependency check if it
+    For each key in d_dependency check if it is avalabie 
-    is avalabie or not
     """
 
     def check_python():
@@ -261,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():
@@ -276,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)
 
@@ -329,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()))
@@ -373,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),
                     ""]
@@ -405,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"""
  _________
@@ -527,3 +548,4 @@ if __name__ == '__main__':
     create_ninja_and_rc(l_installed)
 
     recommendation()
+

install/scripts/build.sh

@@ -7,4 +7,4 @@ mkdir ${BUILD} || exit 1
 tar -zxf Downloads/${TARGET}.tar.gz --strip-components=1 --directory=${BUILD} || exit 1
 _install || exit 1
 rm -rf -- ${BUILD} _build/${TARGET}.log
-exit 0
+exit 0

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/qp_create_ezfio_from_xyz.ml

@@ -17,7 +17,7 @@ let spec =
      ~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 dummy_centers ~threshold ~molecule ~nuclei =
@@ -68,8 +68,21 @@ let dummy_centers ~threshold ~molecule ~nuclei =
     )
      
     
-   
+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 =
 
@@ -345,6 +358,13 @@ let command =
     Command.basic 
     ~summary: "Quantum Package command"
     ~readme:(fun () -> "
+
+=== 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:
@@ -354,7 +374,7 @@ otherwise specific elements can be defined as follows:
 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 d m p xyz_file () ->
        run ?o b c d m p xyz_file )

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)
 ;;
 

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/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/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 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/Hartree_Fock/EZFIO.cfg

@@ -14,7 +14,7 @@ default: 200
 type: Positive_float
 doc: Energy shift on the virtual MOs to improve SCF convergence
 interface: ezfio,provider,ocaml
-default: 0.0
+default: 0.5
 
 [mo_guess_type]
 type: MO_guess

plugins/Hartree_Fock/Fock_matrix.irp.f

@@ -344,30 +344,47 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_ao, (ao_num_align, ao_num) ]
      enddo
    enddo
  else
-   double precision, allocatable :: T(:,:), M(:,:)
+   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))
+   allocate (T(ao_num_align,ao_num),M(ao_num_align,ao_num),stat=ierr)
-   call dgemm('N','N', ao_num,ao_num,ao_num, 1.d0,                   &
+   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
@@ -380,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))
+  allocate (T(ao_num_align,ao_num),M(ao_num_align,ao_num),stat=ierr)
-  call dgemm('N','N', ao_num,ao_num,ao_num, 1.d0,                    &
+  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/NEEDED_CHILDREN_MODULES

@@ -1 +1 @@
-Integrals_Bielec MOGuess
+Integrals_Bielec MOGuess 

plugins/Hartree_Fock/damping_SCF.irp.f

@@ -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,58 +10,103 @@
    integer, allocatable           :: iwork(:)
    double precision, allocatable  :: work(:), F(:,:), S(:,:)
    
-   allocate(F(ao_num_align,ao_num), S(ao_num_align,ao_num) )
-   do j=1,ao_num
-     do i=1,ao_num
-       S(i,j) = ao_overlap(i,j)
-       F(i,j) = Fock_matrix_ao(i,j)
-     enddo
-   enddo
 
-   n = ao_num
+   if (mo_tot_num == ao_num) then
-   lwork = 1+6*n + 2*n*n
+       ! Solve H.C = E.S.C in AO basis set
-   liwork = 3 + 5*n
-   
-   allocate(work(lwork), iwork(liwork) )
 
-   lwork = -1
+       allocate(F(ao_num_align,ao_num), S(ao_num_align,ao_num) )
-   liwork = -1
+       do j=1,ao_num
+         do i=1,ao_num
+           S(i,j) = ao_overlap(i,j)
+           F(i,j) = Fock_matrix_ao(i,j)
+         enddo
+       enddo
 
-   call dsygvd(1,'v','u',ao_num,F,size(F,1),S,size(S,1),&
+       n = ao_num
-     diagonal_Fock_matrix_mo, work, lwork, iwork, liwork, info)
+       lwork = 1+6*n + 2*n*n
-!    call dsygv(1, 'v', 'u',ao_num,F,size(F,1),S,size(S,1),&
+       liwork = 3 + 5*n
-!     diagonal_Fock_matrix_mo, work, lwork, info)
+       
+       allocate(work(lwork), iwork(liwork) )
 
+       lwork = -1
+       liwork = -1
 
+       call dsygvd(1,'v','u',ao_num,F,size(F,1),S,size(S,1),&
+         diagonal_Fock_matrix_mo, work, lwork, iwork, liwork, info)
 
-   if (info /= 0) then
+       if (info /= 0) then
-     print *,  irp_here//' failed : ', info
+         print *,  irp_here//' failed : ', info
-     stop 1
+         stop 1
-   endif
+       endif
-   lwork = int(work(1))
+       lwork = int(work(1))
-   liwork = iwork(1)
+       liwork = iwork(1)
-   deallocate(work,iwork)
+       deallocate(work,iwork)
-   allocate(work(lwork), iwork(liwork) )
+       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),&
+       call dsygvd(1,'v','u',ao_num,F,size(F,1),S,size(S,1),&
-     diagonal_Fock_matrix_mo, work, lwork, iwork, liwork, info)
+         diagonal_Fock_matrix_mo, work, lwork, iwork, liwork, info)
 
-!    call dsygv(1, 'v', 'u',ao_num,F,size(F,1),S,size(S,1),&
+       if (info /= 0) then
-!     diagonal_Fock_matrix_mo, work, lwork, info)
+         print *,  irp_here//' failed : ', info
+         stop 1
+       endif
+       do j=1,mo_tot_num
+         do i=1,ao_num
+           eigenvectors_Fock_matrix_mo(i,j) = F(i,j)
+         enddo
+       enddo
 
-   if (info /= 0) then
+       deallocate(work, iwork, F, S)
-     print *,  irp_here//' failed : ', info
+
-     stop 1
+  else 
-   endif
+
-   do j=1,mo_tot_num
+       ! Solve H.C = E.C in MO basis set
-     do i=1,ao_num
+       
-       eigenvectors_Fock_matrix_mo(i,j) = F(i,j)
+       allocate( F(mo_tot_num_align,mo_tot_num) )
-     enddo
+       do j=1,mo_tot_num
-   enddo
+         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
 
-   deallocate(work, iwork, F, S)
 END_PROVIDER
  
 BEGIN_PROVIDER [double precision, diagonal_Fock_matrix_mo_sum, (mo_tot_num)]

plugins/Hartree_Fock/huckel.irp.f

@@ -13,7 +13,7 @@ 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

plugins/MP2/EZFIO.cfg

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

plugins/Perturbation/perturbation.template.f

@@ -20,14 +20,18 @@ subroutine perturb_buffer_$PERT(i_generator,buffer,buffer_size,e_2_pert_buffer,c
   integer, external              :: connected_to_ref
   logical, external              :: is_in_wavefunction
   
-  integer(bit_kind) :: minilist(Nint,2,N_det_selectors)
+  integer(bit_kind), allocatable :: minilist(:,:,:)
-  integer :: idx_minilist(N_det_selectors), N_minilist
+  integer, allocatable           :: idx_minilist(:)
+  integer                        :: N_minilist
   
-  integer(bit_kind) :: minilist_gen(Nint,2,N_det_generators)
+  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)
@@ -40,6 +44,7 @@ subroutine perturb_buffer_$PERT(i_generator,buffer,buffer_size,e_2_pert_buffer,c
   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
   
@@ -54,8 +59,6 @@ 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(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) 
 
@@ -68,11 +71,13 @@ subroutine perturb_buffer_$PERT(i_generator,buffer,buffer_size,e_2_pert_buffer,c
     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
@@ -81,20 +86,46 @@ 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)
-  do i = 1,buffer_size
+  
-
+  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)
 
     if (c_ref /= 0) then
@@ -105,19 +136,19 @@ subroutine perturb_buffer_by_mono_$PERT(i_generator,buffer,buffer_size,e_2_pert_
       cycle
     endif
     
-   integer :: degree
+    call pt2_$PERT(psi_det_generators(1,1,i_generator),buffer(1,1,i), fock_diag_tmp,        &
-    call pt2_$PERT(buffer(1,1,i),         &
+         c_pert,e_2_pert,H_pert_diag,Nint,N_minilist,n_st,minilist,idx_minilist,N_minilist) 
-        c_pert,e_2_pert,H_pert_diag,Nint,N_det_selectors,n_st)
 
     do k = 1,N_st
-      e_2_pert_buffer(k,i) = e_2_pert(k)
+      e_2_pert_buffer(k,i)   = e_2_pert(k)
-      coef_pert_buffer(k,i) = c_pert(k)
+      coef_pert_buffer(k,i)  = c_pert(k)
-      sum_norm_pert(k) += c_pert(k) * c_pert(k)
+      sum_norm_pert(k)       = sum_norm_pert(k)   + c_pert(k) * c_pert(k)
-      sum_e_2_pert(k) += e_2_pert(k)
+      sum_e_2_pert(k)        = sum_e_2_pert(k)    + e_2_pert(k)
-      sum_H_pert_diag(k) +=  H_pert_diag(k)
+      sum_H_pert_diag(k)     = sum_H_pert_diag(k) + H_pert_diag(k)
     enddo
     
-  enddo
+  enddo 
+  deallocate( minilist, minilist_gen, idx_minilist )
   
 end
 

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

scripts/compilation/qp_create_ninja.py

@@ -36,7 +36,8 @@ except ImportError:
 from qp_path import QP_ROOT, QP_SRC, QP_EZFIO
 
 LIB = "" # join(QP_ROOT, "lib", "rdtsc.o") 
-EZFIO_LIB = join(QP_ROOT, "lib", "libezfio.a") 
+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"""#
@@ -95,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("LIB = {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("")
 
@@ -261,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)]
 

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))

scripts/module/qp_install_module.py

@@ -59,10 +59,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__':
+def main(arguments):
-    arguments = docopt(__doc__)
-
     if arguments["list"]:
 
         if arguments["--installed"]:
@@ -109,12 +147,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>"])
@@ -207,3 +247,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)
+

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/density_matrix.irp.f

@@ -185,9 +185,9 @@ subroutine set_natural_mos
  allocate(tmp(size(one_body_dm_mo,1),size(one_body_dm_mo,2)))
 
  ! Negation to have the occupied MOs first after the diagonalization
- tmp = -one_body_dm_mo
+ tmp = one_body_dm_mo
  label = "Natural"
- call mo_as_eigvectors_of_mo_matrix(tmp,size(tmp,1),size(tmp,2),label)
+ call mo_as_eigvectors_of_mo_matrix(tmp,size(tmp,1),size(tmp,2),label,-1)
  deallocate(tmp)
 
 end

src/Determinants/filter_connected.irp.f

@@ -130,9 +130,7 @@ subroutine filter_connected_i_H_psi0(key1,key2,Nint,sze,idx)
     do i=1,sze
       degree_x2 = popcnt(xor( key1(1,1,i), key2(1,1))) +             &
           popcnt(xor( key1(1,2,i), key2(1,2)))
-      if (degree_x2 > 4) then
+      if (degree_x2 <= 4) then
-        cycle
-      else if(degree_x2 /= 0)then
         idx(l) = i
         l = l+1
       endif
@@ -146,9 +144,7 @@ subroutine filter_connected_i_H_psi0(key1,key2,Nint,sze,idx)
           popcnt(xor( key1(2,1,i), key2(2,1))) +                     &
           popcnt(xor( key1(1,2,i), key2(1,2))) +                     &
           popcnt(xor( key1(2,2,i), key2(2,2)))
-      if (degree_x2 > 4) then
+      if (degree_x2 <= 4) then
-        cycle
-      else if(degree_x2 /= 0)then
         idx(l) = i
         l = l+1
       endif
@@ -164,9 +160,7 @@ subroutine filter_connected_i_H_psi0(key1,key2,Nint,sze,idx)
           popcnt(xor( key1(2,2,i), key2(2,2))) +                     &
           popcnt(xor( key1(3,1,i), key2(3,1))) +                     &
           popcnt(xor( key1(3,2,i), key2(3,2)))
-      if (degree_x2 > 4) then
+      if (degree_x2 <= 4) then
-        cycle
-      else if(degree_x2 /= 0)then
         idx(l) = i
         l = l+1
       endif
@@ -174,6 +168,8 @@ subroutine filter_connected_i_H_psi0(key1,key2,Nint,sze,idx)
     
   else
     
+
+    integer, save :: icount(4) = (/0,0,0,0/)
     !DIR$ LOOP COUNT (1000)
     outer: do i=1,sze
       degree_x2 = 0
@@ -181,21 +177,17 @@ subroutine filter_connected_i_H_psi0(key1,key2,Nint,sze,idx)
       do m=1,Nint
         if ( key1(m,1,i) /=  key2(m,1)) then
           degree_x2 = degree_x2+ popcnt(xor( key1(m,1,i), key2(m,1))) 
-          if (degree_x2 > 4) then
-            cycle outer
-          endif
         endif
         if ( key1(m,2,i) /=  key2(m,2)) then
           degree_x2 = degree_x2+ popcnt(xor( key1(m,2,i), key2(m,2)))
-          if (degree_x2 > 4) then
+        endif
-            cycle outer
+        if (degree_x2 > 4) then
-          endif
+          cycle outer
         endif
       enddo
-      if(degree_x2 /= 0)then
+      idx(l) = i
-          idx(l) = i
+      l = l+1
-          l = l+1
+      icount(3) = icount(3) + 1_8
-      endif
     enddo outer
     
   endif

src/Integrals_Bielec/NEEDED_CHILDREN_MODULES

@@ -1 +1 @@
-Pseudo Bitmask
+Pseudo Bitmask ZMQ

src/Integrals_Bielec/ao_bi_integrals.irp.f

@@ -301,7 +301,7 @@ subroutine compute_ao_bielec_integrals(j,k,l,sze,buffer_value)
   double precision               :: thresh
   thresh = ao_integrals_threshold
   
-  integer                        :: n_centers, i
+  integer                        :: i
   
   if (ao_overlap_abs(j,l) < thresh) then
     buffer_value = 0._integral_kind
@@ -329,6 +329,7 @@ end
 
 BEGIN_PROVIDER [ logical, ao_bielec_integrals_in_map ]
   implicit none
+  use f77_zmq
   use map_module
   BEGIN_DOC
   !  Map of Atomic integrals
@@ -345,9 +346,8 @@ BEGIN_PROVIDER [ logical, ao_bielec_integrals_in_map ]
   integer(key_kind),allocatable  :: buffer_i(:)
   integer,parameter              :: size_buffer = 1024*64
   real(integral_kind),allocatable :: buffer_value(:)
-  integer(omp_lock_kind)         :: lock
   
-  integer                        :: n_integrals, n_centers, thread_num
+  integer                        :: n_integrals, rc
   integer                        :: jl_pairs(2,ao_num*(ao_num+1)/2), kk, m, j1, i1, lmax
   
   integral = ao_bielec_integral(1,1,1,1)
@@ -363,120 +363,61 @@ BEGIN_PROVIDER [ logical, ao_bielec_integrals_in_map ]
     endif
   endif
   
-  kk=1
-  do l = 1, ao_num   ! r2
-    do j = 1, l       ! r2
-      jl_pairs(1,kk) = j
-      jl_pairs(2,kk) = l
-      kk += 1
-    enddo
-  enddo
-  
-  PROVIDE progress_bar
-  call omp_init_lock(lock)
-  lmax = ao_num*(ao_num+1)/2
   print*, 'Providing the AO integrals'
   call wall_time(wall_0)
   call wall_time(wall_1)
   call cpu_time(cpu_1)
-  call start_progress(lmax,'AO integrals (MB)',0.d0)
+
-  !$OMP PARALLEL PRIVATE(i,j,k,l,kk,                                 &
+  integer(ZMQ_PTR)               :: zmq_socket_rep_inproc, zmq_socket_push_inproc 
-      !$OMP integral,buffer_i,buffer_value,n_integrals,              &
+  zmq_socket_rep_inproc = f77_zmq_socket(zmq_context, ZMQ_REP)
-      !$OMP cpu_2,wall_2,i1,j1,thread_num)                           &
+  rc = f77_zmq_bind(zmq_socket_rep_inproc, 'inproc://req_rep')
-      !$OMP DEFAULT(NONE)                                            &
+  if (rc /= 0) then
-      !$OMP SHARED (ao_num, jl_pairs, ao_integrals_map,thresh,       &
+    stop 'Unable to connect zmq_socket_rep_inproc'
-      !$OMP    cpu_1,wall_1,lock, lmax,n_centers,ao_nucl,            &
-      !$OMP    ao_overlap_abs,ao_overlap,abort_here,                 &
-      !$OMP    wall_0,progress_bar,progress_value,                   &
-      !$OMP    ao_bielec_integral_schwartz)
-  
-  allocate(buffer_i(size_buffer))
-  allocate(buffer_value(size_buffer))
-  n_integrals = 0
-!$  thread_num = omp_get_thread_num()
-  
-  !$OMP DO SCHEDULE(dynamic)
-  do kk=1,lmax
-IRP_IF COARRAY
-    if (mod(kk-this_image(),num_images()) /= 0) then
-      cycle
-    endif
-IRP_ENDIF
-    if (abort_here) then
-      cycle
-    endif
-    if (thread_num == 0) then
-      progress_bar(1) = kk
-    endif
-    j = jl_pairs(1,kk)
-    l = jl_pairs(2,kk)
-    j1 = j+ishft(l*l-l,-1)
-    if (ao_overlap_abs(j,l) < thresh) then
-      cycle
-    endif
-    do k = 1, ao_num           ! r1
-      i1 = ishft(k*k-k,-1)
-      if (i1 > j1) then
-        exit
-      endif
-      do i = 1, k
-        i1 += 1
-        if (i1 > j1) then
-          exit
-        endif
-        if (ao_overlap_abs(i,k)*ao_overlap_abs(j,l) < thresh) then
-          cycle
-        endif
-        if (ao_bielec_integral_schwartz(i,k)*ao_bielec_integral_schwartz(j,l) < thresh ) then
-          cycle
-        endif
-        !DIR$ FORCEINLINE
-        integral = ao_bielec_integral(i,k,j,l)
-        if (abs(integral) < thresh) then
-          cycle
-        endif
-        n_integrals += 1
-        !DIR$ FORCEINLINE
-        call bielec_integrals_index(i,j,k,l,buffer_i(n_integrals))
-        buffer_value(n_integrals) = integral
-        if (n_integrals > 1024 ) then
-          if (omp_test_lock(lock)) then
-            call insert_into_ao_integrals_map(n_integrals,buffer_i,buffer_value)
-            n_integrals = 0
-            call omp_unset_lock(lock)
-          endif
-        endif
-        if (n_integrals == size(buffer_i)) then
-          call insert_into_ao_integrals_map(n_integrals,buffer_i,buffer_value)
-          n_integrals = 0
-        endif
-      enddo
-    enddo
-    call wall_time(wall_2)
-    
-    if (thread_num == 0) then
-      if (wall_2 - wall_0 > 1.d0) then
-        wall_0 = wall_2
-        print*, 100.*float(kk)/float(lmax), '% in ',  &
-            wall_2-wall_1, 's', map_mb(ao_integrals_map) ,'MB'
-        progress_value = dble(map_mb(ao_integrals_map))
-      endif
-    endif
-  enddo
-  !$OMP END DO NOWAIT
-  call insert_into_ao_integrals_map(n_integrals,buffer_i,buffer_value)
-  deallocate(buffer_i)
-  deallocate(buffer_value)
-  !$OMP END PARALLEL
-  call omp_destroy_lock(lock)
-  call stop_progress
-  if (abort_here) then
-    stop 'Aborting in AO integrals calculation'
   endif
-IRP_IF COARRAY
+  
-  print*, 'Communicating the map'
+  integer(ZMQ_PTR)               :: thread(0:nproc)
-  call communicate_ao_integrals()
+  external                       :: ao_bielec_integrals_in_map_slave, ao_bielec_integrals_in_map_collector
-IRP_ENDIF COARRAY
+  rc = pthread_create( thread(0), ao_bielec_integrals_in_map_collector )
+  ! Create client threads
+  do i=1,nproc
+    rc = pthread_create( thread(i), ao_bielec_integrals_in_map_slave )
+  enddo
+  
+  character*(64)                 :: message_string
+  
+  do l = ao_num, 1, -1
+    rc = f77_zmq_recv( zmq_socket_rep_inproc, message_string, 64, 0)
+    print *,  l
+    ! TODO : error handling
+    ASSERT (rc >= 0)
+    ASSERT (message == 'get_ao_integrals')
+    rc = f77_zmq_send( zmq_socket_rep_inproc, l, 4, 0)
+  enddo
+  do i=1,nproc
+    rc = f77_zmq_recv( zmq_socket_rep_inproc, message_string, 64, 0)
+    ! TODO : error handling
+    ASSERT (rc >= 0)
+    ASSERT (message == 'get_ao_integrals')
+    rc = f77_zmq_send( zmq_socket_rep_inproc, 0, 4, 0)
+  enddo
+  ! TODO terminer thread(0)
+  
+  rc = f77_zmq_unbind(zmq_socket_rep_inproc, 'inproc://req_rep')
+  do i=1,nproc
+    rc = pthread_join( thread(i) )
+  enddo
+
+  zmq_socket_push_inproc = f77_zmq_socket(zmq_context, ZMQ_PUSH)
+  rc = f77_zmq_connect(zmq_socket_push_inproc, 'inproc://push_pull')
+  if (rc /= 0) then
+    stop 'Unable to connect zmq_socket_push_inproc'
+  endif
+  rc = f77_zmq_send( zmq_socket_push_inproc, -1, 4, ZMQ_SNDMORE)
+  rc = f77_zmq_send( zmq_socket_push_inproc, 0_key_kind, key_kind, ZMQ_SNDMORE)
+  rc = f77_zmq_send( zmq_socket_push_inproc, 0_integral_kind, integral_kind, 0)
+  
+  rc = pthread_join( thread(0) )
+
   print*, 'Sorting the map'
   call map_sort(ao_integrals_map)
   call cpu_time(cpu_2)
@@ -727,6 +668,7 @@ subroutine integrale_new(I_f,a_x,b_x,c_x,d_x,a_y,b_y,c_y,d_y,a_z,b_z,c_z,d_z,p,q
   
   
   implicit none
+  include 'Utils/constants.include.F'
   double precision               :: p,q
   integer                        :: a_x,b_x,c_x,d_x,a_y,b_y,c_y,d_y,a_z,b_z,c_z,d_z
   integer                        :: i, n_iter, n_pt, j
@@ -741,8 +683,9 @@ subroutine integrale_new(I_f,a_x,b_x,c_x,d_x,a_y,b_y,c_y,d_y,a_z,b_z,c_z,d_z,p,q
   p01_1 = 0.5d0/q
   p10_2 = 0.5d0 *  q /(p * q + p * p)
   p01_2 = 0.5d0 *  p /(q * q + q * p)
-  double precision               :: B10(n_pt), B01(n_pt), B00(n_pt)
+  double precision               :: B00(n_pt_max_integrals)
-  double precision               :: t1(n_pt), t2(n_pt)
+  double precision               :: B10(n_pt_max_integrals), B01(n_pt_max_integrals)
+  double precision               :: t1(n_pt_max_integrals), t2(n_pt_max_integrals)
   !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: t1, t2, B10, B01, B00
   ix = a_x+b_x
   jx = c_x+d_x
@@ -797,10 +740,11 @@ recursive subroutine I_x1_new(a,c,B_10,B_01,B_00,res,n_pt)
   !  recursive function involved in the bielectronic integral
   END_DOC
   implicit none
+  include 'Utils/constants.include.F'
   integer, intent(in)            :: a,c,n_pt
-  double precision, intent(in)   :: B_10(n_pt),B_01(n_pt),B_00(n_pt)
+  double precision, intent(in)   :: B_10(n_pt_max_integrals),B_01(n_pt_max_integrals),B_00(n_pt_max_integrals)
-  double precision, intent(out)  :: res(n_pt)
+  double precision, intent(out)  :: res(n_pt_max_integrals)
-  double precision               :: res2(n_pt)
+  double precision               :: res2(n_pt_max_integrals)
   integer                        :: i
   
   if(c<0)then
@@ -832,9 +776,10 @@ recursive subroutine I_x2_new(c,B_10,B_01,B_00,res,n_pt)
   BEGIN_DOC
   !  recursive function involved in the bielectronic integral
   END_DOC
+  include 'Utils/constants.include.F'
   integer, intent(in)            :: c, n_pt
-  double precision, intent(in)   :: B_10(n_pt),B_01(n_pt),B_00(n_pt)
+  double precision, intent(in)   :: B_10(n_pt_max_integrals),B_01(n_pt_max_integrals),B_00(n_pt_max_integrals)
-  double precision, intent(out)  :: res(n_pt)
+  double precision, intent(out)  :: res(n_pt_max_integrals)
   integer                        :: i
   
   if(c==1)then
@@ -1252,3 +1197,57 @@ recursive subroutine I_x2_pol_mult(c,B_10,B_01,B_00,C_00,D_00,d,nd,dim)
 end
 
 
+  
+
+subroutine compute_ao_integrals_jl(j,l,n_integrals,buffer_i,buffer_value)
+  implicit none
+  use map_module
+  BEGIN_DOC
+  !  Parallel client for AO integrals
+  END_DOC
+  
+  integer, intent(in)            :: j,l
+  integer,intent(out)            :: n_integrals
+  integer(key_kind),intent(out)  :: buffer_i(ao_num*ao_num)
+  real(integral_kind),intent(out) :: buffer_value(ao_num*ao_num)
+
+  integer                        :: i,k
+  double precision               :: ao_bielec_integral,cpu_1,cpu_2, wall_1, wall_2
+  double precision               :: integral, wall_0
+  double precision               :: thresh
+  integer                        :: kk, m, j1, i1
+
+  thresh = ao_integrals_threshold
+  
+  n_integrals = 0
+  
+  j1 = j+ishft(l*l-l,-1)
+  do k = 1, ao_num           ! r1
+    i1 = ishft(k*k-k,-1)
+    if (i1 > j1) then
+      exit
+    endif
+    do i = 1, k
+      i1 += 1
+      if (i1 > j1) then
+        exit
+      endif
+      if (ao_overlap_abs(i,k)*ao_overlap_abs(j,l) < thresh) then
+        cycle
+      endif
+      if (ao_bielec_integral_schwartz(i,k)*ao_bielec_integral_schwartz(j,l) < thresh ) then
+        cycle
+      endif
+      !DIR$ FORCEINLINE
+      integral = ao_bielec_integral(i,k,j,l)
+      if (abs(integral) < thresh) then
+        cycle
+      endif
+      n_integrals += 1
+      !DIR$ FORCEINLINE
+      call bielec_integrals_index(i,j,k,l,buffer_i(n_integrals))
+      buffer_value(n_integrals) = integral
+    enddo
+  enddo
+    
+end

src/Integrals_Bielec/ao_bielec_integrals_in_map_slave.irp.f

@@ -0,0 +1,99 @@
+subroutine ao_bielec_integrals_in_map_slave
+  use map_module
+  use f77_zmq
+  implicit none
+  BEGIN_DOC
+! Computes a buffer of integrals
+  END_DOC
+
+  integer                        :: j,l,n_integrals
+  integer                        :: rc
+  character*(8), external        :: zmq_port
+  integer(ZMQ_PTR)               :: zmq_socket_req_inproc, zmq_socket_push_inproc
+  real(integral_kind), allocatable :: buffer_value(:)
+  integer(key_kind), allocatable :: buffer_i(:)
+
+  allocate ( buffer_i(ao_num*ao_num), buffer_value(ao_num*ao_num) )
+
+  ! Sockets
+  zmq_socket_req_inproc = f77_zmq_socket(zmq_context, ZMQ_REQ)
+  rc = f77_zmq_connect(zmq_socket_req_inproc, 'inproc://req_rep')
+  if (rc /= 0) then
+    stop 'Unable to connect zmq_socket_req_inproc'
+  endif
+
+  zmq_socket_push_inproc = f77_zmq_socket(zmq_context, ZMQ_PUSH)
+  rc = f77_zmq_connect(zmq_socket_push_inproc, 'inproc://push_pull')
+  if (rc /= 0) then
+    stop 'Unable to connect zmq_socket_push_inproc'
+  endif
+
+
+
+  rc = f77_zmq_send( zmq_socket_req_inproc, 'get_ao_integrals', 16, 0)
+  rc = f77_zmq_recv( zmq_socket_req_inproc, l, 4, 0)
+
+  do while (l > 0) 
+    rc = f77_zmq_send( zmq_socket_req_inproc, 'get_ao_integrals', 16, 0)
+
+    do j = 1, l
+      if (ao_overlap_abs(j,l) < ao_integrals_threshold) then
+        cycle
+      endif
+      call compute_ao_integrals_jl(j,l,n_integrals,buffer_i,buffer_value) 
+      rc = f77_zmq_send( zmq_socket_push_inproc, n_integrals, 4, ZMQ_SNDMORE)
+      rc = f77_zmq_send( zmq_socket_push_inproc, buffer_i, key_kind*n_integrals, ZMQ_SNDMORE)
+      rc = f77_zmq_send( zmq_socket_push_inproc, buffer_value, integral_kind*n_integrals, 0)
+    enddo
+    rc = f77_zmq_recv( zmq_socket_req_inproc, l, 4, 0)
+  enddo
+
+  deallocate( buffer_i, buffer_value )
+
+  rc = f77_zmq_disconnect(zmq_socket_req_inproc, 'inproc://req_rep')
+end
+
+
+subroutine ao_bielec_integrals_in_map_collector
+  use map_module
+  use f77_zmq
+  implicit none
+  BEGIN_DOC
+! Collects results from the AO integral calculation
+  END_DOC
+
+  integer                        :: j,l,n_integrals
+  integer                        :: rc
+  character*(8), external        :: zmq_port
+  integer(ZMQ_PTR)               :: zmq_socket_pull_inproc
+  real(integral_kind), allocatable :: buffer_value(:)
+  integer(key_kind), allocatable :: buffer_i(:)
+
+  allocate ( buffer_i(ao_num*ao_num), buffer_value(ao_num*ao_num) )
+
+  zmq_socket_pull_inproc = f77_zmq_socket(zmq_context, ZMQ_PULL)
+  rc = f77_zmq_bind(zmq_socket_pull_inproc, 'inproc://push_pull')
+  if (rc /= 0) then
+    stop 'Unable to connect zmq_socket_pull_inproc'
+  endif
+
+  n_integrals = 0
+  do while (n_integrals >= 0)
+    
+    rc = f77_zmq_recv( zmq_socket_pull_inproc, n_integrals, 4, 0)
+    if (n_integrals > -1) then
+      rc = f77_zmq_recv( zmq_socket_pull_inproc, buffer_i, key_kind*n_integrals, 0)
+      rc = f77_zmq_recv( zmq_socket_pull_inproc, buffer_value, integral_kind*n_integrals, 0)
+      call insert_into_ao_integrals_map(n_integrals,buffer_i,buffer_value)
+    else
+      rc = f77_zmq_recv( zmq_socket_pull_inproc, buffer_i, key_kind, 0)
+      rc = f77_zmq_recv( zmq_socket_pull_inproc, buffer_value, integral_kind, 0)
+    endif
+
+  enddo
+
+  rc = f77_zmq_unbind(zmq_socket_pull_inproc, 'inproc://push_pull')
+
+  deallocate( buffer_i, buffer_value )
+end
+

src/Integrals_Bielec/map_integrals.irp.f

@@ -247,8 +247,7 @@ BEGIN_PROVIDER [ type(map_type), mo_integrals_map ]
   print*, 'MO map initialized'
 END_PROVIDER
 
-subroutine insert_into_ao_integrals_map(n_integrals,                 &
+subroutine insert_into_ao_integrals_map(n_integrals,buffer_i, buffer_values)
-      buffer_i, buffer_values)
   use map_module
   implicit none
   BEGIN_DOC
@@ -318,6 +317,7 @@ double precision function get_mo_bielec_integral_schwartz(i,j,k,l,map)
   get_mo_bielec_integral_schwartz = dble(tmp)
 end
 
+
 double precision function mo_bielec_integral(i,j,k,l)
   implicit none
   BEGIN_DOC
@@ -356,36 +356,40 @@ subroutine get_mo_bielec_integrals(j,k,l,sze,out_val,map)
     call map_get_many(map, hash, tmp_val, sze)
     ! Conversion to double precision 
     do i=1,sze
-      out_val(i) = tmp_val(i)
+      out_val(i) = dble(tmp_val(i))
     enddo
   endif
 end
 
-subroutine get_mo_bielec_integrals_existing_ik(j,l,sze,out_array,map)
+subroutine get_mo_bielec_integrals_ij(k,l,sze,out_array,map)
   use map_module
   implicit none
   BEGIN_DOC
   ! Returns multiple integrals <ij|kl> in the MO basis, all
-  ! i(1)j(1) 1/r12 k(2)l(2)
+  ! i(1)j(2) 1/r12 k(1)l(2)
-  ! i for j,k,l fixed.
+  ! i, j for k,l fixed.
   END_DOC
-  integer, intent(in)            :: j,l, sze
+  integer, intent(in)            :: k,l, sze
-  logical, intent(out)           :: out_array(sze,sze)
+  double precision, intent(out)  :: out_array(sze,sze)
   type(map_type), intent(inout)  :: map
-  integer                        :: i,k,kk,ll,m
+  integer                        :: i,j,kk,ll,m
   integer(key_kind),allocatable  :: hash(:)
   integer  ,allocatable          :: pairs(:,:), iorder(:)
+  real(integral_kind), allocatable :: tmp_val(:)
+
   PROVIDE mo_bielec_integrals_in_map
-  allocate (hash(sze*sze), pairs(2,sze*sze),iorder(sze*sze))
+  allocate (hash(sze*sze), pairs(2,sze*sze),iorder(sze*sze), &
+  tmp_val(sze*sze))
   
   kk=0
-  do k=1,sze
+  out_array = 0.d0
+  do j=1,sze
    do i=1,sze
     kk += 1
     !DIR$ FORCEINLINE
     call bielec_integrals_index(i,j,k,l,hash(kk))
     pairs(1,kk) = i
-    pairs(2,kk) = k
+    pairs(2,kk) = j
     iorder(kk) = kk
    enddo
   enddo
@@ -399,16 +403,16 @@ subroutine get_mo_bielec_integrals_existing_ik(j,l,sze,out_array,map)
     call i2radix_sort(hash,iorder,kk,-1)
   endif
 
-  call map_exists_many(mo_integrals_map, hash, kk)
+  call map_get_many(mo_integrals_map, hash, tmp_val, kk)
 
   do ll=1,kk
     m = iorder(ll)
     i=pairs(1,m)
-    k=pairs(2,m)
+    j=pairs(2,m)
-    out_array(i,k) = (hash(ll) /= 0_8)
+    out_array(i,j) = tmp_val(ll)
   enddo  
 
-  deallocate(pairs,hash,iorder)
+  deallocate(pairs,hash,iorder,tmp_val)
 end
 
 integer*8 function get_mo_map_size()
@@ -419,15 +423,6 @@ integer*8 function get_mo_map_size()
   get_mo_map_size = mo_integrals_map % n_elements
 end
 
-subroutine clear_mo_map
-  implicit none
-  BEGIN_DOC
-  ! Frees the memory of the MO map
-  END_DOC
-  call map_deinit(mo_integrals_map)
-  FREE mo_integrals_map
-end
-
 BEGIN_TEMPLATE
 
 subroutine dump_$ao_integrals(filename)

src/Integrals_Bielec/mo_bi_integrals.irp.f

@@ -102,7 +102,7 @@ subroutine add_integrals_to_map(mask_ijkl)
       !$OMP  mo_coef_transp,                                         &
       !$OMP  mo_coef_transp_is_built, list_ijkl,                     &
       !$OMP  mo_coef_is_built, wall_1, abort_here,                   &
-      !$OMP  mo_coef,mo_integrals_threshold,ao_integrals_map,mo_integrals_map,progress_bar,progress_value)
+      !$OMP  mo_coef,mo_integrals_threshold,mo_integrals_map,progress_bar,progress_value)
   n_integrals = 0
   allocate(bielec_tmp_3(mo_tot_num_align, n_j, n_k),                 &
       bielec_tmp_1(mo_tot_num_align),                                &
@@ -315,7 +315,6 @@ IRP_ENDIF
     call ezfio_set_integrals_bielec_disk_access_mo_integrals("Read")
   endif
   
-  
 end
 
 
@@ -504,3 +503,14 @@ BEGIN_PROVIDER [ double precision, mo_bielec_integral_schwartz,(mo_tot_num,mo_to
   enddo
   
 END_PROVIDER
+
+
+subroutine clear_mo_map
+  implicit none
+  BEGIN_DOC
+  ! Frees the memory of the MO map
+  END_DOC
+  call map_deinit(mo_integrals_map)
+  FREE mo_integrals_map mo_bielec_integral_schwartz mo_bielec_integral_jj mo_bielec_integral_jj_anti
+  FREE mo_bielec_integral_jj_exchange mo_bielec_integrals_in_map
+end

src/Integrals_Monoelec/pseudopot.f90

@@ -1585,148 +1585,8 @@ end
        bessel_mod_exp=x**n*accu
        end
 
-!        double precision function bessel_mod(x,n)
-!        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
-!        parameter(NBESSMAX=100)
-!        dimension SI(0:NBESSMAX),DI(0:NBESSMAX)
-!        if(n.lt.0.or.n.gt.NBESSMAX)stop 'pb with argument of bessel_mod'
-!        CALL SPHI(N,X,NBESSMAX,SI,DI)
-!        bessel_mod=si(n)
-!        end
-
-        SUBROUTINE SPHI(N,X,NMAX,SI,DI)
-!
-!       ========================================================
-!       Purpose: Compute modified spherical Bessel functions
-!                of the first kind, in(x) and in'(x)
-!       Input :  x --- Argument of in(x)
-!                n --- Order of in(x) ( n = 0,1,2,... )
-!       Output:  SI(n) --- in(x)
-!                DI(n) --- in'(x)
-!                NM --- Highest order computed
-!       Routines called:
-!                MSTA1 and MSTA2 for computing the starting
-!                point for backward recurrence
-!       ========================================================
-!
-        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
-        DIMENSION SI(0:NMAX),DI(0:NMAX)
-        NM=N
-        IF (DABS(X).LT.1.0D-100) THEN
-           DO 10 K=0,N
-              SI(K)=0.0D0
-10            DI(K)=0.0D0
-           SI(0)=1.0D0
-           DI(1)=0.333333333333333D0
-           RETURN
-        ENDIF
-        SI(0)=DSINH(X)/X
-        SI(1)=-(DSINH(X)/X-DCOSH(X))/X
-        SI0=SI(0)
-        IF (N.GE.2) THEN
-           M=MSTA1(X,200)
-
-           write(34,*)'m=',m
-
-           IF (M.LT.N) THEN
-              NM=M
-           ELSE
-              M=MSTA2(X,N,15)
-           write(34,*)'m=',m
-           ENDIF
-           F0=0.0D0
-           F1=1.0D0-100
-           DO 15 K=M,0,-1
-              F=(2.0D0*K+3.0D0)*F1/X+F0
-              IF (K.LE.NM) SI(K)=F
-              F0=F1
-15            F1=F
-           CS=SI0/F
-           write(34,*)'cs=',cs
-           DO 20 K=0,NM
-20            SI(K)=CS*SI(K)
-        ENDIF
-        DI(0)=SI(1)
-        DO 25 K=1,NM
-25         DI(K)=SI(K-1)-(K+1.0D0)/X*SI(K)
-        RETURN
-        END
 
 
-        INTEGER FUNCTION MSTA1(X,MP)
-!
-!       ===================================================
-!       Purpose: Determine the starting point for backward  
-!                recurrence such that the magnitude of    
-!                Jn(x) at that point is about 10^(-MP)
-!       Input :  x     --- Argument of Jn(x)
-!                MP    --- Value of magnitude
-!       Output:  MSTA1 --- Starting point   
-!       ===================================================
-!
-        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
-        A0=DABS(X)
-        N0=INT(1.1*A0)+1
-        F0=ENVJ(N0,A0)-MP
-        N1=N0+5
-        F1=ENVJ(N1,A0)-MP
-        DO 10 IT=1,20             
-           NN=N1-(N1-N0)/(1.0D0-F0/F1)                  
-           F=ENVJ(NN,A0)-MP
-           IF(ABS(NN-N1).LT.1) GO TO 20
-           N0=N1
-           F0=F1
-           N1=NN
- 10        F1=F
- 20     MSTA1=NN
-        RETURN
-        END
-
-
-        INTEGER FUNCTION MSTA2(X,N,MP)
-!
-!       ===================================================
-!       Purpose: Determine the starting point for backward
-!                recurrence such that all Jn(x) has MP
-!                significant digits
-!       Input :  x  --- Argument of Jn(x)
-!                n  --- Order of Jn(x)
-!                MP --- Significant digit
-!       Output:  MSTA2 --- Starting point
-!       ===================================================
-!
-        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
-        A0=DABS(X)
-        HMP=0.5D0*MP
-        EJN=ENVJ(N,A0)
-        IF (EJN.LE.HMP) THEN
-           OBJ=MP
-           N0=INT(1.1*A0)
-        ELSE
-           OBJ=HMP+EJN
-           N0=N
-        ENDIF
-        F0=ENVJ(N0,A0)-OBJ
-        N1=N0+5
-        F1=ENVJ(N1,A0)-OBJ
-        DO 10 IT=1,20
-           NN=N1-(N1-N0)/(1.0D0-F0/F1)
-           F=ENVJ(NN,A0)-OBJ
-           IF (iABS(NN-N1).LT.1) GO TO 20
-           N0=N1
-           F0=F1
-           N1=NN
-10         F1=F
-20      MSTA2=NN+10
-        RETURN
-        END
-
-        double precision FUNCTION ENVJ(N,X)
-        DOUBLE PRECISION X
-        integer N
-        ENVJ=0.5D0*DLOG10(6.28D0*N)-N*DLOG10(1.36D0*X/N)
-        RETURN
-        END
 
 !c  Computation of real spherical harmonics Ylm(theta,phi)
 !c

src/MOGuess/H_CORE_guess.irp.f

@@ -10,7 +10,7 @@ program H_CORE_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)
   print *,  'save mos'
   call save_mos
   

src/MOGuess/guess_overlap.irp.f

@@ -0,0 +1,15 @@
+program guess_mimi
+  BEGIN_DOC
+! Produce `H_core` MO orbital
+  END_DOC
+  implicit none
+  character*(64)                 :: label
+
+  mo_coef = ao_ortho_lowdin_coef
+  TOUCH mo_coef
+  label = "Guess"
+  call mo_as_eigvectors_of_mo_matrix(ao_overlap,          &
+                                     size(ao_overlap,1),  &
+                                     size(ao_overlap,2),label,-1)
+  call save_mos
+end

src/MOGuess/h_core_guess_routine.irp.f

@@ -9,7 +9,7 @@ subroutine hcore_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)
   print *,  'save mos'
   call save_mos
   SOFT_TOUCH mo_coef mo_label

src/MOGuess/truncate_mos.irp.f

@@ -0,0 +1,10 @@
+program prog_truncate_mo
+  BEGIN_DOC
+! Truncate MO set
+  END_DOC
+  implicit none
+  integer :: n
+  write(*,*) 'Number of MOs to keep'
+  read (*,*) n
+  call save_mos_truncated(n)
+end

src/MO_Basis/utils.irp.f

@@ -21,13 +21,37 @@ subroutine save_mos
   
 end
 
-subroutine mo_as_eigvectors_of_mo_matrix(matrix,n,m,label)
+subroutine save_mos_truncated(n)
   implicit none
-  integer,intent(in)             :: n,m
+  double precision, allocatable  :: buffer(:,:)
+  integer                        :: i,j,n
+  
+  call system('$QP_ROOT/scripts/save_current_mos.sh '//trim(ezfio_filename))
+  
+  call ezfio_set_mo_basis_mo_tot_num(n)
+  call ezfio_set_mo_basis_mo_label(mo_label)
+  call ezfio_set_mo_basis_ao_md5(ao_md5)
+  allocate ( buffer(ao_num,n) )
+  buffer = 0.d0
+  do j = 1, n
+    do i = 1, ao_num
+      buffer(i,j) = mo_coef(i,j)
+    enddo
+  enddo
+  call ezfio_set_mo_basis_mo_coef(buffer)
+  call ezfio_set_mo_basis_mo_occ(mo_occ)
+  deallocate (buffer)
+  
+end
+
+subroutine mo_as_eigvectors_of_mo_matrix(matrix,n,m,label,sign)
+  implicit none
+  integer,intent(in)             :: n,m, sign
   character*(64), intent(in)     :: label
   double precision, intent(in)   :: matrix(n,m)
   
-  double precision, allocatable  :: mo_coef_new(:,:), R(:,:),eigvalues(:)
+  integer :: i,j
+  double precision, allocatable  :: mo_coef_new(:,:), R(:,:),eigvalues(:), A(:,:)
   !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: mo_coef_new, R
   
   call write_time(output_mo_basis)
@@ -35,30 +59,47 @@ subroutine mo_as_eigvectors_of_mo_matrix(matrix,n,m,label)
     print *, irp_here, ': Error : m/= mo_tot_num'
     stop 1
   endif
-  allocate(R(n,m),mo_coef_new(ao_num_align,m),eigvalues(m))
+  allocate(A(n,m),R(n,m),mo_coef_new(ao_num_align,m),eigvalues(m))
+  if (sign == -1) then
+    do j=1,m
+      do i=1,n
+        A(i,j) = -matrix(i,j)
+      enddo
+    enddo
+  else
+    do j=1,m
+      do i=1,n
+        A(i,j) = matrix(i,j)
+      enddo
+    enddo
+  endif
   mo_coef_new = mo_coef
   
-  call lapack_diag(eigvalues,R,matrix,size(matrix,1),size(matrix,2))
+  call lapack_diag(eigvalues,R,A,n,m)
-  integer :: i
   write (output_mo_basis,'(A)'), 'MOs are now **'//trim(label)//'**'
   write (output_mo_basis,'(A)'), ''
   write (output_mo_basis,'(A)'), 'Eigenvalues'
   write (output_mo_basis,'(A)'), '-----------'
   write (output_mo_basis,'(A)'), ''
   write (output_mo_basis,'(A)'), '======== ================'
-  do i = 1, m
+  if (sign == -1) then
-   write (output_mo_basis,'(I8,X,F16.10)'), i,eigvalues(i)
+    do i=1,m
+      eigvalues(i) = -eigvalues(i)
+    enddo
+  endif
+  do i=1,m
+    write (output_mo_basis,'(I8,X,F16.10)'), i,eigvalues(i)
   enddo
   write (output_mo_basis,'(A)'), '======== ================'
   write (output_mo_basis,'(A)'), ''
   
   call dgemm('N','N',ao_num,m,m,1.d0,mo_coef_new,size(mo_coef_new,1),R,size(R,1),0.d0,mo_coef,size(mo_coef,1))
-  deallocate(mo_coef_new,R,eigvalues)
+  deallocate(A,mo_coef_new,R,eigvalues)
   call write_time(output_mo_basis)
   
   mo_label = label
-  SOFT_TOUCH mo_coef mo_label
 end
+
 subroutine mo_as_eigvectors_of_mo_matrix_sort_by_observable(matrix,observable,n,m,label)
   implicit none
   integer,intent(in)             :: n,m

src/Utils/LinearAlgebra.irp.f

@@ -1,3 +1,48 @@
+subroutine svd(A,LDA,U,LDU,D,Vt,LDVt,m,n)
+  implicit none
+  BEGIN_DOC
+  ! Compute A = U.D.Vt
+  !
+  ! LDx : leftmost dimension of x
+  !
+  ! Dimsneion of A is m x n
+  !
+  END_DOC
+  
+  integer, intent(in)             :: LDA, LDU, LDVt, m, n
+  double precision, intent(in)    :: A(LDA,n)
+  double precision, intent(out)   :: U(LDU,n)
+  double precision,intent(out)    :: Vt(LDVt,n)
+  double precision,intent(out)    :: D(n)
+  double precision,allocatable    :: work(:)
+  integer                         :: info, lwork, i, j, k
+  
+  double precision,allocatable    :: A_tmp(:,:)
+  allocate (A_tmp(LDA,n))
+  A_tmp = A
+  
+  ! Find optimal size for temp arrays
+  allocate(work(1))
+  lwork = -1
+  call dgesvd('A','A', n, n, A_tmp, LDA,                             &
+      D, U, LDU, Vt, LDVt, work, lwork, info)
+  lwork = work(1)
+  deallocate(work)
+
+  allocate(work(lwork))
+  call dgesvd('A','A', n, n, A_tmp, LDA,                             &
+      D, U, LDU, Vt, LDVt, work, lwork, info)
+  deallocate(work,A_tmp)
+
+  if (info /= 0) then
+    print *,  info, ': SVD failed'
+    stop
+  endif
+  
+end
+
+
+
 subroutine ortho_lowdin(overlap,LDA,N,C,LDC,m)
   implicit none
   BEGIN_DOC
@@ -29,32 +74,15 @@ subroutine ortho_lowdin(overlap,LDA,N,C,LDC,m)
   !DEC$ ATTRIBUTES ALIGN : 64    :: U, Vt, D, work
   integer                        :: info, lwork, i, j, k
   
-  double precision,allocatable   :: overlap_tmp(:,:)
+  call svd(overlap,lda,U,ldc,D,Vt,lda,m,n)
-  allocate (overlap_tmp(lda,n))
+
-  overlap_tmp = overlap
-  
-  allocate(work(1))
-  lwork = -1
-  call dgesvd('A','A', n, n, overlap_tmp, lda,                       &
-      D, U, ldc, Vt, lda, work, lwork, info)
-  lwork = work(1)
-  deallocate(work)
-  allocate(work(lwork))
-  call dgesvd('A','A', n, n, overlap_tmp, lda,                       &
-      D, U, ldc, Vt, lda, work, lwork, info)
-  deallocate(work,overlap_tmp)
-  if (info /= 0) then
-    print *,  info, ': SVD failed'
-    stop
-  endif
-  
   !$OMP PARALLEL DEFAULT(NONE) &
   !$OMP SHARED(S_half,U,D,Vt,n,C,m) &
   !$OMP PRIVATE(i,j,k)
 
   !$OMP DO
   do i=1,n
-    if ( D(i) < 1.d-6 ) then
+    if ( D(i) < 1.d-12 ) then
       D(i) = 0.d0
     else
       D(i) = 1.d0/dsqrt(D(i))

src/Utils/map_module.f90

@@ -437,97 +437,6 @@ call omp_unset_lock(map%lock)
 
 end
 
-subroutine map_update_verbose(map, key, value, sze, thr)
-  use map_module
-  implicit none
-  type (map_type), intent(inout) :: map
-  integer, intent(in)            :: sze
-  integer(key_kind), intent(inout) :: key(sze)
-  real(integral_kind), intent(inout) :: value(sze)
-  real(integral_kind), intent(in) :: thr
-  
-  integer                        :: i
-  integer(map_size_kind)         :: idx_cache, idx_cache_new
-  integer(cache_map_size_kind)   :: idx
-  integer                        :: sze2
-  integer(cache_key_kind)        :: cache_key
-  integer(map_size_kind)         :: n_elements_temp
-  type (cache_map_type)          :: local_map
-  logical                        :: map_sorted
-! do i = 1, sze
-!  print*,'value in map = ',value(i)
-! enddo
-  
-  sze2 = sze
-  map_sorted = .True.
-  
-  n_elements_temp = 0_8
-  n_elements_temp = n_elements_temp + 1_8
-  do while (sze2>0)
-    i=1
-    do while (i<=sze)
-      if (key(i) /= 0_8) then
-        idx_cache = ishft(key(i),map_shift)
-        if (omp_test_lock(map%map(idx_cache)%lock)) then
-          local_map%key => map%map(idx_cache)%key
-          local_map%value => map%map(idx_cache)%value
-          local_map%sorted = map%map(idx_cache)%sorted
-          local_map%map_size = map%map(idx_cache)%map_size
-          local_map%n_elements = map%map(idx_cache)%n_elements
-          do
-          !DIR$ FORCEINLINE
-          call search_key_big_interval(key(i),local_map%key, local_map%n_elements, idx, 1, local_map%n_elements)
-          if (idx > 0_8) then
-!           print*,'AHAAH'
-!           print*,'local_map%value(idx) = ',local_map%value(idx)
-            local_map%value(idx) = local_map%value(idx) + value(i)
-!           print*,'not a new value !'
-!           print*,'local_map%value(idx) = ',local_map%value(idx)
-          else
-            ! Assert that the map has a proper size
-            if (local_map%n_elements == local_map%map_size) then
-              call cache_map_unique(local_map)
-              call cache_map_reallocate(local_map, local_map%n_elements + local_map%n_elements)
-              call cache_map_shrink(local_map,thr)
-            endif
-            cache_key = iand(key(i),map_mask)
-            local_map%n_elements = local_map%n_elements + 1_8
-            local_map%value(local_map%n_elements) = value(i)
-!           print*,'new value !'
-            local_map%key(local_map%n_elements) = cache_key
-            local_map%sorted = .False.
-            n_elements_temp = n_elements_temp + 1_8
-          endif  ! idx > 0
-          key(i) = 0_8
-          i = i+1
-          sze2 = sze2-1
-          if (i>sze) then
-            i=1
-          endif
-          if ( (ishft(key(i),map_shift) /= idx_cache).or.(key(i)==0_8)) then
-            exit
-          endif
-        enddo
-        map%map(idx_cache)%key => local_map%key
-        map%map(idx_cache)%value => local_map%value
-        map%map(idx_cache)%sorted = local_map%sorted
-        map%map(idx_cache)%n_elements = local_map%n_elements
-        map%map(idx_cache)%map_size = local_map%map_size
-        map_sorted = map_sorted .and. local_map%sorted
-        call omp_unset_lock(map%map(idx_cache)%lock)
-      endif  ! omp_test_lock
-    else
-      i=i+1
-    endif  ! key = 0
-  enddo  ! i
-enddo  ! sze2 > 0
-call omp_set_lock(map%lock)
-map%n_elements = map%n_elements + n_elements_temp
-map%sorted = map%sorted .and. map_sorted
-call omp_unset_lock(map%lock)
-
-end
-
 subroutine map_append(map, key, value, sze)
   use map_module
   implicit none
@@ -587,13 +496,16 @@ subroutine cache_map_get_interval(map, key, value, ibegin, iend, idx)
   integer(cache_map_size_kind), intent(in) :: ibegin, iend
   real(integral_kind), intent(out) :: value
   integer(cache_map_size_kind), intent(inout) :: idx
+  double precision, pointer :: v(:)
+  integer :: i
   
-  call search_key_big_interval(key,map%key, map%n_elements, idx, ibegin, iend)
+!  call search_key_big_interval(key,map%key, map%n_elements, idx, ibegin, iend)
-  if (idx > 0) then
+  call search_key_value_big_interval(key, value, map%key, map%value, map%n_elements, idx, ibegin, iend)
-    value = map%value(idx)
+!  if (idx > 0) then
-  else
+!    value = v(idx)
-    value = 0._integral_kind
+!  else
-  endif
+!    value = 0._integral_kind
+!  endif
 end
 
 
@@ -703,7 +615,7 @@ subroutine search_key_big_interval(key,X,sze,idx,ibegin_in,iend_in)
     
     istep = ishft(iend-ibegin,-1)
     idx = ibegin + istep
-    do while (istep > 32)
+    do while (istep > 16)
       idx = ibegin + istep
       if (cache_key < X(idx)) then
         iend = idx
@@ -740,8 +652,8 @@ subroutine search_key_big_interval(key,X,sze,idx,ibegin_in,iend_in)
       endif
     enddo
     idx = ibegin
-    if (min(iend_in,sze) > ibegin+64) then
+    if (min(iend_in,sze) > ibegin+16) then
-      iend = ibegin+64
+      iend = ibegin+16
       !DIR$ VECTOR ALIGNED
       do while (cache_key > X(idx))
         idx = idx+1
@@ -784,6 +696,126 @@ subroutine search_key_big_interval(key,X,sze,idx,ibegin_in,iend_in)
   
 end
 
+subroutine search_key_value_big_interval(key,value,X,Y,sze,idx,ibegin_in,iend_in)
+  use map_module
+  implicit none
+  integer(cache_map_size_kind), intent(in) :: sze
+  integer(key_kind)           , intent(in) :: key
+  real(integral_kind)         , intent(out) :: value
+  integer(cache_key_kind)     , intent(in) :: X(sze)
+  real(integral_kind)         , intent(in) :: Y(sze)
+  integer(cache_map_size_kind), intent(in) :: ibegin_in, iend_in
+  integer(cache_map_size_kind), intent(out) :: idx
+  
+  integer(cache_map_size_kind)   :: istep, ibegin, iend, i
+  integer(cache_key_kind)        :: cache_key
+  
+  if (sze /= 0) then
+    continue
+  else
+    idx = -1
+    value = 0.d0
+    return
+  endif
+  cache_key = iand(key,map_mask)
+  ibegin = min(ibegin_in,sze)
+  iend   = min(iend_in,sze)
+  if ((cache_key > X(ibegin)) .and. (cache_key < X(iend))) then
+    
+    istep = ishft(iend-ibegin,-1)
+    idx = ibegin + istep
+    do while (istep > 16)
+      idx = ibegin + istep
+      if (cache_key < X(idx)) then
+        iend = idx
+        istep = ishft(idx-ibegin,-1)
+        idx = ibegin + istep
+        if (cache_key < X(idx)) then
+          iend = idx
+          istep = ishft(idx-ibegin,-1)
+          cycle
+        else if (cache_key > X(idx)) then
+          ibegin = idx
+          istep = ishft(iend-idx,-1)
+          cycle
+        else
+          value = Y(idx)
+          return
+        endif
+      else if (cache_key > X(idx)) then
+        ibegin = idx
+        istep = ishft(iend-idx,-1)
+        idx = idx + istep
+        if (cache_key < X(idx)) then
+          iend = idx
+          istep = ishft(idx-ibegin,-1)
+          cycle
+        else if (cache_key > X(idx)) then
+          ibegin = idx
+          istep = ishft(iend-idx,-1)
+          cycle
+        else
+          value = Y(idx)
+          return
+        endif
+      else
+        value = Y(idx)
+        return
+      endif
+    enddo
+    idx = ibegin
+    value = Y(idx)
+    if (min(iend_in,sze) > ibegin+16) then
+      iend = ibegin+16
+      !DIR$ VECTOR ALIGNED
+      do while (cache_key > X(idx))
+        idx = idx+1
+        value = Y(idx)
+      end do
+    else
+      !DIR$ VECTOR ALIGNED
+      do while (cache_key > X(idx))
+        idx = idx+1
+        value = Y(idx)
+        if (idx /= iend) then
+          cycle
+        else
+          exit
+        endif
+      end do
+    endif
+    if (cache_key /= X(idx)) then
+      idx = 1-idx
+      value = 0.d0
+    endif
+    return
+    
+  else
+    
+    if (cache_key < X(ibegin)) then
+      idx = -ibegin
+      value = 0.d0
+      return
+    endif
+    if (cache_key > X(iend)) then
+      idx = -iend
+      value = 0.d0
+      return
+    endif
+    if (cache_key == X(ibegin)) then
+      idx = ibegin
+      value = Y(idx)
+      return
+    endif
+    if (cache_key == X(iend)) then
+      idx = iend
+      value = Y(idx)
+      return
+    endif
+  endif
+  
+end
+
 
 subroutine get_cache_map_n_elements_max(map,n_elements_max)
   use map_module

src/ZMQ/NEEDED_CHILDREN_MODULES

@@ -0,0 +1 @@
+

src/ZMQ/README.rst

@@ -0,0 +1,15 @@
+===
+ZMQ
+===
+
+Socket address : defined as an environment variable : QP_RUN_ADDRESS
+
+
+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.

src/ZMQ/f77_zmq_module.f90

@@ -0,0 +1,4 @@
+module f77_zmq
+ include 'f77_zmq.h'
+end module
+

src/ZMQ/zmq.irp.f

@@ -0,0 +1,105 @@
+use f77_zmq
+
+
+BEGIN_PROVIDER [ integer(ZMQ_PTR), zmq_context ]
+ implicit none
+ BEGIN_DOC
+ ! Context for the ZeroMQ library
+ END_DOC
+ zmq_context = f77_zmq_ctx_new ()
+END_PROVIDER
+
+
+ BEGIN_PROVIDER [ character*(128), qp_run_address ]
+&BEGIN_PROVIDER [ integer, zmq_port_start ]
+ implicit none
+ BEGIN_DOC
+ ! Address of the qp_run socket
+ ! Example : tcp://130.120.229.139:12345
+ END_DOC
+ character*(128) :: buffer
+ call getenv('QP_RUN_ADDRESS',buffer)
+ if (trim(buffer) == '') then
+   stop 'QP_RUN_ADDRESS environment variable not defined'
+ endif
+
+ print *,  trim(buffer)
+ integer :: i
+ do i=len(buffer),1,-1
+   if ( buffer(i:i) == ':') then
+     qp_run_address = trim(buffer(1:i-1))
+     read(buffer(i+1:), *) zmq_port_start
+     exit
+   endif
+ enddo
+END_PROVIDER
+
+   
+function zmq_port(ishift)
+  implicit none
+  integer, intent(in)            :: ishift
+  character*(8)                  :: zmq_port
+  write(zmq_port,'(I8)') zmq_port_start+ishift
+  zmq_port = adjustl(trim(zmq_port))
+end
+
+
+BEGIN_PROVIDER [ integer(ZMQ_PTR), zmq_to_qp_run_socket ]
+ implicit none
+ BEGIN_DOC
+ ! Socket on which the qp_run process replies
+ END_DOC
+ integer                        :: rc
+ zmq_to_qp_run_socket = f77_zmq_socket(zmq_context, ZMQ_REQ)
+ rc = f77_zmq_connect(zmq_to_qp_run_socket, trim(qp_run_address))
+ if (rc /= 0) then
+   stop 'Unable to connect zmq_to_qp_run_socket'
+ endif
+ integer                        :: i
+ i=4
+ rc = f77_zmq_setsockopt(zmq_to_qp_run_socket, ZMQ_SNDTIMEO, 120000, i)
+ if (rc /= 0) then
+   stop 'Unable to set send timout in zmq_to_qp_run_socket'
+ endif
+ rc = f77_zmq_setsockopt(zmq_to_qp_run_socket, ZMQ_RCVTIMEO, 120000, i)
+ if (rc /= 0) then
+   stop 'Unable to set recv timout in zmq_to_qp_run_socket'
+ endif
+END_PROVIDER
+
+BEGIN_PROVIDER [ integer(ZMQ_PTR), zmq_socket_push ]
+  implicit none
+  BEGIN_DOC
+  ! Socket on which to push the results (1)
+  END_DOC
+  integer                        :: rc
+  character*(64)                 :: address
+  character*(8), external        :: zmq_port
+  zmq_socket_push = f77_zmq_socket(zmq_context, ZMQ_PUSH)
+  address = trim(qp_run_address)//':'//zmq_port(1)
+  rc = f77_zmq_connect(zmq_socket_push, trim(address))
+  if (rc /= 0) then
+    stop 'Unable to connect zmq_socket_push'
+  endif
+  
+END_PROVIDER
+
+BEGIN_PROVIDER [ integer(ZMQ_PTR), zmq_socket_pull ]
+  implicit none
+  BEGIN_DOC
+  ! Socket which pulls the results (2)
+  END_DOC
+  integer                        :: rc
+  character*(64)                 :: address
+  character*(8), external        :: zmq_port
+  zmq_socket_pull = f77_zmq_socket(zmq_context, ZMQ_PULL)
+  address = 'tcp://*:'//zmq_port(2)
+  rc = f77_zmq_bind(zmq_socket_pull, trim(address))
+  if (rc /= 0) then
+    stop 'Unable to connect zmq_socket_pull'
+  endif
+  
+END_PROVIDER
+
+
+