Work on tc

devel/dmc_dress/EZFIO.cfg

@@ -4,3 +4,9 @@ doc: Dressing matrix obtained from DMC
 size: (determinants.n_det)
 interface: ezfio, provider
 
+[dmc_delta_htc]
+type: double precision
+doc: Dressing matrix obtained from H_TC
+size: (determinants.n_det)
+interface: ezfio
+

devel/dmc_dress/dress.py

@@ -12,21 +12,26 @@ sys.path.insert(0,QP_PATH)
 from ezfio import ezfio
 
 def read_hamiltonian(inp):
-    text = subprocess.run(["qmcchem", "result", inp], capture_output=True).stdout
+#    text = subprocess.run(["qmcchem", "result", inp], capture_output=True).stdout
+    text = subprocess.run(["cat", "OUT"], capture_output=True).stdout
     inside = None
     h = []
     s = []
     norm = None
     for line in text.splitlines():
-       line = str(line)
-       print (line)
+       line = line.decode("utf-8")
+       if line == "":
+         continue
        if "Psi_norm :" in line:
-           norm = float(line.split()[3])
+           norm = float(line.split()[2])
        if "]" in line:
            inside = None
        if inside == "H":
           data = line.split()
-          h.append(float(data[3]))
+#          h.append(float(data[2]))                            #ave
+          h.append(max(float(data[2]) - float(data[4]),0.))   #diff
+#          if float(data[2]) - float(data[4])>0.: h.append(float(data[2]))  
+#          else: h.append(0.)
        elif "Ci_dress" in line:
            inside = "H"
     h = np.array(h)/norm
@@ -38,3 +43,4 @@ h = read_hamiltonian(sys.argv[1])
 ezfio.set_file(sys.argv[1])
 ezfio.set_dmc_dress_dmc_delta_h(h)
 print(h)
+

devel/dmc_dress/dressing_vector.irp.f

@@ -11,7 +11,7 @@
  double precision, allocatable :: delta(:)
 
  allocate(delta(N_det))
- delta(1:N_det) = dmc_delta_h(1:N_det)
+ delta (1:N_det) = dmc_delta_h(1:N_det) + dmc_delta_htc(1:N_det) 
 
  call dset_order(delta,psi_bilinear_matrix_order_reverse,N_det)
 
@@ -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

devel/trexio/import_trexio_integrals.irp.f

@@ -24,7 +24,8 @@ subroutine run
   double precision, allocatable :: A(:,:)
   double precision, allocatable :: V(:)
   integer         , allocatable :: Vi(:,:)
-  double precision, allocatable :: s
+  double precision              :: s
+  PROVIDE ao_num
 
   f = trexio_open(trexio_filename, 'r', TREXIO_AUTO, rc)
   if (f == 0_8) then
@@ -37,8 +38,9 @@ subroutine run
 
   if (trexio_has_nucleus_repulsion(f) == TREXIO_SUCCESS) then
     rc = trexio_read_nucleus_repulsion(f, s)
+    call trexio_assert(rc, TREXIO_SUCCESS)
     if (rc /= TREXIO_SUCCESS) then
-      print *, irp_here
+      print *, irp_here, rc
       print *, 'Error reading nuclear repulsion'
       stop -1
     endif
@@ -100,20 +102,21 @@ subroutine run
 
   ! AO 2e integrals
   ! ---------------
+  PROVIDE ao_integrals_map
 
-  allocate(buffer_i(ao_num**3), buffer_values(ao_num**3))
-  allocate(Vi(4,ao_num**3), V(ao_num**3))
+  integer*4 :: BUFSIZE
+  BUFSIZE=ao_num**2
+  allocate(buffer_i(BUFSIZE), buffer_values(BUFSIZE))
+  allocate(Vi(4,BUFSIZE), V(BUFSIZE))
 
   integer*8 :: offset, icount
 
+  open(unit=104,file='tmp')
   offset = 0_8
-  icount = 0_8
+  icount = BUFSIZE
   rc = TREXIO_SUCCESS
   do while (icount == size(V))
     rc = trexio_read_ao_2e_int_eri(f, offset, icount, Vi, V)
-    if (rc /= TREXIO_SUCCESS) then
-        exit
-    endif
     do m=1,icount
       i = Vi(1,m)
       j = Vi(2,m)
@@ -122,11 +125,16 @@ subroutine run
       integral = V(m)
       call two_e_integrals_index(i, j, k, l, buffer_i(m) )
       buffer_values(m) = integral
