From 602e9e6fe751a38fd9bb66f5fa741aedfacd4e3d Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 11 May 2020 11:17:03 +0200 Subject: [PATCH] Working on normf --- ocaml/Input_ao_basis.ml | 67 ++++++++++++++++++++++++++++++++------- src/ao_basis/EZFIO.cfg | 12 +++++++ src/ao_basis/aos.irp.f | 37 ++++++++++++++------- src/cipsi/selection.irp.f | 13 ++++---- 4 files changed, 100 insertions(+), 29 deletions(-) diff --git a/ocaml/Input_ao_basis.ml b/ocaml/Input_ao_basis.ml index b0a66b75..f0b3e92c 100644 --- a/ocaml/Input_ao_basis.ml +++ b/ocaml/Input_ao_basis.ml @@ -13,6 +13,8 @@ module Ao_basis : sig ao_coef : AO_coef.t array; ao_expo : AO_expo.t array; ao_cartesian : bool; + ao_normalized : bool; + primitives_normalized : bool; } [@@deriving sexp] ;; val read : unit -> t option @@ -34,6 +36,8 @@ end = struct ao_coef : AO_coef.t array; ao_expo : AO_expo.t array; ao_cartesian : bool; + ao_normalized : bool; + primitives_normalized : bool; } [@@deriving sexp] ;; @@ -107,6 +111,24 @@ end = struct Ezfio.get_ao_basis_ao_cartesian () ;; + let read_ao_normalized () = + if not (Ezfio.has_ao_basis_ao_normalized()) then + get_default "ao_normalized" + |> bool_of_string + |> Ezfio.set_ao_basis_ao_normalized + ; + Ezfio.get_ao_basis_ao_normalized () + ;; + + let read_primitives_normalized () = + if not (Ezfio.has_ao_basis_primitives_normalized()) then + get_default "primitives_normalized" + |> bool_of_string + |> Ezfio.set_ao_basis_primitives_normalized + ; + Ezfio.get_ao_basis_primitives_normalized () + ;; + let to_long_basis b = let ao_num = AO_number.to_int b.ao_num in let gto_array = Array.init (AO_number.to_int b.ao_num) @@ -169,6 +191,8 @@ end = struct ao_coef ; ao_expo ; ao_cartesian ; + ao_normalized ; + primitives_normalized ; } = b in write_md5 b ; @@ -201,6 +225,8 @@ end = struct ~rank:2 ~dim:[| ao_num ; 3 |] ~data:ao_power) ; Ezfio.set_ao_basis_ao_cartesian(ao_cartesian); + Ezfio.set_ao_basis_ao_normalized(ao_normalized); + Ezfio.set_ao_basis_primitives_normalized(primitives_normalized); let ao_coef = Array.to_list ao_coef @@ -233,6 +259,8 @@ end = struct ao_coef = read_ao_coef () ; ao_expo = read_ao_expo () ; ao_cartesian = read_ao_cartesian () ; + ao_normalized = read_ao_normalized () ; + primitives_normalized = read_primitives_normalized () ; } in to_md5 result @@ -340,7 +368,10 @@ end = struct in { ao_basis = name ; ao_num ; ao_prim_num ; ao_prim_num_max ; ao_nucl ; - ao_power ; ao_coef ; ao_expo ; ao_cartesian } + ao_power ; ao_coef ; ao_expo ; ao_cartesian ; + ao_normalized = bool_of_string @@ get_default "ao_normalized"; + primitives_normalized = bool_of_string @@ get_default "primitives_normalized"; + } ;; let reorder b = @@ -394,6 +425,14 @@ Cartesian coordinates (6d,10f,...) :: ao_cartesian = %s +Use normalized primitive functions :: + + primitives_normalized = %s + +Use normalized basis functions :: + + ao_normalized = %s + Basis set (read-only) :: %s @@ -407,6 +446,8 @@ Basis set (read-only) :: " (AO_basis_name.to_string b.ao_basis) (string_of_bool b.ao_cartesian) + (string_of_bool b.primitives_normalized) + (string_of_bool b.ao_normalized) (Basis.to_string short_basis |> String_ext.split ~on:'\n' |> List.map (fun x-> " "^x) @@ -434,16 +475,18 @@ Basis set (read-only) :: let to_string b = Printf.sprintf " -ao_basis = %s -ao_num = %s -ao_prim_num = %s -ao_prim_num_max = %s -ao_nucl = %s -ao_power = %s -ao_coef = %s -ao_expo = %s -ao_cartesian = %s -md5 = %s +ao_basis = %s +ao_num = %s +ao_prim_num = %s +ao_prim_num_max = %s +ao_nucl = %s +ao_power = %s +ao_coef = %s +ao_expo = %s +ao_cartesian = %s +ao_normalized = %s +primitives_normalized = %s +md5 = %s " (AO_basis_name.to_string b.ao_basis) (AO_number.to_string b.ao_num) @@ -459,6 +502,8 @@ md5 = %s (b.ao_expo |> Array.to_list |> List.map AO_expo.to_string |> String.concat ", ") (b.ao_cartesian |> string_of_bool) + (b.ao_normalized |> string_of_bool) + (b.primitives_normalized |> string_of_bool) (to_md5 b |> MD5.to_string ) ;; diff --git a/src/ao_basis/EZFIO.cfg b/src/ao_basis/EZFIO.cfg index c3e2761b..51d726da 100644 --- a/src/ao_basis/EZFIO.cfg +++ b/src/ao_basis/EZFIO.cfg @@ -55,3 +55,15 @@ doc: If |true|, use |AOs| in Cartesian coordinates (6d,10f,...) interface: ezfio, provider default: false +[ao_normalized] +type: logical +doc: Use normalized basis functions +interface: ezfio, provider +default: true + +[primitives_normalized] +type: logical +doc: Use normalized primitive functions +interface: ezfio, provider +default: true + diff --git a/src/ao_basis/aos.irp.f b/src/ao_basis/aos.irp.f index f2bf16f0..e95ac711 100644 --- a/src/ao_basis/aos.irp.f +++ b/src/ao_basis/aos.irp.f @@ -20,25 +20,38 @@ END_PROVIDER C_A(2) = 0.d0 C_A(3) = 0.d0 ao_coef_normalized = 0.d0 + do i=1,ao_num powA(1) = ao_power(i,1) powA(2) = ao_power(i,2) powA(3) = ao_power(i,3) - do j=1,ao_prim_num(i) - call overlap_gaussian_xyz(C_A,C_A,ao_expo(i,j),ao_expo(i,j),powA,powA,overlap_x,overlap_y,overlap_z,norm,nz) - ao_coef_normalized(i,j) = ao_coef(i,j)/sqrt(norm) - enddo + ! Normalization of the primitives + if (primitives_normalized) then + do j=1,ao_prim_num(i) + call overlap_gaussian_xyz(C_A,C_A,ao_expo(i,j),ao_expo(i,j),powA,powA,overlap_x,overlap_y,overlap_z,norm,nz) + ao_coef_normalized(i,j) = ao_coef(i,j)/sqrt(norm) + enddo + else + do j=1,ao_prim_num(i) + ao_coef_normalized(i,j) = ao_coef(i,j) + enddo + endif + ! Normalization of the contracted basis functions - norm = 0.d0 - do j=1,ao_prim_num(i) - do k=1,ao_prim_num(i) - call overlap_gaussian_xyz(C_A,C_A,ao_expo(i,j),ao_expo(i,k),powA,powA,overlap_x,overlap_y,overlap_z,c,nz) - norm = norm+c*ao_coef_normalized(i,j)*ao_coef_normalized(i,k) - enddo - enddo - ao_coef_normalization_factor(i) = 1.d0/sqrt(norm) + if (ao_normalized) then + norm = 0.d0 + do j=1,ao_prim_num(i) + do k=1,ao_prim_num(i) + call overlap_gaussian_xyz(C_A,C_A,ao_expo(i,j),ao_expo(i,k),powA,powA,overlap_x,overlap_y,overlap_z,c,nz) + norm = norm+c*ao_coef_normalized(i,j)*ao_coef_normalized(i,k) + enddo + enddo + ao_coef_normalization_factor(i) = 1.d0/sqrt(norm) + else + ao_coef_normalization_factor(i) = 1.d0 + endif enddo END_PROVIDER diff --git a/src/cipsi/selection.irp.f b/src/cipsi/selection.irp.f index 5237ab94..cf5b9965 100644 --- a/src/cipsi/selection.irp.f +++ b/src/cipsi/selection.irp.f @@ -22,14 +22,14 @@ END_PROVIDER subroutine update_pt2_and_variance_weights(pt2, variance, norm, N_st) implicit none BEGIN_DOC -! Updates the rPT2- and Variance- matching weights. +! Updates the PT2- and Variance- matching weights. END_DOC integer, intent(in) :: N_st double precision, intent(in) :: pt2(N_st) double precision, intent(in) :: variance(N_st) double precision, intent(in) :: norm(N_st) - double precision :: avg, rpt2(N_st), element, dt, x + double precision :: avg, pt2_rpt2(N_st), element, dt, x integer :: k integer, save :: i_iter=0 integer, parameter :: i_itermax = 1 @@ -49,12 +49,13 @@ subroutine update_pt2_and_variance_weights(pt2, variance, norm, N_st) dt = 1.0d0 do k=1,N_st - rpt2(k) = pt2(k)/(1.d0 + norm(k)) + ! PT2 + rPT2 + pt2_rpt2(k) = pt2(k)* (1.d0 + 1.d0/(1.d0 + norm(k))) enddo - avg = sum(rpt2(1:N_st)) / dble(N_st) - 1.d-32 ! Avoid future division by zero + avg = sum(pt2_rpt2(1:N_st)) / dble(N_st) - 1.d-32 ! Avoid future division by zero do k=1,N_st - element = exp(dt*(rpt2(k)/avg -1.d0)) + element = exp(dt*(pt2_rpt2(k)/avg -1.d0)) element = min(2.0d0 , element) element = max(0.5d0 , element) memo_pt2(k,i_iter) = element @@ -71,7 +72,7 @@ subroutine update_pt2_and_variance_weights(pt2, variance, norm, N_st) enddo threshold_davidson_pt2 = min(1.d-6, & - max(threshold_davidson, 1.e-1 * PT2_relative_error * minval(abs(rpt2(1:N_states)))) ) + max(threshold_davidson, 1.e-1 * PT2_relative_error * minval(abs(pt2(1:N_states)))) ) SOFT_TOUCH pt2_match_weight variance_match_weight threshold_davidson_pt2 end