Merge pull request #36 from QuantumPackage/dev-stable

Dev stable

ocaml/Atom.ml

@@ -22,10 +22,15 @@ let of_string ~units s =
     }
   | [ name; x; y; z ] ->
     let e = Element.of_string name in
+    begin
+      try
       { element = e ;
         charge  = Element.to_charge e;
         coord   = Point3d.of_string ~units (String.concat " " [x; y; z])
       }
+      with
+      | err -> (Printf.eprintf "name = \"%s\"\nxyz = (%s,%s,%s)\n%!" name x y z ; raise err)
+    end
   | _ -> raise (AtomError s)
 
 

ocaml/Molecule.ml

@@ -142,13 +142,21 @@ let of_xyz_string
   result
 
 
+let regexp_r = Str.regexp {|
|}
+let regexp_t = Str.regexp {|	|}
+
 
 let of_xyz_file
     ?(charge=(Charge.of_int 0)) ?(multiplicity=(Multiplicity.of_int 1))
     ?(units=Units.Angstrom)
     filename =
   let lines =
-    match Io_ext.input_lines filename with
+    Io_ext.input_lines filename
+    |> List.map (fun s -> Str.global_replace regexp_r "" s)
+    |> List.map (fun s -> Str.global_replace regexp_t " " s)
+  in
+  let lines =
+    match lines with
     | natoms :: title :: rest ->
           let natoms = 
             try
@@ -173,6 +181,8 @@ let of_zmt_file
     ?(units=Units.Angstrom)
     filename =
   Io_ext.read_all filename
+  |> Str.global_replace regexp_r ""
+  |> Str.global_replace regexp_t " "
   |> Zmatrix.of_string
   |> Zmatrix.to_xyz_string
   |> of_xyz_string ~charge ~multiplicity ~units

ocaml/Point3d.ml

@@ -24,7 +24,9 @@ let of_string ~units s =
   let l = s
           |> String_ext.split ~on:' '
           |> List.filter (fun x -> x <> "")
