diff --git a/docs/source/users_guide/qpsh.rst b/docs/source/users_guide/qpsh.rst index 7a1a874c..049ab8e3 100644 --- a/docs/source/users_guide/qpsh.rst +++ b/docs/source/users_guide/qpsh.rst @@ -8,7 +8,7 @@ qpsh :command:`qpsh` is the |qp| shell. It is a Bash shell with all the -required evironment variables loaded, a modified prompt, and the +required environment variables loaded, a modified prompt, and the :ref:`qp` command. diff --git a/ocaml/Input_ao_basis.ml b/ocaml/Input_ao_basis.ml index 506cf069..d9e28e04 100644 --- a/ocaml/Input_ao_basis.ml +++ b/ocaml/Input_ao_basis.ml @@ -58,52 +58,70 @@ end = struct ;; let read_ao_prim_num () = - Ezfio.get_ao_basis_ao_prim_num () - |> Ezfio.flattened_ezfio - |> Array.map AO_prim_number.of_int + if Ezfio.has_ao_basis_ao_prim_num () then + Ezfio.get_ao_basis_ao_prim_num () + |> Ezfio.flattened_ezfio + |> Array.map AO_prim_number.of_int + else + [||] ;; let read_ao_prim_num_max () = - Ezfio.get_ao_basis_ao_prim_num () - |> Ezfio.flattened_ezfio - |> Array.fold_left (fun x y -> if x>y then x else y) 0 - |> AO_prim_number.of_int + if Ezfio.has_ao_basis_ao_prim_num () then + Ezfio.get_ao_basis_ao_prim_num () + |> Ezfio.flattened_ezfio + |> Array.fold_left (fun x y -> if x>y then x else y) 0 + |> AO_prim_number.of_int + else + AO_prim_number.of_int 0 ;; let read_ao_nucl () = - let nmax = Nucl_number.get_max () in - Ezfio.get_ao_basis_ao_nucl () - |> Ezfio.flattened_ezfio - |> Array.map (fun x-> Nucl_number.of_int ~max:nmax x) + if Ezfio.has_ao_basis_ao_nucl () then + let nmax = Nucl_number.get_max () in + Ezfio.get_ao_basis_ao_nucl () + |> Ezfio.flattened_ezfio + |> Array.map (fun x-> Nucl_number.of_int ~max:nmax x) + else + [||] ;; let read_ao_power () = - let x = Ezfio.get_ao_basis_ao_power () in - let dim = x.Ezfio.dim.(0) in - let data = Ezfio.flattened_ezfio x in - let result = Array.init dim (fun x -> "") in - for i=1 to dim - do - if (data.(i-1) > 0) then - result.(i-1) <- result.(i-1)^"x"^(string_of_int data.(i-1)); - if (data.(dim+i-1) > 0) then - result.(i-1) <- result.(i-1)^"y"^(string_of_int data.(dim+i-1)); - if (data.(2*dim+i-1) > 0) then - result.(i-1) <- result.(i-1)^"z"^(string_of_int data.(2*dim+i-1)); - done; - Array.map Angmom.Xyz.of_string result + if Ezfio.has_ao_basis_ao_power () then + let x = Ezfio.get_ao_basis_ao_power () in + let dim = x.Ezfio.dim.(0) in + let data = Ezfio.flattened_ezfio x in + let result = Array.init dim (fun x -> "") in + for i=1 to dim + do + if (data.(i-1) > 0) then + result.(i-1) <- result.(i-1)^"x"^(string_of_int data.(i-1)); + if (data.(dim+i-1) > 0) then + result.(i-1) <- result.(i-1)^"y"^(string_of_int data.(dim+i-1)); + if (data.(2*dim+i-1) > 0) then + result.(i-1) <- result.(i-1)^"z"^(string_of_int data.(2*dim+i-1)); + done; + Array.map Angmom.Xyz.of_string result + else + [||] ;; let read_ao_coef () = - Ezfio.get_ao_basis_ao_coef () - |> Ezfio.flattened_ezfio - |> Array.map AO_coef.of_float + if Ezfio.has_ao_basis_ao_coef () then + Ezfio.get_ao_basis_ao_coef () + |> Ezfio.flattened_ezfio + |> Array.map AO_coef.of_float + else + [||] ;; let read_ao_expo () = - Ezfio.get_ao_basis_ao_expo () - |> Ezfio.flattened_ezfio - |> Array.map AO_expo.of_float + if Ezfio.has_ao_basis_ao_expo () then + Ezfio.get_ao_basis_ao_expo () + |> Ezfio.flattened_ezfio + |> Array.map AO_expo.of_float + else + [||] ;; let read_ao_cartesian () = diff --git a/ocaml/qp_run.ml b/ocaml/qp_run.ml index 0cb862ae..bb886143 100644 --- a/ocaml/qp_run.ml +++ b/ocaml/qp_run.ml @@ -132,6 +132,7 @@ let run slave ?prefix exe ezfio_file = (** Run executable *) let prefix = match prefix with + | Some "gdb" -> "gdb --args " | Some x -> x^" " | None -> "" and exe = diff --git a/plugins/local/cipsi_tc_bi_ortho/stochastic_cipsi.irp.f b/plugins/local/cipsi_tc_bi_ortho/stochastic_cipsi.irp.f index 6b8f3b42..0c61e38f 100644 --- a/plugins/local/cipsi_tc_bi_ortho/stochastic_cipsi.irp.f +++ b/plugins/local/cipsi_tc_bi_ortho/stochastic_cipsi.irp.f @@ -155,6 +155,7 @@ subroutine run_stochastic_cipsi call pt2_alloc(pt2_data_err, N_states) call ZMQ_pt2(E_tc, pt2_data, pt2_data_err, relative_error,0) ! Stochastic PT2 and selection call diagonalize_CI_tc_bi_ortho(ndet, E_tc,norm) + call print_summary_tc(psi_energy_with_nucl_rep, pt2_data, pt2_data_err, N_det, N_configuration, N_states, psi_s2) call pt2_dealloc(pt2_data) call pt2_dealloc(pt2_data_err) diff --git a/plugins/local/jastrow/EZFIO.cfg b/plugins/local/jastrow/EZFIO.cfg index 8fd2d05a..9d4cf431 100644 --- a/plugins/local/jastrow/EZFIO.cfg +++ b/plugins/local/jastrow/EZFIO.cfg @@ -1,7 +1,25 @@ +[log_jpsi] +type: logical +doc: If |true|, the Jpsi is taken as log(1+psi_cor) +interface: ezfio,provider,ocaml +default: False + + +[mu_of_r_tc] +type: character*(32) +doc: type of the mu(r): [ Standard | Erfmu | Erfmugauss ] +interface: ezfio,provider,ocaml +default: Standard + +[mu_of_r_av] +type: logical +doc: If |true|, take the second formula for mu(r) +interface: ezfio,provider,ocaml +default: False [j2e_type] type: character*(32) -doc: type of the 2e-Jastrow: [ None | Mu | Mu_Nu | Mur | Boys | Boys_Handy | Qmckl ] +doc: type of the 2e-Jastrow: [ None | Mu | Mugauss | Mu_Nu | Mur | Murgauss | Bump | Boys | Boys_Handy | Qmckl ] interface: ezfio,provider,ocaml default: Mu diff --git a/plugins/local/non_h_ints_mu/NEED b/plugins/local/non_h_ints_mu/NEED index 5ca1d543..bfc4f311 100644 --- a/plugins/local/non_h_ints_mu/NEED +++ b/plugins/local/non_h_ints_mu/NEED @@ -4,3 +4,4 @@ jastrow ao_tc_eff_map bi_ortho_mos trexio +mu_of_r diff --git a/plugins/local/non_h_ints_mu/deb_deriv_mu.irp.f b/plugins/local/non_h_ints_mu/deb_deriv_mu.irp.f new file mode 100644 index 00000000..1c8d7198 --- /dev/null +++ b/plugins/local/non_h_ints_mu/deb_deriv_mu.irp.f @@ -0,0 +1,28 @@ +program test_j_mu_of_r + implicit none + double precision :: x,mu_min,dmu,mu_max, mu, mu_p, mu_m + double precision :: j_simple,j_p, j_m,numeric_d_mu,d_dx_mu + double precision :: accu + integer :: npt,i + npt = 1000 + mu_min = 0.3d0 + mu_max = 10.d0 + dmu = (mu_max - mu_min)/dble(npt) + x = 0.7d0 + mu = mu_min + do i = 1, npt + call get_deriv_mu_j12(x,mu,d_dx_mu) + mu_p = mu + dmu + mu_m = mu - dmu + j_p = j_simple(x,mu_p) + j_m = j_simple(x,mu_m) + numeric_d_mu = 0.5d0 * (j_p - j_m)/dmu + print*,mu + print*,numeric_d_mu,d_dx_mu,dabs(d_dx_mu-numeric_d_mu) + accu += dabs(d_dx_mu-numeric_d_mu) + mu += dmu + enddo + accu *= dmu + print*,'accu = ',accu +end + diff --git a/plugins/local/non_h_ints_mu/deb_j_bump.irp.f b/plugins/local/non_h_ints_mu/deb_j_bump.irp.f new file mode 100644 index 00000000..82872357 --- /dev/null +++ b/plugins/local/non_h_ints_mu/deb_j_bump.irp.f @@ -0,0 +1,98 @@ +program test_j_mu_of_r + implicit none +! call routine_test_mu_of_r + call routine_test_mu_of_r_tot +end + +subroutine routine_test_mu_of_r_tot + implicit none + integer :: ipoint,k + double precision :: r2(3), weight, dr, r1(3), r1bis(3) + double precision :: accu_grad(3) + double precision :: jast,grad_jast_mu_r1(3),j_bump + double precision :: jast_p,jast_m,num_grad_jast_mu_r1(3) + + dr = 0.00001d0 + r2 = 0.d0 + r2(1) = 0.5d0 + r2(2) = -0.1d0 + r2(3) = 1.0d0 + accu_grad = 0.d0 + do ipoint = 1, n_points_final_grid + r1(1:3) = final_grid_points(1:3,ipoint) + weight = final_weight_at_r_vector(ipoint) +! call grad_j_sum_mu_of_r(r1,r2,jast,grad_jast_mu_r1) + call get_grad_j_bump_mu_of_r(r1,r2,grad_jast_mu_r1) + double precision :: norm,error + norm = 0.D0 + do k = 1, 3 + r1bis= r1 + r1bis(k) += dr + jast_p = j_bump(r1bis,r2,a_boys) + + r1bis= r1 + r1bis(k) -= dr + jast_m = j_bump(r1bis,r2,a_boys) + + num_grad_jast_mu_r1(k) = (jast_p - jast_m)/(2.d0* dr) + norm += num_grad_jast_mu_r1(k)*num_grad_jast_mu_r1(k) + enddo + error = 0.d0 + do k = 1, 3 + error += dabs(grad_jast_mu_r1(k) - num_grad_jast_mu_r1(k)) + enddo + error *= 0.33333333d0 + norm = dsqrt(norm) + if(norm.gt.1.d-05)then + if(dabs(error/norm).gt.dr)then + print*,'/////' + print*,error,norm + print*,grad_jast_mu_r1 + print*,num_grad_jast_mu_r1 + endif + endif + do k = 1,3 + accu_grad(k) += weight * dabs(grad_jast_mu_r1(k) - num_grad_jast_mu_r1(k)) + enddo + enddo + print*,'accu_grad = ' + print*, accu_grad + +end + +subroutine routine_test_mu_of_r + implicit none + integer :: ipoint,k + double precision :: weight, dr, r1(3), r1bis(3),accu_grad(3),num_grad_mu_r1(3) + double precision :: mu_r1,dm_r1, mu_der_r1(3), grad_dm_r1(3) + double precision :: mu_der_rp(3), grad_dm_rp(3),mu_rp + double precision :: mu_der_rm(3), grad_dm_rm(3),mu_rm + + dr = 0.0001d0 + accu_grad = 0.d0 + do ipoint = 1, n_points_final_grid + r1(1:3) = final_grid_points(1:3,ipoint) + weight = final_weight_at_r_vector(ipoint) + call grad_mu_of_r_mean_field(r1,mu_r1,dm_r1, mu_der_r1, grad_dm_r1) + do k = 1, 3 + r1bis= r1 + r1bis(k) += dr + call grad_mu_of_r_mean_field(r1bis,mu_rp, dm_r1, mu_der_rp, grad_dm_r1) + + r1bis= r1 + r1bis(k) -= dr + call grad_mu_of_r_mean_field(r1bis,mu_rm, dm_r1, mu_der_rm, grad_dm_r1) + num_grad_mu_r1(k) = (mu_rp - mu_rm)/(2.d0* dr) +! print*,jast_mu_r1_p,jast_mu_r1_m + enddo + print*,'/////' + print*,mu_der_r1 + print*,num_grad_mu_r1 + do k = 1,3 + accu_grad(k) += weight * dabs(mu_der_r1(k) - num_grad_mu_r1(k)) + enddo + enddo + print*,'accu_grad = ' + print*, accu_grad + +end diff --git a/plugins/local/non_h_ints_mu/deb_j_gauss.irp.f b/plugins/local/non_h_ints_mu/deb_j_gauss.irp.f new file mode 100644 index 00000000..264a6f04 --- /dev/null +++ b/plugins/local/non_h_ints_mu/deb_j_gauss.irp.f @@ -0,0 +1,62 @@ +program test_j_mu_of_r + implicit none +! call routine_test_mu_of_r + call routine_test_mu_of_r_tot +end + +subroutine routine_test_mu_of_r_tot + implicit none + integer :: ipoint,k + double precision :: r2(3), weight, dr, r1(3), r1bis(3) + double precision :: accu_grad(3) + double precision :: jast,grad_jast(3),j_bump,j12_mu + double precision :: jast_p,jast_m,num_grad_jast(3) + + dr = 0.00001d0 + r2 = 0.d0 + r2(1) = 0.5d0 + r2(2) = -0.1d0 + r2(3) = 1.0d0 + accu_grad = 0.d0 + do ipoint = 1, n_points_final_grid + r1(1:3) = final_grid_points(1:3,ipoint) + weight = final_weight_at_r_vector(ipoint) + call grad1_j12_mu(r1, r2, grad_jast) + grad_jast = - grad_jast + double precision :: norm,error + norm = 0.D0 + do k = 1, 3 + r1bis= r1 + r1bis(k) += dr + jast_p = j12_mu(r1bis, r2) + + r1bis= r1 + r1bis(k) -= dr + jast_m = j12_mu(r1bis, r2) + + num_grad_jast(k) = (jast_p - jast_m)/(2.d0* dr) + norm += num_grad_jast(k)*num_grad_jast(k) + enddo + error = 0.d0 + do k = 1, 3 + error += dabs(grad_jast(k) - num_grad_jast(k)) + enddo + error *= 0.33333333d0 + norm = dsqrt(norm) + if(norm.gt.1.d-05)then + if(dabs(error/norm).gt.dr)then + print*,'/////' + print*,error,norm + print*,grad_jast + print*,num_grad_jast + endif + endif + do k = 1,3 + accu_grad(k) += weight * dabs(grad_jast(k) - num_grad_jast(k)) + enddo + enddo + print*,'accu_grad = ' + print*, accu_grad + +end + diff --git a/plugins/local/non_h_ints_mu/deb_j_mu_of_r.irp.f b/plugins/local/non_h_ints_mu/deb_j_mu_of_r.irp.f new file mode 100644 index 00000000..fc29bf50 --- /dev/null +++ b/plugins/local/non_h_ints_mu/deb_j_mu_of_r.irp.f @@ -0,0 +1,97 @@ +program test_j_mu_of_r + implicit none +! call routine_test_mu_of_r + call routine_test_mu_of_r_tot +end + +subroutine routine_test_mu_of_r_tot + implicit none + integer :: ipoint,k + double precision :: r2(3), weight, dr, r1(3), r1bis(3) + double precision :: accu_grad(3) + double precision :: jast,grad_jast_mu_r1(3) + double precision :: jast_p,jast_m,num_grad_jast_mu_r1(3) + + dr = 0.000001d0 + r2 = 0.d0 + r2(1) = 0.5d0 + r2(2) = -0.1d0 + r2(3) = 1.0d0 + accu_grad = 0.d0 + do ipoint = 1, n_points_final_grid + r1(1:3) = final_grid_points(1:3,ipoint) + weight = final_weight_at_r_vector(ipoint) + call grad_j_sum_mu_of_r(r1,r2,jast,grad_jast_mu_r1) + double precision :: norm,error + norm = 0.D0 + do k = 1, 3 + r1bis= r1 + r1bis(k) += dr + call get_j_sum_mu_of_r(r1bis,r2,jast_p) + + r1bis= r1 + r1bis(k) -= dr + call get_j_sum_mu_of_r(r1bis,r2,jast_m) + + num_grad_jast_mu_r1(k) = (jast_p - jast_m)/(2.d0* dr) + norm += num_grad_jast_mu_r1(k)*num_grad_jast_mu_r1(k) + enddo + error = 0.d0 + do k = 1, 3 + error += dabs(grad_jast_mu_r1(k) - num_grad_jast_mu_r1(k)) + enddo + error *= 0.33333333d0 + norm = dsqrt(norm) + if(norm.gt.1.d-05)then + if(dabs(error/norm).gt.10.d0*dr)then + print*,'/////' + print*,error,norm,dabs(error/norm) + print*,grad_jast_mu_r1 + print*,num_grad_jast_mu_r1 + endif + endif + do k = 1,3 + accu_grad(k) += weight * dabs(grad_jast_mu_r1(k) - num_grad_jast_mu_r1(k)) + enddo + enddo + print*,'accu_grad = ' + print*, accu_grad + +end + +subroutine routine_test_mu_of_r + implicit none + integer :: ipoint,k + double precision :: weight, dr, r1(3), r1bis(3),accu_grad(3),num_grad_mu_r1(3) + double precision :: mu_r1,dm_r1, mu_der_r1(3), grad_dm_r1(3) + double precision :: mu_der_rp(3), grad_dm_rp(3),mu_rp + double precision :: mu_der_rm(3), grad_dm_rm(3),mu_rm + + dr = 0.0001d0 + accu_grad = 0.d0 + do ipoint = 1, n_points_final_grid + r1(1:3) = final_grid_points(1:3,ipoint) + weight = final_weight_at_r_vector(ipoint) + call grad_mu_of_r_mean_field(r1,mu_r1,dm_r1, mu_der_r1, grad_dm_r1) + do k = 1, 3 + r1bis= r1 + r1bis(k) += dr + call grad_mu_of_r_mean_field(r1bis,mu_rp, dm_r1, mu_der_rp, grad_dm_r1) + + r1bis= r1 + r1bis(k) -= dr + call grad_mu_of_r_mean_field(r1bis,mu_rm, dm_r1, mu_der_rm, grad_dm_r1) + num_grad_mu_r1(k) = (mu_rp - mu_rm)/(2.d0* dr) +! print*,jast_mu_r1_p,jast_mu_r1_m + enddo + print*,'/////' + print*,mu_der_r1 + print*,num_grad_mu_r1 + do k = 1,3 + accu_grad(k) += weight * dabs(mu_der_r1(k) - num_grad_mu_r1(k)) + enddo + enddo + print*,'accu_grad = ' + print*, accu_grad + +end diff --git a/plugins/local/non_h_ints_mu/deb_j_psi.irp.f b/plugins/local/non_h_ints_mu/deb_j_psi.irp.f new file mode 100644 index 00000000..1b034684 --- /dev/null +++ b/plugins/local/non_h_ints_mu/deb_j_psi.irp.f @@ -0,0 +1,131 @@ +program test_j_mu_of_r + implicit none + call routine_deb_j_psi +! call routine_deb_denom +end + +subroutine routine_deb_j_psi + implicit none + integer :: ipoint,k + double precision :: r2(3), weight, dr, r1(3), r1bis(3) + double precision :: accu_grad(3) + double precision :: jast,grad_jast(3),j_bump,jastrow_psi,grad_jast_bis(3) + double precision :: jast_p,jast_m,num_grad_jast(3) + + dr = 0.00001d0 + r2 = 0.d0 + r2(1) = 0.5d0 + r2(2) = -0.1d0 + r2(3) = 1.0d0 + accu_grad = 0.d0 + do ipoint = 1, n_points_final_grid + r1(1:3) = final_grid_points(1:3,ipoint) + weight = final_weight_at_r_vector(ipoint) + call get_grad_r1_jastrow_psi(r1,r2,grad_jast,jast) +! grad_jast = - grad_jast + double precision :: norm,error + norm = 0.D0 + do k = 1, 3 + r1bis= r1 + r1bis(k) += dr + call get_grad_r1_jastrow_psi(r1bis,r2,grad_jast_bis,jast_p) + + r1bis= r1 + r1bis(k) -= dr + call get_grad_r1_jastrow_psi(r1bis,r2,grad_jast_bis,jast_m) + + num_grad_jast(k) = (jast_p - jast_m)/(2.d0* dr) + norm += num_grad_jast(k)*num_grad_jast(k) + enddo + error = 0.d0 + do k = 1, 3 + error += dabs(grad_jast(k) - num_grad_jast(k)) + enddo + error *= 0.33333333d0 + norm = dsqrt(norm) + if(norm.gt.1.d-05)then + if(dabs(error/norm).gt.dr)then + print*,'/////' + print*,error,norm + print*,grad_jast + print*,num_grad_jast + endif + endif + do k = 1,3 + accu_grad(k) += weight * dabs(grad_jast(k) - num_grad_jast(k)) + enddo + enddo + print*,'accu_grad = ' + print*, accu_grad + +end + + +subroutine routine_deb_denom + implicit none + integer :: ipoint,k,i,j + double precision :: r2(3), weight, dr, r1(3), r1bis(3) + double precision :: accu_grad(3) + double precision :: jast,grad_jast(3),j_bump,jastrow_psi,grad_jast_bis(3) + double precision :: jast_p,jast_m,num_grad_jast(3) + + dr = 0.00001d0 + r2 = 0.d0 + r2(1) = 0.5d0 + r2(2) = -0.1d0 + r2(3) = 1.0d0 + double precision, allocatable :: mos_array_r1(:), mos_array_r2(:) + double precision, allocatable :: mos_grad_array_r1(:,:),mos_grad_array_r2(:,:) + allocate(mos_array_r1(mo_num), mos_array_r2(mo_num)) + allocate(mos_grad_array_r1(3,mo_num), mos_grad_array_r2(3,mo_num)) + do i = 1, 1 + do j = 1, 1 + accu_grad = 0.d0 + call give_all_mos_and_grad_at_r(r2,mos_array_r2,mos_grad_array_r2) + do ipoint = 1, n_points_final_grid + r1(1:3) = final_grid_points(1:3,ipoint) + weight = final_weight_at_r_vector(ipoint) + call give_all_mos_and_grad_at_r(r1,mos_array_r1,mos_grad_array_r1) + call denom_jpsi(i,j,a_boys, mos_array_r1,mos_grad_array_r1,mos_array_r2,jast, grad_jast) + double precision :: norm,error + norm = 0.D0 + do k = 1, 3 + r1bis= r1 + r1bis(k) += dr + call give_all_mos_and_grad_at_r(r1bis,mos_array_r1,mos_grad_array_r1) + call denom_jpsi(i,j,a_boys, mos_array_r1,mos_grad_array_r1,mos_array_r2,jast_p, grad_jast_bis) + + r1bis= r1 + r1bis(k) -= dr + call give_all_mos_and_grad_at_r(r1bis,mos_array_r1,mos_grad_array_r1) + call denom_jpsi(i,j,a_boys, mos_array_r1,mos_grad_array_r1,mos_array_r2,jast_m, grad_jast_bis) + + num_grad_jast(k) = (jast_p - jast_m)/(2.d0* dr) + norm += num_grad_jast(k)*num_grad_jast(k) + enddo + error = 0.d0 + do k = 1, 3 + error += dabs(grad_jast(k) - num_grad_jast(k)) + enddo + error *= 0.33333333d0 + norm = dsqrt(norm) + if(norm.gt.1.d-05)then + if(dabs(error/norm).gt.dr)then + print*,'/////' + print*,error,norm + print*,grad_jast + print*,num_grad_jast + endif + endif + do k = 1,3 + accu_grad(k) += weight * dabs(grad_jast(k) - num_grad_jast(k)) + enddo + enddo + print*,'i,j = ',i,j + print*,'accu_grad = ' + print*, accu_grad + enddo + enddo + +end + diff --git a/plugins/local/non_h_ints_mu/j_bump.irp.f b/plugins/local/non_h_ints_mu/j_bump.irp.f new file mode 100644 index 00000000..1731bb72 --- /dev/null +++ b/plugins/local/non_h_ints_mu/j_bump.irp.f @@ -0,0 +1,90 @@ +double precision function wigner_radius(rho) + implicit none + include 'constants.include.F' + double precision, intent(in) :: rho + wigner_radius = 4.d0 * pi * rho * 0.333333333333d0 + wigner_radius = wigner_radius**(-0.3333333d0) +end + +double precision function j_bump(r1,r2,a) + implicit none + include 'constants.include.F' + double precision, intent(in) :: r1(3),r2(3),a + double precision :: inv_a,factor,x_scaled,scalar + double precision :: r12 + r12 = (r1(1) - r2(1))*(r1(1) - r2(1)) + r12 += (r1(2) - r2(2))*(r1(2) - r2(2)) + r12 += (r1(3) - r2(3))*(r1(3) - r2(3)) + r12 = dsqrt(r12) + inv_a = 1.d0/a + x_scaled = r12*inv_a*inv_sq_pi + x_scaled*= x_scaled + j_bump = 0.5d0 * (r12-a) * dexp(-x_scaled) +end + +subroutine get_grad_j_bump(x,a,grad) + implicit none + BEGIN_DOC + ! gradient of the Jastrow with a bump + ! + ! j(x,a) = 1/2 * (x-a)* exp[-(x/(a*sqrt(pi)))^2] + ! + ! d/dx j(x,a) = 1/(2 pi a^2) * exp[-(x/(a*sqrt(pi)))^2] * (pi a^2 + 2 a x - 2x^2) + END_DOC + include 'constants.include.F' + double precision, intent(in) :: x,a + double precision, intent(out) :: grad + double precision :: inv_a,factor,x_scaled,scalar + inv_a = 1.d0/a + factor = 0.5d0*inv_pi*inv_a*inv_a + x_scaled = x*inv_a*inv_sq_pi + x_scaled*= x_scaled + grad = factor * dexp(-x_scaled) * (pi*a*a + 2.d0 * a*x - 2.d0*x*x) +end + +subroutine get_d_da_j_bump(x,a,d_da) + implicit none + BEGIN_DOC + ! Derivative with respect by to the parameter "a" of the Jastrow with a bump + ! + ! j(x,a) = 1/2 * (x-a)* exp[-(x/(a*sqrt(pi)))^2] + ! + ! d/da j(x,a) = - 1/(pi*a^3) * exp[-(x/(a*sqrt(pi)))^2] * (-2 x^3 + 2 a x^2 + pi a^x3) + END_DOC + include 'constants.include.F' + double precision, intent(in) :: x,a + double precision, intent(out) :: d_da + double precision :: factor, inv_a,x_scaled,scalar + inv_a = 1.d0/a + factor = inv_a*inv_a*inv_a*inv_pi + x_scaled = x*inv_a*inv_sq_pi + x_scaled*= x_scaled + d_da = factor * dexp(-x_scaled) * (-2.d0 * x*x*x + 2.d0*x*x*a+pi*a*a*a) +end + +subroutine get_grad_j_bump_mu_of_r(r1,r2,grad_j_bump) + implicit none + BEGIN_DOC + ! d/dx1 j(x,a(r1,r2)) where j(x,a) is the Jastrow with a bump + ! + ! j(x,a) = 1/2 * (x-a)* exp[-(x/(a*sqrt(pi)))^2] + ! + ! a(r1,r2) = [rho(r1) a(r1) + rho(r2) a(r2)]/[rho(r1) + rho(r2)] + ! + ! d/dx1 j(x,a) = d/dx1 j(x,a(r1,r2)) + END_DOC + double precision, intent(in) :: r1(3),r2(3) + double precision, intent(out):: grad_j_bump(3) + double precision :: r12,r12_vec(3),grad_scal,inv_r12 + r12_vec = r1 - r2 + r12 = (r1(1) - r2(1))*(r1(1) - r2(1)) + r12 += (r1(2) - r2(2))*(r1(2) - r2(2)) + r12 += (r1(3) - r2(3))*(r1(3) - r2(3)) + r12 = dsqrt(r12) + call get_grad_j_bump(r12,a_boys,grad_scal) + if(r12.lt.1.d-10)then + grad_j_bump = 0.d0 + else + grad_j_bump = grad_scal * r12_vec/r12 + endif +end diff --git a/plugins/local/non_h_ints_mu/jast_deriv.irp.f b/plugins/local/non_h_ints_mu/jast_deriv.irp.f index 9a430135..1f97c18a 100644 --- a/plugins/local/non_h_ints_mu/jast_deriv.irp.f +++ b/plugins/local/non_h_ints_mu/jast_deriv.irp.f @@ -31,7 +31,7 @@ grad1_u12_squared_num = 0.d0 if( ((j2e_type .eq. "Mu") .and. (env_type .eq. "None")) .or. & - (j2e_type .eq. "Mur") ) then + (j2e_type .eq. "Mur").or.(j2e_type .eq. "Mugauss") .or. (j2e_type .eq. "Murgauss")) then !$OMP PARALLEL & !$OMP DEFAULT (NONE) & diff --git a/plugins/local/non_h_ints_mu/jast_deriv_mu_of_r.irp.f b/plugins/local/non_h_ints_mu/jast_deriv_mu_of_r.irp.f new file mode 100644 index 00000000..8f5aee0c --- /dev/null +++ b/plugins/local/non_h_ints_mu/jast_deriv_mu_of_r.irp.f @@ -0,0 +1,306 @@ +subroutine get_j_sum_mu_of_r(r1,r2,jast) + implicit none + double precision, intent(in) :: r1(3),r2(3) + double precision, intent(out):: jast + double precision :: mu_r1, dm_r1, grad_mu_r1(3), grad_dm_r1(3), j_mu_r1 + double precision :: mu_r2, dm_r2, grad_mu_r2(3), grad_dm_r2(3), j_mu_r2 + double precision :: j12_mu_input,mu_tot,r12,j_simple + jast = 0.d0 + if(murho_type==0)then +! J(r1,r2) = [rho(r1) * j(mu(r1),r12) + rho(r2) * j(mu(r2),r12)] / [rho(r1) + rho(r2)] + call grad_mu_of_r_mean_field(r1,mu_r1, dm_r1, grad_mu_r1, grad_dm_r1) + call grad_mu_of_r_mean_field(r2,mu_r2, dm_r2, grad_mu_r2, grad_dm_r2) + j_mu_r1 = j12_mu_input(r1, r2, mu_r1) + j_mu_r2 = j12_mu_input(r1, r2, mu_r2) + if(dm_r1 + dm_r2.lt.1.d-7)return + jast = (dm_r1 * j_mu_r1 + dm_r2 * j_mu_r2) / (dm_r1 + dm_r2) + else if(murho_type==1)then +! J(r1,r2) = j(0.5 * (mu(r1)+mu(r2)),r12), MU(r1,r2) = 0.5 *(mu(r1)+mu(r2)) + call grad_mu_of_r_mean_field(r1,mu_r1, dm_r1, grad_mu_r1, grad_dm_r1) + call grad_mu_of_r_mean_field(r2,mu_r2, dm_r2, grad_mu_r2, grad_dm_r2) + mu_tot = 0.5d0 * (mu_r1 + mu_r2) + jast = j12_mu_input(r1, r2, mu_tot) + else if(murho_type==2)then +! MU(r1,r2) = (rho(1) * mu(r1)+ rho(2) * mu(r2))/(rho(1)+rho(2)) +! J(r1,r2) = j(MU(r1,r2),r12) + call grad_mu_of_r_mean_field(r1,mu_r1, dm_r1, grad_mu_r1, grad_dm_r1) + call grad_mu_of_r_mean_field(r2,mu_r2, dm_r2, grad_mu_r2, grad_dm_r2) + double precision :: mu_tmp, dm_tot, dm_tot_inv + dm_tot = dm_r1**a_boys + dm_r2**a_boys ! rho(1)**alpha+rho(2)**alpha + if(dm_tot.lt.1.d-12)then + dm_tot_inv = 1.d+12 + else + dm_tot_inv = 1.d0/dm_tot + endif + mu_tmp = dm_r1**a_boys * mu_r1 + dm_r2**a_boys * mu_r2 !rho(1)**alpha * mu(r1)+ rho(2)**alpha * mu(r2) + mu_tot = nu_erf * mu_tmp*dm_tot_inv ! + r12 = (r1(1) - r2(1)) * (r1(1) - r2(1)) + r12 += (r1(2) - r2(2)) * (r1(2) - r2(2)) + r12 += (r1(3) - r2(3)) * (r1(3) - r2(3)) + r12 = dsqrt(r12) + jast = j_simple(r12,mu_tot) + endif + +end + +subroutine grad_j_sum_mu_of_r(r1,r2,jast,grad_jast) + implicit none + include 'constants.include.F' + BEGIN_DOC + END_DOC + double precision, intent(in) :: r1(3),r2(3) + double precision, intent(out):: jast, grad_jast(3) + jast = 0.d0 + grad_jast = 0.d0 + double precision :: num, denom, grad_num(3), grad_denom(3) + double precision :: j_r1, grad_j_r1(3),j_r2, grad_j_r2(3) + double precision :: dm_r1, grad_dm_r1(3), grad_jmu_r2(3) + double precision :: dm_r2, grad_dm_r2(3),mu_r2, grad_mu_r2(3),mu_r1 + double precision :: j12_mu_input,r12,grad_mu_r1(3),grad_jmu_r1(3) + double precision :: mu_tot,dx,dy,dz,r12_vec(3),d_dmu_j,d_dr12_j + + dx = r1(1) - r2(1) + dy = r1(2) - r2(2) + dz = r1(3) - r2(3) + + r12 = dsqrt(dx * dx + dy * dy + dz * dz) + if(r12.gt.1.d-10)then + r12_vec(1) = dx + r12_vec(2) = dy + r12_vec(3) = dz + r12_vec *= 1.d0/r12 + ! r12_vec = grad_r1 (r12) + else + r12 = 1.d-10 + r12_vec = 0.d0 + endif + + if(murho_type==0)then +! J(r1,r2) = [rho(r1) * j(mu(r1),r12) + rho(r2) * j(mu(r2),r12)] / [rho(r1) + rho(r2)] +! +! = num(r1,r2) / denom(r1,r2) +! +! d/dx1 J(r1,r2) = [denom(r1,r2) X d/dx1 num(r1,r2) - num(r1,r2) X d/dx1 denom(r1,r2) ] / denom(r1,r2)^2 +! +! d/dx1 num(r1,r2) = j(mu(r1),r12)*d/dx1 rho(r1) + rho(r1) * d/dx1 j(mu(r1),r12) +! + rho(r2) d/dx1 j(mu(r2),r12) +! d/dx1 denom(r1,r2) = d/dx1 rho(r1) + call grad_j_mu_of_r_1(r1,r2,j_r1, grad_j_r1,dm_r1, grad_dm_r1) + call grad_mu_of_r_mean_field(r2,mu_r2, dm_r2, grad_mu_r2, grad_dm_r2) + j_r2 = j12_mu_input(r1, r2, mu_r2) ! j(mu(r2),r1,r2) + num = dm_r1 * j_r1 + dm_r2 * j_r2 + denom = dm_r1 + dm_r2 + if(denom.lt.1.d-7)return + jast = num / denom + + grad_denom = grad_dm_r1 + call grad_j12_mu_input(r1, r2, mu_r2, grad_jmu_r2,r12) + grad_num = j_r1 * grad_dm_r1 + dm_r1 * grad_j_r1 + dm_r2 * grad_jmu_r2 + grad_jast = (grad_num * denom - num * grad_denom)/(denom*denom) + else if(murho_type==1)then +! J(r1,r2) = j(0.5 * (mu(r1)+mu(r2)),r12), MU(r1,r2) = 0.5 *(mu(r1)+mu(r2)) +! +! d/dx1 J(r1,r2) = d/dx1 j(MU(r1,r2),r12)|MU=cst +! + d/dMU [j(MU,r12)] +! x d/d(mu(r1)) MU(r1,r2) +! x d/dx1 mu(r1) +! = 0.5 * (1 - erf(MU(r1,r2) *r12))/r12 * (x1 - x2) == grad_jmu_r1 +! + e^{-(r12*MU(r1,r2))^2}/(2 sqrt(pi) * MU(r1,r2)^2) +! x 0.5 +! x d/dx1 mu(r1) + call grad_mu_of_r_mean_field(r1,mu_r1, dm_r1, grad_mu_r1, grad_dm_r1) + call grad_mu_of_r_mean_field(r2,mu_r2, dm_r2, grad_mu_r2, grad_dm_r2) + mu_tot = 0.5d0 * (mu_r1 + mu_r2) + call grad_j12_mu_input(r1, r2, mu_tot, grad_jmu_r1,r12) + grad_jast = grad_jmu_r1 + grad_jast+= dexp(-r12*mu_tot*r12*mu_tot) * inv_sq_pi_2 /(mu_tot* mu_tot) * 0.5d0 * grad_mu_r1 + else if(murho_type==2)then +! MU(r1,r2) = beta * (rho(1)**alpha * mu(r1)+ rho(2)**alpha * mu(r2))/(rho(1)**alpha+rho(2)**alpha) +! J(r1,r2) = j(MU(r1,r2),r12) +! +! d/dx1 J(r1,r2) = d/dx1 j(MU(r1,r2),r12)|MU=cst +! + d/dMU [j(MU,r12)] +! x d/d(mu(r1)) MU(r1,r2) +! x d/dx1 mu(r1) +! = 0.5 * (1 - erf(MU(r1,r2) *r12))/r12 * (x1 - x2) == grad_jmu_r1 +! + 0.5 e^{-(r12*MU(r1,r2))^2}/(2 sqrt(pi) * MU(r1,r2)^2) +! x d/dx1 MU(r1,r2) +! with d/dx1 MU(r1,r2) = beta * {[mu(1) d/dx1 [rho(1)**alpha] + rho(1)**alpha * d/dx1 mu(1)](rho(1)**alpha+rho(2)**alpha) +! - MU(1,2) d/dx1 [rho(1)]**alpha}/(rho(1)**alpha+rho(2)**alpha)^2 +! d/dx1 [rho(1)]**alpha = alpha [rho(1)]**(alpha-1) d/dx1 rho(1) +! + call grad_mu_of_r_mean_field(r1,mu_r1, dm_r1, grad_mu_r1, grad_dm_r1) + call grad_mu_of_r_mean_field(r2,mu_r2, dm_r2, grad_mu_r2, grad_dm_r2) + double precision :: dm_tot,dm_tot_inv,grad_mu_tot(3),mu_tmp,grad_dm_r1_alpha(3),d_dx_j + dm_tot = dm_r1**a_boys + dm_r2**a_boys ! rho(1)**alpha+rho(2)**alpha + grad_dm_r1_alpha = a_boys * dm_r1**(a_boys-1) * grad_dm_r1 + if(dm_tot.lt.1.d-12)then + dm_tot_inv = 1.d+12 + else + dm_tot_inv = 1.d0/dm_tot + endif + mu_tmp = dm_r1**a_boys * mu_r1 + dm_r2**a_boys * mu_r2 !rho(1)**alpha * mu(r1)+ rho(2)**alpha * mu(r2) + mu_tot = nu_erf * mu_tmp*dm_tot_inv ! + grad_mu_tot = ( mu_r1 * grad_dm_r1_alpha + dm_r1**a_boys * grad_mu_r1 ) * dm_tot & + - mu_tmp * grad_dm_r1_alpha + grad_mu_tot *= dm_tot_inv * dm_tot_inv * nu_erf + call get_deriv_r12_j12(r12,mu_tot,d_dr12_j) ! d/dr12 j(MU(r1,r2,r12) + ! d/dx1 j(MU(r1,r2),r12) | MU(r1,r2) = cst + ! d/dr12 j(MU(r1,r2,r12) x d/dx1 r12 + grad_jmu_r1 = d_dr12_j * r12_vec +! call grad_j12_mu_input(r1, r2, mu_tot, grad_jmu_r1,r12) + grad_jast = grad_jmu_r1 + ! d/dMU j(MU(r1,r2),r12) + call get_deriv_mu_j12(r12,mu_tot,d_dmu_j) + grad_jast+= d_dmu_j * grad_mu_tot + else if(murho_type==-1)then +! J(r1,r2) = 0.5 * [j(mu(r1),r12) + j(mu(r2),r12)] +! +! d/dx1 J(r1,r2) = 0.5 * (d/dx1 j(mu(r1),r12) + d/dx1 j(mu(r2),r12)) + call grad_j_mu_of_r_1(r1,r2,j_r1, grad_j_r1,dm_r1, grad_dm_r1) + call grad_mu_of_r_mean_field(r2,mu_r2, dm_r2, grad_mu_r2, grad_dm_r2) + j_r2 = j12_mu_input(r1, r2, mu_r2) ! j(mu(r2),r1,r2) + call grad_j12_mu_input(r1, r2, mu_r2, grad_jmu_r2,r12) + jast = 0.5d0 * (j_r1 + j_r2) + grad_jast = 0.5d0 * (grad_j_r1 + grad_jmu_r2) + + endif + +end + +subroutine grad_j_mu_of_r_1(r1,r2,jast, grad_jast, dm_r1, grad_dm_r1) + implicit none + include 'constants.include.F' + BEGIN_DOC +! grad_r1 of j(mu(r1),r12) + ! + ! + ! d/dx1 j(mu(r1),r12) = exp(-(mu(r1)*r12)**2) /(2 *sqrt(pi) * mu(r1)**2 ) d/dx1 mu(r1) + ! + d/dx1 j(mu(r1),r12) + ! + ! with + ! + ! j(mu,r12) = 1/2 r12 (1 - erf(mu r12)) - 1/2 (sqrt(pi) * mu) e^{-(mu*r12)^2} + ! + ! and d/dx1 j(mu,r12) = 0.5 * (1 - erf(mu *r12))/r12 * (x1 - x2) + ! + ! d/d mu j(mu,r12) = e^{-(r12*mu)^2}/(2 sqrt(pi) * mu^2) + ! + ! here mu(r1) is obtained by MU MEAN FIELD + END_DOC + double precision, intent(in) :: r1(3),r2(3) + double precision, intent(out):: jast, grad_jast(3),dm_r1, grad_dm_r1(3) + double precision :: dx, dy, dz, r12, mu_der(3) + double precision :: mu_tmp, tmp, grad(3), mu_val + jast = 0.d0 + grad = 0.d0 + + dx = r1(1) - r2(1) + dy = r1(2) - r2(2) + dz = r1(3) - r2(3) + r12 = dsqrt(dx * dx + dy * dy + dz * dz) + ! get mu(r1) == mu_val and its gradient d/dx1 mu(r1) == mu_der + call grad_mu_of_r_mean_field(r1,mu_val, dm_r1, mu_der, grad_dm_r1) + mu_tmp = mu_val * r12 + ! evalulation of the jastrow j(mu(r1),r12) + jast = 0.5d0 * r12 * (1.d0 - derf(mu_tmp)) - inv_sq_pi_2 * dexp(-mu_tmp*mu_tmp) / mu_val + + ! tmp = exp(-(mu(r1)*r12)**2) /(2 *sqrt(pi) * mu(r1)**2 ) + tmp = inv_sq_pi_2 * dexp(-mu_tmp*mu_tmp) / (mu_val * mu_val) + ! grad = + grad(1) = tmp * mu_der(1) + grad(2) = tmp * mu_der(2) + grad(3) = tmp * mu_der(3) + + if(r12 .lt. 1d-10) return + tmp = 0.5d0 * (1.d0 - derf(mu_tmp)) / r12 ! d/dx1 j(mu(r1),r12) + grad(1) = grad(1) + tmp * dx + grad(2) = grad(2) + tmp * dy + grad(3) = grad(3) + tmp * dz + + grad_jast = grad +end + +! --- + +double precision function j12_mu_input(r1, r2, mu) + + BEGIN_DOC + ! j(mu,r12) = 1/2 r12 (1 - erf(mu r12)) - 1/2 (sqrt(pi) * mu) e^{-(mu*r12)^2} + END_DOC + include 'constants.include.F' + + implicit none + double precision, intent(in) :: r1(3), r2(3), mu + double precision :: mu_tmp, r12 + + r12 = dsqrt( (r1(1) - r2(1)) * (r1(1) - r2(1)) & + + (r1(2) - r2(2)) * (r1(2) - r2(2)) & + + (r1(3) - r2(3)) * (r1(3) - r2(3)) ) + mu_tmp = mu * r12 + + j12_mu_input = 0.5d0 * r12 * (1.d0 - derf(mu_tmp)) - inv_sq_pi_2 * dexp(-mu_tmp*mu_tmp) / mu + +end + +subroutine grad_j12_mu_input(r1, r2, mu, grad_jmu,r12) + implicit none + BEGIN_DOC + ! grad_jmu = d/dx1 j(mu,r12) assuming mu=cst(r1) + ! + ! = 0.5/r_12 * (x_1 - x_2) * [1 - erf(mu*r12)] + END_DOC + double precision, intent(in) :: r1(3), r2(3), mu + double precision, intent(out):: grad_jmu(3),r12 + double precision :: mu_tmp, dx, dy, dz, grad(3), tmp + grad_jmu = 0.d0 + dx = r1(1) - r2(1) + dy = r1(2) - r2(2) + dz = r1(3) - r2(3) + r12 = dsqrt(dx * dx + dy * dy + dz * dz) + if(r12 .lt. 1d-10) return + mu_tmp = mu * r12 + tmp = 0.5d0 * (1.d0 - derf(mu_tmp)) / r12 ! d/dx1 j(mu(r1),r12) + grad(1) = tmp * dx + grad(2) = tmp * dy + grad(3) = tmp * dz + + grad_jmu = grad +end + +subroutine j12_and_grad_j12_mu_input(r1, r2, mu, jmu, grad_jmu) + implicit none + include 'constants.include.F' + BEGIN_DOC + ! jmu = j(mu,r12) + ! grad_jmu = d/dx1 j(mu,r12) assuming mu=cst(r1) + ! + ! = 0.5/r_12 * (x_1 - x_2) * [1 - erf(mu*r12)] + END_DOC + double precision, intent(in) :: r1(3), r2(3), mu + double precision, intent(out):: grad_jmu(3), jmu + double precision :: mu_tmp, r12, dx, dy, dz, grad(3), tmp + double precision :: erfc_mur12,inv_mu + inv_mu = 1.d0/mu + + grad_jmu = 0.d0 ! initialization when r12 --> 0 + jmu = - inv_sq_pi_2 * inv_mu ! initialization when r12 --> 0 + + dx = r1(1) - r2(1) + dy = r1(2) - r2(2) + dz = r1(3) - r2(3) + r12 = dsqrt(dx * dx + dy * dy + dz * dz) + if(r12 .lt. 1d-10) return + erfc_mur12 = (1.d0 - derf(mu_tmp)) + mu_tmp = mu * r12 + tmp = 0.5d0 * erfc_mur12 / r12 ! d/dx1 j(mu(r1),r12) + grad(1) = tmp * dx + grad(2) = tmp * dy + grad(3) = tmp * dz + + grad_jmu = grad + + jmu= 0.5d0 * r12 * erfc_mur12 - inv_sq_pi_2 * dexp(-mu_tmp*mu_tmp) * inv_mu + + +end diff --git a/plugins/local/non_h_ints_mu/jast_deriv_utils.irp.f b/plugins/local/non_h_ints_mu/jast_deriv_utils.irp.f index 79822508..f13990de 100644 --- a/plugins/local/non_h_ints_mu/jast_deriv_utils.irp.f +++ b/plugins/local/non_h_ints_mu/jast_deriv_utils.irp.f @@ -1,8 +1,73 @@ +subroutine get_deriv_r12_j12(x,mu,d_dx_j) + implicit none + include 'constants.include.F' + BEGIN_DOC + ! d/dr12 j(mu,r12) + END_DOC + double precision, intent(in) :: x,mu + double precision, intent(out) :: d_dx_j + + d_dx_j = 0.d0 + if(x .lt. 1d-10) return + if(j2e_type .eq. "Mu" .or. j2e_type .eq. "Mur") then + d_dx_j = 0.5d0 * (1.d0 - derf(mu * x)) + else if(j2e_type .eq. "Mugauss" .or. j2e_type .eq. "Murgauss" ) then + double precision :: x_tmp + x_tmp = mu * x + ! gradient of j(mu,x) + d_dx_j = 0.5d0 * (1.d0 - derf(x_tmp)) + + ! gradient of gaussian additional term + x_tmp *= alpha_mu_gauss + x_tmp *= x_tmp + d_dx_j += -0.5d0 * mu * c_mu_gauss * x * dexp(-x_tmp) + else + print *, ' Error in get_deriv_r12_j12: Unknown j2e_type = ', j2e_type + stop + endif +end + + +subroutine get_deriv_mu_j12(x,mu,d_d_mu) + implicit none + BEGIN_DOC + ! d/dmu j(mu,r12) + END_DOC + include 'constants.include.F' + double precision, intent(in) :: x,mu + double precision, intent(out) :: d_d_mu + double precision :: x_tmp,inv_mu_2,inv_alpha_2 + + d_d_mu = 0.d0 + if(x .lt. 1d-10) return + x_tmp = x*mu + if(mu.lt.1.d-10) return + inv_mu_2 = mu*mu + inv_mu_2 = 1.d0/inv_mu_2 + if(j2e_type .eq. "Mu" .or. j2e_type .eq. "Mur") then + ! e^{-(r12*mu)^2}/(2 sqrt(pi) * mu^2) + d_d_mu = dexp(-x_tmp*x_tmp) * inv_sq_pi_2 * inv_mu_2 + else if(j2e_type .eq. "Mugauss" .or. j2e_type .eq. "Murgauss" ) then + d_d_mu = dexp(-x_tmp*x_tmp) * inv_sq_pi_2 * inv_mu_2 + inv_alpha_2 = 1.d0/alpha_mu_gauss + inv_alpha_2 *= inv_alpha_2 + x_tmp *= alpha_mu_gauss + x_tmp *= x_tmp + d_d_mu += -0.25d0 * c_mu_gauss*inv_alpha_2*dexp(-x_tmp) * (1.d0 + 2.d0 * x_tmp) * inv_mu_2 + else + print *, ' Error in get_deriv_r12_j12: Unknown j2e_type = ', j2e_type + stop + endif +end + ! --- double precision function j12_mu(r1, r2) + BEGIN_DOC + ! j(mu,r12) = 1/2 r12 (1 - erf(mu r12)) - 1/2 (sqrt(pi) * mu) e^{-(mu*r12)^2} + END_DOC include 'constants.include.F' implicit none @@ -18,6 +83,18 @@ double precision function j12_mu(r1, r2) j12_mu = 0.5d0 * r12 * (1.d0 - derf(mu_tmp)) - inv_sq_pi_2 * dexp(-mu_tmp*mu_tmp) / mu_erf + else if(j2e_type .eq. "Mugauss") then + + r12 = dsqrt( (r1(1) - r2(1)) * (r1(1) - r2(1)) & + + (r1(2) - r2(2)) * (r1(2) - r2(2)) & + + (r1(3) - r2(3)) * (r1(3) - r2(3)) ) + double precision :: r12_tmp + r12_tmp = mu_erf * r12 + + j12_mu = 0.5d0 * r12 * (1.d0 - derf(r12_tmp)) - inv_sq_pi_2 * dexp(-r12_tmp*r12_tmp) / mu_erf + r12_tmp *= alpha_mu_gauss + j12_mu += 0.25d0 * c_mu_gauss / (alpha_mu_gauss*alpha_mu_gauss*mu_erf) * dexp(-r12_tmp*r12_tmp) + else print *, ' Error in j12_mu: Unknown j2e_type = ', j2e_type @@ -57,7 +134,7 @@ subroutine grad1_j12_mu(r1, r2, grad) grad = 0.d0 - if(j2e_type .eq. "Mu") then + if(j2e_type .eq. "Mu".or.j2e_type .eq. "Mugauss") then dx = r1(1) - r2(1) dy = r1(2) - r2(2) @@ -66,31 +143,42 @@ subroutine grad1_j12_mu(r1, r2, grad) r12 = dsqrt(dx * dx + dy * dy + dz * dz) if(r12 .lt. 1d-10) return - tmp = 0.5d0 * (1.d0 - derf(mu_erf * r12)) / r12 + call get_deriv_r12_j12(r12,mu_erf,tmp) +! tmp = 0.5d0 * (1.d0 - derf(mu_erf * r12)) / r12 grad(1) = tmp * dx grad(2) = tmp * dy grad(3) = tmp * dz + grad *= 1.d0/r12 - elseif(j2e_type .eq. "Mur") then + elseif(j2e_type .eq. "Mur" .or. j2e_type .eq. "Murgauss") then + double precision :: jast + call grad_j_sum_mu_of_r(r1,r2,jast,grad) + + elseif(j2e_type .eq. "Bump") then + double precision ::grad_jast(3) + call get_grad_j_bump_mu_of_r(r1,r2,grad_jast) + dx = r1(1) - r2(1) + dy = r1(2) - r2(2) + dz = r1(3) - r2(3) - dx = r1(1) - r2(1) - dy = r1(2) - r2(2) - dz = r1(3) - r2(3) - r12 = dsqrt(dx * dx + dy * dy + dz * dz) + r12 = dsqrt(dx * dx + dy * dy + dz * dz) + if(r12 .lt. 1d-10) then + grad(1) = 0.d0 + grad(2) = 0.d0 + grad(3) = 0.d0 + return + endif - call mu_r_val_and_grad(r1, r2, mu_val, mu_der) - mu_tmp = mu_val * r12 - tmp = inv_sq_pi_2 * dexp(-mu_tmp*mu_tmp) / (mu_val * mu_val) - grad(1) = tmp * mu_der(1) - grad(2) = tmp * mu_der(2) - grad(3) = tmp * mu_der(3) + tmp = 0.5d0 * (1.d0 - derf(mu_erf * r12)) / r12 + + grad(1) = 0.5d0 * tmp * dx + grad(2) = 0.5d0 * tmp * dy + grad(3) = 0.5d0 * tmp * dz + grad(1) += 0.5d0 * grad_jast(1) + grad(2) += 0.5d0 * grad_jast(2) + grad(3) += 0.5d0 * grad_jast(3) - if(r12 .lt. 1d-10) return - tmp = 0.5d0 * (1.d0 - derf(mu_tmp)) / r12 - grad(1) = grad(1) + tmp * dx - grad(2) = grad(2) + tmp * dy - grad(3) = grad(3) + tmp * dz else @@ -369,7 +457,18 @@ end ! --- subroutine mu_r_val_and_grad(r1, r2, mu_val, mu_der) - + BEGIN_DOC +! various flavours of mu(r1,r2) +! depends on essentially the density and other related quantities +! +! change the variable "murho_type" to change type +! +! murho_type == -1 :: mu(r1,r2) = (rho(r1) mu_mf(r1) + rho(r2) mu_mf(r2))/[rho(r1)+rho(r2)] +! +! == 0 :: mu(r1,r2) = (sqrt(rho(r1)) mu_mf(r1) + sqrt(rho(r2)) mu_mf(r2))/[sqrt(rho(r1))+sqrt(rho(r2))] +! +! == -2 :: mu(r1,r2) = 0.5(mu_mf(r1) + mu_mf(r2)) + END_DOC implicit none double precision, intent(in) :: r1(3), r2(3) double precision, intent(out) :: mu_val, mu_der(3) @@ -379,11 +478,50 @@ subroutine mu_r_val_and_grad(r1, r2, mu_val, mu_der) double precision :: rho1, grad_rho1(3),rho2,rho_tot,inv_rho_tot double precision :: f_rho1, f_rho2, d_drho_f_rho1 double precision :: d_dx1_f_rho1(3),d_dx_rho_f_rho(3),nume + double precision :: mu_mf_r1, dm_r1, grad_mu_mf_r1(3), grad_dm_r1(3) + double precision :: mu_mf_r2, dm_r2, grad_mu_mf_r2(3), grad_dm_r2(3) + + double precision :: num, denom, grad_denom(3), grad_num(3) + double precision :: dsqrt_dm_r1 PROVIDE murho_type PROVIDE mu_r_ct mu_erf - if(murho_type .eq. 1) then + if(murho_type .eq. 0) then + call grad_mu_of_r_mean_field(r1,mu_mf_r1, dm_r1, grad_mu_mf_r1, grad_dm_r1) + call grad_mu_of_r_mean_field(r2,mu_mf_r2, dm_r2, grad_mu_mf_r2, grad_dm_r2) + dsqrt_dm_r1 = dsqrt(dm_r1) + denom = (dsqrt_dm_r1 + dsqrt(dm_r2) ) + if(denom.lt.1.d-7)then + mu_val = 1.d+10 + mu_der = 0.d0 + return + endif + num = (dsqrt(dm_r1) * mu_mf_r1 + dsqrt(dm_r2) * mu_mf_r2) + mu_val = num / denom + grad_denom = grad_dm_r1/dsqrt_dm_r1 + grad_num = dsqrt(dm_r1) * grad_mu_mf_r1 + mu_mf_r1 * grad_dm_r1 + mu_der = (grad_num * denom - num * grad_denom)/(denom*denom) + else if(murho_type .eq. -1) then + call grad_mu_of_r_mean_field(r1,mu_mf_r1, dm_r1, grad_mu_mf_r1, grad_dm_r1) + call grad_mu_of_r_mean_field(r2,mu_mf_r2, dm_r2, grad_mu_mf_r2, grad_dm_r2) + denom = (dm_r1 + dm_r2 ) + if(denom.lt.1.d-7)then + mu_val = 1.d+10 + mu_der = 0.d0 + return + endif + num = (dm_r1 * mu_mf_r1 + dm_r2 * mu_mf_r2) + mu_val = num / denom + grad_denom = grad_dm_r1 + grad_num = dm_r1 * grad_mu_mf_r1 + mu_mf_r1 * grad_dm_r1 + mu_der = (grad_num * denom - num * grad_denom)/(denom*denom) + else if(murho_type .eq. -2) then + call grad_mu_of_r_mean_field(r1,mu_mf_r1, dm_r1, grad_mu_mf_r1, grad_dm_r1) + call grad_mu_of_r_mean_field(r2,mu_mf_r2, dm_r2, grad_mu_mf_r2, grad_dm_r2) + mu_val = 0.5d0 * (mu_mf_r1 + mu_mf_r2) + mu_der = 0.5d0 * grad_mu_mf_r1 + else if(murho_type .eq. 1) then ! ! r = 0.5 (r1 + r2) diff --git a/plugins/local/non_h_ints_mu/jast_deriv_utils_vect.irp.f b/plugins/local/non_h_ints_mu/jast_deriv_utils_vect.irp.f index 4fc537c8..d951db93 100644 --- a/plugins/local/non_h_ints_mu/jast_deriv_utils_vect.irp.f +++ b/plugins/local/non_h_ints_mu/jast_deriv_utils_vect.irp.f @@ -33,8 +33,12 @@ subroutine get_grad1_u12_withsq_r1_seq(ipoint, n_grid2, resx, resy, resz, res) r1(2) = final_grid_points(2,ipoint) r1(3) = final_grid_points(3,ipoint) - if( (j2e_type .eq. "Mu") .or. & - (j2e_type .eq. "Mur") .or. & + if( (j2e_type .eq. "Mu") .or. & + (j2e_type .eq. "Mur") .or. & + (j2e_type .eq. "Jpsi") .or. & + (j2e_type .eq. "Mugauss") .or. & + (j2e_type .eq. "Murgauss") .or. & + (j2e_type .eq. "Bump") .or. & (j2e_type .eq. "Boys") ) then if(env_type .eq. "None") then @@ -206,7 +210,43 @@ subroutine grad1_j12_r1_seq(r1, n_grid2, gradx, grady, gradz) gradz(jpoint) = tmp * dz enddo - elseif(j2e_type .eq. "Mur") then + else if(j2e_type .eq. "Mugauss") then + + ! d/dx1 j(mu,r12) = 0.5 * [(1 - erf(mu * r12)) / r12 - mu*c*r12*exp(-(mu*alpha*r12)^2] * (x1 - x2) + + do jpoint = 1, n_points_extra_final_grid ! r2 + + r2(1) = final_grid_points_extra(1,jpoint) + r2(2) = final_grid_points_extra(2,jpoint) + r2(3) = final_grid_points_extra(3,jpoint) + + dx = r1(1) - r2(1) + dy = r1(2) - r2(2) + dz = r1(3) - r2(3) + + r12 = dsqrt(dx * dx + dy * dy + dz * dz) + if(r12 .lt. 1d-10) then + gradx(jpoint) = 0.d0 + grady(jpoint) = 0.d0 + gradz(jpoint) = 0.d0 + cycle + endif + + double precision :: r12_tmp + r12_tmp = mu_erf * r12 + ! gradient of j(mu,r12) + tmp = 0.5d0 * (1.d0 - derf(r12_tmp)) / r12 + ! gradient of gaussian additional term + r12_tmp *= alpha_mu_gauss + r12_tmp *= r12_tmp + tmp += -0.5d0 * mu_erf * c_mu_gauss * r12 * dexp(-r12_tmp)/r12 + + gradx(jpoint) = tmp * dx + grady(jpoint) = tmp * dy + gradz(jpoint) = tmp * dz + enddo + + elseif(j2e_type .eq. "Mur".or.j2e_type .eq. "Murgauss") then ! d/dx1 j(mu(r1,r2),r12) = exp(-(mu(r1,r2)*r12)**2) /(2 *sqrt(pi) * mu(r1,r2)**2 ) d/dx1 mu(r1,r2) ! + 0.5 * (1 - erf(mu(r1,r2) *r12))/r12 * (x1 - x2) @@ -216,31 +256,46 @@ subroutine grad1_j12_r1_seq(r1, n_grid2, gradx, grady, gradz) r2(1) = final_grid_points_extra(1,jpoint) r2(2) = final_grid_points_extra(2,jpoint) r2(3) = final_grid_points_extra(3,jpoint) + double precision :: jast, grad_jast(3) + call grad_j_sum_mu_of_r(r1,r2,jast,grad_jast) + gradx(jpoint) = grad_jast(1) + grady(jpoint) = grad_jast(2) + gradz(jpoint) = grad_jast(3) + enddo + elseif(j2e_type .eq. "Bump") then + + ! d/dx1 jbump(r1,r2) + + do jpoint = 1, n_points_extra_final_grid ! r2 + + r2(1) = final_grid_points_extra(1,jpoint) + r2(2) = final_grid_points_extra(2,jpoint) + r2(3) = final_grid_points_extra(3,jpoint) + call get_grad_j_bump_mu_of_r(r1,r2,grad_jast) + + dx = r1(1) - r2(1) + dy = r1(2) - r2(2) + dz = r1(3) - r2(3) - dx = r1(1) - r2(1) - dy = r1(2) - r2(2) - dz = r1(3) - r2(3) r12 = dsqrt(dx * dx + dy * dy + dz * dz) - - call mu_r_val_and_grad(r1, r2, mu_val, mu_der) - mu_tmp = mu_val * r12 - tmp = inv_sq_pi_2 * dexp(-mu_tmp*mu_tmp) / (mu_val * mu_val) - gradx(jpoint) = tmp * mu_der(1) - grady(jpoint) = tmp * mu_der(2) - gradz(jpoint) = tmp * mu_der(3) - if(r12 .lt. 1d-10) then - gradx(jpoint) = 0.d0 - grady(jpoint) = 0.d0 - gradz(jpoint) = 0.d0 + gradx(jpoint) = 0.d0 + grady(jpoint) = 0.d0 + gradz(jpoint) = 0.d0 cycle endif - tmp = 0.5d0 * (1.d0 - derf(mu_tmp)) / r12 + tmp = 0.5d0 * (1.d0 - derf(mu_erf * r12)) / r12 - gradx(jpoint) = gradx(jpoint) + tmp * dx - grady(jpoint) = grady(jpoint) + tmp * dy - gradz(jpoint) = gradz(jpoint) + tmp * dz + gradx(jpoint) = 0.5d0 * tmp * dx + grady(jpoint) = 0.5d0 * tmp * dy + gradz(jpoint) = 0.5d0 * tmp * dz + gradx(jpoint) += 0.5d0 * grad_jast(1) + grady(jpoint) += 0.5d0 * grad_jast(2) + gradz(jpoint) += 0.5d0 * grad_jast(3) +! gradx(jpoint) = grad_jast(1) +! grady(jpoint) = grad_jast(2) +! gradz(jpoint) = grad_jast(3) enddo elseif(j2e_type .eq. "Boys") then @@ -363,6 +418,17 @@ subroutine grad1_j12_r1_seq(r1, n_grid2, gradx, grady, gradz) enddo ! i_nucl enddo ! jpoint + elseif(j2e_type .eq. "Jpsi") then + double precision :: grad_j_psi_r1(3),jast_psi + do jpoint = 1, n_points_extra_final_grid ! r2 + r2(1) = final_grid_points_extra(1,jpoint) + r2(2) = final_grid_points_extra(2,jpoint) + r2(3) = final_grid_points_extra(3,jpoint) + call get_grad_r1_jastrow_psi(r1,r2,grad_j_psi_r1,jast_psi) + gradx(jpoint) = grad_j_psi_r1(1) + grady(jpoint) = grad_j_psi_r1(2) + gradz(jpoint) = grad_j_psi_r1(3) + enddo else print *, ' Error in grad1_j12_r1_seq: Unknown j2e_type = ', j2e_type @@ -666,8 +732,12 @@ subroutine get_grad1_u12_2e_r1_seq(ipoint, n_grid2, resx, resy, resz) r1(2) = final_grid_points(2,ipoint) r1(3) = final_grid_points(3,ipoint) - if( (j2e_type .eq. "Mu") .or. & - (j2e_type .eq. "Mur") .or. & + if( (j2e_type .eq. "Mu") .or. & + (j2e_type .eq. "Mugauss") .or. & + (j2e_type .eq. "Mur") .or. & + (j2e_type .eq. "Jpsi") .or. & + (j2e_type .eq. "Murgauss") .or. & + (j2e_type .eq. "Bump") .or. & (j2e_type .eq. "Boys") ) then if(env_type .eq. "None") then @@ -784,8 +854,11 @@ subroutine get_u12_2e_r1_seq(ipoint, n_grid2, res) r1(2) = final_grid_points(2,ipoint) r1(3) = final_grid_points(3,ipoint) - if( (j2e_type .eq. "Mu") .or. & - (j2e_type .eq. "Mur") .or. & + if( (j2e_type .eq. "Mu") .or. & + (j2e_type .eq. "Mur") .or. & + (j2e_type .eq. "Mugauss") .or. & + (j2e_type .eq. "Murgauss") .or. & + (j2e_type .eq. "Mugauss") .or. & (j2e_type .eq. "Boys") ) then if(env_type .eq. "None") then diff --git a/plugins/local/non_h_ints_mu/jastrow_psi.irp.f b/plugins/local/non_h_ints_mu/jastrow_psi.irp.f new file mode 100644 index 00000000..4c88e793 --- /dev/null +++ b/plugins/local/non_h_ints_mu/jastrow_psi.irp.f @@ -0,0 +1,124 @@ +BEGIN_PROVIDER [ double precision, c_ij_ab_jastrow, (mo_num, mo_num, elec_alpha_num, elec_beta_num)] + implicit none + integer :: iunit, getUnitAndOpen + c_ij_ab_jastrow = 0.d0 + iunit = getUnitAndOpen(trim(ezfio_work_dir)//'c_ij_ab', 'R') + read(iunit) c_ij_ab_jastrow + close(iunit) + print*,'c_ij_ab_jastrow = ' + integer :: i,j,a,b + do i = 1, elec_beta_num ! r2 + do j = 1, elec_alpha_num ! r1 + do a = elec_beta_num+1, mo_num ! r2 + do b = elec_alpha_num+1, mo_num ! r1 +! print*,b,a,j,i + print*,c_ij_ab_jastrow(b,a,j,i),b,a,j,i + if(dabs(c_ij_ab_jastrow(b,a,j,i)).lt.1.d-12)then + c_ij_ab_jastrow(b,a,j,i) = 0.d0 + endif + enddo + enddo + enddo + enddo +END_PROVIDER + +double precision function jastrow_psi(r1,r2) + implicit none + double precision, intent(in) :: r1(3), r2(3) + integer :: i,j,a,b + double precision, allocatable :: mos_array_r1(:), mos_array_r2(:) + allocate(mos_array_r1(mo_num), mos_array_r2(mo_num)) + call give_all_mos_at_r(r1,mos_array_r1) + call give_all_mos_at_r(r2,mos_array_r2) + double precision :: eps,coef, numerator,denominator + double precision :: phi_i_phi_j + eps = a_boys + jastrow_psi= 0.d0 + do i = 1, elec_beta_num ! r1 + do j = 1, elec_alpha_num ! r2 + phi_i_phi_j = mos_array_r1(i) * mos_array_r2(j) + eps + denominator = 1.d0/phi_i_phi_j + do a = elec_beta_num+1, mo_num ! r1 + do b = elec_alpha_num+1, mo_num ! r2 + coef = c_ij_ab_jastrow(b,a,j,i) + numerator = mos_array_r2(b) * mos_array_r1(a) + jastrow_psi += coef * numerator*denominator + enddo + enddo + enddo + enddo +end + +subroutine get_grad_r1_jastrow_psi(r1,r2,grad_j_psi_r1,jast) + implicit none + double precision, intent(in) :: r1(3), r2(3) + double precision, intent(out):: grad_j_psi_r1(3),jast + integer :: i,j,a,b + double precision, allocatable :: mos_array_r1(:), mos_array_r2(:) + double precision, allocatable :: mos_grad_array_r1(:,:),mos_grad_array_r2(:,:) + double precision :: num_j, denom_j, num_j_grad(3), denom_j_grad(3),delta,coef + double precision :: inv_denom_j + allocate(mos_array_r1(mo_num), mos_array_r2(mo_num)) + allocate(mos_grad_array_r1(3,mo_num), mos_grad_array_r2(3,mo_num)) + delta = a_boys + call give_all_mos_and_grad_at_r(r1,mos_array_r1,mos_grad_array_r1) + call give_all_mos_and_grad_at_r(r2,mos_array_r2,mos_grad_array_r2) + grad_j_psi_r1 = 0.d0 + jast = 0.d0 + do i = 1, elec_beta_num ! r1 + do j = 1, elec_alpha_num ! r2 + call denom_jpsi(i,j,delta,mos_array_r1,mos_grad_array_r1,mos_array_r2,denom_j, denom_j_grad) + inv_denom_j = 1.d0/denom_j + do a = elec_beta_num+1, mo_num ! r1 + do b = elec_alpha_num+1, mo_num ! r2 + call numerator_psi(a,b,mos_array_r1,mos_grad_array_r1,mos_array_r2,num_j, num_j_grad) + coef = c_ij_ab_jastrow(b,a,j,i) + jast += coef * num_j * inv_denom_j + grad_j_psi_r1 += coef * (num_j_grad * denom_j - num_j * denom_j_grad) * inv_denom_j * inv_denom_j + enddo + enddo + enddo + enddo + if(jast.lt.-1.d0.or.dabs(jast).gt.1.d0)then + print*,'pb ! ' + print*,jast + print*,dsqrt(r1(1)**2+r1(2)**2+r1(3)**2),dsqrt(r2(1)**2+r2(2)**2+r2(3)**2) + print*,r1 +! print*,mos_array_r1(1:2) + print*,r2 +! print*,mos_array_r2(1:2) + stop + endif + if(log_jpsi)then + grad_j_psi_r1 = grad_j_psi_r1/(1.d0 + jast) + endif + +end + + +subroutine denom_jpsi(i,j,delta,mos_array_r1,mos_grad_array_r1,mos_array_r2,denom, grad_denom) + implicit none + integer, intent(in) :: i,j + double precision, intent(in) :: mos_array_r1(mo_num),mos_grad_array_r1(3,mo_num),mos_array_r2(mo_num),delta + double precision, intent(out) :: denom, grad_denom(3) + double precision :: coef,phi_i_phi_j,inv_phi_i_phi_j,inv_phi_i_phi_j_2 + phi_i_phi_j = mos_array_r1(i) * mos_array_r2(j) + if(phi_i_phi_j /= 0.d0)then + inv_phi_i_phi_j = 1.d0/phi_i_phi_j + inv_phi_i_phi_j_2 = 1.d0/(phi_i_phi_j * phi_i_phi_j) + else + inv_phi_i_phi_j = huge(1.0) + inv_phi_i_phi_j_2 = huge(1.d0) + endif + denom = phi_i_phi_j + delta * inv_phi_i_phi_j + grad_denom(:) = (1.d0 - delta*inv_phi_i_phi_j_2) * mos_array_r2(j) * mos_grad_array_r1(:,i) +end + +subroutine numerator_psi(a,b,mos_array_r1,mos_grad_array_r1,mos_array_r2,num, grad_num) + implicit none + integer, intent(in) :: a,b + double precision, intent(in) :: mos_array_r1(mo_num),mos_grad_array_r1(3,mo_num),mos_array_r2(mo_num) + double precision, intent(out) :: num, grad_num(3) + num = mos_array_r1(a) * mos_array_r2(b) + grad_num(:) = mos_array_r2(b) * mos_grad_array_r1(:,a) +end diff --git a/plugins/local/non_h_ints_mu/mu_of_r.irp.f b/plugins/local/non_h_ints_mu/mu_of_r.irp.f new file mode 100644 index 00000000..e9cf6c4a --- /dev/null +++ b/plugins/local/non_h_ints_mu/mu_of_r.irp.f @@ -0,0 +1,43 @@ + +subroutine grad_mu_of_r_mean_field(r,mu_mf, dm, grad_mu_mf, grad_dm) + implicit none + BEGIN_DOC + ! returns the value and gradients of the mu(r) mean field, together with the HF density and its gradients. + END_DOC + include 'constants.include.F' + double precision, intent(in) :: r(3) + double precision, intent(out):: grad_mu_mf(3), grad_dm(3) + double precision, intent(out):: mu_mf, dm + double precision :: grad_f_mf_ab(3), grad_two_bod_dens(3),grad_dm_a(3), grad_dm_b(3) + double precision :: f_mf_ab,two_bod_dens, dm_a, dm_b + + double precision :: dist + call get_grad_f_mf_ab(r,grad_f_mf_ab, grad_two_bod_dens,f_mf_ab,two_bod_dens, dm_a, dm_b,grad_dm_a, grad_dm_b) + + dm = dm_a + dm_b + grad_dm(1:3) = grad_dm_a(1:3) + grad_dm_b(1:3) + + if(dabs(two_bod_dens).lt.1.d-10)then + mu_mf = 1.d+10 + grad_mu_mf = 0.d0 + else + if(mu_of_r_tc=="Erfmu")then + mu_mf = 0.3333333333d0 * sqpi * (f_mf_ab/two_bod_dens + 0.25d0) + grad_mu_mf(1:3) = 0.3333333333d0 * sqpi * (grad_f_mf_ab(1:3) * two_bod_dens - f_mf_ab * grad_two_bod_dens(1:3))& + /(two_bod_dens*two_bod_dens) + else if(mu_of_r_tc=="Standard")then + mu_mf = 0.5d0 * sqpi * f_mf_ab/two_bod_dens + grad_mu_mf(1:3) = 0.5d0 * sqpi * (grad_f_mf_ab(1:3) * two_bod_dens - f_mf_ab * grad_two_bod_dens(1:3))& + /(two_bod_dens*two_bod_dens) + else if(mu_of_r_tc=="Erfmugauss")then + mu_mf = (f_mf_ab/two_bod_dens + 0.25d0)/c_mu_gauss_tot + grad_mu_mf(1:3) = 1.d0/c_mu_gauss_tot* (grad_f_mf_ab(1:3) * two_bod_dens - f_mf_ab * grad_two_bod_dens(1:3))& + /(two_bod_dens*two_bod_dens) + else + print*,'Wrong value for mu_of_r_tc !' + stop + endif + endif + +end + diff --git a/src/mu_of_r/mu_of_r_mean_field.irp.f b/plugins/local/non_h_ints_mu/mu_of_r_mean_field.irp.f similarity index 88% rename from src/mu_of_r/mu_of_r_mean_field.irp.f rename to plugins/local/non_h_ints_mu/mu_of_r_mean_field.irp.f index 6abc7e4f..295d58c2 100644 --- a/src/mu_of_r/mu_of_r_mean_field.irp.f +++ b/plugins/local/non_h_ints_mu/mu_of_r_mean_field.irp.f @@ -57,6 +57,9 @@ end subroutine get_grad_f_mf_ab(r,grad_f_mf_ab, grad_two_bod_dens,f_mf_ab,two_bod_dens, dm_a, dm_b,grad_dm_a, grad_dm_b) implicit none + BEGIN_DOC + ! gradient of mu(r) mean field, together with the gradient of the one- and two-body HF density. + END_DOC double precision, intent(in) :: r(3) double precision, intent(out) :: f_mf_ab, two_bod_dens double precision, intent(out) :: grad_two_bod_dens(3), grad_f_mf_ab(3) @@ -146,26 +149,18 @@ subroutine mu_of_r_mean_field(r,mu_mf, dm) endif end -subroutine grad_mu_of_r_mean_field(r,mu_mf, dm, grad_mu_mf, grad_dm) +subroutine mu_of_r_mean_field_tc(r,mu_mf, dm) implicit none - include 'constants.include.F' + include 'constants.include.F' double precision, intent(in) :: r(3) - double precision, intent(out):: grad_mu_mf(3), grad_dm(3) double precision, intent(out):: mu_mf, dm - double precision :: grad_f_mf_ab(3), grad_two_bod_dens(3),grad_dm_a(3), grad_dm_b(3) double precision :: f_mf_ab,two_bod_dens, dm_a, dm_b - call get_grad_f_mf_ab(r,grad_f_mf_ab, grad_two_bod_dens,f_mf_ab,two_bod_dens, dm_a, dm_b,grad_dm_a, grad_dm_b) - + call get_f_mf_ab(r,f_mf_ab,two_bod_dens, dm_a, dm_b) dm = dm_a + dm_b - grad_dm(1:3) = grad_dm_a(1:3) + grad_dm_b(1:3) - if(dabs(two_bod_dens).lt.1.d-10)then mu_mf = 1.d+10 - grad_mu_mf = 0.d0 else - mu_mf = 0.5d0 * sqpi * f_mf_ab/two_bod_dens - grad_mu_mf(1:3) = 0.5d0 * sqpi * (grad_f_mf_ab(1:3) * two_bod_dens - f_mf_ab * grad_two_bod_dens(1:3))& - /(two_bod_dens*two_bod_dens) - endif - + mu_mf = 0.3333333333d0 * sqpi * (f_mf_ab/two_bod_dens + 0.25d0) + endif end + diff --git a/plugins/local/non_h_ints_mu/plot_j_gauss.irp.f b/plugins/local/non_h_ints_mu/plot_j_gauss.irp.f new file mode 100644 index 00000000..a4030d8c --- /dev/null +++ b/plugins/local/non_h_ints_mu/plot_j_gauss.irp.f @@ -0,0 +1,59 @@ +program plot_j_gauss + implicit none + double precision :: xmin, xmax, x, dx + double precision :: mu_min, mu_max, mu, d_mu + double precision :: pot_j_gauss,j_mu_simple,j_gauss_simple,pot_j_mu + double precision, allocatable :: mu_tab(:),j_mu(:),j_mu_gauss(:) + double precision, allocatable :: w_mu(:), w_mu_gauss(:) + + character*(128) :: output + integer :: getUnitAndOpen + integer :: i_unit_output_wee_gauss,i_unit_output_wee_mu + integer :: i_unit_output_j_gauss,i_unit_output_j_mu + output=trim(ezfio_filename)//'.w_ee_mu_gauss' + i_unit_output_wee_gauss = getUnitAndOpen(output,'w') + output=trim(ezfio_filename)//'.w_ee_mu' + i_unit_output_wee_mu = getUnitAndOpen(output,'w') + output=trim(ezfio_filename)//'.j_mu_gauss' + i_unit_output_j_gauss = getUnitAndOpen(output,'w') + output=trim(ezfio_filename)//'.j_mu' + i_unit_output_j_mu = getUnitAndOpen(output,'w') + + integer :: npt, i, j, n_mu + n_mu = 3 + allocate(mu_tab(n_mu),j_mu(n_mu),j_mu_gauss(n_mu),w_mu(n_mu), w_mu_gauss(n_mu)) + mu_min = 0.5d0 + mu_max = 2.d0 + d_mu = (mu_max - mu_min)/dble(n_mu) + mu = mu_min + do i = 1, n_mu + mu_tab(i) = mu + print*,'mu = ',mu + mu += d_mu + enddo + mu_tab(1) = 0.9d0 + mu_tab(2) = 0.95d0 + mu_tab(3) = 1.d0 + + xmin = 0.01d0 + xmax = 10.d0 + npt = 1000 + dx = (xmax - xmin)/dble(npt) + x = xmin + do i = 1, npt + do j = 1, n_mu + mu = mu_tab(j) + w_mu_gauss(j) = pot_j_gauss(x,mu) + w_mu(j) = pot_j_mu(x,mu) + j_mu(j) = j_mu_simple(x,mu) + j_mu_gauss(j) = j_gauss_simple(x,mu) + j_mu(j) + enddo + write(i_unit_output_wee_gauss,'(100(F16.10,X))')x,w_mu_gauss(:) + write(i_unit_output_wee_mu,'(100(F16.10,X))')x,w_mu(:) + write(i_unit_output_j_gauss,'(100(F16.10,X))')x,j_mu_gauss(:) + write(i_unit_output_j_mu,'(100(F16.10,X))')x,j_mu(:) + x += dx + enddo + + +end diff --git a/plugins/local/non_h_ints_mu/plot_mo.irp.f b/plugins/local/non_h_ints_mu/plot_mo.irp.f new file mode 100644 index 00000000..e1ecc783 --- /dev/null +++ b/plugins/local/non_h_ints_mu/plot_mo.irp.f @@ -0,0 +1,19 @@ +program plot_mo + implicit none + integer :: i,npt + double precision :: xmin,xmax,dx,r(3) + double precision,allocatable :: mos_array(:) + allocate(mos_array(mo_num)) + npt = 10000 + xmin =0.d0 + xmax =10.d0 + dx=(xmax-xmin)/dble(npt) + r=0.d0 + r(1) = xmin + do i = 1, npt + call give_all_mos_at_r(r,mos_array) + write(33,'(100(F16.10,X))')r(1),mos_array(1),mos_array(2),mos_array(3) + r(1) += dx + enddo + +end diff --git a/plugins/local/non_h_ints_mu/plot_mu_of_r.irp.f b/plugins/local/non_h_ints_mu/plot_mu_of_r.irp.f index 3a5984bd..4a3ec0d5 100644 --- a/plugins/local/non_h_ints_mu/plot_mu_of_r.irp.f +++ b/plugins/local/non_h_ints_mu/plot_mu_of_r.irp.f @@ -16,15 +16,16 @@ subroutine routine_print integer :: ipoint,nx,i double precision :: xmax,xmin,r(3),dx,sigma double precision :: mu_val, mu_der(3),dm_a,dm_b,grad,grad_dm_a(3), grad_dm_b(3) - xmax = 5.D0 - xmin = -5.D0 + xmax = 3.9D0 + xmin = -3.9D0 nx = 10000 dx = (xmax - xmin)/dble(nx) r = 0.d0 r(1) = xmin do ipoint = 1, nx - call mu_r_val_and_grad(r, r, mu_val, mu_der) - call density_and_grad_alpha_beta(r,dm_a,dm_b, grad_dm_a, grad_dm_b) +! call mu_r_val_and_grad(r, r, mu_val, mu_der) + call grad_mu_of_r_mean_field(r,mu_val, dm_a, mu_der, grad_dm_a) +! call density_and_grad_alpha_beta(r,dm_a,dm_b, grad_dm_a, grad_dm_b) sigma = 0.d0 do i = 1,3 sigma += grad_dm_a(i)**2 @@ -32,7 +33,8 @@ subroutine routine_print sigma=dsqrt(sigma) grad = mu_der(1)**2 + mu_der(2)**2 + mu_der(3)**2 grad = dsqrt(grad) - write(i_unit_output,'(100(F16.7,X))')r(1),mu_val,dm_a+dm_b,grad,sigma/dm_a + print*,r(1),mu_val + write(i_unit_output,'(100(F16.7,X))')r(1),mu_val,dm_a,grad,sigma/dm_a r(1) += dx enddo end diff --git a/plugins/local/non_h_ints_mu/pot_j_gauss.irp.f b/plugins/local/non_h_ints_mu/pot_j_gauss.irp.f new file mode 100644 index 00000000..f9a0a7bc --- /dev/null +++ b/plugins/local/non_h_ints_mu/pot_j_gauss.irp.f @@ -0,0 +1,146 @@ +double precision function j_simple(x,mu) + implicit none + double precision, intent(in) :: x,mu + double precision :: j_mu_simple,j_gauss_simple + if(j2e_type .eq. "Mu".or.j2e_type .eq. "Mur") then + j_simple = j_mu_simple(x,mu) + else if(j2e_type .eq. "Mugauss".or.j2e_type .eq. "Murgauss") then + j_simple = j_gauss_simple(x,mu) + j_mu_simple(x,mu) + endif +end + + +double precision function j_mu_simple(x,mu) + implicit none + double precision, intent(in):: x,mu + include 'constants.include.F' + BEGIN_DOC +! j_mu(mu,x) = 0.5 x (1 - erf(mu x)) - 1/[2 sqrt(pi)mu] exp(-(x*mu)^2) + END_DOC + j_mu_simple = 0.5d0 * x * (1.D0 - derf(mu*x)) - 0.5d0 * inv_sq_pi/mu * dexp(-x*mu*x*mu) + +end + +double precision function j_gauss_simple(x,mu) + implicit none + double precision, intent(in):: x,mu + include 'constants.include.F' + BEGIN_DOC +! j_mu(mu,x) = c/[4 alpha^2 mu] exp(-(alpha * mu * x)^2) +! with c = 27/(8 sqrt(pi)), alpha=3/2 + END_DOC + double precision :: x_tmp + x_tmp = alpha_mu_gauss * mu * x + j_gauss_simple = 0.25d0 * c_mu_gauss / (alpha_mu_gauss*alpha_mu_gauss*mu) * dexp(-x_tmp*x_tmp) + +end + +double precision function j_mu_deriv(x,mu) + implicit none + BEGIN_DOC +! d/dx j_mu(mu,x) = d/dx 0.5 x (1 - erf(mu x)) - 1/[2 sqrt(pi)mu] exp(-(x*mu)^2) +! = 0.5*(1 - erf(mu x)) + END_DOC + include 'constants.include.F' + double precision, intent(in) :: x,mu + j_mu_deriv = 0.5d0 * (1.d0 - derf(mu*x)) +end + +double precision function j_mu_deriv_2(x,mu) + implicit none + BEGIN_DOC +! d^2/dx^2 j_mu(mu,x) = d^2/dx^2 0.5 x (1 - erf(mu x)) - 1/[2 sqrt(pi)mu] exp(-(x*mu)^2) +! = -mu/sqrt(pi) * exp(-(mu x)^2) + END_DOC + include 'constants.include.F' + double precision, intent(in) :: x,mu + j_mu_deriv_2 = - mu * inv_sq_pi * dexp(-x*mu*x*mu) +end + +double precision function j_gauss_deriv(x,mu) + implicit none + include 'constants.include.F' + double precision, intent(in) :: x,mu + BEGIN_DOC +! d/dx j_gauss(mu,x) = d/dx c/[4 alpha^2 mu] exp(-(alpha * mu * x)^2) +! with c = 27/(8 sqrt(pi)), alpha=3/2 +! = -0.5 * mu * c * x * exp(-(alpha * mu * x)^2) + END_DOC + double precision :: x_tmp + x_tmp = alpha_mu_gauss * mu * x + j_gauss_deriv = -0.5d0 * mu * c_mu_gauss * x * exp(-x_tmp*x_tmp) +end + +double precision function j_gauss_deriv_2(x,mu) + implicit none + include 'constants.include.F' + double precision, intent(in) :: x,mu + BEGIN_DOC +! d/dx j_gauss(mu,x) = d/dx c/[4 alpha^2 mu] exp(-(alpha * mu * x)^2) +! with c = 27/(8 sqrt(pi)), alpha=3/2 +! = 0.5 * mu * c * exp(-(alpha * mu * x)^2) * (2 (alpha*mu*x)^2 - 1) + END_DOC + double precision :: x_tmp + x_tmp = alpha_mu_gauss * mu * x + x_tmp = x_tmp * x_tmp + j_gauss_deriv_2 = 0.5d0 * mu * c_mu_gauss * exp(-x_tmp) * (2.d0*x_tmp - 1.d0) +end + +double precision function j_erf_gauss_deriv(x,mu) + implicit none + double precision, intent(in) :: x,mu + BEGIN_DOC +! d/dx (j_gauss(mu,x)+j_mu(mu,x)) + END_DOC + double precision :: j_gauss_deriv,j_mu_deriv + j_erf_gauss_deriv = j_gauss_deriv(x,mu)+j_mu_deriv(x,mu) +end + +double precision function j_erf_gauss_deriv_2(x,mu) + implicit none + double precision, intent(in) :: x,mu + BEGIN_DOC +! d^2/dx^2 (j_gauss(mu,x)+j_mu(mu,x)) + END_DOC + double precision :: j_gauss_deriv_2,j_mu_deriv_2 + j_erf_gauss_deriv_2 = j_gauss_deriv_2(x,mu)+j_mu_deriv_2(x,mu) +end + + +double precision function pot_j_gauss(x,mu) + implicit none + double precision, intent(in) :: x,mu + BEGIN_DOC + ! effective scalar potential associated with the erf_gauss correlation factor + ! + ! 1/x( 1 - 2 * d/dx j_erf_gauss(x,mu)) - d^2/dx^2 j_erf_gauss(x,mu)) - d/dx d/dx (j_erf_gauss(x,mu))^2 + END_DOC + double precision :: j_erf_gauss_deriv_2,j_erf_gauss_deriv + double precision :: deriv_1, deriv_2 + pot_j_gauss = 0.d0 + if(x.ne.0.d0)then + deriv_1 = j_erf_gauss_deriv(x,mu) + deriv_2 = j_erf_gauss_deriv_2(x,mu) + pot_j_gauss = 1.d0/x * (1.d0 - 2.d0 * deriv_1) - deriv_1 * deriv_1 - deriv_2 + endif + +end + +double precision function pot_j_mu(x,mu) + implicit none + double precision, intent(in) :: x,mu + BEGIN_DOC + ! effective scalar potential associated with the correlation factor + ! + ! 1/x( 1 - 2 * d/dx j_erf(x,mu)) - d^2/dx^2 j_erf(x,mu)) - d/dx d/dx (j_erf(x,mu))^2 + END_DOC + double precision :: j_mu_deriv_2,j_mu_deriv + double precision :: deriv_1, deriv_2 + pot_j_mu = 0.d0 + if(x.ne.0.d0)then + deriv_1 = j_mu_deriv(x,mu) + deriv_2 = j_mu_deriv_2(x,mu) + pot_j_mu= 1.d0/x * (1.d0 - 2.d0 * deriv_1) - deriv_1 * deriv_1 - deriv_2 + endif + +end diff --git a/plugins/local/non_h_ints_mu/print_jastrow_psi.irp.f b/plugins/local/non_h_ints_mu/print_jastrow_psi.irp.f new file mode 100644 index 00000000..740743cb --- /dev/null +++ b/plugins/local/non_h_ints_mu/print_jastrow_psi.irp.f @@ -0,0 +1,15 @@ +program print_j_psi + implicit none + integer :: i,j,a,b + do i = 1, elec_beta_num ! r2 + do j = 1, elec_alpha_num ! r1 + do a = elec_beta_num+1, mo_num ! r2 + do b = elec_alpha_num+1, mo_num ! r1 + print*,b,a,j,i + print*,c_ij_ab_jastrow(b,a,j,i) + enddo + enddo + enddo + enddo + +end diff --git a/plugins/local/non_h_ints_mu/tc_integ.irp.f b/plugins/local/non_h_ints_mu/tc_integ.irp.f index ce7ab101..ab88dd41 100644 --- a/plugins/local/non_h_ints_mu/tc_integ.irp.f +++ b/plugins/local/non_h_ints_mu/tc_integ.irp.f @@ -127,8 +127,8 @@ BEGIN_PROVIDER [double precision, int2_grad1_u12_ao, (ao_num, ao_num, n_points_f ! TODO combine 1shot & int2_grad1_u12_ao_num PROVIDE int2_grad1_u12_ao_num int2_grad1_u12_ao = int2_grad1_u12_ao_num - !PROVIDE int2_grad1_u12_ao_num_1shot - !int2_grad1_u12_ao = int2_grad1_u12_ao_num_1shot +! PROVIDE int2_grad1_u12_ao_num_1shot +! int2_grad1_u12_ao = int2_grad1_u12_ao_num_1shot endif elseif(tc_integ_type .eq. "semi-analytic") then diff --git a/plugins/local/non_h_ints_mu/test_mu_of_r_tc.irp.f b/plugins/local/non_h_ints_mu/test_mu_of_r_tc.irp.f new file mode 100644 index 00000000..79e21d1b --- /dev/null +++ b/plugins/local/non_h_ints_mu/test_mu_of_r_tc.irp.f @@ -0,0 +1,157 @@ +program test_mu_of_r_tc + implicit none + BEGIN_DOC +! TODO + END_DOC + ! You specify that you want to avoid any contribution from + ! orbitals coming from core + call test_grad_f_mean_field + call test_grad_mu_mf + call plot_mu_of_r_mf +end + + +subroutine test_grad_f_mean_field + implicit none + integer :: i_point,k + double precision :: weight,r(3) + double precision :: grad_f_mf_ab(3), grad_two_bod_dens(3) + double precision :: grad_dm_a(3), grad_dm_b(3) + double precision :: f_mf_ab,two_bod_dens, dm_a, dm_b + + double precision :: num_grad_f_mf_ab(3), num_grad_two_bod_dens(3) + double precision :: num_grad_dm_a(3), num_grad_dm_b(3) + double precision :: f_mf_ab_p,f_mf_ab_m + double precision :: two_bod_dens_p, two_bod_dens_m + double precision :: dm_a_p, dm_a_m + double precision :: dm_b_p, dm_b_m + double precision :: rbis(3), dr + double precision :: accu_grad_f_mf_ab(3),accu_grad_two_bod_dens(3) + double precision :: accu_grad_dm_a(3),accu_grad_dm_b(3) + double precision :: accu_f_mf_ab, accu_two_bod_dens, accu_dm_a, accu_dm_b + dr = 0.00001d0 + accu_f_mf_ab = 0.d0 + accu_two_bod_dens = 0.d0 + accu_dm_a = 0.d0 + accu_dm_b = 0.d0 + + accu_grad_f_mf_ab = 0.d0 + accu_grad_two_bod_dens = 0.d0 + accu_grad_dm_a = 0.d0 + accu_grad_dm_b = 0.d0 + do i_point = 1, n_points_final_grid + r(1:3) = final_grid_points(1:3,i_point) + weight = final_weight_at_r_vector(i_point) + call get_grad_f_mf_ab(r,grad_f_mf_ab, grad_two_bod_dens,f_mf_ab,two_bod_dens, dm_a, dm_b,grad_dm_a, grad_dm_b) + call get_f_mf_ab(r,f_mf_ab_p,two_bod_dens_p, dm_a_p, dm_b_p) + accu_f_mf_ab += weight * dabs(f_mf_ab - f_mf_ab_p) + accu_two_bod_dens += weight * dabs(two_bod_dens - two_bod_dens_p) + accu_dm_a += weight*dabs(dm_a - dm_a_p) + accu_dm_b += weight*dabs(dm_b - dm_b_p) + do k = 1, 3 + rbis = r + rbis(k) += dr + call get_f_mf_ab(rbis,f_mf_ab_p,two_bod_dens_p, dm_a_p, dm_b_p) + rbis = r + rbis(k) -= dr + call get_f_mf_ab(rbis,f_mf_ab_m,two_bod_dens_m, dm_a_m, dm_b_m) + num_grad_f_mf_ab(k) = (f_mf_ab_p - f_mf_ab_m)/(2.d0*dr) + num_grad_two_bod_dens(k) = (two_bod_dens_p - two_bod_dens_m)/(2.d0*dr) + num_grad_dm_a(k) = (dm_a_p - dm_a_m)/(2.d0*dr) + num_grad_dm_b(k) = (dm_b_p - dm_b_m)/(2.d0*dr) + enddo + do k = 1, 3 + accu_grad_f_mf_ab(k) += weight * dabs(grad_f_mf_ab(k) - num_grad_f_mf_ab(k)) + accu_grad_two_bod_dens(k) += weight * dabs(grad_two_bod_dens(k) - num_grad_two_bod_dens(k)) + accu_grad_dm_a(k) += weight * dabs(grad_dm_a(k) - num_grad_dm_a(k)) + accu_grad_dm_b(k) += weight * dabs(grad_dm_b(k) - num_grad_dm_b(k)) + enddo + enddo + print*,'accu_f_mf_ab = ',accu_f_mf_ab + print*,'accu_two_bod_dens = ',accu_two_bod_dens + print*,'accu_dm_a = ',accu_dm_a + print*,'accu_dm_b = ',accu_dm_b + print*,'accu_grad_f_mf_ab = ' + print*,accu_grad_f_mf_ab + print*,'accu_grad_two_bod_dens = ' + print*,accu_grad_two_bod_dens + print*,'accu_dm_a = ' + print*,accu_grad_dm_a + print*,'accu_dm_b = ' + print*,accu_grad_dm_b + +end + +subroutine test_grad_mu_mf + implicit none + integer :: i_point,k + double precision :: weight,r(3),rbis(3) + double precision :: mu_mf, dm,grad_mu_mf(3), grad_dm(3) + double precision :: mu_mf_p, mu_mf_m, dm_m, dm_p, num_grad_mu_mf(3),dr, num_grad_dm(3) + double precision :: accu_mu, accu_dm, accu_grad_dm(3), accu_grad_mu_mf(3) + dr = 0.00001d0 + accu_grad_mu_mf = 0.d0 + accu_mu = 0.d0 + accu_grad_dm = 0.d0 + accu_dm = 0.d0 + do i_point = 1, n_points_final_grid + r(1:3) = final_grid_points(1:3,i_point) + weight = final_weight_at_r_vector(i_point) + call grad_mu_of_r_mean_field(r,mu_mf, dm, grad_mu_mf, grad_dm) + call mu_of_r_mean_field(r,mu_mf_p, dm_p) + accu_mu += weight*dabs(mu_mf_p - mu_mf) + accu_dm += weight*dabs(dm_p - dm) + do k = 1, 3 + rbis = r + rbis(k) += dr + call mu_of_r_mean_field(rbis,mu_mf_p, dm_p) + rbis = r + rbis(k) -= dr + call mu_of_r_mean_field(rbis,mu_mf_m, dm_m) + + num_grad_mu_mf(k) = (mu_mf_p - mu_mf_m)/(2.d0*dr) + num_grad_dm(k) = (dm_p - dm_m)/(2.d0*dr) + enddo + do k = 1, 3 + accu_grad_dm(k)+= weight *dabs(num_grad_dm(k) - grad_dm(k)) + accu_grad_mu_mf(k)+= weight *dabs(num_grad_mu_mf(k) - grad_mu_mf(k)) + enddo + enddo + print*,'accu_mu = ',accu_mu + print*,'accu_dm = ',accu_dm + print*,'accu_grad_dm = ' + print*, accu_grad_dm + print*,'accu_grad_mu_mf = ' + print*, accu_grad_mu_mf + +end + +subroutine plot_mu_of_r_mf + implicit none + include 'constants.include.F' + integer :: ipoint,npoint + double precision :: dx,r(3),xmax,xmin + double precision :: accu_mu,accu_nelec,mu_mf, dm,mu_mf_tc + character*(128) :: output + integer :: i_unit_output,getUnitAndOpen + output=trim(ezfio_filename)//'.mu_mf' + i_unit_output = getUnitAndOpen(output,'w') + xmax = 5.D0 + xmin = 0.d0 + npoint = 10000 + dx = (xmax - xmin)/dble(npoint) + r = 0.d0 + r(1) = xmin + accu_mu = 0.d0 + accu_nelec = 0.d0 + do ipoint = 1, npoint + call mu_of_r_mean_field(r,mu_mf, dm) + call mu_of_r_mean_field_tc(r,mu_mf_tc, dm) + write(i_unit_output,'(100(F16.10,X))')r(1),mu_mf,mu_mf_tc,dm + accu_mu += mu_mf * dm * r(1)**2*dx*4.D0*pi + accu_nelec += dm * r(1)**2*dx*4.D0*pi + r(1) += dx + enddo + print*,'nelec = ',accu_nelec + print*,'mu average = ',accu_mu/accu_nelec +end diff --git a/plugins/local/tc_int/uninstall b/plugins/local/tc_int/uninstall index e37a5491..c4877879 100755 --- a/plugins/local/tc_int/uninstall +++ b/plugins/local/tc_int/uninstall @@ -9,5 +9,5 @@ then fi rm -rf ${PWD}/CuTC -rm ${QP_ROOT}/lib/libcutcint.so +rm -f ${QP_ROOT}/lib/libcutcint.so diff --git a/plugins/local/tc_keywords/EZFIO.cfg b/plugins/local/tc_keywords/EZFIO.cfg index f3bd75c8..b858fa5b 100644 --- a/plugins/local/tc_keywords/EZFIO.cfg +++ b/plugins/local/tc_keywords/EZFIO.cfg @@ -230,7 +230,7 @@ default: 70 type: character*(32) doc: approach used to evaluate TC integrals [ analytic | numeric | semi-analytic ] interface: ezfio,ocaml,provider -default: semi-analytic +default: numeric [minimize_lr_angles] type: logical diff --git a/scripts/import_champ_jastrow.py b/scripts/import_champ_jastrow.py old mode 100755 new mode 100644 index 489309b7..4bc58ff4 --- a/scripts/import_champ_jastrow.py +++ b/scripts/import_champ_jastrow.py @@ -45,7 +45,7 @@ if __name__ == '__main__': jastrow_file = sys.argv[2] jastrow = import_jastrow(jastrow_file) print (jastrow) - ezfio.set_jastrow_jast_type("Qmckl") + ezfio.set_jastrow_j2e_type("Qmckl") ezfio.set_jastrow_jast_qmckl_type_nucl_num(jastrow['type_num']) charges = ezfio.get_nuclei_nucl_charge() types = {} diff --git a/scripts/qp_import_trexio.py b/scripts/qp_import_trexio.py index 9251a1b0..50806598 100755 --- a/scripts/qp_import_trexio.py +++ b/scripts/qp_import_trexio.py @@ -261,13 +261,10 @@ def write_ezfio(trexio_filename, filename): except: cartesian = True - if not cartesian: - raise TypeError('Only cartesian TREXIO files can be converted') - ao_num = trexio.read_ao_num(trexio_file) ezfio.set_ao_basis_ao_num(ao_num) - if shell_num > 0: + if cartesian and shell_num > 0: ao_shell = trexio.read_ao_shell(trexio_file) at = [ nucl_index[i]+1 for i in ao_shell ] ezfio.set_ao_basis_ao_nucl(at) @@ -330,7 +327,7 @@ def write_ezfio(trexio_filename, filename): print("OK") else: - print("None") + print("None: integrals should be also imported using qp run import_trexio_integrals") # _ diff --git a/src/ao_one_e_ints/one_e_Coul_integrals_cosgtos.irp.f b/src/ao_one_e_ints/one_e_Coul_integrals_cosgtos.irp.f index 7f94f226..88fe4c85 100644 --- a/src/ao_one_e_ints/one_e_Coul_integrals_cosgtos.irp.f +++ b/src/ao_one_e_ints/one_e_Coul_integrals_cosgtos.irp.f @@ -111,8 +111,9 @@ complex*16 function NAI_pol_mult_cosgtos(A_center, B_center, power_A, power_B, a complex*16 :: accu, P_center(3) complex*16 :: d(0:n_pt_in) - complex*16 :: V_n_e_cosgtos - complex*16 :: crint + complex*16, external :: V_n_e_cosgtos + complex*16, external :: crint_2 + complex*16, external :: crint_sum_2 if ( (A_center(1)/=B_center(1)) .or. (A_center(2)/=B_center(2)) .or. (A_center(3)/=B_center(3)) .or. & (A_center(1)/=C_center(1)) .or. (A_center(2)/=C_center(2)) .or. (A_center(3)/=C_center(3)) ) then @@ -158,27 +159,27 @@ complex*16 function NAI_pol_mult_cosgtos(A_center, B_center, power_A, power_B, a n_pt = 2 * ( (power_A(1) + power_B(1)) + (power_A(2) + power_B(2)) + (power_A(3) + power_B(3)) ) if(n_pt == 0) then - NAI_pol_mult_cosgtos = coeff * crint(0, const) + NAI_pol_mult_cosgtos = coeff * crint_2(0, const) return endif - call give_cpolynomial_mult_center_one_e( A_center, B_center, alpha, beta & - , power_A, power_B, C_center, n_pt_in, d, n_pt_out) + call give_cpolynomial_mult_center_one_e(A_center, B_center, alpha, beta, & + power_A, power_B, C_center, n_pt_in, d, n_pt_out) if(n_pt_out < 0) then NAI_pol_mult_cosgtos = (0.d0, 0.d0) return endif - accu = (0.d0, 0.d0) - do i = 0, n_pt_out, 2 - accu += crint(shiftr(i, 1), const) * d(i) - -! print *, shiftr(i, 1), real(const), real(d(i)), real(crint(shiftr(i, 1), const)) - enddo + !accu = (0.d0, 0.d0) + !do i = 0, n_pt_out, 2 + ! accu += crint_2(shiftr(i, 1), const) * d(i) + !enddo + accu = crint_sum_2(n_pt_out, const, d) NAI_pol_mult_cosgtos = accu * coeff -end function NAI_pol_mult_cosgtos + return +end ! --- @@ -312,7 +313,7 @@ subroutine give_cpolynomial_mult_center_one_e( A_center, B_center, alpha, beta & d(i) = d1(i) enddo -end subroutine give_cpolynomial_mult_center_one_e +end ! --- @@ -405,7 +406,7 @@ recursive subroutine I_x1_pol_mult_one_e_cosgtos(a, c, R1x, R1xp, R2x, d, nd, n_ endif -end subroutine I_x1_pol_mult_one_e_cosgtos +end ! --- @@ -467,7 +468,7 @@ recursive subroutine I_x2_pol_mult_one_e_cosgtos(c, R1x, R1xp, R2x, d, nd, dim) endif -end subroutine I_x2_pol_mult_one_e_cosgtos +end ! --- @@ -502,7 +503,7 @@ complex*16 function V_n_e_cosgtos(a_x, a_y, a_z, b_x, b_y, b_z, alpha, beta) * V_theta(a_z + b_z, a_x + b_x + a_y + b_y + 1) endif -end function V_n_e_cosgtos +end ! --- @@ -529,7 +530,7 @@ complex*16 function V_r_cosgtos(n, alpha) V_r_cosgtos = sqpi * fact(n) / fact(shiftr(n, 1)) * (0.5d0/zsqrt(alpha))**(n+1) endif -end function V_r_cosgtos +end ! --- diff --git a/src/ao_two_e_ints/cholesky.irp.f b/src/ao_two_e_ints/cholesky.irp.f index ccaa7239..efafd504 100644 --- a/src/ao_two_e_ints/cholesky.irp.f +++ b/src/ao_two_e_ints/cholesky.irp.f @@ -5,7 +5,7 @@ double precision function get_ao_integ_chol(i,j,k,l) ! i(r1) j(r1) 1/r12 k(r2) l(r2) END_DOC integer, intent(in) :: i,j,k,l - double precision, external :: ddot + double precision, external :: ddot get_ao_integ_chol = ddot(cholesky_ao_num, cholesky_ao_transp(1,i,j), 1, cholesky_ao_transp(1,k,l), 1) end @@ -73,14 +73,13 @@ END_PROVIDER integer, external :: getUnitAndOpen integer :: iunit, ierr - ndim8 = ao_num*ao_num*1_8 + ndim8 = ao_num*ao_num*1_8+1 double precision :: wall0,wall1 - type(c_ptr) :: c_pointer(2) - integer :: fd(2) + type(mmap_type) :: map PROVIDE nproc ao_cholesky_threshold do_direct_integrals qp_max_mem - PROVIDE nucl_coord ao_two_e_integral_schwartz + PROVIDE nucl_coord call set_multiple_levels_omp(.False.) call wall_time(wall0) @@ -143,19 +142,21 @@ END_PROVIDER if (do_direct_integrals) then !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i8) SCHEDULE(dynamic,21) - do i8=ndim8,1,-1 + do i8=ndim8-1,1,-1 D(i8) = ao_two_e_integral(addr1(i8), addr2(i8), & addr1(i8), addr2(i8)) enddo !$OMP END PARALLEL DO else !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i8) SCHEDULE(dynamic,21) - do i8=ndim8,1,-1 + do i8=ndim8-1,1,-1 D(i8) = get_ao_two_e_integral(addr1(i8), addr1(i8), & addr2(i8), addr2(i8), ao_integrals_map) enddo !$OMP END PARALLEL DO endif + ! Just to guarentee termination + D(ndim8) = 0.d0 D_sorted(:) = -D(:) call dsort_noidx_big(D_sorted,ndim8) @@ -179,14 +180,9 @@ END_PROVIDER if (elec_num > 10) then rank_max = min(np,20*elec_num*elec_num) endif - call mmap(trim(ezfio_work_dir)//'cholesky_ao_tmp', (/ ndim8, rank_max /), 8, fd(1), .False., .True., c_pointer(1)) - call c_f_pointer(c_pointer(1), L, (/ ndim8, rank_max /)) - - ! Deleting the file while it is open makes the file invisible on the filesystem, - ! and automatically deleted, even if the program crashes - iunit = getUnitAndOpen(trim(ezfio_work_dir)//'cholesky_ao_tmp', 'R') - close(iunit,status='delete') + call mmap_create_d('', (/ ndim8, rank_max /), .False., .True., map) + L => map%d2 ! 3. N = 0 @@ -205,6 +201,7 @@ END_PROVIDER i = i+1 + block_size = max(N,24) ! Determine nq so that Delta fits in memory @@ -314,9 +311,10 @@ END_PROVIDER ! g. iblock = 0 + do j=1,nq - if ( (Qmax <= Dmin).or.(N+j*1_8 > ndim8) ) exit + if ( (Qmax < Dmin).or.(N+j*1_8 > ndim8) ) exit ! i. rank = N+j @@ -476,7 +474,7 @@ END_PROVIDER enddo !$OMP END PARALLEL DO - call munmap( (/ ndim8, rank_max /), 8, fd(1), c_pointer(1) ) + call mmap_destroy(map) cholesky_ao_num = rank diff --git a/src/ao_two_e_ints/two_e_Coul_integrals_cosgtos.irp.f b/src/ao_two_e_ints/two_e_Coul_integrals_cosgtos.irp.f index ea9ff009..df402ff1 100644 --- a/src/ao_two_e_ints/two_e_Coul_integrals_cosgtos.irp.f +++ b/src/ao_two_e_ints/two_e_Coul_integrals_cosgtos.irp.f @@ -35,11 +35,9 @@ double precision function ao_two_e_integral_cosgtos(i, j, k, l) if(ao_prim_num(i) * ao_prim_num(j) * ao_prim_num(k) * ao_prim_num(l) > 1024) then - !print *, ' with shwartz acc ' ao_two_e_integral_cosgtos = ao_2e_cosgtos_schwartz_accel(i, j, k, l) else - !print *, ' without shwartz acc ' dim1 = n_pt_max_integrals @@ -51,7 +49,6 @@ double precision function ao_two_e_integral_cosgtos(i, j, k, l) ao_two_e_integral_cosgtos = 0.d0 if(num_i /= num_j .or. num_k /= num_l .or. num_j /= num_k) then - !print *, ' not the same center' do p = 1, 3 I_power(p) = ao_power(i,p) @@ -72,72 +69,22 @@ double precision function ao_two_e_integral_cosgtos(i, j, k, l) coef2 = coef1 * ao_coef_norm_ord_transp_cosgtos(q,j) expo2 = ao_expo_ord_transp_cosgtos(q,j) - call give_explicit_cpoly_and_cgaussian( P1_new, P1_center, pp1, fact_p1, iorder_p1 & - , expo1, expo2, I_power, J_power, I_center, J_center, dim1 ) + call give_explicit_cpoly_and_cgaussian(P1_new, P1_center, pp1, fact_p1, iorder_p1, & + expo1, expo2, I_power, J_power, I_center, J_center, dim1) p1_inv = (1.d0,0.d0) / pp1 - call give_explicit_cpoly_and_cgaussian( P2_new, P2_center, pp2, fact_p2, iorder_p2 & - , conjg(expo1), expo2, I_power, J_power, I_center, J_center, dim1 ) + call give_explicit_cpoly_and_cgaussian(P2_new, P2_center, pp2, fact_p2, iorder_p2, & + conjg(expo1), expo2, I_power, J_power, I_center, J_center, dim1) p2_inv = (1.d0,0.d0) / pp2 - call give_explicit_cpoly_and_cgaussian( P3_new, P3_center, pp3, fact_p3, iorder_p3 & - , expo1, conjg(expo2), I_power, J_power, I_center, J_center, dim1 ) + call give_explicit_cpoly_and_cgaussian(P3_new, P3_center, pp3, fact_p3, iorder_p3, & + expo1, conjg(expo2), I_power, J_power, I_center, J_center, dim1) p3_inv = (1.d0,0.d0) / pp3 - call give_explicit_cpoly_and_cgaussian( P4_new, P4_center, pp4, fact_p4, iorder_p4 & - , conjg(expo1), conjg(expo2), I_power, J_power, I_center, J_center, dim1 ) + call give_explicit_cpoly_and_cgaussian(P4_new, P4_center, pp4, fact_p4, iorder_p4, & + conjg(expo1), conjg(expo2), I_power, J_power, I_center, J_center, dim1) p4_inv = (1.d0,0.d0) / pp4 - !integer :: ii - !do ii = 1, 3 - ! print *, 'fact_p1', fact_p1 - ! print *, 'fact_p2', fact_p2 - ! print *, 'fact_p3', fact_p3 - ! print *, 'fact_p4', fact_p4 - ! !print *, pp1, p1_inv - ! !print *, pp2, p2_inv - ! !print *, pp3, p3_inv - ! !print *, pp4, p4_inv - !enddo - ! if( abs(aimag(P1_center(ii))) .gt. 0.d0 ) then - ! print *, ' P_1 is complex !!' - ! print *, P1_center - ! print *, expo1, expo2 - ! print *, conjg(expo1), conjg(expo2) - ! stop - ! endif - ! if( abs(aimag(P2_center(ii))) .gt. 0.d0 ) then - ! print *, ' P_2 is complex !!' - ! print *, P2_center - ! print *, ' old expos:' - ! print *, expo1, expo2 - ! print *, conjg(expo1), conjg(expo2) - ! print *, ' new expo:' - ! print *, pp2, p2_inv - ! print *, ' factor:' - ! print *, fact_p2 - ! print *, ' old centers:' - ! print *, I_center, J_center - ! print *, ' powers:' - ! print *, I_power, J_power - ! stop - ! endif - ! if( abs(aimag(P3_center(ii))) .gt. 0.d0 ) then - ! print *, ' P_3 is complex !!' - ! print *, P3_center - ! print *, expo1, expo2 - ! print *, conjg(expo1), conjg(expo2) - ! stop - ! endif - ! if( abs(aimag(P4_center(ii))) .gt. 0.d0 ) then - ! print *, ' P_4 is complex !!' - ! print *, P4_center - ! print *, expo1, expo2 - ! print *, conjg(expo1), conjg(expo2) - ! stop - ! endif - !enddo - do r = 1, ao_prim_num(k) coef3 = coef2 * ao_coef_norm_ord_transp_cosgtos(r,k) expo3 = ao_expo_ord_transp_cosgtos(r,k) @@ -146,74 +93,47 @@ double precision function ao_two_e_integral_cosgtos(i, j, k, l) coef4 = coef3 * ao_coef_norm_ord_transp_cosgtos(s,l) expo4 = ao_expo_ord_transp_cosgtos(s,l) - call give_explicit_cpoly_and_cgaussian( Q1_new, Q1_center, qq1, fact_q1, iorder_q1 & - , expo3, expo4, K_power, L_power, K_center, L_center, dim1 ) + call give_explicit_cpoly_and_cgaussian(Q1_new, Q1_center, qq1, fact_q1, iorder_q1, & + expo3, expo4, K_power, L_power, K_center, L_center, dim1) q1_inv = (1.d0,0.d0) / qq1 - call give_explicit_cpoly_and_cgaussian( Q2_new, Q2_center, qq2, fact_q2, iorder_q2 & - , conjg(expo3), expo4, K_power, L_power, K_center, L_center, dim1 ) + call give_explicit_cpoly_and_cgaussian(Q2_new, Q2_center, qq2, fact_q2, iorder_q2, & + conjg(expo3), expo4, K_power, L_power, K_center, L_center, dim1) q2_inv = (1.d0,0.d0) / qq2 - !do ii = 1, 3 - ! !print *, qq1, q1_inv - ! !print *, qq2, q2_inv - ! print *, 'fact_q1', fact_q1 - ! print *, 'fact_q2', fact_q2 - !enddo - ! if( abs(aimag(Q1_center(ii))) .gt. 0.d0 ) then - ! print *, ' Q_1 is complex !!' - ! print *, Q1_center - ! print *, expo3, expo4 - ! print *, conjg(expo3), conjg(expo4) - ! stop - ! endif - ! if( abs(aimag(Q2_center(ii))) .gt. 0.d0 ) then - ! print *, ' Q_2 is complex !!' - ! print *, Q2_center - ! print *, expo3, expo4 - ! print *, conjg(expo3), conjg(expo4) - ! stop - ! endif - !enddo + integral1 = general_primitive_integral_cosgtos(dim1, P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1, & + Q1_new, Q1_center, fact_q1, qq1, q1_inv, iorder_q1) + integral2 = general_primitive_integral_cosgtos(dim1, P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1, & + Q2_new, Q2_center, fact_q2, qq2, q2_inv, iorder_q2) - integral1 = general_primitive_integral_cosgtos( dim1, P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1 & - , Q1_new, Q1_center, fact_q1, qq1, q1_inv, iorder_q1 ) + integral3 = general_primitive_integral_cosgtos(dim1, P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2, & + Q1_new, Q1_center, fact_q1, qq1, q1_inv, iorder_q1) - integral2 = general_primitive_integral_cosgtos( dim1, P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1 & - , Q2_new, Q2_center, fact_q2, qq2, q2_inv, iorder_q2 ) + integral4 = general_primitive_integral_cosgtos(dim1, P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2, & + Q2_new, Q2_center, fact_q2, qq2, q2_inv, iorder_q2) - integral3 = general_primitive_integral_cosgtos( dim1, P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2 & - , Q1_new, Q1_center, fact_q1, qq1, q1_inv, iorder_q1 ) + integral5 = general_primitive_integral_cosgtos(dim1, P3_new, P3_center, fact_p3, pp3, p3_inv, iorder_p3, & + Q1_new, Q1_center, fact_q1, qq1, q1_inv, iorder_q1) - integral4 = general_primitive_integral_cosgtos( dim1, P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2 & - , Q2_new, Q2_center, fact_q2, qq2, q2_inv, iorder_q2 ) + integral6 = general_primitive_integral_cosgtos(dim1, P3_new, P3_center, fact_p3, pp3, p3_inv, iorder_p3, & + Q2_new, Q2_center, fact_q2, qq2, q2_inv, iorder_q2) - integral5 = general_primitive_integral_cosgtos( dim1, P3_new, P3_center, fact_p3, pp3, p3_inv, iorder_p3 & - , Q1_new, Q1_center, fact_q1, qq1, q1_inv, iorder_q1 ) + integral7 = general_primitive_integral_cosgtos(dim1, P4_new, P4_center, fact_p4, pp4, p4_inv, iorder_p4, & + Q1_new, Q1_center, fact_q1, qq1, q1_inv, iorder_q1) - integral6 = general_primitive_integral_cosgtos( dim1, P3_new, P3_center, fact_p3, pp3, p3_inv, iorder_p3 & - , Q2_new, Q2_center, fact_q2, qq2, q2_inv, iorder_q2 ) - - integral7 = general_primitive_integral_cosgtos( dim1, P4_new, P4_center, fact_p4, pp4, p4_inv, iorder_p4 & - , Q1_new, Q1_center, fact_q1, qq1, q1_inv, iorder_q1 ) - - integral8 = general_primitive_integral_cosgtos( dim1, P4_new, P4_center, fact_p4, pp4, p4_inv, iorder_p4 & - , Q2_new, Q2_center, fact_q2, qq2, q2_inv, iorder_q2 ) + integral8 = general_primitive_integral_cosgtos(dim1, P4_new, P4_center, fact_p4, pp4, p4_inv, iorder_p4, & + Q2_new, Q2_center, fact_q2, qq2, q2_inv, iorder_q2) integral_tot = integral1 + integral2 + integral3 + integral4 + integral5 + integral6 + integral7 + integral8 - !integral_tot = integral1 - !print*, integral_tot - ao_two_e_integral_cosgtos = ao_two_e_integral_cosgtos + coef4 * 2.d0 * real(integral_tot) enddo ! s - enddo ! r - enddo ! q - enddo ! p + enddo ! r + enddo ! q + enddo ! p else - !print *, ' the same center' do p = 1, 3 I_power(p) = ao_power(i,p) @@ -290,7 +210,7 @@ double precision function ao_two_e_integral_cosgtos(i, j, k, l) endif endif -end function ao_two_e_integral_cosgtos +end ! --- @@ -326,8 +246,8 @@ double precision function ao_2e_cosgtos_schwartz_accel(i, j, k, l) double precision :: thr double precision :: schwartz_ij - complex*16 :: ERI_cosgtos - complex*16 :: general_primitive_integral_cosgtos + complex*16, external :: ERI_cosgtos + complex*16, external :: general_primitive_integral_cosgtos ao_2e_cosgtos_schwartz_accel = 0.d0 @@ -341,7 +261,7 @@ double precision function ao_2e_cosgtos_schwartz_accel(i, j, k, l) thr = ao_integrals_threshold*ao_integrals_threshold - allocate( schwartz_kl(0:ao_prim_num(l),0:ao_prim_num(k)) ) + allocate(schwartz_kl(0:ao_prim_num(l),0:ao_prim_num(k))) if(num_i /= num_j .or. num_k /= num_l .or. num_j /= num_k) then @@ -366,45 +286,45 @@ double precision function ao_2e_cosgtos_schwartz_accel(i, j, k, l) coef2 = coef1 * ao_coef_norm_ord_transp_cosgtos(s,l) * ao_coef_norm_ord_transp_cosgtos(s,l) expo2 = ao_expo_ord_transp_cosgtos(s,l) - call give_explicit_cpoly_and_cgaussian( P1_new, P1_center, pp1, fact_p1, iorder_p1 & - , expo1, expo2, K_power, L_power, K_center, L_center, dim1 ) + call give_explicit_cpoly_and_cgaussian(P1_new, P1_center, pp1, fact_p1, iorder_p1, & + expo1, expo2, K_power, L_power, K_center, L_center, dim1) p1_inv = (1.d0,0.d0) / pp1 - call give_explicit_cpoly_and_cgaussian( P2_new, P2_center, pp2, fact_p2, iorder_p2 & - , conjg(expo1), expo2, K_power, L_power, K_center, L_center, dim1 ) + call give_explicit_cpoly_and_cgaussian(P2_new, P2_center, pp2, fact_p2, iorder_p2, & + conjg(expo1), expo2, K_power, L_power, K_center, L_center, dim1) p2_inv = (1.d0,0.d0) / pp2 - call give_explicit_cpoly_and_cgaussian( P3_new, P3_center, pp3, fact_p3, iorder_p3 & - , expo1, conjg(expo2), K_power, L_power, K_center, L_center, dim1 ) + call give_explicit_cpoly_and_cgaussian(P3_new, P3_center, pp3, fact_p3, iorder_p3, & + expo1, conjg(expo2), K_power, L_power, K_center, L_center, dim1) p3_inv = (1.d0,0.d0) / pp3 - call give_explicit_cpoly_and_cgaussian( P4_new, P4_center, pp4, fact_p4, iorder_p4 & - , conjg(expo1), conjg(expo2), K_power, L_power, K_center, L_center, dim1 ) + call give_explicit_cpoly_and_cgaussian(P4_new, P4_center, pp4, fact_p4, iorder_p4, & + conjg(expo1), conjg(expo2), K_power, L_power, K_center, L_center, dim1) p4_inv = (1.d0,0.d0) / pp4 - integral1 = general_primitive_integral_cosgtos( dim1, P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1 & - , P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1 ) + integral1 = general_primitive_integral_cosgtos(dim1, P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1, & + P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1) - integral2 = general_primitive_integral_cosgtos( dim1, P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1 & - , P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2 ) + integral2 = general_primitive_integral_cosgtos(dim1, P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1, & + P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2) - integral3 = general_primitive_integral_cosgtos( dim1, P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2 & - , P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1 ) + integral3 = general_primitive_integral_cosgtos(dim1, P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2, & + P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1) - integral4 = general_primitive_integral_cosgtos( dim1, P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2 & - , P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2 ) + integral4 = general_primitive_integral_cosgtos(dim1, P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2, & + P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2) - integral5 = general_primitive_integral_cosgtos( dim1, P3_new, P3_center, fact_p3, pp3, p3_inv, iorder_p3 & - , P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1 ) + integral5 = general_primitive_integral_cosgtos(dim1, P3_new, P3_center, fact_p3, pp3, p3_inv, iorder_p3, & + P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1) - integral6 = general_primitive_integral_cosgtos( dim1, P3_new, P3_center, fact_p3, pp3, p3_inv, iorder_p3 & - , P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2 ) + integral6 = general_primitive_integral_cosgtos(dim1, P3_new, P3_center, fact_p3, pp3, p3_inv, iorder_p3, & + P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2) - integral7 = general_primitive_integral_cosgtos( dim1, P4_new, P4_center, fact_p4, pp4, p4_inv, iorder_p4 & - , P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1 ) + integral7 = general_primitive_integral_cosgtos(dim1, P4_new, P4_center, fact_p4, pp4, p4_inv, iorder_p4, & + P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1) - integral8 = general_primitive_integral_cosgtos( dim1, P4_new, P4_center, fact_p4, pp4, p4_inv, iorder_p4 & - , P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2 ) + integral8 = general_primitive_integral_cosgtos(dim1, P4_new, P4_center, fact_p4, pp4, p4_inv, iorder_p4, & + P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2) integral_tot = integral1 + integral2 + integral3 + integral4 + integral5 + integral6 + integral7 + integral8 @@ -544,41 +464,45 @@ double precision function ao_2e_cosgtos_schwartz_accel(i, j, k, l) coef2 = coef1 * ao_coef_norm_ord_transp_cosgtos(s,l) * ao_coef_norm_ord_transp_cosgtos(s,l) expo2 = ao_expo_ord_transp_cosgtos(s,l) - integral1 = ERI_cosgtos( expo1, expo2, expo1, expo2 & - , K_power(1), L_power(1), K_power(1), L_power(1) & - , K_power(2), L_power(2), K_power(2), L_power(2) & - , K_power(3), L_power(3), K_power(3), L_power(3) ) - integral2 = ERI_cosgtos( expo1, expo2, conjg(expo1), expo2 & - , K_power(1), L_power(1), K_power(1), L_power(1) & - , K_power(2), L_power(2), K_power(2), L_power(2) & - , K_power(3), L_power(3), K_power(3), L_power(3) ) + integral1 = ERI_cosgtos(expo1, expo2, expo1, expo2, & + K_power(1), L_power(1), K_power(1), L_power(1), & + K_power(2), L_power(2), K_power(2), L_power(2), & + K_power(3), L_power(3), K_power(3), L_power(3)) - integral3 = ERI_cosgtos( conjg(expo1), expo2, expo1, expo2 & - , K_power(1), L_power(1), K_power(1), L_power(1) & - , K_power(2), L_power(2), K_power(2), L_power(2) & - , K_power(3), L_power(3), K_power(3), L_power(3) ) - integral4 = ERI_cosgtos( conjg(expo1), expo2, conjg(expo1), expo2 & - , K_power(1), L_power(1), K_power(1), L_power(1) & - , K_power(2), L_power(2), K_power(2), L_power(2) & - , K_power(3), L_power(3), K_power(3), L_power(3) ) + integral2 = ERI_cosgtos(expo1, expo2, conjg(expo1), expo2, & + K_power(1), L_power(1), K_power(1), L_power(1), & + K_power(2), L_power(2), K_power(2), L_power(2), & + K_power(3), L_power(3), K_power(3), L_power(3)) - integral5 = ERI_cosgtos( expo1, conjg(expo2), expo1, expo2 & - , K_power(1), L_power(1), K_power(1), L_power(1) & - , K_power(2), L_power(2), K_power(2), L_power(2) & - , K_power(3), L_power(3), K_power(3), L_power(3) ) - integral6 = ERI_cosgtos( expo1, conjg(expo2), conjg(expo1), expo2 & - , K_power(1), L_power(1), K_power(1), L_power(1) & - , K_power(2), L_power(2), K_power(2), L_power(2) & - , K_power(3), L_power(3), K_power(3), L_power(3) ) + integral3 = ERI_cosgtos(conjg(expo1), expo2, expo1, expo2, & + K_power(1), L_power(1), K_power(1), L_power(1), & + K_power(2), L_power(2), K_power(2), L_power(2), & + K_power(3), L_power(3), K_power(3), L_power(3)) - integral7 = ERI_cosgtos( conjg(expo1), conjg(expo2), expo1, expo2 & - , K_power(1), L_power(1), K_power(1), L_power(1) & - , K_power(2), L_power(2), K_power(2), L_power(2) & - , K_power(3), L_power(3), K_power(3), L_power(3) ) - integral8 = ERI_cosgtos( conjg(expo1), conjg(expo2), conjg(expo1), expo2 & - , K_power(1), L_power(1), K_power(1), L_power(1) & - , K_power(2), L_power(2), K_power(2), L_power(2) & - , K_power(3), L_power(3), K_power(3), L_power(3) ) + integral4 = ERI_cosgtos(conjg(expo1), expo2, conjg(expo1), expo2, & + K_power(1), L_power(1), K_power(1), L_power(1), & + K_power(2), L_power(2), K_power(2), L_power(2), & + K_power(3), L_power(3), K_power(3), L_power(3)) + + integral5 = ERI_cosgtos(expo1, conjg(expo2), expo1, expo2, & + K_power(1), L_power(1), K_power(1), L_power(1), & + K_power(2), L_power(2), K_power(2), L_power(2), & + K_power(3), L_power(3), K_power(3), L_power(3)) + + integral6 = ERI_cosgtos(expo1, conjg(expo2), conjg(expo1), expo2, & + K_power(1), L_power(1), K_power(1), L_power(1), & + K_power(2), L_power(2), K_power(2), L_power(2), & + K_power(3), L_power(3), K_power(3), L_power(3)) + + integral7 = ERI_cosgtos(conjg(expo1), conjg(expo2), expo1, expo2, & + K_power(1), L_power(1), K_power(1), L_power(1), & + K_power(2), L_power(2), K_power(2), L_power(2), & + K_power(3), L_power(3), K_power(3), L_power(3)) + + integral8 = ERI_cosgtos(conjg(expo1), conjg(expo2), conjg(expo1), expo2, & + K_power(1), L_power(1), K_power(1), L_power(1), & + K_power(2), L_power(2), K_power(2), L_power(2), & + K_power(3), L_power(3), K_power(3), L_power(3)) integral_tot = integral1 + integral2 + integral3 + integral4 + integral5 + integral6 + integral7 + integral8 @@ -598,45 +522,45 @@ double precision function ao_2e_cosgtos_schwartz_accel(i, j, k, l) coef2 = coef1 * ao_coef_norm_ord_transp_cosgtos(q,j) expo2 = ao_expo_ord_transp_cosgtos(q,j) - integral1 = ERI_cosgtos( expo1, expo2, expo1, expo2 & - , I_power(1), J_power(1), I_power(1), J_power(1) & - , I_power(2), J_power(2), I_power(2), J_power(2) & - , I_power(3), J_power(3), I_power(3), J_power(3) ) + integral1 = ERI_cosgtos(expo1, expo2, expo1, expo2, & + I_power(1), J_power(1), I_power(1), J_power(1), & + I_power(2), J_power(2), I_power(2), J_power(2), & + I_power(3), J_power(3), I_power(3), J_power(3)) - integral2 = ERI_cosgtos( expo1, expo2, conjg(expo1), expo2 & - , I_power(1), J_power(1), I_power(1), J_power(1) & - , I_power(2), J_power(2), I_power(2), J_power(2) & - , I_power(3), J_power(3), I_power(3), J_power(3) ) + integral2 = ERI_cosgtos(expo1, expo2, conjg(expo1), expo2, & + I_power(1), J_power(1), I_power(1), J_power(1), & + I_power(2), J_power(2), I_power(2), J_power(2), & + I_power(3), J_power(3), I_power(3), J_power(3)) - integral3 = ERI_cosgtos( conjg(expo1), expo2, expo1, expo2 & - , I_power(1), J_power(1), I_power(1), J_power(1) & - , I_power(2), J_power(2), I_power(2), J_power(2) & - , I_power(3), J_power(3), I_power(3), J_power(3) ) + integral3 = ERI_cosgtos(conjg(expo1), expo2, expo1, expo2, & + I_power(1), J_power(1), I_power(1), J_power(1), & + I_power(2), J_power(2), I_power(2), J_power(2), & + I_power(3), J_power(3), I_power(3), J_power(3)) - integral4 = ERI_cosgtos( conjg(expo1), expo2, conjg(expo1), expo2 & - , I_power(1), J_power(1), I_power(1), J_power(1) & - , I_power(2), J_power(2), I_power(2), J_power(2) & - , I_power(3), J_power(3), I_power(3), J_power(3) ) + integral4 = ERI_cosgtos(conjg(expo1), expo2, conjg(expo1), expo2, & + I_power(1), J_power(1), I_power(1), J_power(1), & + I_power(2), J_power(2), I_power(2), J_power(2), & + I_power(3), J_power(3), I_power(3), J_power(3)) - integral5 = ERI_cosgtos( expo1, conjg(expo2), expo1, expo2 & - , I_power(1), J_power(1), I_power(1), J_power(1) & - , I_power(2), J_power(2), I_power(2), J_power(2) & - , I_power(3), J_power(3), I_power(3), J_power(3) ) + integral5 = ERI_cosgtos(expo1, conjg(expo2), expo1, expo2, & + I_power(1), J_power(1), I_power(1), J_power(1), & + I_power(2), J_power(2), I_power(2), J_power(2), & + I_power(3), J_power(3), I_power(3), J_power(3)) - integral6 = ERI_cosgtos( expo1, conjg(expo2), conjg(expo1), expo2 & - , I_power(1), J_power(1), I_power(1), J_power(1) & - , I_power(2), J_power(2), I_power(2), J_power(2) & - , I_power(3), J_power(3), I_power(3), J_power(3) ) + integral6 = ERI_cosgtos(expo1, conjg(expo2), conjg(expo1), expo2, & + I_power(1), J_power(1), I_power(1), J_power(1), & + I_power(2), J_power(2), I_power(2), J_power(2), & + I_power(3), J_power(3), I_power(3), J_power(3)) - integral7 = ERI_cosgtos( conjg(expo1), conjg(expo2), expo1, expo2 & - , I_power(1), J_power(1), I_power(1), J_power(1) & - , I_power(2), J_power(2), I_power(2), J_power(2) & - , I_power(3), J_power(3), I_power(3), J_power(3) ) + integral7 = ERI_cosgtos(conjg(expo1), conjg(expo2), expo1, expo2, & + I_power(1), J_power(1), I_power(1), J_power(1), & + I_power(2), J_power(2), I_power(2), J_power(2), & + I_power(3), J_power(3), I_power(3), J_power(3)) - integral8 = ERI_cosgtos( conjg(expo1), conjg(expo2), conjg(expo1), expo2 & - , I_power(1), J_power(1), I_power(1), J_power(1) & - , I_power(2), J_power(2), I_power(2), J_power(2) & - , I_power(3), J_power(3), I_power(3), J_power(3) ) + integral8 = ERI_cosgtos(conjg(expo1), conjg(expo2), conjg(expo1), expo2, & + I_power(1), J_power(1), I_power(1), J_power(1), & + I_power(2), J_power(2), I_power(2), J_power(2), & + I_power(3), J_power(3), I_power(3), J_power(3)) integral_tot = integral1 + integral2 + integral3 + integral4 + integral5 + integral6 + integral7 + integral8 @@ -655,45 +579,45 @@ double precision function ao_2e_cosgtos_schwartz_accel(i, j, k, l) coef4 = coef3 * ao_coef_norm_ord_transp_cosgtos(s,l) expo4 = ao_expo_ord_transp_cosgtos(s,l) - integral1 = ERI_cosgtos( expo1, expo2, expo3, expo4 & - , I_power(1), J_power(1), K_power(1), L_power(1) & - , I_power(2), J_power(2), K_power(2), L_power(2) & - , I_power(3), J_power(3), K_power(3), L_power(3) ) + integral1 = ERI_cosgtos(expo1, expo2, expo3, expo4, & + I_power(1), J_power(1), K_power(1), L_power(1), & + I_power(2), J_power(2), K_power(2), L_power(2), & + I_power(3), J_power(3), K_power(3), L_power(3)) - integral2 = ERI_cosgtos( expo1, expo2, conjg(expo3), expo4 & - , I_power(1), J_power(1), K_power(1), L_power(1) & - , I_power(2), J_power(2), K_power(2), L_power(2) & - , I_power(3), J_power(3), K_power(3), L_power(3) ) + integral2 = ERI_cosgtos(expo1, expo2, conjg(expo3), expo4, & + I_power(1), J_power(1), K_power(1), L_power(1), & + I_power(2), J_power(2), K_power(2), L_power(2), & + I_power(3), J_power(3), K_power(3), L_power(3)) - integral3 = ERI_cosgtos( conjg(expo1), expo2, expo3, expo4 & - , I_power(1), J_power(1), K_power(1), L_power(1) & - , I_power(2), J_power(2), K_power(2), L_power(2) & - , I_power(3), J_power(3), K_power(3), L_power(3) ) + integral3 = ERI_cosgtos(conjg(expo1), expo2, expo3, expo4, & + I_power(1), J_power(1), K_power(1), L_power(1), & + I_power(2), J_power(2), K_power(2), L_power(2), & + I_power(3), J_power(3), K_power(3), L_power(3)) - integral4 = ERI_cosgtos( conjg(expo1), expo2, conjg(expo3), expo4 & - , I_power(1), J_power(1), K_power(1), L_power(1) & - , I_power(2), J_power(2), K_power(2), L_power(2) & - , I_power(3), J_power(3), K_power(3), L_power(3) ) + integral4 = ERI_cosgtos(conjg(expo1), expo2, conjg(expo3), expo4, & + I_power(1), J_power(1), K_power(1), L_power(1), & + I_power(2), J_power(2), K_power(2), L_power(2), & + I_power(3), J_power(3), K_power(3), L_power(3)) - integral5 = ERI_cosgtos( expo1, conjg(expo2), expo3, expo4 & - , I_power(1), J_power(1), K_power(1), L_power(1) & - , I_power(2), J_power(2), K_power(2), L_power(2) & - , I_power(3), J_power(3), K_power(3), L_power(3) ) + integral5 = ERI_cosgtos(expo1, conjg(expo2), expo3, expo4, & + I_power(1), J_power(1), K_power(1), L_power(1), & + I_power(2), J_power(2), K_power(2), L_power(2), & + I_power(3), J_power(3), K_power(3), L_power(3)) - integral6 = ERI_cosgtos( expo1, conjg(expo2), conjg(expo3), expo4 & - , I_power(1), J_power(1), K_power(1), L_power(1) & - , I_power(2), J_power(2), K_power(2), L_power(2) & - , I_power(3), J_power(3), K_power(3), L_power(3) ) + integral6 = ERI_cosgtos(expo1, conjg(expo2), conjg(expo3), expo4, & + I_power(1), J_power(1), K_power(1), L_power(1), & + I_power(2), J_power(2), K_power(2), L_power(2), & + I_power(3), J_power(3), K_power(3), L_power(3)) - integral7 = ERI_cosgtos( conjg(expo1), conjg(expo2), expo3, expo4 & - , I_power(1), J_power(1), K_power(1), L_power(1) & - , I_power(2), J_power(2), K_power(2), L_power(2) & - , I_power(3), J_power(3), K_power(3), L_power(3) ) + integral7 = ERI_cosgtos(conjg(expo1), conjg(expo2), expo3, expo4, & + I_power(1), J_power(1), K_power(1), L_power(1), & + I_power(2), J_power(2), K_power(2), L_power(2), & + I_power(3), J_power(3), K_power(3), L_power(3)) - integral8 = ERI_cosgtos( conjg(expo1), conjg(expo2), conjg(expo3), expo4 & - , I_power(1), J_power(1), K_power(1), L_power(1) & - , I_power(2), J_power(2), K_power(2), L_power(2) & - , I_power(3), J_power(3), K_power(3), L_power(3) ) + integral8 = ERI_cosgtos(conjg(expo1), conjg(expo2), conjg(expo3), expo4, & + I_power(1), J_power(1), K_power(1), L_power(1), & + I_power(2), J_power(2), K_power(2), L_power(2), & + I_power(3), J_power(3), K_power(3), L_power(3)) integral_tot = integral1 + integral2 + integral3 + integral4 + integral5 + integral6 + integral7 + integral8 @@ -707,11 +631,11 @@ double precision function ao_2e_cosgtos_schwartz_accel(i, j, k, l) deallocate(schwartz_kl) -end function ao_2e_cosgtos_schwartz_accel +end ! --- -BEGIN_PROVIDER [ double precision, ao_2e_cosgtos_schwartz, (ao_num,ao_num)] +BEGIN_PROVIDER [double precision, ao_2e_cosgtos_schwartz, (ao_num, ao_num)] BEGIN_DOC ! Needed to compute Schwartz inequalities @@ -739,8 +663,8 @@ END_PROVIDER ! --- -complex*16 function general_primitive_integral_cosgtos( dim, P_new, P_center, fact_p, p, p_inv, iorder_p & - , Q_new, Q_center, fact_q, q, q_inv, iorder_q ) +complex*16 function general_primitive_integral_cosgtos(dim, P_new, P_center, fact_p, p, p_inv, iorder_p, & + Q_new, Q_center, fact_q, q, q_inv, iorder_q) BEGIN_DOC ! @@ -765,7 +689,7 @@ complex*16 function general_primitive_integral_cosgtos( dim, P_new, P_center, fa complex*16 :: dx(0:max_dim), Ix_pol(0:max_dim), dy(0:max_dim), Iy_pol(0:max_dim), dz(0:max_dim), Iz_pol(0:max_dim) complex*16 :: d1(0:max_dim), d_poly(0:max_dim) - complex*16 :: crint_sum + complex*16 :: crint_sum_2 !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: dx, Ix_pol, dy, Iy_pol, dz, Iz_pol @@ -912,13 +836,11 @@ complex*16 function general_primitive_integral_cosgtos( dim, P_new, P_center, fa !DIR$ FORCEINLINE call multiply_cpoly(d_poly, n_pt_tmp, Iz_pol, n_Iz, d1, n_pt_out) - accu = crint_sum(n_pt_out, const, d1) -! print *, n_pt_out, real(d1(0:n_pt_out)) -! print *, real(accu) + accu = crint_sum_2(n_pt_out, const, d1) general_primitive_integral_cosgtos = fact_p * fact_q * accu * pi_5_2 * p_inv * q_inv / sq_ppq -end function general_primitive_integral_cosgtos +end ! --- @@ -994,7 +916,7 @@ complex*16 function ERI_cosgtos(alpha, beta, delta, gama, a_x, b_x, c_x, d_x, a_ ERI_cosgtos = I_f * coeff -end function ERI_cosgtos +end ! --- @@ -1076,7 +998,7 @@ subroutine integrale_new_cosgtos(I_f, a_x, b_x, c_x, d_x, a_y, b_y, c_y, d_y, a_ I_f += gauleg_w(i, j) * t1(i) enddo -end subroutine integrale_new_cosgtos +end ! --- @@ -1123,7 +1045,7 @@ recursive subroutine I_x1_new_cosgtos(a, c, B_10, B_01, B_00, res, n_pt) endif -end subroutine I_x1_new_cosgtos +end ! --- @@ -1163,7 +1085,7 @@ recursive subroutine I_x2_new_cosgtos(c, B_10, B_01, B_00, res, n_pt) endif -end subroutine I_x2_new_cosgtos +end ! --- @@ -1234,7 +1156,7 @@ subroutine give_cpolynom_mult_center_x( P_center, Q_center, a_x, d_x, p, q, n_pt return endif -end subroutine give_cpolynom_mult_center_x +end ! --- @@ -1277,7 +1199,7 @@ subroutine I_x1_pol_mult_cosgtos(a, c, B_10, B_01, B_00, C_00, D_00, d, nd, n_pt endif -end subroutine I_x1_pol_mult_cosgtos +end ! --- @@ -1370,7 +1292,7 @@ recursive subroutine I_x1_pol_mult_recurs_cosgtos(a, c, B_10, B_01, B_00, C_00, !DIR$ FORCEINLINE call multiply_cpoly(Y, ny, C_00, 2, d, nd) -end subroutine I_x1_pol_mult_recurs_cosgtos +end ! --- @@ -1426,7 +1348,7 @@ recursive subroutine I_x1_pol_mult_a1_cosgtos(c,B_10,B_01,B_00,C_00,D_00,d,nd,n_ !DIR$ FORCEINLINE call multiply_cpoly(Y, ny, C_00, 2, d, nd) -end subroutine I_x1_pol_mult_a1_cosgtos +end ! --- @@ -1490,7 +1412,7 @@ recursive subroutine I_x1_pol_mult_a2_cosgtos(c, B_10, B_01, B_00, C_00, D_00, d !DIR$ FORCEINLINE call multiply_cpoly(Y, ny, C_00, 2, d, nd) -end subroutine I_x1_pol_mult_a2_cosgtos +end ! --- @@ -1575,7 +1497,7 @@ recursive subroutine I_x2_pol_mult_cosgtos(c, B_10, B_01, B_00, C_00, D_00, d, n end select -end subroutine I_x2_pol_mult_cosgtos +end ! --- diff --git a/src/ao_two_e_ints/two_e_integrals.irp.f b/src/ao_two_e_ints/two_e_integrals.irp.f index d12f3d45..7b4a2e5a 100644 --- a/src/ao_two_e_ints/two_e_integrals.irp.f +++ b/src/ao_two_e_ints/two_e_integrals.irp.f @@ -1,6 +1,22 @@ ! --- +logical function do_schwartz_accel(i,j,k,l) + implicit none + BEGIN_DOC + ! If true, use Schwatrz to accelerate direct integral calculation + END_DOC + integer, intent(in) :: i, j, k, l + if (do_ao_cholesky) then + do_schwartz_accel = .False. + else + do_schwartz_accel = (ao_prim_num(i) * ao_prim_num(j) * & + ao_prim_num(k) * ao_prim_num(l) > 1024 ) + endif + +end function + + double precision function ao_two_e_integral(i, j, k, l) BEGIN_DOC @@ -25,6 +41,7 @@ double precision function ao_two_e_integral(i, j, k, l) double precision, external :: ao_two_e_integral_cosgtos double precision, external :: ao_two_e_integral_schwartz_accel + logical, external :: do_schwartz_accel if(use_cosgtos) then !print *, ' use_cosgtos for ao_two_e_integral ?', use_cosgtos @@ -35,7 +52,7 @@ double precision function ao_two_e_integral(i, j, k, l) ao_two_e_integral = ao_two_e_integral_erf(i, j, k, l) - else if (ao_prim_num(i) * ao_prim_num(j) * ao_prim_num(k) * ao_prim_num(l) > 1024 ) then + else if (do_schwartz_accel(i,j,k,l)) then ao_two_e_integral = ao_two_e_integral_schwartz_accel(i,j,k,l) diff --git a/src/ccsd/ccsd_space_orb_sub.irp.f b/src/ccsd/ccsd_space_orb_sub.irp.f index d8131a9c..e4907f22 100644 --- a/src/ccsd/ccsd_space_orb_sub.irp.f +++ b/src/ccsd/ccsd_space_orb_sub.irp.f @@ -154,14 +154,14 @@ subroutine run_ccsd_space_orb allocate(all_err(nO*nV+nO*nO*nV*(nV*1_8),cc_diis_depth), all_t(nO*nV+nO*nO*nV*(nV*1_8),cc_diis_depth)) !$OMP PARALLEL PRIVATE(i,j) DEFAULT(SHARED) + !$OMP DO COLLAPSE(2) do j=1,cc_diis_depth - !$OMP DO do i=1, size(all_err,1) all_err(i,j) = 0d0 all_t(i,j) = 0d0 enddo - !$OMP END DO NOWAIT enddo + !$OMP END DO NOWAIT !$OMP END PARALLEL endif @@ -237,6 +237,7 @@ subroutine run_ccsd_space_orb call update_t2(nO,nV,cc_space_f_o,cc_space_f_v,r2%f,t2%f) else print*,'Unkown cc_method_method: '//cc_update_method + call abort endif call update_tau_space(nO,nV,t1%f,t1,t2,tau) diff --git a/src/davidson/diagonalization_hs2_dressed.irp.f b/src/davidson/diagonalization_hs2_dressed.irp.f index 191e0021..d299f982 100644 --- a/src/davidson/diagonalization_hs2_dressed.irp.f +++ b/src/davidson/diagonalization_hs2_dressed.irp.f @@ -223,12 +223,11 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_ exit endif - if (itermax > 4) then - itermax = itermax - 1 - else if (m==1.and.disk_based_davidson) then + if (disk_based_davidson) then m=0 disk_based = .True. - itermax = 6 + else if (itermax > 4) then + itermax = itermax - 1 else nproc_target = nproc_target - 1 endif @@ -267,14 +266,12 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_ if (disk_based) then ! Create memory-mapped files for W and S - type(c_ptr) :: ptr_w, ptr_s - integer :: fd_s, fd_w - call mmap(trim(ezfio_work_dir)//'davidson_w', (/int(sze,8),int(N_st_diag*itermax,8)/),& - 8, fd_w, .False., .True., ptr_w) - call mmap(trim(ezfio_work_dir)//'davidson_s', (/int(sze,8),int(N_st_diag*itermax,8)/),& - 4, fd_s, .False., .True., ptr_s) - call c_f_pointer(ptr_w, w, (/sze,N_st_diag*itermax/)) - call c_f_pointer(ptr_s, s, (/sze,N_st_diag*itermax/)) + type(mmap_type) :: map_s, map_w + + call mmap_create_d('', (/ 1_8*sze, 1_8*N_st_diag*itermax /), .False., .True., map_w) + call mmap_create_s('', (/ 1_8*sze, 1_8*N_st_diag*itermax /), .False., .True., map_s) + w => map_w%d2 + s => map_s%s2 else allocate(W(sze,N_st_diag*itermax), S(sze,N_st_diag*itermax)) endif @@ -755,13 +752,8 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_ if (disk_based)then ! Remove temp files - integer, external :: getUnitAndOpen - call munmap( (/int(sze,8),int(N_st_diag*itermax,8)/), 8, fd_w, ptr_w ) - fd_w = getUnitAndOpen(trim(ezfio_work_dir)//'davidson_w','r') - close(fd_w,status='delete') - call munmap( (/int(sze,8),int(N_st_diag*itermax,8)/), 4, fd_s, ptr_s ) - fd_s = getUnitAndOpen(trim(ezfio_work_dir)//'davidson_s','r') - close(fd_s,status='delete') + call mmap_destroy(map_w) + call mmap_destroy(map_s) else deallocate(W,S) endif @@ -774,6 +766,7 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_ lambda & ) FREE nthreads_davidson + end diff --git a/src/davidson/diagonalize_ci.irp.f b/src/davidson/diagonalize_ci.irp.f index 6b852905..764309ed 100644 --- a/src/davidson/diagonalize_ci.irp.f +++ b/src/davidson/diagonalize_ci.irp.f @@ -330,6 +330,10 @@ END_PROVIDER deallocate(eigenvectors,eigenvalues) endif +! ! Dominant determinants for each states +! call print_dominant_det(psi_det,CI_eigenvectors,N_det,N_states,N_int) +! call wf_overlap(psi_det,psi_coef,N_states,N_det,psi_det,CI_eigenvectors,N_states,N_det) + END_PROVIDER subroutine diagonalize_CI diff --git a/src/davidson/u0_h_u0.irp.f b/src/davidson/u0_h_u0.irp.f index 7ef154a3..808bbb5d 100644 --- a/src/davidson/u0_h_u0.irp.f +++ b/src/davidson/u0_h_u0.irp.f @@ -179,10 +179,12 @@ subroutine H_u_0_nstates_openmp_work_$N_int(v_t,u_t,N_st,sze,istart,iend,ishift, ! ! compute_singles = (mem+rss > qp_max_mem) ! -! if (.not.compute_singles) then -! provide singles_beta_csc -! endif -compute_singles=.True. + compute_singles=.True. + + if (.not.compute_singles) then + provide singles_alpha_csc singles_beta_csc + endif + maxab = max(N_det_alpha_unique, N_det_beta_unique)+1 @@ -287,8 +289,7 @@ compute_singles=.True. tmp_det2(1:$N_int,2) = psi_det_beta_unique(1:$N_int, lcol) -!--- -! if (compute_singles) then + if (compute_singles) then l_a = psi_bilinear_matrix_columns_loc(lcol) ASSERT (l_a <= N_det) @@ -311,69 +312,67 @@ compute_singles=.True. buffer, idx, tmp_det(1,1), j, & singles_a, n_singles_a ) -!----- -! else -! -! ! Search for singles -! -!call cpu_time(time0) -! ! Right boundary -! l_a = psi_bilinear_matrix_columns_loc(lcol+1)-1 -! ASSERT (l_a <= N_det) -! do j=1,psi_bilinear_matrix_columns_loc(lcol+1) - psi_bilinear_matrix_columns_loc(lcol) -! lrow = psi_bilinear_matrix_rows(l_a) -! ASSERT (lrow <= N_det_alpha_unique) -! -! left = singles_alpha_csc_idx(krow) -! right_max = -1_8 -! right = singles_alpha_csc_idx(krow+1) -! do while (right-left>0_8) -! k8 = shiftr(right+left,1) -! if (singles_alpha_csc(k8) > lrow) then -! right = k8 -! else if (singles_alpha_csc(k8) < lrow) then -! left = k8 + 1_8 -! else -! right_max = k8+1_8 -! exit -! endif -! enddo -! if (right_max > 0_8) exit -! l_a = l_a-1 -! enddo -! if (right_max < 0_8) right_max = singles_alpha_csc_idx(krow) -! -! ! Search -! n_singles_a = 0 -! l_a = psi_bilinear_matrix_columns_loc(lcol) -! ASSERT (l_a <= N_det) -! -! last_found = singles_alpha_csc_idx(krow) -! do j=1,psi_bilinear_matrix_columns_loc(lcol+1) - psi_bilinear_matrix_columns_loc(lcol) -! lrow = psi_bilinear_matrix_rows(l_a) -! ASSERT (lrow <= N_det_alpha_unique) -! -! left = last_found -! right = right_max -! do while (right-left>0_8) -! k8 = shiftr(right+left,1) -! if (singles_alpha_csc(k8) > lrow) then -! right = k8 -! else if (singles_alpha_csc(k8) < lrow) then -! left = k8 + 1_8 -! else -! n_singles_a += 1 -! singles_a(n_singles_a) = l_a -! last_found = k8+1_8 -! exit -! endif -! enddo -! l_a = l_a+1 -! enddo -! j = j-1 -! -! endif -!----- + else + + ! Search for singles + + ! Right boundary + l_a = psi_bilinear_matrix_columns_loc(lcol+1)-1 + ASSERT (l_a <= N_det) + do j=1,psi_bilinear_matrix_columns_loc(lcol+1) - psi_bilinear_matrix_columns_loc(lcol) + lrow = psi_bilinear_matrix_rows(l_a) + ASSERT (lrow <= N_det_alpha_unique) + + left = singles_alpha_csc_idx(krow) + right_max = -1_8 + right = singles_alpha_csc_idx(krow+1) + do while (right-left>0_8) + k8 = shiftr(right+left,1) + if (singles_alpha_csc(k8) > lrow) then + right = k8 + else if (singles_alpha_csc(k8) < lrow) then + left = k8 + 1_8 + else + right_max = k8+1_8 + exit + endif + enddo + if (right_max > 0_8) exit + l_a = l_a-1 + enddo + if (right_max < 0_8) right_max = singles_alpha_csc_idx(krow) + + ! Search + n_singles_a = 0 + l_a = psi_bilinear_matrix_columns_loc(lcol) + ASSERT (l_a <= N_det) + + last_found = singles_alpha_csc_idx(krow) + do j=1,psi_bilinear_matrix_columns_loc(lcol+1) - psi_bilinear_matrix_columns_loc(lcol) + lrow = psi_bilinear_matrix_rows(l_a) + ASSERT (lrow <= N_det_alpha_unique) + + left = last_found + right = right_max + do while (right-left>0_8) + k8 = shiftr(right+left,1) + if (singles_alpha_csc(k8) > lrow) then + right = k8 + else if (singles_alpha_csc(k8) < lrow) then + left = k8 + 1_8 + else + n_singles_a += 1 + singles_a(n_singles_a) = l_a + last_found = k8+1_8 + exit + endif + enddo + l_a = l_a+1 + enddo + j = j-1 + + endif + ! Loop over alpha singles ! ----------------------- diff --git a/src/davidson/u0_hs2_u0.irp.f b/src/davidson/u0_hs2_u0.irp.f index 38fb56bd..3afe4ec6 100644 --- a/src/davidson/u0_hs2_u0.irp.f +++ b/src/davidson/u0_hs2_u0.irp.f @@ -218,10 +218,13 @@ subroutine H_S2_u_0_nstates_openmp_work_$N_int(v_t,s_t,u_t,N_st,sze,istart,iend, ! ! compute_singles = (mem+rss > qp_max_mem) ! -! if (.not.compute_singles) then -! provide singles_beta_csc -! endif -compute_singles=.True. + compute_singles=.True. + + if (.not.compute_singles) then + provide singles_alpha_csc + provide singles_beta_csc + endif + maxab = max(N_det_alpha_unique, N_det_beta_unique)+1 @@ -314,6 +317,7 @@ compute_singles=.True. singles_b(n_singles_b) = singles_beta_csc(k8) enddo endif + endif kcol_prev = kcol @@ -326,8 +330,7 @@ compute_singles=.True. tmp_det2(1:$N_int,2) = psi_det_beta_unique(1:$N_int, lcol) -!--- -! if (compute_singles) then + if (compute_singles) then l_a = psi_bilinear_matrix_columns_loc(lcol) ASSERT (l_a <= N_det) @@ -352,69 +355,66 @@ compute_singles=.True. buffer, idx, tmp_det(1,1), j, & singles_a, n_singles_a ) -!----- -! else -! -! ! Search for singles -! -!call cpu_time(time0) -! ! Right boundary -! l_a = psi_bilinear_matrix_columns_loc(lcol+1)-1 -! ASSERT (l_a <= N_det) -! do j=1,psi_bilinear_matrix_columns_loc(lcol+1) - psi_bilinear_matrix_columns_loc(lcol) -! lrow = psi_bilinear_matrix_rows(l_a) -! ASSERT (lrow <= N_det_alpha_unique) -! -! left = singles_alpha_csc_idx(krow) -! right_max = -1_8 -! right = singles_alpha_csc_idx(krow+1) -! do while (right-left>0_8) -! k8 = shiftr(right+left,1) -! if (singles_alpha_csc(k8) > lrow) then -! right = k8 -! else if (singles_alpha_csc(k8) < lrow) then -! left = k8 + 1_8 -! else -! right_max = k8+1_8 -! exit -! endif -! enddo -! if (right_max > 0_8) exit -! l_a = l_a-1 -! enddo -! if (right_max < 0_8) right_max = singles_alpha_csc_idx(krow) -! -! ! Search -! n_singles_a = 0 -! l_a = psi_bilinear_matrix_columns_loc(lcol) -! ASSERT (l_a <= N_det) -! -! last_found = singles_alpha_csc_idx(krow) -! do j=1,psi_bilinear_matrix_columns_loc(lcol+1) - psi_bilinear_matrix_columns_loc(lcol) -! lrow = psi_bilinear_matrix_rows(l_a) -! ASSERT (lrow <= N_det_alpha_unique) -! -! left = last_found -! right = right_max -! do while (right-left>0_8) -! k8 = shiftr(right+left,1) -! if (singles_alpha_csc(k8) > lrow) then -! right = k8 -! else if (singles_alpha_csc(k8) < lrow) then -! left = k8 + 1_8 -! else -! n_singles_a += 1 -! singles_a(n_singles_a) = l_a -! last_found = k8+1_8 -! exit -! endif -! enddo -! l_a = l_a+1 -! enddo -! j = j-1 -! -! endif -!----- + else + + ! Search for singles + + ! Right boundary + l_a = psi_bilinear_matrix_columns_loc(lcol+1)-1 + ASSERT (l_a <= N_det) + do j=1,psi_bilinear_matrix_columns_loc(lcol+1) - psi_bilinear_matrix_columns_loc(lcol) + lrow = psi_bilinear_matrix_rows(l_a) + ASSERT (lrow <= N_det_alpha_unique) + + left = singles_alpha_csc_idx(krow) + right_max = -1_8 + right = singles_alpha_csc_idx(krow+1) + do while (right-left>0_8) + k8 = shiftr(right+left,1) + if (singles_alpha_csc(k8) > lrow) then + right = k8 + else if (singles_alpha_csc(k8) < lrow) then + left = k8 + 1_8 + else + right_max = k8+1_8 + exit + endif + enddo + if (right_max > 0_8) exit + l_a = l_a-1 + enddo + if (right_max < 0_8) right_max = singles_alpha_csc_idx(krow) + + ! Search + n_singles_a = 0 + l_a = psi_bilinear_matrix_columns_loc(lcol) + ASSERT (l_a <= N_det) + + last_found = singles_alpha_csc_idx(krow) + do j=1,psi_bilinear_matrix_columns_loc(lcol+1) - psi_bilinear_matrix_columns_loc(lcol) + lrow = psi_bilinear_matrix_rows(l_a) + ASSERT (lrow <= N_det_alpha_unique) + + left = last_found + right = right_max + do while (right-left>0_8) + k8 = shiftr(right+left,1) + if (singles_alpha_csc(k8) > lrow) then + right = k8 + else if (singles_alpha_csc(k8) < lrow) then + left = k8 + 1_8 + else + n_singles_a += 1 + singles_a(n_singles_a) = l_a + last_found = k8+1_8 + exit + endif + enddo + l_a = l_a+1 + enddo + j = j-1 + + endif ! Loop over alpha singles ! ----------------------- diff --git a/src/davidson_keywords/EZFIO.cfg b/src/davidson_keywords/EZFIO.cfg index 6337b96f..81b7f949 100644 --- a/src/davidson_keywords/EZFIO.cfg +++ b/src/davidson_keywords/EZFIO.cfg @@ -48,7 +48,7 @@ default: false [distributed_davidson] type: logical -doc: If |true|, use the distributed algorithm -default: True +doc: If |true|, use the distributed algorithm. If you plan to run multi-node calculations, set this to true before running. +default: False interface: ezfio,provider,ocaml diff --git a/src/determinants/generate_cas_space.irp.f b/src/determinants/generate_cas_space.irp.f index 47a2ca30..05201c74 100644 --- a/src/determinants/generate_cas_space.irp.f +++ b/src/determinants/generate_cas_space.irp.f @@ -33,7 +33,7 @@ subroutine generate_cas_space print *, 'CAS(', n_alpha_act+n_beta_act, ', ', n_act_orb, ')' print *, '' - n_det_alpha_unique = binom_int(n_act_orb, n_alpha_act) + n_det_alpha_unique = int(binom_int(n_act_orb, n_alpha_act),4) TOUCH n_det_alpha_unique n = n_alpha_act @@ -56,7 +56,7 @@ subroutine generate_cas_space u = ior(t1,t2) enddo - n_det_beta_unique = binom_int(n_act_orb, n_beta_act) + n_det_beta_unique = int(binom_int(n_act_orb, n_beta_act),4) TOUCH n_det_beta_unique n = n_beta_act diff --git a/src/determinants/ref_bitmask.irp.f b/src/determinants/ref_bitmask.irp.f index 4e029ceb..18fa2396 100644 --- a/src/determinants/ref_bitmask.irp.f +++ b/src/determinants/ref_bitmask.irp.f @@ -30,31 +30,30 @@ ref_bitmask_energy += mo_one_e_integrals(occ(i,1),occ(i,1)) + mo_one_e_integrals(occ(i,2),occ(i,2)) ref_bitmask_kinetic_energy += mo_kinetic_integrals(occ(i,1),occ(i,1)) + mo_kinetic_integrals(occ(i,2),occ(i,2)) ref_bitmask_n_e_energy += mo_integrals_n_e(occ(i,1),occ(i,1)) + mo_integrals_n_e(occ(i,2),occ(i,2)) + do j = i+1, elec_alpha_num + ref_bitmask_two_e_energy += mo_two_e_integrals_jj_anti(occ(j,1),occ(i,1)) + ref_bitmask_energy += mo_two_e_integrals_jj_anti(occ(j,1),occ(i,1)) + enddo + do j= 1, elec_alpha_num + ref_bitmask_two_e_energy += mo_two_e_integrals_jj(occ(j,1),occ(i,2)) + ref_bitmask_energy += mo_two_e_integrals_jj(occ(j,1),occ(i,2)) + enddo + do j = i+1, elec_beta_num + ref_bitmask_two_e_energy += mo_two_e_integrals_jj_anti(occ(j,2),occ(i,2)) + ref_bitmask_energy += mo_two_e_integrals_jj_anti(occ(j,2),occ(i,2)) + enddo enddo do i = elec_beta_num+1,elec_alpha_num ref_bitmask_energy += mo_one_e_integrals(occ(i,1),occ(i,1)) ref_bitmask_kinetic_energy += mo_kinetic_integrals(occ(i,1),occ(i,1)) ref_bitmask_n_e_energy += mo_integrals_n_e(occ(i,1),occ(i,1)) - enddo - - do j= 1, elec_alpha_num - do i = j+1, elec_alpha_num - ref_bitmask_two_e_energy += mo_two_e_integrals_jj_anti(occ(i,1),occ(j,1)) - ref_bitmask_energy += mo_two_e_integrals_jj_anti(occ(i,1),occ(j,1)) + do j = i+1, elec_alpha_num + ref_bitmask_two_e_energy += mo_two_e_integrals_jj_anti(occ(j,1),occ(i,1)) + ref_bitmask_energy += mo_two_e_integrals_jj_anti(occ(j,1),occ(i,1)) enddo enddo - do j= 1, elec_beta_num - do i = j+1, elec_beta_num - ref_bitmask_two_e_energy += mo_two_e_integrals_jj_anti(occ(i,2),occ(j,2)) - ref_bitmask_energy += mo_two_e_integrals_jj_anti(occ(i,2),occ(j,2)) - enddo - do i= 1, elec_alpha_num - ref_bitmask_two_e_energy += mo_two_e_integrals_jj(occ(i,1),occ(j,2)) - ref_bitmask_energy += mo_two_e_integrals_jj(occ(i,1),occ(j,2)) - enddo - enddo ref_bitmask_one_e_energy = ref_bitmask_kinetic_energy + ref_bitmask_n_e_energy ref_bitmask_energy_ab = 0.d0 diff --git a/src/determinants/spindeterminants.irp.f b/src/determinants/spindeterminants.irp.f index dd55e112..87c5d360 100644 --- a/src/determinants/spindeterminants.irp.f +++ b/src/determinants/spindeterminants.irp.f @@ -910,6 +910,8 @@ subroutine copy_psi_bilinear_to_psi(psi, isize) end +use mmap_module + BEGIN_PROVIDER [ integer*8, singles_alpha_csc_idx, (N_det_alpha_unique+1) ] &BEGIN_PROVIDER [ integer*8, singles_alpha_csc_size ] implicit none @@ -925,12 +927,11 @@ end idx0(i) = i enddo - !$OMP PARALLEL DEFAULT(NONE) & - !$OMP SHARED(N_det_alpha_unique, psi_det_alpha_unique, & - !$OMP idx0, N_int, singles_alpha_csc, & - !$OMP elec_alpha_num, mo_num, singles_alpha_csc_idx) & + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP SHARED(N_det_alpha_unique, psi_det_alpha_unique, & + !$OMP idx0, N_int, singles_alpha_csc_idx) & !$OMP PRIVATE(i,s,j) - allocate (s(elec_alpha_num * (mo_num-elec_alpha_num) )) + allocate (s(N_det_alpha_unique)) !$OMP DO SCHEDULE(static,64) do i=1, N_det_alpha_unique call get_all_spin_singles( & @@ -966,7 +967,7 @@ BEGIN_PROVIDER [ integer, singles_alpha_csc, (singles_alpha_csc_size) ] !$OMP PARALLEL DO DEFAULT(NONE) & !$OMP SHARED(N_det_alpha_unique, psi_det_alpha_unique, & !$OMP idx0, N_int, singles_alpha_csc, singles_alpha_csc_idx)& - !$OMP PRIVATE(i,k) SCHEDULE(static,1) + !$OMP PRIVATE(i,k) SCHEDULE(static) do i=1, N_det_alpha_unique call get_all_spin_singles( & psi_det_alpha_unique, idx0, psi_det_alpha_unique(1,i), N_int,& @@ -978,7 +979,36 @@ BEGIN_PROVIDER [ integer, singles_alpha_csc, (singles_alpha_csc_size) ] END_PROVIDER +BEGIN_PROVIDER [ type(mmap_type), singles_alpha_csc_map ] + implicit none + BEGIN_DOC + ! Indices of all single excitations + END_DOC + integer :: i, k + integer, allocatable :: idx0(:) + call mmap_create_i('', (/ 1_8*singles_alpha_csc_size /), & + .False., .False., singles_alpha_csc_map) + + allocate (idx0(N_det_alpha_unique)) + do i=1, N_det_alpha_unique + idx0(i) = i + enddo + + !$OMP PARALLEL DO DEFAULT(NONE) & + !$OMP SHARED(N_det_alpha_unique, psi_det_alpha_unique, & + !$OMP idx0, N_int, singles_alpha_csc_map, singles_alpha_csc_idx)& + !$OMP PRIVATE(i,k) SCHEDULE(static) + do i=1, N_det_alpha_unique + call get_all_spin_singles( & + psi_det_alpha_unique, idx0, psi_det_alpha_unique(1,i), N_int, N_det_alpha_unique, & + singles_alpha_csc_map%i1(singles_alpha_csc_idx(i):singles_alpha_csc_idx(i)+N_det_alpha_unique-1),& + k) + enddo + !$OMP END PARALLEL DO + deallocate(idx0) + +END_PROVIDER BEGIN_PROVIDER [ integer*8, singles_beta_csc_idx, (N_det_beta_unique+1) ] @@ -996,13 +1026,12 @@ END_PROVIDER idx0(i) = i enddo - !$OMP PARALLEL DEFAULT(NONE) & + !$OMP PARALLEL DEFAULT(NONE) & !$OMP SHARED(N_det_beta_unique, psi_det_beta_unique, & - !$OMP idx0, N_int, singles_beta_csc, & - !$OMP elec_beta_num, mo_num, singles_beta_csc_idx) & + !$OMP idx0, N_int, singles_beta_csc_idx) & !$OMP PRIVATE(i,s,j) - allocate (s(elec_beta_num*(mo_num-elec_beta_num))) - !$OMP DO SCHEDULE(static,1) + allocate (s(N_det_beta_unique)) + !$OMP DO SCHEDULE(static) do i=1, N_det_beta_unique call get_all_spin_singles( & psi_det_beta_unique, idx0, psi_det_beta_unique(1,i), N_int,& @@ -1037,7 +1066,7 @@ BEGIN_PROVIDER [ integer, singles_beta_csc, (singles_beta_csc_size) ] !$OMP PARALLEL DO DEFAULT(NONE) & !$OMP SHARED(N_det_beta_unique, psi_det_beta_unique, & !$OMP idx0, N_int, singles_beta_csc, singles_beta_csc_idx)& - !$OMP PRIVATE(i,k) SCHEDULE(static,64) + !$OMP PRIVATE(i,k) SCHEDULE(static) do i=1, N_det_beta_unique call get_all_spin_singles( & psi_det_beta_unique, idx0, psi_det_beta_unique(1,i), N_int,& @@ -1049,6 +1078,37 @@ BEGIN_PROVIDER [ integer, singles_beta_csc, (singles_beta_csc_size) ] END_PROVIDER +BEGIN_PROVIDER [ type(mmap_type), singles_beta_csc_map ] + implicit none + BEGIN_DOC + ! Indices of all single excitations + END_DOC + integer :: i, k + integer, allocatable :: idx0(:) + + call mmap_create_i('', (/ 1_8*singles_beta_csc_size /), & + .False., .False., singles_beta_csc_map) + + allocate (idx0(N_det_beta_unique)) + do i=1, N_det_beta_unique + idx0(i) = i + enddo + + !$OMP PARALLEL DO DEFAULT(NONE) & + !$OMP SHARED(N_det_beta_unique, psi_det_beta_unique, & + !$OMP idx0, N_int, singles_beta_csc_map, singles_beta_csc_idx)& + !$OMP PRIVATE(i,k) SCHEDULE(static) + do i=1, N_det_beta_unique + call get_all_spin_singles( & + psi_det_beta_unique, idx0, psi_det_beta_unique(1,i), N_int, N_det_beta_unique, & + singles_beta_csc_map%i1(singles_beta_csc_idx(i):singles_beta_csc_idx(i)+N_det_beta_unique-1),& + k) + enddo + !$OMP END PARALLEL DO + deallocate(idx0) + +END_PROVIDER + @@ -1111,16 +1171,16 @@ subroutine get_all_spin_singles_1(buffer, idx, spindet, size_buffer, singles, n_ integer :: i integer(bit_kind) :: v integer :: degree - integer :: add_single(0:64) = (/ 0, 0, 1, 0, 0, (0, i=1,60) /) include 'utils/constants.include.F' - n_singles = 1 + n_singles = 0 do i=1,size_buffer degree = popcnt(xor( spindet, buffer(i) )) - singles(n_singles) = idx(i) - n_singles = n_singles+add_single(degree) + if (degree == 2) then + n_singles = n_singles+1 + singles(n_singles) = idx(i) + endif enddo - n_singles = n_singles-1 end @@ -1142,15 +1202,15 @@ subroutine get_all_spin_doubles_1(buffer, idx, spindet, size_buffer, doubles, n_ integer :: i include 'utils/constants.include.F' integer :: degree - integer :: add_double(0:64) = (/ 0, 0, 0, 0, 1, (0, i=1,60) /) - n_doubles = 1 + n_doubles = 0 do i=1,size_buffer degree = popcnt(xor( spindet, buffer(i) )) - doubles(n_doubles) = idx(i) - n_doubles = n_doubles+add_double(degree) + if (degree == 4) then + n_doubles = n_doubles+1 + doubles(n_doubles) = idx(i) + endif enddo - n_doubles = n_doubles-1 end @@ -1181,8 +1241,8 @@ subroutine get_all_spin_singles_and_doubles_$N_int(buffer, idx, spindet, size_bu integer(bit_kind) :: xorvec($N_int) integer :: degree - n_singles = 1 - n_doubles = 1 + n_singles = 0 + n_doubles = 0 do i=1,size_buffer do k=1,$N_int @@ -1196,16 +1256,14 @@ subroutine get_all_spin_singles_and_doubles_$N_int(buffer, idx, spindet, size_bu enddo if ( degree == 4 ) then - doubles(n_doubles) = idx(i) n_doubles = n_doubles+1 + doubles(n_doubles) = idx(i) else if ( degree == 2 ) then - singles(n_singles) = idx(i) n_singles = n_singles+1 + singles(n_singles) = idx(i) endif enddo - n_singles = n_singles-1 - n_doubles = n_doubles-1 end @@ -1230,7 +1288,7 @@ subroutine get_all_spin_singles_$N_int(buffer, idx, spindet, size_buffer, single integer(bit_kind) :: xorvec($N_int) integer :: degree - n_singles = 1 + n_singles = 0 do i=1,size_buffer do k=1,$N_int @@ -1247,11 +1305,10 @@ subroutine get_all_spin_singles_$N_int(buffer, idx, spindet, size_buffer, single cycle endif - singles(n_singles) = idx(i) n_singles = n_singles+1 + singles(n_singles) = idx(i) enddo - n_singles = n_singles-1 end @@ -1275,7 +1332,7 @@ subroutine get_all_spin_doubles_$N_int(buffer, idx, spindet, size_buffer, double include 'utils/constants.include.F' integer(bit_kind) :: xorvec($N_int) - n_doubles = 1 + n_doubles = 0 do i=1,size_buffer do k=1,$N_int @@ -1292,13 +1349,11 @@ subroutine get_all_spin_doubles_$N_int(buffer, idx, spindet, size_buffer, double cycle endif - doubles(n_doubles) = idx(i) n_doubles = n_doubles+1 + doubles(n_doubles) = idx(i) enddo - n_doubles = n_doubles-1 - end SUBST [ N_int ] diff --git a/src/ezfio_files/ezfio.irp.f b/src/ezfio_files/ezfio.irp.f index 02f45571..ef19d551 100644 --- a/src/ezfio_files/ezfio.irp.f +++ b/src/ezfio_files/ezfio.irp.f @@ -60,3 +60,16 @@ BEGIN_PROVIDER [ character*(1024), ezfio_work_dir ] ezfio_work_dir = trim(ezfio_filename)//'/work/' END_PROVIDER +BEGIN_PROVIDER [ character*(1024), ezfio_work_dir_pid ] + use c_functions + implicit none + BEGIN_DOC + ! EZFIO/work/pid_ + END_DOC + character*(32) :: pid_str + integer :: getpid + + write(pid_str,*) getpid() + ezfio_work_dir_pid = trim(ezfio_work_dir)//'/'//trim(pid_str)//'_' +END_PROVIDER + diff --git a/src/gpu/gpu_module.F90 b/src/gpu/gpu_module.F90 index 6050075f..90b82667 100644 --- a/src/gpu/gpu_module.F90 +++ b/src/gpu/gpu_module.F90 @@ -143,7 +143,7 @@ module gpu b, ldb, c, ldc) bind(C, name='gpu_dgeam') import type(c_ptr), value, intent(in) :: handle - character(c_char), intent(in), value :: transa, transb + character(c_char), intent(in) :: transa, transb integer(c_int64_t), intent(in), value :: m, n, lda, ldb, ldc real(c_double), intent(in) :: alpha, beta type(c_ptr), value :: a, b, c @@ -153,7 +153,7 @@ module gpu b, ldb, c, ldc) bind(C, name='gpu_sgeam') import type(c_ptr), value, intent(in) :: handle - character(c_char), intent(in), value :: transa, transb + character(c_char), intent(in) :: transa, transb integer(c_int64_t), intent(in), value :: m, n, lda, ldb, ldc real(c_float), intent(in) :: alpha, beta real(c_float) :: a, b, c @@ -194,7 +194,7 @@ module gpu b, ldb, beta, c, ldc) bind(C, name='gpu_sgemm') import type(c_ptr), value, intent(in) :: handle - character(c_char), intent(in), value :: transa, transb + character(c_char), intent(in) :: transa, transb integer(c_int64_t), intent(in), value :: m, n, k, lda, ldb, ldc real(c_float), intent(in) :: alpha, beta real(c_float) :: a, b, c @@ -268,8 +268,12 @@ module gpu implicit none type(gpu_double1), intent(inout) :: ptr integer, intent(in) :: s + integer*8 :: s_8, n - call gpu_allocate_c(ptr%c, s*8_8) + s_8 = s + n = s_8 * 8_8 + + call gpu_allocate_c(ptr%c, n) call c_f_pointer(ptr%c, ptr%f, (/ s /)) end subroutine @@ -277,8 +281,13 @@ module gpu implicit none type(gpu_double2), intent(inout) :: ptr integer, intent(in) :: s1, s2 + integer*8 :: s1_8, s2_8, n - call gpu_allocate_c(ptr%c, s1*s2*8_8) + s1_8 = s1 + s2_8 = s2 + n = s1_8 * s2_8 * 8_8 + + call gpu_allocate_c(ptr%c, n) call c_f_pointer(ptr%c, ptr%f, (/ s1, s2 /)) end subroutine @@ -286,8 +295,14 @@ module gpu implicit none type(gpu_double3), intent(inout) :: ptr integer, intent(in) :: s1, s2, s3 + integer*8 :: s1_8, s2_8, s3_8, n - call gpu_allocate_c(ptr%c, s1*s2*s3*8_8) + s1_8 = s1 + s2_8 = s2 + s3_8 = s3 + n = s1_8 * s2_8 * s3_8 * 8_8 + + call gpu_allocate_c(ptr%c, n) call c_f_pointer(ptr%c, ptr%f, (/ s1, s2, s3 /)) end subroutine @@ -295,8 +310,15 @@ module gpu implicit none type(gpu_double4), intent(inout) :: ptr integer, intent(in) :: s1, s2, s3, s4 + integer*8 :: s1_8, s2_8, s3_8, s4_8, n - call gpu_allocate_c(ptr%c, s1*s2*s3*s4*8_8) + s1_8 = s1 + s2_8 = s2 + s3_8 = s3 + s4_8 = s4 + n = s1_8 * s2_8 * s3_8 * s4_8 * 8_8 + + call gpu_allocate_c(ptr%c, n) call c_f_pointer(ptr%c, ptr%f, (/ s1, s2, s3, s4 /)) end subroutine @@ -304,8 +326,16 @@ module gpu implicit none type(gpu_double5), intent(inout) :: ptr integer, intent(in) :: s1, s2, s3, s4, s5 + integer*8 :: s1_8, s2_8, s3_8, s4_8, s5_8, n - call gpu_allocate_c(ptr%c, s1*s2*s3*s4*s5*8_8) + s1_8 = s1 + s2_8 = s2 + s3_8 = s3 + s4_8 = s4 + s5_8 = s5 + n = s1_8 * s2_8 * s3_8 * s4_8 * s5_8 * 8_8 + + call gpu_allocate_c(ptr%c, n) call c_f_pointer(ptr%c, ptr%f, (/ s1, s2, s3, s4, s5 /)) end subroutine @@ -313,8 +343,17 @@ module gpu implicit none type(gpu_double6), intent(inout) :: ptr integer, intent(in) :: s1, s2, s3, s4, s5, s6 + integer*8 :: s1_8, s2_8, s3_8, s4_8, s5_8, s6_8, n - call gpu_allocate_c(ptr%c, s1*s2*s3*s4*s5*s6*8_8) + s1_8 = s1 + s2_8 = s2 + s3_8 = s3 + s4_8 = s4 + s5_8 = s5 + s6_8 = s6 + n = s1_8 * s2_8 * s3_8 * s4_8 * s5_8 * s6_8 * 8_8 + + call gpu_allocate_c(ptr%c, n) call c_f_pointer(ptr%c, ptr%f, (/ s1, s2, s3, s4, s5, s6 /)) end subroutine diff --git a/src/hartree_fock/deb_ao_2e_int.irp.f b/src/hartree_fock/deb_ao_2e_int.irp.f new file mode 100644 index 00000000..469eb654 --- /dev/null +++ b/src/hartree_fock/deb_ao_2e_int.irp.f @@ -0,0 +1,188 @@ + +program deb_ao_2e_int + + !call check_ao_two_e_integral_cosgtos() + call check_crint1() + !call check_crint2() + +end + +! --- + +subroutine check_ao_two_e_integral_cosgtos() + + implicit none + + integer :: i, j, k, l + double precision :: tmp1, tmp2 + double precision :: acc, nrm, dif + + double precision, external :: ao_two_e_integral + double precision, external :: ao_two_e_integral_cosgtos + + acc = 0.d0 + nrm = 0.d0 + + i = 1 + j = 6 + k = 1 + l = 16 +! do i = 1, ao_num +! do k = 1, ao_num +! do j = 1, ao_num +! do l = 1, ao_num + + tmp1 = ao_two_e_integral (i, j, k, l) + tmp2 = ao_two_e_integral_cosgtos(i, j, k, l) + + dif = dabs(tmp1 - tmp2) + if(dif .gt. 1d-12) then + print*, ' error on:', i, j, k, l + print*, tmp1, tmp2, dif + stop + endif + +! acc += dif +! nrm += dabs(tmp1) +! enddo +! enddo +! enddo +! enddo + + print *, ' acc (%) = ', dif * 100.d0 / nrm + +end + +! --- + +subroutine check_crint1() + + implicit none + integer :: i, n, i_rho + double precision :: dif_thr + double precision :: dif_re, dif_im, acc_re, nrm_re, acc_im, nrm_im + complex*16 :: rho_test(1:10) = (/ (1d-12, 0.d0), & + (+1d-9, +1d-6), & + (-1d-6, -1d-5), & + (+1d-3, -1d-2), & + (-1d-1, +1d-1), & + (+1d-0, +1d-1), & + (-1d+1, +1d+1), & + (+1d+2, +1d+1), & + (-1d+3, +1d+2), & + (+1d+4, +1d+4) /) + complex*16 :: rho + complex*16 :: int_an, int_nm + double precision, external :: rint + complex*16, external :: crint_1, crint_2, crint_quad + + n = 10 + dif_thr = 1d-7 + + do i_rho = 8, 10 + !do i_rho = 7, 7 + + !rho = (-10.d0, 0.1d0) + !rho = (+10.d0, 0.1d0) + rho = rho_test(i_rho) + print*, "rho = ", real(rho), aimag(rho) + + acc_re = 0.d0 + nrm_re = 0.d0 + acc_im = 0.d0 + nrm_im = 0.d0 + do i = 0, n + !int_an = crint_1 (i, rho) + int_an = crint_2 (i, rho) + int_nm = crint_quad(i, rho) + + dif_re = dabs(real(int_an) - real(int_nm)) + dif_im = dabs(aimag(int_an) - aimag(int_nm)) + + if((dif_re .gt. dif_thr) .or. (dif_im .gt. dif_thr)) then + print*, ' error on i =', i + print*, real(int_an), real(int_nm), dif_re + print*, aimag(int_an), aimag(int_nm), dif_im + !print*, rint(i, real(rho)) + print*, crint_1(i, rho) + !print*, crint_2(i, rho) + stop + endif + acc_re += dif_re + nrm_re += dabs(real(int_nm)) + acc_im += dif_im + nrm_im += dabs(aimag(int_nm)) + enddo + + print*, "accuracy on real part (%):", 100.d0 * acc_re / (nrm_re+1d-15) + print*, "accuracy on imag part (%):", 100.d0 * acc_im / (nrm_im+1d-15) + enddo + +end + +! --- + +subroutine check_crint2() + + implicit none + + integer :: i, n, i_rho + double precision :: dif_thr + double precision :: dif_re, dif_im, acc_re, nrm_re, acc_im, nrm_im + complex*16 :: rho_test(1:10) = (/ (1d-12, 0.d0), & + (+1d-9, +1d-6), & + (-1d-6, -1d-5), & + (+1d-3, -1d-2), & + (-1d-1, +1d-1), & + (+1d-0, +1d-1), & + (-1d+1, +1d+1), & + (+1d+2, +1d+1), & + (-1d+3, +1d+2), & + (+1d+4, +1d+4) /) + complex*16 :: rho + complex*16 :: int_an, int_nm + complex*16, external :: crint_1, crint_2 + + n = 30 + dif_thr = 1d-12 + + do i_rho = 1, 10 + rho = rho_test(i_rho) + print*, "rho = ", real(rho), aimag(rho) + + acc_re = 0.d0 + nrm_re = 0.d0 + acc_im = 0.d0 + nrm_im = 0.d0 + do i = 0, n + int_an = crint_1(i, rho) + int_nm = crint_2(i, rho) + + dif_re = dabs(real(int_an) - real(int_nm)) + !if(dif_re .gt. dif_thr) then + ! print*, ' error in real part:', i + ! print*, real(int_an), real(int_nm), dif_re + ! stop + !endif + acc_re += dif_re + nrm_re += dabs(real(int_nm)) + + dif_im = dabs(aimag(int_an) - aimag(int_nm)) + !if(dif_im .gt. dif_thr) then + ! print*, ' error in imag part:', i + ! print*, aimag(int_an), aimag(int_nm), dif_im + ! stop + !endif + acc_im += dif_im + nrm_im += dabs(aimag(int_nm)) + enddo + + print*, "accuracy on real part (%):", 100.d0 * acc_re / (nrm_re+1d-15) + print*, "accuracy on imag part (%):", 100.d0 * acc_im / (nrm_im+1d-15) + enddo + +end + +! --- + + diff --git a/src/hartree_fock/hf_energy.irp.f b/src/hartree_fock/hf_energy.irp.f index a0f9f897..862d64b0 100644 --- a/src/hartree_fock/hf_energy.irp.f +++ b/src/hartree_fock/hf_energy.irp.f @@ -19,16 +19,41 @@ END_PROVIDER ! Hartree-Fock energy containing the nuclear repulsion, and its one- and two-body components. END_DOC integer :: i,j - HF_energy = nuclear_repulsion + double precision :: tmp1, tmp2 + HF_energy = 0.d0 HF_two_electron_energy = 0.d0 HF_one_electron_energy = 0.d0 do j=1,ao_num do i=1,ao_num - HF_two_electron_energy += 0.5d0 * ( ao_two_e_integral_alpha(i,j) * SCF_density_matrix_ao_alpha(i,j) & - +ao_two_e_integral_beta(i,j) * SCF_density_matrix_ao_beta(i,j) ) - HF_one_electron_energy += ao_one_e_integrals(i,j) * (SCF_density_matrix_ao_alpha(i,j) + SCF_density_matrix_ao_beta (i,j) ) + tmp1 = 0.5d0 * ( ao_two_e_integral_alpha(i,j) * SCF_density_matrix_ao_alpha(i,j) & + +ao_two_e_integral_beta (i,j) * SCF_density_matrix_ao_beta (i,j) ) + tmp2 = ao_one_e_integrals(i,j) * (SCF_density_matrix_ao_alpha(i,j) + SCF_density_matrix_ao_beta (i,j) ) + HF_two_electron_energy += tmp1 + HF_one_electron_energy += tmp2 + HF_energy += tmp1 + tmp2 + enddo + enddo + HF_energy += nuclear_repulsion +END_PROVIDER + + + BEGIN_PROVIDER [ double precision, HF_kinetic_energy] +&BEGIN_PROVIDER [ double precision, HF_n_e_energy] + implicit none + BEGIN_DOC + ! Hartree-Fock energy containing the nuclear repulsion, and its one- and two-body components. + END_DOC + integer :: i,j + double precision :: tmp1, tmp2 + HF_n_e_energy = 0.d0 + HF_kinetic_energy = 0.d0 + do j=1,ao_num + do i=1,ao_num + tmp1 = ao_integrals_n_e(i,j) * (SCF_density_matrix_ao_alpha(i,j) + SCF_density_matrix_ao_beta (i,j) ) + tmp2 = ao_kinetic_integrals(i,j) * (SCF_density_matrix_ao_alpha(i,j) + SCF_density_matrix_ao_beta (i,j) ) + HF_n_e_energy += tmp1 + HF_kinetic_energy += tmp2 enddo enddo - HF_energy += HF_two_electron_energy + HF_one_electron_energy END_PROVIDER diff --git a/src/hartree_fock/print_scf_int.irp.f b/src/hartree_fock/print_scf_int.irp.f deleted file mode 100644 index ee7590f6..00000000 --- a/src/hartree_fock/print_scf_int.irp.f +++ /dev/null @@ -1,114 +0,0 @@ - -program print_scf_int - - call main() - -end - -subroutine main() - - implicit none - integer :: i, j, k, l - - print *, " Hcore:" - do j = 1, ao_num - do i = 1, ao_num - print *, i, j, ao_one_e_integrals(i,j) - enddo - enddo - - print *, " P:" - do j = 1, ao_num - do i = 1, ao_num - print *, i, j, SCF_density_matrix_ao_alpha(i,j) - enddo - enddo - - - double precision :: integ, density_a, density_b, density - double precision :: J_scf(ao_num, ao_num) - double precision :: K_scf(ao_num, ao_num) - - - double precision, external :: get_ao_two_e_integral - PROVIDE ao_integrals_map - - print *, " J:" - !do j = 1, ao_num - ! do l = 1, ao_num - ! do i = 1, ao_num - ! do k = 1, ao_num - ! ! < 1:k, 2:l | 1:i, 2:j > - ! print *, '< k l | i j >', k, l, i, j - ! print *, get_ao_two_e_integral(i, j, k, l, ao_integrals_map) - ! enddo - ! enddo - ! enddo - !enddo - - !do k = 1, ao_num - ! do i = 1, ao_num - ! do j = 1, ao_num - ! do l = 1, ao_num - ! ! ( 1:k, 1:i | 2:l, 2:j ) - ! print *, '(k i | l j)', k, i, l, j - ! print *, get_ao_two_e_integral(l, j, k, i, ao_integrals_map) - ! enddo - ! enddo - ! print *, '' - ! enddo - !enddo - - J_scf = 0.d0 - K_scf = 0.d0 - do i = 1, ao_num - do k = 1, ao_num - do j = 1, ao_num - do l = 1, ao_num - - density_a = SCF_density_matrix_ao_alpha(l,j) - density_b = SCF_density_matrix_ao_beta (l,j) - density = density_a + density_b - - integ = get_ao_two_e_integral(l, j, k, i, ao_integrals_map) - J_scf(k,i) += density * integ - integ = get_ao_two_e_integral(l, i, k, j, ao_integrals_map) - K_scf(k,i) -= density_a * integ - enddo - enddo - enddo - enddo - - print *, 'J x P' - do i = 1, ao_num - do k = 1, ao_num - print *, k, i, J_scf(k,i) - enddo - enddo - - print *, '' - print *, 'K x P' - do i = 1, ao_num - do k = 1, ao_num - print *, k, i, K_scf(k,i) - enddo - enddo - - print *, '' - print *, 'F in AO' - do i = 1, ao_num - do k = 1, ao_num - print *, k, i, Fock_matrix_ao(k,i) - enddo - enddo - - print *, '' - print *, 'F in MO' - do i = 1, ao_num - do k = 1, ao_num - print *, k, i, 2.d0 * Fock_matrix_mo_alpha(k,i) - enddo - enddo - -end - diff --git a/src/mo_basis/mos.irp.f b/src/mo_basis/mos.irp.f index 57ebb533..1ed859ee 100644 --- a/src/mo_basis/mos.irp.f +++ b/src/mo_basis/mos.irp.f @@ -277,7 +277,7 @@ subroutine ao_to_mo(A_ao,LDA_ao,A_mo,LDA_mo) T, ao_num, & 0.d0, A_mo, size(A_mo,1)) - call restore_symmetry(mo_num,mo_num,A_mo,size(A_mo,1),1.d-12) + call restore_symmetry(mo_num,mo_num,A_mo,size(A_mo,1),1.d-15) deallocate(T) end diff --git a/src/mo_one_e_ints/mo_one_e_ints.irp.f b/src/mo_one_e_ints/mo_one_e_ints.irp.f index 926ac7bd..fba55beb 100644 --- a/src/mo_one_e_ints/mo_one_e_ints.irp.f +++ b/src/mo_one_e_ints/mo_one_e_ints.irp.f @@ -18,6 +18,6 @@ BEGIN_PROVIDER [ double precision, mo_one_e_integrals,(mo_num,mo_num)] call ezfio_set_mo_one_e_ints_mo_one_e_integrals(mo_one_e_integrals) print *, 'MO one-e integrals written to disk' ENDIF - call nullify_small_elements(mo_num,mo_num,mo_one_e_integrals,size(mo_one_e_integrals,1),1.d-10) + call nullify_small_elements(mo_num,mo_num,mo_one_e_integrals,size(mo_one_e_integrals,1),1.d-15) END_PROVIDER diff --git a/src/mo_two_e_ints/mo_bi_integrals.irp.f b/src/mo_two_e_ints/mo_bi_integrals.irp.f index 549bbed2..6ee23f4a 100644 --- a/src/mo_two_e_ints/mo_bi_integrals.irp.f +++ b/src/mo_two_e_ints/mo_bi_integrals.irp.f @@ -70,6 +70,10 @@ BEGIN_PROVIDER [ logical, mo_two_e_integrals_in_map ] else call add_integrals_to_map(full_ijkl_bitmask_4) endif + double precision, external :: map_mb + print*,'Molecular integrals provided:' + print*,' Size of MO map ', map_mb(mo_integrals_map) ,'MB' + print*,' Number of MO integrals: ', mo_map_size endif call wall_time(wall_2) @@ -78,10 +82,6 @@ BEGIN_PROVIDER [ logical, mo_two_e_integrals_in_map ] integer*8 :: get_mo_map_size, mo_map_size mo_map_size = get_mo_map_size() - double precision, external :: map_mb - print*,'Molecular integrals provided:' - print*,' Size of MO map ', map_mb(mo_integrals_map) ,'MB' - print*,' Number of MO integrals: ', mo_map_size print*,' cpu time :',cpu_2 - cpu_1, 's' print*,' wall time :',wall_2 - wall_1, 's ( x ', (cpu_2-cpu_1)/(wall_2-wall_1), ')' diff --git a/src/mol_properties/print_properties.irp.f b/src/mol_properties/print_properties.irp.f index af413a88..c392c148 100644 --- a/src/mol_properties/print_properties.irp.f +++ b/src/mol_properties/print_properties.irp.f @@ -100,7 +100,7 @@ subroutine print_transition_dipole_moment dip_str = d_x**2 + d_y**2 + d_z**2 d = multi_s_dipole_moment(istate,jstate) f = 2d0/3d0 * d * d * dabs(ci_energy_no_diag(istate) - ci_energy_no_diag(jstate)) - write(*,'(I4,I4,A4,I3,6(F12.6))') (istate-1), (jstate-1), ' ->', (istate-1), d_x, d_y, d_z, d, dip_str, f + write(*,'(I4,I4,A4,I3,6(F12.6))') (jstate -1) * (2*N_states-jstate)/2 + istate - jstate, (jstate-1), ' ->', (istate-1), d_x, d_y, d_z, d, dip_str, f enddo enddo @@ -117,7 +117,7 @@ subroutine print_transition_dipole_moment dip_str = d_x**2 + d_y**2 + d_z**2 f = 2d0/3d0 * d * d * dabs(ci_energy_no_diag(istate) - ci_energy_no_diag(jstate)) d = multi_s_dipole_moment(istate,jstate) * au_to_D - write(*,'(I4,I4,A4,I3,6(F12.6))') (istate-1), (jstate-1), ' ->', (istate-1), d_x, d_y, d_z, d, dip_str, f + write(*,'(I4,I4,A4,I3,6(F12.6))') (jstate -1) * (2*N_states-jstate)/2 + istate - jstate, (jstate-1), ' ->', (istate-1), d_x, d_y, d_z, d, dip_str, f enddo enddo print*,'==============================================' @@ -181,10 +181,9 @@ subroutine print_oscillator_strength ! Mixed gauge f_m = 2d0/3d0 * d * v - write(*,'(A19,I3,A9,F10.6,A5,F7.1,A10,F9.6,A6,F9.6,A6,F9.6,A8,F7.3)') ' # Transition n.', (istate-1), ': Excit.=', dabs((ci_energy_no_diag(istate) - ci_energy_no_diag(jstate)))*ha_to_ev, & + write(*,'(A19,I3,A9,F10.6,A5,F7.1,A10,F9.6,A6,F9.6,A6,F9.6,A8,F7.3)') ' # Transition n.', (jstate -1) * (2*N_states-jstate)/2 + istate - jstate, ': Excit.=', dabs((ci_energy_no_diag(istate) - ci_energy_no_diag(jstate)))*ha_to_ev, & ' eV ( ',dabs((ci_energy_no_diag(istate) - ci_energy_no_diag(jstate)))*Ha_to_nm,' nm), f_l=',f_l, ', f_v=', f_v, ', f_m=', f_m, ', =', s2_values(istate) !write(*,'(I4,I4,A4,I3,A6,F6.1,A6,F6.1)') (istate-1), (jstate-1), ' ->', (istate-1), ', %T1=', percent_exc(2,istate), ', %T2=',percent_exc(3,istate) - enddo enddo diff --git a/src/mu_of_r/test_proj_op.irp.f b/src/mu_of_r/test_proj_op.irp.f index cf53c772..94052b18 100644 --- a/src/mu_of_r/test_proj_op.irp.f +++ b/src/mu_of_r/test_proj_op.irp.f @@ -12,15 +12,12 @@ program projected_operators mu_of_r_potential = "cas_full" touch mu_of_r_potential print*,'Using Valence Only functions' -! call test_f_HF_valence_ab -! call routine_full_mos -! call test_f_ii_valence_ab -! call test_f_ia_valence_ab -! call test_f_ii_ia_aa_valence_ab -! call test -! call test_f_mean_field -! call test_grad_f_mean_field - call test_grad_mu_mf + call test_f_HF_valence_ab + call routine_full_mos + call test_f_ii_valence_ab + call test_f_ia_valence_ab + call test_f_ii_ia_aa_valence_ab + call test end @@ -39,138 +36,3 @@ subroutine test end -subroutine test_f_mean_field - implicit none - integer :: i_point - double precision :: weight,r(3) - double precision :: ref_f, new_f, accu_f - double precision :: ref_two_dens, new_two_dens, accu_two_dens, dm_a, dm_b - accu_f = 0.d0 - accu_two_dens = 0.d0 - do i_point = 1, n_points_final_grid - r(1:3) = final_grid_points(1:3,i_point) - weight = final_weight_at_r_vector(i_point) - call get_f_mf_ab(r,new_f,new_two_dens, dm_a, dm_b) - call f_HF_valence_ab(r,r,ref_f,ref_two_dens) - accu_f += weight * dabs(new_f- ref_f) - accu_two_dens += weight * dabs(new_two_dens - ref_two_dens) - enddo - print*,'accu_f = ',accu_f - print*,'accu_two_dens = ',accu_two_dens - -end - -subroutine test_grad_f_mean_field - implicit none - integer :: i_point,k - double precision :: weight,r(3) - double precision :: grad_f_mf_ab(3), grad_two_bod_dens(3) - double precision :: grad_dm_a(3), grad_dm_b(3) - double precision :: f_mf_ab,two_bod_dens, dm_a, dm_b - - double precision :: num_grad_f_mf_ab(3), num_grad_two_bod_dens(3) - double precision :: num_grad_dm_a(3), num_grad_dm_b(3) - double precision :: f_mf_ab_p,f_mf_ab_m - double precision :: two_bod_dens_p, two_bod_dens_m - double precision :: dm_a_p, dm_a_m - double precision :: dm_b_p, dm_b_m - double precision :: rbis(3), dr - double precision :: accu_grad_f_mf_ab(3),accu_grad_two_bod_dens(3) - double precision :: accu_grad_dm_a(3),accu_grad_dm_b(3) - double precision :: accu_f_mf_ab, accu_two_bod_dens, accu_dm_a, accu_dm_b - dr = 0.00001d0 - accu_f_mf_ab = 0.d0 - accu_two_bod_dens = 0.d0 - accu_dm_a = 0.d0 - accu_dm_b = 0.d0 - - accu_grad_f_mf_ab = 0.d0 - accu_grad_two_bod_dens = 0.d0 - accu_grad_dm_a = 0.d0 - accu_grad_dm_b = 0.d0 - do i_point = 1, n_points_final_grid - r(1:3) = final_grid_points(1:3,i_point) - weight = final_weight_at_r_vector(i_point) - call get_grad_f_mf_ab(r,grad_f_mf_ab, grad_two_bod_dens,f_mf_ab,two_bod_dens, dm_a, dm_b,grad_dm_a, grad_dm_b) - call get_f_mf_ab(r,f_mf_ab_p,two_bod_dens_p, dm_a_p, dm_b_p) - accu_f_mf_ab += weight * dabs(f_mf_ab - f_mf_ab_p) - accu_two_bod_dens += weight * dabs(two_bod_dens - two_bod_dens_p) - accu_dm_a += weight*dabs(dm_a - dm_a_p) - accu_dm_b += weight*dabs(dm_b - dm_b_p) - do k = 1, 3 - rbis = r - rbis(k) += dr - call get_f_mf_ab(rbis,f_mf_ab_p,two_bod_dens_p, dm_a_p, dm_b_p) - rbis = r - rbis(k) -= dr - call get_f_mf_ab(rbis,f_mf_ab_m,two_bod_dens_m, dm_a_m, dm_b_m) - num_grad_f_mf_ab(k) = (f_mf_ab_p - f_mf_ab_m)/(2.d0*dr) - num_grad_two_bod_dens(k) = (two_bod_dens_p - two_bod_dens_m)/(2.d0*dr) - num_grad_dm_a(k) = (dm_a_p - dm_a_m)/(2.d0*dr) - num_grad_dm_b(k) = (dm_b_p - dm_b_m)/(2.d0*dr) - enddo - do k = 1, 3 - accu_grad_f_mf_ab(k) += weight * dabs(grad_f_mf_ab(k) - num_grad_f_mf_ab(k)) - accu_grad_two_bod_dens(k) += weight * dabs(grad_two_bod_dens(k) - num_grad_two_bod_dens(k)) - accu_grad_dm_a(k) += weight * dabs(grad_dm_a(k) - num_grad_dm_a(k)) - accu_grad_dm_b(k) += weight * dabs(grad_dm_b(k) - num_grad_dm_b(k)) - enddo - enddo - print*,'accu_f_mf_ab = ',accu_f_mf_ab - print*,'accu_two_bod_dens = ',accu_two_bod_dens - print*,'accu_dm_a = ',accu_dm_a - print*,'accu_dm_b = ',accu_dm_b - print*,'accu_grad_f_mf_ab = ' - print*,accu_grad_f_mf_ab - print*,'accu_grad_two_bod_dens = ' - print*,accu_grad_two_bod_dens - print*,'accu_dm_a = ' - print*,accu_grad_dm_a - print*,'accu_dm_b = ' - print*,accu_grad_dm_b - -end - -subroutine test_grad_mu_mf - implicit none - integer :: i_point,k - double precision :: weight,r(3),rbis(3) - double precision :: mu_mf, dm,grad_mu_mf(3), grad_dm(3) - double precision :: mu_mf_p, mu_mf_m, dm_m, dm_p, num_grad_mu_mf(3),dr, num_grad_dm(3) - double precision :: accu_mu, accu_dm, accu_grad_dm(3), accu_grad_mu_mf(3) - dr = 0.00001d0 - accu_grad_mu_mf = 0.d0 - accu_mu = 0.d0 - accu_grad_dm = 0.d0 - accu_dm = 0.d0 - do i_point = 1, n_points_final_grid - r(1:3) = final_grid_points(1:3,i_point) - weight = final_weight_at_r_vector(i_point) - call grad_mu_of_r_mean_field(r,mu_mf, dm, grad_mu_mf, grad_dm) - call mu_of_r_mean_field(r,mu_mf_p, dm_p) - accu_mu += weight*dabs(mu_mf_p - mu_mf) - accu_dm += weight*dabs(dm_p - dm) - do k = 1, 3 - rbis = r - rbis(k) += dr - call mu_of_r_mean_field(rbis,mu_mf_p, dm_p) - rbis = r - rbis(k) -= dr - call mu_of_r_mean_field(rbis,mu_mf_m, dm_m) - - num_grad_mu_mf(k) = (mu_mf_p - mu_mf_m)/(2.d0*dr) - num_grad_dm(k) = (dm_p - dm_m)/(2.d0*dr) - enddo - do k = 1, 3 - accu_grad_dm(k)+= weight *dabs(num_grad_dm(k) - grad_dm(k)) - accu_grad_mu_mf(k)+= weight *dabs(num_grad_mu_mf(k) - grad_mu_mf(k)) - enddo - enddo - print*,'accu_mu = ',accu_mu - print*,'accu_dm = ',accu_dm - print*,'accu_grad_dm = ' - print*, accu_grad_dm - print*,'accu_grad_mu_mf = ' - print*, accu_grad_mu_mf - -end diff --git a/src/scf_utils/EZFIO.cfg b/src/scf_utils/EZFIO.cfg index 7b950b14..1ca97217 100644 --- a/src/scf_utils/EZFIO.cfg +++ b/src/scf_utils/EZFIO.cfg @@ -45,6 +45,12 @@ type: double precision doc: Calculated HF energy interface: ezfio +[do_mom] +type: logical +doc: If true, this will run a MOM calculation. The overlap will be computed at each step with respect to the initial MOs. After an initial Hartree-Fock calculation, the guess can be created by swapping molecular orbitals through the qp run swap_mos command. +interface: ezfio,provider,ocaml +default: False + [frozen_orb_scf] type: logical doc: If true, leave untouched all the orbitals defined as core and optimize all the orbitals defined as active with qp_set_mo_class diff --git a/src/scf_utils/fock_matrix.irp.f b/src/scf_utils/fock_matrix.irp.f index 269a441b..051da9eb 100644 --- a/src/scf_utils/fock_matrix.irp.f +++ b/src/scf_utils/fock_matrix.irp.f @@ -253,17 +253,18 @@ BEGIN_PROVIDER [ double precision, SCF_energy ] BEGIN_DOC ! Hartree-Fock energy END_DOC - SCF_energy = nuclear_repulsion - integer :: i,j + + SCF_energy = 0.d0 + do j=1,ao_num do i=1,ao_num - SCF_energy += 0.5d0 * ( & + SCF_energy += & (ao_one_e_integrals(i,j) + Fock_matrix_ao_alpha(i,j) ) * SCF_density_matrix_ao_alpha(i,j) +& - (ao_one_e_integrals(i,j) + Fock_matrix_ao_beta (i,j) ) * SCF_density_matrix_ao_beta (i,j) ) + (ao_one_e_integrals(i,j) + Fock_matrix_ao_beta (i,j) ) * SCF_density_matrix_ao_beta (i,j) enddo enddo - SCF_energy += extra_e_contrib_density + SCF_energy = 0.5d0 * SCF_energy + extra_e_contrib_density + nuclear_repulsion END_PROVIDER diff --git a/src/scf_utils/reorder_mo_max_overlap.irp.f b/src/scf_utils/reorder_mo_max_overlap.irp.f new file mode 100644 index 00000000..29f16735 --- /dev/null +++ b/src/scf_utils/reorder_mo_max_overlap.irp.f @@ -0,0 +1,96 @@ + subroutine reorder_mo_max_overlap + implicit none + BEGIN_DOC + ! routines that compute the projection of each MO of the current `mo_coef` on the space spanned by the occupied orbitals of `mo_coef_begin_iteration` + END_DOC + integer :: i,j,k,l + double precision, allocatable :: overlap(:,:) + double precision, allocatable :: proj(:) + integer, allocatable :: iorder(:) + double precision, allocatable :: mo_coef_tmp(:,:) + double precision, allocatable :: tmp(:,:) + allocate(overlap(mo_num,mo_num),proj(mo_num),iorder(mo_num),mo_coef_tmp(ao_num,mo_num),tmp(mo_num,ao_num)) + + overlap(:,:) = 0d0 + mo_coef_tmp(:,:) = 0d0 + proj(:) = 0d0 + iorder(:) = 0d0 + tmp(:,:) = 0d0 + + ! These matrix products compute the overlap bewteen the initial and the current MOs + call dgemm('T','N', mo_num, ao_num, ao_num, 1.d0, & + mo_coef_begin_iteration, size(mo_coef_begin_iteration,1), & + ao_overlap, size(ao_overlap,1), 0.d0, & + tmp, size(tmp,1)) + + call dgemm('N','N', mo_num, mo_num, ao_num, 1.d0, & + tmp, size(tmp,1), & + mo_coef, size(mo_coef, 1), 0.d0, & + overlap, size(overlap,1) ) + + + ! for each orbital compute the best overlap + do i = 1, mo_num + iorder(i) = i ! initialize the iorder list as we need it to sort later + do j = 1, elec_alpha_num + proj(i) += overlap(j,i)*overlap(j,i) ! compute the projection of current orbital i on the occupied space of the initial orbitals + enddo + proj(i) = dsqrt(proj(i)) + enddo + ! sort the list of projection to find the mos with the largest overlap + call dsort(proj(:),iorder(:),mo_num) + ! reorder orbitals according to projection + do i=1,mo_num + mo_coef_tmp(:,i) = mo_coef(:,iorder(mo_num+1-i)) + enddo + + ! update the orbitals + mo_coef(:,:) = mo_coef_tmp(:,:) + + ! if the determinant is open-shell we need to make sure that the singly occupied orbital correspond to the initial ones + if (elec_alpha_num > elec_beta_num) then + double precision, allocatable :: overlap_alpha(:,:) + double precision, allocatable :: proj_alpha(:) + integer, allocatable :: iorder_alpha(:) + allocate(overlap_alpha(mo_num,elec_alpha_num),proj_alpha(elec_alpha_num),iorder_alpha(elec_alpha_num)) + overlap_alpha(:,:) = 0d0 + mo_coef_tmp(:,:) = 0d0 + proj_alpha(:) = 0d0 + iorder_alpha(:) = 0d0 + tmp(:,:) = 0d0 + ! These matrix products compute the overlap bewteen the initial and the current MOs + call dgemm('T','N', mo_num, ao_num, ao_num, 1.d0, & + mo_coef_begin_iteration, size(mo_coef_begin_iteration,1), & + ao_overlap, size(ao_overlap,1), 0.d0, & + tmp, size(tmp,1)) + + call dgemm('N','N', mo_num, elec_alpha_num, ao_num, 1.d0, & + tmp, size(tmp,1), & + mo_coef, size(mo_coef, 1), 0.d0, & + overlap_alpha, size(overlap_alpha,1) ) + + do i = 1, elec_alpha_num + iorder_alpha(i) = i ! initialize the iorder list as we need it to sort later + do j = 1, elec_beta_num + proj_alpha(i) += overlap_alpha(j,i)*overlap_alpha(j,i) ! compute the projection of current orbital i on the beta occupied space of the initial orbitals + enddo + proj_alpha(i) = dsqrt(proj_alpha(i)) + enddo + + ! sort the list of projection to find the mos with the largest overlap + call dsort(proj_alpha(:),iorder_alpha(:),elec_alpha_num) + ! reorder orbitals according to projection + do i=1,elec_alpha_num + mo_coef_tmp(:,i) = mo_coef(:,iorder_alpha(elec_alpha_num+1-i)) + enddo + do i=1,elec_alpha_num + mo_coef(:,i) = mo_coef_tmp(:,i) + enddo + + deallocate(overlap_alpha, proj_alpha, iorder_alpha) + endif + + deallocate(overlap, proj, iorder, mo_coef_tmp, tmp) + + end + diff --git a/src/scf_utils/roothaan_hall_scf.irp.f b/src/scf_utils/roothaan_hall_scf.irp.f index e0fe5319..947917af 100644 --- a/src/scf_utils/roothaan_hall_scf.irp.f +++ b/src/scf_utils/roothaan_hall_scf.irp.f @@ -51,6 +51,11 @@ END_DOC ! PROVIDE FPS_SPF_matrix_AO Fock_matrix_AO + ! Initialize MO to run IMOM + if(do_mom)then + call initialize_mo_coef_begin_iteration + endif + converged = .False. do while ( .not.converged .and. (iteration_SCF < n_it_SCF_max) ) @@ -88,16 +93,17 @@ END_DOC Fock_matrix_AO_beta = Fock_matrix_AO*0.5d0 TOUCH Fock_matrix_AO_alpha Fock_matrix_AO_beta - endif - + endif MO_coef = eigenvectors_Fock_matrix_MO + if(do_mom)then + call reorder_mo_max_overlap + endif if(frozen_orb_scf)then - call reorder_core_orb - call initialize_mo_coef_begin_iteration + call reorder_core_orb + call initialize_mo_coef_begin_iteration endif TOUCH MO_coef - ! Calculate error vectors max_error_DIIS = maxval(Abs(FPS_SPF_Matrix_MO)) @@ -106,41 +112,46 @@ END_DOC energy_SCF = SCF_energy Delta_Energy_SCF = energy_SCF - energy_SCF_previous - if ( (SCF_algorithm == 'DIIS').and.(Delta_Energy_SCF > 0.d0) ) then + if ( (SCF_algorithm == 'DIIS').and.(Delta_Energy_SCF > 0.d0).and.(.not.do_mom) ) then Fock_matrix_AO(1:ao_num,1:ao_num) = Fock_matrix_DIIS (1:ao_num,1:ao_num,index_dim_DIIS) Fock_matrix_AO_alpha = Fock_matrix_AO*0.5d0 Fock_matrix_AO_beta = Fock_matrix_AO*0.5d0 TOUCH Fock_matrix_AO_alpha Fock_matrix_AO_beta endif - double precision :: level_shift_save - level_shift_save = level_shift - mo_coef_save(1:ao_num,1:mo_num) = mo_coef(1:ao_num,1:mo_num) - do while (Delta_energy_SCF > 0.d0) - mo_coef(1:ao_num,1:mo_num) = mo_coef_save - if (level_shift <= .1d0) then - level_shift = 1.d0 - else - level_shift = level_shift * 3.0d0 - endif - TOUCH mo_coef level_shift - mo_coef(1:ao_num,1:mo_num) = eigenvectors_Fock_matrix_MO(1:ao_num,1:mo_num) - if(frozen_orb_scf)then - call reorder_core_orb - call initialize_mo_coef_begin_iteration - endif - TOUCH mo_coef - Delta_Energy_SCF = SCF_energy - energy_SCF_previous - energy_SCF = SCF_energy - if (level_shift-level_shift_save > 40.d0) then - level_shift = level_shift_save * 4.d0 - SOFT_TOUCH level_shift - exit - endif - dim_DIIS=0 - enddo - level_shift = level_shift * 0.5d0 - SOFT_TOUCH level_shift + if (.not.do_mom) then + double precision :: level_shift_save + level_shift_save = level_shift + mo_coef_save(1:ao_num,1:mo_num) = mo_coef(1:ao_num,1:mo_num) + do while (Delta_energy_SCF > 0.d0) + mo_coef(1:ao_num,1:mo_num) = mo_coef_save + if (level_shift <= .1d0) then + level_shift = 1.d0 + else + level_shift = level_shift * 3.0d0 + endif + TOUCH mo_coef level_shift + mo_coef(1:ao_num,1:mo_num) = eigenvectors_Fock_matrix_MO(1:ao_num,1:mo_num) + if(do_mom)then + call reorder_mo_max_overlap + endif + if(frozen_orb_scf)then + call reorder_core_orb + call initialize_mo_coef_begin_iteration + endif + TOUCH mo_coef + Delta_Energy_SCF = SCF_energy - energy_SCF_previous + energy_SCF = SCF_energy + if (level_shift-level_shift_save > 40.d0) then + level_shift = level_shift_save * 4.d0 + SOFT_TOUCH level_shift + exit + endif + dim_DIIS=0 + enddo + level_shift = level_shift * 0.5d0 + SOFT_TOUCH level_shift + endif energy_SCF_previous = energy_SCF converged = ( (max_error_DIIS <= threshold_DIIS_nonzero) .and. & @@ -205,7 +216,7 @@ END_DOC if(.not.frozen_orb_scf)then call mo_as_eigvectors_of_mo_matrix(Fock_matrix_mo,size(Fock_matrix_mo,1), & - size(Fock_matrix_mo,2),mo_label,1,.true.) + size(Fock_matrix_mo,2),mo_label,1,.true.) call restore_symmetry(ao_num, mo_num, mo_coef, size(mo_coef,1), 1.d-10) call orthonormalize_mos endif @@ -228,6 +239,9 @@ END_DOC i = j+1 enddo + if(do_mom)then + call reorder_mo_max_overlap + endif call save_mos diff --git a/src/tools/print_energy.irp.f b/src/tools/print_energy.irp.f index 0e67828e..485a0910 100644 --- a/src/tools/print_energy.irp.f +++ b/src/tools/print_energy.irp.f @@ -15,5 +15,21 @@ end subroutine run implicit none call print_mol_properties - print *, psi_energy + nuclear_repulsion + print *, psi_energy + nuclear_repulsion + call print_energy_components + print *, 'E(HF) = ', HF_energy + print *, 'E(CI) = ', psi_energy + nuclear_repulsion + print *, '' + print *, 'E_kin(CI) = ', ref_bitmask_kinetic_energy + print *, 'E_kin(HF) = ', HF_kinetic_energy + print *, '' + print *, 'E_ne (CI) = ', ref_bitmask_n_e_energy + print *, 'E_ne (HF) = ', HF_n_e_energy + print *, '' + print *, 'E_1e (CI) = ', ref_bitmask_one_e_energy + print *, 'E_1e (HF) = ', HF_one_electron_energy + print *, '' + print *, 'E_2e (CI) = ', ref_bitmask_two_e_energy + print *, 'E_2e (HF) = ', HF_two_electron_energy + end diff --git a/src/tools/truncate_wf.irp.f b/src/tools/truncate_wf.irp.f index 6c66c8ec..1556e7d8 100644 --- a/src/tools/truncate_wf.irp.f +++ b/src/tools/truncate_wf.irp.f @@ -56,9 +56,14 @@ subroutine routine_s2 double precision :: accu(N_states) print *, 'Weights of the CFG' - do i=1,N_det + integer :: step + + step = max(1,N_det/100) + do i=1,N_det-1,step print *, i, real(weight_configuration(det_to_configuration(i),:)), real(sum(weight_configuration(det_to_configuration(i),:))) enddo + i=N_det + print *, i, real(weight_configuration(det_to_configuration(i),:)), real(sum(weight_configuration(det_to_configuration(i),:))) print*, 'Min weight of the configuration?' read(5,*) wmin diff --git a/src/trexio/import_trexio_integrals.irp.f b/src/trexio/import_trexio_integrals.irp.f index 556ed7bc..405c0863 100644 --- a/src/trexio/import_trexio_integrals.irp.f +++ b/src/trexio/import_trexio_integrals.irp.f @@ -46,6 +46,8 @@ subroutine run(f) double precision, allocatable :: tmp(:,:,:) integer*8 :: offset, icount + integer :: k_num + integer, external :: getUnitAndOpen if (trexio_has_nucleus_repulsion(f) == TREXIO_SUCCESS) then @@ -163,46 +165,48 @@ subroutine run(f) deallocate(Vi, V, tmp) print *, 'Cholesky AO integrals read from TREXIO file' - endif - rc = trexio_has_ao_2e_int_eri(f) - if (rc /= TREXIO_HAS_NOT) then - PROVIDE ao_integrals_map + else - BUFSIZE=ao_num**2 - allocate(buffer_i(BUFSIZE), buffer_values(BUFSIZE)) - allocate(Vi(4,BUFSIZE), V(BUFSIZE)) + rc = trexio_has_ao_2e_int_eri(f) + if (rc /= TREXIO_HAS_NOT) then + PROVIDE ao_integrals_map - offset = 0_8 - icount = BUFSIZE - rc = TREXIO_SUCCESS - do while (icount == size(V)) - rc = trexio_read_ao_2e_int_eri(f, offset, icount, Vi, V) - do m=1,icount - i = Vi(1,m) - j = Vi(2,m) - k = Vi(3,m) - l = Vi(4,m) - integral = V(m) - call two_e_integrals_index(i, j, k, l, buffer_i(m) ) - buffer_values(m) = integral - enddo - call insert_into_ao_integrals_map(int(icount,4),buffer_i,buffer_values) - offset = offset + icount - if (rc /= TREXIO_SUCCESS) then - exit - endif - end do - n_integrals = offset + BUFSIZE=ao_num**2 + allocate(buffer_i(BUFSIZE), buffer_values(BUFSIZE)) + allocate(Vi(4,BUFSIZE), V(BUFSIZE)) - call map_sort(ao_integrals_map) - call map_unique(ao_integrals_map) + offset = 0_8 + icount = BUFSIZE + rc = TREXIO_SUCCESS + do while (icount == size(V)) + rc = trexio_read_ao_2e_int_eri(f, offset, icount, Vi, V) + do m=1,icount + i = Vi(1,m) + j = Vi(2,m) + k = Vi(3,m) + l = Vi(4,m) + integral = V(m) + call two_e_integrals_index(i, j, k, l, buffer_i(m) ) + buffer_values(m) = integral + enddo + call insert_into_ao_integrals_map(int(icount,4),buffer_i,buffer_values) + offset = offset + icount + if (rc /= TREXIO_SUCCESS) then + exit + endif + end do + n_integrals = offset - call map_save_to_disk(trim(ezfio_filename)//'/work/ao_ints',ao_integrals_map) - call ezfio_set_ao_two_e_ints_io_ao_two_e_integrals('Read') + call map_sort(ao_integrals_map) + call map_unique(ao_integrals_map) - deallocate(buffer_i, buffer_values, Vi, V) - print *, 'AO integrals read from TREXIO file' + call map_save_to_disk(trim(ezfio_filename)//'/work/ao_ints',ao_integrals_map) + call ezfio_set_ao_two_e_ints_io_ao_two_e_integrals('Read') + + deallocate(buffer_i, buffer_values, Vi, V) + print *, 'AO integrals read from TREXIO file' + endif endif else print *, 'AO integrals not found in TREXIO file' @@ -270,44 +274,47 @@ subroutine run(f) deallocate(Vi, V, tmp) print *, 'Cholesky MO integrals read from TREXIO file' - endif - rc = trexio_has_mo_2e_int_eri(f) - if (rc /= TREXIO_HAS_NOT) then - BUFSIZE=mo_num**2 - allocate(buffer_i(BUFSIZE), buffer_values(BUFSIZE)) - allocate(Vi(4,BUFSIZE), V(BUFSIZE)) + else + + rc = trexio_has_mo_2e_int_eri(f) + if (rc /= TREXIO_HAS_NOT) then + BUFSIZE=mo_num**2 + allocate(buffer_i(BUFSIZE), buffer_values(BUFSIZE)) + allocate(Vi(4,BUFSIZE), V(BUFSIZE)) - offset = 0_8 - icount = BUFSIZE - rc = TREXIO_SUCCESS - do while (icount == size(V)) - rc = trexio_read_mo_2e_int_eri(f, offset, icount, Vi, V) - do m=1,icount - i = Vi(1,m) - j = Vi(2,m) - k = Vi(3,m) - l = Vi(4,m) - integral = V(m) - call two_e_integrals_index(i, j, k, l, buffer_i(m) ) - buffer_values(m) = integral - enddo - call map_append(mo_integrals_map, buffer_i, buffer_values, int(icount,4)) - offset = offset + icount - if (rc /= TREXIO_SUCCESS) then - exit - endif - end do - n_integrals = offset + offset = 0_8 + icount = BUFSIZE + rc = TREXIO_SUCCESS + do while (icount == size(V)) + rc = trexio_read_mo_2e_int_eri(f, offset, icount, Vi, V) + do m=1,icount + i = Vi(1,m) + j = Vi(2,m) + k = Vi(3,m) + l = Vi(4,m) + integral = V(m) + call two_e_integrals_index(i, j, k, l, buffer_i(m) ) + buffer_values(m) = integral + enddo + call map_append(mo_integrals_map, buffer_i, buffer_values, int(icount,4)) + offset = offset + icount + if (rc /= TREXIO_SUCCESS) then + exit + endif + end do + n_integrals = offset - call map_sort(mo_integrals_map) - call map_unique(mo_integrals_map) + call map_sort(mo_integrals_map) + call map_unique(mo_integrals_map) + + call map_save_to_disk(trim(ezfio_filename)//'/work/mo_ints',mo_integrals_map) + call ezfio_set_mo_two_e_ints_io_mo_two_e_integrals('Read') + deallocate(buffer_i, buffer_values, Vi, V) + print *, 'MO integrals read from TREXIO file' + endif - call map_save_to_disk(trim(ezfio_filename)//'/work/mo_ints',mo_integrals_map) - call ezfio_set_mo_two_e_ints_io_mo_two_e_integrals('Read') - deallocate(buffer_i, buffer_values, Vi, V) - print *, 'MO integrals read from TREXIO file' endif else diff --git a/src/utils/cgtos_utils.irp.f b/src/utils/cgtos_utils.irp.f index a820d5f2..9c25edc2 100644 --- a/src/utils/cgtos_utils.irp.f +++ b/src/utils/cgtos_utils.irp.f @@ -56,7 +56,7 @@ subroutine give_explicit_cpoly_and_cgaussian_x(P_new, P_center, p, fact_k, iorde call multiply_cpoly(P_a(0), a, P_b(0), b, P_new(0), n_new) iorder = a + b -end subroutine give_explicit_cpoly_and_cgaussian_x +end ! --- @@ -141,7 +141,7 @@ subroutine give_explicit_cpoly_and_cgaussian(P_new, P_center, p, fact_k, iorder, !DIR$ FORCEINLINE call multiply_cpoly(P_a(0,3), a(3), P_b(0,3), b(3), P_new(0,3), n_new) -end subroutine give_explicit_cpoly_and_cgaussian +end ! --- @@ -249,7 +249,7 @@ subroutine cgaussian_product(a, xa, b, xb, k, p, xp) xp(2) = ( a * xa(2) + b * xb(2) ) * p_inv xp(3) = ( a * xa(3) + b * xb(3) ) * p_inv -end subroutine cgaussian_product +end ! --- @@ -290,7 +290,7 @@ subroutine cgaussian_product_x(a, xa, b, xb, k, p, xp) k = zexp(-k) xp = (a*xa + b*xb) * p_inv -end subroutine cgaussian_product_x +end ! --- @@ -338,7 +338,7 @@ subroutine multiply_cpoly(b, nb, c, nc, d, nd) exit enddo -end subroutine multiply_cpoly +end ! --- @@ -373,7 +373,7 @@ subroutine add_cpoly(b, nb, c, nc, d, nd) if(nd < 0) exit enddo -end subroutine add_cpoly +end ! --- @@ -413,7 +413,7 @@ subroutine add_cpoly_multiply(b, nb, cst, d, nd) endif -end subroutine add_cpoly_multiply +end ! --- @@ -475,7 +475,7 @@ subroutine recentered_cpoly2(P_A, x_A, x_P, a, P_B, x_B, x_Q, b) P_B(i) = binom_func(b,b-i) * pows_b(b-i) enddo -end subroutine recentered_cpoly2 +end ! --- @@ -514,267 +514,7 @@ complex*16 function Fc_integral(n, inv_sq_p) Fc_integral = sqpi * 0.5d0**n * inv_sq_p**dble(n+1) * fact(n) / fact(shiftr(n, 1)) -end function Fc_integral - -! --- - -complex*16 function crint(n, rho) - - implicit none - include 'constants.include.F' - - integer, intent(in) :: n - complex*16, intent(in) :: rho - - integer :: i, mmax - double precision :: rho_mod, rho_re, rho_im - double precision :: sq_rho_re, sq_rho_im - double precision :: n_tmp - complex*16 :: sq_rho, rho_inv, rho_exp - - complex*16 :: crint_smallz, cpx_erf - - rho_re = REAL (rho) - rho_im = AIMAG(rho) - rho_mod = dsqrt(rho_re*rho_re + rho_im*rho_im) - - if(rho_mod < 10.d0) then - ! small z - - if(rho_mod .lt. 1.d-10) then - crint = 1.d0 / dble(n + n + 1) - else - crint = crint_smallz(n, rho) - endif - - else - ! large z - - if(rho_mod .gt. 40.d0) then - - n_tmp = dble(n) + 0.5d0 - crint = 0.5d0 * gamma(n_tmp) / (rho**n_tmp) - - else - - ! get \sqrt(rho) - sq_rho_re = sq_op5 * dsqrt(rho_re + rho_mod) - sq_rho_im = 0.5d0 * rho_im / sq_rho_re - sq_rho = sq_rho_re + (0.d0, 1.d0) * sq_rho_im - - rho_exp = 0.5d0 * zexp(-rho) - rho_inv = (1.d0, 0.d0) / rho - - crint = 0.5d0 * sqpi * cpx_erf(sq_rho_re, sq_rho_im) / sq_rho - mmax = n - if(mmax .gt. 0) then - do i = 0, mmax-1 - crint = ((dble(i) + 0.5d0) * crint - rho_exp) * rho_inv - enddo - endif - - ! *** - - endif - - endif - -! print *, n, real(rho), real(crint) - -end function crint - -! --- - -complex*16 function crint_sum(n_pt_out, rho, d1) - - implicit none - include 'constants.include.F' - - integer, intent(in) :: n_pt_out - complex*16, intent(in) :: rho, d1(0:n_pt_out) - - integer :: n, i, mmax - double precision :: rho_mod, rho_re, rho_im - double precision :: sq_rho_re, sq_rho_im - complex*16 :: sq_rho, F0 - complex*16 :: rho_tmp, rho_inv, rho_exp - complex*16, allocatable :: Fm(:) - - complex*16 :: crint_smallz, cpx_erf - - rho_re = REAL (rho) - rho_im = AIMAG(rho) - rho_mod = dsqrt(rho_re*rho_re + rho_im*rho_im) - - if(rho_mod < 10.d0) then - ! small z - - if(rho_mod .lt. 1.d-10) then - -! print *, ' 111' -! print *, ' rho = ', rho - - crint_sum = d1(0) -! print *, 0, 1 - - do i = 2, n_pt_out, 2 - - n = shiftr(i, 1) - crint_sum = crint_sum + d1(i) / dble(n+n+1) - -! print *, n, 1.d0 / dble(n+n+1) - enddo - - ! *** - - else - -! print *, ' 222' -! print *, ' rho = ', real(rho) -! if(abs(aimag(rho)) .gt. 1d-15) then -! print *, ' complex rho', rho -! stop -! endif - - crint_sum = d1(0) * crint_smallz(0, rho) - -! print *, 0, real(d1(0)), real(crint_smallz(0, rho)) -! if(abs(aimag(d1(0))) .gt. 1d-15) then -! print *, ' complex d1(0)', d1(0) -! stop -! endif - - do i = 2, n_pt_out, 2 - n = shiftr(i, 1) - crint_sum = crint_sum + d1(i) * crint_smallz(n, rho) - -! print *, n, real(d1(i)), real(crint_smallz(n, rho)) -! if(abs(aimag(d1(i))) .gt. 1d-15) then -! print *, ' complex d1(i)', i, d1(i) -! stop -! endif - - enddo - -! print *, 'sum = ', real(crint_sum) -! if(abs(aimag(crint_sum)) .gt. 1d-15) then -! print *, ' complex crint_sum', crint_sum -! stop -! endif - - ! *** - - endif - - else - ! large z - - if(rho_mod .gt. 40.d0) then - -! print *, ' 333' -! print *, ' rho = ', rho - - rho_inv = (1.d0, 0.d0) / rho - rho_tmp = 0.5d0 * sqpi * zsqrt(rho_inv) - crint_sum = rho_tmp * d1(0) -! print *, 0, rho_tmp - - do i = 2, n_pt_out, 2 - n = shiftr(i, 1) - rho_tmp = rho_tmp * (dble(n) + 0.5d0) * rho_inv - crint_sum = crint_sum + rho_tmp * d1(i) -! print *, n, rho_tmp - enddo - - ! *** - - else - -! print *, ' 444' -! print *, ' rho = ', rho - - ! get \sqrt(rho) - sq_rho_re = sq_op5 * dsqrt(rho_re + rho_mod) - sq_rho_im = 0.5d0 * rho_im / sq_rho_re - sq_rho = sq_rho_re + (0.d0, 1.d0) * sq_rho_im - !sq_rho = zsqrt(rho) - - - F0 = 0.5d0 * sqpi * cpx_erf(sq_rho_re, sq_rho_im) / sq_rho - crint_sum = F0 * d1(0) -! print *, 0, F0 - - rho_exp = 0.5d0 * zexp(-rho) - rho_inv = (1.d0, 0.d0) / rho - - mmax = shiftr(n_pt_out, 1) - if(mmax .gt. 0) then - - allocate( Fm(mmax) ) - Fm(1:mmax) = (0.d0, 0.d0) - - do n = 0, mmax-1 - F0 = ((dble(n) + 0.5d0) * F0 - rho_exp) * rho_inv - Fm(n+1) = F0 -! print *, n, F0 - enddo - - do i = 2, n_pt_out, 2 - n = shiftr(i, 1) - crint_sum = crint_sum + Fm(n) * d1(i) - enddo - deallocate(Fm) - endif - - ! *** - - endif - - endif - -end function crint_sum - -! --- - -complex*16 function crint_smallz(n, rho) - - BEGIN_DOC - ! Standard version of rint - END_DOC - - implicit none - integer, intent(in) :: n - complex*16, intent(in) :: rho - - integer, parameter :: kmax = 40 - double precision, parameter :: eps = 1.d-13 - - integer :: k - double precision :: delta_mod - complex*16 :: rho_k, ct, delta_k - - ct = 0.5d0 * zexp(-rho) * gamma(dble(n) + 0.5d0) - rho_k = (1.d0, 0.d0) - crint_smallz = ct * rho_k / gamma(dble(n) + 1.5d0) - - do k = 1, kmax - - rho_k = rho_k * rho - delta_k = ct * rho_k / gamma(dble(n+k) + 1.5d0) - crint_smallz = crint_smallz + delta_k - - delta_mod = dsqrt(REAL(delta_k)*REAL(delta_k) + AIMAG(delta_k)*AIMAG(delta_k)) - if(delta_mod .lt. eps) return - enddo - - if(delta_mod > eps) then - write(*,*) ' pb in crint_smallz !' - write(*,*) ' n, rho = ', n, rho - write(*,*) ' delta_mod = ', delta_mod - stop 1 - endif - -end function crint_smallz +end ! --- diff --git a/src/utils/constants.include.F b/src/utils/constants.include.F index 7b01f888..830b71a1 100644 --- a/src/utils/constants.include.F +++ b/src/utils/constants.include.F @@ -9,6 +9,9 @@ double precision, parameter :: pi_5_2 = 34.9868366552d0 double precision, parameter :: dfour_pi = 4.d0*dacos(-1.d0) double precision, parameter :: dtwo_pi = 2.d0*dacos(-1.d0) double precision, parameter :: inv_sq_pi = 1.d0/dsqrt(dacos(-1.d0)) +double precision, parameter :: c_mu_gauss = 27.d0/(8.d0*dsqrt(dacos(-1.d0))) +double precision, parameter :: c_mu_gauss_tot = 1.5d0*27.d0/(8.d0*dsqrt(dacos(-1.d0)))+3.d0/dsqrt(dacos(-1.d0)) +double precision, parameter :: alpha_mu_gauss = 1.5d0 double precision, parameter :: inv_sq_pi_2 = 0.5d0/dsqrt(dacos(-1.d0)) double precision, parameter :: thresh = 1.d-15 double precision, parameter :: cx_lda = -0.73855876638202234d0 diff --git a/src/utils/cpx_boys.irp.f b/src/utils/cpx_boys.irp.f new file mode 100644 index 00000000..9ffcc817 --- /dev/null +++ b/src/utils/cpx_boys.irp.f @@ -0,0 +1,543 @@ + +! --- + +complex*16 function crint_1(n, rho) + + implicit none + include 'constants.include.F' + + integer, intent(in) :: n + complex*16, intent(in) :: rho + + integer :: i, mmax + double precision :: rho_mod, rho_re, rho_im + double precision :: sq_rho_re, sq_rho_im + double precision :: n_tmp + complex*16 :: sq_rho, rho_inv, rho_exp + + complex*16 :: crint_smallz, cpx_erf_1 + complex*16 :: cpx_erf_2 + + rho_re = real (rho) + rho_im = aimag(rho) + rho_mod = dsqrt(rho_re*rho_re + rho_im*rho_im) + + if(rho_mod < 10.d0) then + ! small z + if(rho_mod .lt. 1.d-15) then + crint_1 = 1.d0 / dble(n + n + 1) + else + crint_1 = crint_smallz(n, rho) + endif + + else + ! large z + + if(rho_mod .gt. 40.d0) then + + n_tmp = dble(n) + 0.5d0 + crint_1 = 0.5d0 * gamma(n_tmp) / (rho**n_tmp) + + else + + ! get \sqrt(rho) + sq_rho_re = sq_op5 * dsqrt(rho_re + rho_mod) + sq_rho_im = 0.5d0 * rho_im / sq_rho_re + sq_rho = sq_rho_re + (0.d0, 1.d0) * sq_rho_im + + rho_exp = 0.5d0 * zexp(-rho) + rho_inv = (1.d0, 0.d0) / rho + + !print*, sq_rho_re, sq_rho_im + !print*, cpx_erf_1(sq_rho_re, sq_rho_im) + !print*, cpx_erf_2(sq_rho_re, sq_rho_im) + + crint_1 = 0.5d0 * sqpi * cpx_erf_1(sq_rho_re, sq_rho_im) / sq_rho + mmax = n + if(mmax .gt. 0) then + do i = 0, mmax-1 + crint_1 = ((dble(i) + 0.5d0) * crint_1 - rho_exp) * rho_inv + enddo + endif + + endif + + endif + +end + +! --- + +complex*16 function crint_quad(n, rho) + + implicit none + + integer, intent(in) :: n + complex*16, intent(in) :: rho + + integer :: i_quad, n_quad + double precision :: tmp_inv, tmp + + n_quad = 1000000000 + tmp_inv = 1.d0 / dble(n_quad) + + !crint_quad = 0.5d0 * zexp(-rho) + !do i_quad = 1, n_quad - 1 + ! tmp = tmp_inv * dble(i_quad) + ! tmp = tmp * tmp + ! crint_quad += zexp(-rho*tmp) * tmp**n + !enddo + !crint_quad = crint_quad * tmp_inv + + !crint_quad = 0.5d0 * zexp(-rho) + !do i_quad = 1, n_quad - 1 + ! tmp = tmp_inv * dble(i_quad) + ! crint_quad += zexp(-rho*tmp) * tmp**n / dsqrt(tmp) + !enddo + !crint_quad = crint_quad * 0.5d0 * tmp_inv + + ! Composite Boole's Rule + crint_quad = 7.d0 * zexp(-rho) + do i_quad = 1, n_quad - 1 + tmp = tmp_inv * dble(i_quad) + tmp = tmp * tmp + if(modulo(i_quad, 4) .eq. 0) then + crint_quad += 14.d0 * zexp(-rho*tmp) * tmp**n + else if(modulo(i_quad, 2) .eq. 0) then + crint_quad += 12.d0 * zexp(-rho*tmp) * tmp**n + else + crint_quad += 32.d0 * zexp(-rho*tmp) * tmp**n + endif + enddo + crint_quad = crint_quad * 2.d0 * tmp_inv / 45.d0 + + ! Composite Simpson's 3/8 rule + !crint_quad = zexp(-rho) + !do i_quad = 1, n_quad - 1 + ! tmp = tmp_inv * dble(i_quad) + ! tmp = tmp * tmp + ! if(modulo(i_quad, 3) .eq. 0) then + ! crint_quad += 2.d0 * zexp(-rho*tmp) * tmp**n + ! else + ! crint_quad += 3.d0 * zexp(-rho*tmp) * tmp**n + ! endif + !enddo + !crint_quad = crint_quad * 3.d0 * tmp_inv / 8.d0 + +end + +! --- + +complex*16 function crint_sum_1(n_pt_out, rho, d1) + + implicit none + include 'constants.include.F' + + integer, intent(in) :: n_pt_out + complex*16, intent(in) :: rho, d1(0:n_pt_out) + + integer :: n, i, mmax + double precision :: rho_mod, rho_re, rho_im + double precision :: sq_rho_re, sq_rho_im + complex*16 :: sq_rho, F0 + complex*16 :: rho_tmp, rho_inv, rho_exp + complex*16, allocatable :: Fm(:) + + complex*16 :: crint_smallz, cpx_erf_1 + + + rho_re = real (rho) + rho_im = aimag(rho) + rho_mod = dsqrt(rho_re*rho_re + rho_im*rho_im) + +! ! debug +! double precision :: d1_real(0:n_pt_out) +! double precision :: rint_sum +! do i = 0, n_pt_out +! d1_real(i) = real(d1(i)) +! enddo +! crint_sum_1 = rint_sum(n_pt_out, rho_re, d1_real) +! return + + if(rho_mod < 10.d0) then + ! small z + + if(rho_mod .lt. 1.d-15) then + + crint_sum_1 = d1(0) + do i = 2, n_pt_out, 2 + n = shiftr(i, 1) + crint_sum_1 = crint_sum_1 + d1(i) / dble(n+n+1) + enddo + + else + + crint_sum_1 = d1(0) * crint_smallz(0, rho) + do i = 2, n_pt_out, 2 + n = shiftr(i, 1) + crint_sum_1 = crint_sum_1 + d1(i) * crint_smallz(n, rho) + enddo + + endif + + else + ! large z + + if(rho_mod .gt. 40.d0) then + + rho_inv = (1.d0, 0.d0) / rho + rho_tmp = 0.5d0 * sqpi * zsqrt(rho_inv) + + crint_sum_1 = rho_tmp * d1(0) + do i = 2, n_pt_out, 2 + n = shiftr(i, 1) + rho_tmp = rho_tmp * (dble(n) + 0.5d0) * rho_inv + crint_sum_1 = crint_sum_1 + rho_tmp * d1(i) + enddo + + else + + ! get \sqrt(rho) + sq_rho_re = sq_op5 * dsqrt(rho_re + rho_mod) + sq_rho_im = 0.5d0 * rho_im / sq_rho_re + sq_rho = sq_rho_re + (0.d0, 1.d0) * sq_rho_im + + F0 = 0.5d0 * sqpi * cpx_erf_1(sq_rho_re, sq_rho_im) / sq_rho + crint_sum_1 = F0 * d1(0) + + rho_exp = 0.5d0 * zexp(-rho) + rho_inv = (1.d0, 0.d0) / rho + + mmax = shiftr(n_pt_out, 1) + if(mmax .gt. 0) then + + allocate(Fm(mmax)) + Fm(1:mmax) = (0.d0, 0.d0) + + do n = 0, mmax-1 + F0 = ((dble(n) + 0.5d0) * F0 - rho_exp) * rho_inv + Fm(n+1) = F0 + enddo + + do i = 2, n_pt_out, 2 + n = shiftr(i, 1) + crint_sum_1 = crint_sum_1 + Fm(n) * d1(i) + enddo + + deallocate(Fm) + endif + + endif ! rho_mod + endif ! rho_mod + +end + +! --- + +complex*16 function crint_smallz(n, rho) + + BEGIN_DOC + ! Standard version of rint + END_DOC + + implicit none + integer, intent(in) :: n + complex*16, intent(in) :: rho + + integer, parameter :: kmax = 40 + double precision, parameter :: eps = 1.d-13 + + integer :: k + double precision :: delta_mod + complex*16 :: rho_k, ct, delta_k + + ct = 0.5d0 * zexp(-rho) * gamma(dble(n) + 0.5d0) + crint_smallz = ct / gamma(dble(n) + 1.5d0) + + rho_k = (1.d0, 0.d0) + do k = 1, kmax + + rho_k = rho_k * rho + delta_k = ct * rho_k / gamma(dble(n+k) + 1.5d0) + crint_smallz = crint_smallz + delta_k + + delta_mod = dsqrt(real(delta_k)*real(delta_k) + aimag(delta_k)*aimag(delta_k)) + if(delta_mod .lt. eps) return + enddo + + if(delta_mod > eps) then + write(*,*) ' pb in crint_smallz !' + write(*,*) ' n, rho = ', n, rho + write(*,*) ' delta_mod = ', delta_mod + !stop 1 + endif + +end + +! --- + +complex*16 function crint_2(n, rho) + + implicit none + + integer, intent(in) :: n + complex*16, intent(in) :: rho + + double precision :: tmp + complex*16 :: rho2 + complex*16 :: vals(0:n) + complex*16, external :: crint_smallz + + if(abs(rho) < 10.d0) then + + if(abs(rho) .lt. 1d-6) then + tmp = 2.d0 * dble(n) + rho2 = rho * rho + crint_2 = 1.d0 / (tmp + 1.d0) & + - rho / (tmp + 3.d0) & + + 0.5d0 * rho2 / (tmp + 5.d0) & + - 0.16666666666666666d0 * rho * rho2 / (tmp + 7.d0) + else + crint_2 = crint_smallz(n, rho) + endif + + else + + if(real(rho) .ge. 0.d0) then + call zboysfun(n, rho, vals) + crint_2 = vals(n) + else + call zboysfunnrp(n, rho, vals) + crint_2 = vals(n) * zexp(-rho) + endif + + endif + + return +end + +! --- + +subroutine zboysfun(n_max, x, vals) + + BEGIN_DOC + ! + ! Computes values of the Boys function for n = 0, 1, ..., n_max + ! for a complex valued argument + ! + ! Input: x --- argument, complex*16, Re(x) >= 0 + ! Output: vals --- values of the Boys function, n = 0, 1, ..., n_max + ! + END_DOC + + implicit none + + integer, intent(in) :: n_max + complex*16, intent(in) :: x + complex*16, intent(out) :: vals(0:n_max) + + integer :: n + complex*16 :: yy, x_inv + + call zboysfun00(x, vals(0)) + + yy = 0.5d0 * zexp(-x) + x_inv = (1.d0, 0.d0) / x + do n = 1, n_max + vals(n) = ((dble(n) - 0.5d0) * vals(n-1) - yy) * x_inv + enddo + + return +end + +! --- + +subroutine zboysfunnrp(n_max, x, vals) + + BEGIN_DOC + ! + ! Computes values of e^z F(n,z) for n = 0, 1, ..., n_max + ! (where F(n,z) are the Boys functions) + ! for a complex valued argument WITH NEGATIVE REAL PART + ! + ! Input: x --- argument, complex *16 Re(x)<=0 + ! Output: vals --- values of e^z F(n,z), n = 0, 1, ..., n_max + ! + END_DOC + + implicit none + + integer, intent(in) :: n_max + complex*16, intent(in) :: x + complex*16, intent(out) :: vals(0:n_max) + + integer :: n + complex*16 :: x_inv + + call zboysfun00nrp(x, vals(0)) + + x_inv = (1.d0, 0.d0) / x + do n = 1, n_max + vals(n) = ((dble(n) - 0.5d0) * vals(n-1) - 0.5d0) * x_inv + enddo + + return +end + +! --- + +complex*16 function crint_sum_2(n_pt_out, rho, d1) + + implicit none + + integer, intent(in) :: n_pt_out + complex*16, intent(in) :: rho, d1(0:n_pt_out) + + integer :: n, i + integer :: n_max + + complex*16, allocatable :: vals(:) + + !complex*16, external :: crint_2 + !crint_sum_2 = (0.d0, 0.d0) + !do i = 0, n_pt_out, 2 + ! n = shiftr(i, 1) + ! crint_sum_2 = crint_sum_2 + d1(i) * crint_2(n, rho) + !enddo + + n_max = shiftr(n_pt_out, 1) + + allocate(vals(0:n_max)) + call crint_2_vec(n_max, rho, vals) + + crint_sum_2 = d1(0) * vals(0) + do i = 2, n_pt_out, 2 + n = shiftr(i, 1) + crint_sum_2 += d1(i) * vals(n) + enddo + + deallocate(vals) + + return +end + +! --- + +subroutine crint_2_vec(n_max, rho, vals) + + implicit none + + integer, intent(in) :: n_max + complex*16, intent(in) :: rho + complex*16, intent(out) :: vals(0:n_max) + + integer :: n + double precision :: tmp, abs_rho + complex*16 :: rho2, rho3, erho + + + abs_rho = abs(rho) + + if(abs_rho < 10.d0) then + + if(abs_rho .lt. 1d-6) then + + ! use finite expansion for very small rho + + ! rho^2 / 2 + rho2 = 0.5d0 * rho * rho + ! rho^3 / 6 + rho3 = 0.3333333333333333d0 * rho * rho2 + + vals(0) = 1.d0 - 0.3333333333333333d0 * rho + 0.2d0 * rho2 - 0.14285714285714285d0 * rho3 + do n = 1, n_max + tmp = 2.d0 * dble(n) + vals(n) = 1.d0 / (tmp + 1.d0) - rho / (tmp + 3.d0) & + + rho2 / (tmp + 5.d0) - rho3 / (tmp + 7.d0) + enddo + + else + + call crint_smallz_vec(n_max, rho, vals) + + endif + + else + + if(real(rho) .ge. 0.d0) then + + call zboysfun(n_max, rho, vals) + + else + + call zboysfunnrp(n_max, rho, vals) + erho = zexp(-rho) + do n = 0, n_max + vals(n) = vals(n) * erho + enddo + + endif + + endif + + return +end + +! --- + +subroutine crint_smallz_vec(n_max, rho, vals) + + BEGIN_DOC + ! Standard version of rint + END_DOC + + implicit none + integer, intent(in) :: n_max + complex*16, intent(in) :: rho + complex*16, intent(out) :: vals(0:n_max) + + integer, parameter :: kmax = 40 + double precision, parameter :: eps = 1.d-13 + + integer :: k, n + complex*16 :: ct, delta_k + complex*16 :: rhoe + complex*16, allocatable :: rho_k(:) + + + allocate(rho_k(0:kmax)) + + rho_k(0) = (1.d0, 0.d0) + do k = 1, kmax + rho_k(k) = rho_k(k-1) * rho + enddo + + rhoe = 0.5d0 * zexp(-rho) + + do n = 0, n_max + + ct = rhoe * gamma(dble(n) + 0.5d0) + vals(n) = ct / gamma(dble(n) + 1.5d0) + + do k = 1, kmax + delta_k = ct * rho_k(k) / gamma(dble(n+k) + 1.5d0) + vals(n) += delta_k + if(abs(delta_k) .lt. eps) then + exit + endif + enddo + + !if(abs(delta_k) > eps) then + ! write(*,*) ' pb in crint_smallz_vec !' + ! write(*,*) ' n, rho = ', n, rho + ! write(*,*) ' |delta_k| = ', abs(delta_k) + !endif + enddo + + deallocate(rho_k) + + return +end + +! --- + diff --git a/src/utils/cpx_erf.irp.f b/src/utils/cpx_erf.irp.f index 61f81055..1c5fa61d 100644 --- a/src/utils/cpx_erf.irp.f +++ b/src/utils/cpx_erf.irp.f @@ -1,7 +1,7 @@ ! --- -complex*16 function cpx_erf(x, y) +complex*16 function cpx_erf_1(x, y) BEGIN_DOC ! @@ -25,25 +25,25 @@ complex*16 function cpx_erf(x, y) if(yabs .lt. 1.d-15) then - cpx_erf = (1.d0, 0.d0) * derf(x) + cpx_erf_1 = (1.d0, 0.d0) * derf(x) return else erf_tmp1 = (1.d0, 0.d0) * derf(x) erf_tmp2 = erf_E(x, yabs) + erf_F(x, yabs) - erf_tmp3 = zexp(-(0.d0, 2.d0) * x * yabs) * ( erf_G(x, yabs) + erf_H(x, yabs) ) + erf_tmp3 = zexp(-(0.d0, 2.d0) * x * yabs) * (erf_G(x, yabs) + erf_H(x, yabs)) erf_tot = erf_tmp1 + erf_tmp2 - erf_tmp3 endif if(y .gt. 0.d0) then - cpx_erf = erf_tot + cpx_erf_1 = erf_tot else - cpx_erf = CONJG(erf_tot) + cpx_erf_1 = conjg(erf_tot) endif -end function cpx_erf +end ! --- @@ -54,7 +54,7 @@ complex*16 function erf_E(x, yabs) double precision, intent(in) :: x, yabs - if( (dabs(x).gt.6.d0) .or. (x==0.d0) ) then + if((dabs(x).gt.6.d0) .or. (x==0.d0)) then erf_E = (0.d0, 0.d0) return endif @@ -70,7 +70,7 @@ complex*16 function erf_E(x, yabs) endif -end function erf_E +end ! --- @@ -109,7 +109,7 @@ double precision function erf_F(x, yabs) endif -end function erf_F +end ! --- @@ -149,7 +149,7 @@ complex*16 function erf_G(x, yabs) enddo -end function erf_G +end ! --- @@ -172,7 +172,7 @@ complex*16 function erf_H(x, yabs) endif - if( (dabs(x) .lt. 10d0) .and. (yabs .lt. 6.1d0) ) then + if((dabs(x) .lt. 10d0) .and. (yabs .lt. 6.1d0)) then x2 = x * x ct = 0.5d0 * inv_pi @@ -186,7 +186,7 @@ complex*16 function erf_H(x, yabs) tmp2 = dexp(-tmp1-idble*yabs) * (x + (0.d0, 1.d0)*tmp0) / tmp1 erf_H = erf_H + tmp2 - tmp_mod = dsqrt(REAL(tmp2)*REAL(tmp2) + AIMAG(tmp2)*AIMAG(tmp2)) + tmp_mod = dsqrt(real(tmp2)*real(tmp2) + aimag(tmp2)*aimag(tmp2)) if(tmp_mod .lt. 1d-15) exit enddo erf_H = ct * erf_H @@ -197,8 +197,394 @@ complex*16 function erf_H(x, yabs) endif -end function erf_H +end ! --- +complex*16 function cpx_erf_2(x, y) + + BEGIN_DOC + ! + ! compute erf(z) for z = x + i y + ! + ! Beylkin & Sharma, J. Chem. Phys. 155, 174117 (2021) + ! https://doi.org/10.1063/5.0062444 + ! + END_DOC + + implicit none + + double precision, intent(in) :: x, y + + double precision :: yabs + complex*16 :: z + + yabs = dabs(y) + + if(yabs .lt. 1.d-15) then + + cpx_erf_2 = (1.d0, 0.d0) * derf(x) + return + + else + + z = x + (0.d0, 1.d0) * y + + if(x .ge. 0.d0) then + call zboysfun00(z, cpx_erf_2) + else + call zboysfun00nrp(z, cpx_erf_2) + cpx_erf_2 = cpx_erf_2 * zexp(-z) + endif + + endif + + return +end + +! --- + +subroutine zboysfun00(z, val) + + BEGIN_DOC + ! + ! Computes values of the Boys function for n=0 + ! for a complex valued argument + ! + ! Input: z --- argument, complex*16, Real(z) >= 0 + ! Output: val --- value of the Boys function n=0 + ! + ! Beylkin & Sharma, J. Chem. Phys. 155, 174117 (2021) + ! https://doi.org/10.1063/5.0062444 + ! + END_DOC + + implicit none + + double precision, parameter :: asymcoef(1:7) = (/ -0.499999999999999799d0, & + 0.249999999999993161d0, & + -0.374999999999766599d0, & + 0.937499999992027020d0, & + -3.28124999972738868d0, & + 14.7656249906697030d0, & + -81.2109371803307752d0 /) + + double precision, parameter :: taylcoef(0:10) = (/ 1.0d0, & + -0.333333333333333333d0, & + 0.1d0, & + -0.238095238095238095d-01, & + 0.462962962962962963d-02, & + -0.757575757575757576d-03, & + 0.106837606837606838d-03, & + -0.132275132275132275d-04, & + 1.458916900093370682d-06, & + -1.450385222315046877d-07, & + 1.3122532963802805073d-08 /) + + double precision, parameter :: sqpio2 = 0.886226925452758014d0 + + double precision, parameter :: pp(1:22) = (/ 0.001477878263796956477d0, & + 0.013317276413725817441d0, & + 0.037063591452052541530d0, & + 0.072752512422882761543d0, & + 0.120236941228785688896d0, & + 0.179574293958937717967d0, & + 0.253534046984087292596d0, & + 0.350388652780721927513d0, & + 0.482109575931276669313d0, & + 0.663028993158374107103d0, & + 0.911814736856590885929d0, & + 1.2539502287919293d0, & + 1.7244634233573395d0, & + 2.3715248262781863d0, & + 3.2613796996078355d0, & + 4.485130169059591d0, & + 6.168062135122484d0, & + 8.48247187231787d0, & + 11.665305486296793d0, & + 16.042417132288328d0, & + 22.06192951814709d0, & + 30.340112094708307d0 /) + + double precision, parameter :: ff(1:22) = (/ 0.0866431027201416556d0, & + 0.0857720608434394764d0, & + 0.0839350436829178814d0, & + 0.0809661970413229146d0, & + 0.0769089548492978618d0, & + 0.0731552078711821626d0, & + 0.0726950035163157228d0, & + 0.0752842556089304050d0, & + 0.0770943953645196145d0, & + 0.0754250625677530441d0, & + 0.0689686192650315305d0, & + 0.05744480422143023d0, & + 0.04208199434694545d0, & + 0.025838539448223282d0, & + 0.012445024157255563d0, & + 0.004292541592599837d0, & + 0.0009354342987735969d0, & + 0.10840885466502504d-03, & + 5.271867966761674d-06, & + 7.765974039750418d-08, & + 2.2138172422680093d-10, & + 6.594161760037707d-14 /) + + complex*16, intent(in) :: z + complex*16, intent(out) :: val + + integer :: k + complex*16 :: z1, zz, y + + zz = zexp(-z) + + if(abs(z) .ge. 100.0d0) then + + ! large |z| + z1 = 1.0d0 / zsqrt(z) + y = 1.0d0 / z + val = asymcoef(7) + do k = 6, 1, -1 + val = val * y + asymcoef(k) + enddo + val = zz * val * y + z1 * sqpio2 + + else if(abs(z) .le. 0.35d0) then + + ! small |z| + val = taylcoef(10) * (1.d0, 0.d0) + do k = 9, 0, -1 + val = val * z + taylcoef(k) + enddo + + else + + ! intermediate |z| + val = sqpio2 / zsqrt(z) - 0.5d0 * zz * sum(ff(1:22)/(z+pp(1:22))) + !val = (0.d0, 0.d0) + !do k = 1, 22 + ! val += ff(k) / (z + pp(k)) + !enddo + !val = sqpio2 / zsqrt(z) - 0.5d0 * zz * val + + endif + + return +end + +! --- + +subroutine zboysfun00nrp(z, val) + + BEGIN_DOC + ! + ! Computes values of the exp(z) F(0,z) + ! (where F(0,z) is the Boys function) + ! for a complex valued argument with Real(z)<=0 + ! + ! Input: z --- argument, complex*16, !!! Real(z)<=0 !!! + ! Output: val --- value of the function !!! exp(z) F(0,z) !!!, where F(0,z) is the Boys function + ! + ! Beylkin & Sharma, J. Chem. Phys. 155, 174117 (2021) + ! https://doi.org/10.1063/5.0062444 + ! + END_DOC + + implicit none + + double precision, parameter :: asymcoef(1:7) = (/ -0.499999999999999799d0, & + 0.249999999999993161d0, & + -0.374999999999766599d0, & + 0.937499999992027020d0, & + -3.28124999972738868d0, & + 14.7656249906697030d0, & + -81.2109371803307752d0 /) + + double precision, parameter :: taylcoef(0:10) = (/ 1.0d0, & + -0.333333333333333333d0, & + 0.1d0, & + -0.238095238095238095d-01, & + 0.462962962962962963d-02, & + -0.757575757575757576d-03, & + 0.106837606837606838d-03, & + -0.132275132275132275d-04, & + 1.458916900093370682d-06, & + -1.450385222315046877d-07, & + 1.3122532963802805073d-08 /) + + double precision, parameter :: tol = 1.0d-03 + double precision, parameter :: sqpio2 = 0.886226925452758014d0 ! sqrt(pi)/2 + double precision, parameter :: pi = 3.14159265358979324d0 + double precision, parameter :: etmax = 25.7903399171930621d0 + double precision, parameter :: etmax1 = 26.7903399171930621d0 + complex*16, parameter :: ima = (0.d0, 1.d0) + + double precision, parameter :: pp(1:16) = (/ 0.005299532504175031d0, & + 0.0277124884633837d0, & + 0.06718439880608407d0, & + 0.12229779582249845d0, & + 0.19106187779867811d0, & + 0.27099161117138637d0, & + 0.35919822461037054d0, & + 0.45249374508118123d0, & + 0.5475062549188188d0, & + 0.6408017753896295d0, & + 0.7290083888286136d0, & + 0.8089381222013219d0, & + 0.8777022041775016d0, & + 0.9328156011939159d0, & + 0.9722875115366163d0, & + 0.994700467495825d0 /) + + double precision, parameter :: ww(1:16) = (/ 0.013576229705876844d0, & + 0.03112676196932382d0, & + 0.04757925584124612d0, & + 0.062314485627766904d0, & + 0.07479799440828848d0, & + 0.08457825969750153d0, & + 0.09130170752246194d0, & + 0.0947253052275344d0, & + 0.0947253052275344d0, & + 0.09130170752246194d0, & + 0.08457825969750153d0, & + 0.07479799440828848d0, & + 0.062314485627766904d0, & + 0.04757925584124612d0, & + 0.03112676196932382d0, & + 0.013576229705876844d0 /) + + double precision, parameter :: qq (1:16) = (/ 0.0007243228510223928d0, & + 0.01980651726441906d0, & + 0.11641097769229371d0, & + 0.38573968881461146d0, & + 0.9414671037609641d0, & + 1.8939510935716377d0, & + 3.3275564293459383d0, & + 5.280587297262129d0, & + 7.730992222360452d0, & + 10.590207725831563d0, & + 13.706359477128965d0, & + 16.876705473663804d0, & + 19.867876155236257d0, & + 22.441333930203022d0, & + 24.380717439613566d0, & + 25.51771075067431d0 /) + + + double precision, parameter :: qq1 (1:16) = (/ 0.0007524078957852004d0,& + 0.020574499281252233d0, & + 0.12092472113522865d0, & + 0.40069643967765295d0, & + 0.9779717449089211d0, & + 1.9673875468969015d0, & + 3.4565797939091802d0, & + 5.485337886599723d0, & + 8.030755321535683d0, & + 11.000834641174064d0, & + 14.237812708111456d0, & + 17.531086359214406d0, & + 20.6382373144543d0, & + 23.31147887603379d0, & + 25.326060444703632d0, & + 26.507139770710722d0 /) + + double precision, parameter :: uu(1:16) = (/ 0.9992759394074501d0, & + 0.9803883431758104d0, & + 0.8901093330366746d0, & + 0.6799475005849274d0, & + 0.3900551639790145d0, & + 0.15047608763371934d0, & + 0.0358806749968974d0, & + 0.005089440900100864d0, & + 0.00043900830706867264d0, & + 0.000025161192619824898d0, & + 1.1153308427285078d-6, & + 4.68317018372038d-8, & + 2.3522908467181876d-9, & + 1.7941242138648815d-10, & + 2.5798173021885247d-11, & + 8.27559122014575d-12 /) + + + double precision, parameter :: uu1(1:16) = (/ 0.999247875092057d0, & + 0.979635711599488d0, & + 0.8861006617341018d0, & + 0.6698533710831932d0, & + 0.3760730980014839d0, & + 0.13982165701683388d0, & + 0.031537442321301304d0, & + 0.004147133581658446d0, & + 0.0003253024081883165d0, & + 0.000016687766678889653d0, & + 6.555359391864376d-7, & + 2.4341421258295026d-8, & + 1.0887481200652014d-9, & + 7.51542178140961d-11, & + 1.002378402152542d-11, & + 3.0767730761654096d-12 /) + + complex*16, intent(in) :: z + complex*16, intent(out) :: val + + integer :: k + complex*16 :: z1, zz, y, zsum, tmp, zt, q, p + + zz = zexp(z) + + if(abs(z) .ge. 100.0d0) then + ! large |z| + z1 = 1.0d0 / zsqrt(z) + y = 1.0d0 / z + val = asymcoef(7) + do k = 6, 1, -1 + val = val * y + asymcoef(k) + enddo + val = val * y + z1 * sqpio2 * zz + return + endif + + if(abs(z) .le. 0.35d0) then + ! small |z| + val = taylcoef(10) * (1.d0, 0.d0) + do k = 9, 0, -1 + val = val * z + taylcoef(k) + enddo + val = val * zz + return + endif + + if(abs(etmax+z) .ge. 0.5d0) then + ! intermediate |z| + zsum = (0.d0, 0.d0) + do k = 1, 16 + if(abs(z + qq(k)) .ge. tol) then + zsum = zsum + ww(k) * (zz - uu(k)) / (qq(k) + z) + else + q = z + qq(k) + p = 1.0d0 - 0.5d0*q + q*q/6.0d0 - q*q*q/24.0d0 + q*q*q*q/120.0d0 + zsum = zsum + ww(k) * p *zz + endif + enddo + zt = ima * sqrt(z / etmax) + tmp = 0.5d0 * ima * log((1.0d0 - zt) / (1.0d0 + zt)) + val = sqrt(etmax) * zsum / sqrt(pi) + zz * tmp / sqrt(pi*z) + else + zsum = (0.d0, 0.d0) + do k = 1, 16 + if(abs(z + qq1(k)) .ge. tol) then + zsum = zsum + ww(k) * (zz - uu1(k)) / (qq1(k) + z) + else + q = z + qq1(k) + p = 1.0d0 - 0.5d0*q + q*q/6.0d0 - q*q*q/24.0d0 + q*q*q*q/120.0d0 + zsum = zsum + ww(k) * p * zz + endif + enddo + zt = ima * zsqrt(z / etmax1) + tmp = 0.5d0 * ima * log((1.0d0 - zt) / (1.0d0 + zt)) + val = dsqrt(etmax1) * zsum / dsqrt(pi) + zz * tmp / zsqrt(pi*z) + endif + + return +end + +! --- diff --git a/src/utils/mmap.f90 b/src/utils/mmap.f90 index e342b422..4ac32233 100644 --- a/src/utils/mmap.f90 +++ b/src/utils/mmap.f90 @@ -2,6 +2,34 @@ module mmap_module use iso_c_binding + type mmap_type + type(c_ptr) :: ptr ! Pointer to the data + character*(128) :: filename ! Name of the file + integer*8 :: length ! Size of the array in bytes + integer :: fd ! File descriptor + + ! Pointers to data + integer, pointer :: i1(:) + integer, pointer :: i2(:,:) + integer, pointer :: i3(:,:,:) + integer, pointer :: i4(:,:,:,:) + + integer*8, pointer :: i81(:) + integer*8, pointer :: i82(:,:) + integer*8, pointer :: i83(:,:,:) + integer*8, pointer :: i84(:,:,:,:) + + double precision, pointer :: d1(:) + double precision, pointer :: d2(:,:) + double precision, pointer :: d3(:,:,:) + double precision, pointer :: d4(:,:,:,:) + + real, pointer :: s1(:) + real, pointer :: s2(:,:) + real, pointer :: s3(:,:,:) + real, pointer :: s4(:,:,:,:) + end type mmap_type + interface ! File descriptors @@ -82,7 +110,7 @@ module mmap_module length = length * shape(i) enddo fd_ = fd - call c_munmap_fortran( length, fd_, map) + call c_munmap_fortran(length, fd_, map) end subroutine subroutine msync(shape, bytes, fd, map) @@ -106,6 +134,200 @@ module mmap_module call c_msync_fortran( length, fd_, map) end subroutine + + ! Functions for the mmap_type + + subroutine mmap_create(filename, shape, bytes, read_only, single_node, map) + implicit none + character*(*), intent(in) :: filename ! Name of the mapped file + integer*8, intent(in) :: shape(:) ! Shape of the array to map + integer, intent(in) :: bytes ! Number of bytes per element + logical, intent(in) :: read_only ! If true, mmap is read-only + logical, intent(in) :: single_node! If true, mmap is on a single node + type(mmap_type), intent(out) :: map ! mmap + + integer :: i + logical :: temporary + + temporary = ( trim(filename) == '' ) + + if (.not.temporary) then + map%filename = filename + else + call getenv('EZFIO_FILE', map%filename) + map%filename = trim(map%filename) // '/work/tmpfile' + endif + + map%length = int(bytes,8) + do i=1,size(shape) + map%length = map%length * shape(i) + enddo + call mmap(map%filename, & + shape, & + bytes, & + map%fd, & + read_only, & + single_node, & + map%ptr) + + if (temporary) then + ! Deleting the file while it is open makes the file invisible on the filesystem, + ! and automatically deleted, even if the program crashes + open(UNIT=47, FILE=trim(map%filename), STATUS='OLD') + close(47,STATUS='DELETE') + endif + + map%d1 => NULL() + map%d2 => NULL() + map%d3 => NULL() + map%d4 => NULL() + map%s1 => NULL() + map%s2 => NULL() + map%s3 => NULL() + map%s4 => NULL() + map%i1 => NULL() + map%i2 => NULL() + map%i3 => NULL() + map%i4 => NULL() + map%i81 => NULL() + map%i82 => NULL() + map%i83 => NULL() + map%i84 => NULL() + + end + + subroutine mmap_create_d(filename, shape, read_only, single_node, map) + implicit none + character*(*), intent(in) :: filename ! Name of the mapped file + integer*8, intent(in) :: shape(:) ! Shape of the array to map + logical, intent(in) :: read_only ! If true, mmap is read-only + logical, intent(in) :: single_node! If true, mmap is on a single node + type(mmap_type), intent(out) :: map ! mmap + + call mmap_create(filename, shape, 8, read_only, single_node, map) + + select case (size(shape)) + case (1) + call c_f_pointer(map%ptr, map%d1, shape) + case (2) + call c_f_pointer(map%ptr, map%d2, shape) + case (3) + call c_f_pointer(map%ptr, map%d3, shape) + case (4) + call c_f_pointer(map%ptr, map%d4, shape) + case default + stop 'mmap: dimension not implemented' + end select + end subroutine + + subroutine mmap_create_s(filename, shape, read_only, single_node, map) + implicit none + character*(*), intent(in) :: filename ! Name of the mapped file + integer*8, intent(in) :: shape(:) ! Shape of the array to map + logical, intent(in) :: read_only ! If true, mmap is read-only + logical, intent(in) :: single_node! If true, mmap is on a single node + type(mmap_type), intent(out) :: map ! mmap + + call mmap_create(filename, shape, 4, read_only, single_node, map) + + select case (size(shape)) + case (1) + call c_f_pointer(map%ptr, map%s1, shape) + case (2) + call c_f_pointer(map%ptr, map%s2, shape) + case (3) + call c_f_pointer(map%ptr, map%s3, shape) + case (4) + call c_f_pointer(map%ptr, map%s4, shape) + case default + stop 'mmap: dimension not implemented' + end select + end subroutine + + subroutine mmap_create_i(filename, shape, read_only, single_node, map) + implicit none + character*(*), intent(in) :: filename ! Name of the mapped file + integer*8, intent(in) :: shape(:) ! Shape of the array to map + logical, intent(in) :: read_only ! If true, mmap is read-only + logical, intent(in) :: single_node! If true, mmap is on a single node + type(mmap_type), intent(out) :: map ! mmap + + call mmap_create(filename, shape, 4, read_only, single_node, map) + + select case (size(shape)) + case (1) + call c_f_pointer(map%ptr, map%i1, shape) + case (2) + call c_f_pointer(map%ptr, map%i2, shape) + case (3) + call c_f_pointer(map%ptr, map%i3, shape) + case (4) + call c_f_pointer(map%ptr, map%i4, shape) + case default + stop 'mmap: dimension not implemented' + end select + end subroutine + + subroutine mmap_create_i8(filename, shape, read_only, single_node, map) + implicit none + character*(*), intent(in) :: filename ! Name of the mapped file + integer*8, intent(in) :: shape(:) ! Shape of the array to map + logical, intent(in) :: read_only ! If true, mmap is read-only + logical, intent(in) :: single_node! If true, mmap is on a single node + type(mmap_type), intent(out) :: map ! mmap + + call mmap_create(filename, shape, 8, read_only, single_node, map) + + select case (size(shape)) + case (1) + call c_f_pointer(map%ptr, map%i81, shape) + case (2) + call c_f_pointer(map%ptr, map%i82, shape) + case (3) + call c_f_pointer(map%ptr, map%i83, shape) + case (4) + call c_f_pointer(map%ptr, map%i84, shape) + case default + stop 'mmap: dimension not implemented' + end select + end subroutine + + subroutine mmap_destroy(map) + implicit none + type(mmap_type), intent(inout) :: map + + call c_munmap_fortran(map%length, map%fd, map%ptr) + + map%ptr = C_NULL_PTR + map%filename = '' + map%length = 0 + map%fd = 0 + map%s1 => NULL() + map%s2 => NULL() + map%s3 => NULL() + map%s4 => NULL() + map%d1 => NULL() + map%d2 => NULL() + map%d3 => NULL() + map%d4 => NULL() + map%i1 => NULL() + map%i2 => NULL() + map%i3 => NULL() + map%i4 => NULL() + map%i81 => NULL() + map%i82 => NULL() + map%i83 => NULL() + map%i84 => NULL() + end subroutine + + + subroutine mmap_sync(map) + implicit none + type(mmap_type), intent(inout) :: map + + call c_msync_fortran(map%length, map%fd, map%ptr) + end subroutine + end module mmap_module diff --git a/src/utils_cc/diis.irp.f b/src/utils_cc/diis.irp.f index fe771373..a64a454f 100644 --- a/src/utils_cc/diis.irp.f +++ b/src/utils_cc/diis.irp.f @@ -53,10 +53,10 @@ subroutine diis_cc(all_err,all_t,sze,m,iter,t) !$OMP END PARALLEL do i = 1, m_iter - B(i,m_iter+1) = -1 + B(i,m_iter+1) = -1.d0 enddo do j = 1, m_iter - B(m_iter+1,j) = -1 + B(m_iter+1,j) = -1.d0 enddo ! Debug !print*,'B' @@ -493,7 +493,7 @@ subroutine update_t_ccsd_diis_v3(nO,nV,nb_iter,f_o,f_v,r1,r2,t1,t2,all_err,all_t do i = 1, nO*nV tmp(i) = t1(i) enddo - !$OMP END DO NOWAIT + !$OMP END DO !$OMP DO do i = 1, nO*nO*nV*nV tmp(i+nO*nV) = t2(i) @@ -515,7 +515,7 @@ subroutine update_t_ccsd_diis_v3(nO,nV,nb_iter,f_o,f_v,r1,r2,t1,t2,all_err,all_t do i = 1, nO*nV t1(i) = tmp(i) enddo - !$OMP END DO NOWAIT + !$OMP END DO !$OMP DO do i = 1, nO*nO*nV*nV t2(i) = tmp(i+nO*nV)