diff --git a/plugins/local/tc_scf/routines_rotates.irp.f b/plugins/local/tc_scf/routines_rotates.irp.f index cc825429..c42e846e 100644 --- a/plugins/local/tc_scf/routines_rotates.irp.f +++ b/plugins/local/tc_scf/routines_rotates.irp.f @@ -103,7 +103,7 @@ subroutine routine_save_rotated_mos(thr_deg, good_angles) double precision, allocatable :: stmp(:,:), T(:,:), Snew(:,:), smat2(:,:) double precision, allocatable :: mo_l_coef_tmp(:,:), mo_r_coef_tmp(:,:), mo_l_coef_new(:,:) - E_thr = 1d-8 + E_thr = 1d-04 E_old = TC_HF_energy allocate(mo_l_coef_old(ao_num,mo_num), mo_r_coef_old(ao_num,mo_num)) mo_r_coef_old = mo_r_coef @@ -164,10 +164,42 @@ subroutine routine_save_rotated_mos(thr_deg, good_angles) allocate(mo_r_coef_tmp(ao_num,n_degen), mo_l_coef_tmp(ao_num,n_degen), mo_l_coef_new(ao_num,n_degen)) allocate(T(n_degen,n_degen), Snew(n_degen,n_degen)) - do j = 1, n_degen - mo_r_coef_tmp(1:ao_num,j) = mo_r_coef_new(1:ao_num,list_degen(i,j)) - mo_l_coef_tmp(1:ao_num,j) = mo_l_coef(1:ao_num,list_degen(i,j)) - enddo + print*,'Right orbitals before' + do j = 1, n_degen + write(*,'(100(F16.10,X))') mo_r_coef_new(1:ao_num,list_degen(i,j)) + enddo + print*,'Left orbitals before' + do j = 1, n_degen + write(*,'(100(F16.10,X))')mo_l_coef(1:ao_num,list_degen(i,j)) + enddo + if(angle_left_right(list_degen(i,1)).gt.80.d0.and.n_degen==2)then + integer :: i_list, j_list + i_list = list_degen(i,1) + j_list = list_degen(i,2) + print*,'Huge angle !!! == ',angle_left_right(list_degen(i,1)),angle_left_right(list_degen(i,2)) + print*,'i_list = ',i_list + print*,'i_list = ',j_list + print*,'Swapping left/right orbitals' + call print_strong_overlap(i_list, j_list) + mo_r_coef_tmp(1:ao_num,1) = mo_r_coef_new(1:ao_num,i_list) + mo_r_coef_tmp(1:ao_num,2) = mo_l_coef(1:ao_num,i_list) + mo_l_coef_tmp(1:ao_num,1) = mo_l_coef(1:ao_num,j_list) + mo_l_coef_tmp(1:ao_num,2) = mo_r_coef_new(1:ao_num,j_list) + else + do j = 1, n_degen + print*,'i_list = ',list_degen(i,j) + mo_r_coef_tmp(1:ao_num,j) = mo_r_coef_new(1:ao_num,list_degen(i,j)) + mo_l_coef_tmp(1:ao_num,j) = mo_l_coef(1:ao_num,list_degen(i,j)) + enddo + endif + print*,'Right orbitals ' + do j = 1, n_degen + write(*,'(100(F16.10,X))')mo_r_coef_tmp(1:ao_num,j) + enddo + print*,'Left orbitals ' + do j = 1, n_degen + write(*,'(100(F16.10,X))')mo_l_coef_tmp(1:ao_num,j) + enddo ! Orthogonalization of right functions print *, ' Orthogonalization of RIGHT functions' print *, ' ------------------------------------' @@ -445,3 +477,31 @@ subroutine sort_by_tc_fock end + +subroutine print_strong_overlap(i_list, j_list) + implicit none + integer, intent(in) :: i_list,j_list + double precision :: o_i, o_j,o_ij + double precision :: s_mat_r(2,2),s_mat_l(2,2) + o_i = dsqrt(overlap_mo_r(i_list, i_list)) + o_j = dsqrt(overlap_mo_r(j_list, j_list)) + o_ij = overlap_mo_r(j_list, i_list) + s_mat_r(1,1) = o_i*o_i + s_mat_r(2,1) = o_ij/(o_i * o_j) + s_mat_r(2,2) = o_j*o_j + s_mat_r(1,2) = s_mat_r(2,1) + print*,'Right overlap matrix ' + write(*,'(2(F10.5,X))')s_mat_r(1:2,1) + write(*,'(2(F10.5,X))')s_mat_r(1:2,2) + o_i = dsqrt(overlap_mo_l(i_list, i_list)) + o_j = dsqrt(overlap_mo_l(j_list, j_list)) + o_ij = overlap_mo_l(j_list, i_list) + s_mat_l(1,1) = o_i*o_i + s_mat_l(2,1) = o_ij/(o_i * o_j) + s_mat_l(2,2) = o_j*o_j + s_mat_l(1,2) = s_mat_l(2,1) + print*,'Left overlap matrix ' + write(*,'(2(F10.5,X))')s_mat_l(1:2,1) + write(*,'(2(F10.5,X))')s_mat_l(1:2,2) + +end diff --git a/src/mo_one_e_ints/spread_dipole_mo.irp.f b/src/mo_one_e_ints/spread_dipole_mo.irp.f index e4484433..b0a7198b 100644 --- a/src/mo_one_e_ints/spread_dipole_mo.irp.f +++ b/src/mo_one_e_ints/spread_dipole_mo.irp.f @@ -58,3 +58,21 @@ END_PROVIDER ) END_PROVIDER + BEGIN_PROVIDER [double precision, mo_spread_centered_x, (mo_num, mo_num) ] +&BEGIN_PROVIDER [double precision, mo_spread_centered_y, (mo_num, mo_num) ] +&BEGIN_PROVIDER [double precision, mo_spread_centered_z, (mo_num, mo_num) ] + BEGIN_DOC + ! array of the integrals of MO_i * (x^2 - ^2) MO_j = MO_i x^2 MO_j - (MO_i x MO_j)^2 + ! array of the integrals of MO_i * (y^2 - ^2) MO_j = MO_i y^2 MO_j - (MO_i y MO_j)^2 + ! array of the integrals of MO_i * (z^2 - ^2) MO_j = MO_i z^2 MO_j - (MO_i z MO_j)^2 + END_DOC + implicit none + integer :: i,j + do i = 1, mo_num + do j = 1, mo_num + mo_spread_centered_x(j,i) = mo_spread_x(j,i) - mo_dipole_x(j,i)**2 + mo_spread_centered_y(j,i) = mo_spread_y(j,i) - mo_dipole_y(j,i)**2 + mo_spread_centered_z(j,i) = mo_spread_z(j,i) - mo_dipole_z(j,i)**2 + enddo + enddo +END_PROVIDER diff --git a/src/utils/linear_algebra.irp.f b/src/utils/linear_algebra.irp.f index 175beff3..26e390b7 100644 --- a/src/utils/linear_algebra.irp.f +++ b/src/utils/linear_algebra.irp.f @@ -1920,8 +1920,12 @@ subroutine exp_matrix(X,n,exp_X) call get_A_squared(X,n,r2_mat) call lapack_diagd(eigvalues,eigvectors,r2_mat,n,n) eigvalues=-eigvalues + do i = 1,n + ! t = dsqrt(t^2) where t^2 are eigenvalues of X^2 + eigvalues(i) = dsqrt(eigvalues(i)) + enddo - if(.False.)then + if(.false.)then !!! For debugging and following the book intermediate ! rebuilding the matrix : X^2 = -W t^2 W^T as in 3.1.30 ! matrix_tmp1 = W t^2 @@ -1932,14 +1936,16 @@ subroutine exp_matrix(X,n,exp_X) enddo eigvalues_mat=0.d0 do i = 1,n - ! t = dsqrt(t^2) where t^2 are eigenvalues of X^2 - eigvalues(i) = dsqrt(eigvalues(i)) eigvalues_mat(i,i) = eigvalues(i)*eigvalues(i) enddo call dgemm('N','N',n,n,n,1.d0,eigvectors,size(eigvectors,1), & eigvalues_mat,size(eigvalues_mat,1),0.d0,matrix_tmp1,size(matrix_tmp1,1)) call dgemm('N','T',n,n,n,-1.d0,matrix_tmp1,size(matrix_tmp1,1), & eigvectors,size(eigvectors,1),0.d0,matrix_tmp2,size(matrix_tmp2,1)) + print*,'r2_mat = ' + do i = 1, n + write(*,'(100(F16.10,X))')r2_mat(:,i) + enddo print*,'r2_mat new = ' do i = 1, n write(*,'(100(F16.10,X))')matrix_tmp2(:,i) @@ -1964,7 +1970,8 @@ subroutine exp_matrix(X,n,exp_X) if(dabs(eigvalues(i)).gt.1.d-4)then eigvalues_mat(i,i) = dsin(eigvalues(i))/eigvalues(i) else ! Taylor development of sin(x)/x near x=0 = 1 - x^2/6 - eigvalues_mat(i,i) = 1.d0 - eigvalues(i)*eigvalues(i)*c_1_3*0.5d0 + eigvalues_mat(i,i) = 1.d0 - eigvalues(i)*eigvalues(i)*c_1_3*0.5d0 & + + eigvalues(i)*eigvalues(i)*eigvalues(i)*eigvalues(i)*c_1_3*0.025d0 endif enddo ! matrix_tmp1 = W t^-1 sin(t)