diff --git a/.merlin b/.merlin new file mode 100644 index 0000000..05d34eb --- /dev/null +++ b/.merlin @@ -0,0 +1,7 @@ +PKG str unix bigarray lacaml +S . +S Nuclei +S Utils +S Basis +S HartreeFock +B _build/** diff --git a/Basis/TwoElectronRR.ml b/Basis/TwoElectronRR.ml index 2abedb3..ea5f70e 100644 --- a/Basis/TwoElectronRR.ml +++ b/Basis/TwoElectronRR.ml @@ -41,30 +41,9 @@ let rec hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d) | _ -> Z in - (* - if debug then begin - Printf.printf "\n---- %d %d %d %d ----\n" totAngMom_a totAngMom_b totAngMom_c totAngMom_d; - let (x,y,z) = angMom_a in Printf.printf "%d %d %d\n" x y z; - let (x,y,z) = angMom_b in Printf.printf "%d %d %d\n" x y z; - let (x,y,z) = angMom_c in Printf.printf "%d %d %d\n" x y z; - let (x,y,z) = angMom_d in Printf.printf "%d %d %d\n" x y z; - Printf.printf "%f %f %f %f\n%f %f %f\n%f %f %f\n%f %f %f\n" expo_b expo_d - expo_inv_p expo_inv_q - (get X center_ab) (get Y center_ab) (get Z center_ab) - (get X center_cd) (get Y center_cd) (get Z center_cd) - (get X center_pq) (get Y center_pq) (get Z center_pq) - end; - *) (** Vertical recurrence relations *) let rec vrr0 angMom_a = - (* - if debug then - begin - let (x,y,z) = angMom_a in - Printf.printf "vrr0: %d : %d %d %d\n" angMom_a.tot x y z - end; - *) match angMom_a.tot with | 0 -> zero_m_array @@ -78,59 +57,45 @@ let rec hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d) let am = Powers.decr xyz angMom_a in let amxyz = Powers.get xyz am in - if amxyz < 0 then empty else - let f1 = expo_inv_p *. (Coordinate.get xyz center_pq) - and f2 = expo_b *. expo_inv_p *. (Coordinate.get xyz center_ab) - in - let result = Array.create_float maxsze in - if amxyz < 1 then - begin - let v1 = - vrr0 am - in - for m=0 to maxm-1 do - result.(m) <- f1 *. v1.(m+1) -. f2 *. v1.(m) - done; - result.(maxm) <- -. f2 *. v1.(maxm) - end - else - begin - let amm = Powers.decr xyz am in - let v3 = vrr0 amm in - let v1 = vrr0 am in - let f3 = (float_of_int amxyz) *. expo_inv_p *. 0.5 in - for m=0 to maxm-1 do - result.(m) <- f1 *. v1.(m+1) -. f2 *. v1.(m) - +. f3 *. (v3.(m) +. expo_inv_p *. v3.(m+1)) - done; - result.(maxm) <- f3 *. v3.(maxm) - end; - result + if amxyz >= 0 then + let f1 = expo_inv_p *. (Coordinate.get xyz center_pq) + and f2 = expo_b *. expo_inv_p *. (Coordinate.get xyz center_ab) + in + let result = Array.create_float maxsze in + if amxyz < 1 then + begin + let v1 = + vrr0 am + in + for m=0 to maxm-1 do + result.(m) <- f1 *. v1.(m+1) -. f2 *. v1.(m) + done; + result.(maxm) <- -. f2 *. v1.(maxm) + end + else + begin + let amm = Powers.decr xyz am in + let v3 = vrr0 amm in + let v1 = vrr0 am in + let f3 = (float_of_int amxyz) *. expo_inv_p *. 0.5 in + for m=0 to maxm-1 do + result.(m) <- f1 *. v1.(m+1) -. f2 *. v1.(m) + +. f3 *. (v3.(m) +. expo_inv_p *. v3.(m+1)) + done; + result.(maxm) <- f3 *. v3.(maxm) + end; + result + else + empty in Zmap.add map_1d key result; result and vrr angMom_a angMom_c = - (* - if debug then - begin - let angMom_ax, angMom_ay, angMom_az = angMom_a in - let angMom_cx, angMom_cy, angMom_cz = angMom_c in - Printf.printf "vrr : %d %d : %d %d %d %d %d %d\n" angMom_a.tot angMom_c.tot - angMom_ax angMom_ay angMom_az angMom_cx angMom_cy angMom_cz - end; - *) - match (angMom_a.tot, angMom_c.tot) with - | (i,0) -> if (i>0) then - vrr0 angMom_a - (* - OneElectronRR.hvrr_one_e (angMom_a, angMom_b) (angMom_a.tot, angMom_b.tot) - (maxm, zero_m_array) (expo_b) (expo_inv_p) (center_ab, center_pq, center_ab) - map_1d - *) - else zero_m_array + | (i,0) -> if (i>0) then vrr0 angMom_a + else zero_m_array | (_,_) -> let key = Zkey.of_powers (Zkey.Six (angMom_a, angMom_c)) in @@ -142,58 +107,59 @@ let rec hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d) let cmxyz = Powers.get xyz cm in let axyz = Powers.get xyz angMom_a in - if cmxyz < 0 then empty - else - let f1 = - -. expo_d *. expo_inv_q *. (Coordinate.get xyz center_cd) - and f2 = - expo_inv_q *. (Coordinate.get xyz center_pq) - in - let result = Array.make maxsze 0. in - if ( (abs_float f1 > cutoff) || (abs_float f2 > cutoff) ) then - begin - let v1 = - vrr angMom_a cm - in - for m=0 to maxm-1 do - result.(m) <- f1 *. v1.(m) -. f2 *. v1.(m+1) ; - done; - result.(maxm) <- f1 *. v1.(maxm) ; - end; - if cmxyz > 0 then - begin - let f3 = - (float_of_int cmxyz) *. expo_inv_q *. 0.5 - in - if (abs_float f3 > cutoff) || - (abs_float (f3 *. expo_inv_q) > cutoff) then + if cmxyz >= 0 then + let f1 = + -. expo_d *. expo_inv_q *. (Coordinate.get xyz center_cd) + and f2 = + expo_inv_q *. (Coordinate.get xyz center_pq) + in + let result = Array.make maxsze 0. in + if ( (abs_float f1 > cutoff) || (abs_float f2 > cutoff) ) then begin - let v3 = - let cmm = Powers.decr xyz cm in - vrr angMom_a cmm + let v1 = + vrr angMom_a cm in for m=0 to maxm-1 do - result.(m) <- result.(m) +. - f3 *. (v3.(m) +. expo_inv_q *. v3.(m+1)) + result.(m) <- f1 *. v1.(m) -. f2 *. v1.(m+1) ; done; - result.(maxm) <- result.(maxm) +. f3 *. v3.(maxm) - end - end; - if (axyz > 0) && (cmxyz >= 0) then - begin - let am = Powers.decr xyz angMom_a in - let f5 = - (float_of_int axyz) *. expo_inv_p *. expo_inv_q *. 0.5 - in - if (abs_float f5 > cutoff) then - let v5 = - vrr am cm - in - for m=0 to maxm-1 do - result.(m) <- result.(m) -. f5 *. v5.(m+1) - done - end; - result + result.(maxm) <- f1 *. v1.(maxm) ; + end; + if cmxyz > 0 then + begin + let f3 = + (float_of_int cmxyz) *. expo_inv_q *. 0.5 + in + if (abs_float f3 > cutoff) || + (abs_float (f3 *. expo_inv_q) > cutoff) then + begin + let v3 = + let cmm = Powers.decr xyz cm in + vrr angMom_a cmm + in + for m=0 to maxm-1 do + result.(m) <- result.(m) +. + f3 *. (v3.(m) +. expo_inv_q *. v3.(m+1)) + done; + result.(maxm) <- result.(maxm) +. f3 *. v3.(maxm) + end + end; + if axyz > 0 && cmxyz >= 0 then + begin + let am = Powers.decr xyz angMom_a in + let f5 = + (float_of_int axyz) *. expo_inv_p *. expo_inv_q *. 0.5 + in + if (abs_float f5 > cutoff) then + let v5 = + vrr am cm + in + for m=0 to maxm-1 do + result.(m) <- result.(m) -. f5 *. v5.(m+1) + done + end; + result + else + empty in Zmap.add map_2d key result; result @@ -204,74 +170,47 @@ let rec hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d) (** Horizontal recurrence relations *) and hrr0 angMom_a angMom_b angMom_c = - (* - if debug then - begin - let angMom_ax, angMom_ay, angMom_az = angMom_a - and angMom_bx, angMom_by, angMom_bz = angMom_b - and angMom_cx, angMom_cy, angMom_cz = angMom_c in - Printf.printf "hrr0: %d %d %d : %d %d %d %d %d %d %d %d %d\n" - angMom_ax angMom_ay angMom_az - angMom_bx angMom_by angMom_bz - angMom_cx angMom_cy angMom_cz - end; - *) - match angMom_b.tot with | 0 -> (vrr angMom_a angMom_c).(0) | 1 -> let xyz = get_xyz angMom_b in let ap = Powers.incr xyz angMom_a in - let v1 = - vrr ap angMom_c - in - let f2 = - (Coordinate.get xyz center_ab) - in + let v1 = vrr ap angMom_c in + let f2 = Coordinate.get xyz center_ab in if (abs_float f2 < cutoff) then v1.(0) else - let v2 = - vrr angMom_a angMom_c - in - v1.(0) +. f2 *. v2.(0) + let v2 = vrr angMom_a angMom_c in + v1.(0) +. f2 *. v2.(0) | _ -> let xyz = get_xyz angMom_b in let bxyz = Powers.get xyz angMom_b in - if (bxyz < 1) then 0. else - let ap = Powers.incr xyz angMom_a in - let bm = Powers.decr xyz angMom_b in - let h1 = - hrr0 ap bm angMom_c - in - let f2 = - (Coordinate.get xyz center_ab) - in - if (abs_float f2 < cutoff) then h1 else - let h2 = - hrr0 angMom_a bm angMom_c - in - h1 +. f2 *. h2 + if bxyz > 0 then + let ap = Powers.incr xyz angMom_a in + let bm = Powers.decr xyz angMom_b in + let h1 = hrr0 ap bm angMom_c in + let f2 = Coordinate.get xyz center_ab in + if abs_float f2 < cutoff then h1 else + let h2 = hrr0 angMom_a bm angMom_c in + h1 +. f2 *. h2 + else 0. and hrr angMom_a angMom_b angMom_c angMom_d = match (angMom_b.tot, angMom_d.tot) with - | (_,0) -> if (angMom_b.tot = 0) then - (vrr angMom_a angMom_c).(0) + | (_,0) -> + if (angMom_b.tot = 0) then + (vrr angMom_a angMom_c).(0) else - hrr0 angMom_a angMom_b angMom_c + hrr0 angMom_a angMom_b angMom_c | (_,_) -> let xyz = get_xyz angMom_d in let cp = Powers.incr xyz angMom_c in let dm = Powers.decr xyz angMom_d in - let h1 = - hrr angMom_a angMom_b cp dm - in + let h1 = hrr angMom_a angMom_b cp dm in let f2 = Coordinate.get xyz center_cd in - if (abs_float f2 < cutoff) then h1 else - let h2 = - hrr angMom_a angMom_b angMom_c dm - in - h1 +. f2 *. h2 + if abs_float f2 < cutoff then h1 else + let h2 = hrr angMom_a angMom_b angMom_c dm in + h1 +. f2 *. h2 in hrr angMom_a angMom_b angMom_c angMom_d diff --git a/Basis/TwoElectronRRVectorized.ml b/Basis/TwoElectronRRVectorized.ml index 066316d..e6e154c 100644 --- a/Basis/TwoElectronRRVectorized.ml +++ b/Basis/TwoElectronRRVectorized.ml @@ -159,17 +159,10 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) | None -> None | Some v1 -> begin - let result = Array.make_matrix np nq 0. in - for l=0 to np-1 do - for k=0 to nq-1 do - result.(l).(k) <- v1.(l).(k) *. f1.(k) - done - done; - Some ( - Array.init np (fun l -> - let v1_l = v1.(l) in - Array.init nq (fun k -> v1_l.(k) *. f1.(k)) - )) + Some (Array.init np (fun l -> + let v1_l = v1.(l) in + Array.mapi (fun k f1k -> v1_l.(k) *. f1k) f1 + ) ) end else None in diff --git a/Makefile b/Makefile index b9fec0c..b7b083a 100644 --- a/Makefile +++ b/Makefile @@ -1,6 +1,6 @@ .NOPARALLEL: -INCLUDE_DIRS=Nuclei,Utils,Basis +INCLUDE_DIRS=Nuclei,Utils,Basis,HartreeFock LIBS= PKGS= OCAMLCFLAGS="-g -warn-error A" diff --git a/Nuclei/Charge.mli b/Nuclei/Charge.mli index 83b09ad..a54e81a 100644 --- a/Nuclei/Charge.mli +++ b/Nuclei/Charge.mli @@ -1,4 +1,4 @@ -type t = float +type t = private float (** Float conversion functions *) val to_float : t -> float diff --git a/Nuclei/Nuclei.ml b/Nuclei/Nuclei.ml index cab3cad..51fea16 100644 --- a/Nuclei/Nuclei.ml +++ b/Nuclei/Nuclei.ml @@ -1,3 +1,4 @@ +open Util open Xyz_ast type t = (Element.t * Coordinate.t) array @@ -70,3 +71,18 @@ let of_filename ~filename = of_xyz_file filename +let repulsion ~nuclei = + let get_charge e = + Element.to_charge e + |> Charge.to_float + in + Array.fold_left ( fun accu (e1, coord1) -> + accu +. + Array.fold_left (fun accu (e2, coord2) -> + let r = Coordinate.(norm (coord1 |- coord2)) in + if r > 0. then + accu +. 0.5 *. (get_charge e2) *. (get_charge e1) /. r + else accu + ) 0. nuclei + ) 0. nuclei + diff --git a/Simulation.ml b/Simulation.ml index d979da5..f87723c 100644 --- a/Simulation.ml +++ b/Simulation.ml @@ -5,6 +5,7 @@ type t = { eN_ints : NucInt.t lazy_t; kin_ints : KinInt.t lazy_t; ee_ints : ERI.t lazy_t; + nuclear_repulsion : float; } let make ~nuclei ~basis = @@ -14,6 +15,7 @@ let make ~nuclei ~basis = eN_ints = lazy (NucInt.of_basis_nuclei basis nuclei); kin_ints = lazy (KinInt.of_basis basis); ee_ints = lazy (ERI.of_basis basis); + nuclear_repulsion = Nuclei.repulsion nuclei; } diff --git a/Utils/Util.ml b/Utils/Util.ml index ed7782c..b04c0f2 100644 --- a/Utils/Util.ml +++ b/Utils/Util.ml @@ -144,3 +144,11 @@ let boys_function ~maxm t = result end + + +let array_sum a = + Array.fold_left ( +. ) 0. a + +let array_product a = + Array.fold_left ( *. ) 0. a + diff --git a/run_integrals.ml b/run_integrals.ml index c165659..33efa5d 100644 --- a/run_integrals.ml +++ b/run_integrals.ml @@ -36,6 +36,8 @@ let run ~out = in print_endline @@ Nuclei.to_string s.Simulation.nuclei; + print_endline "Nuclear repulsion : "; + print_float s.Simulation.nuclear_repulsion; print_newline (); print_endline @@ Basis.to_string s.Simulation.basis; let overlap = Lazy.force s.Simulation.overlap in