Merge branch 'dev' of github.com:QuantumPackage/qp2 into dev

README.md

@@ -1,6 +1,7 @@
 # Quantum Package 2.2
 
-<img src="https://raw.githubusercontent.com/QuantumPackage/qp2/master/data/qp2.png" width="250">
+<!--- img src="https://raw.githubusercontent.com/QuantumPackage/qp2/master/data/qp2.png" width="250" --->
+<img src="https://trex-coe.eu/sites/default/files/styles/responsive_no_crop/public/2021-12/Risorsa%2014_0.png"  width="250">
 
 
 [![DOI](https://zenodo.org/badge/167513335.svg)](https://zenodo.org/badge/latestdoi/167513335)
@@ -44,9 +45,19 @@ https://arxiv.org/abs/1902.08154
 
 # Credits
 
+* [TREX Center of Excellence](https://trex-coe.eu)
+* [ERC PTEROSOR](https://lcpq.github.io/PTEROSOR)
 * [CNRS](http://www.cnrs.fr)
 * [Laboratoire de Chimie et Physique Quantiques](http://lcpq.ups-tlse.fr)
 * [Laboratoire de Chimie Théorique](http://www.lct.jussieu.fr)
 * [Argonne  Leadership Computing  Facility](http://alcf.anl.gov)
 * [CALMIP](https://www.calmip.univ-toulouse.fr)
 
+
+------------------------------
+
+<img src="https://lcpq.github.io/PTEROSOR/img/ERC.png" width="200" />
+
+[TREX: Targeting Real Chemical Accuracy at the Exascale](https://trex-coe.eu) project has received funding from the European Union’s Horizon 2020 - Research and Innovation program - under grant agreement no. 952165. The content of this document does not represent the opinion of the European Union, and the European Union is not responsible for any use that might be made of such content.
+
+[PTEROSOR](https://lcpq.github.io/PTEROSOR) project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Grant agreement No. 863481).

bin/qp_convert_output_to_ezfio

@@ -146,6 +146,17 @@ def write_ezfio(res, filename):
     ezfio.set_ao_basis_ao_nucl(at)
     ezfio.set_ao_basis_ao_prim_num(num_prim)
     ezfio.set_ao_basis_ao_power(power_x + power_y + power_z)
+    try:
+      normf = res.normf
+      if normf == 0:
+        ezfio.set_ao_basis_ao_normalized(True)
+      elif normf == 1:
+        ezfio.set_ao_basis_ao_normalized(False)
+      else:
+        print("BUG in NORMF")
+        sys.exit(0)
+    except AttributeError:
+      ezfio.set_ao_basis_ao_normalized(True)
 
     # ~#~#~#~#~#~#~ #
     # P a r s i n g #

etc/qp.rc

@@ -80,8 +80,6 @@ function qp()
       if [[ -d $NAME ]] ; then
         [[ -d $EZFIO_FILE ]] && ezfio unset_file
         ezfio set_file $NAME
-      else
-        qp_create_ezfio -h | more
       fi
       unset _ARGS
       ;;

ocaml/Command_line.ml

@@ -1,3 +1,5 @@
+exception Error of string
+
 type short_opt     = char
 type long_opt      = string 
 type optional      = Mandatory | Optional 
@@ -181,15 +183,16 @@ let set_specs specs_in =
     Getopt.parse_cmdline cmd_specs (fun x -> anon_args := !anon_args @ [x]);
 
     if show_help () then
-        (help () ; exit 0);
+        help ()
+    else
+        (* Check that all mandatory arguments are set *)
+        List.filter (fun x -> x.short <> ' ' && x.opt = Mandatory) !specs
+        |> List.iter (fun x -> 
+            match get x.long with
+            | Some _ -> ()
+            | None -> raise (Error ("--"^x.long^" option is missing."))
+            )
 
-    (* Check that all mandatory arguments are set *)
-    List.filter (fun x -> x.short <> ' ' && x.opt = Mandatory) !specs
-    |> List.iter (fun x -> 
-        match get x.long with
-        | Some _ -> ()
-        | None -> failwith ("Error: --"^x.long^" option is missing.")
-        )
 ;;
 
 

ocaml/Command_line.mli

@@ -59,6 +59,8 @@ let () =
 *)
 
 
+exception Error of string
+
 type short_opt = char
 
 type long_opt = string

ocaml/Molecule.ml

@@ -101,7 +101,7 @@ let to_string_general ~f m =
   |> String.concat "\n"
 
 let to_string =
-  to_string_general ~f:(fun x -> Atom.to_string Units.Angstrom x)
+  to_string_general ~f:(fun x -> Atom.to_string ~units:Units.Angstrom x)
 
 let to_xyz =
   to_string_general ~f:Atom.to_xyz
@@ -113,7 +113,7 @@ let of_xyz_string
     s =
   let l = String_ext.split s ~on:'\n'
        |> List.filter (fun x -> x <> "")
-       |> list_map (fun x -> Atom.of_string units x)
+       |> list_map (fun x -> Atom.of_string ~units x)
   in
   let ne = ( get_charge {
         nuclei=l ;

ocaml/qp_create_ezfio.ml

@@ -677,6 +677,7 @@ let run ?o b au c d m p cart xyz_file =
 
 let () =
 
+  try (
 
   let open Command_line in
   begin
@@ -734,7 +735,7 @@ If a file with the same name as the basis set exists, this file will be read.  O
 
   let basis =
     match Command_line.get "basis" with
-    | None -> assert false
+    | None -> ""
     | Some x -> x
   in
 
@@ -773,10 +774,14 @@ If a file with the same name as the basis set exists, this file will be read.  O
 
   let xyz_filename =
     match Command_line.anon_args () with
-    | [x] -> x
+    | []  -> failwith "input file is missing"
-    | _ -> (Command_line.help () ; failwith "input file is missing")
+    | x::_ -> x
   in
 
   run ?o:output basis au charge dummy multiplicity pseudo cart xyz_filename
+  )
+  with
+  | Failure txt  -> Printf.eprintf "Fatal error: %s\n%!" txt
+  | Command_line.Error txt  -> Printf.eprintf "Command line error: %s\n%!" txt
 
 

ocaml/qp_run.ml

@@ -110,7 +110,7 @@ let run slave ?prefix exe ezfio_file =
   let task_thread =
      let thread =
       Thread.create ( fun () ->
-         TaskServer.run port_number )
+         TaskServer.run ~port:port_number )
      in
      thread ();
   in

src/cisd/cisd.irp.f

@@ -47,6 +47,37 @@ program cisd
   PROVIDE N_states
   read_wf = .False.
   SOFT_TOUCH read_wf
+
+  integer :: i,k
+
+  if(pseudo_sym)then
+   call H_apply_cisd_sym
+  else
+   call H_apply_cisd
+  endif
+  double precision :: r1, r2
+  double precision, allocatable :: U_csf(:,:)
+
+  allocate(U_csf(N_csf,N_states))
+  U_csf = 0.d0
+  U_csf(1,1) = 1.d0
+  do k=2,N_states
+    do i=1,N_csf
+        call random_number(r1)
+        call random_number(r2)
+        r1 = dsqrt(-2.d0*dlog(r1))
+        r2 = dacos(-1.d0)*2.d0*r2
+        U_csf(i,k) = r1*dcos(r2)
+    enddo
+    U_csf(k,k) = U_csf(k,k) +100.d0
+  enddo
+  do k=1,N_states
+    call normalize(U_csf(1,k),N_csf)
+  enddo
+  call convertWFfromCSFtoDET(N_states,U_csf(1,1),psi_coef(1,1))
+  deallocate(U_csf)
+  SOFT_TOUCH psi_coef
+
   call run
 end
 
@@ -56,13 +87,7 @@ subroutine run
   double precision               :: cisdq(N_states), delta_e
   double precision,external      :: diag_h_mat_elem
 
-  if(pseudo_sym)then
-   call H_apply_cisd_sym
-  else
-   call H_apply_cisd
-  endif
   psi_coef = ci_eigenvectors
-  SOFT_TOUCH psi_coef
   call save_wavefunction_truncated(save_threshold)
   call ezfio_set_cisd_energy(CI_energy)
 

src/csf/cfgCI_utils.c

@@ -68,10 +68,16 @@ void getBFIndexList(int NSOMO, int *BF1, int *IdxListBF1){
                     break;
                 }
             }
-            BFcopy[Iidx] = -1;
+            if(countN1 <= 0){
-            BFcopy[Jidx] = -1;
+                BFcopy[Iidx] = -1;
-            IdxListBF1[Jidx] = Iidx;
+                IdxListBF1[Iidx] = Iidx;
-            IdxListBF1[Iidx] = Jidx;
+            }
+            else{
+                BFcopy[Iidx] = -1;
+                BFcopy[Jidx] = -1;
+                IdxListBF1[Jidx] = Iidx;
+                IdxListBF1[Iidx] = Jidx;
+            }
         }
     }
 
@@ -328,10 +334,21 @@ void convertCSFtoDetBasis(int64_t Isomo, int MS, int rowsmax, int colsmax, doubl
                  Get Overlap
     ************************************/
     // Fill matrix
+
+    int rowsbftodetI, colsbftodetI;
+
+    /***********************************
+         Get BFtoDeterminant Matrix
+    ************************************/
+
+    printf("In convertcsftodet\n");
+    convertBFtoDetBasis(Isomo, MS, &bftodetmatrixI, &rowsbftodetI, &colsbftodetI);
+
     int rowsI = 0;
     int colsI = 0;
 
-    getOverlapMatrix(Isomo, MS, &overlapMatrixI, &rowsI, &colsI, &NSOMO);
+    //getOverlapMatrix(Isomo, MS, &overlapMatrixI, &rowsI, &colsI, &NSOMO);
+    getOverlapMatrix_withDet(bftodetmatrixI, rowsbftodetI, colsbftodetI, Isomo, MS, &overlapMatrixI, &rowsI, &colsI, &NSOMO);
 
 
     /***********************************
@@ -342,14 +359,6 @@ void convertCSFtoDetBasis(int64_t Isomo, int MS, int rowsmax, int colsmax, doubl
 
     gramSchmidt(overlapMatrixI, rowsI, colsI, orthoMatrixI);
 
-    /***********************************
-         Get BFtoDeterminant Matrix
-    ************************************/
-
-    int rowsbftodetI, colsbftodetI;
-
-    convertBFtoDetBasis(Isomo, MS, &bftodetmatrixI, &rowsbftodetI, &colsbftodetI);
-
     /***********************************
          Get Final CSF to Det Matrix
     ************************************/
@@ -1297,12 +1306,16 @@ void getbftodetfunction(Tree *dettree, int NSOMO, int MS, int *BF1, double *rowv
     double fac = 1.0;
     for(int i = 0; i < NSOMO; i++)
         donepq[i] = 0.0;
+    for(int i=0;i<npairs;++i){
+      for(int j=0;j<NSOMO;++j)
+        detslist[i*NSOMO + j]=0;
+    }
 
     for(int i = 0; i < NSOMO; i++){
         idxp = BF1[i];
         idxq = BF1[idxp];
         // Do one pair only once
-        if(donepq[idxp] > 0.0 || donepq[idxq] > 0.0) continue;
+        if(donepq[idxp] > 0.0 || donepq[idxq] > 0.0 || idxp == idxq) continue;
         fac *= 2.0;
         donepq[idxp] = 1.0;
         donepq[idxq] = 1.0;
@@ -1319,15 +1332,19 @@ void getbftodetfunction(Tree *dettree, int NSOMO, int MS, int *BF1, double *rowv
         }
         shft /= 2;
     }
-
+    
     // Now get the addresses
     int inpdet[NSOMO];
     int phase_cfg_to_qp=1;
     int addr = -1;
     for(int i = 0; i < npairs; i++){
-        for(int j = 0; j < NSOMO; j++)
+        for(int j = 0; j < NSOMO; j++) {
             inpdet[j] = detslist[i*NSOMO + j];
+            printf(" %d ",inpdet[j]);
+        }
+        printf("\n");
         findAddofDetDriver(dettree, NSOMO, inpdet, &addr);
+        printf("(%d) - addr  = %d\n",i,addr);
         // Calculate the phase for cfg to QP2 conversion
         //get_phase_cfg_to_qp_inpList(inpdet, NSOMO, &phase_cfg_to_qp);
         //rowvec[addr] = 1.0 * phaselist[i]*phase_cfg_to_qp/sqrt(fac);
@@ -1416,7 +1433,6 @@ void convertBFtoDetBasis(int64_t Isomo, int MS, double **bftodetmatrixptr, int *
         getIthBFDriver(&bftree, NSOMO, addI, BF1);
         getBFIndexList(NSOMO, BF1, IdxListBF1);
 
-
         // Get ith row
         getbftodetfunction(&dettree, NSOMO, MS, IdxListBF1, rowvec);
 
@@ -1425,6 +1441,11 @@ void convertBFtoDetBasis(int64_t Isomo, int MS, double **bftodetmatrixptr, int *
 
         for(int k=0;k<ndets;k++)
             rowvec[k]=0.0;
+
+        for(int j=0;j<NSOMO;++j){
+          BF1[j]=0;
+          IdxListBF1[j]=0;
+        }
     }
 
     // Garbage collection
@@ -1680,7 +1701,6 @@ void getApqIJMatrixDriverArrayInp(int64_t Isomo, int64_t Jsomo, int32_t orbp, in
 
     gramSchmidt(overlapMatrixJ, rowsJ, colsJ, orthoMatrixJ);
 
-
     int rowsA = 0;
     int colsA = 0;
 

src/csf/tree_utils.c

@@ -17,27 +17,31 @@ void buildTree(Tree *bftree,
     if(isomo > NSOMOMax || icpl < 0 || izeros > zeromax ) return;
 
     // If we find a valid BF assign its address
-    if(isomo == NSOMOMax){
+    if(isomo == NSOMOMax && icpl == MSmax){
         (*inode)->addr = bftree->rootNode->addr;
         bftree->rootNode->addr += 1;
         return;
     }
 
     // Call 0 branch
-    if(((*inode)->C0) == NULL && izeros+1 <= zeromax){
+    if(izeros+1 <= zeromax){
-        ((*inode)->C0) = malloc(sizeof(Node));
+      if(((*inode)->C0) == NULL){
-        (*(*inode)->C0) = (Node){ .C0 = NULL, .C1 = NULL, .PREV = *inode, .addr = -1, .cpl = 0, .iSOMO = isomo };
+          ((*inode)->C0) = malloc(sizeof(Node));
-        buildTree(bftree, &(*inode)->C0, isomo+1, izeros+1, icpl+1, NSOMOMax, MSmax);
+          (*(*inode)->C0) = (Node){ .C0 = NULL, .C1 = NULL, .PREV = *inode, .addr = -1, .cpl = 0, .iSOMO = isomo };
+          buildTree(bftree, &(*inode)->C0, isomo+1, izeros+1, icpl+1, NSOMOMax, MSmax);
+      }
+      else buildTree(bftree, &(*inode)->C0, isomo+1, izeros+1, icpl+1, NSOMOMax, MSmax);
     }
-    else buildTree(bftree, &(*inode)->C0, isomo+1, izeros+1, icpl+1, NSOMOMax, MSmax);
 
     // Call 1 branch
-    if(((*inode)->C1) == NULL && icpl-1 >= 0){
+    if(icpl-1 >=0){
-        ((*inode)->C1) = malloc(sizeof(Node));
+      if(((*inode)->C1) == NULL){
-        (*(*inode)->C1) = (Node){ .C0 = NULL, .C1 = NULL, .PREV = *inode, .addr = -1, .cpl = 1, .iSOMO = isomo };
+          ((*inode)->C1) = malloc(sizeof(Node));
-        buildTree(bftree, &(*inode)->C1, isomo+1, izeros+0, icpl-1, NSOMOMax, MSmax);
+          (*(*inode)->C1) = (Node){ .C0 = NULL, .C1 = NULL, .PREV = *inode, .addr = -1, .cpl = 1, .iSOMO = isomo };
+          buildTree(bftree, &(*inode)->C1, isomo+1, izeros+0, icpl-1, NSOMOMax, MSmax);
+      }
+      else buildTree(bftree, &(*inode)->C1, isomo+1, izeros+0, icpl-1, NSOMOMax, MSmax);
     }
-    else buildTree(bftree, &(*inode)->C1, isomo+1, izeros+0, icpl-1, NSOMOMax, MSmax);
 
     return;
 }

src/davidson/diagonalization_hcsf_dressed.irp.f

@@ -124,7 +124,7 @@ subroutine davidson_diag_csf_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,sze_csf,N
     stop -1
   endif
 
-  itermax = max(2,min(davidson_sze_max, sze/N_st_diag))+1
+  itermax = max(2,min(davidson_sze_max, sze_csf/N_st_diag))+1
   itertot = 0
 
   if (state_following) then
@@ -263,29 +263,20 @@ subroutine davidson_diag_csf_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,sze_csf,N
   ! ===================
 
   converged = .False.
-
+  call convertWFfromDETtoCSF(N_st_diag,u_in(1,1),U_csf(1,1))
   do k=N_st+1,N_st_diag
-    do i=1,sze
+    do i=1,sze_csf
         call random_number(r1)
         call random_number(r2)
         r1 = dsqrt(-2.d0*dlog(r1))
         r2 = dtwo_pi*r2
-        u_in(i,k) = r1*dcos(r2) * u_in(i,k-N_st)
+        U_csf(i,k) = r1*dcos(r2) * u_csf(i,k-N_st)
     enddo
-    u_in(k,k) = u_in(k,k) + 10.d0
+    U_csf(k,k) = u_csf(k,k) + 10.d0
   enddo
   do k=1,N_st_diag
-    call normalize(u_in(1,k),sze)
+    call normalize(U_csf(1,k),sze_csf)
   enddo
-
-  do k=1,N_st_diag
-    do i=1,sze
-      U(i,k) = u_in(i,k)
-    enddo
-  enddo
-
-  ! Make random verctors eigenstates of S2
-  call convertWFfromDETtoCSF(N_st_diag,U(1,1),U_csf(1,1))
   call convertWFfromCSFtoDET(N_st_diag,U_csf(1,1),U(1,1))
 
   do while (.not.converged)

src/davidson/diagonalize_ci.irp.f

@@ -1,3 +1,13 @@
+BEGIN_PROVIDER [ character*(3), sigma_vector_algorithm ]
+ implicit none
+ BEGIN_DOC
+ ! If 'det', use <Psi_det|H|Psi_det> in Davidson
+ !
+ ! If 'cfg', use <Psi_csf|H|Psi_csf> in Davidson
+ END_DOC
+ sigma_vector_algorithm = 'det'
+END_PROVIDER
+
 BEGIN_PROVIDER [ double precision, CI_energy, (N_states_diag) ]
   implicit none
   BEGIN_DOC
@@ -61,9 +71,18 @@ END_PROVIDER
    if (diag_algorithm == "Davidson") then
 
      if (do_csf) then
-       call davidson_diag_H_csf(psi_det,CI_eigenvectors, &
+       if (sigma_vector_algorithm == 'det') then
-         size(CI_eigenvectors,1),CI_electronic_energy,               &
+         call davidson_diag_H_csf(psi_det,CI_eigenvectors, &
-         N_det,N_csf,min(N_det,N_states),min(N_det,N_states_diag),N_int,0,converged)
+           size(CI_eigenvectors,1),CI_electronic_energy,               &
+           N_det,N_csf,min(N_det,N_states),min(N_det,N_states_diag),N_int,0,converged)
+!       else if (sigma_vector_algorithm == 'cfg') then
+!       call davidson_diag_H_csf(psi_det,CI_eigenvectors, &
+!         size(CI_eigenvectors,1),CI_electronic_energy,               &
+!         N_det,N_csf,min(N_det,N_states),min(N_det,N_states_diag),N_int,0,converged)
+!       else
+!         print *, irp_here
+!         stop 'bug'
+       endif
      else
        call davidson_diag_HS2(psi_det,CI_eigenvectors, CI_s2, &
          size(CI_eigenvectors,1),CI_electronic_energy,               &

src/determinants/determinants.irp.f

@@ -77,7 +77,7 @@ BEGIN_PROVIDER [ integer, psi_det_size ]
   END_DOC
   PROVIDE ezfio_filename
   logical                        :: exists
-  psi_det_size = 1
+  psi_det_size = N_states
   PROVIDE mpi_master
   if (read_wf) then
     if (mpi_master) then
@@ -85,7 +85,7 @@ BEGIN_PROVIDER [ integer, psi_det_size ]
       if (exists) then
         call ezfio_get_determinants_n_det(psi_det_size)
       else
-        psi_det_size = 1
+        psi_det_size = N_states
       endif
       call write_int(6,psi_det_size,'Dimension of the psi arrays')
     endif

src/mo_two_e_ints/mo_bi_integrals.irp.f

@@ -130,7 +130,6 @@ subroutine four_idx_dgemm
   real(integral_kind), allocatable :: buffer_value(:)
   size_buffer = min(ao_num*ao_num*ao_num,16000000)
 
-  print *, 'Storing'
   !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(i,j,k,l,buffer_value,buffer_i,n_integrals)
   allocate ( buffer_i(size_buffer), buffer_value(size_buffer) )
 
@@ -164,7 +163,6 @@ subroutine four_idx_dgemm
 
   deallocate (a1)
 
-  print *, 'Unique'
   call map_unique(mo_integrals_map)
 
   integer*8                      :: get_mo_map_size, mo_map_size

src/tools/truncate_wf.irp.f

@@ -54,11 +54,23 @@ subroutine routine_s2
  double precision, allocatable :: psi_coef_tmp(:,:)
  integer :: i,j,k
  double precision :: accu(N_states)
+ integer :: weights(0:16), ix
+ double precision :: x
 
- print *,  'Weights of the CFG'
+ weights(:) = 0
  do i=1,N_det
-   print *,  i, real(weight_configuration(det_to_configuration(i),:)), real(sum(weight_configuration(det_to_configuration(i),:)))
+   x = -dlog(1.d-32+sum(weight_configuration(det_to_configuration(i),:)))/dlog(10.d0)
+   ix = min(int(x), 16)
+   weights(ix) += 1
  enddo
+
+ print *,  'Histogram of the weights of the CFG'
+ do i=0,15
+   print *, '  10^{-', i, '}   ', weights(i) 
+ end do
+ print *, '< 10^{-', 15, '}   ', weights(16) 
+
+
  print*, 'Min weight of the configuration?'
  read(5,*) wmin