diff --git a/ocaml/Atom.ml b/ocaml/Atom.ml index d02b20d8..49e788e8 100644 --- a/ocaml/Atom.ml +++ b/ocaml/Atom.ml @@ -22,10 +22,15 @@ let of_string ~units s = } | [ name; x; y; z ] -> let e = Element.of_string name in - { element = e ; - charge = Element.to_charge e; - coord = Point3d.of_string ~units (String.concat " " [x; y; z]) - } + 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) diff --git a/ocaml/Molecule.ml b/ocaml/Molecule.ml index 603244c8..3771b6f9 100644 --- a/ocaml/Molecule.ml +++ b/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 diff --git a/ocaml/Point3d.ml b/ocaml/Point3d.ml index 57b02bfe..4df375e4 100644 --- a/ocaml/Point3d.ml +++ b/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 ; diff --git a/plugins/local/non_h_ints_mu/NEED b/plugins/local/non_h_ints_mu/NEED index 48c1c24b..5ca1d543 100644 --- a/plugins/local/non_h_ints_mu/NEED +++ b/plugins/local/non_h_ints_mu/NEED @@ -3,3 +3,4 @@ hamiltonian jastrow ao_tc_eff_map bi_ortho_mos +trexio diff --git a/plugins/local/non_h_ints_mu/deb_aos.irp.f b/plugins/local/non_h_ints_mu/deb_aos.irp.f index a84e1b91..4012f47c 100644 --- a/plugins/local/non_h_ints_mu/deb_aos.irp.f +++ b/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 - enddo - enddo - - do ipoint = 1, n_points_final_grid - r(:) = final_grid_points(:,ipoint) - do i = 1, mo_num - mo_val = mos_in_r_array (i,ipoint) - mo_der(:) = mos_grad_in_r_array(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(2010, '(5(f15.7, 3X))') mo_val, mo_der, mo_lap + write(111, '(5(f15.7, 3X))') ao_val, ao_der, ao_lap enddo enddo + do ipoint = 1, n_points_final_grid + do i = 1, ao_num + 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) + write(222, '(5(f15.7, 3X))') ao_val, ao_der, ao_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 diff --git a/plugins/local/non_h_ints_mu/qmckl.irp.f b/plugins/local/non_h_ints_mu/qmckl.irp.f index 1df80457..4d419e24 100644 --- a/plugins/local/non_h_ints_mu/qmckl.irp.f +++ b/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 + + diff --git a/plugins/local/tc_keywords/EZFIO.cfg b/plugins/local/tc_keywords/EZFIO.cfg index 33b9db57..1e89eaa4 100644 --- a/plugins/local/tc_keywords/EZFIO.cfg +++ b/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 diff --git a/scripts/get_fci_tc_conv.sh b/scripts/get_fci_tc_conv.sh index 643f3ac0..f0c99baf 100755 --- a/scripts/get_fci_tc_conv.sh +++ b/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 diff --git a/src/mo_optimization/gradient_list_opt.irp.f b/src/mo_optimization/gradient_list_opt.irp.f index 9b7228c7..9331c80f 100644 --- a/src/mo_optimization/gradient_list_opt.irp.f +++ b/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. diff --git a/src/mo_optimization/gradient_opt.irp.f b/src/mo_optimization/gradient_opt.irp.f index 25be6b5a..10d42b35 100644 --- a/src/mo_optimization/gradient_opt.irp.f +++ b/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. diff --git a/src/scf_utils/roothaan_hall_scf.irp.f b/src/scf_utils/roothaan_hall_scf.irp.f index 3f5c8549..e0fe5319 100644 --- a/src/scf_utils/roothaan_hall_scf.irp.f +++ b/src/scf_utils/roothaan_hall_scf.irp.f @@ -217,7 +217,7 @@ END_DOC do while (i