dressing multi states

devel/dmc_dress/EZFIO.cfg

@@ -1,7 +1,7 @@
 [dmc_delta_h]
 type: double precision
 doc: Dressing matrix obtained from DMC
-size: (determinants.n_det)
+size: (determinants.n_det,determinants.n_states)
 interface: ezfio, provider
 
 [la]

devel/dmc_dress/dmc_dress_1state.irp.f

@@ -0,0 +1,36 @@
+program dmc_dress_1state
+ implicit none
+ BEGIN_DOC
+! Program that extracts the lowest states of the Hamiltonian dressed by the QMC
+! dressing vector stored in :option:`dmc_dressing dmc_delta_h`
+!
+ END_DOC
+ read_wf = .True.
+ touch read_wf
+! call pre
+ call routine
+ call save_wavefunction_general(N_det,N_states,psi_det_sorted,size(psi_coef_sorted,1),psi_coef_sorted)
+end
+
+subroutine pre
+ implicit none
+ double precision, allocatable :: left(:,:), right(:,:), tmp(:,:), res(:,:)
+ integer :: i
+ allocate (left(1,1:N_det), right(1:N_det,1), tmp(1:N_det,1), res(1,1))
+ left(1,1:N_det) = psi_coef(1:N_det,1)
+ right(1:N_det,1) = psi_coef(1:N_det,1)
+
+ tmp(1:N_det,1:1) = matmul(h_matrix_dressed(1:N_det,1:N_det), right(1:N_det,1:1))
+ res(1:1,1:1) = matmul(left(1:1,1:N_det), tmp(1:N_det,1:1))
+ print *,  'E_in = ', res(1,1)
+ do i=1,N_det
+ print *,  'HPsi/c0 = ', tmp(i,1)/psi_coef(i,1)  
+ enddo
+end
+subroutine routine
+ implicit none
+ psi_coef(1:N_det,1) = ci_eigenvectors_dressed(1:N_det,1)
+ print*,'N_det = ',N_det
+ print *, 'E = ', ci_energy_dressed(1) + nuclear_repulsion
+ SOFT_TOUCH psi_coef
+end

devel/dmc_dress/dmc_dress_multistate.irp.f

@@ -0,0 +1,37 @@
+program dmc_dress_multistate
+
+  BEGIN_DOC
+  ! Program that extracts the lowest states of the Hamiltonian dressed by the QMC
+  ! dressing vector stored in :option:`dmc_dressing dmc_delta_h`
+  !
+ END_DOC
+
+  implicit none
+
+  read_wf = .True.
+  touch read_wf
+ 
+  call routine()
+  call save_wavefunction_general(N_det, N_states, psi_det_sorted, size(psi_coef_sorted, 1), psi_coef_sorted)
+
+end
+
+! ---
+
+subroutine routine
+
+  implicit none
+  integer :: k
+
+  do k = 1, N_states
+    print *, ' state:', k
+    psi_coef(1:N_det,k) = ci_eigenvectors_nonsym_dressed(1:N_det,k)
+    print *, 'E = ', ci_energy_nonsym_dressed(k) + nuclear_repulsion
+  enddo
+
+  SOFT_TOUCH psi_coef
+
+end
+
+! ---
+

devel/dmc_dress/dressing_vector.irp.f

@@ -1,33 +1,76 @@
+
+! ---
+
  BEGIN_PROVIDER [ double precision, dressing_column_h, (N_det,N_states) ]
 &BEGIN_PROVIDER [ double precision, dressing_column_s, (N_det,N_states) ]
- implicit none
- BEGIN_DOC
- ! \Delta_{state-specific}. \Psi
- ! Diagonal element is divided by 2 because Delta = D + D^t
- END_DOC
 
- integer :: i,ii,k,j, l
- double precision :: f, tmp
- double precision, allocatable :: delta(:)
+  BEGIN_DOC
+  ! \Delta_{state-specific}. \Psi
+  ! Diagonal element is divided by 2 because Delta = D + D^t
+  END_DOC
 
- allocate(delta(N_det))
- delta(1:N_det) = dmc_delta_h(1:N_det)
+  implicit none
+  integer :: i,ii,k,j, l
+  double precision :: f, tmp
+  double precision, allocatable :: delta(:)
 
- call dset_order(delta,psi_bilinear_matrix_order_reverse,N_det)
+  allocate(delta(N_det))
+  delta(1:N_det) = dmc_delta_h(1:N_det,1)
 
- dressing_column_h(:,:) = 0.d0
- dressing_column_s(:,:) = 0.d0
+  call dset_order(delta,psi_bilinear_matrix_order_reverse,N_det)
 
- l = dressed_column_idx(1)
- do j = 1, n_det
-   if (j == l) cycle
-   dressing_column_h(j,1)  = delta(j) 
-   dressing_column_h(l,1) -= psi_coef(j,1) * delta(j) / psi_coef(l,1)
- enddo
- dressing_column_h(l,1) += delta(l) 
- dressing_column_h(l,1) *= 0.5d0
+  dressing_column_h(:,:) = 0.d0
+  dressing_column_s(:,:) = 0.d0
+
+  l = dressed_column_idx(1)
+  do j = 1, n_det
+    if(j == l) cycle
+    dressing_column_h(j,1)  = delta(j) 
+    dressing_column_h(l,1) -= psi_coef(j,1) * delta(j) / psi_coef(l,1)
+  enddo
+  dressing_column_h(l,1) += delta(l) 
+  dressing_column_h(l,1) *= 0.5d0
+
+  deallocate(delta)
 
 END_PROVIDER
 
+! ---
 
+BEGIN_PROVIDER [ double precision, dressing_delta, (N_det, N_states) ]
+
+  BEGIN_DOC
+  ! 
+  ! dressing_delta is:
+  ! [\delta_K]_I = < I | \tilde{H} - H | \Phi_K >
+  !
+  END_DOC
+
+  implicit none
+  integer                       :: i, j, k
+  double precision, allocatable :: delta(:,:)
+
+  dressing_delta(1:N_det,1:N_states) = 0.d0
+  
+  allocate(delta(N_det,N_states))
+
+  do k = 1, N_states
+
+    do j = 1, N_det
+      delta(j,k) = dmc_delta_h(j,k)
+    enddo
+
+    call dset_order(delta(1:N_det,k), psi_bilinear_matrix_order_reverse, N_det)
+    
+    do j = 1, N_det
+      dressing_delta(j,k) = delta(j,k) 
+    enddo
+
+  enddo
+
+  deallocate(delta)
+
+END_PROVIDER
+
+! ---