Merge branch 'master' of gitlab.com:scemama/qp_plugins_scemama

devel/trexio/EZFIO.cfg

@@ -10,3 +10,21 @@ doc: Name of the exported TREXIO file
 interface: ezfio, ocaml, provider
 default: None
 
+[export_rdm]
+type: logical
+doc: If True, export two-body reduced density matrix
+interface: ezfio, ocaml, provider
+default: False
+
+[export_ao_ints]
+type: logical
+doc: If True, export two-electron integrals in AO basis
+interface: ezfio, ocaml, provider
+default: False
+
+[export_mo_ints]
+type: logical
+doc: If True, export two-electron integrals in MO basis
+interface: ezfio, ocaml, provider
+default: True
+

devel/trexio/LIB

@@ -1,2 +1,2 @@
--L/home/scemama/TREX/trexio/_install/lib -ltrexio
+-ltrexio
 

devel/trexio/NEED

@@ -4,3 +4,4 @@ mo_one_e_ints
 mo_two_e_ints
 ao_two_e_ints
 ao_one_e_ints
+two_body_rdm

devel/trexio/export_trexio.irp.f

@@ -1,4 +1,11 @@
 program export_trexio
+  implicit none
+  read_wf = .True.
+  SOFT_TOUCH read_wf
+  call run
+end
+
+subroutine run
   use trexio
   implicit none
   BEGIN_DOC
@@ -7,6 +14,7 @@ program export_trexio
 
   integer(8)                     :: f ! TREXIO file handle
   integer                        :: rc
+  double precision, allocatable  :: factor(:)
 
   print *, 'TREXIO file : '//trim(trexio_filename)
   print *, ''
@@ -30,6 +38,8 @@ program export_trexio
 ! Electrons
 ! ---------
 
+  print *, 'Electrons'
+
   rc = trexio_write_electron_up_num(f, elec_alpha_num)
   call check_success(rc)
 
@@ -40,6 +50,8 @@ program export_trexio
 ! Nuclei
 ! ------
 
+  print *, 'Nuclei'
+
   rc = trexio_write_nucleus_num(f, nucl_num)
   call check_success(rc)
 
@@ -52,10 +64,15 @@ program export_trexio
   rc = trexio_write_nucleus_label(f, nucl_label, 32)
   call check_success(rc)
 
+  rc = trexio_write_nucleus_repulsion(f, nuclear_repulsion)
+  call check_success(rc)
+
 
 ! Pseudo-potentials
 ! -----------------
 
+  print *, 'ECP'
+
   double precision, allocatable :: tmp_double(:,:)
   integer, allocatable          :: tmp_int(:,:)
 
@@ -94,39 +111,36 @@ program export_trexio
 ! Basis
 ! -----
 
+  print *, 'Basis'
+
+
   rc = trexio_write_basis_type(f, 'Gaussian', len('Gaussian'))
   call check_success(rc)
 
-  rc = trexio_write_basis_num(f, shell_num)
-  call check_success(rc)
-
-  rc = trexio_write_basis_nucleus_shell_num(f, nucleus_shell_num)
-  call check_success(rc)
-
-  rc = trexio_write_basis_nucleus_index(f, basis_nucleus_index)
-  call check_success(rc)
-
-  rc = trexio_write_basis_shell_ang_mom(f, shell_ang_mom)
-  call check_success(rc)
-
   rc = trexio_write_basis_prim_num(f, prim_num)
   call check_success(rc)
 
-  rc = trexio_write_basis_shell_prim_num(f, shell_prim_num)
-  call check_success(rc)
+   rc = trexio_write_basis_shell_num(f, shell_num)
+   call check_success(rc)
 
-  double precision, allocatable :: factor(:)
-  allocate(factor(shell_num))
-  if (ao_normalized) then
-    factor(1:shell_num) = shell_normalization_factor(1:shell_num)
-  else
-    factor(1:shell_num) = 1.d0
-  endif
-  rc = trexio_write_basis_shell_factor(f, factor)
-  call check_success(rc)
-  deallocate(factor)
+   rc = trexio_write_basis_nucleus_index(f, basis_nucleus_index)
+   call check_success(rc)
 
