Merge branch 'master' of github.com:LCPQ/quantum_package

Conflicts:
	src/Dets/determinants.irp.f
	src/NEEDED_MODULES

ocaml/Input_determinants.ml

@@ -508,11 +508,18 @@ psi_det                = %s
     let psi_coef = 
       let rec read_coefs accu = function
       | [] -> List.rev accu
+      | ""::""::tail -> read_coefs accu tail
       | ""::c::tail -> 
+          let c =
+            Float.of_string c
+            |> Det_coef.of_float 
+          in
           read_coefs (c::accu) tail
       | _::tail -> read_coefs accu tail
       in
-      let a = read_coefs [] dets
+      let a =
+        read_coefs [] dets
+        |> List.map ~f:(fun x -> Det_coef.to_string x)
         |> String.concat ~sep:" "
       in
       "(psi_coef ("^a^"))"

ocaml/qp_edit.ml

@@ -2,16 +2,13 @@ open Qputils;;
 open Qptypes;;
 open Core.Std;;
 
-let file_header filename = Printf.sprintf
+(** Interactive editing of the input.
-"
-==================================================================
-                       Quantum Package
-==================================================================
 
-Editing file `%s`
+@author A. Scemama
+*)
 
-" filename
 
+(** Keywords used to define input sections *)
 type keyword = 
 | Ao_basis
 | Bielec_integrals
@@ -24,6 +21,7 @@ type keyword =
 | Nuclei
 ;;
 
+
 let keyword_to_string = function
 | Ao_basis         -> "AO basis"
 | Bielec_integrals -> "Two electron integrals"
@@ -36,12 +34,30 @@ let keyword_to_string = function
 | Nuclei           -> "Molecule"
 ;;
 
+
+
+(** Create the header of the temporary file *)
+let file_header filename = 
+  Printf.sprintf  "
+==================================================================
+                       Quantum Package
+==================================================================
+
+Editing file `%s`
+
+" filename
+;;
+  
+
+(** Creates the header of a section *)
 let make_header kw =
   let s = keyword_to_string kw in
   let l = String.length s in
   "\n\n"^s^"\n"^(String.init l ~f:(fun _ -> '='))^"\n\n"
 ;;
 
+
+(** Returns the rst string of section [s] *)
 let get s = 
   let header = (make_header s) in
   let f (read,to_rst) = 
@@ -79,6 +95,8 @@ let get s =
   rst
 ;;
 
+
+(** Applies the changes from the string [str] corresponding to section [s] *)
 let set str s = 
   let header = (make_header s) in
   match String.substr_index ~pos:0 ~pattern:header str with
@@ -103,7 +121,7 @@ let set str s =
         | None -> ()
         with
         | _ -> (Printf.eprintf "Info: Read error in %s\n%!"
-               (keyword_to_string s))
+               (keyword_to_string s); ignore (of_rst str) )
       in
       let open Input in
         match s with
@@ -120,6 +138,7 @@ let set str s =
 ;;
 
 
+(** Creates the temporary file for interactive editing *)
 let create_temp_file ezfio_filename fields =
   let temp_filename  = Filename.temp_file "qp_edit_" ".rst" in
   begin
@@ -132,6 +151,8 @@ let create_temp_file ezfio_filename fields =
   ; temp_filename
 ;;
 
+
+
 let run ezfio_filename =
 
   (* Open EZFIO *)
@@ -173,8 +194,8 @@ let run ezfio_filename =
     | Some editor -> editor
     | None -> "vi"
   in
-  let command = Printf.sprintf "%s %s" editor temp_filename in
+  Printf.sprintf "%s %s" editor temp_filename 
-  Sys.command_exn command;
+  |> Sys.command_exn ;
 
   (* Re-read the temp file *)
   let temp_string  =
@@ -188,6 +209,24 @@ let run ezfio_filename =
 ;;
 
 
