From 4b0021d47f4bddd6049b350c4ffb54e4976eea73 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 31 May 2017 01:12:06 +0200 Subject: [PATCH] Fixed CASSD --- ocaml/Gamess.ml | 457 ++++++++++++++++++++++++ plugins/Selectors_CASSD/selectors.irp.f | 10 +- 2 files changed, 462 insertions(+), 5 deletions(-) create mode 100644 ocaml/Gamess.ml diff --git a/ocaml/Gamess.ml b/ocaml/Gamess.ml new file mode 100644 index 00000000..3b21df51 --- /dev/null +++ b/ocaml/Gamess.ml @@ -0,0 +1,457 @@ +(** CONTRL *) +type scftyp_t = RHF | ROHF | MCSCF | NONE +let string_of_scftyp = function +| RHF -> "RHF" +| ROHF -> "ROHF" +| MCSCF -> "MCSCF" +| NONE -> "NONE" + +type contrl = +{ scftyp: scftyp_t ; + maxit: int; + ispher: int; + icharg: int; + mult: int; + mplevl: int; +} + +let string_of_contrl c = + Printf.sprintf " $CONTRL + EXETYP=RUN COORD=UNIQUE UNITS=ANGS + RUNTYP=ENERGY SCFTYP=%s CITYP=NONE + MAXIT=%d + ISPHER=%d + MULT=%d + ICHARG=%d + MPLEVL=%d + $END" + (string_of_scftyp c.scftyp) + c.maxit c.ispher c.mult c.icharg c.mplevl + +let make_contrl ?(maxit=100) ?(ispher=1) ?(mplevl=0) ~mult ~charge scftyp = + { scftyp ; maxit ; ispher ; mult ; icharg=charge ; mplevl } + + +(** Vec *) +type vec_t = +| Canonical of string +| Natural of string + +let read_mos guide filename = + let text = + let ic = open_in filename in + let n = in_channel_length ic in + let s = Bytes.create n in + really_input ic s 0 n; + close_in ic; + s + in + + let re_vec = + Str.regexp " \\$VEC *\n" + and re_natural = + Str.regexp guide + and re_end = + Str.regexp " \\$END *\n" + and re_eol = + Str.regexp "\n" + in + let i = + Str.search_forward re_natural text 0 + in + let start = + Str.search_forward re_vec text i + in + let i = + Str.search_forward re_end text start + in + let finish = + Str.search_forward re_eol text i + in + String.sub text start (finish-start) + +let read_until_found f tries = + let result = + List.fold_left (fun accu x -> + match accu with + | Some mos -> Some mos + | None -> + begin + try + Some (read_mos x f) + with Not_found -> + None + end + ) None tries + in + match result with + | Some mos -> mos + | None -> raise Not_found + +let read_natural_mos f = + let tries = [ + "--- NATURAL ORBITALS OF MCSCF ---" ; + "MP2 NATURAL ORBITALS" ] + in + read_until_found f tries + +let read_canonical_mos f = + let tries = [ + "--- OPTIMIZED MCSCF MO-S ---" ; + "--- CLOSED SHELL ORBITALS ---" ; + "--- OPEN SHELL ORBITALS ---" + ] + in + read_until_found f tries + +let string_of_vec = function +| Natural filename -> read_natural_mos filename +| Canonical filename -> read_canonical_mos filename + +(** GUESS *) +type guess_t = +| Huckel +| Hcore +| Canonical of (int*string) +| Natural of (int*string) + +let guess_of_string s = + match String.lowercase s with + | "huckel" -> Huckel + | "hcore" -> Hcore + | _ -> raise (Invalid_argument "Bad MO guess") + +let string_of_guess g = + [ + " $GUESS\n" ; " GUESS=" ; + begin + match g with + | Hcore -> "HCORE\n" + | Huckel -> "HUCKEL\n" + | Canonical (norb,_) | Natural (norb,_) -> Printf.sprintf "MOREAD\n NORB=%d\n" norb + end + ; " $END" ; + match g with + | Hcore + | Huckel -> "" + | Natural (_,filename) -> "\n\n"^(string_of_vec (Natural filename)) + | Canonical (_,filename) ->"\n\n"^(string_of_vec (Canonical filename)) + ] |> String.concat "" + + +(** BASIS *) +let string_of_basis = + Printf.sprintf " $BASIS + GBASIS=%s + $END" + + +(** DATA *) +type coord_t = +| Atom of Element.t +| Diatomic_homo of (Element.t*float) +| Diatomic of (Element.t*Element.t*float) +| Xyz of (Element.t*float*float*float) list + + +type data_t = +{ sym: Sym.t ; + title: string; + xyz: string; + nucl_charge: int; +} + +let data_of_atom ele = + let atom = + Element.to_string ele + in + let charge = + Element.to_charge ele + |> Charge.to_int + in + { sym=Sym.D4h ; + title=Printf.sprintf "%s" atom ; + xyz=Printf.sprintf "%s %d.0 0. 0. 0." atom charge ; + nucl_charge = charge + } + +let data_of_diatomic_homo ele r = + assert (r > 0.); + let atom = + Element.to_string ele + in + let charge = + Element.to_charge ele + |> Charge.to_int + in + { sym=Sym.D4h ; + title=Printf.sprintf "%s2" atom ; + xyz=Printf.sprintf "%s %d.0 0. 0. %f" atom charge (-.r *. 0.5) ; + nucl_charge = 2*charge + } + +let data_of_diatomic ele1 ele2 r = + assert (r > 0.); + let atom1, atom2 = + Element.to_string ele1, + Element.to_string ele2 + in + let charge1, charge2 = + Charge.to_int @@ Element.to_charge ele1, + Charge.to_int @@ Element.to_charge ele2 + in + { sym=Sym.C4v ; + title=Printf.sprintf "%s%s" atom1 atom2 ; + xyz=Printf.sprintf "%s %d.0 0. 0. 0.\n%s %d.0 0. 0. %f" + atom1 charge1 atom2 charge2 r ; + nucl_charge = charge1 + charge2 + } + +let data_of_xyz l = + { sym = Sym.C1 ; + title = "..." ; + xyz = String.concat "\n" ( + List.map (fun (e,x,y,z) -> Printf.sprintf "%s %f %f %f %f" + (Element.to_string e) (Element.to_charge e) + x y z) l ) ; + nucl_charge = List.fold_left (fun accu (e,_,_,_) -> + accu + (int_of_float @@ Element.to_charge e) ) 0 l + } + +let make_data = function +| Atom ele -> data_of_atom ele +| Diatomic_homo (ele,r) -> data_of_diatomic_homo ele r +| Diatomic (ele1,ele2,r) -> data_of_diatomic ele1 ele2 r +| Xyz l -> data_of_xyz l + +let string_of_data d = + String.concat "\n" [ " $DATA" ; + d.title ; + Sym.to_data d.sym ; + ] ^ d.xyz ^ "\n $END" + + +(** GUGDM *) +type gugdm2_t = int + +let string_of_gugdm2 = function +| 1 -> "" +| i when i<1 -> raise (Invalid_argument "Nstates must be > 0") +| i -> + let s = + Array.make i "1." + |> Array.to_list + |> String.concat "," + in + Printf.sprintf " + $GUGDM2 + WSTATE(1)=%s + $END +" s + + +type gugdia_t = +{ nstate : int ; + itermx : int ; +} + +let string_of_gugdia g = + Printf.sprintf " + $GUGDIA + PRTTOL=0.0001 + NSTATE=%d + ITERMX=%d + $END +" g.nstate g.itermx + + +let make_gugdia ?(itermx=500) nstate = + assert (nstate > 0); + assert (itermx > 1); + { nstate ; itermx } + + +(** MCSCF *) +type mcscf_t = FULLNR | SOSCF | FOCAS + +let string_of_mcscf m = + " $MCSCF\n" ^ + begin + match m with + | FOCAS -> " FOCAS=.T. SOSCF=.F. FULLNR=.F." + | SOSCF -> " FOCAS=.F. SOSCF=.T. FULLNR=.F." + | FULLNR -> " FOCAS=.F. SOSCF=.F. FULLNR=.T." + end ^ " + CISTEP=GUGA EKT=.F. QUAD=.F. JACOBI=.f. + MAXIT=1000 + $END" + + +type drt_t = +{ nmcc: int ; + ndoc: int ; + nalp: int ; + nval: int ; + istsym: int; +} + + +let make_drt ?(istsym=1) n_elec_alpha n_elec_beta n_e n_act = + let n_elec_tot = + n_elec_alpha + n_elec_beta + in + let nmcc = + (n_elec_tot - n_e)/2 + in + let ndoc = + n_elec_beta - nmcc + in + let nalp = + (n_elec_alpha - nmcc - ndoc) + in + let nval = + n_act - ndoc - nalp + in + { nmcc ; ndoc ; nalp ; nval ; istsym } + +let string_of_drt drt sym = + Printf.sprintf " $DRT + NMCC=%d + NDOC=%d + NALP=%d + NVAL=%d + NEXT=0 + ISTSYM=%d + FORS=.TRUE. + GROUP=C1 + MXNINT= 600000 + NPRT=2 + $END" + drt.nmcc drt.ndoc drt.nalp drt.nval drt.istsym + +(** MP2 *) +let string_of_mp2 = " $MP2 + MP2PRP=.TRUE. + $END" + + +(** Computation *) +type computation = HF | MP2 | CAS of (int*int) + +type system = +{ mult: int ; charge: int ; basis: string ; coord: coord_t } + +let n_elec system = + let data = + make_data system.coord + in + data.nucl_charge - system.charge + +let n_elec_alpha_beta system = + let n = + n_elec system + and m = + system.mult + in + let alpha = + (n+m-1)/2 + in + let beta = + n - alpha + in + (alpha, beta) + + +let create_single_det_input ~mp2 ~guess ?(vecfile="") s = + let scftyp = + match s.mult with + | 1 -> RHF + | _ -> ROHF + and mult = s.mult + and charge = s.charge + and n_elec_alpha, _ = + n_elec_alpha_beta s + and mplevl = + if mp2 then 2 else 0 + in + [ + make_contrl ~mult ~charge ~mplevl scftyp + |> string_of_contrl + ; + begin + match vecfile with + | "" -> string_of_guess guess + | vecfile -> string_of_guess (Canonical (n_elec_alpha, vecfile)) + end + ; + string_of_basis s.basis + ; + if mp2 then + string_of_mp2 + else + "" + ; + make_data s.coord + |> string_of_data + ] |> String.concat "\n\n" + + +let create_hf_input ~guess = + create_single_det_input ~mp2:false ~guess + +let create_mp2_input ~guess = + create_single_det_input ~mp2:true ~guess + + +let create_cas_input ?(vecfile="") ~guess ~nstate s n_e n_a = + let scftyp = MCSCF + and mult = s.mult + and charge = s.charge + in + let n_elec_alpha, n_elec_beta = + n_elec_alpha_beta s + in + let drt = + make_drt n_elec_alpha n_elec_beta n_e n_a + in + let data = + make_data s.coord + in + [ + make_contrl ~mult ~charge scftyp + |> string_of_contrl + ; + begin + match vecfile with + | "" -> string_of_guess guess + | vecfile -> + let norb = + drt.nmcc + drt.ndoc + drt.nval + drt.nalp + in + try + string_of_guess (Natural (norb, vecfile)) + with Not_found -> + string_of_guess (Canonical (norb, vecfile)) + end + ; + string_of_basis s.basis + ; + string_of_mcscf FULLNR + ; + string_of_drt drt data.sym + ; + make_gugdia nstate + |> string_of_gugdia + ; + string_of_gugdm2 nstate + ; + string_of_data data + ] |> String.concat "\n\n" + + +let create_input ?(vecfile="") ?(guess=Huckel) ~system ~nstate = function +| HF -> create_hf_input ~vecfile ~guess system +| MP2 -> create_mp2_input ~vecfile ~guess system +| CAS (n_e,n_a) -> create_cas_input ~vecfile ~nstate ~guess system n_e n_a + + diff --git a/plugins/Selectors_CASSD/selectors.irp.f b/plugins/Selectors_CASSD/selectors.irp.f index 693f56e4..4023d8d6 100644 --- a/plugins/Selectors_CASSD/selectors.irp.f +++ b/plugins/Selectors_CASSD/selectors.irp.f @@ -23,15 +23,15 @@ END_PROVIDER integer :: i, k, l, m logical :: good - do i=1,N_det_generators + do i=1,N_det do k=1,N_int - psi_selectors(k,1,i) = psi_det_generators(k,1,i) - psi_selectors(k,2,i) = psi_det_generators(k,2,i) + psi_selectors(k,1,i) = psi_det_sorted(k,1,i) + psi_selectors(k,2,i) = psi_det_sorted(k,2,i) enddo enddo do k=1,N_states - do i=1,N_det_selectors - psi_selectors_coef(i,k) = psi_coef_generators(i,k) + do i=1,N_det + psi_selectors_coef(i,k) = psi_coef_sorted(i,k) enddo enddo