-  rc = trexio_write_basis_shell_prim_index(f, shell_prim_index)
+   rc = trexio_write_basis_shell_ang_mom(f, shell_ang_mom)
+   call check_success(rc)
+
+   allocate(factor(shell_num))
+   if (ao_normalized) then
+     factor(1:shell_num) = shell_normalization_factor(1:shell_num)
+   else
+     factor(1:shell_num) = 1.d0
+   endif
+   rc = trexio_write_basis_shell_factor(f, factor)
+   call check_success(rc)
+
+   deallocate(factor)
+
+  rc = trexio_write_basis_shell_index(f, shell_index)
   call check_success(rc)
 
   rc = trexio_write_basis_exponent(f, prim_expo)
@@ -146,10 +160,11 @@ program export_trexio
   deallocate(factor)
 
 
-
 ! Atomic orbitals
 ! ---------------
 
+  print *, 'AOs'
+
   rc = trexio_write_ao_num(f, ao_num)
   call check_success(rc)
 
@@ -168,7 +183,6 @@ program export_trexio
   C_A(3) = 0.d0
 
   allocate(factor(ao_num))
-  print *, ao_first_of_shell
   if (ao_normalized) then
     do i=1,ao_num
       l = ao_first_of_shell(ao_shell(i))
@@ -184,6 +198,8 @@ program export_trexio
 ! One-e AO integrals
 ! ------------------
 
+  print *, 'AO integrals'
+
   rc = trexio_write_ao_1e_int_overlap(f,ao_overlap)
   call check_success(rc)
 
@@ -204,10 +220,59 @@ program export_trexio
   rc = trexio_write_ao_1e_int_core_hamiltonian(f,ao_one_e_integrals)
   call check_success(rc)
 
+! Two-e AO integrals
+! ------------------
+
+  if (export_ao_ints) then
+    PROVIDE ao_two_e_integrals_in_map
+
+    integer(8), parameter :: BUFSIZE=10000_8
+    double precision :: eri_buffer(BUFSIZE), integral
+    integer(4) :: eri_index(4,BUFSIZE)
+    integer(8) :: icount, offset
+
+    double precision, external :: get_ao_two_e_integral
+
+
+    icount = 0_8
+    offset = 0_8
+    do l=1,ao_num
+      do k=1,ao_num
+        do j=l,ao_num
+          do i=k,ao_num
+            if (i==j .and. k<l) cycle
+            if (i<j) cycle
+            integral = get_ao_two_e_integral(i,j,k,l,ao_integrals_map)
+            if (integral == 0.d0) cycle
+            icount += 1_8
+            eri_buffer(icount) = integral
+            eri_index(1,icount) = i
+            eri_index(2,icount) = j
+            eri_index(3,icount) = k
+            eri_index(4,icount) = l
+            if (icount == BUFSIZE) then
+              rc = trexio_write_ao_2e_int_eri(f, offset, icount, eri_index, eri_buffer)
+              call check_success(rc)
+              icount = 0_8
+              offset += icount
+            end if
+          end do
+        end do
+      end do
+    end do
+
+    if (icount >= 0_8) then
+      rc = trexio_write_ao_2e_int_eri(f, offset, icount, eri_index, eri_buffer)
+      call check_success(rc)
+    end if
+  end if
+
 
 ! Molecular orbitals
 ! ------------------
 
+  print *, 'MOs'
+
 !  rc = trexio_write_mo_type(f, mo_label)
 !  call check_success(rc)
 
@@ -221,6 +286,8 @@ program export_trexio
 ! One-e MO integrals
 ! ------------------
 
+  print *, 'MO integrals'
+
   rc = trexio_write_mo_1e_int_kinetic(f,mo_kinetic_integrals)
   call check_success(rc)
 
@@ -234,22 +301,103 @@ program export_trexio
     rc = trexio_write_mo_1e_int_ecp_non_local(f,mo_pseudo_integrals_non_local)
     call check_success(rc)
   endif
-!
+
   rc = trexio_write_mo_1e_int_core_hamiltonian(f,one_e_dm_mo)
   call check_success(rc)
 
+! Two-e MO integrals
+! ------------------
 
-! RDM
-! ----
+  if (export_mo_ints) then
+    PROVIDE mo_two_e_integrals_in_map
 
