Auto-generation of libint required files

install/scripts/install_ocaml.sh

@@ -5,7 +5,7 @@ QP_ROOT=$PWD
 cd -
 
 # Normal installation
-PACKAGES="core cryptokit ocamlfind sexplib ZMQ"
+PACKAGES="core cryptokit ocamlfind sexplib pa_sexp_conv ZMQ"
 
 # Needed for ZeroMQ
 export C_INCLUDE_PATH="${QP_ROOT}"/lib:"${C_INCLUDE_PATH}"

ocaml/Atom.ml

@@ -1,4 +1,4 @@
-open Core.Std;;
+open Core.Std
 
 exception AtomError of string
 
@@ -27,12 +27,18 @@ let of_string ~units s =
       coord   = Point3d.of_string ~units (String.concat [x; y; z] ~sep:" ")
     }
   | _ -> raise (AtomError s)
-;;
+
   
 let to_string ~units a =
   [ Element.to_string a.element ;
     Charge.to_string  a.charge ;
     Point3d.to_string ~units a.coord ]
   |> String.concat ~sep:"   "
-;;
+
+
+let to_xyz a =
+  Printf.sprintf "%-3s  %s" 
+    (Element.to_string a.element)
+    (Point3d.to_string ~units:Units.Angstrom a.coord)
+
 

ocaml/Atom.mli

@@ -7,3 +7,4 @@ val sexp_of_t : t -> Sexplib.Sexp.t
 
 val of_string : units:Units.units -> string -> t
 val to_string : units:Units.units -> t -> string
+val to_xyz    : t -> string

ocaml/Basis.ml

@@ -35,11 +35,11 @@ let read_element in_channel at_number element =
   read in_channel at_number 
 
 
-let to_string b =
+
+let to_string_general ~fmt ~atom_sep b =
   let new_nucleus n = 
     Printf.sprintf "Atom %d" n
   in
-  
   let rec do_work accu current_nucleus = function
   | [] -> List.rev accu
   | (g,n)::tail -> 
@@ -47,15 +47,29 @@ let to_string b =
     in
     let accu = 
        if (n <> current_nucleus) then
-         (new_nucleus n)::""::accu
+         (new_nucleus n)::atom_sep::accu
        else
          accu
     in
-    do_work ((Gto.to_string g)::accu) n tail
+    do_work ((Gto.to_string ~fmt g)::accu) n tail
   in
   do_work [new_nucleus 1] 1 b
   |> String.concat ~sep:"\n"
 
+let to_string_gamess =
+    to_string_general ~fmt:Gto.Gamess ~atom_sep:""
+
+let to_string_gaussian b =
+  String.concat ~sep:"\n" [ "****" ;
+    to_string_general ~fmt:Gto.Gaussian ~atom_sep:"****" b ;
+    "****"
+  ]
+
+let to_string ?(fmt=Gto.Gamess) =
+  match fmt with
+  | Gto.Gamess   -> to_string_gamess
+  | Gto.Gaussian -> to_string_gaussian
+
 
 include To_md5
 let to_md5 = to_md5 sexp_of_t

ocaml/Basis.mli

@@ -14,7 +14,7 @@ val read_element :
   in_channel -> Nucl_number.t -> Element.t -> (Gto.t * Nucl_number.t) list
 
 (** Convert the basis to a string *)
-val to_string :  (Gto.t * Nucl_number.t) list -> string
+val to_string : ?fmt:Gto.fmt -> (Gto.t * Nucl_number.t) list -> string
 
 (** Convert the basis to an MD5 hash *)
 val to_md5 :  (Gto.t * Nucl_number.t) list -> MD5.t

ocaml/Gto.ml

@@ -4,6 +4,10 @@ open Qptypes
 exception GTO_Read_Failure of string
 exception End_Of_Basis
 