+(** Create a backup file in case of an exception *)
+let create_backup ezfio_filename =
+  Printf.sprintf "
+    rm -f %s/backup.tgz ;
+    tar -zcf .backup.tgz %s && mv .backup.tgz %s/backup.tgz
+  "
+    ezfio_filename ezfio_filename ezfio_filename
+  |> Sys.command_exn
+;;
+
+
+(** Restore the backup file when an exception occuprs *)
+let restore_backup ezfio_filename =
+  Printf.sprintf "tar -zxf %s/backup.tgz"
+    ezfio_filename 
+  |> Sys.command_exn
+;;
+
 
 let spec =
   let open Command.Spec in
@@ -216,11 +255,34 @@ Edit input data
        with
        | _ msg -> print_string ("\n\nError\n\n"^msg^"\n\n")
     *)
-    (fun ezfio_file () -> run ezfio_file)
+    (fun ezfio_file () ->
+       try
+          run ezfio_file ;
+          (* create_backup ezfio_file; *)
+       with
+       | Failure exc 
+       | Invalid_argument exc as e -> 
+         begin
+           Printf.eprintf "=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-\n\n";
+           Printf.eprintf "%s\n\n" exc;
+           Printf.eprintf "=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-\n\n";
+           (* restore_backup ezfio_file; *) 
+           raise e
+         end
+       | Assert_failure (file, line, ch) as e ->
+         begin
+           Printf.eprintf "=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-\n\n";
+           Printf.eprintf "Assert error in file $QPACKAGE_ROOT/ocaml/%s, line %d, character %d\n\n" file line ch;
+           Printf.eprintf "=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-\n\n";
+           (* restore_backup ezfio_file; *)
+           raise e
+         end
+    )
 ;;
 
 let () =
-    Command.run command
+  Command.run command;
+  exit 0
 ;;
 
 

scripts/install_ocaml.sh

@@ -6,6 +6,11 @@
 QPACKAGE_ROOT=${PWD}
 PACKAGES="core cryptokit"
 
+function asksure() {
+  echo -n "Are you sure (Y/N)? "
+  return $retval
+}
+
 if [[ -f quantum_package.rc ]]
 then
   source quantum_package.rc
@@ -14,7 +19,18 @@ make -C ocaml Qptypes.ml &> /dev/null
 if [[ $? -ne 0 ]]
 then
 
-  rm -rf -- ${HOME}/ocamlbrew
+  if [[ -d ${HOME}/ocamlbrew ]]
+  then
+    echo "Remove directory ${HOME}/ocamlbrew? [Y/n]"
+    while read -r -n 1 -s answer; do
+      if [[ $answer = [YyNn] ]]; then
+         [[ $answer = [Yy] ]] && rm -rf -- ${HOME}/ocamlbrew
+         [[ $answer = [Nn] ]] && exit 1
+        break
+      fi
+    done
+    
+  fi
   scripts/fetch_from_web.py "https://raw.github.com/hcarty/ocamlbrew/master/ocamlbrew-install" ocamlbrew-install.sh 
   cat < ocamlbrew-install.sh | env OCAMLBREW_FLAGS="-r" bash | tee ocamlbrew_install.log
   grep "source " ocamlbrew_install.log | grep "etc/ocamlbrew.bashrc"  >> quantum_package.rc

src/Dets/determinants.irp.f

@@ -213,7 +213,6 @@ END_PROVIDER
 END_PROVIDER 
 
 
-
 BEGIN_PROVIDER [ double precision, psi_coef, (psi_det_size,N_states_diag) ]
   implicit none
   BEGIN_DOC

src/Dets/diagonalize_CI.irp.f

@@ -61,22 +61,35 @@ END_PROVIDER
     call lapack_diag(eigenvalues,eigenvectors,                       &
         H_matrix_all_dets,size(H_matrix_all_dets,1),N_det)
     CI_electronic_energy(:) = 0.d0
+    do i=1,N_det
+       CI_eigenvectors(i,1) = eigenvectors(i,1)
+    enddo
     integer :: i_state
     double precision :: s2
