diff --git a/data/ezfio_defaults b/data/ezfio_defaults index 9d20a5b7..5da4d813 100644 --- a/data/ezfio_defaults +++ b/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 diff --git a/ocaml/Input_hartree_fock.ml b/ocaml/Input_hartree_fock.ml index 845fb070..82110647 100644 --- a/ocaml/Input_hartree_fock.ml +++ b/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) - (Threshold.to_string b.thresh_scf) + (MO_guess.to_string b.guess) + (Strictly_positive_int.to_string b.n_it_scf_max) + (Threshold.to_string b.thresh_scf) |> Rst_string.of_string ;; diff --git a/ocaml/qp_create_ezfio_from_xyz.ml b/ocaml/qp_create_ezfio_from_xyz.ml index 1eb3d59b..dfc9b4bd 100644 --- a/ocaml/qp_create_ezfio_from_xyz.ml +++ b/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 = diff --git a/ocaml/qptypes_generator.ml b/ocaml/qptypes_generator.ml index e47213ff..f5a62356 100644 --- a/ocaml/qptypes_generator.ml +++ b/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 diff --git a/scripts/get_basis.sh b/scripts/get_basis.sh index d01324fe..d6182799 100755 --- a/scripts/get_basis.sh +++ b/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="" diff --git a/src/Hartree_Fock/Huckel_guess.irp.f b/src/Hartree_Fock/Huckel_guess.irp.f new file mode 100644 index 00000000..ed264005 --- /dev/null +++ b/src/Hartree_Fock/Huckel_guess.irp.f @@ -0,0 +1,6 @@ +program guess + implicit none + character*(64) :: label + call huckel_guess + +end diff --git a/src/Hartree_Fock/README.rst b/src/Hartree_Fock/README.rst index 8f468717..8060f092 100644 --- a/src/Hartree_Fock/README.rst +++ b/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 ============== diff --git a/src/Hartree_Fock/SCF.irp.f b/src/Hartree_Fock/SCF.irp.f index b9dcd80b..33e1ac6c 100644 --- a/src/Hartree_Fock/SCF.irp.f +++ b/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 - TOUCH mo_coef - mo_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),mo_label) - SOFT_TOUCH mo_coef mo_label + if (mo_guess_type == "HCore") then + mo_coef = ao_ortho_lowdin_coef + TOUCH mo_coef + mo_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),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 diff --git a/src/Hartree_Fock/damping_SCF.irp.f b/src/Hartree_Fock/damping_SCF.irp.f index 4c180c16..f0f59080 100644 --- a/src/Hartree_Fock/damping_SCF.irp.f +++ b/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 ezfio_set_hartree_fock_energy(HF_energy) + call write_double(output_hartree_fock, E_min, 'Hartree-Fock energy') + call ezfio_set_hartree_fock_energy(E_min) call write_time(output_hartree_fock) diff --git a/src/Hartree_Fock/hartree_fock.ezfio_config b/src/Hartree_Fock/hartree_fock.ezfio_config index 07030bc9..3909f652 100644 --- a/src/Hartree_Fock/hartree_fock.ezfio_config +++ b/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) diff --git a/src/Hartree_Fock/huckel.irp.f b/src/Hartree_Fock/huckel.irp.f new file mode 100644 index 00000000..4ea8d93f --- /dev/null +++ b/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 diff --git a/src/Hartree_Fock/options.irp.f b/src/Hartree_Fock/options.irp.f index 20b96b7c..22f86d6c 100644 --- a/src/Hartree_Fock/options.irp.f +++ b/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 diff --git a/src/MOGuess/H_CORE_guess.irp.f b/src/MOGuess/H_CORE_guess.irp.f index 17b53aaa..2e10dd94 100644 --- a/src/MOGuess/H_CORE_guess.irp.f +++ b/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