Added Davidson in MRCC

ocaml/qp_create_ezfio_from_xyz.ml

@@ -14,7 +14,7 @@ let spec =
   +> flag "m" (optional_with_default 1 int)
      ~doc:"int Spin multiplicity (2S+1) of the molecule. Default is 1."
   +> flag "p" no_arg
-     ~doc:"Using pseudo."
+     ~doc:"Using pseudopotentials"
   +> anon ("xyz_file" %: string)
 ;;
 

scripts/cache_compile.py

@@ -0,0 +1,62 @@
+#!/usr/bin/env python
+
+import os
+import sys
+import shelve
+import hashlib
+import re
+
+r = re.compile(ur'-c\s+(\S+\.[fF]90)\s+-o\s+(\S+\.o)')
+p = re.compile(ur'-I IRPF90_temp/\S*\s+')
+mod = re.compile(ur'module\s+(?P<mod>\S+).+end\s?module\s+(?P=mod)?', re.MULTILINE | re.IGNORECASE)
+
+TMPDIR="/tmp/qp_compiler/"
+
+def main():
+    # Create temp directory
+    if "qp_compiler" not in os.listdir("/tmp"):
+      os.mkdir("/tmp/qp_compiler/")
+
+    line = sys.argv[1:]
+    command = " ".join(line)
+    command_clean = p.sub('',command)
+
+    try:
+      match = r.search(command_clean)
+      input  = match.group(1)
+      output = match.group(2)
+    except:
+      os.system(command)
+      return
+    m = hashlib.md5()
+
+    # Fread : read input
+    with open(input,'r') as file:
+      fread = file.read()
+      m.update( " ".join( [ command, fread ] ))
+
+    # Md5 Key containing command + content of Fread
+    key = TMPDIR+m.hexdigest()
+    try:
+        # Try to return the content of the .o file
+        with open(key,'r') as file:
+            result = file.read()
+    except IOError:
+        # Compile the file -> .o
+        os.system(command)
+        # Read the .o
+        with open(output,'r') as file:
+            result = file.read()
+        # Copy the .o in database
+        if not mod.search(fread.replace('\n',' ')):
+            with open(key,'w') as file:
+                file.write(result)
+        else:
+            print input+' -> module'
+    else:
+        # Write the .o file
+        with open(output,'w') as file:
+            file.write(result)
+
+if __name__ == '__main__':
+    main()

scripts/install/install_m4.sh

@@ -16,6 +16,7 @@ fi
 
 cd ${QPACKAGE_ROOT}
 
+rm -f l${QPACKAGE_ROOT}/bin/m4
 if [[ -z ${M4} ]] 
 then 
     rm -f -- bin/m4

scripts/pseudo/put_pseudo_in_ezfio.py

