mirror of https://gitlab.com/scemama/QCaml.git
Nuclear repulsion
This commit is contained in:
parent
34edda6317
commit
0455da3b46
|
|
@ -0,0 +1,7 @@
|
|||
PKG str unix bigarray lacaml
|
||||
S .
|
||||
S Nuclei
|
||||
S Utils
|
||||
S Basis
|
||||
S HartreeFock
|
||||
B _build/**
|
||||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
2
Makefile
2
Makefile
|
|
@ -1,6 +1,6 @@
|
|||
.NOPARALLEL:
|
||||
|
||||
INCLUDE_DIRS=Nuclei,Utils,Basis
|
||||
INCLUDE_DIRS=Nuclei,Utils,Basis,HartreeFock
|
||||
LIBS=
|
||||
PKGS=
|
||||
OCAMLCFLAGS="-g -warn-error A"
|
||||
|
|
|
|||
|
|
@ -1,4 +1,4 @@
|
|||
type t = float
|
||||
type t = private float
|
||||
|
||||
(** Float conversion functions *)
|
||||
val to_float : t -> float
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
||||
|
|
|
|||
|
|
@ -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;
|
||||
}
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
Loading…
Reference in New Issue