+  write(104,'(4(I5,X),2D22.15)') i,j,k,l, integral
     enddo
     call insert_into_ao_integrals_map(int(icount,4),buffer_i,buffer_values)
     offset = offset + icount
+    if (rc /= TREXIO_SUCCESS) then
+        exit
+    endif
   end do
   n_integrals = offset
+  close(104)
 
   call map_sort(ao_integrals_map)
   call map_unique(ao_integrals_map)

devel/trexio/install

@@ -10,4 +10,14 @@ fi
 
 pkg-config --libs trexio > LIB
 
+scripts_list="qp_import_trexio.py"
+
+# Destroy ONLY the symbolic link for the scripts to be used in the
+# ${QP_ROOT}/scripts/ directory.
+for i in $scripts_list
+do
+  find ${QP_ROOT}/scripts/$i -type l -delete
+done
+
+# Create symlink in scripts
 ln --symbolic ${PWD}/qp_import_trexio.py $QP_ROOT/scripts

devel/trexio/qp_import_trexio.py

@@ -68,21 +68,13 @@ def write_ezfio(trexio_filename, filename):
         trexio_file = trexio.File(trexio_filename,mode='r',back_end=trexio.TREXIO_HDF5)
 
     basis_type = trexio.read_basis_type(trexio_file)
-    if basis_type.lower() != "gaussian":
-        raise TypeError
+#    if basis_type.lower() != "gaussian":
+#        raise TypeError
 
     ezfio.set_file(filename)
+    ezfio.set_trexio_trexio_file(trexio_filename)
 
 
-    print("Electrons\t...\t", end=' ')
-
-    num_alpha = trexio.read_electron_up_num(trexio_file)
-    num_beta  = trexio.read_electron_dn_num(trexio_file)
-    ezfio.set_electrons_elec_alpha_num(num_alpha)
-    ezfio.set_electrons_elec_beta_num(num_beta)
-
-    print("OK")
-
 
     print("Nuclei\t\t...\t", end=' ')
 
@@ -106,6 +98,23 @@ def write_ezfio(trexio_filename, filename):
     print("OK")
 
 
+    print("Electrons\t...\t", end=' ')
+
+    try:
+        num_beta = trexio.read_electron_dn_num(trexio_file)
+    except:
+        num_beta = sum(charge)//2
+
+    try:
+        num_alpha = trexio.read_electron_up_num(trexio_file)
+    except:
+        num_alpha = sum(charge) - num_beta
+
+    ezfio.set_electrons_elec_alpha_num(num_alpha)
+    ezfio.set_electrons_elec_beta_num(num_beta)
+
+    print("OK")
+
     print("Basis\t\t...\t", end=' ')
 
     ezfio.set_basis_basis("Read from TREXIO")
@@ -229,6 +238,7 @@ def write_ezfio(trexio_filename, filename):
             coef.append(coefficient[j])
             expo.append(exponent[j])
 
+    ezfio.set_ao_basis_ao_prim_num_max(prim_num_max)
     ezfio.set_ao_basis_ao_coef(coef)
     ezfio.set_ao_basis_ao_expo(expo)
     ezfio.set_ao_basis_ao_basis("Read from TREXIO")
@@ -265,17 +275,21 @@ def write_ezfio(trexio_filename, filename):
     except trexio.Error:
       pass
 
-    MoMatrix = trexio.read_mo_coefficient(trexio_file)
-    mo_num = trexio.read_mo_num(trexio_file)
+    try:
+      mo_num = trexio.read_mo_num(trexio_file)
+      ezfio.set_mo_basis_mo_num(mo_num)
 
-    ezfio.set_mo_basis_mo_num(mo_num)
-    ezfio.set_mo_basis_mo_coef(MoMatrix)
-    mo_occ = [ 0. for i in range(mo_num) ]
-    for i in range(num_alpha):
-       mo_occ[i] += 1.
-    for i in range(num_beta):
-       mo_occ[i] += 1.
-    ezfio.set_mo_basis_mo_occ(mo_occ)
+      MoMatrix = trexio.read_mo_coefficient(trexio_file)
+      ezfio.set_mo_basis_mo_coef(MoMatrix)
+
+      mo_occ = [ 0. for i in range(mo_num) ]
+      for i in range(num_alpha):
+         mo_occ[i] += 1.
+      for i in range(num_beta):
+         mo_occ[i] += 1.
+      ezfio.set_mo_basis_mo_occ(mo_occ)
+    except:
+      pass
 
     print("OK")
 

devel/trexio/uninstall

@@ -8,7 +8,7 @@ then
   exit -1
 fi
 