-    j=0
     i_state = 0
-    do while(i_state.lt.min(N_states_diag,N_det))
+    do j=1,N_det
-      j+=1
       call get_s2_u0(psi_det,eigenvectors(1,j),N_det,N_det,s2)
-      if(dabs(s2-expected_s2).le.0.1d0)then
+      print *, 'j = ',j,s2, expected_s2
+      if(dabs(s2-expected_s2).le.0.3d0)then
        i_state += 1
+       print *,  'i_state = ',i_state
        do i=1,N_det
          CI_eigenvectors(i,i_state) = eigenvectors(i,j)
        enddo
        CI_electronic_energy(i_state) = eigenvalues(j)
        CI_eigenvectors_s2(i_state) = s2
       endif
+      if (i_state.ge.N_states_diag) then
+        exit
+      endif
     enddo
+!    if(i_state < min(N_states_diag,N_det))then
+!     print *, 'pb with the number of states'
+!     print *, 'i_state = ',i_state
+!     print *, 'N_states_diag ',N_states_diag
+!     print *,'stopping ...'
+!     stop
+!    endif
     deallocate(eigenvectors,eigenvalues)
   endif
   

src/Utils/one_e_integration.irp.f

@@ -31,45 +31,6 @@ double precision function overlap_gaussian_x(A_center,B_center,alpha,beta,power_
   overlap_gaussian_x*= fact_p
 end
 
-subroutine test(alpha,beta,gama,a,b,A_center,B_center,Nucl_center,overlap_x,overlap_y,overlap_z,overlap)
-   implicit none
-  include 'constants.F'
- integer, intent(in)            :: a(3),b(3)         ! powers : (x-xa)**a_x = (x-A(1))**a(1)
- double precision, intent(in)   :: alpha, beta, gama ! exponents
- double precision, intent(in)   :: A_center(3)       ! A center
- double precision, intent(in)   :: B_center (3)      ! B center
- double precision, intent(in)   :: Nucl_center(3)    ! B center
- double precision, intent(out)  :: overlap_x,overlap_y,overlap_z,overlap
- integer :: i,j
- double precision :: dx,Lx,nx,x(3)
- nx = 100000000
- Lx = 25.d0
- dx = dble(Lx/nx)
- overlap_x = 0.d0
- overlap_y = 0.d0
- overlap_z = 0.d0
- x(1) = -12.5d0
- x(2) = -12.5d0
- x(3) = -12.5d0
- do i = 1,nx 
-  overlap_x += (x(1) - A_center(1))**a(1) * (x(1) - B_center(1))**b(1) &
-              * dexp(-alpha*(x(1) - A_center(1))**2) * dexp(-beta*(x(1) - B_center(1))**2) * dexp(-gama*(x(1) - Nucl_center(1))**2) 
-
-  overlap_y += (x(2) - A_center(2))**a(2) * (x(2) - B_center(2))**b(2) &
-              * dexp(-alpha*(x(2) - A_center(2))**2) * dexp(-beta*(x(2) - B_center(2))**2) * dexp(-gama*(x(2) - Nucl_center(2))**2) 
-  overlap_z += (x(3) - A_center(3))**a(3) * (x(3) - B_center(3))**b(3) &
-              * dexp(-alpha*(x(3) - A_center(3))**2) * dexp(-beta*(x(3) - B_center(3))**2) * dexp(-gama*(x(3) - Nucl_center(3))**2) 
-  x(1) += dx
-  x(2) += dx
-  x(3) += dx
- enddo
- overlap_x = overlap_x * dx
- overlap_y = overlap_y * dx
- overlap_z = overlap_z * dx
- overlap = overlap_x * overlap_y * overlap_z
-
-end
-
 
 subroutine overlap_A_B_C(dim,alpha,beta,gama,a,b,A_center,B_center,Nucl_center,overlap)
  implicit none