-!  rc = trexio_write_rdm_one_e(f,one_e_dm_mo)
-!  call check_success(rc)
-!
-!  rc = trexio_write_rdm_one_e_up(f,one_e_dm_mo_alpha_average)
-!  call check_success(rc)
-!
-!  rc = trexio_write_rdm_one_e_dn(f,one_e_dm_mo_beta_average)
-!  call check_success(rc)
+    double precision, external :: mo_two_e_integral
+
+
+    icount = 0_8
+    offset = 0_8
+    do l=1,mo_num
+      do k=1,mo_num
+        do j=l,mo_num
+          do i=k,mo_num
+            if (i==j .and. k<l) cycle
+            if (i<j) cycle
+            integral = mo_two_e_integral(i,j,k,l)
+            if (integral == 0.d0) cycle
+            icount += 1_8
+            eri_buffer(icount) = integral
+            eri_index(1,icount) = i
+            eri_index(2,icount) = j
+            eri_index(3,icount) = k
+            eri_index(4,icount) = l
+            if (icount == BUFSIZE) then
+              rc = trexio_write_mo_2e_int_eri(f, offset, icount, eri_index, eri_buffer)
+              call check_success(rc)
+              icount = 0_8
+              offset += icount
+            end if
+          end do
+        end do
+      end do
+    end do
+
+    if (icount >= 0_8) then
+      rc = trexio_write_mo_2e_int_eri(f, offset, icount, eri_index, eri_buffer)
+      call check_success(rc)
+    end if
+  end if
+
+
+
+! One-e RDM
+! ---------
+
+  rc = trexio_write_rdm_1e(f,one_e_dm_mo)
+  call check_success(rc)
+
+  rc = trexio_write_rdm_1e_up(f,one_e_dm_mo_alpha_average)
+  call check_success(rc)
+
+  rc = trexio_write_rdm_1e_dn(f,one_e_dm_mo_beta_average)
+  call check_success(rc)
+
+
+! Two-e RDM
+! ---------
+
+  if (export_rdm) then
+    PROVIDE two_e_dm_mo
+
+    icount = 0_8
+    offset = 0_8
+    do l=1,mo_num
+      do k=1,mo_num
+        do j=1,mo_num
+          do i=1,mo_num
+            integral = two_e_dm_mo(i,j,k,l)
+            if (integral == 0.d0) cycle
+            icount += 1_8
+            eri_buffer(icount) = integral
+            eri_index(1,icount) = i
+            eri_index(2,icount) = j
+            eri_index(3,icount) = k
+            eri_index(4,icount) = l
+            if (icount == BUFSIZE) then
+              rc = trexio_write_rdm_2e(f, offset, icount, eri_index, eri_buffer)
+              call check_success(rc)
+              icount = 0_8
+              offset += icount
+            end if
+          end do
+        end do
+      end do
+    end do
+
+    if (icount >= 0_8) then
+      rc = trexio_write_rdm_2e(f, offset, icount, eri_index, eri_buffer)
+      call check_success(rc)
+    end if
+  end if
 
 
 ! ------------------------------------------------------------------------------

stable/champ/NEED

@@ -1 +1 @@
-determinants fci
+determinants fci csf

stable/champ/save_for_champ.irp.f

@@ -1,4 +1,4 @@
-program qmcpack
+program save_for_camp
   implicit none
   BEGIN_DOC
 ! Generates a file for CHAMP
