Initial tests

devel/dmc_dress/EZFIO.cfg

@@ -8,5 +8,5 @@ interface: ezfio, provider
 type: double precision
 doc: Dressing matrix obtained from H_TC
 size: (determinants.n_det)
-interface: ezfio, provider
+interface: ezfio
 

devel/dmc_dress/dressing_vector.irp.f

@@ -31,3 +31,38 @@ END_PROVIDER
 
 
 
+BEGIN_PROVIDER [ double precision, dmc_delta_htc , (n_det) ]
+  implicit none
+  BEGIN_DOC
+! Dressing matrix obtained from H_TC
+  END_DOC
+
+  logical                        :: has
+  PROVIDE ezfio_filename
+  if (mpi_master) then
+    if (size(dmc_delta_htc) == 0) return
+
+    call ezfio_has_dmc_dress_dmc_delta_htc(has)
+    if (has) then
+      write(6,'(A)') '.. >>>>> [ IO READ: dmc_delta_htc ] <<<<< ..'
+      call ezfio_get_dmc_dress_dmc_delta_htc(dmc_delta_htc)
+    else
+      dmc_delta_htc(:) = 0.d0
+    endif
+  endif
+  IRP_IF MPI_DEBUG
+    print *,  irp_here, mpi_rank
+    call MPI_BARRIER(MPI_COMM_WORLD, ierr)
+  IRP_ENDIF
+  IRP_IF MPI
+    include 'mpif.h'
+    integer :: ierr
+    call MPI_BCAST( dmc_delta_htc, (n_det), MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
+    if (ierr /= MPI_SUCCESS) then
+      stop 'Unable to read dmc_delta_htc with MPI'
+    endif
+  IRP_ENDIF
+
+  call write_time(6)
+
+END_PROVIDER

fnmf/EZFIO.cfg

@@ -0,0 +1,17 @@
+[E_dmc]
+type: double precision
+doc: DMC energy
+interface: ezfio, provider
+
+[dmc_h]
+type: double precision
+doc: Dressing matrix obtained from DMC
+size: (determinants.n_det)
+interface: ezfio, provider
+
+[dmc_s]
+type: double precision
+doc: Dressing matrix obtained from H_TC
+size: (determinants.n_det)
+interface: ezfio
+

fnmf/NEED

@@ -0,0 +1,2 @@
+determinants
+davidson_undressed

fnmf/README.rst

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

fnmf/dmc_data.irp.f

@@ -0,0 +1,70 @@
+ BEGIN_PROVIDER [ double precision, h_dmc_row , (N_det) ]
+&BEGIN_PROVIDER [ double precision, s_dmc_row , (N_det) ]
+ implicit none
+ BEGIN_DOC
+ ! Data sampled with QMC=Chem
+ END_DOC
+
+! h_dmc_row(:) = h_dmc(:)
+! s_dmc_row(:) = s_dmc(:)
+! call dset_order(h_dmc_row,psi_bilinear_matrix_order_reverse,N_det)
+! call dset_order(s_dmc_row,psi_bilinear_matrix_order_reverse,N_det)
+
+  integer :: i
+  do i=1,N_det
+    s_dmc_row(i) = psi_coef(i,1)
+    call i_h_psi(psi_det(1,1,i), psi_det, psi_coef, N_int, N_det, &
+         N_det, 1, h_dmc_row(i) )
+  enddo
+
+END_PROVIDER
+
+BEGIN_PROVIDER [ integer, mat_size ]
+ implicit none
+ BEGIN_DOC
+! Size of the matrices
+ END_DOC
+ mat_size = N_det+1
+END_PROVIDER
+
+BEGIN_PROVIDER [ double precision, H_dmc_mat, (mat_size, mat_size) ]
+ implicit none
+ BEGIN_DOC
+! Hamiltonian extended with DMC data
+ END_DOC
+ integer :: i,j
+ do j=1,N_det
+   do i=1,N_det
+      call i_h_j(psi_det(1,1,i), psi_det(1,1,j), N_int, H_dmc_mat(i,j))
+   enddo
+ enddo
+
+ do i=1,N_det
+   call i_h_psi(psi_det(1,1,i), psi_det, psi_coef, N_int, N_det, &
+        N_det, 1, H_dmc_mat(i,N_det+1) )
+   H_dmc_mat(N_det+1,i) = h_dmc_row(i)
+ enddo
+ H_dmc_mat(mat_size,mat_size) = E_dmc
+
+END_PROVIDER
+
+BEGIN_PROVIDER [ double precision, S_dmc_mat, (mat_size, mat_size) ]
+ implicit none
+ BEGIN_DOC
+! Overlap matrix extended with DMC data
+ END_DOC
+ integer :: i,j
+ S_dmc_mat = 0.d0
+ do i=1,mat_size
+   S_dmc_mat(i,i) = 1.d0
+ enddo
+
+ do i=1,N_det
+   S_dmc_mat(i,N_det+1) = psi_coef(i,1)
+   S_dmc_mat(N_det+1,i) = S_dmc_row(i)
+ enddo
+
+END_PROVIDER
+
+
+

fnmf/fnmf.irp.f

@@ -0,0 +1,27 @@
+program fnmf
+  implicit none
+  BEGIN_DOC
+! TODO : Put the documentation of the program here
+  END_DOC
+  read_wf = .True.
+  TOUCH read_wf
+  call run
+end
+
+
+subroutine run
+ implicit none
+ integer :: i, n_real
+ double precision, allocatable :: beta(:), vr(:,:), vl(:,:)
+
+ allocate(beta(mat_size), vr(mat_size,mat_size), vl(mat_size, mat_size))
+ call lapack_g_non_sym_real(mat_size, H_dmc_mat, size(H_dmc_mat,1), &
+  S_dmc_mat, size(S_dmc_mat,1), beta, n_real, vl, size(vl,1), vr, size(vr,1))
+
+  print *, 'EV     VR    VL'
+  print *, '---------------------------'
+  do i=1,mat_size
+    print '(3(F16.12,X))', beta(i), vr(i,1), vl(i,1)
+  enddo
+end
+

fnmf/non_hermit.irp.f

@@ -0,0 +1,55 @@
+subroutine lapack_g_non_sym_real(n, H, LDH, S, LDS, beta, &
+      n_real, vl, LDVL, vr, LDVR)
+  implicit none
+  integer, intent(in) :: n, LDH, LDS, LDVL, LDVR
+  double precision, intent(in) :: H(LDH,n), S(LDS,n)
+  double precision, intent(out) :: VL(LDVL,n), VR(LDVR,n), beta(n)
+  integer, intent(out) :: n_real
+
+  integer :: lwork, info, i,j 
+  double precision, allocatable :: work(:)
+  double precision, allocatable :: alphar(:), alphai(:), vltmp(:,:), vrtmp(:,:)
+  integer, allocatable :: iorder(:)
+
+
+  lwork = -1
+  allocate(work(1), alphar(n), alphai(n), vltmp(n,n), vrtmp(n,n))
+
+  call dggev('V', 'V', n, H, size(H,1), S, size(S,1), alphar, alphai, beta, &
+             vltmp, size(vltmp,1), vrtmp, size(vrtmp,1), work, lwork, info)
+
+  lwork = int(work(1))
+  deallocate(work)
+  allocate(work(lwork))
+  call dggev('V', 'V', n, H, size(H,1), S, size(S,1), alphar, alphai, beta, &
+             vltmp, size(vltmp,1), vrtmp, size(vrtmp,1), work, lwork, info)
+
+  deallocate(work)
+  if (info /= 0) then
+    stop 'DGGEV Diagonalization failed'
+  endif
+
+  allocate(iorder(n))
+  n_real = 0
+  do i=1,n
+    iorder(i) = i
+    if (dabs(alphai(i)) < 1.d-10) then 
+       n_real += 1
+       alphar(i) = dble(huge(1.0))
+    endif
+  enddo
+
+  beta(:) = alphar(:)/beta(:)
+  call dsort(beta, iorder, n)
+
+  do i=1,n_real
+    do j=1,n
+      vr(j,i) = vrtmp(j,iorder(i))
+      vl(j,i) = vltmp(j,iorder(i))
+    end do
+  end do
+
+  deallocate(vrtmp, vltmp, iorder)
+
+end
+