-          |> list_map float_of_string
+          |> list_map (fun x ->
+                 try float_of_string x with
+                 | Failure msg -> (Printf.eprintf "Bad string: \"%s\"\n%!" x ; failwith msg) )
           |> Array.of_list
   in
   { x = l.(0) *. f ;

plugins/local/non_h_ints_mu/NEED

@@ -3,3 +3,4 @@ hamiltonian
 jastrow
 ao_tc_eff_map
 bi_ortho_mos
+trexio

plugins/local/non_h_ints_mu/deb_aos.irp.f

@@ -35,32 +35,59 @@ subroutine print_aos()
 
   PROVIDE final_grid_points aos_in_r_array aos_grad_in_r_array aos_lapl_in_r_array
 
-  write(1000, *) n_points_final_grid
   do ipoint = 1, n_points_final_grid
     r(:) = final_grid_points(:,ipoint)
     write(1000, '(3(f15.7, 3X))') r
   enddo
   
+double precision :: accu_vgl(5)
+double precision :: accu_vgl_nrm(5)
+
   do ipoint = 1, n_points_final_grid
-    r(:) = final_grid_points(:,ipoint)
     do i = 1, ao_num
       ao_val    = aos_in_r_array     (i,ipoint)
       ao_der(:) = aos_grad_in_r_array(i,ipoint,:)
       ao_lap    = aos_lapl_in_r_array(1,i,ipoint) + aos_lapl_in_r_array(2,i,ipoint) + aos_lapl_in_r_array(3,i,ipoint)
-      write(1010, '(5(f15.7, 3X))') ao_val, ao_der, ao_lap
+      write(111, '(5(f15.7, 3X))') ao_val, ao_der, ao_lap
     enddo
   enddo
  
   do ipoint = 1, n_points_final_grid
-    r(:) = final_grid_points(:,ipoint)
+    do i = 1, ao_num
-    do i = 1, mo_num
+      ao_val    = aos_in_r_array_qmckl     (i,ipoint)
-      mo_val    = mos_in_r_array     (i,ipoint)
+      ao_der(:) = aos_grad_in_r_array_qmckl(i,ipoint,:)
-      mo_der(:) = mos_grad_in_r_array(i,ipoint,:)
+      ao_lap    = aos_lapl_in_r_array_qmckl(i,ipoint)
-      mo_lap    = mos_lapl_in_r_array(i,ipoint,1) + mos_lapl_in_r_array(i,ipoint,2) + mos_lapl_in_r_array(i,ipoint,3)
+      write(222, '(5(f15.7, 3X))') ao_val, ao_der, ao_lap
-      write(2010, '(5(f15.7, 3X))') mo_val, mo_der, mo_lap
     enddo
   enddo
  
+  accu_vgl = 0.d0
+  accu_vgl_nrm = 0.d0
+  do ipoint = 1, n_points_final_grid
+    do i = 1, ao_num
+      ao_val    = aos_in_r_array     (i,ipoint)
+      ao_der(:) = aos_grad_in_r_array(i,ipoint,:)
+      ao_lap    = aos_lapl_in_r_array(1,i,ipoint) + aos_lapl_in_r_array(2,i,ipoint) + aos_lapl_in_r_array(3,i,ipoint)
+      accu_vgl_nrm(1) += dabs(ao_val)
+      accu_vgl_nrm(2) += dabs(ao_der(1))
+      accu_vgl_nrm(3) += dabs(ao_der(2))
+      accu_vgl_nrm(4) += dabs(ao_der(3))
+      accu_vgl_nrm(5) += dabs(ao_lap)
+
+      ao_val    -= aos_in_r_array_qmckl     (i,ipoint)
+      ao_der(:) -= aos_grad_in_r_array_qmckl(i,ipoint,:)
+      ao_lap    -= aos_lapl_in_r_array_qmckl(i,ipoint)
+      accu_vgl(1) += dabs(ao_val)
+      accu_vgl(2) += dabs(ao_der(1))
+      accu_vgl(3) += dabs(ao_der(2))
+      accu_vgl(4) += dabs(ao_der(3))
+      accu_vgl(5) += dabs(ao_lap)
+    enddo
+
+  enddo
+  accu_vgl(:) *= 1.d0 / accu_vgl_nrm(:)
+  print *, accu_vgl
+
   return
 end
 

plugins/local/non_h_ints_mu/qmckl.irp.f

@@ -75,3 +75,107 @@ BEGIN_PROVIDER [ integer*8, qmckl_ctx_jastrow ]
   endif
 
 END_PROVIDER
+
+
+ BEGIN_PROVIDER [ double precision, aos_in_r_array_qmckl, (ao_num,n_points_final_grid)]
+&BEGIN_PROVIDER [ double precision, aos_grad_in_r_array_qmckl, (ao_num,n_points_final_grid,3)]
+&BEGIN_PROVIDER [ double precision, aos_lapl_in_r_array_qmckl, (ao_num, n_points_final_grid)]
+ implicit none
+ BEGIN_DOC
+ ! AOS computed with qmckl
+ END_DOC
+ use qmckl
+
+ integer*8 :: qmckl_ctx
+ integer(qmckl_exit_code) :: rc
+
+ qmckl_ctx = qmckl_context_create()
+
+ rc = qmckl_trexio_read(qmckl_ctx, trexio_file, 1_8*len(trim(trexio_filename)))
+ if (rc /= QMCKL_SUCCESS) then
+   print *, irp_here, 'qmckl error in read_trexio'
+   rc = qmckl_check(qmckl_ctx, rc)
+   stop -1
+ endif
+
+ rc = qmckl_set_point(qmckl_ctx, 'N', n_points_final_grid*1_8, final_grid_points, n_points_final_grid*3_8)
+ if (rc /= QMCKL_SUCCESS) then
+   print *, irp_here, 'qmckl error in set_electron_point'
+   rc = qmckl_check(qmckl_ctx, rc)
+   stop -1
+ endif
+
+ double precision, allocatable :: vgl(:,:,:)
+ allocate( vgl(ao_num,5,n_points_final_grid))
+ rc = qmckl_get_ao_basis_ao_vgl_inplace(qmckl_ctx, vgl, n_points_final_grid*ao_num*5_8)
+ if (rc /= QMCKL_SUCCESS) then
+   print *, irp_here, 'qmckl error in get_ao_vgl'
+   rc = qmckl_check(qmckl_ctx, rc)
+   stop -1
+ endif
+
+ integer :: i,k
+ do k=1,n_points_final_grid
+   do i=1,ao_num
+     aos_in_r_array_qmckl(i,k) = vgl(i,1,k)
+     aos_grad_in_r_array_qmckl(i,k,1) = vgl(i,2,k)
+     aos_grad_in_r_array_qmckl(i,k,2) = vgl(i,3,k)
+     aos_grad_in_r_array_qmckl(i,k,3) = vgl(i,4,k)
+     aos_lapl_in_r_array_qmckl(i,k) = vgl(i,5,k)
+   enddo
+ enddo
+
+END_PROVIDER
+
+
+ BEGIN_PROVIDER [ double precision, mos_in_r_array_qmckl, (mo_num,n_points_final_grid)]
+&BEGIN_PROVIDER [ double precision, mos_grad_in_r_array_qmckl, (mo_num,n_points_final_grid,3)]
+&BEGIN_PROVIDER [ double precision, mos_lapl_in_r_array_qmckl, (mo_num, n_points_final_grid)]
+ implicit none
+ BEGIN_DOC
+ ! moS computed with qmckl
+ END_DOC
+ use qmckl
+
+ integer*8 :: qmckl_ctx
+ integer(qmckl_exit_code) :: rc
+
+ qmckl_ctx = qmckl_context_create()
+
+ rc = qmckl_trexio_read(qmckl_ctx, trexio_file, 1_8*len(trim(trexio_filename)))
+ if (rc /= QMCKL_SUCCESS) then
+   print *, irp_here, 'qmckl error in read_trexio'
+   rc = qmckl_check(qmckl_ctx, rc)
+   stop -1
+ endif
+
+ rc = qmckl_set_point(qmckl_ctx, 'N', n_points_final_grid*1_8, final_grid_points, n_points_final_grid*3_8)
+ if (rc /= QMCKL_SUCCESS) then
+   print *, irp_here, 'qmckl error in set_electron_point'
+   rc = qmckl_check(qmckl_ctx, rc)
+   stop -1
+ endif
+
+ double precision, allocatable :: vgl(:,:,:)
+ allocate( vgl(mo_num,5,n_points_final_grid))
+ rc = qmckl_get_mo_basis_mo_vgl_inplace(qmckl_ctx, vgl, n_points_final_grid*mo_num*5_8)
+ if (rc /= QMCKL_SUCCESS) then
+   print *, irp_here, 'qmckl error in get_mo_vgl'
+   rc = qmckl_check(qmckl_ctx, rc)
+   stop -1
+ endif
+
+ integer :: i,k
+ do k=1,n_points_final_grid
+   do i=1,mo_num
+     mos_in_r_array_qmckl(i,k) = vgl(i,1,k)
+     mos_grad_in_r_array_qmckl(i,k,1) = vgl(i,2,k)
+     mos_grad_in_r_array_qmckl(i,k,2) = vgl(i,3,k)
+     mos_grad_in_r_array_qmckl(i,k,3) = vgl(i,4,k)
+     mos_lapl_in_r_array_qmckl(i,k) = vgl(i,5,k)
+   enddo
+ enddo
+
+END_PROVIDER
+
+

plugins/local/tc_keywords/EZFIO.cfg

@@ -100,7 +100,7 @@ doc: If |true|, the states are re-ordered to match the input states
 default: False
 interface: ezfio,provider,ocaml
 
-[symetric_fock_tc]
+[symmetric_fock_tc]
 type: logical
 doc: If |true|, using F+F^t as Fock TC
 interface: ezfio,provider,ocaml

scripts/get_fci_tc_conv.sh

@@ -1,2 +1,2 @@
 file=$1
-grep "Ndet,E,E+PT2,E+RPT2,|PT2|=" $file | cut -d "=" -f 2 > ${file}.conv_fci_tc
+grep "Ndet,E,E+PT2,pt2_minus,pt2_plus,pt2_abs=" $file | cut -d "=" -f 2 > ${file}.conv_fci_tc

src/mo_optimization/gradient_list_opt.irp.f

@@ -319,7 +319,7 @@ call omp_set_max_active_levels(4)
 ! \end{equation}
 
 ! We need a vector to use the gradient. Here the gradient is a
-! antisymetric matrix so we can transform it in a vector of length
+! antisymmetric matrix so we can transform it in a vector of length
 ! mo_num*(mo_num-1)/2.
 
 ! Here we do these two things at the same time.

src/mo_optimization/gradient_opt.irp.f

@@ -284,7 +284,7 @@ call omp_set_max_active_levels(4)
 ! \end{equation}
 
 ! We need a vector to use the gradient. Here the gradient is a
-! antisymetric matrix so we can transform it in a vector of length
+! antisymmetric matrix so we can transform it in a vector of length
 ! mo_num*(mo_num-1)/2.
 
 ! Here we do these two things at the same time.

src/scf_utils/roothaan_hall_scf.irp.f

@@ -217,7 +217,7 @@ END_DOC
   do while (i<mo_num)
     j=i
     m=1
-    do while ( (j<mo_num).and.(fock_matrix_diag_mo(j+1)-fock_matrix_diag_mo(i) < 1.d-8) )
+    do while ( (j+1<mo_num).and.(fock_matrix_diag_mo(j+1)-fock_matrix_diag_mo(i) < 1.d-8) )
       j += 1
       m += 1
     enddo

src/utils/util.irp.f

@@ -327,12 +327,12 @@ subroutine wall_time(t)
 end
 
 BEGIN_PROVIDER [ integer, nproc ]
+  use omp_lib
   implicit none
   BEGIN_DOC
   ! Number of current OpenMP threads
   END_DOC
 
-  integer, external              :: omp_get_num_threads
   nproc = 1
   !$OMP PARALLEL
   !$OMP MASTER