Add nuclear repulsion

devel/mpn/.gitignore

@@ -0,0 +1,59 @@
+IRPF90_temp/
+IRPF90_man/
+build.ninja
+irpf90.make
+ezfio_interface.irp.f
+irpf90_entities
+tags
+Makefile
+ao_basis
+ao_one_e_ints
+ao_two_e_erf_ints
+ao_two_e_ints
+aux_quantities
+becke_numerical_grid
+bitmask
+cis
+cisd
+cipsi
+davidson
+davidson_dressed
+davidson_undressed
+density_for_dft
+determinants
+dft_keywords
+dft_utils_in_r
+dft_utils_one_e
+dft_utils_two_body
+dressing
+dummy
+electrons
+ezfio_files
+fci
+generators_cas
+generators_full
+hartree_fock
+iterations
+kohn_sham
+kohn_sham_rs
+mo_basis
+mo_guess
+mo_one_e_ints
+mo_two_e_erf_ints
+mo_two_e_ints
+mpi
+mrpt_utils
+nuclei
+perturbation
+pseudo
+psiref_cas
+psiref_utils
+scf_utils
+selectors_cassd
+selectors_full
+selectors_utils
+single_ref_method
+slave
+tools
+utils
+zmq

devel/mpn/EZFIO.cfg

@@ -0,0 +1,12 @@
+[energy]
+type: double precision
+doc: Calculated Selected |FCI| energy
+interface: ezfio
+size: (determinants.n_states)
+
+[mp_order]
+type: integer
+doc: Max order of MPn
+interface: ezfio, provider, ocaml
+default: 4
+

devel/mpn/NEED

@@ -0,0 +1,3 @@
+davidson_undressed
+hartree_fock
+determinants

devel/mpn/README.rst

@@ -0,0 +1,4 @@
+===
+mpn
+===
+

devel/mpn/energies.irp.f

@@ -0,0 +1,23 @@
+BEGIN_PROVIDER [ double precision, energy_det_i, (N_det) ]
+ implicit none
+ BEGIN_DOC
+ ! Fock Energy of determinant |I> (sum of epsilon_i) 
+ END_DOC
+ integer :: i, k, n
+ integer :: list(elec_alpha_num)
+
+ do k=1,N_det
+   call bitstring_to_list(psi_det(1,1,k), list, n, N_int)
+   energy_det_i(k) = nuclear_repulsion
+   do i=1,n
+     energy_det_i(k) += fock_matrix_diag_mo(list(i))
+   enddo
+   call bitstring_to_list(psi_det(1,2,k), list, n, N_int)
+   do i=1,n
+     energy_det_i(k) += fock_matrix_diag_mo(list(i))
+   enddo
+ enddo
+END_PROVIDER
+
+
+

devel/mpn/mpn.irp.f

