10
1
mirror of https://gitlab.com/scemama/QCaml.git synced 2024-06-14 09:15:19 +02:00
QCaml/SCF/Fock.ml

66 lines
1.5 KiB
OCaml
Raw Normal View History

2018-02-23 18:44:31 +01:00
open Lacaml.D
open Simulation
2018-03-13 18:24:00 +01:00
open Constants
2018-02-23 18:44:31 +01:00
2018-03-26 16:47:08 +02:00
2018-05-30 18:07:05 +02:00
type t =
{
fock : Mat.t ;
core : Mat.t ;
coulomb : Mat.t ;
exchange : Mat.t ;
}
2018-02-23 18:44:31 +01:00
let make ~density simulation =
let m_P = density
and m_T = Lazy.force simulation.kin_ints
and m_V = Lazy.force simulation.eN_ints
and m_G = Lazy.force simulation.ee_ints
in
let nBas = Mat.dim1 m_T
in
2018-05-30 18:07:05 +02:00
let m_Hc = Mat.add m_T m_V
2018-05-31 18:50:34 +02:00
and m_J = Array.make_matrix nBas nBas 0.
and m_K = Array.make_matrix nBas nBas 0.
2018-05-30 18:07:05 +02:00
in
2018-02-23 18:44:31 +01:00
for sigma = 1 to nBas do
for nu = 1 to nBas do
2018-05-31 18:50:34 +02:00
let m_Jnu = m_J.(nu-1) in
2018-02-23 18:44:31 +01:00
for lambda = 1 to nBas do
let p = m_P.{lambda,sigma} in
2018-03-13 18:24:00 +01:00
if abs_float p > epsilon then
for mu = 1 to nu do
2018-05-31 18:50:34 +02:00
m_Jnu.(mu-1) <- m_Jnu.(mu-1) +. p *.
2018-05-30 18:07:05 +02:00
ERI.get_phys m_G mu lambda nu sigma
done
done
done
done;
for nu = 1 to nBas do
2018-05-31 18:50:34 +02:00
let m_Knu = m_K.(nu-1) in
2018-05-30 18:07:05 +02:00
for sigma = 1 to nBas do
for lambda = 1 to nBas do
let p = 0.5 *. m_P.{lambda,sigma} in
if abs_float p > epsilon then
for mu = 1 to nu do
2018-05-31 18:50:34 +02:00
m_Knu.(mu-1) <- m_Knu.(mu-1) -. p *.
2018-05-30 18:07:05 +02:00
ERI.get_phys m_G mu lambda sigma nu
2018-03-13 18:24:00 +01:00
done
2018-02-23 18:44:31 +01:00
done
done
done;
for nu = 1 to nBas do
2018-05-31 18:50:34 +02:00
for mu = 1 to nu-1 do
m_J.(mu-1).(nu-1) <- m_J.(nu-1).(mu-1);
m_K.(mu-1).(nu-1) <- m_K.(nu-1).(mu-1);
2018-02-23 18:44:31 +01:00
done
done;
2018-05-31 18:50:34 +02:00
let m_J = Mat.of_array m_J
and m_K = Mat.of_array m_K
in
2018-05-30 18:07:05 +02:00
{ fock = Mat.add m_Hc @@ Mat.add m_J m_K ;
core = m_Hc ; coulomb = m_J ; exchange = m_K }
2018-02-23 18:44:31 +01:00