diff --git a/ocaml/Molecule.ml b/ocaml/Molecule.ml index 68148288..27ac4bd0 100644 --- a/ocaml/Molecule.ml +++ b/ocaml/Molecule.ml @@ -157,7 +157,7 @@ let of_xyz_file | _ -> raise XYZError end; String.concat "\n" rest - | _ -> failwith ("Problem in xyz file "^filename) + | _ -> raise XYZError in of_xyz_string ~charge:charge ~multiplicity:multiplicity ~units:units lines diff --git a/plugins/MRPT/jm_mrpt2.irp.f b/plugins/MRPT/jm_mrpt2.irp.f index 89deef2e..df9a1dbc 100644 --- a/plugins/MRPT/jm_mrpt2.irp.f +++ b/plugins/MRPT/jm_mrpt2.irp.f @@ -11,7 +11,7 @@ end subroutine routine_3 implicit none !provide fock_virt_total_spin_trace - provide delta_ij + provide delta_ij_mrpt print *, 'N_det = ', N_det print *, 'N_states = ', N_states diff --git a/plugins/MRPT_Utils/mrpt_utils.irp.f b/plugins/MRPT_Utils/mrpt_utils.irp.f index 83f087bc..34d26127 100644 --- a/plugins/MRPT_Utils/mrpt_utils.irp.f +++ b/plugins/MRPT_Utils/mrpt_utils.irp.f @@ -1,5 +1,5 @@ - BEGIN_PROVIDER [ double precision, delta_ij, (N_det,N_det,N_states) ] + BEGIN_PROVIDER [ double precision, delta_ij_mrpt, (N_det,N_det,N_states) ] &BEGIN_PROVIDER [ double precision, second_order_pt_new, (N_states) ] &BEGIN_PROVIDER [ double precision, second_order_pt_new_1h, (N_states) ] &BEGIN_PROVIDER [ double precision, second_order_pt_new_1p, (N_states) ] @@ -19,7 +19,7 @@ double precision, allocatable :: delta_ij_tmp(:,:,:) - delta_ij = 0.d0 + delta_ij_mrpt = 0.d0 allocate (delta_ij_tmp(N_det,N_det,N_states)) @@ -32,7 +32,7 @@ do i = 1, N_det do j = 1, N_det accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state) - delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state) + delta_ij_mrpt(j,i,i_state) += delta_ij_tmp(j,i,i_state) enddo enddo second_order_pt_new_1h(i_state) = accu(i_state) @@ -47,7 +47,7 @@ do i = 1, N_det do j = 1, N_det accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state) - delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state) + delta_ij_mrpt(j,i,i_state) += delta_ij_tmp(j,i,i_state) enddo enddo second_order_pt_new_1p(i_state) = accu(i_state) @@ -63,7 +63,7 @@ do i = 1, N_det do j = 1, N_det accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state) - delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state) + delta_ij_mrpt(j,i,i_state) += delta_ij_tmp(j,i,i_state) enddo enddo second_order_pt_new_1h1p(i_state) = accu(i_state) @@ -79,7 +79,7 @@ do i = 1, N_det do j = 1, N_det accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state) - delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state) + delta_ij_mrpt(j,i,i_state) += delta_ij_tmp(j,i,i_state) enddo enddo second_order_pt_new_1h1p(i_state) = accu(i_state) @@ -95,7 +95,7 @@ do i = 1, N_det do j = 1, N_det accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state) - delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state) + delta_ij_mrpt(j,i,i_state) += delta_ij_tmp(j,i,i_state) enddo enddo second_order_pt_new_2h(i_state) = accu(i_state) @@ -110,7 +110,7 @@ do i = 1, N_det do j = 1, N_det accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state) - delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state) + delta_ij_mrpt(j,i,i_state) += delta_ij_tmp(j,i,i_state) enddo enddo second_order_pt_new_2p(i_state) = accu(i_state) @@ -126,7 +126,7 @@ do i = 1, N_det do j = 1, N_det accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state) - delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state) + delta_ij_mrpt(j,i,i_state) += delta_ij_tmp(j,i,i_state) enddo enddo second_order_pt_new_1h2p(i_state) = accu(i_state) @@ -142,7 +142,7 @@ do i = 1, N_det do j = 1, N_det accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state) - delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state) + delta_ij_mrpt(j,i,i_state) += delta_ij_tmp(j,i,i_state) enddo enddo second_order_pt_new_2h1p(i_state) = accu(i_state) @@ -157,7 +157,7 @@ !do i = 1, N_det ! do j = 1, N_det ! accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state) -! delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state) +! delta_ij_mrpt(j,i,i_state) += delta_ij_tmp(j,i,i_state) ! enddo !enddo !second_order_pt_new_2h2p(i_state) = accu(i_state) @@ -168,7 +168,7 @@ call give_2h2p(contrib_2h2p) do i_state = 1, N_states do i = 1, N_det - delta_ij(i,i,i_state) += contrib_2h2p(i_state) + delta_ij_mrpt(i,i,i_state) += contrib_2h2p(i_state) enddo second_order_pt_new_2h2p(i_state) = contrib_2h2p(i_state) enddo @@ -179,9 +179,9 @@ accu = 0.d0 do i_state = 1, N_states do i = 1, N_det -! write(*,'(1000(F16.10,x))')delta_ij(i,:,:) +! write(*,'(1000(F16.10,x))')delta_ij_mrpt(i,:,:) do j = i_state, N_det - accu(i_state) += delta_ij(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state) + accu(i_state) += delta_ij_mrpt(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state) enddo enddo second_order_pt_new(i_state) = accu(i_state) @@ -199,7 +199,7 @@ END_PROVIDER do i_state = 1, N_states do i = 1,N_det do j = 1,N_det - Hmatrix_dressed_pt2_new(j,i,i_state) = H_matrix_all_dets(j,i) + delta_ij(j,i,i_state) + Hmatrix_dressed_pt2_new(j,i,i_state) = H_matrix_all_dets(j,i) + delta_ij_mrpt(j,i,i_state) enddo enddo enddo @@ -214,7 +214,7 @@ END_PROVIDER do i = 1,N_det do j = i,N_det Hmatrix_dressed_pt2_new_symmetrized(j,i,i_state) = H_matrix_all_dets(j,i) & - + 0.5d0 * ( delta_ij(j,i,i_state) + delta_ij(i,j,i_state) ) + + 0.5d0 * ( delta_ij_mrpt(j,i,i_state) + delta_ij_mrpt(i,j,i_state) ) Hmatrix_dressed_pt2_new_symmetrized(i,j,i_state) = Hmatrix_dressed_pt2_new_symmetrized(j,i,i_state) enddo enddo diff --git a/src/Determinants/EZFIO.cfg b/src/Determinants/EZFIO.cfg index 9d0512f4..6bf6faff 100644 --- a/src/Determinants/EZFIO.cfg +++ b/src/Determinants/EZFIO.cfg @@ -2,7 +2,7 @@ type: Det_number_max doc: Max number of determinants in the wave function interface: ezfio,provider,ocaml -default: 10000 +default: 1000000 [N_det_max_property] type: Det_number_max @@ -14,7 +14,7 @@ default: 10000 type: Det_number_max doc: Maximum number of determinants diagonalized by Jacobi interface: ezfio,provider,ocaml -default: 1000 +default: 2000 [N_states] type: States_number @@ -50,7 +50,7 @@ default: 0.99 type: Threshold doc: Thresholds on selectors (fraction of the norm) interface: ezfio,provider,ocaml -default: 0.999 +default: 0.995 [n_int] interface: ezfio diff --git a/src/MO_Basis/mos.irp.f b/src/MO_Basis/mos.irp.f index 8962dd00..09bde5c0 100644 --- a/src/MO_Basis/mos.irp.f +++ b/src/MO_Basis/mos.irp.f @@ -153,32 +153,39 @@ BEGIN_PROVIDER [ double precision, mo_occ, (mo_tot_num) ] endif END_PROVIDER + subroutine ao_to_mo(A_ao,LDA_ao,A_mo,LDA_mo) implicit none BEGIN_DOC ! Transform A from the AO basis to the MO basis ! - ! C.A_ao.Ct + ! (S.C)t.A_ao.S.C END_DOC integer, intent(in) :: LDA_ao,LDA_mo double precision, intent(in) :: A_ao(LDA_ao) double precision, intent(out) :: A_mo(LDA_mo) - double precision, allocatable :: T(:,:) + double precision, allocatable :: T(:,:), SC(:,:) + allocate ( SC(ao_num_align,mo_tot_num) ) allocate ( T(ao_num_align,mo_tot_num) ) !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: T - call dgemm('N','N', ao_num, mo_tot_num, ao_num, & - 1.d0, A_ao,LDA_ao, & + call dgemm('N','N', ao_num, ao_num, mo_tot_num, & + 1.d0, ao_overlap,size(ao_overlap,1), & mo_coef, size(mo_coef,1), & - 0.d0, T, ao_num_align) + 0.d0, SC, ao_num_align) - call dgemm('T','N', mo_tot_num, mo_tot_num, ao_num, & - 1.d0, mo_coef,size(mo_coef,1), & - T, ao_num_align, & - 0.d0, A_mo, LDA_mo) + call dgemm('T','N', ao_num, mo_tot_num, ao_num, & + 1.d0, SC, size(SC,1), & + A_ao, size(A_ao,1), & + 0.d0, T, size(T,1)) - deallocate(T) + call dgemm('N','N', ao_num, mo_tot_num, ao_num, & + 1.d0, T,size(T,1), & + SC, size(SC,1), & + 0.d0, A_mo, size(A_mo,1)) + + deallocate(T,SC) end subroutine mo_to_ao(A_mo,LDA_mo,A_ao,LDA_ao) @@ -186,39 +193,7 @@ subroutine mo_to_ao(A_mo,LDA_mo,A_ao,LDA_ao) BEGIN_DOC ! Transform A from the MO basis to the AO basis ! - ! (S.C).A_mo.(S.C)t - END_DOC - integer, intent(in) :: LDA_ao,LDA_mo - double precision, intent(in) :: A_mo(LDA_mo) - double precision, intent(out) :: A_ao(LDA_ao) - double precision, allocatable :: T(:,:), SC(:,:) - - allocate ( SC(ao_num_align,mo_tot_num) ) - allocate ( T(mo_tot_num_align,ao_num) ) - !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: T - - call dgemm('N','N', ao_num, mo_tot_num, ao_num, & - 1.d0, ao_overlap,size(ao_overlap,1), & - mo_coef, size(mo_coef,1), & - 0.d0, SC, ao_num_align) - - call dgemm('N','T', mo_tot_num, ao_num, mo_tot_num, & - 1.d0, A_mo,LDA_mo, & - SC, size(SC,1), & - 0.d0, T, mo_tot_num_align) - - call dgemm('N','N', ao_num, ao_num, mo_tot_num, & - 1.d0, SC,size(SC,1), & - T, mo_tot_num_align, & - 0.d0, A_ao, LDA_ao) - - deallocate(T,SC) -end - -subroutine mo_to_ao_no_overlap(A_mo,LDA_mo,A_ao,LDA_ao) - implicit none - BEGIN_DOC - ! Transform A from the MO basis to the S^-1 AO basis + ! C.A_mo.Ct END_DOC integer, intent(in) :: LDA_ao,LDA_mo double precision, intent(in) :: A_mo(LDA_mo) @@ -231,16 +206,17 @@ subroutine mo_to_ao_no_overlap(A_mo,LDA_mo,A_ao,LDA_ao) call dgemm('N','T', mo_tot_num, ao_num, mo_tot_num, & 1.d0, A_mo,LDA_mo, & mo_coef, size(mo_coef,1), & - 0.d0, T, mo_tot_num_align) + 0.d0, T, size(T,1)) call dgemm('N','N', ao_num, ao_num, mo_tot_num, & 1.d0, mo_coef,size(mo_coef,1), & - T, mo_tot_num_align, & - 0.d0, A_ao, LDA_ao) + T, size(T,1), & + 0.d0, A_ao, size(A_ao,1)) deallocate(T) end + subroutine mix_mo_jk(j,k) implicit none integer, intent(in) :: j,k