@@ -0,0 +1,100 @@
+program mpn
+  implicit none
+  BEGIN_DOC
+! TODO : Put the documentation of the program here
+  END_DOC
+  integer :: i, k, l
+  double precision, allocatable :: c_pert(:,:)
+  double precision, allocatable :: e_pert(:)
+  double precision, allocatable :: hc(:), s2(:)
+
+  n_states_diag = 1
+  TOUCH n_states_diag
+  call generate_fci_space
+  allocate(c_pert(N_det,0:mp_order))
+  allocate(s2(N_det))
+  allocate(e_pert(mp_order+1))
+  e_pert = 0.d0
+  c_pert(:,:) = 0.d0
+  c_pert(1,0) = 1.d0
+
+  e_pert(1) = nuclear_repulsion
+  do k=1,mp_order
+    ! H_ij C^(k-1)
+    call h_s2_u_0_nstates_zmq(c_pert(1,k),s2,c_pert(1,k-1),1,N_det)
+    e_pert(k) += c_pert(1,k)
+    print *, k, e_pert(k), sum(e_pert) 
+
+    c_pert(1,k) = 0.d0
+    c_pert(:,k) = -c_pert(:,k)
+    do l=1,k
+      do i=2,N_det
+        c_pert(i,k) = c_pert(i,k) + e_pert(l) * c_pert(i,k-l)
+      enddo
+    enddo
+    do i=2,N_det
+      c_pert(i,k) = c_pert(i,k) + energy_det_i(i) * c_pert(i,k-1)
+    enddo
+    do i=2,N_det
+      c_pert(i,k) = c_pert(i,k) / (energy_det_i(i) - energy_det_i(1))
+    enddo
+
+  enddo
+
+
+end
+
+subroutine generate_fci_space
+  use bitmasks
+  implicit none
+  integer :: i, sze
+  integer(bit_kind) :: o(N_int,2)
+
+  if (mo_num > 64) then
+     stop 'No more than 64 MOs'
+  endif
+  o(:,1) = full_ijkl_bitmask(:)
+  o(:,2) = 0_bit_kind
+
+  call configuration_to_dets_size(o,n_det_alpha_unique,elec_alpha_num,N_int)
+  TOUCH n_det_alpha_unique
+
+  integer :: k,n,m, t, t1, t2
+  integer(bit_kind) :: u
+  k=0
+  n = elec_alpha_num
+  m = mo_num - n
+  u = shiftl(1_bit_kind,n) -1
+  do while (u < shiftl(1_bit_kind,n+m))
+    k = k+1
+    psi_det_alpha_unique(1,k) = u
+    t = ior(u,u-1)
+    t1 = t+1
+    t2 = shiftr((iand(not(t),t1)-1), trailz(u)+1)
+    u = ior(t1,t2)
+  enddo
+
+
+  call configuration_to_dets_size(o,n_det_beta_unique,elec_beta_num,N_int)
+  TOUCH n_det_beta_unique
+
+  k=0
+  n = elec_beta_num
+  m = mo_num - n
+  u = shiftl(1_bit_kind,n) -1
+  do while (u < shiftl(1_bit_kind,n+m))
+    k = k+1
+    psi_det_beta_unique(1,k) = u
+    t = ior(u,u-1)
+    t1 = t+1
+    t2 = shiftr((iand(not(t),t1)-1), trailz(u)+1)
+    u = ior(t1,t2)
+  enddo
+
+  call generate_all_alpha_beta_det_products
+
+  print *, N_det
+
+
+end
+

devel/svdwf/svdwf.irp.f

@@ -12,36 +12,53 @@ subroutine run
   implicit none
   double precision, allocatable :: U(:,:), Vt(:,:), D(:), A(:,:)
   integer :: i, j, k, p, q
+  double precision :: entropy
 
   allocate( A (n_det_alpha_unique, n_det_beta_unique), &
             U (n_det_alpha_unique, n_det_alpha_unique), &
             Vt(n_det_beta_unique, n_det_beta_unique),  &
             D(max(n_det_beta_unique,n_det_alpha_unique)) )
 
-  A = 0.D0
+  do j=1,n_det_beta_unique
+    do i=1,n_det_alpha_unique
+      A(i,j)  = 0.D0
+    enddo
+  enddo
+
   do k=1,N_det
     i = psi_bilinear_matrix_rows(k)
     j = psi_bilinear_matrix_columns(k)
     A(i,j) = psi_bilinear_matrix_values(k,1)
   enddo
 
-  call randomized_svd(A, size(A,1), &
+  if (N_det == 1) then
+    D(1) = 1.d0
+    U(1,1) = 1.d0
+    Vt(1,1) = 1.d0
+  else
+    call randomized_svd(A, size(A,1), &
            U, size(U,1), D, Vt, size(Vt,1), n_det_alpha_unique, n_det_beta_unique, &
-           6,1000)
+           6,min(n_det_beta_unique,1000))
+  endif
 
+  entropy = 0.d0
   do i=1,n_det_beta_unique
     print *, i, real(D(i)), real(D(i)**2), real(sum(D(1:i)**2))
+    entropy -= D(i) * dlog(D(i))
     if (D(i) < 1.d-15) then
       k = i
       exit
     endif
   enddo
   print *, 'threshold: ', 2.858 * D(k/2)
-  do i=1,n_det_alpha_unique
-    print '(I6,4(X,F12.8))', i, U(i,1:4)
-  enddo
-  print *, ''
-  do i=1,n_det_beta_unique
-    print '(I6,4(X,F12.8))', i, Vt(1:4,i)
-  enddo
+  print *, 'Entropy  : ', entropy
+
+!  do i=1,n_det_alpha_unique
+!    print '(I6,4(X,F12.8))', i, U(i,1:4)
+!  enddo
+!  print *, ''
+!  do i=1,n_det_beta_unique
+!    print '(I6,4(X,F12.8))', i, Vt(1:4,i)
+!  enddo
+
 end

stable/amplitudes/EZFIO.cfg.save

@@ -1,7 +0,0 @@
-[t1_amplitudes]
-type: double precision
-doc: Amplitudes for the single-excitation operator
-interface: ezfio,provider
-size: (mo_basis.mo_num,mo_basis.mo_num)
-
-