@@ -23,6 +23,7 @@ program qmcpack
   call system('rm '//trim(ezfio_filename)//'/mo_basis/ao_md5')
   call system('$QP_ROOT/src/champ/qp_convert.py '//trim(ezfio_filename))
 
+  call write_champ_csf()
 !  integer :: iunit
 !  integer, external :: getUnitAndOpen
 !  iunit = getUnitAndOpen(trim(ezfio_filename)//'.H','w')

stable/champ/write_champ_csf.irp.f

@@ -0,0 +1,99 @@
+subroutine write_champ_csf
+  implicit none
+
+  integer, parameter             :: istate=1
+  character*(2048)               :: format
+  integer                        :: i, n_elements, j, k
+  integer, allocatable           :: list(:,:)
+
+  integer                        :: startdet, enddet, iunit
+  integer                        :: ndetI, bfIcfg, s
+  double precision               :: phasedet
+  integer                        :: idx
+
+  open(newunit=iunit, file='det.'//ezfio_filename)
+
+  write(iunit, '(A, X, I4,  X, A, X, I4)') '&electrons nelec', elec_num, 'nup', elec_alpha_num
+
+
+  ! Determinants
+  write(iunit, '(A, I10)') 'determinants', N_det
+
+  write(format,*) '(', N_det, '(F12.8, X))'
+  write(iunit, format) psi_coef(1:N_det,istate)
+
+  write(format,*) '( ', elec_alpha_num, '(I4,X), 2X, ', elec_beta_num, '(I4,X))'
+  allocate ( list(bit_kind_size,2) )
+  do i=1, N_det
+   call bitstring_to_list( psi_det(1,1,i), list(1,1), n_elements, N_int)
+   call bitstring_to_list( psi_det(1,2,i), list(1,2), n_elements, N_int)
+   write(iunit,format) list(1:elec_alpha_num,1), list(1:elec_beta_num,2)
+  end do
+  write(iunit, '(A)') 'end'
+
+  ! CSF
+  write(iunit, '(A, I10, X, I3)') 'csf', N_csf, N_states
+
+  write(format,*) '(', N_csf, '(F12.8, X))'
+  do i=1,N_states
+    write(iunit, format) psi_csf_coef(1:N_csf,i)
+  end do
+
+  write(iunit, '(A)') 'end'
+
+  ! CSF map
+  phasedet = 1.0d0
+
+  ndetI = 0
+  do i=1,N_configuration
+    startdet = psi_configuration_to_psi_det(1,i)
+    enddet = psi_configuration_to_psi_det(2,i)
+    s = 0
+    do k=1,N_int
+      if (psi_configuration(k,1,i) == 0_bit_kind) cycle
+      s = s + popcnt(psi_configuration(k,1,i))
+    enddo
+    bfIcfg = max(1,nint((binom(s,(s+1)/2)-binom(s,((s+1)/2)+1))))
+    do k=1,bfIcfg
+      do j = startdet, enddet
+        if (DetToCSFTransformationMatrix(s,k,j-startdet+1) == 0.d0) cycle
+        ndetI += 1
+      end do
+    end do
+  end do
+
+  write(iunit, '(A)') 'csfmap'
+  write(iunit, '(I10, I10, I10)') n_csf, n_det, ndetI
+
+  do i=1,N_configuration
+    startdet = psi_configuration_to_psi_det(1,i)
+    enddet = psi_configuration_to_psi_det(2,i)
+
+    s = 0
+    do k=1,N_int
+      if (psi_configuration(k,1,i) == 0_bit_kind) cycle
+      s = s + popcnt(psi_configuration(k,1,i))
+    enddo
+    bfIcfg = max(1,nint((binom(s,(s+1)/2)-binom(s,((s+1)/2)+1))))
+
+    do k=1,bfIcfg
+      ndetI = 0
+      do j = startdet, enddet
+        if (DetToCSFTransformationMatrix(s,k,j-startdet+1) == 0.d0) cycle
+        ndetI += 1
+      end do
+
+      write(iunit, '(I4)') ndetI
+      do j = startdet, enddet
+        if (DetToCSFTransformationMatrix(s,k,j-startdet+1) == 0.d0) cycle
+        idx = psi_configuration_to_psi_det_data(j)
+        call get_phase_qp_to_cfg(psi_det(1,1,idx), psi_det(1,2,idx), phasedet)
+        write(iunit, '(I10, F18.12)') idx, DetToCSFTransformationMatrix(s,k,j-startdet+1)*phasedet
+      end do
+    end do
+
+  end do
+
+  write(iunit, '(A)') 'end'
+
+end

stable/utilities/truncate_wf.irp.f

@@ -59,7 +59,7 @@ subroutine routine_s2
  do i=1,N_det
    print *,  i, real(weight_configuration(det_to_configuration(i),:)), real(sum(weight_configuration(det_to_configuration(i),:)))
  enddo
- print*, 'Min weight of the occupation pattern ?'
+ print*, 'Min weight of the configuration?'
  read(5,*) wmin
 
  ndet_max = 0