diff --git a/src/ao_many_one_e_ints/grad_lapl_jmu_modif.irp.f b/src/ao_many_one_e_ints/grad_lapl_jmu_modif.irp.f index 9d34e1d7..18d71134 100644 --- a/src/ao_many_one_e_ints/grad_lapl_jmu_modif.irp.f +++ b/src/ao_many_one_e_ints/grad_lapl_jmu_modif.irp.f @@ -195,7 +195,7 @@ END_PROVIDER ! --- -BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b, (ao_num, ao_num, n_points_final_grid)] +BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_fit, (ao_num, ao_num, n_points_final_grid)] BEGIN_DOC ! @@ -212,14 +212,14 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b, (ao_num, ao_num, n_points_ double precision, external :: overlap_gauss_r12_ao_with1s - print*, ' providing v_ij_u_cst_mu_j1b ...' + print*, ' providing v_ij_u_cst_mu_j1b_fit ...' call wall_time(wall0) provide mu_erf final_grid_points j1b_pen PROVIDE ng_fit_jast expo_gauss_j_mu_x coef_gauss_j_mu_x PROVIDE List_all_comb_b2_size List_all_comb_b2_coef List_all_comb_b2_expo List_all_comb_b2_cent - v_ij_u_cst_mu_j1b = 0.d0 + v_ij_u_cst_mu_j1b_fit = 0.d0 !$OMP PARALLEL DEFAULT (NONE) & !$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center, & @@ -228,7 +228,7 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b, (ao_num, ao_num, n_points_ !$OMP final_grid_points, ng_fit_jast, & !$OMP expo_gauss_j_mu_x, coef_gauss_j_mu_x, & !$OMP List_all_comb_b2_coef, List_all_comb_b2_expo, & - !$OMP List_all_comb_b2_cent, v_ij_u_cst_mu_j1b) + !$OMP List_all_comb_b2_cent, v_ij_u_cst_mu_j1b_fit) !$OMP DO do ipoint = 1, n_points_final_grid r(1) = final_grid_points(1,ipoint) @@ -278,7 +278,7 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b, (ao_num, ao_num, n_points_ enddo - v_ij_u_cst_mu_j1b(j,i,ipoint) = tmp + v_ij_u_cst_mu_j1b_fit(j,i,ipoint) = tmp enddo enddo enddo @@ -288,13 +288,13 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b, (ao_num, ao_num, n_points_ do ipoint = 1, n_points_final_grid do i = 2, ao_num do j = 1, i-1 - v_ij_u_cst_mu_j1b(j,i,ipoint) = v_ij_u_cst_mu_j1b(i,j,ipoint) + v_ij_u_cst_mu_j1b_fit(j,i,ipoint) = v_ij_u_cst_mu_j1b_fit(i,j,ipoint) enddo enddo enddo call wall_time(wall1) - print*, ' wall time for v_ij_u_cst_mu_j1b', wall1 - wall0 + print*, ' wall time for v_ij_u_cst_mu_j1b_fit', wall1 - wall0 END_PROVIDER @@ -327,7 +327,6 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_an, (ao_num, ao_num, n_poin call wall_time(wall0) provide mu_erf final_grid_points j1b_pen - PROVIDE ng_fit_jast expo_gauss_j_mu_x coef_gauss_j_mu_x PROVIDE List_all_comb_b2_size List_all_comb_b2_coef List_all_comb_b2_expo List_all_comb_b2_cent ct = inv_sq_pi_2 / mu_erf @@ -340,7 +339,6 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_an, (ao_num, ao_num, n_poin !$OMP int_e2, int_c3, int_e3) & !$OMP SHARED (n_points_final_grid, ao_num, List_all_comb_b2_size, & !$OMP final_grid_points, mu_erf, ct, & - !$OMP expo_gauss_j_mu_x, coef_gauss_j_mu_x, & !$OMP List_all_comb_b2_coef, List_all_comb_b2_expo, & !$OMP List_all_comb_b2_cent, v_ij_u_cst_mu_j1b_an) !$OMP DO diff --git a/src/bi_ort_ints/bi_ort_ints.irp.f b/src/bi_ort_ints/bi_ort_ints.irp.f index 2ff96326..cac46b18 100644 --- a/src/bi_ort_ints/bi_ort_ints.irp.f +++ b/src/bi_ort_ints/bi_ort_ints.irp.f @@ -9,10 +9,9 @@ program bi_ort_ints implicit none my_grid_becke = .True. - !my_n_pt_r_grid = 10 - !my_n_pt_a_grid = 14 - my_n_pt_r_grid = 30 - my_n_pt_a_grid = 50 + PROVIDE tc_grid1_a tc_grid1_r + my_n_pt_r_grid = tc_grid1_r + my_n_pt_a_grid = tc_grid1_a touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid ! call test_3e diff --git a/src/bi_ort_ints/total_twoe_pot.irp.f b/src/bi_ort_ints/total_twoe_pot.irp.f index c1bacbd0..f03e8a34 100644 --- a/src/bi_ort_ints/total_twoe_pot.irp.f +++ b/src/bi_ort_ints/total_twoe_pot.irp.f @@ -261,51 +261,55 @@ END_PROVIDER ! --- - BEGIN_PROVIDER [ double precision, mo_bi_ortho_tc_two_e_jj, (mo_num,mo_num) ] -&BEGIN_PROVIDER [ double precision, mo_bi_ortho_tc_two_e_jj_exchange, (mo_num,mo_num) ] -&BEGIN_PROVIDER [ double precision, mo_bi_ortho_tc_two_e_jj_anti, (mo_num,mo_num) ] - implicit none + BEGIN_PROVIDER [ double precision, mo_bi_ortho_tc_two_e_jj, (mo_num,mo_num)] +&BEGIN_PROVIDER [ double precision, mo_bi_ortho_tc_two_e_jj_exchange, (mo_num,mo_num)] +&BEGIN_PROVIDER [ double precision, mo_bi_ortho_tc_two_e_jj_anti, (mo_num,mo_num)] + BEGIN_DOC ! mo_bi_ortho_tc_two_e_jj(i,j) = J_ij = ! mo_bi_ortho_tc_two_e_jj_exchange(i,j) = K_ij = ! mo_bi_ortho_tc_two_e_jj_anti(i,j) = J_ij - K_ij END_DOC - integer :: i,j - double precision :: get_two_e_integral + implicit none + integer :: i, j - mo_bi_ortho_tc_two_e_jj = 0.d0 + mo_bi_ortho_tc_two_e_jj = 0.d0 mo_bi_ortho_tc_two_e_jj_exchange = 0.d0 - do i=1,mo_num - do j=1,mo_num - mo_bi_ortho_tc_two_e_jj(i,j) = mo_bi_ortho_tc_two_e(j,i,j,i) + do i = 1, mo_num + do j = 1, mo_num + mo_bi_ortho_tc_two_e_jj(i,j) = mo_bi_ortho_tc_two_e(j,i,j,i) mo_bi_ortho_tc_two_e_jj_exchange(i,j) = mo_bi_ortho_tc_two_e(i,j,j,i) - mo_bi_ortho_tc_two_e_jj_anti(i,j) = mo_bi_ortho_tc_two_e_jj(i,j) - mo_bi_ortho_tc_two_e_jj_exchange(i,j) + mo_bi_ortho_tc_two_e_jj_anti(i,j) = mo_bi_ortho_tc_two_e_jj(i,j) - mo_bi_ortho_tc_two_e_jj_exchange(i,j) enddo enddo END_PROVIDER - BEGIN_PROVIDER [double precision, tc_2e_3idx_coulomb_integrals, (mo_num,mo_num, mo_num)] -&BEGIN_PROVIDER [double precision, tc_2e_3idx_exchange_integrals,(mo_num,mo_num, mo_num)] - implicit none - BEGIN_DOC - ! tc_2e_3idx_coulomb_integrals(j,k,i) = - ! - ! tc_2e_3idx_exchange_integrals(j,k,i) = - END_DOC - integer :: i,j,k,l - double precision :: get_two_e_integral - double precision :: integral +! --- - do i = 1, mo_num - do k = 1, mo_num - do j = 1, mo_num - tc_2e_3idx_coulomb_integrals(j, k,i) = mo_bi_ortho_tc_two_e(j ,k ,j ,i ) - tc_2e_3idx_exchange_integrals(j,k,i) = mo_bi_ortho_tc_two_e(k ,j ,j ,i ) - enddo + BEGIN_PROVIDER [double precision, tc_2e_3idx_coulomb_integrals , (mo_num,mo_num,mo_num)] +&BEGIN_PROVIDER [double precision, tc_2e_3idx_exchange_integrals, (mo_num,mo_num,mo_num)] + + BEGIN_DOC + ! tc_2e_3idx_coulomb_integrals (j,k,i) = + ! tc_2e_3idx_exchange_integrals(j,k,i) = + END_DOC + + implicit none + integer :: i, j, k + + do i = 1, mo_num + do k = 1, mo_num + do j = 1, mo_num + tc_2e_3idx_coulomb_integrals(j, k,i) = mo_bi_ortho_tc_two_e(j ,k ,j ,i ) + tc_2e_3idx_exchange_integrals(j,k,i) = mo_bi_ortho_tc_two_e(k ,j ,j ,i ) + enddo + enddo enddo - enddo END_PROVIDER + +! --- + diff --git a/src/cipsi_tc_bi_ortho/stochastic_cipsi.irp.f b/src/cipsi_tc_bi_ortho/stochastic_cipsi.irp.f index a06f28e9..66d82964 100644 --- a/src/cipsi_tc_bi_ortho/stochastic_cipsi.irp.f +++ b/src/cipsi_tc_bi_ortho/stochastic_cipsi.irp.f @@ -1,47 +1,54 @@ + +! --- + subroutine run_stochastic_cipsi + + BEGIN_DOC + ! Selected Full Configuration Interaction with Stochastic selection and PT2. + END_DOC + use selection_types implicit none - BEGIN_DOC -! Selected Full Configuration Interaction with Stochastic selection and PT2. - END_DOC - integer :: i,j,k,ndet - double precision, allocatable :: zeros(:) - integer :: to_select - type(pt2_type) :: pt2_data, pt2_data_err - logical, external :: qp_stop - logical :: print_pt2 + integer :: i, j, k, ndet + integer :: to_select + logical :: print_pt2 + logical :: has + type(pt2_type) :: pt2_data, pt2_data_err + double precision :: rss + double precision :: correlation_energy_ratio, E_denom, E_tc, norm + double precision :: hf_energy_ref + double precision :: relative_error + double precision, allocatable :: ept2(:), pt1(:), extrap_energy(:) + double precision, allocatable :: zeros(:) - double precision :: rss - double precision, external :: memory_of_double - double precision :: correlation_energy_ratio,E_denom,E_tc,norm - double precision, allocatable :: ept2(:), pt1(:),extrap_energy(:) + logical, external :: qp_stop + double precision, external :: memory_of_double + + PROVIDE mo_l_coef mo_r_coef PROVIDE H_apply_buffer_allocated distributed_davidson - print*,'Diagonal elements of the Fock matrix ' + print*, ' Diagonal elements of the Fock matrix ' do i = 1, mo_num - write(*,*)i,Fock_matrix_tc_mo_tot(i,i) + write(*,*) i, Fock_matrix_tc_mo_tot(i,i) enddo + N_iter = 1 threshold_generators = 1.d0 SOFT_TOUCH threshold_generators rss = memory_of_double(N_states)*4.d0 - call check_mem(rss,irp_here) + call check_mem(rss, irp_here) - allocate (zeros(N_states)) + allocate(zeros(N_states)) call pt2_alloc(pt2_data, N_states) call pt2_alloc(pt2_data_err, N_states) - double precision :: hf_energy_ref - logical :: has - double precision :: relative_error + relative_error = PT2_relative_error - relative_error=PT2_relative_error - - zeros = 0.d0 - pt2_data % pt2 = -huge(1.e0) - pt2_data % rpt2 = -huge(1.e0) - pt2_data % overlap= 0.d0 + zeros = 0.d0 + pt2_data % pt2 = -huge(1.e0) + pt2_data % rpt2 = -huge(1.e0) + pt2_data % overlap = 0.d0 pt2_data % variance = huge(1.e0) !!!! WARNING !!!! SEEMS TO BE PROBLEM WTH make_s2_eigenfunction !!!! THE DETERMINANTS CAN APPEAR TWICE IN THE WFT DURING SELECTION @@ -49,7 +56,7 @@ subroutine run_stochastic_cipsi ! call make_s2_eigenfunction ! endif print_pt2 = .False. - call diagonalize_CI_tc_bi_ortho(ndet, E_tc,norm,pt2_data,print_pt2) + call diagonalize_CI_tc_bi_ortho(ndet, E_tc, norm, pt2_data, print_pt2) ! call routine_save_right @@ -74,14 +81,12 @@ subroutine run_stochastic_cipsi ! soft_touch thresh_it_dav print_pt2 = .True. - do while ( & - (N_det < N_det_max) .and. & - (maxval(abs(pt2_data % pt2(1:N_states))) > pt2_max) & - ) - print*,'maxval(abs(pt2_data % pt2(1:N_states)))',maxval(abs(pt2_data % pt2(1:N_states))) - print*,pt2_max - write(*,'(A)') '--------------------------------------------------------------------------------' + do while( (N_det < N_det_max) .and. & + (maxval(abs(pt2_data % pt2(1:N_states))) > pt2_max)) + print*,'maxval(abs(pt2_data % pt2(1:N_states)))',maxval(abs(pt2_data % pt2(1:N_states))) + print*,pt2_max + write(*,'(A)') '--------------------------------------------------------------------------------' to_select = int(sqrt(dble(N_states))*dble(N_det)*selection_factor) to_select = max(N_states_diag, to_select) @@ -94,8 +99,7 @@ subroutine run_stochastic_cipsi call ZMQ_pt2(E_denom, pt2_data, pt2_data_err, relative_error,to_select) ! Stochastic PT2 and selection ! stop - call print_summary(psi_energy_with_nucl_rep, & - pt2_data, pt2_data_err, N_det,N_configuration,N_states,psi_s2) + call print_summary(psi_energy_with_nucl_rep, pt2_data, pt2_data_err, N_det, N_configuration, N_states, psi_s2) call save_energy(psi_energy_with_nucl_rep, pt2_data % pt2) @@ -109,13 +113,13 @@ subroutine run_stochastic_cipsi ! Add selected determinants call copy_H_apply_buffer_to_wf_tc() - PROVIDE psi_l_coef_bi_ortho psi_r_coef_bi_ortho - PROVIDE psi_det - PROVIDE psi_det_sorted_tc + PROVIDE psi_l_coef_bi_ortho psi_r_coef_bi_ortho + PROVIDE psi_det + PROVIDE psi_det_sorted_tc ept2(N_iter-1) = E_tc + nuclear_repulsion + (pt2_data % pt2(1))/norm - pt1(N_iter-1) = dsqrt(pt2_data % overlap(1,1)) - call diagonalize_CI_tc_bi_ortho(ndet, E_tc,norm,pt2_data,print_pt2) + pt1(N_iter-1) = dsqrt(pt2_data % overlap(1,1)) + call diagonalize_CI_tc_bi_ortho(ndet, E_tc, norm, pt2_data, print_pt2) ! stop if (qp_stop()) exit enddo diff --git a/src/fci_tc_bi/diagonalize_ci.irp.f b/src/fci_tc_bi/diagonalize_ci.irp.f index 8b59785a..6c8f3431 100644 --- a/src/fci_tc_bi/diagonalize_ci.irp.f +++ b/src/fci_tc_bi/diagonalize_ci.irp.f @@ -17,6 +17,8 @@ subroutine diagonalize_CI_tc_bi_ortho(ndet, E_tc, norm, pt2_data, print_pt2) integer :: i, j double precision :: pt2_tmp, pt1_norm, rpt2_tmp, abs_pt2 + PROVIDE mo_l_coef mo_r_coef + pt2_tmp = pt2_data % pt2(1) abs_pt2 = pt2_data % variance(1) pt1_norm = pt2_data % overlap(1,1) @@ -60,7 +62,7 @@ subroutine diagonalize_CI_tc_bi_ortho(ndet, E_tc, norm, pt2_data, print_pt2) psi_coef(i,j) = dabs(psi_l_coef_bi_ortho(i,j) * psi_r_coef_bi_ortho(i,j)) enddo enddo - SOFT_TOUCH eigval_left_tc_bi_orth eigval_right_tc_bi_orth leigvec_tc_bi_orth reigvec_tc_bi_orth norm_ground_left_right_bi_orth + SOFT_TOUCH eigval_left_tc_bi_orth eigval_right_tc_bi_orth leigvec_tc_bi_orth reigvec_tc_bi_orth norm_ground_left_right_bi_orth SOFT_TOUCH psi_l_coef_bi_ortho psi_r_coef_bi_ortho psi_coef psi_energy psi_s2 call save_tc_bi_ortho_wavefunction() diff --git a/src/fci_tc_bi/fci_tc_bi_ortho.irp.f b/src/fci_tc_bi/fci_tc_bi_ortho.irp.f index cf3361bf..1c1c0411 100644 --- a/src/fci_tc_bi/fci_tc_bi_ortho.irp.f +++ b/src/fci_tc_bi/fci_tc_bi_ortho.irp.f @@ -1,5 +1,8 @@ -program fci - implicit none + +! --- + +program fci_tc_bi + BEGIN_DOC ! Selected Full Configuration Interaction with stochastic selection ! and PT2. @@ -36,10 +39,12 @@ program fci ! END_DOC + implicit none my_grid_becke = .True. - my_n_pt_r_grid = 30 - my_n_pt_a_grid = 50 + PROVIDE tc_grid1_a tc_grid1_r + my_n_pt_r_grid = tc_grid1_r + my_n_pt_a_grid = tc_grid1_a touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid pruning = -1.d0 @@ -62,18 +67,18 @@ subroutine run_cipsi_tc() PROVIDE psi_det psi_coef mo_bi_ortho_tc_two_e mo_bi_ortho_tc_one_e - if(elec_alpha_num+elec_beta_num .ge. 3) then + if((elec_alpha_num+elec_beta_num) .ge. 3) then if(three_body_h_tc) then call provide_all_three_ints_bi_ortho() endif endif - FREE int2_grad1_u12_ao - FREE int2_grad1_u12_bimo_transp int2_grad1_u12_ao_transp + FREE int2_grad1_u12_ao int2_grad1_u12_ao_t int2_grad1_u12_ao_transp + FREE int2_grad1_u12_bimo_transp write(json_unit,json_array_open_fmt) 'fci_tc' - if (do_pt2) then + if(do_pt2) then call run_stochastic_cipsi() else call run_cipsi() @@ -88,14 +93,14 @@ subroutine run_cipsi_tc() PROVIDE mo_bi_ortho_tc_one_e mo_bi_ortho_tc_two_e pt2_min_parallel_tasks - if(elec_alpha_num+elec_beta_num .ge. 3) then + if((elec_alpha_num+elec_beta_num) .ge. 3) then if(three_body_h_tc) then call provide_all_three_ints_bi_ortho() endif endif - FREE int2_grad1_u12_ao - FREE int2_grad1_u12_bimo_transp int2_grad1_u12_ao_transp + FREE int2_grad1_u12_ao int2_grad1_u12_ao_t int2_grad1_u12_ao_transp + FREE int2_grad1_u12_bimo_transp call run_slave_cipsi diff --git a/src/fci_tc_bi/pt2_tc.irp.f b/src/fci_tc_bi/pt2_tc.irp.f index 96a54825..390042bf 100644 --- a/src/fci_tc_bi/pt2_tc.irp.f +++ b/src/fci_tc_bi/pt2_tc.irp.f @@ -1,31 +1,42 @@ + +! --- + program tc_pt2_prog + implicit none + my_grid_becke = .True. - my_n_pt_r_grid = 30 - my_n_pt_a_grid = 50 + PROVIDE tc_grid1_a tc_grid1_r + my_n_pt_r_grid = tc_grid1_r + my_n_pt_a_grid = tc_grid1_a touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid + pruning = -1.d0 touch pruning + ! pt2_relative_error = 0.01d0 ! touch pt2_relative_error - call run_pt2_tc + call run_pt2_tc() end +! --- -subroutine run_pt2_tc +subroutine run_pt2_tc() - implicit none + implicit none PROVIDE psi_det psi_coef mo_bi_ortho_tc_two_e mo_bi_ortho_tc_one_e - if(elec_alpha_num+elec_beta_num.ge.3)then + + if(elec_alpha_num+elec_beta_num.ge.3) then if(three_body_h_tc)then - call provide_all_three_ints_bi_ortho + call provide_all_three_ints_bi_ortho() endif endif - ! --- - - call tc_pt2 + call tc_pt2() end + +! --- + diff --git a/src/non_h_ints_mu/debug_fit.irp.f b/src/non_h_ints_mu/debug_fit.irp.f index 146028d5..05d2db68 100644 --- a/src/non_h_ints_mu/debug_fit.irp.f +++ b/src/non_h_ints_mu/debug_fit.irp.f @@ -6,13 +6,9 @@ program debug_fit implicit none my_grid_becke = .True. - - my_n_pt_r_grid = 30 - my_n_pt_a_grid = 50 - !my_n_pt_r_grid = 100 - !my_n_pt_a_grid = 170 - !my_n_pt_r_grid = 150 - !my_n_pt_a_grid = 194 + PROVIDE tc_grid1_a tc_grid1_r + my_n_pt_r_grid = tc_grid1_r + my_n_pt_a_grid = tc_grid1_a touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid PROVIDE mu_erf j1b_pen diff --git a/src/non_h_ints_mu/debug_integ_jmu_modif.irp.f b/src/non_h_ints_mu/debug_integ_jmu_modif.irp.f index f3e93360..b9e8df25 100644 --- a/src/non_h_ints_mu/debug_integ_jmu_modif.irp.f +++ b/src/non_h_ints_mu/debug_integ_jmu_modif.irp.f @@ -6,13 +6,9 @@ program debug_integ_jmu_modif implicit none my_grid_becke = .True. - - !my_n_pt_r_grid = 30 - !my_n_pt_a_grid = 50 - !my_n_pt_r_grid = 100 - !my_n_pt_a_grid = 170 - my_n_pt_r_grid = 150 - my_n_pt_a_grid = 194 + PROVIDE tc_grid1_a tc_grid1_r + my_n_pt_r_grid = tc_grid1_r + my_n_pt_a_grid = tc_grid1_a touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid PROVIDE mu_erf j1b_pen @@ -48,22 +44,21 @@ subroutine test_v_ij_u_cst_mu_j1b() print*, ' test_v_ij_u_cst_mu_j1b ...' - PROVIDE v_ij_u_cst_mu_j1b + PROVIDE v_ij_u_cst_mu_j1b_fit eps_ij = 1d-3 acc_tot = 0.d0 normalz = 0.d0 - !do ipoint = 1, 10 do ipoint = 1, n_points_final_grid do j = 1, ao_num do i = 1, ao_num - i_exc = v_ij_u_cst_mu_j1b(i,j,ipoint) - i_num = num_v_ij_u_cst_mu_j1b(i,j,ipoint) + i_exc = v_ij_u_cst_mu_j1b_fit(i,j,ipoint) + i_num = num_v_ij_u_cst_mu_j1b (i,j,ipoint) acc_ij = dabs(i_exc - i_num) if(acc_ij .gt. eps_ij) then - print *, ' problem in v_ij_u_cst_mu_j1b on', i, j, ipoint + print *, ' problem in v_ij_u_cst_mu_j1b_fit on', i, j, ipoint print *, ' analyt integ = ', i_exc print *, ' numeri integ = ', i_num print *, ' diff = ', acc_ij diff --git a/src/non_h_ints_mu/new_grad_tc.irp.f b/src/non_h_ints_mu/new_grad_tc.irp.f index 499ffe9d..dc76431d 100644 --- a/src/non_h_ints_mu/new_grad_tc.irp.f +++ b/src/non_h_ints_mu/new_grad_tc.irp.f @@ -1,68 +1,3 @@ -! --- - -!BEGIN_PROVIDER [ double precision, int1_grad2_u12_ao, (3, ao_num, ao_num, n_points_final_grid)] -! -! BEGIN_DOC -! ! -! ! int1_grad2_u12_ao(:,i,j,ipoint) = \int dr1 [-1 * \grad_r2 J(r1,r2)] \phi_i(r1) \phi_j(r1) -! ! -! ! where r1 = r(ipoint) -! ! -! ! if J(r1,r2) = u12: -! ! -! ! int1_grad2_u12_ao(:,i,j,ipoint) = +0.5 x \int dr1 [-(r1 - r2) (erf(mu * r12)-1)r_12] \phi_i(r1) \phi_j(r1) -! ! = -0.5 * [ v_ij_erf_rk_cst_mu(i,j,ipoint) * r(:) - x_v_ij_erf_rk_cst_mu(i,j,ipoint,:) ] -! ! = -int2_grad1_u12_ao(i,j,ipoint,:) -! ! -! ! if J(r1,r2) = u12 x v1 x v2 -! ! -! ! int1_grad2_u12_ao(:,i,j,ipoint) = v2 x [ 0.5 x \int dr1 [-(r1 - r2) (erf(mu * r12)-1)r_12] v1 \phi_i(r1) \phi_j(r1) ] -! ! - \grad_2 v2 x [ \int dr1 u12 v1 \phi_i(r1) \phi_j(r1) ] -! ! = -0.5 v_1b(ipoint) * v_ij_erf_rk_cst_mu_j1b(i,j,ipoint) * r(:) -! ! + 0.5 v_1b(ipoint) * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,:) -! ! - v_1b_grad[:,ipoint] * v_ij_u_cst_mu_j1b(i,j,ipoint) -! ! -! ! -! END_DOC -! -! implicit none -! integer :: ipoint, i, j -! double precision :: x, y, z, tmp_x, tmp_y, tmp_z, tmp0, tmp1, tmp2 -! -! PROVIDE j1b_type -! -! if(j1b_type .eq. 3) then -! -! do ipoint = 1, n_points_final_grid -! x = final_grid_points(1,ipoint) -! y = final_grid_points(2,ipoint) -! z = final_grid_points(3,ipoint) -! -! tmp0 = 0.5d0 * v_1b(ipoint) -! tmp_x = v_1b_grad(1,ipoint) -! tmp_y = v_1b_grad(2,ipoint) -! tmp_z = v_1b_grad(3,ipoint) -! -! do j = 1, ao_num -! do i = 1, ao_num -! -! tmp1 = tmp0 * v_ij_erf_rk_cst_mu_j1b(i,j,ipoint) -! tmp2 = v_ij_u_cst_mu_j1b(i,j,ipoint) -! -! int1_grad2_u12_ao(1,i,j,ipoint) = -tmp1 * x + tmp0 * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,1) - tmp2 * tmp_x -! int1_grad2_u12_ao(2,i,j,ipoint) = -tmp1 * y + tmp0 * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,2) - tmp2 * tmp_y -! int1_grad2_u12_ao(3,i,j,ipoint) = -tmp1 * z + tmp0 * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,3) - tmp2 * tmp_z -! enddo -! enddo -! enddo -! -! else -! -! int1_grad2_u12_ao = -1.d0 * int2_grad1_u12_ao -! -! endif -! -!END_PROVIDER ! --- @@ -98,22 +33,14 @@ BEGIN_PROVIDER [double precision, tc_grad_and_lapl_ao_loop, (ao_num, ao_num, ao_ weight1 = 0.5d0 * final_weight_at_r_vector(ipoint) do i = 1, ao_num - !ao_i_r = weight1 * aos_in_r_array_transp (ipoint,i) - !ao_i_dx = weight1 * aos_grad_in_r_array_transp_bis(ipoint,i,1) - !ao_i_dy = weight1 * aos_grad_in_r_array_transp_bis(ipoint,i,2) - !ao_i_dz = weight1 * aos_grad_in_r_array_transp_bis(ipoint,i,3) ao_i_r = weight1 * aos_in_r_array (i,ipoint) ao_i_dx = weight1 * aos_grad_in_r_array(i,ipoint,1) ao_i_dy = weight1 * aos_grad_in_r_array(i,ipoint,2) ao_i_dz = weight1 * aos_grad_in_r_array(i,ipoint,3) do k = 1, ao_num - !ao_k_r = aos_in_r_array_transp(ipoint,k) ao_k_r = aos_in_r_array(k,ipoint) - !tmp_x = ao_k_r * ao_i_dx - ao_i_r * aos_grad_in_r_array_transp_bis(ipoint,k,1) - !tmp_y = ao_k_r * ao_i_dy - ao_i_r * aos_grad_in_r_array_transp_bis(ipoint,k,2) - !tmp_z = ao_k_r * ao_i_dz - ao_i_r * aos_grad_in_r_array_transp_bis(ipoint,k,3) tmp_x = ao_k_r * ao_i_dx - ao_i_r * aos_grad_in_r_array(k,ipoint,1) tmp_y = ao_k_r * ao_i_dy - ao_i_r * aos_grad_in_r_array(k,ipoint,2) tmp_z = ao_k_r * ao_i_dz - ao_i_r * aos_grad_in_r_array(k,ipoint,3) @@ -134,44 +61,11 @@ BEGIN_PROVIDER [double precision, tc_grad_and_lapl_ao_loop, (ao_num, ao_num, ao_ ! --- - !do ipoint = 1, n_points_final_grid - ! weight1 = 0.5d0 * final_weight_at_r_vector(ipoint) - - ! do l = 1, ao_num - ! ao_l_r = weight1 * aos_in_r_array_transp (ipoint,l) - ! ao_l_dx = weight1 * aos_grad_in_r_array_transp_bis(ipoint,l,1) - ! ao_l_dy = weight1 * aos_grad_in_r_array_transp_bis(ipoint,l,2) - ! ao_l_dz = weight1 * aos_grad_in_r_array_transp_bis(ipoint,l,3) - - ! do j = 1, ao_num - ! ao_j_r = aos_in_r_array_transp(ipoint,j) - - ! tmp_x = ao_j_r * ao_l_dx - ao_l_r * aos_grad_in_r_array_transp_bis(ipoint,j,1) - ! tmp_y = ao_j_r * ao_l_dy - ao_l_r * aos_grad_in_r_array_transp_bis(ipoint,j,2) - ! tmp_z = ao_j_r * ao_l_dz - ao_l_r * aos_grad_in_r_array_transp_bis(ipoint,j,3) - - ! do i = 1, ao_num - ! do k = 1, ao_num - - ! contrib_x = int2_grad1_u12_ao(k,i,ipoint,1) * tmp_x - ! contrib_y = int2_grad1_u12_ao(k,i,ipoint,2) * tmp_y - ! contrib_z = int2_grad1_u12_ao(k,i,ipoint,3) * tmp_z - - ! ac_mat(k,i,l,j) += contrib_x + contrib_y + contrib_z - ! enddo - ! enddo - ! enddo - ! enddo - !enddo - - ! --- - do j = 1, ao_num do l = 1, ao_num do i = 1, ao_num do k = 1, ao_num tc_grad_and_lapl_ao_loop(k,i,l,j) = ac_mat(k,i,l,j) + ac_mat(l,j,k,i) - !tc_grad_and_lapl_ao_loop(k,i,l,j) = ac_mat(k,i,l,j) enddo enddo enddo diff --git a/src/non_h_ints_mu/test_non_h_ints.irp.f b/src/non_h_ints_mu/test_non_h_ints.irp.f index a6e0a311..aff53c2d 100644 --- a/src/non_h_ints_mu/test_non_h_ints.irp.f +++ b/src/non_h_ints_mu/test_non_h_ints.irp.f @@ -1,19 +1,18 @@ + +! --- + program test_non_h implicit none my_grid_becke = .True. - my_n_pt_r_grid = 50 - my_n_pt_a_grid = 74 - !my_n_pt_r_grid = 400 - !my_n_pt_a_grid = 974 - -! my_n_pt_r_grid = 10 ! small grid for quick debug -! my_n_pt_a_grid = 26 ! small grid for quick debug + PROVIDE tc_grid1_a tc_grid1_r + my_n_pt_r_grid = tc_grid1_r + my_n_pt_a_grid = tc_grid1_a touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid - !call routine_grad_squared - !call routine_fit + !call routine_grad_squared() + !call routine_fit() call test_ipp() end diff --git a/src/tc_bi_ortho/compute_deltamu_right.irp.f b/src/tc_bi_ortho/compute_deltamu_right.irp.f index 7ca2c890..ab9dc093 100644 --- a/src/tc_bi_ortho/compute_deltamu_right.irp.f +++ b/src/tc_bi_ortho/compute_deltamu_right.irp.f @@ -3,8 +3,9 @@ program compute_deltamu_right implicit none my_grid_becke = .True. - my_n_pt_r_grid = 30 - my_n_pt_a_grid = 50 + PROVIDE tc_grid1_a tc_grid1_r + my_n_pt_r_grid = tc_grid1_r + my_n_pt_a_grid = tc_grid1_a touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid read_wf = .True. diff --git a/src/tc_bi_ortho/print_tc_dump.irp.f b/src/tc_bi_ortho/print_tc_dump.irp.f index 0a7e08d2..868de444 100644 --- a/src/tc_bi_ortho/print_tc_dump.irp.f +++ b/src/tc_bi_ortho/print_tc_dump.irp.f @@ -6,10 +6,9 @@ program tc_bi_ortho implicit none my_grid_becke = .True. - my_n_pt_r_grid = 30 - my_n_pt_a_grid = 50 - !my_n_pt_r_grid = 100 - !my_n_pt_a_grid = 170 + PROVIDE tc_grid1_a tc_grid1_r + my_n_pt_r_grid = tc_grid1_r + my_n_pt_a_grid = tc_grid1_a touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid call ERI_dump() diff --git a/src/tc_bi_ortho/print_tc_energy.irp.f b/src/tc_bi_ortho/print_tc_energy.irp.f index 522f4cd7..7bca72a1 100644 --- a/src/tc_bi_ortho/print_tc_energy.irp.f +++ b/src/tc_bi_ortho/print_tc_energy.irp.f @@ -7,16 +7,12 @@ program print_tc_energy implicit none print *, 'Hello world' + my_grid_becke = .True. - - !my_n_pt_r_grid = 30 - !my_n_pt_a_grid = 50 - - my_n_pt_r_grid = 100 - my_n_pt_a_grid = 170 - - !my_n_pt_r_grid = 100 - !my_n_pt_a_grid = 266 + PROVIDE tc_grid1_a tc_grid1_r + my_n_pt_r_grid = tc_grid1_r + my_n_pt_a_grid = tc_grid1_a + touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid read_wf = .True. touch read_wf @@ -24,8 +20,6 @@ program print_tc_energy PROVIDE j1b_type print*, 'j1b_type = ', j1b_type - touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid - call write_tc_energy() end diff --git a/src/tc_bi_ortho/print_tc_spin_dens.irp.f b/src/tc_bi_ortho/print_tc_spin_dens.irp.f index 8308140d..c7da5bc8 100644 --- a/src/tc_bi_ortho/print_tc_spin_dens.irp.f +++ b/src/tc_bi_ortho/print_tc_spin_dens.irp.f @@ -1,16 +1,26 @@ + +! --- + program test_spin_dens - implicit none + BEGIN_DOC -! TODO : Reads psi_det in the EZFIO folder and prints out the left- and right-eigenvectors together with the energy. Saves the left-right wave functions at the end. + ! TODO : Reads psi_det in the EZFIO folder and prints out the left- and right-eigenvectors together with the energy. Saves the left-right wave functions at the end. END_DOC + + implicit none + print *, 'Hello world' + my_grid_becke = .True. - my_n_pt_r_grid = 30 - my_n_pt_a_grid = 50 + PROVIDE tc_grid1_a tc_grid1_r + my_n_pt_r_grid = tc_grid1_r + my_n_pt_a_grid = tc_grid1_a + touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid + read_wf = .True. touch read_wf - touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid - call tc_print_mulliken_sd -! call test + + call tc_print_mulliken_sd() + !call test end diff --git a/src/tc_bi_ortho/print_tc_var.irp.f b/src/tc_bi_ortho/print_tc_var.irp.f index fa0a4363..bec34f18 100644 --- a/src/tc_bi_ortho/print_tc_var.irp.f +++ b/src/tc_bi_ortho/print_tc_var.irp.f @@ -7,12 +7,15 @@ program print_tc_var implicit none print *, 'Hello world' + my_grid_becke = .True. - my_n_pt_r_grid = 30 - my_n_pt_a_grid = 50 + PROVIDE tc_grid1_a tc_grid1_r + my_n_pt_r_grid = tc_grid1_r + my_n_pt_a_grid = tc_grid1_a + touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid + read_wf = .True. touch read_wf - touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid call write_tc_var() diff --git a/src/tc_bi_ortho/print_tc_wf.irp.f b/src/tc_bi_ortho/print_tc_wf.irp.f index 0c4198a9..c755485b 100644 --- a/src/tc_bi_ortho/print_tc_wf.irp.f +++ b/src/tc_bi_ortho/print_tc_wf.irp.f @@ -1,20 +1,31 @@ + +! --- + program print_tc_bi_ortho - implicit none + BEGIN_DOC -! TODO : Put the documentation of the program here + ! TODO : Put the documentation of the program here END_DOC + + implicit none + print *, 'Hello world' + my_grid_becke = .True. - my_n_pt_r_grid = 30 - my_n_pt_a_grid = 50 + PROVIDE tc_grid1_a tc_grid1_r + my_n_pt_r_grid = tc_grid1_r + my_n_pt_a_grid = tc_grid1_a + touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid + read_wf = .True. touch read_wf - touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid + ! if(three_body_h_tc)then ! call provide_all_three_ints_bi_ortho ! endif ! call routine - call write_l_r_wf + call write_l_r_wf + end subroutine write_l_r_wf diff --git a/src/tc_bi_ortho/pt2_tc_cisd.irp.f b/src/tc_bi_ortho/pt2_tc_cisd.irp.f index 9cb9a600..8940a4f6 100644 --- a/src/tc_bi_ortho/pt2_tc_cisd.irp.f +++ b/src/tc_bi_ortho/pt2_tc_cisd.irp.f @@ -7,12 +7,16 @@ program pt2_tc_cisd ! END_DOC + implicit none + my_grid_becke = .True. - my_n_pt_r_grid = 30 - my_n_pt_a_grid = 50 + PROVIDE tc_grid1_a tc_grid1_r + my_n_pt_r_grid = tc_grid1_r + my_n_pt_a_grid = tc_grid1_a + touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid + read_wf = .True. touch read_wf - touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid print*, ' nb of states = ', N_states print*, ' nb of det = ', N_det diff --git a/src/tc_bi_ortho/save_tc_bi_ortho_nat.irp.f b/src/tc_bi_ortho/save_tc_bi_ortho_nat.irp.f index 8b6eb1d1..47ade8df 100644 --- a/src/tc_bi_ortho/save_tc_bi_ortho_nat.irp.f +++ b/src/tc_bi_ortho/save_tc_bi_ortho_nat.irp.f @@ -1,35 +1,59 @@ - program tc_natorb_bi_ortho - implicit none - BEGIN_DOC - ! TODO : Put the documentation of the program here - END_DOC - print *, 'Hello world' - my_grid_becke = .True. - my_n_pt_r_grid = 30 - my_n_pt_a_grid = 50 - read_wf = .True. - touch read_wf - touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid - call print_energy_and_mos - call save_tc_natorb -! call minimize_tc_orb_angles - end - - subroutine save_tc_natorb + +! --- + +program tc_natorb_bi_ortho + + BEGIN_DOC + ! TODO : Put the documentation of the program here + END_DOC + implicit none + + print *, 'Hello world' + + my_grid_becke = .True. + PROVIDE tc_grid1_a tc_grid1_r + my_n_pt_r_grid = tc_grid1_r + my_n_pt_a_grid = tc_grid1_a + touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid + + read_wf = .True. + touch read_wf + + call print_energy_and_mos() + call save_tc_natorb() + !call minimize_tc_orb_angles() + +end + +! --- + +subroutine save_tc_natorb() + + implicit none + print*,'Saving the natorbs ' + provide natorb_tc_leigvec_ao natorb_tc_reigvec_ao + call ezfio_set_bi_ortho_mos_mo_l_coef(natorb_tc_leigvec_ao) call ezfio_set_bi_ortho_mos_mo_r_coef(natorb_tc_reigvec_ao) - call save_ref_determinant_nstates_1 + call save_ref_determinant_nstates_1() call ezfio_set_determinants_read_wf(.False.) - end + +end + +! --- - subroutine save_ref_determinant_nstates_1 - implicit none +subroutine save_ref_determinant_nstates_1() + use bitmasks - double precision :: buffer(1,N_states) + implicit none + double precision :: buffer(1,N_states) + buffer = 0.d0 buffer(1,1) = 1.d0 - call save_wavefunction_general(1,1,ref_bitmask,1,buffer) - end + call save_wavefunction_general(1, 1, ref_bitmask, 1, buffer) + +end + diff --git a/src/tc_bi_ortho/select_dets_bi_ortho.irp.f b/src/tc_bi_ortho/select_dets_bi_ortho.irp.f index e6bf3d6e..e923548a 100644 --- a/src/tc_bi_ortho/select_dets_bi_ortho.irp.f +++ b/src/tc_bi_ortho/select_dets_bi_ortho.irp.f @@ -1,15 +1,24 @@ -program tc_bi_ortho - implicit none + +! --- + +program select_dets_bi_ortho() + BEGIN_DOC -! TODO : Put the documentation of the program here + ! TODO : Put the documentation of the program here END_DOC + + implicit none + print *, 'Hello world' + my_grid_becke = .True. - my_n_pt_r_grid = 30 - my_n_pt_a_grid = 50 + PROVIDE tc_grid1_a tc_grid1_r + my_n_pt_r_grid = tc_grid1_r + my_n_pt_a_grid = tc_grid1_a + touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid + read_wf = .True. touch read_wf - touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid !!!!!!!!!!!!!!! WARNING NO 3-BODY !!!!!!!!!!!!!!! WARNING NO 3-BODY @@ -22,6 +31,8 @@ program tc_bi_ortho ! call test end +! --- + subroutine routine_test implicit none use bitmasks ! you need to include the bitmasks_module.f90 features @@ -57,5 +68,7 @@ subroutine routine_test enddo call save_wavefunction_general(n_good,n_states,dets,n_good,coef_new) - end + +! --- + diff --git a/src/tc_bi_ortho/slater_tc_3e_slow.irp.f b/src/tc_bi_ortho/slater_tc_3e_slow.irp.f index e6e78534..3a2e3606 100644 --- a/src/tc_bi_ortho/slater_tc_3e_slow.irp.f +++ b/src/tc_bi_ortho/slater_tc_3e_slow.irp.f @@ -1,4 +1,6 @@ +! --- + subroutine diag_htilde_three_body_ints_bi_ort_slow(Nint, key_i, hthree) BEGIN_DOC @@ -22,12 +24,12 @@ subroutine diag_htilde_three_body_ints_bi_ort_slow(Nint, key_i, hthree) if(core_tc_op) then do i = 1, Nint - key_i_core(i,1) = xor(key_i(i,1),core_bitmask(i,1)) - key_i_core(i,2) = xor(key_i(i,2),core_bitmask(i,2)) + key_i_core(i,1) = xor(key_i(i,1), core_bitmask(i,1)) + key_i_core(i,2) = xor(key_i(i,2), core_bitmask(i,2)) enddo - call bitstring_to_list_ab(key_i_core,occ,Ne,Nint) + call bitstring_to_list_ab(key_i_core, occ, Ne, Nint) else - call bitstring_to_list_ab(key_i,occ,Ne,Nint) + call bitstring_to_list_ab(key_i, occ, Ne, Nint) endif hthree = 0.d0 diff --git a/src/tc_bi_ortho/slater_tc_opt.irp.f b/src/tc_bi_ortho/slater_tc_opt.irp.f index d1c1d45d..002f870e 100644 --- a/src/tc_bi_ortho/slater_tc_opt.irp.f +++ b/src/tc_bi_ortho/slater_tc_opt.irp.f @@ -47,13 +47,19 @@ subroutine htilde_mu_mat_opt_bi_ortho_tot(key_j, key_i, Nint, htot) END_DOC use bitmasks - integer, intent(in) :: Nint - integer(bit_kind), intent(in) :: key_i(Nint,2), key_j(Nint,2) - double precision, intent(out) :: htot - double precision :: hmono, htwoe, hthree - call htilde_mu_mat_opt_bi_ortho(key_j, key_i, Nint, hmono, htwoe, hthree, htot) + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: key_i(Nint,2), key_j(Nint,2) + double precision, intent(out) :: htot + double precision :: hmono, htwoe, hthree + + call htilde_mu_mat_opt_bi_ortho(key_j, key_i, Nint, hmono, htwoe, hthree, htot) + end + +! --- + subroutine htilde_mu_mat_opt_bi_ortho(key_j, key_i, Nint, hmono, htwoe, hthree, htot) + BEGIN_DOC ! ! where |key_j> is developed on the LEFT basis and |key_i> is developed on the RIGHT basis @@ -81,11 +87,11 @@ subroutine htilde_mu_mat_opt_bi_ortho(key_j, key_i, Nint, hmono, htwoe, hthree, call get_excitation_degree(key_i, key_j, degree, Nint) if(degree.gt.2) return - if(degree == 0)then + if(degree == 0) then call diag_htilde_mu_mat_fock_bi_ortho (Nint, key_i, hmono, htwoe, hthree, htot) - else if (degree == 1)then - call single_htilde_mu_mat_fock_bi_ortho(Nint,key_j, key_i , hmono, htwoe, hthree, htot) - else if(degree == 2)then + else if (degree == 1) then + call single_htilde_mu_mat_fock_bi_ortho(Nint, key_j, key_i , hmono, htwoe, hthree, htot) + else if(degree == 2) then call double_htilde_mu_mat_fock_bi_ortho(Nint, key_j, key_i, hmono, htwoe, hthree, htot) endif diff --git a/src/tc_bi_ortho/slater_tc_opt_diag.irp.f b/src/tc_bi_ortho/slater_tc_opt_diag.irp.f index 5c156a4d..d95c87b1 100644 --- a/src/tc_bi_ortho/slater_tc_opt_diag.irp.f +++ b/src/tc_bi_ortho/slater_tc_opt_diag.irp.f @@ -111,7 +111,7 @@ end ! --- -subroutine ac_tc_operator(iorb,ispin,key,hmono,htwoe,hthree,Nint,na,nb) +subroutine ac_tc_operator(iorb, ispin, key, hmono, htwoe, hthree, Nint, na, nb) BEGIN_DOC ! Routine that computes one- and two-body energy corresponding @@ -127,17 +127,17 @@ subroutine ac_tc_operator(iorb,ispin,key,hmono,htwoe,hthree,Nint,na,nb) use bitmasks implicit none - integer, intent(in) :: iorb, ispin, Nint - integer, intent(inout) :: na, nb + integer, intent(in) :: iorb, ispin, Nint + integer, intent(inout) :: na, nb integer(bit_kind), intent(inout) :: key(Nint,2) - double precision, intent(inout) :: hmono,htwoe,hthree + double precision, intent(inout) :: hmono, htwoe, hthree - integer :: occ(Nint*bit_kind_size,2) - integer :: other_spin - integer :: k,l,i,jj,mm,j,m - double precision :: direct_int, exchange_int + integer :: occ(Nint*bit_kind_size,2) + integer :: other_spin + integer :: k, l, i, jj, mm, j, m + integer :: tmp(2) + double precision :: direct_int, exchange_int - if (iorb < 1) then print *, irp_here, ': iorb < 1' print *, iorb, mo_num @@ -153,7 +153,6 @@ subroutine ac_tc_operator(iorb,ispin,key,hmono,htwoe,hthree,Nint,na,nb) ASSERT (ispin < 3) ASSERT (Nint > 0) - integer :: tmp(2) !DIR$ FORCEINLINE call bitstring_to_list_ab(key, occ, tmp, Nint) ASSERT (tmp(1) == elec_alpha_num) @@ -169,50 +168,54 @@ subroutine ac_tc_operator(iorb,ispin,key,hmono,htwoe,hthree,Nint,na,nb) hmono = hmono + mo_bi_ortho_tc_one_e(iorb,iorb) ! Same spin - do i=1,na + do i = 1, na htwoe = htwoe + mo_bi_ortho_tc_two_e_jj_anti(occ(i,ispin),iorb) enddo ! Opposite spin - do i=1,nb + do i = 1, nb htwoe = htwoe + mo_bi_ortho_tc_two_e_jj(occ(i,other_spin),iorb) enddo - if(three_body_h_tc.and.elec_num.gt.2.and.three_e_3_idx_term)then - !!!!! 3-e part - !! same-spin/same-spin - do j = 1, na - jj = occ(j,ispin) - do m = j+1, na - mm = occ(m,ispin) - hthree += three_e_diag_parrallel_spin_prov(mm,jj,iorb) + if(three_body_h_tc .and. (elec_num.gt.2) .and. three_e_3_idx_term) then + + !!!!! 3-e part + !! same-spin/same-spin + do j = 1, na + jj = occ(j,ispin) + do m = j+1, na + mm = occ(m,ispin) + hthree += three_e_diag_parrallel_spin_prov(mm,jj,iorb) + enddo enddo - enddo - !! same-spin/oposite-spin - do j = 1, na - jj = occ(j,ispin) - do m = 1, nb - mm = occ(m,other_spin) - direct_int = three_e_3_idx_direct_bi_ort(mm,jj,iorb) ! USES 3-IDX TENSOR - exchange_int = three_e_3_idx_exch12_bi_ort(mm,jj,iorb) ! USES 3-IDX TENSOR - hthree += direct_int - exchange_int + !! same-spin/oposite-spin + do j = 1, na + jj = occ(j,ispin) + do m = 1, nb + mm = occ(m,other_spin) + direct_int = three_e_3_idx_direct_bi_ort(mm,jj,iorb) ! USES 3-IDX TENSOR + exchange_int = three_e_3_idx_exch12_bi_ort(mm,jj,iorb) ! USES 3-IDX TENSOR + hthree += direct_int - exchange_int + enddo enddo - enddo - !! oposite-spin/opposite-spin + !! oposite-spin/opposite-spin do j = 1, nb - jj = occ(j,other_spin) - do m = j+1, nb - mm = occ(m,other_spin) - direct_int = three_e_3_idx_direct_bi_ort(mm,jj,iorb) ! USES 3-IDX TENSOR - exchange_int = three_e_3_idx_exch23_bi_ort(mm,jj,iorb) ! USES 3-IDX TENSOR - hthree += direct_int - exchange_int - enddo + jj = occ(j,other_spin) + do m = j+1, nb + mm = occ(m,other_spin) + direct_int = three_e_3_idx_direct_bi_ort(mm,jj,iorb) ! USES 3-IDX TENSOR + exchange_int = three_e_3_idx_exch23_bi_ort(mm,jj,iorb) ! USES 3-IDX TENSOR + hthree += direct_int - exchange_int + enddo enddo endif na = na+1 + end +! --- + subroutine a_tc_operator(iorb,ispin,key,hmono,htwoe,hthree,Nint,na,nb) use bitmasks implicit none diff --git a/src/tc_bi_ortho/slater_tc_slow.irp.f b/src/tc_bi_ortho/slater_tc_slow.irp.f index 083d4ba3..d78540e7 100644 --- a/src/tc_bi_ortho/slater_tc_slow.irp.f +++ b/src/tc_bi_ortho/slater_tc_slow.irp.f @@ -21,7 +21,7 @@ subroutine htilde_mu_mat_bi_ortho_tot_slow(key_j, key_i, Nint, htot) integer :: degree call get_excitation_degree(key_j, key_i, degree, Nint) - if(degree.gt.2)then + if(degree.gt.2) then htot = 0.d0 else call htilde_mu_mat_bi_ortho_slow(key_j, key_i, Nint, hmono, htwoe, hthree, htot) @@ -60,22 +60,22 @@ subroutine htilde_mu_mat_bi_ortho_slow(key_j, key_i, Nint, hmono, htwoe, hthree, call get_excitation_degree(key_i, key_j, degree, Nint) if(degree.gt.2) return - if(degree == 0)then + if(degree == 0) then call diag_htilde_mu_mat_bi_ortho_slow(Nint, key_i, hmono, htwoe, htot) - else if (degree == 1)then + else if (degree == 1) then call single_htilde_mu_mat_bi_ortho_slow(Nint, key_j, key_i, hmono, htwoe, htot) - else if(degree == 2)then + else if(degree == 2) then call double_htilde_mu_mat_bi_ortho_slow(Nint, key_j, key_i, hmono, htwoe, htot) endif if(three_body_h_tc) then if(degree == 2) then - if(.not.double_normal_ord.and.elec_num.gt.2.and.three_e_5_idx_term) then + if((.not.double_normal_ord) .and. (elec_num .gt. 2) .and. three_e_5_idx_term) then call double_htilde_three_body_ints_bi_ort_slow(Nint, key_j, key_i, hthree) endif - else if(degree == 1.and.elec_num.gt.2.and.three_e_4_idx_term) then + else if((degree == 1) .and. (elec_num .gt. 2) .and. three_e_4_idx_term) then call single_htilde_three_body_ints_bi_ort_slow(Nint, key_j, key_i, hthree) - else if(degree == 0.and.elec_num.gt.2.and.three_e_3_idx_term) then + else if((degree == 0) .and. (elec_num .gt. 2) .and. three_e_3_idx_term) then call diag_htilde_three_body_ints_bi_ort_slow(Nint, key_i, hthree) endif endif diff --git a/src/tc_bi_ortho/tc_bi_ortho.irp.f b/src/tc_bi_ortho/tc_bi_ortho.irp.f index f69684c0..2887c7be 100644 --- a/src/tc_bi_ortho/tc_bi_ortho.irp.f +++ b/src/tc_bi_ortho/tc_bi_ortho.irp.f @@ -8,11 +8,13 @@ program tc_bi_ortho END_DOC my_grid_becke = .True. - my_n_pt_r_grid = 30 - my_n_pt_a_grid = 50 + PROVIDE tc_grid1_a tc_grid1_r + my_n_pt_r_grid = tc_grid1_r + my_n_pt_a_grid = tc_grid1_a + touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid + read_wf = .True. touch read_wf - touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid print*, ' nb of states = ', N_states print*, ' nb of det = ', N_det @@ -20,22 +22,29 @@ program tc_bi_ortho call routine_diag() call write_tc_energy() call save_tc_bi_ortho_wavefunction() -end - -subroutine test - implicit none - integer :: i,j - double precision :: hmono,htwoe,hthree,htot - use bitmasks - print*,'reading the wave function ' - do i = 1, N_det - call debug_det(psi_det(1,1,i),N_int) - print*,i,psi_l_coef_bi_ortho(i,1)*psi_r_coef_bi_ortho(i,1) - print*,i,psi_l_coef_bi_ortho(i,1),psi_r_coef_bi_ortho(i,1) - enddo end +! --- + +subroutine test() + + use bitmasks + implicit none + integer :: i, j + double precision :: hmono, htwoe, hthree, htot + + print*, 'reading the wave function ' + do i = 1, N_det + call debug_det(psi_det(1,1,i), N_int) + print*, i, psi_l_coef_bi_ortho(i,1)*psi_r_coef_bi_ortho(i,1) + print*, i, psi_l_coef_bi_ortho(i,1),psi_r_coef_bi_ortho(i,1) + enddo + +end + +! --- + subroutine routine_diag() implicit none diff --git a/src/tc_bi_ortho/tc_bi_ortho_prop.irp.f b/src/tc_bi_ortho/tc_bi_ortho_prop.irp.f index 28f122ee..9168fb3d 100644 --- a/src/tc_bi_ortho/tc_bi_ortho_prop.irp.f +++ b/src/tc_bi_ortho/tc_bi_ortho_prop.irp.f @@ -1,19 +1,32 @@ + +! --- + program tc_bi_ortho_prop - implicit none + BEGIN_DOC -! TODO : Put the documentation of the program here + ! TODO : Put the documentation of the program here END_DOC + + implicit none + print *, 'Hello world' + my_grid_becke = .True. - my_n_pt_r_grid = 30 - my_n_pt_a_grid = 50 + PROVIDE tc_grid1_a tc_grid1_r + my_n_pt_r_grid = tc_grid1_r + my_n_pt_a_grid = tc_grid1_a + touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid + read_wf = .True. touch read_wf - touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid -! call routine_diag - call test + + !call routine_diag + call test + end +! --- + subroutine test implicit none integer :: i diff --git a/src/tc_bi_ortho/tc_cisd_sc2.irp.f b/src/tc_bi_ortho/tc_cisd_sc2.irp.f index 0fb9f524..d4c8c55d 100644 --- a/src/tc_bi_ortho/tc_cisd_sc2.irp.f +++ b/src/tc_bi_ortho/tc_cisd_sc2.irp.f @@ -1,20 +1,32 @@ -program tc_bi_ortho - implicit none + +! --- + +program tc_cisd_sc2 + BEGIN_DOC -! TODO : Put the documentation of the program here + ! TODO : Put the documentation of the program here END_DOC + + implicit none + print *, 'Hello world' + my_grid_becke = .True. - my_n_pt_r_grid = 30 - my_n_pt_a_grid = 50 + PROVIDE tc_grid1_a tc_grid1_r + my_n_pt_r_grid = tc_grid1_r + my_n_pt_a_grid = tc_grid1_a + touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid + read_wf = .True. touch read_wf - touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid call test + end -subroutine test +! --- + +subroutine test() implicit none ! double precision, allocatable :: dressing_dets(:),e_corr_dets(:) ! allocate(dressing_dets(N_det),e_corr_dets(N_det)) diff --git a/src/tc_bi_ortho/tc_h_eigvectors.irp.f b/src/tc_bi_ortho/tc_h_eigvectors.irp.f index fa946d6a..f027c38f 100644 --- a/src/tc_bi_ortho/tc_h_eigvectors.irp.f +++ b/src/tc_bi_ortho/tc_h_eigvectors.irp.f @@ -1,42 +1,56 @@ + +! --- + use bitmasks - BEGIN_PROVIDER [ integer, index_HF_psi_det] - implicit none - integer :: i,degree - do i = 1, N_det - call get_excitation_degree(HF_bitmask,psi_det(1,1,i),degree,N_int) - if(degree == 0)then - index_HF_psi_det = i - exit - endif - enddo - END_PROVIDER +! --- + +BEGIN_PROVIDER [integer, index_HF_psi_det] -subroutine diagonalize_CI_tc implicit none + integer :: i, degree + + do i = 1, N_det + call get_excitation_degree(HF_bitmask, psi_det(1,1,i), degree, N_int) + if(degree == 0) then + index_HF_psi_det = i + exit + endif + enddo + +END_PROVIDER + +! --- + +subroutine diagonalize_CI_tc() + BEGIN_DOC -! Replace the coefficients of the |CI| states by the coefficients of the -! eigenstates of the |CI| matrix. + ! Replace the coefficients of the |CI| states by the coefficients of the + ! eigenstates of the |CI| matrix. END_DOC - integer :: i,j - do j=1,N_states - do i=1,N_det + + implicit none + integer :: i, j + + do j = 1, N_states + do i = 1, N_det psi_l_coef_bi_ortho(i,j) = leigvec_tc_bi_orth(i,j) psi_r_coef_bi_ortho(i,j) = reigvec_tc_bi_orth(i,j) enddo enddo SOFT_TOUCH psi_l_coef_bi_ortho psi_r_coef_bi_ortho + end +! --- - - BEGIN_PROVIDER [double precision, eigval_right_tc_bi_orth, (N_states)] -&BEGIN_PROVIDER [double precision, eigval_left_tc_bi_orth, (N_states)] -&BEGIN_PROVIDER [double precision, reigvec_tc_bi_orth, (N_det,N_states)] -&BEGIN_PROVIDER [double precision, leigvec_tc_bi_orth, (N_det,N_states)] -&BEGIN_PROVIDER [double precision, s2_eigvec_tc_bi_orth, (N_states)] -&BEGIN_PROVIDER [double precision, norm_ground_left_right_bi_orth ] + BEGIN_PROVIDER [double precision, eigval_right_tc_bi_orth, (N_states) ] +&BEGIN_PROVIDER [double precision, eigval_left_tc_bi_orth , (N_states) ] +&BEGIN_PROVIDER [double precision, reigvec_tc_bi_orth , (N_det,N_states)] +&BEGIN_PROVIDER [double precision, leigvec_tc_bi_orth , (N_det,N_states)] +&BEGIN_PROVIDER [double precision, s2_eigvec_tc_bi_orth , (N_states) ] +&BEGIN_PROVIDER [double precision, norm_ground_left_right_bi_orth ] BEGIN_DOC ! eigenvalues, right and left eigenvectors of the transcorrelated Hamiltonian on the BI-ORTHO basis @@ -44,29 +58,29 @@ end implicit none integer :: i, idx_dress, j, istate, k + integer :: i_good_state, i_other_state, i_state + integer :: n_real_tc_bi_orth_eigval_right, igood_r, igood_l logical :: converged, dagger - integer :: n_real_tc_bi_orth_eigval_right,igood_r,igood_l - double precision, allocatable :: reigvec_tc_bi_orth_tmp(:,:),leigvec_tc_bi_orth_tmp(:,:),eigval_right_tmp(:) + double precision, parameter :: alpha = 0.1d0 + integer, allocatable :: index_good_state_array(:) + integer, allocatable :: iorder(:) + logical, allocatable :: good_state_array(:) + double precision, allocatable :: reigvec_tc_bi_orth_tmp(:,:), leigvec_tc_bi_orth_tmp(:,:),eigval_right_tmp(:) double precision, allocatable :: s2_values_tmp(:), H_prime(:,:), expect_e(:) - double precision, parameter :: alpha = 0.1d0 - integer :: i_good_state,i_other_state, i_state - integer, allocatable :: index_good_state_array(:) - logical, allocatable :: good_state_array(:) - double precision, allocatable :: coef_hf_r(:),coef_hf_l(:) - double precision, allocatable :: Stmp(:,:) - integer, allocatable :: iorder(:) + double precision, allocatable :: coef_hf_r(:),coef_hf_l(:) + double precision, allocatable :: Stmp(:,:) PROVIDE N_det N_int - if(n_det .le. N_det_max_full) then + if(N_det .le. N_det_max_full) then - allocate(reigvec_tc_bi_orth_tmp(N_det,N_det),leigvec_tc_bi_orth_tmp(N_det,N_det),eigval_right_tmp(N_det),expect_e(N_det)) - allocate (H_prime(N_det,N_det),s2_values_tmp(N_det)) + allocate(reigvec_tc_bi_orth_tmp(N_det,N_det), leigvec_tc_bi_orth_tmp(N_det,N_det), eigval_right_tmp(N_det), expect_e(N_det)) + allocate(H_prime(N_det,N_det), s2_values_tmp(N_det)) H_prime(1:N_det,1:N_det) = htilde_matrix_elmt_bi_ortho(1:N_det,1:N_det) if(s2_eig) then H_prime(1:N_det,1:N_det) += alpha * S2_matrix_all_dets(1:N_det,1:N_det) - do j=1,N_det + do j = 1, N_det H_prime(j,j) = H_prime(j,j) - alpha*expected_s2 enddo endif diff --git a/src/tc_bi_ortho/tc_hmat.irp.f b/src/tc_bi_ortho/tc_hmat.irp.f index ec072531..4c21b5bd 100644 --- a/src/tc_bi_ortho/tc_hmat.irp.f +++ b/src/tc_bi_ortho/tc_hmat.irp.f @@ -31,7 +31,9 @@ END_PROVIDER - BEGIN_PROVIDER [double precision, htilde_matrix_elmt_bi_ortho_tranp, (N_det,N_det)] +! --- + +BEGIN_PROVIDER [double precision, htilde_matrix_elmt_bi_ortho_tranp, (N_det,N_det)] implicit none integer ::i,j do i = 1, N_det diff --git a/src/tc_bi_ortho/tc_som.irp.f b/src/tc_bi_ortho/tc_som.irp.f index a7e4d09e..427508d2 100644 --- a/src/tc_bi_ortho/tc_som.irp.f +++ b/src/tc_bi_ortho/tc_som.irp.f @@ -12,10 +12,9 @@ program tc_som print *, ' do not forget to do tc-scf first' my_grid_becke = .True. - my_n_pt_r_grid = 30 - my_n_pt_a_grid = 50 -! my_n_pt_r_grid = 10 ! small grid for quick debug -! my_n_pt_a_grid = 26 ! small grid for quick debug + PROVIDE tc_grid1_a tc_grid1_r + my_n_pt_r_grid = tc_grid1_r + my_n_pt_a_grid = tc_grid1_a touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid PROVIDE mu_erf diff --git a/src/tc_bi_ortho/test_natorb.irp.f b/src/tc_bi_ortho/test_natorb.irp.f index 54c9a827..5b8801f7 100644 --- a/src/tc_bi_ortho/test_natorb.irp.f +++ b/src/tc_bi_ortho/test_natorb.irp.f @@ -1,21 +1,34 @@ + +! --- + program test_natorb - implicit none + BEGIN_DOC -! TODO : Reads psi_det in the EZFIO folder and prints out the left- and right-eigenvectors together with the energy. Saves the left-right wave functions at the end. + ! TODO : Reads psi_det in the EZFIO folder and prints out the left- and right-eigenvectors together with the energy. Saves the left-right wave functions at the end. END_DOC + + implicit none + print *, 'Hello world' + my_grid_becke = .True. - my_n_pt_r_grid = 30 - my_n_pt_a_grid = 50 + PROVIDE tc_grid1_a tc_grid1_r + my_n_pt_r_grid = tc_grid1_r + my_n_pt_a_grid = tc_grid1_a + touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid + read_wf = .True. touch read_wf - touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid - call routine -! call test + + call routine() + ! call test() end -subroutine routine +! --- + +subroutine routine() + implicit none double precision, allocatable :: fock_diag(:),eigval(:),leigvec(:,:),reigvec(:,:),mat_ref(:,:) allocate(eigval(mo_num),leigvec(mo_num,mo_num),reigvec(mo_num,mo_num),fock_diag(mo_num),mat_ref(mo_num, mo_num)) diff --git a/src/tc_bi_ortho/test_normal_order.irp.f b/src/tc_bi_ortho/test_normal_order.irp.f index cb0c355c..e805eecb 100644 --- a/src/tc_bi_ortho/test_normal_order.irp.f +++ b/src/tc_bi_ortho/test_normal_order.irp.f @@ -1,19 +1,32 @@ + +! --- + program test_normal_order - implicit none + BEGIN_DOC -! TODO : Put the documentation of the program here + ! TODO : Put the documentation of the program here END_DOC + + implicit none + print *, 'Hello world' + my_grid_becke = .True. - my_n_pt_r_grid = 30 - my_n_pt_a_grid = 50 + PROVIDE tc_grid1_a tc_grid1_r + my_n_pt_r_grid = tc_grid1_r + my_n_pt_a_grid = tc_grid1_a + touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid + read_wf = .True. touch read_wf - touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid - call provide_all_three_ints_bi_ortho - call test + + call provide_all_three_ints_bi_ortho() + call test() + end +! --- + subroutine test implicit none use bitmasks ! you need to include the bitmasks_module.f90 features diff --git a/src/tc_bi_ortho/test_s2_tc.irp.f b/src/tc_bi_ortho/test_s2_tc.irp.f index 1f7bdfda..b398507a 100644 --- a/src/tc_bi_ortho/test_s2_tc.irp.f +++ b/src/tc_bi_ortho/test_s2_tc.irp.f @@ -1,14 +1,22 @@ + +! --- + program test_tc - implicit none - read_wf = .True. - my_grid_becke = .True. - my_n_pt_r_grid = 30 - my_n_pt_a_grid = 50 - read_wf = .True. - touch read_wf - touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid - call routine_test_s2 - call routine_test_s2_davidson + + implicit none + + my_grid_becke = .True. + PROVIDE tc_grid1_a tc_grid1_r + my_n_pt_r_grid = tc_grid1_r + my_n_pt_a_grid = tc_grid1_a + touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid + + read_wf = .True. + touch read_wf + + call routine_test_s2 + call routine_test_s2_davidson + end subroutine routine_test_s2 diff --git a/src/tc_bi_ortho/test_tc_bi_ortho.irp.f b/src/tc_bi_ortho/test_tc_bi_ortho.irp.f index 902f7295..b6beb65b 100644 --- a/src/tc_bi_ortho/test_tc_bi_ortho.irp.f +++ b/src/tc_bi_ortho/test_tc_bi_ortho.irp.f @@ -1,15 +1,24 @@ + +! --- + program tc_bi_ortho - implicit none + BEGIN_DOC -! TODO : Put the documentation of the program here + ! TODO : Put the documentation of the program here END_DOC + + implicit none + print *, 'Hello world' + my_grid_becke = .True. - my_n_pt_r_grid = 30 - my_n_pt_a_grid = 50 + PROVIDE tc_grid1_a tc_grid1_r + my_n_pt_r_grid = tc_grid1_r + my_n_pt_a_grid = tc_grid1_a + touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid + read_wf = .True. touch read_wf - touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid ! call test_h_u0 ! call test_slater_tc_opt diff --git a/src/tc_bi_ortho/test_tc_fock.irp.f b/src/tc_bi_ortho/test_tc_fock.irp.f index b7de067f..182c03d7 100644 --- a/src/tc_bi_ortho/test_tc_fock.irp.f +++ b/src/tc_bi_ortho/test_tc_fock.irp.f @@ -1,22 +1,32 @@ + +! --- + program test_tc_fock - implicit none + BEGIN_DOC -! TODO : Put the documentation of the program here + ! TODO : Put the documentation of the program here END_DOC + + implicit none + print *, 'Hello world' + my_grid_becke = .True. - my_n_pt_r_grid = 30 - my_n_pt_a_grid = 50 + PROVIDE tc_grid1_a tc_grid1_r + my_n_pt_r_grid = tc_grid1_r + my_n_pt_a_grid = tc_grid1_a + touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid + read_wf = .True. touch read_wf - touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid !call routine_1 !call routine_2 ! call routine_3() ! call test_3e - call routine_tot + call routine_tot + end ! --- diff --git a/src/tc_keywords/EZFIO.cfg b/src/tc_keywords/EZFIO.cfg index ea1503c3..29c238e1 100644 --- a/src/tc_keywords/EZFIO.cfg +++ b/src/tc_keywords/EZFIO.cfg @@ -262,3 +262,16 @@ doc: If |true|, use Manu IPP interface: ezfio,provider,ocaml default: True +[tc_grid1_a] +type: integer +doc: size of angular grid over r1 +interface: ezfio,provider,ocaml +default: 50 + +[tc_grid1_r] +type: integer +doc: size of radial grid over r1 +interface: ezfio,provider,ocaml +default: 30 + + diff --git a/src/tc_scf/combine_lr_tcscf.irp.f b/src/tc_scf/combine_lr_tcscf.irp.f index b257f4a5..a22614ba 100644 --- a/src/tc_scf/combine_lr_tcscf.irp.f +++ b/src/tc_scf/combine_lr_tcscf.irp.f @@ -10,8 +10,9 @@ program combine_lr_tcscf implicit none my_grid_becke = .True. - my_n_pt_r_grid = 30 - my_n_pt_a_grid = 50 + PROVIDE tc_grid1_a tc_grid1_r + my_n_pt_r_grid = tc_grid1_r + my_n_pt_a_grid = tc_grid1_a touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid bi_ortho = .True. diff --git a/src/tc_scf/minimize_tc_angles.irp.f b/src/tc_scf/minimize_tc_angles.irp.f index 5d3ff7f0..c7752930 100644 --- a/src/tc_scf/minimize_tc_angles.irp.f +++ b/src/tc_scf/minimize_tc_angles.irp.f @@ -1,17 +1,26 @@ -program print_angles - implicit none - BEGIN_DOC -! program that minimizes the angle between left- and right-orbitals when degeneracies are found in the TC-Fock matrix - END_DOC + +! --- + +program minimize_tc_angles + + BEGIN_DOC + ! program that minimizes the angle between left- and right-orbitals when degeneracies are found in the TC-Fock matrix + END_DOC + + implicit none + my_grid_becke = .True. - my_n_pt_r_grid = 30 - my_n_pt_a_grid = 50 + PROVIDE tc_grid1_a tc_grid1_r + my_n_pt_r_grid = tc_grid1_r + my_n_pt_a_grid = tc_grid1_a touch my_n_pt_r_grid my_n_pt_a_grid -! call sort_by_tc_fock + + ! call sort_by_tc_fock ! TODO ! check if rotations of orbitals affect the TC energy ! and refuse the step call minimize_tc_orb_angles + end diff --git a/src/tc_scf/molden_lr_mos.irp.f b/src/tc_scf/molden_lr_mos.irp.f index e12fcd1c..b86009ee 100644 --- a/src/tc_scf/molden_lr_mos.irp.f +++ b/src/tc_scf/molden_lr_mos.irp.f @@ -11,10 +11,9 @@ program molden_lr_mos print *, 'starting ...' my_grid_becke = .True. - my_n_pt_r_grid = 30 - my_n_pt_a_grid = 50 -! my_n_pt_r_grid = 10 ! small grid for quick debug -! my_n_pt_a_grid = 26 ! small grid for quick debug + PROVIDE tc_grid1_a tc_grid1_r + my_n_pt_r_grid = tc_grid1_r + my_n_pt_a_grid = tc_grid1_a touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid !call molden_lr diff --git a/src/tc_scf/print_fit_param.irp.f b/src/tc_scf/print_fit_param.irp.f index f8bcfa7f..e62f0dde 100644 --- a/src/tc_scf/print_fit_param.irp.f +++ b/src/tc_scf/print_fit_param.irp.f @@ -7,10 +7,9 @@ program print_fit_param implicit none my_grid_becke = .True. - my_n_pt_r_grid = 30 - my_n_pt_a_grid = 50 -! my_n_pt_r_grid = 10 ! small grid for quick debug -! my_n_pt_a_grid = 26 ! small grid for quick debug + PROVIDE tc_grid1_a tc_grid1_r + my_n_pt_r_grid = tc_grid1_r + my_n_pt_a_grid = tc_grid1_a touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid !call create_guess diff --git a/src/tc_scf/print_tcscf_energy.irp.f b/src/tc_scf/print_tcscf_energy.irp.f index 96512cb7..05b8df23 100644 --- a/src/tc_scf/print_tcscf_energy.irp.f +++ b/src/tc_scf/print_tcscf_energy.irp.f @@ -8,16 +8,9 @@ program print_tcscf_energy print *, 'Hello world' my_grid_becke = .True. - - !my_n_pt_r_grid = 30 - !my_n_pt_a_grid = 50 - - my_n_pt_r_grid = 100 - my_n_pt_a_grid = 170 - - !my_n_pt_r_grid = 100 - !my_n_pt_a_grid = 266 - + PROVIDE tc_grid1_a tc_grid1_r + my_n_pt_r_grid = tc_grid1_r + my_n_pt_a_grid = tc_grid1_a touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid call main() diff --git a/src/tc_scf/rotate_tcscf_orbitals.irp.f b/src/tc_scf/rotate_tcscf_orbitals.irp.f index 2567faf0..0f2663e5 100644 --- a/src/tc_scf/rotate_tcscf_orbitals.irp.f +++ b/src/tc_scf/rotate_tcscf_orbitals.irp.f @@ -10,8 +10,9 @@ program rotate_tcscf_orbitals implicit none my_grid_becke = .True. - my_n_pt_r_grid = 30 - my_n_pt_a_grid = 50 + PROVIDE tc_grid1_a tc_grid1_r + my_n_pt_r_grid = tc_grid1_r + my_n_pt_a_grid = tc_grid1_a touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid bi_ortho = .True. diff --git a/src/tc_scf/tc_petermann_factor.irp.f b/src/tc_scf/tc_petermann_factor.irp.f index d3722098..2e9c67e2 100644 --- a/src/tc_scf/tc_petermann_factor.irp.f +++ b/src/tc_scf/tc_petermann_factor.irp.f @@ -10,10 +10,9 @@ program tc_petermann_factor implicit none my_grid_becke = .True. - my_n_pt_r_grid = 30 - my_n_pt_a_grid = 50 -! my_n_pt_r_grid = 10 ! small grid for quick debug -! my_n_pt_a_grid = 26 ! small grid for quick debug + PROVIDE tc_grid1_a tc_grid1_r + my_n_pt_r_grid = tc_grid1_r + my_n_pt_a_grid = tc_grid1_a touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid call main() diff --git a/src/tc_scf/tc_scf.irp.f b/src/tc_scf/tc_scf.irp.f index 2f2d803f..e4c38741 100644 --- a/src/tc_scf/tc_scf.irp.f +++ b/src/tc_scf/tc_scf.irp.f @@ -14,14 +14,10 @@ program tc_scf my_grid_becke = .True. - !my_n_pt_r_grid = 30 - !my_n_pt_a_grid = 50 + PROVIDE tc_grid1_a tc_grid1_r + my_n_pt_r_grid = tc_grid1_r + my_n_pt_a_grid = tc_grid1_a - my_n_pt_r_grid = 100 - my_n_pt_a_grid = 170 - -! my_n_pt_r_grid = 10 ! small grid for quick debug -! my_n_pt_a_grid = 26 ! small grid for quick debug touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid PROVIDE mu_erf diff --git a/src/tc_scf/test_int.irp.f b/src/tc_scf/test_int.irp.f index b9287d58..649d0f3e 100644 --- a/src/tc_scf/test_int.irp.f +++ b/src/tc_scf/test_int.irp.f @@ -9,10 +9,9 @@ program test_ints print *, ' starting test_ints ...' my_grid_becke = .True. - my_n_pt_r_grid = 30 - my_n_pt_a_grid = 50 -! my_n_pt_r_grid = 15 ! small grid for quick debug -! my_n_pt_a_grid = 26 ! small grid for quick debug + PROVIDE tc_grid1_a tc_grid1_r + my_n_pt_r_grid = tc_grid1_r + my_n_pt_a_grid = tc_grid1_a touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid my_extra_grid_becke = .True. @@ -280,7 +279,7 @@ subroutine routine_v_ij_u_cst_mu_j1b_test do i = 1, ao_num do j = 1, ao_num array(j,i,l,k) += v_ij_u_cst_mu_j1b_test(j,i,ipoint) * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight - array_ref(j,i,l,k) += v_ij_u_cst_mu_j1b(j,i,ipoint) * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight + array_ref(j,i,l,k) += v_ij_u_cst_mu_j1b_fit (j,i,ipoint) * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight enddo enddo enddo @@ -506,7 +505,7 @@ subroutine routine_v_ij_u_cst_mu_j1b do i = 1, ao_num do j = 1, ao_num array(j,i,l,k) += v_ij_u_cst_mu_j1b_test(j,i,ipoint) * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight - array_ref(j,i,l,k) += v_ij_u_cst_mu_j1b(j,i,ipoint) * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight + array_ref(j,i,l,k) += v_ij_u_cst_mu_j1b_fit (j,i,ipoint) * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight enddo enddo enddo