LCPQ/QCaml
10
1
Fork 0
mirror of https://gitlab.com/scemama/QCaml.git synced 2024-11-09 07:33:40 +01:00
Code Releases Activity

Changed sign of exchange matrix

Browse Source
This commit is contained in:
Anthony Scemama 2020-04-20 02:06:49 +02:00
parent d7be3ac86e
commit 0094b36c86
2 changed files with 14 additions and 14 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

24
SCF/Fock.ml
View File

@ -53,11 +53,11 @@ let make_rhf ~density ?(threshold=Constants.epsilon) ao_basis =
in
if (integral <> 0.) then begin
m_Jnu.(mu-1) <- m_Jnu.(mu-1) +. pJ *. integral;
m_Ksigma.(mu-1) <- m_Ksigma.(mu-1) -. pK *. integral
m_Ksigma.(mu-1) <- m_Ksigma.(mu-1) +. pK *. integral
end
done;
for mu = nu+1 to sigma do
m_Ksigma.(mu-1) <- m_Ksigma.(mu-1) -. pK *.
m_Ksigma.(mu-1) <- m_Ksigma.(mu-1) +. pK *.
ERI.get_phys m_G mu lambda nu sigma
done
end
@ -69,7 +69,7 @@ let make_rhf ~density ?(threshold=Constants.epsilon) ao_basis =
in
if (integral <> 0.) then begin
m_Jnu.(mu-1) <- m_Jnu.(mu-1) +. pJ *. integral;
m_Ksigma.(mu-1) <- m_Ksigma.(mu-1) -. pK *. integral
m_Ksigma.(mu-1) <- m_Ksigma.(mu-1) +. pK *. integral
end
done;
for mu = sigma+1 to nu do
@ -80,7 +80,7 @@ let make_rhf ~density ?(threshold=Constants.epsilon) ao_basis =
| (false, true , _) ->
for mu = 1 to sigma do
m_Ksigma.(mu-1) <-
m_Ksigma.(mu-1) -. pK *. ERI.get_phys m_G mu lambda nu sigma
m_Ksigma.(mu-1) +. pK *. ERI.get_phys m_G mu lambda nu sigma
done
| (true , false, _) ->
for mu = 1 to nu do
@ -103,7 +103,7 @@ let make_rhf ~density ?(threshold=Constants.epsilon) ao_basis =
let m_J = Mat.of_array m_J
and m_K = Mat.of_array m_K
in
{ fock = Mat.add m_Hc (Mat.add m_J m_K) ;
{ fock = Mat.add m_Hc (Mat.sub m_J m_K) ;
core = m_Hc ; coulomb = m_J ; exchange = m_K }
@ -141,11 +141,11 @@ let make_uhf ~density_same ~density_other ?(threshold=Constants.epsilon) ao_basi
in
if (integral <> 0.) then begin
m_Jnu.(mu-1) <- m_Jnu.(mu-1) +. pJ *. integral;
m_Ksigma.(mu-1) <- m_Ksigma.(mu-1) -. pK *. integral
m_Ksigma.(mu-1) <- m_Ksigma.(mu-1) +. pK *. integral
end
done;
for mu = nu+1 to sigma do
m_Ksigma.(mu-1) <- m_Ksigma.(mu-1) -. pK *.
m_Ksigma.(mu-1) <- m_Ksigma.(mu-1) +. pK *.
ERI.get_phys m_G mu lambda nu sigma
done
end
@ -157,7 +157,7 @@ let make_uhf ~density_same ~density_other ?(threshold=Constants.epsilon) ao_basi
in
if (integral <> 0.) then begin
m_Jnu.(mu-1) <- m_Jnu.(mu-1) +. pJ *. integral;
m_Ksigma.(mu-1) <- m_Ksigma.(mu-1) -. pK *. integral
m_Ksigma.(mu-1) <- m_Ksigma.(mu-1) +. pK *. integral
end
done;
for mu = sigma+1 to nu do
@ -168,7 +168,7 @@ let make_uhf ~density_same ~density_other ?(threshold=Constants.epsilon) ao_basi
| (false, true , _) ->
for mu = 1 to sigma do
m_Ksigma.(mu-1) <-
m_Ksigma.(mu-1) -. pK *. ERI.get_phys m_G mu lambda nu sigma
m_Ksigma.(mu-1) +. pK *. ERI.get_phys m_G mu lambda nu sigma
done
| (true , false, _) ->
for mu = 1 to nu do
@ -191,7 +191,7 @@ let make_uhf ~density_same ~density_other ?(threshold=Constants.epsilon) ao_basi
let m_J = Mat.of_array m_J
and m_K = Mat.of_array m_K
in
{ fock = Mat.add m_Hc (Mat.add m_J m_K) ;
{ fock = Mat.add m_Hc (Mat.sub m_J m_K) ;
core = m_Hc ; coulomb = m_J ; exchange = m_K }
@ -207,7 +207,7 @@ let op ~f f1 f2 =
and m_K = f f1.exchange f2.exchange
in
{
fock = Mat.add m_Hc (Mat.add m_J m_K);
fock = Mat.add m_Hc (Mat.sub m_J m_K);
core = m_Hc;
coulomb = m_J;
exchange = m_K;
@ -226,7 +226,7 @@ let scale alpha f1 =
Mat.scal alpha m_J;
Mat.scal alpha m_K;
{
fock = Mat.add m_Hc (Mat.add m_J m_K);
fock = Mat.add m_Hc (Mat.sub m_J m_K);
core = m_Hc;
coulomb = m_J;
exchange = m_K;

4
SCF/HartreeFock.ml
View File

@ -171,7 +171,7 @@ let exchange_energy t =
| RHF -> let m_P = of_some data.density in
let fock = of_some data.fock in
let m_K = Fock.exchange fock in
0.5 *. Mat.gemm_trace m_P m_K
-. 0.5 *. Mat.gemm_trace m_P m_K
| ROHF -> let m_P_a = of_some data.density_a in
let m_P_b = of_some data.density_b in
@ -179,7 +179,7 @@ let exchange_energy t =
let fock_b = of_some data.fock_b in
let m_K_a = Fock.exchange fock_a in
let m_K_b = Fock.exchange fock_b in
0.5 *. ( (Mat.gemm_trace m_P_a m_K_a) +. (Mat.gemm_trace m_P_b m_K_b) )
-. 0.5 *. ( (Mat.gemm_trace m_P_a m_K_a) +. (Mat.gemm_trace m_P_b m_K_b) )
| _ -> failwith "Not implemented"