Added Huckel guess and fixed compatibility with new EMSL

data/ezfio_defaults

@@ -45,6 +45,7 @@ cassd
 hartree_fock
   n_it_scf_max               200
   thresh_scf                 1.e-10
+  guess                      "Huckel"
 
 cisd_selected
   n_det_max_cisd         10000

ocaml/Input_hartree_fock.ml

@@ -6,6 +6,7 @@ module Hartree_fock : sig
   type t = 
     { n_it_scf_max       : Strictly_positive_int.t;
       thresh_scf         : Threshold.t;
+      guess              : MO_guess.t;
     } with sexp
   ;;
   val read  : unit -> t option
@@ -17,6 +18,7 @@ end = struct
   type t = 
     { n_it_scf_max       : Strictly_positive_int.t;
       thresh_scf         : Threshold.t;
+      guess              : MO_guess.t;
     } with sexp
   ;;
 
@@ -52,34 +54,57 @@ end = struct
     |> Ezfio.set_hartree_fock_thresh_scf
   ;;
 
+  let read_guess () = 
+    if not (Ezfio.has_hartree_fock_guess ()) then
+      get_default "guess"
+      |> String.strip ~drop:(fun x -> x = '"')
+      |> Ezfio.set_hartree_fock_guess
+    ;
+    Ezfio.get_hartree_fock_guess ()
+    |> MO_guess.of_string
+  ;;
+
+  let write_guess guess =
+    MO_guess.to_string guess
+    |> Ezfio.set_hartree_fock_guess
+  ;;
 
   let read () = 
     Some
     { n_it_scf_max       = read_n_it_scf_max ();
       thresh_scf         = read_thresh_scf ();
+      guess              = read_guess ();
     }
   ;;
 
 
   let write { n_it_scf_max ;
               thresh_scf ;
+              guess      ;
             } =
     write_n_it_scf_max n_it_scf_max;
-    write_thresh_scf   thresh_scf
+    write_thresh_scf thresh_scf; 
+    write_guess  guess
   ;;
 
 
   let to_string b =
     Printf.sprintf "
 n_it_scf_max = %s
-thresh_scf = %s
+thresh_scf  = %s
+guess       = %s
 "
       (Strictly_positive_int.to_string b.n_it_scf_max)
       (Threshold.to_string b.thresh_scf)
+      (MO_guess.to_string b.guess)
   ;;
 
   let to_rst b =
     Printf.sprintf "
+Type of MO guess [ Huckel | HCore ] ::
+
+  guess = %s
+
 Max number of SCF iterations ::
 
   n_it_scf_max = %s
@@ -89,8 +114,9 @@ SCF convergence criterion (on energy) ::
   thresh_scf = %s
 
 "
-        (Strictly_positive_int.to_string b.n_it_scf_max)
+      (MO_guess.to_string b.guess)
-        (Threshold.to_string b.thresh_scf)
+      (Strictly_positive_int.to_string b.n_it_scf_max)
+      (Threshold.to_string b.thresh_scf)
     |> Rst_string.of_string
   ;;
 

ocaml/qp_create_ezfio_from_xyz.ml

@@ -60,20 +60,8 @@ let run ?o b c m xyz_file =
       | None -> (* Principal basis *)
         let basis = elem_and_basis_name in
         let command =
-          let rec apply accu = function
-          | [] -> accu
-          | atom::rest -> 
-            let new_accu =
-              accu ^ " " ^ (Element.to_string atom)
-            in
-            apply new_accu rest
-          in
-          let accu = 
             Qpackage.root ^ "/scripts/get_basis.sh \"" ^ temp_filename 
               ^ "\" \"" ^ basis ^"\""
-          in
-          List.map nuclei ~f:(fun x -> x.Atom.element)
-          |> apply accu 
         in
         begin
           let filename = 

ocaml/qptypes_generator.ml

@@ -138,7 +138,32 @@ let input_ezfio = "
 ;;
 
 let untouched = "
