diff --git a/plugins/local/bi_ort_ints/no_dressing.irp.f b/plugins/local/bi_ort_ints/no_dressing.irp.f index 721ac0f8..fd2c6285 100644 --- a/plugins/local/bi_ort_ints/no_dressing.irp.f +++ b/plugins/local/bi_ort_ints/no_dressing.irp.f @@ -336,9 +336,6 @@ BEGIN_PROVIDER [double precision, noL_0e] double precision, allocatable :: tmp_M(:,:), tmp_S(:), tmp_O(:), tmp_J(:,:) double precision, allocatable :: tmp_M_priv(:,:), tmp_S_priv(:), tmp_O_priv(:), tmp_J_priv(:,:) - call wall_time(t0) - print*, " Providing noL_0e ..." - if(elec_alpha_num .eq. elec_beta_num) then allocate(tmp(elec_beta_num)) @@ -713,11 +710,6 @@ BEGIN_PROVIDER [double precision, noL_0e] endif - call wall_time(t1) - print*, " Wall time for noL_0e (min) = ", (t1 - t0)/60.d0 - - print*, " noL_0e = ", noL_0e - END_PROVIDER ! --- 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 33563102..db06e835 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 @@ -4,7 +4,7 @@ subroutine get_grad1_u12_withsq_r1_seq(ipoint, n_grid2, resx, resy, resz, res) BEGIN_DOC - ! + ! ! grad_1 u(r1,r2) ! ! we use grid for r1 and extra_grid for r2 @@ -167,7 +167,7 @@ subroutine grad1_j12_r1_seq(r1, n_grid2, gradx, grady, gradz) integer :: jpoint integer :: i_nucl, p, mpA, npA, opA double precision :: r2(3) - double precision :: dx, dy, dz, r12, tmp, r12_inv + double precision :: dx, dy, dz, r12, tmp double precision :: mu_val, mu_tmp, mu_der(3) double precision :: rn(3), f1A, grad1_f1A(3), f2A, grad2_f2A(3), g12, grad1_g12(3) double precision :: tmp1, tmp2 @@ -181,7 +181,7 @@ subroutine grad1_j12_r1_seq(r1, n_grid2, gradx, grady, gradz) ! d/dy1 j(mu,r12) = 0.5 * [(1 - erf(mu * r12)) / r12] * (y1 - y2) ! d/dz1 j(mu,r12) = 0.5 * [(1 - erf(mu * r12)) / r12] * (z1 - z2) - do jpoint = 1, n_points_extra_final_grid ! 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) @@ -191,19 +191,15 @@ subroutine grad1_j12_r1_seq(r1, n_grid2, gradx, grady, gradz) dy = r1(2) - r2(2) dz = r1(3) - r2(3) - r12 = dx * dx + dy * dy + dz * dz - - if(r12 .lt. 1d-20) then - gradx(jpoint) = 0.d0 - grady(jpoint) = 0.d0 - gradz(jpoint) = 0.d0 + 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 - r12_inv = 1.d0/dsqrt(r12) - r12 = r12*r12_inv - - tmp = 0.5d0 * (1.d0 - derf(mu_erf * r12)) * r12_inv + tmp = 0.5d0 * (1.d0 - derf(mu_erf * r12)) / r12 gradx(jpoint) = tmp * dx grady(jpoint) = tmp * dy @@ -212,10 +208,10 @@ subroutine grad1_j12_r1_seq(r1, n_grid2, gradx, grady, gradz) elseif(j2e_type .eq. "Mur") 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) + ! 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) - do jpoint = 1, n_points_extra_final_grid ! 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) @@ -224,29 +220,23 @@ subroutine grad1_j12_r1_seq(r1, n_grid2, gradx, grady, gradz) dx = r1(1) - r2(1) dy = r1(2) - r2(2) dz = r1(3) - r2(3) + r12 = dsqrt(dx * dx + dy * dy + dz * dz) - r12 = 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-20) then + if(r12 .lt. 1d-10) then gradx(jpoint) = 0.d0 grady(jpoint) = 0.d0 gradz(jpoint) = 0.d0 cycle endif - r12_inv = 1.d0/dsqrt(r12) - r12 = r12*r12_inv - - 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) - - tmp = 0.5d0 * (1.d0 - derf(mu_tmp)) * r12_inv + tmp = 0.5d0 * (1.d0 - derf(mu_tmp)) / r12 gradx(jpoint) = gradx(jpoint) + tmp * dx grady(jpoint) = grady(jpoint) + tmp * dy @@ -264,7 +254,7 @@ subroutine grad1_j12_r1_seq(r1, n_grid2, gradx, grady, gradz) PROVIDE a_boys - do jpoint = 1, n_points_extra_final_grid ! 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) @@ -273,17 +263,14 @@ subroutine grad1_j12_r1_seq(r1, n_grid2, gradx, grady, gradz) dx = r1(1) - r2(1) dy = r1(2) - r2(2) dz = r1(3) - r2(3) - r12 = dx * dx + dy * dy + dz * dz - + 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 + gradx(jpoint) = 0.d0 + grady(jpoint) = 0.d0 + gradz(jpoint) = 0.d0 cycle endif - r12 = dsqrt(r12) - tmp = 1.d0 + a_boys * r12 tmp = 0.5d0 / (r12 * tmp * tmp) @@ -294,13 +281,16 @@ subroutine grad1_j12_r1_seq(r1, n_grid2, gradx, grady, gradz) elseif(j2e_type .eq. "Boys_Handy") then - integer :: powmax - powmax = max(maxval(jBH_m),maxval(jBH_n)) - + integer :: powmax1, powmax, powmax2 double precision, allocatable :: f1A_power(:), f2A_power(:), double_p(:), g12_power(:) - allocate (f1A_power(-1:powmax), f2A_power(-1:powmax), g12_power(-1:powmax), double_p(0:powmax)) - do p=0,powmax + powmax1 = max(maxval(jBH_m), maxval(jBH_n)) + powmax2 = maxval(jBH_o) + powmax = max(powmax1, powmax2) + + allocate(f1A_power(-1:powmax), f2A_power(-1:powmax), g12_power(-1:powmax), double_p(0:powmax)) + + do p = 0, powmax double_p(p) = dble(p) enddo @@ -318,11 +308,10 @@ subroutine grad1_j12_r1_seq(r1, n_grid2, gradx, grady, gradz) r2(2) = final_grid_points_extra(2,jpoint) r2(3) = final_grid_points_extra(3,jpoint) - gradx(jpoint) = 0.d0 - grady(jpoint) = 0.d0 - gradz(jpoint) = 0.d0 - - do i_nucl = 1, nucl_num + gradx(jpoint) = 0.d0 + grady(jpoint) = 0.d0 + gradz(jpoint) = 0.d0 + do i_nucl = 1, nucl_num rn(1) = nucl_coord(i_nucl,1) rn(2) = nucl_coord(i_nucl,2) @@ -332,23 +321,15 @@ subroutine grad1_j12_r1_seq(r1, n_grid2, gradx, grady, gradz) call jBH_elem_fct_grad(jBH_en(i_nucl), r2, rn, f2A, grad2_f2A) call jBH_elem_fct_grad(jBH_ee(i_nucl), r1, r2, g12, grad1_g12) - ! Compute powers of f1A and f2A - - do p = 1, maxval(jBH_m(:,i_nucl)) + do p = 1, powmax1 f1A_power(p) = f1A_power(p-1) * f1A - enddo - - do p = 1, maxval(jBH_n(:,i_nucl)) f2A_power(p) = f2A_power(p-1) * f2A enddo - - do p = 1, maxval(jBH_o(:,i_nucl)) + do p = 1, powmax2 g12_power(p) = g12_power(p-1) * g12 enddo - - do p = 1, jBH_size mpA = jBH_m(p,i_nucl) npA = jBH_n(p,i_nucl) @@ -358,27 +339,22 @@ subroutine grad1_j12_r1_seq(r1, n_grid2, gradx, grady, gradz) tmp = tmp * 0.5d0 endif -!TODO : Powers to optimize here - -! tmp1 = 0.d0 -! if(mpA .gt. 0) then -! tmp1 = tmp1 + dble(mpA) * f1A**(mpA-1) * f2A**npA -! endif -! if(npA .gt. 0) then -! tmp1 = tmp1 + dble(npA) * f1A**(npA-1) * f2A**mpA -! endif -! tmp1 = tmp1 * g12**(opA) -! -! tmp2 = 0.d0 -! if(opA .gt. 0) then -! tmp2 = tmp2 + dble(opA) * g12**(opA-1) * (f1A**(mpA) * f2A**(npA) + f1A**(npA) * f2A**(mpA)) -! endif - tmp1 = double_p(mpA) * f1A_power(mpA-1) * f2A_power(npA) + double_p(npA) * f1A_power(npA-1) * f2A_power(mpA) tmp1 = tmp1 * g12_power(opA) - tmp2 = double_p(opA) * g12_power(opA-1) * (f1A_power(mpA) * f2A_power(npA) + f1A_power(npA) * f2A_power(mpA)) + !tmp1 = 0.d0 + !if(mpA .gt. 0) then + ! tmp1 = tmp1 + dble(mpA) * f1A**dble(mpA-1) * f2A**dble(npA) + !endif + !if(npA .gt. 0) then + ! tmp1 = tmp1 + dble(npA) * f1A**dble(npA-1) * f2A**dble(mpA) + !endif + !tmp1 = tmp1 * g12**dble(opA) + !tmp2 = 0.d0 + !if(opA .gt. 0) then + ! tmp2 = tmp2 + dble(opA) * g12**dble(opA-1) * (f1A**dble(mpA) * f2A**dble(npA) + f1A**dble(npA) * f2A**dble(mpA)) + !endif gradx(jpoint) = gradx(jpoint) + tmp * (tmp1 * grad1_f1A(1) + tmp2 * grad1_g12(1)) grady(jpoint) = grady(jpoint) + tmp * (tmp1 * grad1_f1A(2) + tmp2 * grad1_g12(2)) @@ -418,10 +394,10 @@ subroutine grad1_jmu_r1_seq(mu, r1, n_grid2, gradx, grady, gradz) integer :: jpoint double precision :: r2(3) - double precision :: dx, dy, dz, r12, r12_inv, tmp + double precision :: dx, dy, dz, r12, tmp - do jpoint = 1, n_points_extra_final_grid ! 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) @@ -431,19 +407,15 @@ subroutine grad1_jmu_r1_seq(mu, r1, n_grid2, gradx, grady, gradz) dy = r1(2) - r2(2) dz = r1(3) - r2(3) - r12 = dx * dx + dy * dy + dz * dz - - if(r12 .lt. 1d-20) then - gradx(jpoint) = 0.d0 - grady(jpoint) = 0.d0 - gradz(jpoint) = 0.d0 + 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 - r12_inv = 1.d0 / dsqrt(r12) - r12 = r12 * r12_inv - - tmp = 0.5d0 * (1.d0 - derf(mu * r12)) * r12_inv + tmp = 0.5d0 * (1.d0 - derf(mu * r12)) / r12 gradx(jpoint) = tmp * dx grady(jpoint) = tmp * dy @@ -467,7 +439,7 @@ subroutine j12_r1_seq(r1, n_grid2, res) integer :: jpoint double precision :: r2(3) double precision :: dx, dy, dz - double precision :: mu_tmp, r12, mu_erf_inv + double precision :: mu_tmp, r12 PROVIDE final_grid_points_extra @@ -475,21 +447,20 @@ subroutine j12_r1_seq(r1, n_grid2, res) PROVIDE mu_erf - mu_erf_inv = 1.d0 / mu_erf - do jpoint = 1, n_points_extra_final_grid ! 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) - + dx = r1(1) - r2(1) dy = r1(2) - r2(2) dz = r1(3) - r2(3) r12 = dsqrt(dx * dx + dy * dy + dz * dz) mu_tmp = mu_erf * r12 - - res(jpoint) = 0.5d0 * r12 * (1.d0 - derf(mu_tmp)) - inv_sq_pi_2 * dexp(-mu_tmp*mu_tmp) * mu_erf_inv + + res(jpoint) = 0.5d0 * r12 * (1.d0 - derf(mu_tmp)) - inv_sq_pi_2 * dexp(-mu_tmp*mu_tmp) / mu_erf enddo elseif(j2e_type .eq. "Boys") then @@ -498,7 +469,7 @@ subroutine j12_r1_seq(r1, n_grid2, res) PROVIDE a_boys - do jpoint = 1, n_points_extra_final_grid ! 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) @@ -540,19 +511,19 @@ subroutine jmu_r1_seq(mu, r1, n_grid2, res) tmp1 = inv_sq_pi_2 / mu - do jpoint = 1, n_points_extra_final_grid ! 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) - + dx = r1(1) - r2(1) dy = r1(2) - r2(2) dz = r1(3) - r2(3) r12 = dsqrt(dx * dx + dy * dy + dz * dz) tmp2 = mu * r12 - + res(jpoint) = 0.5d0 * r12 * (1.d0 - derf(tmp2)) - tmp1 * dexp(-tmp2*tmp2) enddo @@ -579,7 +550,7 @@ subroutine env_nucl_r1_seq(n_grid2, res) res = 1.d0 - do jpoint = 1, n_points_extra_final_grid ! r2 + do jpoint = 1, n_points_extra_final_grid ! r2 r(1) = final_grid_points_extra(1,jpoint) r(2) = final_grid_points_extra(2,jpoint) r(3) = final_grid_points_extra(3,jpoint) @@ -598,7 +569,7 @@ subroutine env_nucl_r1_seq(n_grid2, res) res = 1.d0 - do jpoint = 1, n_points_extra_final_grid ! r2 + do jpoint = 1, n_points_extra_final_grid ! r2 r(1) = final_grid_points_extra(1,jpoint) r(2) = final_grid_points_extra(2,jpoint) r(3) = final_grid_points_extra(3,jpoint) @@ -618,7 +589,7 @@ subroutine env_nucl_r1_seq(n_grid2, res) res = 1.d0 - do jpoint = 1, n_points_extra_final_grid ! r2 + do jpoint = 1, n_points_extra_final_grid ! r2 r(1) = final_grid_points_extra(1,jpoint) r(2) = final_grid_points_extra(2,jpoint) r(3) = final_grid_points_extra(3,jpoint) @@ -636,7 +607,7 @@ subroutine env_nucl_r1_seq(n_grid2, res) res = 1.d0 - do jpoint = 1, n_points_extra_final_grid ! r2 + do jpoint = 1, n_points_extra_final_grid ! r2 r(1) = final_grid_points_extra(1,jpoint) r(2) = final_grid_points_extra(2,jpoint) r(3) = final_grid_points_extra(3,jpoint) @@ -666,7 +637,7 @@ end subroutine get_grad1_u12_2e_r1_seq(ipoint, n_grid2, resx, resy, resz) BEGIN_DOC - ! + ! ! grad_1 u_2e(r1,r2) ! ! we use grid for r1 and extra_grid for r2 @@ -786,7 +757,7 @@ end subroutine get_u12_2e_r1_seq(ipoint, n_grid2, res) BEGIN_DOC - ! + ! ! u_2e(r1,r2) ! ! we use grid for r1 and extra_grid for r2 @@ -909,7 +880,7 @@ subroutine jBH_elem_fct_grad(alpha, r1, r2, fct, grad1_fct) endif return -end +end ! --- diff --git a/plugins/local/slater_tc/tc_hmat.irp.f b/plugins/local/slater_tc/tc_hmat.irp.f index abec410d..cc780364 100644 --- a/plugins/local/slater_tc/tc_hmat.irp.f +++ b/plugins/local/slater_tc/tc_hmat.irp.f @@ -22,6 +22,7 @@ BEGIN_PROVIDER [double precision, htilde_matrix_elmt_bi_ortho, (N_det,N_det)] if(noL_standard) then PROVIDE noL_0e + print*, "noL_0e =", noL_0e PROVIDE noL_1e PROVIDE noL_2e endif diff --git a/plugins/local/tc_bi_ortho/print_tc_energy.irp.f b/plugins/local/tc_bi_ortho/print_tc_energy.irp.f index 1fa0c6d9..979d792b 100644 --- a/plugins/local/tc_bi_ortho/print_tc_energy.irp.f +++ b/plugins/local/tc_bi_ortho/print_tc_energy.irp.f @@ -9,15 +9,6 @@ program print_tc_energy read_wf = .True. touch read_wf - PROVIDE j2e_type - PROVIDE j1e_type - PROVIDE env_type - - print *, ' j2e_type = ', j2e_type - print *, ' j1e_type = ', j1e_type - print *, ' env_type = ', env_type - - my_grid_becke = .True. PROVIDE tc_grid1_a tc_grid1_r my_n_pt_r_grid = tc_grid1_r @@ -38,6 +29,24 @@ program print_tc_energy call write_int(6, my_n_pt_a_extra_grid, 'angular internal grid over') endif + call main() + +end + +! --- + +subroutine main() + + implicit none + + PROVIDE j2e_type + PROVIDE j1e_type + PROVIDE env_type + + print *, ' j2e_type = ', j2e_type + print *, ' j1e_type = ', j1e_type + print *, ' env_type = ', env_type + call write_tc_energy() end diff --git a/plugins/local/tc_scf/tc_scf.irp.f b/plugins/local/tc_scf/tc_scf.irp.f index f099b90e..83da03ec 100644 --- a/plugins/local/tc_scf/tc_scf.irp.f +++ b/plugins/local/tc_scf/tc_scf.irp.f @@ -7,15 +7,6 @@ program tc_scf END_DOC implicit none - integer :: i - logical :: good_angles - - print *, ' TC-SCF with:' - print *, ' j2e_type = ', j2e_type - print *, ' j1e_type = ', j1e_type - print *, ' env_type = ', env_type - - write(json_unit,json_array_open_fmt) 'tc-scf' my_grid_becke = .True. PROVIDE tc_grid1_a tc_grid1_r @@ -37,6 +28,26 @@ program tc_scf call write_int(6, my_n_pt_a_extra_grid, 'angular internal grid over') endif + call main() + +end + +! --- + +subroutine main() + + implicit none + + integer :: i + logical :: good_angles + + print *, ' TC-SCF with:' + print *, ' j2e_type = ', j2e_type + print *, ' j1e_type = ', j1e_type + print *, ' env_type = ', env_type + + write(json_unit,json_array_open_fmt) 'tc-scf' + call rh_tcscf_diis() PROVIDE Fock_matrix_tc_diag_mo_tot @@ -84,7 +95,7 @@ subroutine create_guess() SOFT_TOUCH mo_label endif -end subroutine create_guess +end ! ---