From 347e918a4a4dacc96238a52b9cbd324c51d5d1a3 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 4 Mar 2019 17:40:50 +0100 Subject: [PATCH] Add print_ci_vector in tools (#11) * Fixed energies of non-expected s2 * Moved diag_algorithm in Davdison * Fixed travis * Added print_ci_vector * Documentation * Cleaned qp_set_mo_class.ml * Removed Core in taskserver --- ocaml/TaskServer.ml | 91 +++++++++++++---------- ocaml/TaskServer.mli | 2 +- ocaml/_tags | 3 + ocaml/c_bindings.c | 69 +++++++++++++++++ ocaml/myocamlbuild.ml | 1 + ocaml/qp_print_basis.ml | 7 +- ocaml/qp_set_mo_class.ml | 79 ++++++-------------- src/determinants/s2.irp.f | 2 +- src/scf_utils/fock_matrix.irp.f | 48 ++++++------ src/scf_utils/scf_density_matrix_ao.irp.f | 6 +- src/tools/print_ci_vectors.irp.f | 35 +++++++++ 11 files changed, 213 insertions(+), 130 deletions(-) create mode 100644 ocaml/c_bindings.c create mode 100644 src/tools/print_ci_vectors.irp.f diff --git a/ocaml/TaskServer.ml b/ocaml/TaskServer.ml index 76349844..0d457e38 100644 --- a/ocaml/TaskServer.ml +++ b/ocaml/TaskServer.ml @@ -1,8 +1,5 @@ -open Core open Qptypes -module StringHashtbl = Hashtbl.Make(String) - type pub_state = | Waiting | Running of string @@ -29,15 +26,15 @@ type t = progress_bar : Progress_bar.t option ; running : bool; accepting_clients : bool; - data : string StringHashtbl.t; + data : (string, string) Hashtbl.t; } let debug_env = - match Sys.getenv "QP_TASK_DEBUG" with - | Some x -> x <> "" - | None -> false + try + Sys.getenv "QP_TASK_DEBUG"; true + with Not_found -> false let debug str = @@ -64,7 +61,7 @@ let bind_socket ~socket_type ~socket ~port = Zmq.Socket.bind socket @@ Printf.sprintf "tcp://*:%d" port; loop (-1) with - | Unix.Unix_error _ -> (Time.pause @@ Time.Span.of_sec 1. ; loop (i-1) ) + | Unix.Unix_error _ -> (Unix.sleep 1 ; loop (i-1) ) | other_exception -> raise other_exception in loop 60 @@ -77,28 +74,34 @@ let hostname = lazy ( ) +external get_ipv4_address_for_interface : string -> string = + "get_ipv4_address_for_interface" ;; + let ip_address = lazy ( - match Sys.getenv "QP_NIC" with + let interface = + try Some (Sys.getenv "QP_NIC") + with Not_found -> None + in + match interface with | None -> begin try - Lazy.force hostname - |> Unix.Inet_addr.of_string_or_getbyname - |> Unix.Inet_addr.to_string + let host = + Lazy.force hostname + |> Unix.gethostbyname + in + Unix.string_of_inet_addr host.h_addr_list.(0); with | Unix.Unix_error _ -> failwith "Unable to find IP address from host name." end | Some interface -> - begin - try - ok_exn Linux_ext.get_ipv4_address_for_interface interface - with - | Unix.Unix_error _ -> - Lazy.force hostname - |> Unix.Inet_addr.of_string_or_getbyname - |> Unix.Inet_addr.to_string - end + let result = get_ipv4_address_for_interface interface in + if String.sub result 0 5 = "error" then + Printf.sprintf "Unable to use network interface %s" interface + |> failwith + else + result ) @@ -209,7 +212,7 @@ let end_job msg program_state rep_socket pair_socket = address_inproc = None; running = true; accepting_clients = false; - data = StringHashtbl.create (); + data = Hashtbl.create 23; } and wait n = @@ -335,8 +338,10 @@ let del_task msg program_state rep_socket = and success () = let queue = - List.fold ~f:(fun queue task_id -> Queuing_system.del_task ~task_id queue) - ~init:program_state.queue task_ids + List.fold_left + (fun queue task_id -> Queuing_system.del_task ~task_id queue) + program_state.queue + task_ids in let accepting_clients = (Queuing_system.number_of_queued queue > Queuing_system.number_of_clients queue) @@ -382,11 +387,12 @@ let add_task msg program_state rep_socket = in let result = - let new_queue, new_bar = - List.fold ~f:(fun (queue, bar) task -> - Queuing_system.add_task ~task queue, - increment_progress_bar bar) - ~init:(program_state.queue, program_state.progress_bar) tasks + let new_queue, new_bar = + List.fold_left (fun (queue, bar) task -> + Queuing_system.add_task ~task queue, + increment_progress_bar bar) + (program_state.queue, program_state.progress_bar) + tasks in { program_state with queue = new_queue; @@ -547,10 +553,11 @@ let task_done msg program_state rep_socket = and success () = let new_queue, new_bar = - List.fold ~f:(fun (queue, bar) task_id -> + List.fold_left (fun (queue, bar) task_id -> Queuing_system.end_task ~task_id ~client_id queue, increment_progress_bar bar) - ~init:(program_state.queue, program_state.progress_bar) task_ids + (program_state.queue, program_state.progress_bar) + task_ids in let accepting_clients = @@ -593,7 +600,7 @@ let put_data msg rest_of_msg program_state rep_socket = in let success () = - StringHashtbl.set program_state.data ~key ~data:value ; + Hashtbl.add program_state.data key value ; Message.PutDataReply (Message.PutDataReply_msg.create ()) |> Message.to_string |> Zmq.Socket.send rep_socket; @@ -623,9 +630,8 @@ let get_data msg program_state rep_socket = let success () = let value = - match StringHashtbl.find program_state.data key with - | Some value -> value - | None -> "\000" + try Hashtbl.find program_state.data key with + | Not_found -> "\000" in Message.GetDataReply (Message.GetDataReply_msg.create ~value) |> Message.to_string_list @@ -677,13 +683,16 @@ let abort program_state rep_socket = aux [] queue 1 in let queue = - List.fold ~f:(fun queue task_id -> - Queuing_system.end_task ~task_id ~client_id queue) - ~init:queue tasks + List.fold_left + (fun queue task_id -> Queuing_system.end_task ~task_id ~client_id queue) + queue + tasks in let queue = - List.fold ~f:(fun queue task_id -> Queuing_system.del_task ~task_id queue) - ~init:queue tasks + List.fold_left + (fun queue task_id -> Queuing_system.del_task ~task_id queue) + queue + tasks in let queue = Queuing_system.del_client ~client_id queue @@ -777,7 +786,7 @@ let run ~port = address_inproc = None; progress_bar = None ; accepting_clients = false; - data = StringHashtbl.create (); + data = Hashtbl.create 23; } in diff --git a/ocaml/TaskServer.mli b/ocaml/TaskServer.mli index e3801423..e4bd87cf 100644 --- a/ocaml/TaskServer.mli +++ b/ocaml/TaskServer.mli @@ -7,7 +7,7 @@ type t = progress_bar : Progress_bar.t option ; running : bool; accepting_clients : bool; - data : (string, string) Core.Hashtbl.t ; + data : (string, string) Hashtbl.t ; } diff --git a/ocaml/_tags b/ocaml/_tags index 42843d25..1ed06ebb 100644 --- a/ocaml/_tags +++ b/ocaml/_tags @@ -1,3 +1,6 @@ true: package(core,cryptokit,zmq,str,ppx_sexp_conv,ppx_deriving,getopt) true: thread false: profile +<*byte> : linkdep(c_bindings.o), custom +<*.native>: linkdep(c_bindings.o) + diff --git a/ocaml/c_bindings.c b/ocaml/c_bindings.c new file mode 100644 index 00000000..8fb5c9b2 --- /dev/null +++ b/ocaml/c_bindings.c @@ -0,0 +1,69 @@ +#include +#include +#include +#include +#include + +#include + + + +/* Adapted from + https://github.com/monadbobo/ocaml-core/blob/master/base/core/lib/linux_ext_stubs.c +*/ + +#include +#include +#include +#include +#include +#include + +CAMLprim value get_ipv4_address_for_interface(value v_interface) +{ + CAMLparam1(v_interface); + struct ifreq ifr; + int fd = -1; + value res; + char* error = NULL; + + memset(&ifr, 0, sizeof(ifr)); + ifr.ifr_addr.sa_family = AF_INET; + /* [ifr] is already initialized to zero, so it doesn't matter if the + incoming string is too long, and [strncpy] fails to add a \0. */ + strncpy(ifr.ifr_name, String_val(v_interface), IFNAMSIZ - 1); + + caml_enter_blocking_section(); + fd = socket(AF_INET, SOCK_DGRAM, 0); + + if (fd == -1) + error = "error: couldn't allocate socket"; + else { + if (ioctl(fd, SIOCGIFADDR, &ifr) < 0) + error = "error: ioctl(fd, SIOCGIFADDR, ...) failed"; + + (void) close(fd); + } + + caml_leave_blocking_section(); + + if (error == NULL) { + /* This is weird but doing the usual casting causes errors when using + * the new gcc on CentOS 6. This solution was picked up on Red Hat's + * bugzilla or something. It also works to memcpy a sockaddr into + * a sockaddr_in. This is faster hopefully. + */ + union { + struct sockaddr sa; + struct sockaddr_in sain; + } u; + u.sa = ifr.ifr_addr; + res = caml_copy_string(inet_ntoa(u.sain.sin_addr)); + } + else + res = caml_copy_string(error); + CAMLreturn(res); + +} + + diff --git a/ocaml/myocamlbuild.ml b/ocaml/myocamlbuild.ml index 8282e794..d0909c46 100644 --- a/ocaml/myocamlbuild.ml +++ b/ocaml/myocamlbuild.ml @@ -7,6 +7,7 @@ dispatch begin function | After_rules -> begin flag ["ocaml";"compile";"native";"gprof"] (S [ A "-p"]); + pdep ["link"] "linkdep" (fun param -> [param]); end | _ -> () end diff --git a/ocaml/qp_print_basis.ml b/ocaml/qp_print_basis.ml index 2b21bf62..74c36761 100644 --- a/ocaml/qp_print_basis.ml +++ b/ocaml/qp_print_basis.ml @@ -1,11 +1,10 @@ -open Core open Qptypes let basis () = let ezfio_filename = Sys.argv.(1) in - if (not (Sys.file_exists_exn ezfio_filename)) then + if (not (Sys.file_exists ezfio_filename)) then failwith "Error reading EZFIO file"; Ezfio.set_file ezfio_filename; let basis = @@ -22,7 +21,7 @@ let mo () = let ezfio_filename = Sys.argv.(1) in - if (not (Sys.file_exists_exn ezfio_filename)) then + if (not (Sys.file_exists ezfio_filename)) then failwith "Error reading EZFIO file"; Ezfio.set_file ezfio_filename; let mo_coef = @@ -39,7 +38,7 @@ let psi_det () = let ezfio_filename = Sys.argv.(1) in - if (not (Sys.file_exists_exn ezfio_filename)) then + if (not (Sys.file_exists ezfio_filename)) then failwith "Error reading EZFIO file"; Ezfio.set_file ezfio_filename; let psi_det = diff --git a/ocaml/qp_set_mo_class.ml b/ocaml/qp_set_mo_class.ml index fe3348d9..e806082c 100644 --- a/ocaml/qp_set_mo_class.ml +++ b/ocaml/qp_set_mo_class.ml @@ -1,6 +1,5 @@ open Qputils open Qptypes -open Core (* * Command-line arguments @@ -46,7 +45,7 @@ let set ~core ~inact ~act ~virt ~del = let mo_class = - Array.init mo_num ~f:(fun i -> None) + Array.init mo_num (fun i -> None) in (* Check input data *) @@ -113,7 +112,8 @@ let set ~core ~inact ~act ~virt ~del = and av = Excitation.create_single act virt in let single_excitations = [ ia ; aa ; av ] - |> List.map ~f:Excitation.(fun x -> + |> List.map (fun x -> + let open Excitation in match x with | Single (x,y) -> ( MO_class.to_bitlist n_int (Hole.to_mo_class x), @@ -128,7 +128,8 @@ let set ~core ~inact ~act ~virt ~del = Excitation.double_of_singles aa aa ; Excitation.double_of_singles aa av ; Excitation.double_of_singles av av ] - |> List.map ~f:Excitation.(fun x -> + |> List.map (fun x -> + let open Excitation in match x with | Single _ -> assert false | Double (x,y,z,t) -> @@ -146,19 +147,20 @@ let set ~core ~inact ~act ~virt ~del = and extract_hole2 (_,_,h,_) = h and extract_particle2 (_,_,_,p) = p in + let init = Bitlist.zero n_int in let result = [ - List.map ~f:extract_hole single_excitations - |> List.fold ~init:(Bitlist.zero n_int) ~f:Bitlist.or_operator ; - List.map ~f:extract_particle single_excitations - |> List.fold ~init:(Bitlist.zero n_int) ~f:Bitlist.or_operator ; - List.map ~f:extract_hole1 double_excitations - |> List.fold ~init:(Bitlist.zero n_int) ~f:Bitlist.or_operator ; - List.map ~f:extract_particle1 double_excitations - |> List.fold ~init:(Bitlist.zero n_int) ~f:Bitlist.or_operator ; - List.map ~f:extract_hole2 double_excitations - |> List.fold ~init:(Bitlist.zero n_int) ~f:Bitlist.or_operator ; - List.map ~f:extract_particle2 double_excitations - |> List.fold ~init:(Bitlist.zero n_int) ~f:Bitlist.or_operator ; + List.map extract_hole single_excitations + |> List.fold_left Bitlist.or_operator init; + List.map extract_particle single_excitations + |> List.fold_left Bitlist.or_operator init; + List.map extract_hole1 double_excitations + |> List.fold_left Bitlist.or_operator init; + List.map extract_particle1 double_excitations + |> List.fold_left Bitlist.or_operator init; + List.map extract_hole2 double_excitations + |> List.fold_left Bitlist.or_operator init; + List.map extract_particle2 double_excitations + |> List.fold_left Bitlist.or_operator init; ] in @@ -167,10 +169,11 @@ let set ~core ~inact ~act ~virt ~del = *) (* Write masks *) - let result = List.map ~f:(fun x -> - let y = Bitlist.to_int64_list x in y@y ) + let result = + List.map (fun x -> + let y = Bitlist.to_int64_list x in y@y ) result - |> List.concat + |> List.concat in Ezfio.set_bitmasks_n_int (N_int_number.to_int n_int); @@ -194,7 +197,7 @@ let set ~core ~inact ~act ~virt ~del = let data = Array.to_list mo_class - |> List.map ~f:(fun x -> match x with + |> List.map (fun x -> match x with |None -> assert false | Some x -> MO_class.to_string x ) @@ -276,42 +279,6 @@ let run ~q ?(core="[]") ?(inact="[]") ?(act="[]") ?(virt="[]") ?(del="[]") ezfio set ~core ~inact ~act ~virt ~del -let ezfio_file = - let failure filename = - eprintf "'%s' is not an EZFIO file.\n%!" filename; - exit 1 - in - Command.Spec.Arg_type.create - (fun filename -> - match Sys.is_directory filename with - | `Yes -> - begin - match Sys.is_file (filename ^ "/.version") with - | `Yes -> filename - | _ -> failure filename - end - | _ -> failure filename - ) - - -let default range = - let failure filename = - eprintf "'%s' is not a regular file.\n%!" filename; - exit 1 - in - Command.Spec.Arg_type.create - (fun filename -> - match Sys.is_directory filename with - | `Yes -> - begin - match Sys.is_file (filename^"/.version") with - | `Yes -> filename - | _ -> failure filename - end - | _ -> failure filename - ) - - let () = let open Command_line in diff --git a/src/determinants/s2.irp.f b/src/determinants/s2.irp.f index 52c111b2..391d0073 100644 --- a/src/determinants/s2.irp.f +++ b/src/determinants/s2.irp.f @@ -33,7 +33,7 @@ subroutine get_s2(key_i,key_j,Nint,s2) implicit none use bitmasks BEGIN_DOC - ! Returns + ! Returns $\langle S^2 \rangle - S_z^2 S_z$ END_DOC integer, intent(in) :: Nint integer(bit_kind), intent(in) :: key_i(Nint,2) diff --git a/src/scf_utils/fock_matrix.irp.f b/src/scf_utils/fock_matrix.irp.f index 4c65f829..fc9eaadd 100644 --- a/src/scf_utils/fock_matrix.irp.f +++ b/src/scf_utils/fock_matrix.irp.f @@ -24,33 +24,33 @@ do j=1,elec_beta_num ! F-K - do i=1,elec_beta_num + do i=1,elec_beta_num !CC Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j))& - (Fock_matrix_mo_beta(i,j) - Fock_matrix_mo_alpha(i,j)) enddo ! F+K/2 - do i=elec_beta_num+1,elec_alpha_num + do i=elec_beta_num+1,elec_alpha_num !CA Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j))& + 0.5d0*(Fock_matrix_mo_beta(i,j) - Fock_matrix_mo_alpha(i,j)) enddo ! F - do i=elec_alpha_num+1, mo_num + do i=elec_alpha_num+1, mo_num !CV Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j)) enddo enddo do j=elec_beta_num+1,elec_alpha_num ! F+K/2 - do i=1,elec_beta_num + do i=1,elec_beta_num !AC Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j))& + 0.5d0*(Fock_matrix_mo_beta(i,j) - Fock_matrix_mo_alpha(i,j)) enddo ! F - do i=elec_beta_num+1,elec_alpha_num + do i=elec_beta_num+1,elec_alpha_num !AA Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j)) enddo ! F-K/2 - do i=elec_alpha_num+1, mo_num + do i=elec_alpha_num+1, mo_num !AV Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j))& - 0.5d0*(Fock_matrix_mo_beta(i,j) - Fock_matrix_mo_alpha(i,j)) enddo @@ -58,16 +58,16 @@ do j=elec_alpha_num+1, mo_num ! F - do i=1,elec_beta_num + do i=1,elec_beta_num !VC Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j)) enddo ! F-K/2 - do i=elec_beta_num+1,elec_alpha_num + do i=elec_beta_num+1,elec_alpha_num !VA Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j))& - 0.5d0*(Fock_matrix_mo_beta(i,j) - Fock_matrix_mo_alpha(i,j)) enddo ! F+K - do i=elec_alpha_num+1,mo_num + do i=elec_alpha_num+1,mo_num !VV Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j)) & + (Fock_matrix_mo_beta(i,j) - Fock_matrix_mo_alpha(i,j)) enddo @@ -123,22 +123,22 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_ao, (ao_num, ao_num) ] END_DOC if(frozen_orb_scf)then - call mo_to_ao(Fock_matrix_mo,size(Fock_matrix_mo,1), & - Fock_matrix_ao,size(Fock_matrix_ao,1)) - else - if ( (elec_alpha_num == elec_beta_num).and. & - (level_shift == 0.) ) & - then - integer :: i,j - do j=1,ao_num - do i=1,ao_num - Fock_matrix_ao(i,j) = Fock_matrix_ao_alpha(i,j) - enddo - enddo - else - call mo_to_ao(Fock_matrix_mo,size(Fock_matrix_mo,1), & + call mo_to_ao(Fock_matrix_mo,size(Fock_matrix_mo,1), & Fock_matrix_ao,size(Fock_matrix_ao,1)) - endif + else + if ( (elec_alpha_num == elec_beta_num).and. & + (level_shift == 0.) ) & + then + integer :: i,j + do j=1,ao_num + do i=1,ao_num + Fock_matrix_ao(i,j) = Fock_matrix_ao_alpha(i,j) + enddo + enddo + else + call mo_to_ao(Fock_matrix_mo,size(Fock_matrix_mo,1), & + Fock_matrix_ao,size(Fock_matrix_ao,1)) + endif endif END_PROVIDER diff --git a/src/scf_utils/scf_density_matrix_ao.irp.f b/src/scf_utils/scf_density_matrix_ao.irp.f index 6008c8f9..55fa8e7c 100644 --- a/src/scf_utils/scf_density_matrix_ao.irp.f +++ b/src/scf_utils/scf_density_matrix_ao.irp.f @@ -1,7 +1,7 @@ BEGIN_PROVIDER [double precision, SCF_density_matrix_ao_alpha, (ao_num,ao_num) ] implicit none BEGIN_DOC - ! S^{-1}.P_alpha.S^{-1} + ! $C.C^t$ over $\alpha$ MOs END_DOC call dgemm('N','T',ao_num,ao_num,elec_alpha_num,1.d0, & @@ -14,7 +14,7 @@ END_PROVIDER BEGIN_PROVIDER [ double precision, SCF_density_matrix_ao_beta, (ao_num,ao_num) ] implicit none BEGIN_DOC - ! S^{-1}.P_beta.S^{-1} + ! $C.C^t$ over $\beta$ MOs END_DOC call dgemm('N','T',ao_num,ao_num,elec_beta_num,1.d0, & @@ -27,7 +27,7 @@ END_PROVIDER BEGIN_PROVIDER [ double precision, SCF_density_matrix_ao, (ao_num,ao_num) ] implicit none BEGIN_DOC - ! S^{-1}.P.S^{-1} where P = C.C^t + ! Sum of $\alpha$ and $\beta$ density matrices END_DOC ASSERT (size(SCF_density_matrix_ao,1) == size(SCF_density_matrix_ao_alpha,1)) if (elec_alpha_num== elec_beta_num) then diff --git a/src/tools/print_ci_vectors.irp.f b/src/tools/print_ci_vectors.irp.f new file mode 100644 index 00000000..1682f734 --- /dev/null +++ b/src/tools/print_ci_vectors.irp.f @@ -0,0 +1,35 @@ +program print_wf + implicit none + BEGIN_DOC + ! Print the ground state wave function stored in the |EZFIO| directory + ! in the intermediate normalization. + ! + ! It also prints a lot of information regarding the excitation + ! operators from the reference determinant ! and a first-order + ! perturbative analysis of the wave function. + ! + ! If the wave function strongly deviates from the first-order analysis, + ! something funny is going on :) + END_DOC + + + ! this has to be done in order to be sure that N_det, psi_det and + ! psi_coef are the wave function stored in the |EZFIO| directory. + read_wf = .True. + touch read_wf + call routine +end + +subroutine routine + implicit none + integer :: i,k + integer :: degree + do i = 1, N_det + print *, 'Determinant ', i + call debug_det(psi_det(1,1,i),N_int) + print '(4E20.12,X)', (psi_coef(i,k), k=1,N_states) + print *, '' + print *, '' + enddo + +end