+module MO_guess : sig
+  type t with sexp
+  val to_string : t -> string
+  val of_string : string -> t
+end = struct
+  type t = 
+  | Huckel
+  | HCore
+  with sexp
+
+  let to_string = function
+  | Huckel -> \"Huckel\"
+  | HCore  -> \"HCore\"
+
+  let of_string  s = 
+    let s = 
+      String.lowercase s
+    in
+    match s with
+    | \"huckel\" -> Huckel
+    | \"hcore\"  -> HCore
+    | _ -> failwith (\"Wrong Guess type : \"^s)
+
+end
 "
+;;
 
 let template = format_of_string "
 module %s : sig

scripts/get_basis.sh

@@ -13,7 +13,9 @@ export PYTHONPATH="${EMSL_API_ROOT}":${PYTHONPATH}
 
 tmpfile="$1"
 shift
-basis="$1"
+
+# Case insensitive basis in input
+basis=$( ${EMSL_API_ROOT}/EMSL_api.py list_basis | cut -d "'" -f 2 | grep -i "^${1}\$")
 shift
 atoms=""
 

src/Hartree_Fock/Huckel_guess.irp.f

@@ -0,0 +1,6 @@
+program guess
+  implicit none
+  character*(64)                 :: label
+  call huckel_guess
+  
+end

src/Hartree_Fock/README.rst

@@ -2,6 +2,7 @@
 Hartree-Fock Module
 ===================
 
+From the 140 molecules of the G2 set, only LiO, ONa don't converge well.
 
 Needed Modules
 ==============

src/Hartree_Fock/SCF.irp.f

@@ -8,17 +8,24 @@ end
 subroutine create_guess
   implicit none
   BEGIN_DOC
-! Create an H_core guess if no MOs are present in the EZFIO directory
+! Create an MO guess if no MOs are present in the EZFIO directory
   END_DOC
   logical                        :: exists
   PROVIDE ezfio_filename
   call ezfio_has_mo_basis_mo_coef(exists)
   if (.not.exists) then
-    mo_coef = ao_ortho_lowdin_coef
+    if (mo_guess_type == "HCore") then
-    TOUCH mo_coef
+      mo_coef = ao_ortho_lowdin_coef
-    mo_label = 'Guess'
+      TOUCH mo_coef
-    call mo_as_eigvectors_of_mo_matrix(mo_mono_elec_integral,size(mo_mono_elec_integral,1),size(mo_mono_elec_integral,2),mo_label)
+      mo_label = 'Guess'
-    SOFT_TOUCH mo_coef mo_label
+      call mo_as_eigvectors_of_mo_matrix(mo_mono_elec_integral,size(mo_mono_elec_integral,1),size(mo_mono_elec_integral,2),mo_label)
+      SOFT_TOUCH mo_coef mo_label
+    else if (mo_guess_type == "Huckel") then
+      call huckel_guess
+    else
+      print *,  'Unrecognized MO guess type : '//mo_guess_type
+      stop 1
+    endif
   endif
 end
 

src/Hartree_Fock/damping_SCF.irp.f

@@ -71,7 +71,7 @@ subroutine damping_SCF
     delta_alpha = D_new_alpha - D_alpha
     delta_beta  = D_new_beta  - D_beta 
     
-    lambda = 0.5d0
+    lambda = 1.d0
     E_half = 0.d0
     do while (E_half > E)
       HF_density_matrix_ao_alpha = D_alpha + lambda * delta_alpha
@@ -118,8 +118,8 @@ subroutine damping_SCF
   
   call mo_as_eigvectors_of_mo_matrix(Fock_matrix_mo,size(Fock_matrix_mo,1),size(Fock_matrix_mo,2),mo_label)
 
-  call write_double(output_hartree_fock, HF_energy, 'Hartree-Fock energy')
+  call write_double(output_hartree_fock, E_min, 'Hartree-Fock energy')
-  call ezfio_set_hartree_fock_energy(HF_energy)
+  call ezfio_set_hartree_fock_energy(E_min)
 
   call write_time(output_hartree_fock)
 

src/Hartree_Fock/hartree_fock.ezfio_config

@@ -2,4 +2,5 @@ hartree_fock
   thresh_scf       double precision
   n_it_scf_max     integer
   energy           double precision
+  guess            character*(32)
 

src/Hartree_Fock/huckel.irp.f

@@ -0,0 +1,34 @@
+subroutine huckel_guess
+  implicit none
+  BEGIN_DOC
+! Build the MOs using the extended Huckel model
+  END_DOC
+  integer                        :: i,j
+  double precision               :: tmp_matrix(ao_num_align,ao_num),accu
+  double precision               :: c
+  character*(64)                 :: label
+
+  mo_coef = ao_ortho_lowdin_coef
+  TOUCH mo_coef
+  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)
+  TOUCH mo_coef
+
+  c = 0.5d0 * 1.75d0
+  do j=1,ao_num
+    do i=1,ao_num
+      if (i/=j) then
+        Fock_matrix_ao(i,j) = c*ao_overlap(i,j)*(ao_mono_elec_integral(i,i) + &
+           ao_mono_elec_integral(j,j))
+      else
+        Fock_matrix_ao(i,j) = Fock_matrix_alpha_ao(i,j)
+      endif
+    enddo
+  enddo
+  TOUCH Fock_matrix_ao
+  mo_coef = eigenvectors_fock_matrix_mo
+  SOFT_TOUCH mo_coef
+  call save_mos
+
+end

src/Hartree_Fock/options.irp.f

@@ -18,6 +18,15 @@ T.set_ezfio_name( "n_it_scf_max" )
 T.set_output    ( "output_Hartree_Fock" )
 print T
 
+T = EZFIO_Provider()
+T.set_type      ( "character*(32)" )
+T.set_name      ( "mo_guess_type" )
+T.set_doc       ( "Initial MO guess. Can be [ Huckel | HCore ]" )
+T.set_ezfio_dir ( "Hartree_Fock" )
+T.set_ezfio_name( "guess" )
+T.set_output    ( "output_Hartree_Fock" )
+print T
+
 
 END_SHELL
 

src/MOGuess/H_CORE_guess.irp.f

@@ -4,7 +4,8 @@ program H_CORE_guess
   mo_coef = ao_ortho_lowdin_coef
   TOUCH mo_coef
   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)
+  call mo_as_eigvectors_of_mo_matrix(mo_mono_elec_integral,          &
+      size(mo_mono_elec_integral,1),size(mo_mono_elec_integral,2),label)
   print *,  'save mos'
   call save_mos