+type fmt =
+| Gamess
+| Gaussian
+
 type t =
 { sym  : Symmetry.t ;
   lc   : ((Primitive.t * AO_coef.t) list)
@@ -68,8 +72,8 @@ let read_one in_channel =
 
 
 
-(** Transform the gto to a string *)
-let to_string { sym = sym ; lc = lc } =
+(** Write the GTO in Gamess format *)
+let to_string_gamess { sym = sym ; lc = lc } =
   let result = 
     Printf.sprintf "%s %3d" (Symmetry.to_string sym) (List.length lc)
   in
@@ -88,3 +92,30 @@ let to_string { sym = sym ; lc = lc } =
   |> String.concat ~sep:"\n"
 
 
+(** Write the GTO in Gaussian format *)
+let to_string_gaussian { sym = sym ; lc = lc } =
+  let result = 
+    Printf.sprintf "%s %3d   1.00" (Symmetry.to_string sym) (List.length lc)
+  in
+  let rec do_work accu i = function
+  | [] -> List.rev accu
+  | (p,c)::tail -> 
+    let p = AO_expo.to_float p.Primitive.expo
+    and c = AO_coef.to_float c
+    in
+    let result = 
+      Printf.sprintf "%15.7f       %15.7f" p c
+    in
+    do_work (result::accu) (i+1) tail
+  in
+  (do_work [result] 1 lc)
+  |> String.concat ~sep:"\n"
+
+
+(** Transform the gto to a string *)
+let to_string ?(fmt=Gamess) =
+  match fmt with
+  | Gamess   -> to_string_gamess
+  | Gaussian -> to_string_gaussian
+
+

ocaml/Gto.mli

@@ -1,5 +1,9 @@
 exception GTO_Read_Failure of string
 exception End_Of_Basis
+type fmt =
+| Gamess
+| Gaussian
+
 type t =
   { sym : Symmetry.t ;
     lc  : (Primitive.t * Qptypes.AO_coef.t) list;
@@ -13,4 +17,4 @@ val of_prim_coef_list :
 val read_one : in_channel -> t
 
 (** Convert to string for printing *)
-val to_string : t -> string
+val to_string : ?fmt:fmt -> t -> string

ocaml/Input_ao_basis.ml

@@ -17,6 +17,7 @@ module Ao_basis : sig
   ;;
   val read : unit -> t option
   val to_string : t -> string
+  val to_basis  : t -> Basis.t
   val write  : t -> unit
   val to_md5 : t -> MD5.t
   val to_rst : t -> Rst_string.t

ocaml/Input_nuclei.ml

@@ -13,6 +13,7 @@ module Nuclei : sig
   val read  : unit -> t option
   val write : t -> unit 
   val to_string : t -> string
+  val to_atom_list : t -> Atom.t list
   val to_rst : t -> Rst_string.t
   val of_rst : Rst_string.t -> t option
 end = struct
@@ -134,6 +135,22 @@ end = struct
   ;;         
 
 
+  let to_atom_list b =
+     let rec loop accu (coord, charge, label) = function 
+        | -1 -> accu 
+        | i ->  
+            let atom = 
+               { Atom.element = label.(i) ; 
+                 Atom.charge  = charge.(i) ; 
+                 Atom.coord   = coord.(i) ; 
+               } 
+            in 
+            loop (atom::accu)  (coord, charge, label) (i-1) 
+     in 
+     loop [] (b.nucl_coord, b.nucl_charge, b.nucl_label) 
+          ( (Nucl_number.to_int b.nucl_num) - 1) 
+  ;;
+
   let to_string b =
     Printf.sprintf "
 nucl_num         = %s

ocaml/Libint.ml

@@ -0,0 +1,85 @@
+open Core.Std
+open Qptypes
+
+let molecule = lazy (
+    let atom_list = 
+        match Input.Nuclei.read () with
+        | Some data ->  Input.Nuclei.to_atom_list data
+        | None -> failwith "No coordinate found"
+    and data = 
+        match Input.Electrons.read () with
+        | Some data -> data
+        | None -> failwith "No electrons found"
+    in
+    { Molecule.nuclei = atom_list;
+      Molecule.elec_alpha = data.Input.Electrons.elec_alpha_num;
+      Molecule.elec_beta = data.Input.Electrons.elec_beta_num
+    }
+)
+
+let write_xyz_file libint_dir =
+    assert (Sys.is_directory_exn libint_dir);
+    let filename = 
+       Filename.concat libint_dir "xyz" 
+    and molecule = 
+        Lazy.force molecule
+    in
+    Out_channel.with_file  filename ~f:(fun oc -> 
+        Molecule.to_xyz molecule
+        |> Out_channel.output_string oc
+    )
+
+let write_basis_file libint_dir = 
+    assert (Sys.is_directory_exn libint_dir);
+    let filename = 
+       Filename.concat libint_dir "basis.g94" 
+    and molecule = 
+        Lazy.force molecule
+    in
+    let text = 
+       let rec substitute accu i = function
+       | ele :: tail -> 
+         let pattern = 
+           Printf.sprintf "Atom %d" i
+         in
+         let new_string = 
+            String.substr_replace_first accu ~pattern ~with_:(Printf.sprintf "%s  0" ele)
+         in
+         substitute new_string (i+1) tail
+       | [] -> accu
+       in
+       let accu = 
+          let b = 
+            match Input.Ao_basis.read () with
+            | Some data -> Input.Ao_basis.to_basis data
+            | None -> failwith "No AO basis"
+          in
+          Basis.to_string ~fmt:Gto.Gaussian b
+       and atom_names =
+          List.map molecule.Molecule.nuclei ~f:(fun x -> Element.to_string x.Atom.element)
+       in
+       substitute accu 1 atom_names
+    in
+    Out_channel.with_file  filename ~f:(fun oc -> 
+        Out_channel.output_string oc text
+    )
+    
+
+let write_files ezfio_filename =
+    assert (Sys.is_directory_exn ezfio_filename);
+
+    let libint_dir =
+      Filename.concat ezfio_filename "libint" 
+    in
+
+    let () =
+      match Sys.is_directory libint_dir with
+      | `Yes -> ()
+      | `No -> Unix.mkdir libint_dir
+      | `Unknown -> failwith ("Unable to tell if "^libint_dir^" exists.")
+    in
+
+    write_xyz_file libint_dir;
+    write_basis_file libint_dir
+
+

ocaml/Molecule.ml

@@ -85,7 +85,7 @@ let name m =
   String.concat (result)
 
 
-let to_string m =
+let to_string_general ~f m = 
   let { nuclei  ; elec_alpha ; elec_beta } = m
   in
   let n =
@@ -94,10 +94,15 @@ let to_string m =
   let title =
      name m
   in
-  [ Int.to_string n ; title ] @ 
-   (List.map ~f:(fun x -> Atom.to_string Units.Angstrom x) nuclei)
+  [ Int.to_string n ; title ] @  (List.map ~f nuclei)
   |> String.concat ~sep:"\n"
 
+let to_string =
+  to_string_general ~f:(fun x -> Atom.to_string Units.Angstrom x)
+
+let to_xyz = 
+  to_string_general ~f:Atom.to_xyz 
+
 
 let of_xyz_string
     ?(charge=(Charge.of_int 0)) ?(multiplicity=(Multiplicity.of_int 1))

ocaml/Molecule.mli

@@ -20,6 +20,7 @@ val name : t -> string
 
 (** Conversion for printing *)
 val to_string : t -> string
+val to_xyz    : t -> string
 
 
 (** Creates a molecule from an xyz file *)

ocaml/TaskServer.mli

@@ -0,0 +1,84 @@
+type t =
+{
+    queue           : Queuing_system.t ;
+    state           : Message.State.t option ;
+    address_tcp     : Address.Tcp.t option ;
+    address_inproc  : Address.Inproc.t option ;
+    psi             : Message.Psi.t option;
+    progress_bar    : Progress_bar.t option ;
+    running         : bool;
+}
+
+
+(** {1} Debugging *)
+
+(** Fetch the QP_TASK_DEBUG environment variable *)
+val debug_env : bool
+
+(** Print a debug message *)
+val debug : string -> unit
+
+(** {1} ZMQ *)
+
+(** ZeroMQ context *)
+val zmq_context : ZMQ.Context.t
+
+(** Bind a ZMQ socket *)
+val bind_socket :
+  socket_type:string -> socket:'a ZMQ.Socket.t -> address:string -> unit
+
+(** Name of the host on which the server runs *)
+val hostname : string lazy_t
+
+(** IP address of the current host *)
+val ip_address : string lazy_t
+
+(** Standard messages *)
+val reply_ok : [> `Req ] ZMQ.Socket.t -> unit
+val reply_wrong_state : [> `Req ] ZMQ.Socket.t -> unit
+
+(** Stop server *)
+val stop : port:int -> unit
+
+(** {1} Server functions *)
+
+(** Create a new job *)
+val new_job : Message.Newjob_msg.t -> t -> [> `Req ] ZMQ.Socket.t -> t
+
+(** Finish a running job *)
+val end_job : Message.Endjob_msg.t -> t -> [> `Req ] ZMQ.Socket.t -> t
+
+(** Connect a client *)
+val connect: Message.Connect_msg.t -> t -> [> `Req ] ZMQ.Socket.t -> t
+
+(** Disconnect a client *)
+val disconnect: Message.Disconnect_msg.t -> t -> [> `Req ] ZMQ.Socket.t -> t
+
+(** Add a task to the pool *)
+val add_task: Message.AddTask_msg.t -> t -> [> `Req ] ZMQ.Socket.t -> t
+
+(** Mark the task as done by the client *)
+val task_done: Message.TaskDone_msg.t -> t -> [> `Req ] ZMQ.Socket.t -> t
+
+(** Delete a task when it has been pulled by the collector *)
+val del_task: Message.DelTask_msg.t -> t -> [> `Req ] ZMQ.Socket.t -> t
+
+(** The client get a new task to execute *)
+val get_task: Message.GetTask_msg.t -> t -> [> `Req ] ZMQ.Socket.t -> t
+
+(** Terminate server *)
+val terminate : t -> [> `Req ] ZMQ.Socket.t -> t
+
+(** Put a wave function in the task server *)
+val put_psi :
+   Message.PutPsi_msg.t -> string list -> t -> [> `Req ] ZMQ.Socket.t -> t
+
+(** Get the wave function stored in the task server *)
+val get_psi : Message.GetPsi_msg.t -> t -> [> `Req ] ZMQ.Socket.t -> t
+
+(** Reply an Error message *)
+val error : string -> t -> [> `Req ] ZMQ.Socket.t -> t
+
+(** Run server *)
+val run : port:int -> unit
+

ocaml/_tags

@@ -1,3 +1,3 @@
-true: package(core,sexplib.syntax,cryptokit,ZMQ)
+true: package(core,sexplib,pa_sexp_conv,cryptokit,ZMQ)
 true: thread
 false: profile

plugins/Hartree_Fock/debug_libint.irp.f

@@ -1,20 +1,10 @@
  program debug_libint
-   use libint_module
 
    implicit none
    double precision :: ao_bielec_integral
 
-   integer, allocatable :: s2bf(:)
    double precision, allocatable :: buffer_int(:)
 
-   call init_libint(trim(ezfio_filename)//char(0))
-
-   integer :: nb_shell_f
-   nb_shell_f = get_nb_shell()
-
-   allocate(s2bf(2*nb_shell_f))
-   call map_shell_to_basis_function_interval(2*nb_shell_f,s2bf)
-
    integer :: s1, s2,s3,s4
    integer :: bf1,bf2,bf3,bf4
    integer :: bf1_begin,bf2_begin,bf3_begin,bf4_begin
@@ -26,30 +16,30 @@
    ! Loop over the shell !
    ! =================== !
 
-   do s1 = 1, nb_shell_f
+   do s1 = 1, shell_num
 
-      print*, s1, "/", nb_shell_f
+      print*, s1, "/", shell_num
 
-      bf1_begin = s2bf(2*s1-1)
-      bf1_end = s2bf(2*s1)
+      bf1_begin = shell_idx(1,s1)
+      bf1_end = shell_idx(2,s1)
       n1 = 1 + bf1_end - bf1_begin
    
-      do s2 = 1, nb_shell_f
+      do s2 = 1, shell_num
    
-         bf2_begin = s2bf(2*s2-1)
-         bf2_end = s2bf(2*s2)
+         bf2_begin = shell_idx(1,s2)
+         bf2_end = shell_idx(2,s2)
          n2 = 1 + bf2_end - bf2_begin
    
-         do s3 = 1, nb_shell_f
+         do s3 = 1, shell_num
      
-            bf3_begin = s2bf(2*s3-1)
-            bf3_end = s2bf(2*s3)
+            bf3_begin = shell_idx(1,s3)
+            bf3_end = shell_idx(2,s3)
             n3 = 1 + bf3_end - bf3_begin
    
-            do s4 = 1, nb_shell_f
+            do s4 = 1, shell_num
    
-               bf4_begin = s2bf(2*s4-1)
-               bf4_end = s2bf(2*s4)     
+               bf4_begin = shell_idx(1,s4)
+               bf4_end = shell_idx(2,s4)     
                n4 = 1 + bf4_end - bf4_begin
    
                ! ========================== !
@@ -93,19 +83,19 @@
                            libint = buffer_int(f1234) * norm
 
                            !Verify with the manu's one
-!                           ref = ao_bielec_integral(bf1,bf2,bf3,bf4)
-!   
-!                           if ( (ABS(ABS(ref) - ABS(libint)).GE.1e-6) )  THEN
-!
-!                              print*, bf1,bf2,bf3,bf4
-!                              print*,"r", ref
-!                              print*,"l", libint
-!                              print*,"r/l", ref/libint
-!                              print*,"l/r", libint/ref
-!                              print*,"n", norm
-!
-!                              call exit(1)
-!                           end if
+                           ref = ao_bielec_integral(bf1,bf2,bf3,bf4)
+   
+                           if ( (ABS(ABS(ref) - ABS(libint)).GE.1e-6) )  THEN
+
+                              print*, bf1,bf2,bf3,bf4
+                              print*,"r", ref
+                              print*,"l", libint
+                              print*,"r/l", ref/libint
+                              print*,"l/r", libint/ref
+                              print*,"n", norm
+
+                              call exit(1)
+                           end if
 
                          enddo
                      enddo
@@ -121,4 +111,4 @@
    enddo
    
    call finalize_libint()
- end debug_libint
+ end 

scripts/compilation/qp_create_ninja.py

@@ -38,7 +38,7 @@ from qp_path import QP_ROOT, QP_SRC, QP_EZFIO
 LIB = "" # join(QP_ROOT, "lib", "rdtsc.o") 
 EZFIO_LIB = join(QP_ROOT, "lib", "libezfio_irp.a") 
 ZMQ_LIB = join(QP_ROOT, "lib", "libf77zmq.a") + " "  + join(QP_ROOT, "lib", "libzmq.a") + " -lstdc++ -lrt"
-INT_LIB = join(QP_ROOT, "libint","lib", ".libs", "libint2.a")
+INT_LIB = join(QP_ROOT, "install", "libint","lib", ".libs", "libint2.a")
 ROOT_BUILD_NINJA = join(QP_ROOT, "config", "build.ninja")
 
 header = r"""#

scripts/ezfio_interface/qp_edit_template

@@ -6,7 +6,7 @@ open Core.Std;;
 
 WARNING 
 This file is autogenerad by
-`${{QP_ROOT}}/script/ezfio_interface/ei_handler.py`
+`${{QP_ROOT}}/scripts/ezfio_interface/ei_handler.py`
 *)
 
 
@@ -120,7 +120,7 @@ let set str s =
         | Nuclei           -> write Nuclei.(of_rst, write) s
         | Ao_basis         -> () (* TODO *)
         | Mo_basis         -> () (* TODO *)
-    end 
+    end
 ;;
 
 
@@ -169,7 +169,9 @@ let run check_only ezfio_filename =
   in
 
   (* Create the temp file *)
-  let temp_filename = create_temp_file ezfio_filename tasks in
+  let temp_filename =
+      create_temp_file ezfio_filename tasks
+  in
 
   (* Open the temp file with external editor *)
   let editor = 
@@ -191,9 +193,10 @@ let run check_only ezfio_filename =
       In_channel.input_all in_channel) 
   in
   List.iter ~f:(fun x -> set temp_string x) tasks;
+  Libint.write_files (!Ezfio.ezfio_filename);
 
   (* Remove temp_file *)
-  Sys.remove temp_filename;
+  Sys.remove temp_filename
 ;;
 
 

src/Integrals_Bielec/integral_bielec.cc

@@ -18,7 +18,7 @@
  */
 
 #if __cplusplus <= 199711L
-# error "Hartree-Fock test requires C++11 support"
+# error "C++11 support is required"
 #endif
 
 // standard C++ headers
@@ -244,4 +244,4 @@ void compute_ao_bielec_integrals_shell(int s1, int s2, int s3, int s4, int sze,
 
   for(auto i=0; i!=sze; i++)
     values[i] = buf_1234[i];
-  };
+};

src/Integrals_Monoelec/shell.irp.f

@@ -0,0 +1,30 @@
+use libint_module
+BEGIN_PROVIDER [ logical, libint_initialized ]
+ implicit none
+ BEGIN_DOC
+ ! if true, libint is initialized
+ END_DOC
+ call init_libint(trim(ezfio_filename)//char((0)))
+END_PROVIDER
+
+BEGIN_PROVIDER [ integer, shell_num ]
+ implicit none
+ BEGIN_DOC
+ ! Number of shells
+ END_DOC
+ PROVIDE libint_initialized
+ shell_num = get_nb_shell()
+
+END_PROVIDER
+
+
+BEGIN_PROVIDER [ integer, shell_idx, (2,shell_num) ]
+ implicit none
+ BEGIN_DOC
+ ! Contains the 1st and last AO index in each shell
+ END_DOC
+ PROVIDE libint_initialized
+ call map_shell_to_basis_function_interval(2*shell_num,shell_idx)
+END_PROVIDER
+
+