@@ -186,7 +186,7 @@ def add_zero(array, size, type):
 
 def make_it_square(matrix, dim, type=float):
     """
-    matix the matrix to squate
+    matix the matrix to square
     dim array  [lmax, kmax]
     type the null value you want
         [[[28.59107316], [19.37583724]], [[50.25646328]]]

setup_environment.sh

@@ -28,7 +28,7 @@ EOF
 
 
 source quantum_package.rc
-
+mkdir -p install_logs
 echo "${BLUE}===== Installing IRPF90 ===== ${BLACK}"
 ${QPACKAGE_ROOT}/scripts/install/install_irpf90.sh     | tee ${QPACKAGE_ROOT}/install_logs/install_irpf90.log
 if [[ ! -d ${QPACKAGE_ROOT}/irpf90 ]]  ||  [[ ! -x ${QPACKAGE_ROOT}/bin/irpf90 ]] ||  [[ ! -x ${QPACKAGE_ROOT}/bin/irpman ]]

src/Bitmask/bitmasks_module.f90

@@ -8,4 +8,4 @@ module bitmasks
   integer, parameter             :: d_part2 = 4
   integer, parameter             :: s_hole  = 5
   integer, parameter             :: s_part  = 6
-end module
+end module bitmasks

src/DensityFit/NEEDED_CHILDREN_MODULES

@@ -1 +1 @@
-AOs
+AOs Pseudo

src/DensityFit/README.rst

@@ -10,6 +10,51 @@ Documentation
 .. Do not edit this section. It was auto-generated from the
 .. NEEDED_MODULES file.
 
+`aux_basis_coef <http://github.com/LCPQ/quantum_package/tree/master/src/DensityFit/aux_basis.irp.f#L94>`_
+  Exponents and coefficients of the auxiliary basis
+
+`aux_basis_coef_transp <http://github.com/LCPQ/quantum_package/tree/master/src/DensityFit/aux_basis.irp.f#L37>`_
+  Exponents of the auxiliary basis
+
+`aux_basis_expo <http://github.com/LCPQ/quantum_package/tree/master/src/DensityFit/aux_basis.irp.f#L93>`_
+  Exponents and coefficients of the auxiliary basis
+
+`aux_basis_expo_transp <http://github.com/LCPQ/quantum_package/tree/master/src/DensityFit/aux_basis.irp.f#L36>`_
+  Exponents of the auxiliary basis
+
+`aux_basis_idx <http://github.com/LCPQ/quantum_package/tree/master/src/DensityFit/aux_basis.irp.f#L20>`_
+  aux_basis_idx(k) -> i,j
+
+`aux_basis_nucl <http://github.com/LCPQ/quantum_package/tree/master/src/DensityFit/aux_basis.irp.f#L40>`_
+  Exponents of the auxiliary basis
+
+`aux_basis_num <http://github.com/LCPQ/quantum_package/tree/master/src/DensityFit/aux_basis.irp.f#L2>`_
+  Number of auxiliary basis functions
+
+`aux_basis_num_8 <http://github.com/LCPQ/quantum_package/tree/master/src/DensityFit/aux_basis.irp.f#L3>`_
+  Number of auxiliary basis functions
+
+`aux_basis_num_sqrt <http://github.com/LCPQ/quantum_package/tree/master/src/DensityFit/aux_basis.irp.f#L1>`_
+  Number of auxiliary basis functions
+
+`aux_basis_overlap_matrix <http://github.com/LCPQ/quantum_package/tree/master/src/DensityFit/aux_basis.irp.f#L69>`_
+  Auxiliary basis set
+
+`aux_basis_power <http://github.com/LCPQ/quantum_package/tree/master/src/DensityFit/aux_basis.irp.f#L39>`_
+  Exponents of the auxiliary basis
+
+`aux_basis_prim_num <http://github.com/LCPQ/quantum_package/tree/master/src/DensityFit/aux_basis.irp.f#L38>`_
+  Exponents of the auxiliary basis
+
+`aux_basis_prim_num_max <http://github.com/LCPQ/quantum_package/tree/master/src/DensityFit/aux_basis.irp.f#L111>`_
+  = ao_prim_num_max
+
+`save_aux_basis <http://github.com/LCPQ/quantum_package/tree/master/src/DensityFit/aux_basis.irp.f#L120>`_
+  Undocumented
+
+`aux_basis_four_overlap <http://github.com/LCPQ/quantum_package/tree/master/src/DensityFit/overlap.irp.f#L1>`_
+  \int \chi_i(r) \chi_j(r) \chi_k(r) \chi_l(r) dr
+
 
 
 Needed Modules
@@ -21,4 +66,5 @@ Needed Modules
 .. image:: tree_dependancy.png
 
 * `AOs <http://github.com/LCPQ/quantum_package/tree/master/src/AOs>`_
+* `Pseudo <http://github.com/LCPQ/quantum_package/tree/master/src/Pseudo>`_
 

src/DensityFit/aux_basis.ezfio_config

@@ -0,0 +1,12 @@
+aux_basis
+  aux_basis_num             integer
+  aux_basis_num_sqrt        integer
+  aux_basis_idx             integer       (2,aux_basis_aux_basis_num)
+  aux_basis_prim_num        integer       (aux_basis_aux_basis_num_sqrt)
+  aux_basis_nucl            integer       (aux_basis_aux_basis_num_sqrt)
+  aux_basis_power           integer       (aux_basis_aux_basis_num_sqrt,3)
+  aux_basis_prim_num_max    integer       = maxval(aux_basis_aux_basis_prim_num)
+  aux_basis_coef            double precision (aux_basis_aux_basis_num_sqrt,aux_basis_aux_basis_prim_num_max)
+  aux_basis_expo            double precision (aux_basis_aux_basis_num_sqrt,aux_basis_aux_basis_prim_num_max)
+
+

src/DensityFit/aux_basis.irp.f

@@ -1,11 +1,19 @@
- BEGIN_PROVIDER [ integer, aux_basis_num ]
+ BEGIN_PROVIDER [ integer, aux_basis_num_sqrt ]
+&BEGIN_PROVIDER [ integer, aux_basis_num ]
 &BEGIN_PROVIDER [ integer, aux_basis_num_8 ]
  implicit none
  BEGIN_DOC
  ! Number of auxiliary basis functions
  END_DOC
  integer :: align_double
- aux_basis_num = ao_num * (ao_num+1)/2
+
+ if (do_pseudo) then
+   aux_basis_num_sqrt = ao_num + ao_pseudo_num
+!   aux_basis_num_sqrt = ao_num 
+ else
+ endif
+
+ aux_basis_num = aux_basis_num_sqrt * (aux_basis_num_sqrt+1)/2 
  aux_basis_num_8 = align_double(aux_basis_num)
 END_PROVIDER
 
@@ -16,7 +24,7 @@ BEGIN_PROVIDER [ integer, aux_basis_idx, (2,aux_basis_num) ]
  END_DOC
  integer :: i,j,k
  k=0
- do j=1,ao_num
+ do j=1,aux_basis_num_sqrt
    do i=1,j
      k = k+1
      aux_basis_idx(1,k) = i
@@ -25,6 +33,38 @@ BEGIN_PROVIDER [ integer, aux_basis_idx, (2,aux_basis_num) ]
  enddo
 END_PROVIDER
 
+ BEGIN_PROVIDER [ double precision, aux_basis_expo_transp, (ao_prim_num_max_align,aux_basis_num_sqrt) ]
+&BEGIN_PROVIDER [ double precision, aux_basis_coef_transp, (ao_prim_num_max_align,aux_basis_num_sqrt) ]
+&BEGIN_PROVIDER [ integer, aux_basis_prim_num, (aux_basis_num_sqrt) ]
+&BEGIN_PROVIDER [ integer, aux_basis_power, (aux_basis_num_sqrt,3) ]
+&BEGIN_PROVIDER [ integer, aux_basis_nucl, (aux_basis_num_sqrt) ]
+ implicit none
+ BEGIN_DOC
+ ! Exponents of the auxiliary basis
+ END_DOC
+ integer                        :: i,j
+ do j=1,ao_num
+   do i=1,ao_prim_num_max
+     aux_basis_expo_transp(i,j) = ao_expo_ordered_transp(i,j)
+     aux_basis_coef_transp(i,j) = ao_coef_normalized_ordered_transp(i,j)
+   enddo
+ enddo
+ do i=1,ao_num
+   aux_basis_prim_num(i) = ao_prim_num(i)
+   aux_basis_nucl(i) = ao_nucl(i)
+   aux_basis_power(i,1:3) = ao_power(i,1:3)
+ enddo
+
+ do j=1,ao_pseudo_num
+   aux_basis_expo_transp(1,ao_num+j) = 0.5d0*pseudo_ao_expo(j)
+   aux_basis_coef_transp(1,ao_num+j) = 1.d0
+   aux_basis_power(ao_num+j,1:3) = 0
+   aux_basis_prim_num(ao_num+j) = 1
+   aux_basis_nucl(ao_num+j) = pseudo_ao_nucl(j)
+ enddo
+  
+END_PROVIDER
+
 
 BEGIN_PROVIDER [ double precision, aux_basis_overlap_matrix, (aux_basis_num_8,aux_basis_num) ]
  implicit none
@@ -32,9 +72,9 @@ BEGIN_PROVIDER [ double precision, aux_basis_overlap_matrix, (aux_basis_num_8,au
 ! Auxiliary basis set
  END_DOC
  integer :: m,n,i,j,k,l
- double precision :: ao_four_overlap
+ double precision :: aux_basis_four_overlap
 
- aux_basis_overlap_matrix(1,1) = ao_four_overlap(1,1,1,1)
+ aux_basis_overlap_matrix(1,1) = aux_basis_four_overlap(1,1,1,1)
  !$OMP PARALLEL DO PRIVATE(i,j,k,l,m,n) SCHEDULE(GUIDED)
  do m=1,aux_basis_num
    i = aux_basis_idx(1,m)
@@ -42,7 +82,7 @@ BEGIN_PROVIDER [ double precision, aux_basis_overlap_matrix, (aux_basis_num_8,au
    do n=1,m
      k = aux_basis_idx(1,n)
      l = aux_basis_idx(2,n)
-     aux_basis_overlap_matrix(m,n) = ao_four_overlap(i,j,k,l)
+     aux_basis_overlap_matrix(m,n) = aux_basis_four_overlap(i,j,k,l)
      aux_basis_overlap_matrix(n,m) = aux_basis_overlap_matrix(m,n) 
    enddo
  enddo
@@ -50,5 +90,41 @@ BEGIN_PROVIDER [ double precision, aux_basis_overlap_matrix, (aux_basis_num_8,au
 
 END_PROVIDER
 
+ BEGIN_PROVIDER [ double precision, aux_basis_expo, (aux_basis_num_sqrt,aux_basis_prim_num_max) ]
+&BEGIN_PROVIDER [ double precision, aux_basis_coef, (aux_basis_num_sqrt,aux_basis_prim_num_max) ]
+ implicit none
+ BEGIN_DOC
+ ! Exponents and coefficients of the auxiliary basis
+ END_DOC
+ integer :: i,j
+ aux_basis_expo = 0.d0
+ aux_basis_coef = 0.d0
+ do j=1,aux_basis_num_sqrt
+   do i=1,aux_basis_prim_num(j)
+     aux_basis_expo(j,i) = aux_basis_expo_transp(i,j)
+     aux_basis_coef(j,i) = aux_basis_coef_transp(i,j)
+   enddo
+ enddo
+
+END_PROVIDER
+
+BEGIN_PROVIDER [ integer, aux_basis_prim_num_max ]
+ implicit none
+ BEGIN_DOC
+ ! = ao_prim_num_max
+ END_DOC
+ aux_basis_prim_num_max = ao_prim_num_max
+END_PROVIDER
 
 
+subroutine save_aux_basis
+ implicit none
+ call ezfio_set_aux_basis_aux_basis_num(aux_basis_num)
+ call ezfio_set_aux_basis_aux_basis_num_sqrt(aux_basis_num_sqrt)
+ call ezfio_set_aux_basis_aux_basis_idx(aux_basis_idx)
+ call ezfio_set_aux_basis_aux_basis_prim_num(aux_basis_prim_num)
+ call ezfio_set_aux_basis_aux_basis_nucl(aux_basis_nucl)
+ call ezfio_set_aux_basis_aux_basis_power(aux_basis_power)
+ call ezfio_set_aux_basis_aux_basis_coef(aux_basis_coef)
+ call ezfio_set_aux_basis_aux_basis_expo(aux_basis_expo)
+end

src/DensityFit/overlap.irp.f

@@ -1,4 +1,4 @@
-double precision function ao_four_overlap(i,j,k,l)
+double precision function aux_basis_four_overlap(i,j,k,l)
   implicit none
   BEGIN_DOC
 ! \int \chi_i(r) \chi_j(r) \chi_k(r) \chi_l(r) dr
@@ -15,44 +15,44 @@ double precision function ao_four_overlap(i,j,k,l)
   
   dim1 = n_pt_max_integrals
   
-  num_i = ao_nucl(i)
-  num_j = ao_nucl(j)
-  num_k = ao_nucl(k)
-  num_l = ao_nucl(l)
-  ao_four_overlap = 0.d0
+  num_i = aux_basis_nucl(i)
+  num_j = aux_basis_nucl(j)
+  num_k = aux_basis_nucl(k)
+  num_l = aux_basis_nucl(l)
+  aux_basis_four_overlap = 0.d0
   
   do p = 1, 3
-    I_power(p) = ao_power(i,p)
-    J_power(p) = ao_power(j,p)
-    K_power(p) = ao_power(k,p)
-    L_power(p) = ao_power(l,p)
+    I_power(p) = aux_basis_power(i,p)
+    J_power(p) = aux_basis_power(j,p)
+    K_power(p) = aux_basis_power(k,p)
+    L_power(p) = aux_basis_power(l,p)
     I_center(p) = nucl_coord(num_i,p)
     J_center(p) = nucl_coord(num_j,p)
     K_center(p) = nucl_coord(num_k,p)
     L_center(p) = nucl_coord(num_l,p)
   enddo
   
-  do p = 1, ao_prim_num(i)
+  do p = 1, aux_basis_prim_num(i)
     double precision               :: coef1
-    coef1 = ao_coef_normalized_ordered_transp(p,i)
-    do q = 1, ao_prim_num(j)
+    coef1 = aux_basis_coef_transp(p,i)
+    do q = 1, aux_basis_prim_num(j)
       call give_explicit_poly_and_gaussian(P_new,P_center,pp,fact_p,iorder_p,&
-          ao_expo_ordered_transp(p,i),ao_expo_ordered_transp(q,j),                 &
+          aux_basis_expo_transp(p,i),aux_basis_expo_transp(q,j),                 &
           I_power,J_power,I_center,J_center,dim1)
       double precision               :: coef2
-      coef2 = coef1*ao_coef_normalized_ordered_transp(q,j)*fact_p
-      do r = 1, ao_prim_num(k)
+      coef2 = coef1*aux_basis_coef_transp(q,j)*fact_p
+      do r = 1, aux_basis_prim_num(k)
         double precision               :: coef3
-        coef3 = coef2*ao_coef_normalized_ordered_transp(r,k)
-        do s = 1, ao_prim_num(l)
+        coef3 = coef2*aux_basis_coef_transp(r,k)
+        do s = 1, aux_basis_prim_num(l)
           double precision               :: general_primitive_integral
-          call give_explicit_poly_and_gaussian(Q_new,Q_center,qq,fact_q,iorder_q,&
-              ao_expo_ordered_transp(r,k),ao_expo_ordered_transp(s,l),             &
+          call give_explicit_poly_and_gaussian(Q_new,Q_center,qq,fact_q,iorder_q,    &
+              aux_basis_expo_transp(r,k),aux_basis_expo_transp(s,l),             &
               K_power,L_power,K_center,L_center,dim1)
           double precision               :: coef4
-          coef4 = coef3*ao_coef_normalized_ordered_transp(s,l)*fact_q
+          coef4 = coef3*aux_basis_coef_transp(s,l)*fact_q
           call overlap_gaussian_xyz(P_center,Q_center,pp,qq,iorder_p,iorder_q,overlap_x,overlap_y,overlap_z,overlap,dim1)
-          ao_four_overlap +=  coef4 * overlap 
+          aux_basis_four_overlap +=  coef4 * overlap 
         enddo ! s
       enddo  ! r
     enddo   ! q

src/Determinants/README.rst

@@ -532,12 +532,6 @@ Documentation
 `save_casino <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/save_for_casino.irp.f#L1>`_
   Undocumented
 
-`save_dets_qmcchem <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/save_for_qmcchem.irp.f#L1>`_
-  Undocumented
-
-`save_for_qmc <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/save_for_qmcchem.irp.f#L46>`_
-  Undocumented
-
 `save_natorb <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/save_natorb.irp.f#L1>`_
   Undocumented
 

src/Determinants/save_for_qmcchem.irp.f

@@ -1,51 +0,0 @@
-subroutine save_dets_qmcchem
- use bitmasks
- implicit none
- character :: c(mo_tot_num)
- integer :: i,k
-
- integer, allocatable :: occ(:,:,:), occ_tmp(:,:)
- !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: occ, occ_tmp
-
- call ezfio_set_determinants_det_num(N_det)
- call ezfio_set_determinants_det_coef(psi_coef_sorted(1,1))
-
- allocate (occ(elec_alpha_num,N_det,2))
- ! OMP PARALLEL DEFAULT(NONE) &
- ! OMP PRIVATE(occ_tmp,i,k)&
- ! OMP SHARED(N_det,psi_det_sorted,elec_alpha_num, &
- ! OMP   occ,elec_beta_num,N_int)
- allocate (occ_tmp(N_int*bit_kind_size,2))
- occ_tmp = 0
- ! OMP DO 
- do i=1,N_det
-  call bitstring_to_list(psi_det_sorted(1,1,i), occ_tmp(1,1), elec_alpha_num, N_int )
-  call bitstring_to_list(psi_det_sorted(1,2,i), occ_tmp(1,2), elec_beta_num, N_int )
-  do k=1,elec_alpha_num
-    occ(k,i,1) = occ_tmp(k,1)
-    occ(k,i,2) = occ_tmp(k,2)
-  enddo
- enddo
- ! OMP END DO
- deallocate(occ_tmp)
- ! OMP END PARALLEL
- call ezfio_set_determinants_det_occ(occ)
- call write_int(output_determinants,N_det,'Determinants saved for QMC')
- deallocate(occ)
- open(unit=31,file=trim(ezfio_filename)//'/mo_basis/mo_classif')
- write(31,'(I1)') 1
- write(31,*) mo_tot_num
- do i=1,mo_tot_num
-   write(31,'(A)') 'a'
- enddo
- close(31)
- call system('gzip -f '//trim(ezfio_filename)//'/mo_basis/mo_classif')
-
-end
-
-program save_for_qmc
- read_wf = .True.
- TOUCH read_wf
-!  call save_dets_qmcchem
- call write_spindeterminants
-end

src/Integrals_Monoelec/README.rst

@@ -97,18 +97,12 @@ Documentation
 `ao_pseudo_integral <http://github.com/LCPQ/quantum_package/tree/master/src/Integrals_Monoelec/pot_ao_pseudo_ints.irp.f#L1>`_
   Pseudo-potential
 
-`ao_pseudo_integral_local <http://github.com/LCPQ/quantum_package/tree/master/src/Integrals_Monoelec/pot_ao_pseudo_ints.irp.f#L13>`_
+`ao_pseudo_integral_local <http://github.com/LCPQ/quantum_package/tree/master/src/Integrals_Monoelec/pot_ao_pseudo_ints.irp.f#L15>`_
   Local pseudo-potential
 
-`ao_pseudo_integral_non_local <http://github.com/LCPQ/quantum_package/tree/master/src/Integrals_Monoelec/pot_ao_pseudo_ints.irp.f#L119>`_
+`ao_pseudo_integral_non_local <http://github.com/LCPQ/quantum_package/tree/master/src/Integrals_Monoelec/pot_ao_pseudo_ints.irp.f#L121>`_
   Local pseudo-potential
 
-`pseudo_grid <http://github.com/LCPQ/quantum_package/tree/master/src/Integrals_Monoelec/pot_ao_pseudo_ints.irp.f#L223>`_
-  Grid points for f(|r-r_A|) = \int Y_{lm}^{C} (|r-r_C|, \Omega_C) \chi_i^{A} (r-r_A) d\Omega_C
-  .br
-  <img src="http://latex.codecogs.com/gif.latex?f(|r-r_A|)&space;=&space;\int&space;Y_{lm}^{C}&space;(|r-r_C|,&space;\Omega_C)&space;\chi_i^{A}&space;(r-r_A)&space;d\Omega_C"
-  title="f(|r-r_A|) = \int Y_{lm}^{C} (|r-r_C|, \Omega_C) \chi_i^{A} (r-r_A) d\Omega_C" />
-
 `mo_nucl_elec_integral <http://github.com/LCPQ/quantum_package/tree/master/src/Integrals_Monoelec/pot_mo_ints.irp.f#L1>`_
   interaction nuclear electron on the MO basis
 

src/Integrals_Monoelec/pot_ao_pseudo_ints.irp.f

@@ -4,7 +4,9 @@ BEGIN_PROVIDER [ double precision, ao_pseudo_integral, (ao_num_align,ao_num)]
 ! Pseudo-potential
   END_DOC
   if (do_pseudo) then
-    ao_pseudo_integral = ao_pseudo_integral_local !+ ao_pseudo_integral_non_local 
+    ao_pseudo_integral = ao_pseudo_integral_local + ao_pseudo_integral_non_local 
+!    ao_pseudo_integral = ao_pseudo_integral_local 
+!    ao_pseudo_integral = ao_pseudo_integral_non_local 
   else
     ao_pseudo_integral = 0.d0
   endif
@@ -220,67 +222,5 @@ END_PROVIDER
 
  END_PROVIDER
 
-BEGIN_PROVIDER [ double precision, pseudo_grid, (pseudo_grid_size,ao_num,-pseudo_lmax:pseudo_lmax,0:pseudo_lmax,nucl_num)  ]
- implicit none
- BEGIN_DOC
-! Grid points for f(|r-r_A|) = \int Y_{lm}^{C} (|r-r_C|, \Omega_C) \chi_i^{A} (r-r_A) d\Omega_C
-!
-! <img src="http://latex.codecogs.com/gif.latex?f(|r-r_A|)&space;=&space;\int&space;Y_{lm}^{C}&space;(|r-r_C|,&space;\Omega_C)&space;\chi_i^{A}&space;(r-r_A)&space;d\Omega_C"
-! title="f(|r-r_A|) = \int Y_{lm}^{C} (|r-r_C|, \Omega_C) \chi_i^{A} (r-r_A) d\Omega_C" />
- END_DOC
- ! l,m    : Y(l,m) parameters
- ! c(3)   : pseudopotential center
- ! a(3)   : Atomic Orbital center
- ! n_a(3) : Powers of x,y,z in the Atomic Orbital
- ! g_a    : Atomic Orbital exponent      
- ! r      : Distance between the Atomic Orbital center and the considered point
- double precision, external :: ylm_orb
- integer :: n_a(3)
- double precision :: a(3), c(3), g_a
- integer :: i,j,k,l,m,n,p
- double precision :: r(pseudo_grid_size), dr, Ulc
- double precision, parameter :: rmax= 10.d0
-
- dr = rmax/dble(pseudo_grid_size)
- r(1) = 0.d0
- do j=2,pseudo_grid_size
-   r(j) = r(j-1) + dr
- enddo
-
- pseudo_grid = 0.d0
- do k=1,nucl_num
-   c(1:3) = nucl_coord(k,1:3)
-   do l=0,pseudo_lmax
-     do i=1,ao_num
-       a(1:3) = nucl_coord(ao_nucl(i),1:3)
-       n_a(1:3) = ao_power(i,1:3)
-       do j=1,pseudo_grid_size
-         do p=1,ao_prim_num(i)
-           g_a = ao_expo_ordered_transp(p,i)
-           do m=-l,l
-             double precision               :: y
-             !            y = ylm_orb(l,m,c,a,n_a,g_a,r(j))
-             !            if (y /= 0.d0) then
-             !              print *, 'y = ', y
-             !              print *, 'l = ', l
-             !              print *, 'm = ', m
-             !              print *, 'c = ', c
-             !              print *, 'a = ', a
-             !              print *, 'n = ', n_a
-             !              print *, 'g = ', g_a
-             !              print *, 'r = ', r(j)
-             !              print *,  ''
-             !            endif
-             y = ylm_orb(l,m,c,a,n_a,g_a,r(j))
-             pseudo_grid(j,i,m,l,k) = pseudo_grid(j,i,m,l,k) +           &
-                 ao_coef_normalized_ordered_transp(p,i)*y
-           enddo
-         enddo
-       enddo
-     enddo
-   enddo
- enddo
-
-END_PROVIDER
 
 

src/Integrals_Monoelec/pot_mo_pseudo_ints.irp.f

@@ -28,7 +28,6 @@ BEGIN_PROVIDER [double precision, mo_pseudo_integral, (mo_tot_num_align,mo_tot_n
     enddo
   enddo
   !$OMP END PARALLEL DO
-  call ezfio_set_pseudo_pseudo_matrix(mo_pseudo_integral(1:mo_tot_num,1:mo_tot_num))
 END_PROVIDER
 
 

src/Integrals_Monoelec/pseudopot.f90

@@ -1485,7 +1485,6 @@ end
         a=bessel_mod_exp(n,x)
         return
        endif
-       print *,  n,x
        if(n.eq.0)a=dsinh(x)/x
        if(n.eq.1)a=(x*dcosh(x)-dsinh(x))/x**2
        if(n.ge.2)a=bessel_mod_recur(n-2,x)-(2*n-1)/x*bessel_mod_recur(n-1,x)

src/MRCC/README.rst

@@ -53,7 +53,7 @@ Documentation
 `ci_electronic_energy_dressed <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC/mrcc_utils.irp.f#L78>`_
   Eigenvectors/values of the CI matrix
 
-`ci_energy_dressed <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC/mrcc_utils.irp.f#L144>`_
+`ci_energy_dressed <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC/mrcc_utils.irp.f#L145>`_
   N_states lowest eigenvalues of the dressed CI matrix
 
 `delta_ij <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC/mrcc_utils.irp.f#L43>`_
@@ -62,7 +62,7 @@ Documentation
 `delta_ij_non_cas <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC/mrcc_utils.irp.f#L34>`_
   Dressing matrix in SD basis
 
-`diagonalize_ci_dressed <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC/mrcc_utils.irp.f#L159>`_
+`diagonalize_ci_dressed <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC/mrcc_utils.irp.f#L160>`_
   Replace the coefficients of the CI states by the coefficients of the
   eigenstates of the CI matrix
 

src/MRCC/mrcc_utils.irp.f

@@ -92,9 +92,10 @@ END_PROVIDER
   
   if (diag_algorithm == "Davidson") then
     
-    stop 'use Lapack'
-!    call davidson_diag(psi_det,CI_eigenvectors_dressed,CI_electronic_energy_dressed, &
-!        size(CI_eigenvectors_dressed,1),N_det,N_states_diag,N_int,output_determinants)
+    integer :: istate
+    istate = 1
+    call davidson_diag_mrcc(psi_det,CI_eigenvectors_dressed,CI_electronic_energy_dressed, &
+        size(CI_eigenvectors_dressed,1),N_det,N_states_diag,N_int,output_determinants,istate)
     
   else if (diag_algorithm == "Lapack") then
     

src/Makefile.config.example

@@ -4,6 +4,7 @@ DEBUG  = 0
 
 IRPF90_FLAGS+= --align=32
 FC     = ifort -g 
+#FC     = cache_compile.py ifort -g    # Accelerates compilation
 FCFLAGS= 
 FCFLAGS+= -xHost
 #FCFLAGS+= -xAVX

src/NEEDED_MODULES

@@ -1 +1 @@
-AOs Integrals_Bielec Bitmask CAS_SD CID CID_SC2_selected CID_selected CIS CISD CISD_SC2_selected CISD_selected DDCI_selected Determinants Electrons Ezfio_files FCIdump Full_CI Generators_CAS Generators_full Hartree_Fock MOGuess Molden Integrals_Monoelec MOs MP2 MRCC Nuclei Pseudo Selectors_full Utils DensityFit
+AOs Integrals_Bielec Bitmask CAS_SD CID CID_SC2_selected CID_selected CIS CISD CISD_SC2_selected CISD_selected DDCI_selected Determinants Electrons Ezfio_files FCIdump Full_CI Generators_CAS Generators_full Hartree_Fock MOGuess Molden Integrals_Monoelec MOs MP2 MRCC Nuclei Pseudo Selectors_full Utils DensityFit QmcChem

src/Nuclei/NEEDED_CHILDREN_MODULES

@@ -1 +1 @@
-Ezfio_files
+Ezfio_files Utils

src/Nuclei/README.rst

@@ -15,6 +15,7 @@ Needed Modules
 .. image:: tree_dependancy.png
 
 * `Ezfio_files <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files>`_
+* `Utils <http://github.com/LCPQ/quantum_package/tree/master/src/Utils>`_
 
 Documentation
 =============

src/Pseudo/EZFIO.cfg

@@ -57,10 +57,16 @@ default: False
 
 [pseudo_grid_size]
 type: integer
-doc: Size of the QMC grid 
+doc: Nb of points of the QMC grid 
 interface: input
 default: 100
 
+[pseudo_grid_rmax]
+type: double precision
+doc: R_maxof the QMC grid 
+interface: input
+default: 4.0
+
 [pseudo_grid]
 type: double precision
 doc: QMC grid 
@@ -71,6 +77,6 @@ size: (pseudo.pseudo_grid_size,ao_basis.ao_num,-pseudo.pseudo_lmax:pseudo.pseudo
 type: double precision
 doc: QMC grid 
 interface: output
-size: (mo_basis.mo_tot_num,mo_basis.mo_tot_num)
+size: (aux_basis.aux_basis_num_sqrt,aux_basis.aux_basis_num_sqrt)
 
 

src/Pseudo/NEEDED_CHILDREN_MODULES

@@ -1 +1 @@
-Ezfio_files
+Ezfio_files Nuclei

src/Pseudo/README.rst

@@ -19,4 +19,5 @@ Needed Modules
 .. image:: tree_dependancy.png
 
 * `Ezfio_files <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files>`_
+* `Nuclei <http://github.com/LCPQ/quantum_package/tree/master/src/Nuclei>`_
 

src/QmcChem/Makefile

@@ -0,0 +1,6 @@
+# Define here all new external source files and objects.Don't forget to prefix the
+# object files with IRPF90_temp/
+SRC=
+OBJ=
+
+include $(QPACKAGE_ROOT)/src/Makefile.common

src/QmcChem/NEEDED_CHILDREN_MODULES

@@ -0,0 +1 @@
+Determinants DensityFit

src/QmcChem/README.rst

@@ -0,0 +1,23 @@
+==============
+QmcChem Module
+==============
+
+Documentation
+=============
+
+.. Do not edit this section. It was auto-generated from the
+.. NEEDED_MODULES file.
+
+
+
+Needed Modules
+==============
+
+.. Do not edit this section. It was auto-generated from the
+.. NEEDED_MODULES file.
+
+.. image:: tree_dependancy.png
+
+* `Determinants <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants>`_
+* `DensityFit <http://github.com/LCPQ/quantum_package/tree/master/src/DensityFit>`_
+

src/QmcChem/pot_ao_pseudo_ints.irp.f

@@ -0,0 +1,275 @@
+BEGIN_PROVIDER [ double precision, aux_pseudo_integral, (aux_basis_num_sqrt,aux_basis_num_sqrt)]
+  implicit none
+  BEGIN_DOC
+! Pseudo-potential
+  END_DOC
+  if (do_pseudo) then
+!    aux_pseudo_integral = aux_pseudo_integral_local + aux_pseudo_integral_non_local 
+!    aux_pseudo_integral = aux_pseudo_integral_local 
+    aux_pseudo_integral = aux_pseudo_integral_non_local 
+  else
+    aux_pseudo_integral = 0.d0
+  endif
+END_PROVIDER
+
+ BEGIN_PROVIDER [ double precision, aux_pseudo_integral_local, (aux_basis_num_sqrt,aux_basis_num_sqrt)]
+  implicit none
+  BEGIN_DOC
+! Local pseudo-potential
+  END_DOC
+ double precision  :: alpha, beta, gama, delta
+ integer           :: num_A,num_B
+ double precision  :: A_center(3),B_center(3),C_center(3)
+ integer           :: power_A(3),power_B(3)
+ integer           :: i,j,k,l,n_pt_in,m
+ double precision  ::  Vloc, Vpseudo
+
+ double precision  :: cpu_1, cpu_2, wall_1, wall_2, wall_0
+ integer           :: thread_num
+
+  aux_pseudo_integral_local = 0.d0
+
+  !! Dump array 
+  integer, allocatable ::  n_k_dump(:)
+  double precision, allocatable ::  v_k_dump(:), dz_k_dump(:) 
+ 
+  allocate(n_k_dump(1:pseudo_klocmax), v_k_dump(1:pseudo_klocmax), dz_k_dump(1:pseudo_klocmax)) 
+
+
+  !  _                
+  ! /   _. |  _     | 
+  ! \_ (_| | (_ |_| | 
+  !                   
+
+  print*, 'Providing the nuclear electron pseudo integrals '
+
+  call wall_time(wall_1)
+  call cpu_time(cpu_1)
+
+  !$OMP PARALLEL &
+  !$OMP DEFAULT (NONE) &
+  !$OMP PRIVATE (i,j,k,l,m,alpha,beta,A_center,B_center,C_center,power_A,power_B, &
+  !$OMP          num_A,num_B,Z,c,n_pt_in, &
+  !$OMP          v_k_dump,n_k_dump, dz_k_dump, &
+  !$OMP          wall_0,wall_2,thread_num) & 
+  !$OMP SHARED (aux_basis_num_sqrt,aux_basis_prim_num,aux_basis_expo_transp,aux_basis_power,aux_basis_nucl,nucl_coord,aux_basis_coef_transp, &
+  !$OMP         aux_pseudo_integral_local,nucl_num,nucl_charge, &
+  !$OMP         pseudo_klocmax,pseudo_lmax,pseudo_kmax,pseudo_v_k,pseudo_n_k, pseudo_dz_k, &
+  !$OMP         wall_1)
+  
+  !$OMP DO SCHEDULE (guided)
+
+  do j = 1, aux_basis_num_sqrt
+
+   num_A = aux_basis_nucl(j)
+   power_A(1:3)= aux_basis_power(j,1:3)
+   A_center(1:3) = nucl_coord(num_A,1:3)
+
+  do i = 1, aux_basis_num_sqrt
+
+   num_B = aux_basis_nucl(i)
+   power_B(1:3)= aux_basis_power(i,1:3)
+   B_center(1:3) = nucl_coord(num_B,1:3)
+
+    do l=1,aux_basis_prim_num(j)
+     alpha = aux_basis_expo_transp(l,j)
+  
+    do m=1,aux_basis_prim_num(i)
+      beta = aux_basis_expo_transp(m,i)
+      double precision :: c
+      c = 0.d0
+       
+      do  k = 1, nucl_num
+        double precision :: Z
+        Z = nucl_charge(k)
+  
+        C_center(1:3) = nucl_coord(k,1:3)
+  
+        v_k_dump = pseudo_v_k(k,1:pseudo_klocmax)
+        n_k_dump = pseudo_n_k(k,1:pseudo_klocmax)
+        dz_k_dump =  pseudo_dz_k(k,1:pseudo_klocmax)
+
+        c = c + Vloc(pseudo_klocmax, v_k_dump,n_k_dump, dz_k_dump, &
+                     A_center,power_A,alpha,B_center,power_B,beta,C_center)
+  
+      enddo
+      aux_pseudo_integral_local(i,j) = aux_pseudo_integral_local(i,j) + &
+                                   aux_basis_coef_transp(l,j)*aux_basis_coef_transp(m,i)*c
+     enddo
+     enddo
+  enddo
+
+    call wall_time(wall_2)
+    if (thread_num == 0) then
+      if (wall_2 - wall_0 > 1.d0) then
+        wall_0 = wall_2
+        print*, 100.*float(j)/float(aux_basis_num_sqrt), '% in ',  &
+                                 wall_2-wall_1, 's'
+      endif
+    endif
+  enddo
+
+ !$OMP END DO
+ !$OMP END PARALLEL
+
+
+  deallocate(n_k_dump,v_k_dump, dz_k_dump)
+
+ END_PROVIDER
+
+
+ BEGIN_PROVIDER [ double precision, aux_pseudo_integral_non_local, (aux_basis_num_sqrt,aux_basis_num_sqrt)]
+  implicit none
+  BEGIN_DOC
+! Local pseudo-potential
+  END_DOC
+ double precision  :: alpha, beta, gama, delta
+ integer           :: num_A,num_B
+ double precision  :: A_center(3),B_center(3),C_center(3)
+ integer           :: power_A(3),power_B(3)
+ integer           :: i,j,k,l,n_pt_in,m
+ double precision  ::  Vloc, Vpseudo
+
+ double precision  :: cpu_1, cpu_2, wall_1, wall_2, wall_0
+ integer           :: thread_num
+
+  aux_pseudo_integral_non_local = 0.d0
+
+  !! Dump array 
+  integer, allocatable ::  n_kl_dump(:,:)
+  double precision, allocatable ::  v_kl_dump(:,:), dz_kl_dump(:,:) 
+ 
+  allocate(n_kl_dump(pseudo_kmax,0:pseudo_lmax), v_kl_dump(pseudo_kmax,0:pseudo_lmax), dz_kl_dump(pseudo_kmax,0:pseudo_lmax)) 
+
+  !  _                
+  ! /   _. |  _     | 
+  ! \_ (_| | (_ |_| | 
+  !                   
+
+  print*, 'Providing the nuclear electron pseudo integrals '
+
+  call wall_time(wall_1)
+  call cpu_time(cpu_1)
+
+  !$OMP PARALLEL &
+  !$OMP DEFAULT (NONE) &
+  !$OMP PRIVATE (i,j,k,l,m,alpha,beta,A_center,B_center,C_center,power_A,power_B, &
+  !$OMP          num_A,num_B,Z,c,n_pt_in, &
+  !$OMP          n_kl_dump, v_kl_dump, dz_kl_dump, &
+  !$OMP          wall_0,wall_2,thread_num) & 
+  !$OMP SHARED (aux_basis_num_sqrt,aux_basis_prim_num,aux_basis_expo_transp,aux_basis_power,aux_basis_nucl,nucl_coord,aux_basis_coef_transp, &
+  !$OMP         aux_pseudo_integral_non_local,nucl_num,nucl_charge, &
+  !$OMP         pseudo_klocmax,pseudo_lmax,pseudo_kmax,pseudo_n_kl, pseudo_v_kl, pseudo_dz_kl, &
+  !$OMP         wall_1)
+  
+  !$OMP DO SCHEDULE (guided)
+
+  do j = 1, aux_basis_num_sqrt
+
+   num_A = aux_basis_nucl(j)
+   power_A(1:3)= aux_basis_power(j,1:3)
+   A_center(1:3) = nucl_coord(num_A,1:3)
+
+  do i = 1, aux_basis_num_sqrt
+
+   num_B = aux_basis_nucl(i)
+   power_B(1:3)= aux_basis_power(i,1:3)
+   B_center(1:3) = nucl_coord(num_B,1:3)
+
+    do l=1,aux_basis_prim_num(j)
+     alpha = aux_basis_expo_transp(l,j)
+  
+    do m=1,aux_basis_prim_num(i)
+      beta = aux_basis_expo_transp(m,i)
+      double precision :: c
+      c = 0.d0
+       
+      do  k = 1, nucl_num
+        double precision :: Z
+        Z = nucl_charge(k)
+  
+        C_center(1:3) = nucl_coord(k,1:3)
+  
+        n_kl_dump = pseudo_n_kl(k,1:pseudo_kmax,0:pseudo_lmax)
+        v_kl_dump = pseudo_v_kl(k,1:pseudo_kmax,0:pseudo_lmax)
+        dz_kl_dump = pseudo_dz_kl(k,1:pseudo_kmax,0:pseudo_lmax)
+
+        c = c + Vpseudo(pseudo_lmax,pseudo_kmax,v_kl_dump,n_kl_dump,dz_kl_dump,A_center,power_A,alpha,B_center,power_B,beta,C_center)
+  
+      enddo
+      aux_pseudo_integral_non_local(i,j) = aux_pseudo_integral_non_local(i,j) + &
+                                   aux_basis_coef_transp(l,j)*aux_basis_coef_transp(m,i)*c
+     enddo
+     enddo
+  enddo
+
+    call wall_time(wall_2)
+    if (thread_num == 0) then
+      if (wall_2 - wall_0 > 1.d0) then
+        wall_0 = wall_2
+        print*, 100.*float(j)/float(aux_basis_num_sqrt), '% in ',  &
+                                 wall_2-wall_1, 's'
+      endif
+    endif
+  enddo
+
+ !$OMP END DO
+ !$OMP END PARALLEL
+
+
+  deallocate(n_kl_dump,v_kl_dump, dz_kl_dump)
+
+
+ END_PROVIDER
+
+BEGIN_PROVIDER [ double precision, pseudo_grid, (pseudo_grid_size,ao_num,-pseudo_lmax:pseudo_lmax,0:pseudo_lmax,nucl_num)  ]
+ implicit none
+ BEGIN_DOC
+! Grid points for f(|r-r_A|) = \int Y_{lm}^{C} (|r-r_C|, \Omega_C) \chi_i^{A} (r-r_A) d\Omega_C
+!
+! <img src="http://latex.codecogs.com/gif.latex?f(|r-r_A|)&space;=&space;\int&space;Y_{lm}^{C}&space;(|r-r_C|,&space;\Omega_C)&space;\chi_i^{A}&space;(r-r_A)&space;d\Omega_C"
+! title="f(|r-r_A|) = \int Y_{lm}^{C} (|r-r_C|, \Omega_C) \chi_i^{A} (r-r_A) d\Omega_C" />
+ END_DOC
+ ! l,m    : Y(l,m) parameters
+ ! c(3)   : pseudopotential center
+ ! a(3)   : Atomic Orbital center
+ ! n_a(3) : Powers of x,y,z in the Atomic Orbital
+ ! g_a    : Atomic Orbital exponent      
+ ! r      : Distance between the Atomic Orbital center and the considered point
+ double precision, external :: ylm_orb
+ integer :: n_a(3)
+ double precision :: a(3), c(3), g_a
+ integer :: i,j,k,l,m,n,p
+ double precision :: r(pseudo_grid_size), dr, Ulc
+ double precision :: y
+
+ dr = pseudo_grid_rmax/dble(pseudo_grid_size)
+ r(1) = 0.d0
+ do j=2,pseudo_grid_size
+   r(j) = r(j-1) + dr
+ enddo
+
+ pseudo_grid = 0.d0
+ do k=1,nucl_num
+   c(1:3) = nucl_coord(k,1:3)
+   do l=0,pseudo_lmax
+     do i=1,ao_num
+       a(1:3) = nucl_coord(ao_nucl(i),1:3)
+       n_a(1:3) = ao_power(i,1:3)
+       do j=1,pseudo_grid_size
+         do p=1,ao_prim_num(i)
+           g_a = ao_expo_ordered_transp(p,i)
+           do m=-l,l
+             y = ylm_orb(l,m,c,a,n_a,g_a,r(j))
+             pseudo_grid(j,i,m,l,k) = pseudo_grid(j,i,m,l,k) +           &
+                 ao_coef_normalized_ordered_transp(p,i)*y
+           enddo
+         enddo
+       enddo
+     enddo
+   enddo
+ enddo
+
+END_PROVIDER
+
+

src/QmcChem/pseudo.irp.f

@@ -0,0 +1,170 @@
+subroutine write_pseudopotential
+  implicit none
+  BEGIN_DOC
+!  Write the pseudo_potential into the EZFIO file
+  END_DOC
+!  call ezfio_set_pseudo_pseudo_matrix(pseudo_matrix)
+  call ezfio_set_pseudo_pseudo_grid(pseudo_grid)
+end
+
+
+BEGIN_PROVIDER [ double precision, pseudo_matrix, (aux_basis_num_sqrt,aux_basis_num_sqrt) ]
+ implicit none
+ BEGIN_DOC
+ ! Pseudo-potential expressed in the basis of ao products
+ END_DOC
+
+ integer :: i,j,k,l
+ integer :: info, m,n, lwork, lda, ldu, ldvt
+ integer, allocatable :: iwork(:)
+ character :: jobz
+ double precision, allocatable :: a(:,:),work(:)
+
+ double precision,allocatable :: U(:,:)
+ double precision,allocatable :: Vt(:,:)
+ double precision,allocatable :: S(:), B(:)
+
+
+
+ jobz = 'A'
+ m = aux_basis_num
+ n = aux_basis_num
+ lda = size(aux_basis_overlap_matrix,1)
+ ldu = lda
+ ldvt = lda
+ lwork = -1
+
+! allocate (A(lda,n), U(ldu,n), Vt(ldvt,n), S(n), work(1), b(n), iwork(8*n))
+ allocate (A(lda,n), U(ldu,n), Vt(ldvt,n), S(n), work(1), b(n),iwork(1))
+
+ work(1) = 1
+ do i=1,n
+   do j=1,n
+     A(i,j) = aux_basis_overlap_matrix(i,j)
+   enddo
+ enddo
+
+! call dgesdd(jobz, m, n, A, lda, s, u, ldu, vt, ldvt, work, lwork, iwork, info)
+ call dgesvd(jobz, jobz, m, n, a, lda, s, u, ldu, vt, ldvt, work, lwork, info)
+ lwork = int(work(1))
+ deallocate(work)
+
+ print *,  'Fitting pseudo-potentials'
+
+ allocate(work(lwork))
+! call dgesdd(jobz, m, n, A, lda, s, u, ldu, vt, ldvt, work, lwork, iwork, info)
+ call dgesvd(jobz, jobz, m, n, a, lda, s, u, ldu, vt, ldvt, work, lwork, info)
+ deallocate(work)
+
+ do i=1,n
+ print *,  i, s(i)
+ enddo
+
+ do k=1,n
+   if (s(k) < 1.d-1) then
+     s(k) = 0.d0
+   else
+     s(k) = 1.d0/s(k)
+   endif
+   do m=1,n
+     Vt(m,k) = S(m) * Vt(m,k) 
+   enddo
+ enddo
+ call dgemm('N','N',n,n,n,1.d0,U,lda,Vt,ldvt,0.d0,A,lda)
+! do k=1,n
+!  do l=1,n
+!   A(k,l) = 0.d0
+!   do m=1,n
+!    A(k,l) = A(k,l) + U(k,m) * Vt(m,l)
+!   enddo
+! enddo
+
+ do k=1,n
+   i = aux_basis_idx(1,k)
+   j = aux_basis_idx(2,k)
+   b(k) = aux_pseudo_integral(i,j)
+ enddo
+
+ do k=1,n
+   S(k) = 0.d0
+ enddo
+ do l=1,n
+   do k=1,n
+     S(k) = S(k) + A(k,l) * b(l)
+   enddo
+ enddo
+
+ do k=1,aux_basis_num
+   i = aux_basis_idx(1,k)
+   j = aux_basis_idx(2,k)
+   pseudo_matrix(i,j) = S(k)
+   pseudo_matrix(j,i) = S(k)
+ enddo
+ deallocate(a,b,s,iwork,u,vt)
+
+print *,  'Done'
+ if (info /= 0) then
+   print *,  info
+   stop 'pseudo fit failed'
+ endif
+END_PROVIDER
+
+
+
+
+
+!BEGIN_PROVIDER [ double precision, pseudo_matrix, (ao_num,ao_num) ]
+! implicit none
+! BEGIN_DOC
+! ! Pseudo-potential expressed in the basis of ao products
+! END_DOC
+!
+! integer :: i,j,k
+! integer :: info, n, lwork, lda, ldb, nrhs
+! character :: uplo
+! integer, allocatable :: ipiv(:)
+! double precision, allocatable :: a(:,:),work(:), b(:)
+!
+! uplo = 'L'
+! n = aux_basis_num
+! nrhs = 1
+! lda = size(aux_basis_overlap_matrix,1)
+! ldb = n
+! lwork = -1
+!
+! print *,  'Fitting pseudo-potentials'
+! allocate(work(1),a(lda,n),ipiv(n),b(n))
+! work(1) = 1
+! do i=1,n
+!   do j=1,n
+!     a(i,j) = aux_basis_overlap_matrix(i,j)
+!   enddo
+! enddo
+!
+! do k=1,n
+!   i = aux_basis_idx(1,k)
+!   j = aux_basis_idx(2,k)
+!   b(k) = ao_pseudo_integral(i,j)
+! enddo
+! call dsysv( uplo, n, nrhs, a, lda, ipiv, b, ldb, work, lwork, info )
+! lwork = int(work(1))
+! deallocate(work)
+!
+! allocate(work(lwork))
+! call dsysv( uplo, n, nrhs, a, lda, ipiv, b, ldb, work, lwork, info )
+! deallocate(work,ipiv)
+! do k=1,aux_basis_num
+!   i = aux_basis_idx(1,k)
+!   j = aux_basis_idx(2,k)
+!   pseudo_matrix(i,j) = b(k)
+!   pseudo_matrix(j,i) = b(k)
+! enddo
+! deallocate(a,b)
+!
+!print *,  'Done'
+! if (info /= 0) then
+!   print *,  info
+!   stop 'pseudo fit failed'
+! endif
+!END_PROVIDER
+

src/QmcChem/save_for_qmcchem.irp.f

@@ -0,0 +1,8 @@
+program save_for_qmc
+ read_wf = .True.
+ TOUCH read_wf
+!  call save_dets_qmcchem
+ call write_spindeterminants 
+! call write_pseudopotential
+! call save_aux_basis
+end