From 0094b36c86993685e79ce6bbe2e97c49a0836b78 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 20 Apr 2020 02:06:49 +0200 Subject: [PATCH] Changed sign of exchange matrix --- SCF/Fock.ml | 24 ++++++++++++------------ SCF/HartreeFock.ml | 4 ++-- 2 files changed, 14 insertions(+), 14 deletions(-) diff --git a/SCF/Fock.ml b/SCF/Fock.ml index e8986fd..2508f43 100644 --- a/SCF/Fock.ml +++ b/SCF/Fock.ml @@ -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; diff --git a/SCF/HartreeFock.ml b/SCF/HartreeFock.ml index 0bc46ed..e5bb4fc 100644 --- a/SCF/HartreeFock.ml +++ b/SCF/HartreeFock.ml @@ -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"