-scripts_list="qp_import_trexio.py "
+scripts_list="qp_import_trexio.py"
 
 # Destroy ONLY the symbolic link for the scripts to be used in the
 # ${QP_ROOT}/scripts/ directory.

fnmf/EZFIO.cfg

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

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 - nuclear_repulsion
+
+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,76 @@
+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(:,:)
+ double precision, allocatable :: htmp(:,:), stmp(:,:)
+
+ allocate(Htmp(N_det,N_det), Stmp(N_det,N_det))
+
+ htmp(:, :) = H_dmc_mat(1:N_det, 1:N_det)
+ stmp(:, :) = S_dmc_mat(1:N_det, 1:N_det)
+
+ htmp(1, 1) = H_dmc_mat(N_det+1, N_det+1)
+ stmp(1, 1) = S_dmc_mat(N_det+1, N_det+1)
+
+ htmp(2:N_det, 1) = H_dmc_mat(2:N_det, N_det+1)
+ stmp(2:N_det, 1) = S_dmc_mat(2:N_det, N_det+1)
+
+ htmp(1, 2:N_det) = H_dmc_mat(N_det+1, 2:N_det)
+ stmp(1, 2:N_det) = S_dmc_mat(N_det+1, 2:N_det)
+
+
+ allocate(beta(N_det), vr(N_det,N_det), vl(N_det, N_det))
+
+ call lapack_g_non_sym_real(N_det, Htmp, size(Htmp,1), &
+  Stmp, size(Stmp,1), beta, n_real, vl, size(vl,1), vr, size(vr,1))
+
+
+  vl(:,1) = psi_coef(:,1)
+
+  psi_coef(:,1) *= vr(1,1)
+  do i=2,N_det
+    psi_coef(i,1) += vr(i,1)
+  end do
+  double precision, external :: dnrm2, ddot
+  double precision :: norm
+
+  norm = dnrm2(N_det, psi_coef, 1)
+  psi_coef(:,1) *= 1.d0/norm
+  SOFT_TOUCH psi_coef
+
+  print *, '-------------------------------------------------'
+  print *, ' psi_old  psi_new'
+  print *, '-------------------------------------------------'
+  do i=1,N_det
+    print '(2(E16.8,X))', vl(i,1), psi_coef(i,1) 
+  enddo
+  print *, '-------------------------------------------------'
+
+  print *,  ''
+  print *, 'Overlap        : ', ddot(N_det, psi_coef, 1, vl, 1)
+  print *, 'DMC energy     : ', htmp(1,1) + nuclear_repulsion
+  print *, 'Updated energy : ', beta(1) + nuclear_repulsion
+
+  print *,  'H'
+  do i=1,N_det
+    print *,  htmp(i, 1:N_det)
+  enddo
+  print *,  'S'
+  do i=1,N_det
+    print *,  stmp(i, 1:N_det)
+  enddo
+  call save_wavefunction
+ 
+end
+

fnmf/non_hermit.irp.f

@@ -0,0 +1,74 @@
+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(:), Htmp(:,:), Stmp(:,:)
+  double precision, allocatable :: alphar(:), alphai(:), vltmp(:,:), vrtmp(:,:)
+  integer, allocatable :: iorder(:), list_good(:)
+
+
+  lwork = -1
+  allocate(work(1), alphar(n), alphai(n), vltmp(n,n), vrtmp(n,n), &
+  Htmp(n,n), Stmp(n,n), list_good(n))
+
+  Htmp(1:n,1:n) = H(1:n,1:n)
+  Stmp(1:n,1:n) = S(1:n,1:n)
+  call dggev('V', 'V', n, Htmp, size(Htmp,1), Stmp, size(Stmp,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, Htmp, size(Htmp,1), Stmp, size(Stmp,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
+    if (dabs(alphai(i)) < 1.d-10) then 
+       n_real += 1
+       list_good(n_real) = i
+    else
+       alphar(i) = dble(huge(1.0))
+    endif
+  enddo
+
+  do i=1,n
+    if (dabs(beta(i)/(alphar(i)+1.d-12)) < 1.d-10) then
+      beta(i) = 0.d0
+    else
+      beta(i) = alphar(i)/beta(i)
+    endif
+  enddo
+
+  do i=1,n_real
+    iorder(i) = i
+    do j=1,n
+      vrtmp(j,i) = vrtmp(j,list_good(i))
+      vltmp(j,i) = vltmp(j,list_good(i))
+    end do
+  end do
+
+  call dsort(beta, iorder, n_real)
+
+  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, Htmp, Stmp)
+
+end
+