diff --git a/src/bi_ort_ints/semi_num_ints_mo.irp.f b/src/bi_ort_ints/semi_num_ints_mo.irp.f index 33f512cf..4762c25e 100644 --- a/src/bi_ort_ints/semi_num_ints_mo.irp.f +++ b/src/bi_ort_ints/semi_num_ints_mo.irp.f @@ -122,35 +122,40 @@ BEGIN_PROVIDER [ double precision, int2_grad1_u12_ao_transp, (ao_num, ao_num, 3, END_PROVIDER +! --- + BEGIN_PROVIDER [ double precision, int2_grad1_u12_bimo_transp, (mo_num, mo_num, 3, n_points_final_grid)] implicit none integer :: ipoint + double precision :: wall0, wall1 - print*,'providing int2_grad1_u12_bimo_transp' - double precision :: wall0, wall1 - call wall_time(wall0) - !$OMP PARALLEL & - !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (ipoint) & - !$OMP SHARED (n_points_final_grid,int2_grad1_u12_ao_transp,int2_grad1_u12_bimo_transp) - !$OMP DO SCHEDULE (dynamic) - do ipoint = 1, n_points_final_grid - call ao_to_mo_bi_ortho( int2_grad1_u12_ao_transp (1,1,1,ipoint), size(int2_grad1_u12_ao_transp , 1) & - , int2_grad1_u12_bimo_transp(1,1,1,ipoint), size(int2_grad1_u12_bimo_transp, 1) ) - call ao_to_mo_bi_ortho( int2_grad1_u12_ao_transp (1,1,2,ipoint), size(int2_grad1_u12_ao_transp , 1) & - , int2_grad1_u12_bimo_transp(1,1,2,ipoint), size(int2_grad1_u12_bimo_transp, 1) ) - call ao_to_mo_bi_ortho( int2_grad1_u12_ao_transp (1,1,3,ipoint), size(int2_grad1_u12_ao_transp , 1) & - , int2_grad1_u12_bimo_transp(1,1,3,ipoint), size(int2_grad1_u12_bimo_transp, 1) ) - enddo - !$OMP END DO - !$OMP END PARALLEL - call wall_time(wall1) - print*,'Wall time for providing int2_grad1_u12_bimo_transp',wall1 - wall0 + !print *, ' providing int2_grad1_u12_bimo_transp' + + call wall_time(wall0) + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (ipoint) & + !$OMP SHARED (n_points_final_grid,int2_grad1_u12_ao_transp,int2_grad1_u12_bimo_transp) + !$OMP DO SCHEDULE (dynamic) + do ipoint = 1, n_points_final_grid + call ao_to_mo_bi_ortho( int2_grad1_u12_ao_transp (1,1,1,ipoint), size(int2_grad1_u12_ao_transp , 1) & + , int2_grad1_u12_bimo_transp(1,1,1,ipoint), size(int2_grad1_u12_bimo_transp, 1) ) + call ao_to_mo_bi_ortho( int2_grad1_u12_ao_transp (1,1,2,ipoint), size(int2_grad1_u12_ao_transp , 1) & + , int2_grad1_u12_bimo_transp(1,1,2,ipoint), size(int2_grad1_u12_bimo_transp, 1) ) + call ao_to_mo_bi_ortho( int2_grad1_u12_ao_transp (1,1,3,ipoint), size(int2_grad1_u12_ao_transp , 1) & + , int2_grad1_u12_bimo_transp(1,1,3,ipoint), size(int2_grad1_u12_bimo_transp, 1) ) + enddo + !$OMP END DO + !$OMP END PARALLEL + + call wall_time(wall1) + !print *, ' Wall time for providing int2_grad1_u12_bimo_transp',wall1 - wall0 END_PROVIDER ! --- + BEGIN_PROVIDER [ double precision, int2_grad1_u12_bimo_t, (n_points_final_grid,3, mo_num, mo_num )] implicit none integer :: i, j, ipoint diff --git a/src/bi_ortho_mos/mos_rl.irp.f b/src/bi_ortho_mos/mos_rl.irp.f index 9e3ed358..d51999fc 100644 --- a/src/bi_ortho_mos/mos_rl.irp.f +++ b/src/bi_ortho_mos/mos_rl.irp.f @@ -61,7 +61,7 @@ subroutine mo_to_ao_bi_ortho(A_mo, LDA_mo, A_ao, LDA_ao) , 0.d0, tmp_1, size(tmp_1, 1) ) ! (ao_overlap x mo_r_coef) x A_mo - allocate( tmp_1(ao_num,mo_num) ) + allocate( tmp_2(ao_num,mo_num) ) call dgemm( 'N', 'N', ao_num, mo_num, mo_num, 1.d0 & , tmp_1, size(tmp_1, 1), A_mo, LDA_mo & , 0.d0, tmp_2, size(tmp_2, 1) ) @@ -73,7 +73,7 @@ subroutine mo_to_ao_bi_ortho(A_mo, LDA_mo, A_ao, LDA_ao) , 0.d0, tmp_1, size(tmp_1, 1) ) ! (ao_overlap x mo_r_coef) x A_mo x (ao_overlap x mo_l_coef).T - call dgemm( 'N', 'T', ao_num, mo_num, mo_num, 1.d0 & + call dgemm( 'N', 'T', ao_num, ao_num, mo_num, 1.d0 & , tmp_2, size(tmp_2, 1), tmp_1, size(tmp_1, 1) & , 0.d0, A_ao, LDA_ao ) diff --git a/src/non_hermit_dav/biorthog.irp.f b/src/non_hermit_dav/biorthog.irp.f index 926a20f1..89e5b4f4 100644 --- a/src/non_hermit_dav/biorthog.irp.f +++ b/src/non_hermit_dav/biorthog.irp.f @@ -283,16 +283,16 @@ subroutine non_hrmt_bieig(n, A, thr_d, thr_nd, leigvec, reigvec, n_real_eigv, ei ! ------------------------------------------------------------------------------------- ! - print *, ' ' - print *, ' Computing the left/right eigenvectors ...' - print *, ' ' + !print *, ' ' + !print *, ' Computing the left/right eigenvectors ...' + !print *, ' ' - allocate( WR(n), WI(n), VL(n,n), VR(n,n) ) + allocate(WR(n), WI(n), VL(n,n), VR(n,n)) - print *, ' fock matrix' - do i = 1, n - write(*, '(1000(F16.10,X))') A(i,:) - enddo + !print *, ' fock matrix' + !do i = 1, n + ! write(*, '(1000(F16.10,X))') A(i,:) + !enddo !thr_cut = 1.d-15 !call cancel_small_elmts(A, n, thr_cut) @@ -301,11 +301,11 @@ subroutine non_hrmt_bieig(n, A, thr_d, thr_nd, leigvec, reigvec, n_real_eigv, ei call lapack_diag_non_sym(n, A, WR, WI, VL, VR) !call lapack_diag_non_sym_new(n, A, WR, WI, VL, VR) - print *, ' ' - print *, ' eigenvalues' - do i = 1, n - write(*, '(1000(F16.10,X))') WR(i), WI(i) - enddo + !print *, ' ' + !print *, ' eigenvalues' + !do i = 1, n + ! write(*, '(1000(F16.10,X))') WR(i), WI(i) + !enddo !print *, ' right eigenvect bef' !do i = 1, n ! write(*, '(1000(F16.10,X))') VR(:,i) @@ -328,9 +328,10 @@ subroutine non_hrmt_bieig(n, A, thr_d, thr_nd, leigvec, reigvec, n_real_eigv, ei ! track & sort the real eigenvalues n_good = 0 - thr = 1.d-3 + !thr = 100d0 + thr = Im_thresh_tcscf do i = 1, n - print*, 'Re(i) + Im(i)', WR(i), WI(i) + !print*, 'Re(i) + Im(i)', WR(i), WI(i) if(dabs(WI(i)) .lt. thr) then n_good += 1 else @@ -402,23 +403,24 @@ subroutine non_hrmt_bieig(n, A, thr_d, thr_nd, leigvec, reigvec, n_real_eigv, ei if( (accu_nd .lt. thr_nd) .and. (dabs(accu_d-dble(n_real_eigv))/dble(n_real_eigv) .lt. thr_d) ) then - print *, ' lapack vectors are normalized and bi-orthogonalized' + !print *, ' lapack vectors are normalized and bi-orthogonalized' deallocate(S) return - elseif( (accu_nd .lt. thr_nd) .and. (dabs(accu_d-dble(n_real_eigv))/dble(n_real_eigv) .gt. thr_d) ) then + ! accu_nd is modified after adding the normalization + !elseif( (accu_nd .lt. thr_nd) .and. (dabs(accu_d-dble(n_real_eigv))/dble(n_real_eigv) .gt. thr_d) ) then - print *, ' lapack vectors are not normalized but bi-orthogonalized' - call check_biorthog_binormalize(n, n_real_eigv, leigvec, reigvec, thr_d, thr_nd, .true.) + ! print *, ' lapack vectors are not normalized but bi-orthogonalized' + ! call check_biorthog_binormalize(n, n_real_eigv, leigvec, reigvec, thr_d, thr_nd, .true.) - call check_EIGVEC(n, n, A, eigval, leigvec, reigvec, thr_diag, thr_norm, .true.) + ! call check_EIGVEC(n, n, A, eigval, leigvec, reigvec, thr_diag, thr_norm, .true.) - deallocate(S) - return + ! deallocate(S) + ! return else - print *, ' lapack vectors are not normalized neither bi-orthogonalized' + !print *, ' lapack vectors are not normalized neither bi-orthogonalized' ! --- diff --git a/src/non_hermit_dav/lapack_diag_non_hermit.irp.f b/src/non_hermit_dav/lapack_diag_non_hermit.irp.f index 53c62ce8..0d652af4 100644 --- a/src/non_hermit_dav/lapack_diag_non_hermit.irp.f +++ b/src/non_hermit_dav/lapack_diag_non_hermit.irp.f @@ -930,7 +930,7 @@ subroutine check_EIGVEC(n, m, A, eigval, leigvec, reigvec, thr_diag, thr_norm, s tmp_abs = tmp_abs + tmp V_nrm = V_nrm + U_nrm - write(*,'(I4,X,(100(F25.16,X)))')j,eigval(j), tmp, U_nrm + !write(*,'(I4,X,(100(F25.16,X)))') j,eigval(j), tmp, U_nrm enddo @@ -973,7 +973,7 @@ subroutine check_EIGVEC(n, m, A, eigval, leigvec, reigvec, thr_diag, thr_norm, s tmp_abs = tmp_abs + tmp V_nrm = V_nrm + U_nrm - write(*,'(I4,X,(100(F25.16,X)))')j,eigval(j), tmp, U_nrm + !write(*,'(I4,X,(100(F25.16,X)))') j,eigval(j), tmp, U_nrm enddo @@ -1082,7 +1082,7 @@ subroutine impose_weighted_orthog_svd(n, m, W, C) double precision, allocatable :: S(:,:), tmp(:,:) double precision, allocatable :: U(:,:), Vt(:,:), D(:) - print *, ' apply SVD to orthogonalize & normalize weighted vectors' + !print *, ' apply SVD to orthogonalize & normalize weighted vectors' ! --- @@ -1097,10 +1097,10 @@ subroutine impose_weighted_orthog_svd(n, m, W, C) , 0.d0, S, size(S, 1) ) deallocate(tmp) - print *, ' overlap bef SVD: ' - do i = 1, m - write(*, '(1000(F16.10,X))') S(i,:) - enddo + !print *, ' overlap bef SVD: ' + !do i = 1, m + ! write(*, '(1000(F16.10,X))') S(i,:) + !enddo ! --- @@ -1160,10 +1160,10 @@ subroutine impose_weighted_orthog_svd(n, m, W, C) , 0.d0, S, size(S, 1) ) deallocate(tmp) - print *, ' overlap aft SVD: ' - do i = 1, m - write(*, '(1000(F16.10,X))') S(i,:) - enddo + !print *, ' overlap aft SVD: ' + !do i = 1, m + ! write(*, '(1000(F16.10,X))') S(i,:) + !enddo deallocate(S) @@ -1185,7 +1185,7 @@ subroutine impose_orthog_svd(n, m, C) double precision, allocatable :: S(:,:), tmp(:,:) double precision, allocatable :: U(:,:), Vt(:,:), D(:) - print *, ' apply SVD to orthogonalize & normalize vectors' + !print *, ' apply SVD to orthogonalize & normalize vectors' ! --- @@ -1196,10 +1196,10 @@ subroutine impose_orthog_svd(n, m, C) , C, size(C, 1), C, size(C, 1) & , 0.d0, S, size(S, 1) ) - print *, ' eigenvec overlap bef SVD: ' - do i = 1, m - write(*, '(1000(F16.10,X))') S(i,:) - enddo + !print *, ' eigenvec overlap bef SVD: ' + !do i = 1, m + ! write(*, '(1000(F16.10,X))') S(i,:) + !enddo ! --- @@ -1224,6 +1224,7 @@ subroutine impose_orthog_svd(n, m, C) if(num_linear_dependencies > 0) then write(*,*) ' linear dependencies = ', num_linear_dependencies write(*,*) ' m = ', m + write(*,*) ' try with Graham-Schmidt' stop endif @@ -1256,10 +1257,10 @@ subroutine impose_orthog_svd(n, m, C) , C, size(C, 1), C, size(C, 1) & , 0.d0, S, size(S, 1) ) - print *, ' eigenvec overlap aft SVD: ' - do i = 1, m - write(*, '(1000(F16.10,X))') S(i,:) - enddo + !print *, ' eigenvec overlap aft SVD: ' + !do i = 1, m + ! write(*, '(1000(F16.10,X))') S(i,:) + !enddo deallocate(S) @@ -1296,10 +1297,10 @@ subroutine impose_orthog_svd_overlap(n, m, C, overlap) , 0.d0, S, size(S, 1) ) deallocate(Stmp) - print *, ' eigenvec overlap bef SVD: ' - do i = 1, m - write(*, '(1000(F16.10,X))') S(i,:) - enddo + !print *, ' eigenvec overlap bef SVD: ' + !do i = 1, m + ! write(*, '(1000(F16.10,X))') S(i,:) + !enddo ! --- @@ -1358,10 +1359,10 @@ subroutine impose_orthog_svd_overlap(n, m, C, overlap) , 0.d0, S, size(S, 1) ) deallocate(Stmp) - print *, ' eigenvec overlap aft SVD: ' - do i = 1, m - write(*, '(1000(F16.10,X))') S(i,:) - enddo + !print *, ' eigenvec overlap aft SVD: ' + !do i = 1, m + ! write(*, '(1000(F16.10,X))') S(i,:) + !enddo deallocate(S) end subroutine impose_orthog_svd_overlap @@ -1528,11 +1529,11 @@ subroutine impose_orthog_degen_eigvec(n, e0, C0) enddo - do i = 1, n - if(deg_num(i).gt.1) then - print *, ' degen on', i, deg_num(i) - endif - enddo + !do i = 1, n + ! if(deg_num(i) .gt. 1) then + ! print *, ' degen on', i, deg_num(i) + ! endif + !enddo ! --- @@ -1677,7 +1678,7 @@ subroutine check_biorthog_binormalize(n, m, Vl, Vr, thr_d, thr_nd, stop_ifnot) double precision :: accu_d, accu_nd, s_tmp double precision, allocatable :: S(:,:) - print *, ' check bi-orthonormality' + !print *, ' check bi-orthonormality' ! --- @@ -1714,15 +1715,19 @@ subroutine check_biorthog_binormalize(n, m, Vl, Vr, thr_d, thr_nd, stop_ifnot) enddo enddo accu_nd = dsqrt(accu_nd) / dble(m) - print*, ' diag acc: ', accu_d - print*, ' nondiag acc: ', accu_nd + !print*, ' diag acc bef = ', accu_d + !print*, ' nondiag acc bef = ', accu_nd ! --- if( (accu_nd .lt. thr_nd) .and. (dabs(accu_d-dble(m))/dble(m) .gt. thr_d) ) then do i = 1, m - print *, i, S(i,i) + if(S(i,i) <= 0.d0) then + print *, ' overap negative' + print *, i, S(i,i) + exit + endif if(dabs(S(i,i) - 1.d0) .gt. thr_d) then s_tmp = 1.d0 / dsqrt(S(i,i)) do j = 1, n @@ -1757,8 +1762,8 @@ subroutine check_biorthog_binormalize(n, m, Vl, Vr, thr_d, thr_nd, stop_ifnot) enddo enddo accu_nd = dsqrt(accu_nd) / dble(m) - print *, ' diag acc: ', accu_d - print *, ' nondiag acc: ', accu_nd + !print *, ' diag acc aft = ', accu_d + !print *, ' nondiag acc aft = ', accu_nd deallocate(S) @@ -1801,10 +1806,10 @@ subroutine check_weighted_biorthog(n, m, W, Vl, Vr, thr_d, thr_nd, accu_d, accu_ , 0.d0, S, size(S, 1) ) deallocate(tmp) - print *, ' overlap matrix:' - do i = 1, m - write(*,'(1000(F16.10,X))') S(i,:) - enddo + !print *, ' overlap matrix:' + !do i = 1, m + ! write(*,'(1000(F16.10,X))') S(i,:) + !enddo accu_d = 0.d0 accu_nd = 0.d0 @@ -1852,17 +1857,18 @@ subroutine check_biorthog(n, m, Vl, Vr, accu_d, accu_nd, S, thr_d, thr_nd, stop_ integer :: i, j double precision, allocatable :: SS(:,:) - print *, ' check bi-orthogonality' + !print *, ' check bi-orthogonality' ! --- call dgemm( 'T', 'N', m, m, n, 1.d0 & , Vl, size(Vl, 1), Vr, size(Vr, 1) & , 0.d0, S, size(S, 1) ) - print *, ' overlap matrix:' - do i = 1, m - write(*,'(1000(F16.10,X))') S(i,:) - enddo + + !print *, ' overlap matrix:' + !do i = 1, m + ! write(*,'(1000(F16.10,X))') S(i,:) + !enddo accu_d = 0.d0 accu_nd = 0.d0 @@ -1877,12 +1883,12 @@ subroutine check_biorthog(n, m, Vl, Vr, accu_d, accu_nd, S, thr_d, thr_nd, stop_ enddo accu_nd = dsqrt(accu_nd) / dble(m) - print *, ' accu_nd = ', accu_nd - print *, ' accu_d = ', dabs(accu_d-dble(m))/dble(m) + !print *, ' accu_nd = ', accu_nd + !print *, ' accu_d = ', dabs(accu_d-dble(m))/dble(m) ! --- - if( stop_ifnot .and. ((accu_nd .gt. thr_nd) .or. dabs(accu_d-dble(m))/dble(m) .gt. thr_d) ) then + if(stop_ifnot .and. ((accu_nd .gt. thr_nd) .or. dabs(accu_d-dble(m))/dble(m) .gt. thr_d)) then print *, ' non bi-orthogonal vectors !' print *, ' accu_nd = ', accu_nd print *, ' accu_d = ', dabs(accu_d-dble(m))/dble(m) @@ -1912,12 +1918,12 @@ subroutine check_orthog(n, m, V, accu_d, accu_nd, S) , V, size(V, 1), V, size(V, 1) & , 0.d0, S, size(S, 1) ) - print *, '' - print *, ' overlap matrix:' - do i = 1, m - write(*,'(1000(F16.10,X))') S(i,:) - enddo - print *, '' + !print *, '' + !print *, ' overlap matrix:' + !do i = 1, m + ! write(*,'(1000(F16.10,X))') S(i,:) + !enddo + !print *, '' accu_d = 0.d0 accu_nd = 0.d0 @@ -1981,11 +1987,11 @@ subroutine impose_biorthog_degen_eigvec(n, e0, L0, R0) enddo enddo - do i = 1, n - if(deg_num(i).gt.1) then - print *, ' degen on', i, deg_num(i), e0(i) - endif - enddo + !do i = 1, n + ! if(deg_num(i) .gt. 1) then + ! print *, ' degen on', i, deg_num(i), e0(i) + ! endif + !enddo ! --- @@ -2181,11 +2187,11 @@ subroutine impose_unique_biorthog_degen_eigvec(n, thr_d, thr_nd, e0, C0, W0, L0, enddo enddo - do i = 1, n - if(deg_num(i).gt.1) then - print *, ' degen on', i, deg_num(i) - endif - enddo + !do i = 1, n + ! if(deg_num(i) .gt. 1) then + ! print *, ' degen on', i, deg_num(i) + ! endif + !enddo ! --- @@ -2414,10 +2420,10 @@ subroutine impose_biorthog_svd(n, m, L, R) , L, size(L, 1), R, size(R, 1) & , 0.d0, S, size(S, 1) ) - print *, ' overlap bef SVD: ' - do i = 1, m - write(*, '(1000(F16.10,X))') S(i,:) - enddo + !print *, ' overlap bef SVD: ' + !do i = 1, m + ! write(*, '(1000(F16.10,X))') S(i,:) + !enddo ! --- @@ -2489,10 +2495,11 @@ subroutine impose_biorthog_svd(n, m, L, R) , L, size(L, 1), R, size(R, 1) & , 0.d0, S, size(S, 1) ) - print *, ' overlap aft SVD: ' - do i = 1, m - write(*, '(1000(F16.10,X))') S(i,:) - enddo + !print *, ' overlap aft SVD: ' + !do i = 1, m + ! write(*, '(1000(F16.10,X))') S(i,:) + !enddo + deallocate(S) ! --- @@ -2806,10 +2813,10 @@ subroutine impose_weighted_biorthog_svd(n, m, overlap, L, R) , 0.d0, S, size(S, 1) ) deallocate(Stmp) - print *, ' overlap bef SVD: ' - do i = 1, m - write(*, '(1000(F25.16,X))') S(i,:) - enddo + !print *, ' overlap bef SVD: ' + !do i = 1, m + ! write(*, '(1000(F25.16,X))') S(i,:) + !enddo ! --- @@ -2886,10 +2893,11 @@ subroutine impose_weighted_biorthog_svd(n, m, overlap, L, R) , 0.d0, S, size(S, 1) ) deallocate(Stmp) - print *, ' overlap aft SVD with overlap: ' - do i = 1, m - write(*, '(1000(F16.10,X))') S(i,:) - enddo + !print *, ' overlap aft SVD with overlap: ' + !do i = 1, m + ! write(*, '(1000(F16.10,X))') S(i,:) + !enddo + deallocate(S) return diff --git a/src/scf_utils/roothaan_hall_scf.irp.f b/src/scf_utils/roothaan_hall_scf.irp.f index 3b9eaeb4..56a1ed8e 100644 --- a/src/scf_utils/roothaan_hall_scf.irp.f +++ b/src/scf_utils/roothaan_hall_scf.irp.f @@ -29,11 +29,11 @@ END_DOC call write_time(6) - print*,'Energy of the guess = ',SCF_energy + print*,'energy of the guess = ',SCF_energy write(6,'(A4, 1X, A16, 1X, A16, 1X, A16, 1X, A16)') & '====','================','================','================','================' write(6,'(A4, 1X, A16, 1X, A16, 1X, A16, 1X, A16)') & - ' N ', 'Energy ', 'Energy diff ', 'DIIS error ', 'Level shift ' + ' N ', 'energy ', 'energy diff ', 'DIIS error ', 'Level shift ' write(6,'(A4, 1X, A16, 1X, A16, 1X, A16, 1X, A16)') & '====','================','================','================','================' @@ -69,9 +69,9 @@ END_DOC if ( (scf_algorithm == 'DIIS').and.(dabs(Delta_energy_SCF) > 1.d-6) ) then ! Store Fock and error matrices at each iteration + index_dim_DIIS = mod(dim_DIIS-1,max_dim_DIIS)+1 do j=1,ao_num do i=1,ao_num - index_dim_DIIS = mod(dim_DIIS-1,max_dim_DIIS)+1 Fock_matrix_DIIS (i,j,index_dim_DIIS) = Fock_matrix_AO(i,j) error_matrix_DIIS(i,j,index_dim_DIIS) = FPS_SPF_matrix_AO(i,j) enddo @@ -106,8 +106,8 @@ END_DOC ! SCF energy energy_SCF = SCF_energy - Delta_Energy_SCF = energy_SCF - energy_SCF_previous - if ( (SCF_algorithm == 'DIIS').and.(Delta_Energy_SCF > 0.d0) ) then + Delta_energy_SCF = energy_SCF - energy_SCF_previous + if ( (SCF_algorithm == 'DIIS').and.(Delta_energy_SCF > 0.d0) ) then Fock_matrix_AO(1:ao_num,1:ao_num) = Fock_matrix_DIIS (1:ao_num,1:ao_num,index_dim_DIIS) Fock_matrix_AO_alpha = Fock_matrix_AO*0.5d0 Fock_matrix_AO_beta = Fock_matrix_AO*0.5d0 @@ -131,15 +131,17 @@ END_DOC call initialize_mo_coef_begin_iteration endif TOUCH mo_coef - Delta_Energy_SCF = SCF_energy - energy_SCF_previous + Delta_energy_SCF = SCF_energy - energy_SCF_previous energy_SCF = SCF_energy if (level_shift-level_shift_save > 40.d0) then level_shift = level_shift_save * 4.d0 SOFT_TOUCH level_shift exit endif + dim_DIIS=0 enddo + level_shift = level_shift * 0.5d0 SOFT_TOUCH level_shift energy_SCF_previous = energy_SCF @@ -175,7 +177,7 @@ END_DOC call save_mos endif - call write_double(6, Energy_SCF, 'SCF energy') + call write_double(6, energy_SCF, 'SCF energy') call write_time(6) diff --git a/src/tc_keywords/EZFIO.cfg b/src/tc_keywords/EZFIO.cfg index 8db73f9a..9df73154 100644 --- a/src/tc_keywords/EZFIO.cfg +++ b/src/tc_keywords/EZFIO.cfg @@ -86,13 +86,13 @@ default: False type: Threshold doc: Threshold on the convergence of the Hartree Fock energy. interface: ezfio,provider,ocaml -default: 1.e-10 +default: 1.e-12 [n_it_tcscf_max] type: Strictly_positive_int doc: Maximum number of SCF iterations interface: ezfio,provider,ocaml -default: 500 +default: 100 [selection_tc] type: integer @@ -160,3 +160,9 @@ doc: Type of TCSCF algorithm used. Possible choices are [Simple | DIIS] interface: ezfio,provider,ocaml default: DIIS +[im_thresh_tcscf] +type: Threshold +doc: Thresholds on the Imag part of energy +interface: ezfio,provider,ocaml +default: 1.e-7 + diff --git a/src/tc_scf/diago_bi_ort_tcfock.irp.f b/src/tc_scf/diago_bi_ort_tcfock.irp.f index 29ca0efe..9c571f8a 100644 --- a/src/tc_scf/diago_bi_ort_tcfock.irp.f +++ b/src/tc_scf/diago_bi_ort_tcfock.irp.f @@ -76,6 +76,8 @@ , fock_tc_reigvec_mo, size(fock_tc_reigvec_mo, 1) & , 0.d0, overlap_fock_tc_eigvec_mo, size(overlap_fock_tc_eigvec_mo, 1) ) + ! --- + accu_d = 0.d0 accu_nd = 0.d0 do i = 1, mo_num @@ -92,22 +94,24 @@ endif enddo enddo - accu_nd = dsqrt(accu_nd)/accu_d - + accu_nd = dsqrt(accu_nd) / accu_d if(accu_nd .gt. thr_nd) then print *, ' bi-orthog failed' - print*,'accu_nd MO = ', accu_nd, thr_nd - print*,'overlap_fock_tc_eigvec_mo = ' + print *, ' accu_nd MO = ', accu_nd, thr_nd + print *, ' overlap_fock_tc_eigvec_mo = ' do i = 1, mo_num write(*,'(100(F16.10,X))') overlap_fock_tc_eigvec_mo(i,:) enddo - stop + stop endif - if( dabs(accu_d - dble(mo_num))/dble(mo_num) .gt. thr_d ) then - print *, 'mo_num = ', mo_num - print *, 'accu_d MO = ', accu_d, thr_d - print *, 'normalizing vectors ...' + ! --- + + if(dabs(accu_d - dble(mo_num))/dble(mo_num) .gt. thr_d) then + + print *, ' mo_num = ', mo_num + print *, ' accu_d MO = ', accu_d, thr_d + print *, ' normalizing vectors ...' do i = 1, mo_num norm = dsqrt(dabs(overlap_fock_tc_eigvec_mo(i,i))) if(norm .gt. thr_d) then @@ -117,12 +121,43 @@ enddo endif enddo + call dgemm( "T", "N", mo_num, mo_num, mo_num, 1.d0 & , fock_tc_leigvec_mo, size(fock_tc_leigvec_mo, 1) & , fock_tc_reigvec_mo, size(fock_tc_reigvec_mo, 1) & , 0.d0, overlap_fock_tc_eigvec_mo, size(overlap_fock_tc_eigvec_mo, 1) ) + + accu_d = 0.d0 + accu_nd = 0.d0 + do i = 1, mo_num + do k = 1, mo_num + if(i==k) then + accu_tmp = overlap_fock_tc_eigvec_mo(k,i) + accu_d += dabs(accu_tmp) + else + accu_tmp = overlap_fock_tc_eigvec_mo(k,i) + accu_nd += accu_tmp * accu_tmp + if(dabs(overlap_fock_tc_eigvec_mo(k,i)) .gt. thr_nd)then + print *, 'k,i', k, i, overlap_fock_tc_eigvec_mo(k,i) + endif + endif + enddo + enddo + accu_nd = dsqrt(accu_nd) / accu_d + if(accu_nd .gt. thr_nd) then + print *, ' bi-orthog failed' + print *, ' accu_nd MO = ', accu_nd, thr_nd + print *, ' overlap_fock_tc_eigvec_mo = ' + do i = 1, mo_num + write(*,'(100(F16.10,X))') overlap_fock_tc_eigvec_mo(i,:) + enddo + stop + endif + endif + ! --- + END_PROVIDER ! --- diff --git a/src/tc_scf/diis_tcscf.irp.f b/src/tc_scf/diis_tcscf.irp.f index cf339175..ff1077f5 100644 --- a/src/tc_scf/diis_tcscf.irp.f +++ b/src/tc_scf/diis_tcscf.irp.f @@ -27,6 +27,7 @@ BEGIN_PROVIDER [double precision, Q_alpha, (ao_num, ao_num) ] implicit none + Q_alpha = 0.d0 call dgemm( 'N', 'T', ao_num, ao_num, elec_alpha_num, 1.d0 & , mo_r_coef, size(mo_r_coef, 1), mo_l_coef, size(mo_l_coef, 1) & , 0.d0, Q_alpha, size(Q_alpha, 1) ) @@ -47,6 +48,7 @@ BEGIN_PROVIDER [ double precision, Q_beta, (ao_num, ao_num) ] implicit none + Q_beta = 0.d0 call dgemm( 'N', 'T', ao_num, ao_num, elec_beta_num, 1.d0 & , mo_r_coef, size(mo_r_coef, 1), mo_l_coef, size(mo_l_coef, 1) & , 0.d0, Q_beta, size(Q_beta, 1) ) @@ -113,15 +115,18 @@ BEGIN_PROVIDER [double precision, FQS_SQF_ao, (ao_num, ao_num)] , 0.d0, FQS_SQF_ao, size(FQS_SQF_ao, 1) ) ! S x Q + tmp = 0.d0 call dgemm( 'N', 'N', ao_num, ao_num, ao_num, 1.d0 & , ao_overlap, size(ao_overlap, 1), Q_matrix, size(Q_matrix, 1) & , 0.d0, tmp, size(tmp, 1) ) - ! F x P x S - S x P x F + ! F x Q x S - S x Q x F call dgemm( 'N', 'N', ao_num, ao_num, ao_num, -1.d0 & , tmp, size(tmp, 1), Fock_matrix_tc_ao_tot, size(Fock_matrix_tc_ao_tot, 1) & , 1.d0, FQS_SQF_ao, size(FQS_SQF_ao, 1) ) + deallocate(tmp) + END_PROVIDER ! --- diff --git a/src/tc_scf/fock_tc.irp.f b/src/tc_scf/fock_tc.irp.f index b9611836..c3642a7e 100644 --- a/src/tc_scf/fock_tc.irp.f +++ b/src/tc_scf/fock_tc.irp.f @@ -74,68 +74,61 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_tc_ao_beta, (ao_num, ao_num)] + two_e_tc_non_hermit_integral_beta END_PROVIDER -! --- - -!BEGIN_PROVIDER [ double precision, Fock_matrix_tc_ao_tot, (ao_num, ao_num) ] -! implicit none -! BEGIN_DOC -! ! Total alpha+beta TC Fock matrix : h_c + Two-e^TC terms on the AO basis -! END_DOC -! Fock_matrix_tc_ao_tot = 0.5d0 * (Fock_matrix_tc_ao_alpha + Fock_matrix_tc_ao_beta) -!END_PROVIDER ! --- BEGIN_PROVIDER [ double precision, Fock_matrix_tc_mo_alpha, (mo_num, mo_num) ] - implicit none + BEGIN_DOC - ! Total alpha TC Fock matrix : h_c + Two-e^TC terms on the MO basis + ! Total alpha TC Fock matrix : h_c + Two-e^TC terms on the MO basis END_DOC - if(bi_ortho)then - call ao_to_mo_bi_ortho( Fock_matrix_tc_ao_alpha, size(Fock_matrix_tc_ao_alpha, 1) & - , Fock_matrix_tc_mo_alpha, size(Fock_matrix_tc_mo_alpha, 1) ) - if(three_body_h_tc)then - Fock_matrix_tc_mo_alpha += fock_a_tot_3e_bi_orth - endif + + implicit none + + if(bi_ortho) then + + call ao_to_mo_bi_ortho( Fock_matrix_tc_ao_alpha, size(Fock_matrix_tc_ao_alpha, 1) & + , Fock_matrix_tc_mo_alpha, size(Fock_matrix_tc_mo_alpha, 1) ) + if(three_body_h_tc) then + Fock_matrix_tc_mo_alpha += fock_a_tot_3e_bi_orth + endif + else - call ao_to_mo( Fock_matrix_tc_ao_alpha, size(Fock_matrix_tc_ao_alpha, 1) & + call ao_to_mo( Fock_matrix_tc_ao_alpha, size(Fock_matrix_tc_ao_alpha, 1) & , Fock_matrix_tc_mo_alpha, size(Fock_matrix_tc_mo_alpha, 1) ) + endif + END_PROVIDER ! --- BEGIN_PROVIDER [ double precision, Fock_matrix_tc_mo_beta, (mo_num,mo_num) ] - implicit none + BEGIN_DOC - ! Total beta TC Fock matrix : h_c + Two-e^TC terms on the MO basis + ! Total beta TC Fock matrix : h_c + Two-e^TC terms on the MO basis END_DOC - if(bi_ortho)then + + implicit none + + if(bi_ortho) then + call ao_to_mo_bi_ortho( Fock_matrix_tc_ao_beta, size(Fock_matrix_tc_ao_beta, 1) & , Fock_matrix_tc_mo_beta, size(Fock_matrix_tc_mo_beta, 1) ) - if(three_body_h_tc)then - Fock_matrix_tc_mo_beta += fock_b_tot_3e_bi_orth - endif + + if(three_body_h_tc) then + Fock_matrix_tc_mo_beta += fock_b_tot_3e_bi_orth + endif + else - call ao_to_mo( Fock_matrix_tc_ao_beta, size(Fock_matrix_tc_ao_beta, 1) & + + call ao_to_mo( Fock_matrix_tc_ao_beta, size(Fock_matrix_tc_ao_beta, 1) & , Fock_matrix_tc_mo_beta, size(Fock_matrix_tc_mo_beta, 1) ) + endif + END_PROVIDER -! --- - -!BEGIN_PROVIDER [ double precision, Fock_matrix_tc_mo_tot, (mo_num, mo_num)] -! implicit none -! BEGIN_DOC -! ! Total alpha+beta TC Fock matrix : h_c + Two-e^TC terms on the MO basis -! END_DOC -! Fock_matrix_tc_mo_tot = 0.5d0 * (Fock_matrix_tc_mo_alpha + Fock_matrix_tc_mo_beta) -! if(three_body_h_tc) then -! Fock_matrix_tc_mo_tot += fock_3_mat -! endif -! !call restore_symmetry(mo_num, mo_num, Fock_matrix_tc_mo_tot, mo_num, 1.d-10) -!END_PROVIDER - ! --- BEGIN_PROVIDER [ double precision, grad_non_hermit_left] diff --git a/src/tc_scf/fock_three.irp.f b/src/tc_scf/fock_three.irp.f index f73a5049..35b6aac6 100644 --- a/src/tc_scf/fock_three.irp.f +++ b/src/tc_scf/fock_three.irp.f @@ -70,52 +70,72 @@ subroutine give_fock_ia_three_e_total(i,a,contrib) end +! --- + BEGIN_PROVIDER [double precision, diag_three_elem_hf] - implicit none - integer :: i,j,k,ipoint,mm - double precision :: contrib,weight,four_third,one_third,two_third,exchange_int_231 - print*,'providing diag_three_elem_hf' - if(.not.three_body_h_tc)then - diag_three_elem_hf = 0.d0 - else - if(.not.bi_ortho)then - one_third = 1.d0/3.d0 - two_third = 2.d0/3.d0 - four_third = 4.d0/3.d0 - diag_three_elem_hf = 0.d0 - do i = 1, elec_beta_num - do j = 1, elec_beta_num - do k = 1, elec_beta_num - call give_integrals_3_body(k,j,i,j,i,k,exchange_int_231) - diag_three_elem_hf += two_third * exchange_int_231 - enddo - enddo - enddo - do mm = 1, 3 - do ipoint = 1, n_points_final_grid - weight = final_weight_at_r_vector(ipoint) - contrib = 3.d0 * fock_3_w_kk_sum(ipoint,mm) * fock_3_rho_beta(ipoint) * fock_3_w_kk_sum(ipoint,mm) & - -2.d0 * fock_3_w_kl_mo_k_mo_l(ipoint,mm) * fock_3_w_kk_sum(ipoint,mm) & - -1.d0 * fock_3_rho_beta(ipoint) * fock_3_w_kl_w_kl(ipoint,mm) - contrib *= four_third - contrib += -two_third * fock_3_rho_beta(ipoint) * fock_3_w_kl_w_kl(ipoint,mm) & - - four_third * fock_3_w_kk_sum(ipoint,mm) * fock_3_w_kl_mo_k_mo_l(ipoint,mm) - diag_three_elem_hf += weight * contrib - enddo - enddo - diag_three_elem_hf = - diag_three_elem_hf + + implicit none + integer :: i, j, k, ipoint, mm + double precision :: contrib, weight, four_third, one_third, two_third, exchange_int_231 + double precision :: integral_aaa, hthree, integral_aab, integral_abb, integral_bbb + + !print *, ' providing diag_three_elem_hf' + + if(.not. three_body_h_tc) then + + diag_three_elem_hf = 0.d0 + else - double precision :: integral_aaa,hthree, integral_aab,integral_abb,integral_bbb - provide mo_l_coef mo_r_coef - call give_aaa_contrib(integral_aaa) - call give_aab_contrib(integral_aab) - call give_abb_contrib(integral_abb) - call give_bbb_contrib(integral_bbb) - diag_three_elem_hf = integral_aaa + integral_aab + integral_abb + integral_bbb + + if(.not. bi_ortho) then + + ! --- + + one_third = 1.d0/3.d0 + two_third = 2.d0/3.d0 + four_third = 4.d0/3.d0 + diag_three_elem_hf = 0.d0 + do i = 1, elec_beta_num + do j = 1, elec_beta_num + do k = 1, elec_beta_num + call give_integrals_3_body(k, j, i, j, i, k,exchange_int_231) + diag_three_elem_hf += two_third * exchange_int_231 + enddo + enddo + enddo + do mm = 1, 3 + do ipoint = 1, n_points_final_grid + weight = final_weight_at_r_vector(ipoint) + contrib = 3.d0 * fock_3_w_kk_sum(ipoint,mm) * fock_3_rho_beta(ipoint) * fock_3_w_kk_sum(ipoint,mm) & + - 2.d0 * fock_3_w_kl_mo_k_mo_l(ipoint,mm) * fock_3_w_kk_sum(ipoint,mm) & + - 1.d0 * fock_3_rho_beta(ipoint) * fock_3_w_kl_w_kl(ipoint,mm) + contrib *= four_third + contrib += -two_third * fock_3_rho_beta(ipoint) * fock_3_w_kl_w_kl(ipoint,mm) & + -four_third * fock_3_w_kk_sum(ipoint,mm) * fock_3_w_kl_mo_k_mo_l(ipoint,mm) + diag_three_elem_hf += weight * contrib + enddo + enddo + + diag_three_elem_hf = - diag_three_elem_hf + + ! --- + + else + + provide mo_l_coef mo_r_coef + call give_aaa_contrib(integral_aaa) + call give_aab_contrib(integral_aab) + call give_abb_contrib(integral_abb) + call give_bbb_contrib(integral_bbb) + diag_three_elem_hf = integral_aaa + integral_aab + integral_abb + integral_bbb + + endif + endif - endif + END_PROVIDER +! --- BEGIN_PROVIDER [ double precision, fock_3_mat_a_op_sh, (mo_num, mo_num)] implicit none diff --git a/src/tc_scf/rh_tcscf.irp.f b/src/tc_scf/rh_tcscf.irp.f index dc7a34fc..597c3e67 100644 --- a/src/tc_scf/rh_tcscf.irp.f +++ b/src/tc_scf/rh_tcscf.irp.f @@ -16,7 +16,8 @@ subroutine rh_tcscf() double precision :: energy_TCSCF_previous, delta_energy_TCSCF double precision :: gradie_TCSCF_previous, delta_gradie_TCSCF double precision :: max_error_DIIS_TCSCF - double precision :: level_shift_TCSCF_save + double precision :: level_shift_save + double precision :: delta_energy_tmp, delta_gradie_tmp double precision, allocatable :: F_DIIS(:,:,:), e_DIIS(:,:,:) double precision, allocatable :: mo_r_coef_save(:,:), mo_l_coef_save(:,:) @@ -60,8 +61,8 @@ subroutine rh_tcscf() PROVIDE FQS_SQF_ao Fock_matrix_tc_ao_tot do while( (max_error_DIIS_TCSCF > threshold_DIIS_nonzero_TCSCF) .or. & - (dabs(delta_energy_TCSCF) > thresh_TCSCF) .or. & - (dabs(delta_gradie_TCSCF) > dsqrt(thresh_TCSCF)) ) + !(dabs(delta_energy_TCSCF) > thresh_TCSCF) .or. & + (dabs(gradie_TCSCF_previous) > dsqrt(thresh_TCSCF)) ) iteration_TCSCF += 1 if(iteration_TCSCF > n_it_TCSCF_max) then @@ -69,11 +70,6 @@ subroutine rh_tcscf() exit endif - ! TODO - !if(frozen_orb_scf) then - ! call initialize_mo_coef_begin_iteration - !endif - ! current size of the DIIS space dim_DIIS = min(dim_DIIS+1, max_dim_DIIS_TCSCF) @@ -91,13 +87,19 @@ subroutine rh_tcscf() enddo ! Compute the extrapolated Fock matrix - call extrapolate_TC_Fock_matrix( e_DIIS, F_DIIS & + call extrapolate_TC_Fock_matrix( e_DIIS, F_DIIS & , Fock_matrix_tc_ao_tot, size(Fock_matrix_tc_ao_tot, 1) & , iteration_TCSCF, dim_DIIS ) Fock_matrix_tc_ao_alpha = 0.5d0 * Fock_matrix_tc_ao_tot Fock_matrix_tc_ao_beta = 0.5d0 * Fock_matrix_tc_ao_tot - TOUCH Fock_matrix_tc_ao_alpha Fock_matrix_tc_ao_beta + !TOUCH Fock_matrix_tc_ao_alpha Fock_matrix_tc_ao_beta + + call ao_to_mo_bi_ortho( Fock_matrix_tc_ao_alpha, size(Fock_matrix_tc_ao_alpha, 1) & + , Fock_matrix_tc_mo_alpha, size(Fock_matrix_tc_mo_alpha, 1) ) + call ao_to_mo_bi_ortho( Fock_matrix_tc_ao_beta , size(Fock_matrix_tc_ao_beta , 1) & + , Fock_matrix_tc_mo_beta , size(Fock_matrix_tc_mo_beta , 1) ) + TOUCH Fock_matrix_tc_mo_alpha Fock_matrix_tc_mo_beta endif @@ -109,15 +111,54 @@ subroutine rh_tcscf() ! --- - ! TODO - !if(frozen_orb_scf) then - ! call reorder_core_orb - ! call initialize_mo_coef_begin_iteration - !endif - ! calculate error vectors max_error_DIIS_TCSCF = maxval(abs(FQS_SQF_mo)) + ! --- + + delta_energy_tmp = TC_HF_energy - energy_TCSCF_previous + delta_gradie_tmp = grad_non_hermit - gradie_TCSCF_previous + + ! --- + + do while((dabs(delta_energy_tmp) > 0.1d0) .and. (iteration_TCSCF > 1)) +! print *, ' very big step : ', delta_energy_tmp +! print *, ' TC level shift = ', level_shift_TCSCF + + mo_l_coef(1:ao_num,1:mo_num) = mo_l_coef_save(1:ao_num,1:mo_num) + mo_r_coef(1:ao_num,1:mo_num) = mo_r_coef_save(1:ao_num,1:mo_num) + + if(level_shift_TCSCF <= .1d0) then + level_shift_TCSCF = 1.d0 + else + level_shift_TCSCF = level_shift_TCSCF * 3.0d0 + endif + TOUCH mo_l_coef mo_r_coef level_shift_TCSCF + + mo_l_coef(1:ao_num,1:mo_num) = fock_tc_leigvec_ao(1:ao_num,1:mo_num) + mo_r_coef(1:ao_num,1:mo_num) = fock_tc_reigvec_ao(1:ao_num,1:mo_num) + TOUCH mo_l_coef mo_r_coef + + delta_energy_tmp = TC_HF_energy - energy_TCSCF_previous + + if(level_shift_TCSCF - level_shift_save > 40.d0) then + level_shift_TCSCF = level_shift_save * 4.d0 + SOFT_TOUCH level_shift_TCSCF + exit + endif + + dim_DIIS = 0 + enddo +! print *, ' very big step : ', delta_energy_tmp +! print *, ' TC level shift = ', level_shift_TCSCF + + ! --- + + level_shift_TCSCF = 0.d0 + !level_shift_TCSCF = level_shift_TCSCF * 0.5d0 + SOFT_TOUCH level_shift_TCSCF + + gradie_TCSCF = grad_non_hermit energy_TCSCF = TC_HF_energy energy_TCSCF_1e = TC_HF_one_e_energy energy_TCSCF_2e = TC_HF_two_e_energy @@ -125,78 +166,17 @@ subroutine rh_tcscf() if(three_body_h_tc) then energy_TCSCF_3e = diag_three_elem_hf endif - gradie_TCSCF = grad_non_hermit delta_energy_TCSCF = energy_TCSCF - energy_TCSCF_previous delta_gradie_TCSCF = gradie_TCSCF - gradie_TCSCF_previous - if((TCSCF_algorithm == 'DIIS') .and. (delta_gradie_TCSCF > 0.d0)) then - Fock_matrix_tc_ao_tot(1:ao_num,1:ao_num) = F_DIIS(1:ao_num,1:ao_num,index_dim_DIIS) - Fock_matrix_tc_ao_alpha = 0.5d0 * Fock_matrix_tc_ao_tot - Fock_matrix_tc_ao_beta = 0.5d0 * Fock_matrix_tc_ao_tot - TOUCH Fock_matrix_tc_ao_alpha Fock_matrix_tc_ao_beta - endif - - ! --- - - level_shift_TCSCF_save = level_shift_TCSCF - - mo_r_coef_save(1:ao_num,1:mo_num) = mo_r_coef(1:ao_num,1:mo_num) - mo_l_coef_save(1:ao_num,1:mo_num) = mo_l_coef(1:ao_num,1:mo_num) - - do while(delta_gradie_TCSCF > 0.d0) - - mo_r_coef(1:ao_num,1:mo_num) = mo_r_coef_save(1:ao_num,1:mo_num) - mo_l_coef(1:ao_num,1:mo_num) = mo_l_coef_save(1:ao_num,1:mo_num) - - if(level_shift_TCSCF <= .1d0) then - level_shift_TCSCF = 1.d0 - else - level_shift_TCSCF = level_shift_TCSCF * 3.0d0 - endif - TOUCH mo_r_coef mo_l_coef level_shift_TCSCF - - mo_l_coef(1:ao_num,1:mo_num) = fock_tc_leigvec_ao(1:ao_num,1:mo_num) - mo_r_coef(1:ao_num,1:mo_num) = fock_tc_reigvec_ao(1:ao_num,1:mo_num) - - !if(frozen_orb_scf) then - ! call reorder_core_orb - ! call initialize_mo_coef_begin_iteration - !endif - TOUCH mo_l_coef mo_r_coef - - energy_TCSCF = TC_HF_energy - energy_TCSCF_1e = TC_HF_one_e_energy - energy_TCSCF_2e = TC_HF_two_e_energy - energy_TCSCF_3e = 0.d0 - if(three_body_h_tc) then - energy_TCSCF_3e = diag_three_elem_hf - endif - gradie_TCSCF = grad_non_hermit - delta_energy_TCSCF = energy_TCSCF - energy_TCSCF_previous - delta_gradie_TCSCF = gradie_TCSCF - gradie_TCSCF_previous - - if(level_shift_TCSCF - level_shift_TCSCF_save > 40.d0) then - level_shift_TCSCF = level_shift_TCSCF_save * 4.d0 - SOFT_TOUCH level_shift_TCSCF - exit - endif - - dim_DIIS = 0 - enddo - - ! --- - - level_shift_TCSCF = level_shift_TCSCF * 0.5d0 - SOFT_TOUCH level_shift_TCSCF - energy_TCSCF_previous = energy_TCSCF - energy_TCSCF_1e = TC_HF_one_e_energy - energy_TCSCF_2e = TC_HF_two_e_energy - energy_TCSCF_3e = 0.d0 - if(three_body_h_tc) then - energy_TCSCF_3e = diag_three_elem_hf - endif - gradie_TCSCF_previous = grad_non_hermit + gradie_TCSCF_previous = gradie_TCSCF + + + level_shift_save = level_shift_TCSCF + mo_l_coef_save(1:ao_num,1:mo_num) = mo_l_coef(1:ao_num,1:mo_num) + mo_r_coef_save(1:ao_num,1:mo_num) = mo_r_coef(1:ao_num,1:mo_num) + print *, ' iteration = ', iteration_TCSCF print *, ' total TC energy = ', energy_TCSCF @@ -204,36 +184,25 @@ subroutine rh_tcscf() print *, ' 2-e TC energy = ', energy_TCSCF_2e print *, ' 3-e TC energy = ', energy_TCSCF_3e print *, ' |delta TC energy| = ', delta_energy_TCSCF + print *, ' TC gradient = ', gradie_TCSCF print *, ' delta TC gradient = ', delta_gradie_TCSCF print *, ' max TC DIIS error = ', max_error_DIIS_TCSCF print *, ' TC DIIS dim = ', dim_DIIS print *, ' TC level shift = ', level_shift_TCSCF + print *, ' ' - if(delta_gradie_TCSCF < 0.d0) then - call ezfio_set_bi_ortho_mos_mo_l_coef(mo_l_coef) - call ezfio_set_bi_ortho_mos_mo_r_coef(mo_r_coef) - call ezfio_set_tc_scf_bitc_energy(energy_TCSCF) - endif + call ezfio_set_bi_ortho_mos_mo_l_coef(mo_l_coef) + call ezfio_set_bi_ortho_mos_mo_r_coef(mo_r_coef) if(qp_stop()) exit enddo ! --- - !if(iteration_TCSCF < n_it_TCSCF_max) then - ! mo_label = 'Canonical' - !endif - - !if(.not.frozen_orb_scf) then - ! call mo_as_eigvectors_of_mo_matrix(Fock_matrix_mo, size(Fock_matrix_mo,1), size(Fock_matrix_mo, 2), mo_label, 1, .true.) - ! call restore_symmetry(ao_num, mo_num, mo_coef, size(mo_coef, 1), 1.d-10) - ! call orthonormalize_mos - ! call save_mos - !endif - !call write_double(6, energy_TCSCF, 'TCSCF energy') - call write_time(6) + deallocate(mo_r_coef_save, mo_l_coef_save, F_DIIS, e_DIIS) + end ! --- diff --git a/src/tc_scf/routines_rotates.irp.f b/src/tc_scf/routines_rotates.irp.f index 42925e41..15264768 100644 --- a/src/tc_scf/routines_rotates.irp.f +++ b/src/tc_scf/routines_rotates.irp.f @@ -116,7 +116,7 @@ subroutine routine_save_rotated_mos(thr_deg, good_angles) print *, ' ------------------------------------' call orthog_functions(ao_num, n_degen, mo_l_coef_tmp, ao_overlap) - print *, ' Overlap lef-right ' + print *, ' Overlap left-right ' call build_s_matrix(ao_num, n_degen, mo_r_coef_tmp, mo_l_coef_tmp, ao_overlap, stmp) do j = 1, n_degen write(*,'(100(F8.4,X))') stmp(:,j) diff --git a/src/tc_scf/tc_scf.irp.f b/src/tc_scf/tc_scf.irp.f index 2b751e50..283ec2ae 100644 --- a/src/tc_scf/tc_scf.irp.f +++ b/src/tc_scf/tc_scf.irp.f @@ -15,8 +15,8 @@ program tc_scf ! 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 - !call create_guess - !call orthonormalize_mos + call create_guess() + call orthonormalize_mos() PROVIDE tcscf_algorithm if(tcscf_algorithm == 'DIIS') then @@ -42,7 +42,8 @@ subroutine create_guess logical :: exists PROVIDE ezfio_filename - call ezfio_has_mo_basis_mo_coef(exists) + !call ezfio_has_mo_basis_mo_coef(exists) + exists = .false. if (.not.exists) then mo_label = 'Guess' @@ -106,7 +107,7 @@ subroutine simple_tcscf() else - print*,'grad_hermit = ',grad_hermit + print *, ' grad_hermit = ', grad_hermit call save_good_hermit_tc_eigvectors TOUCH mo_coef call save_mos @@ -117,58 +118,70 @@ subroutine simple_tcscf() if(bi_ortho) then - !do while( it .lt. n_it_tcscf_max .and. (e_delta .gt. dsqrt(thresh_tcscf)) ) - !do while( it .lt. n_it_tcscf_max .and. (e_delta .gt. thresh_tcscf) ) - !do while( it .lt. n_it_tcscf_max .and. (rho_delta .gt. thresh_tcscf) ) - do while( it .lt. n_it_tcscf_max .and. (grad_non_hermit_right.gt. dsqrt(thresh_tcscf)) ) + !do while(e_delta .gt. dsqrt(thresh_tcscf)) ) + !do while(e_delta .gt. thresh_tcscf) ) + !do while(rho_delta .gt. thresh_tcscf) ) + !do while(grad_non_hermit_right .gt. dsqrt(thresh_tcscf)) + do while(grad_non_hermit .gt. dsqrt(thresh_tcscf)) it += 1 - print*,'iteration = ', it - print*,'***' - print*,'TC HF total energy = ', TC_HF_energy - print*,'TC HF 1 e energy = ', TC_HF_one_e_energy - print*,'TC HF 2 non hermit = ', TC_HF_two_e_energy - if(three_body_h_tc)then - print*,'TC HF 3 body = ', diag_three_elem_hf + if(it > n_it_tcscf_max) then + print *, ' max of TCSCF iterations is reached ', n_it_TCSCF_max + exit endif - print*,'***' - e_delta = dabs( TC_HF_energy - e_save ) - print*, 'it, delta E = ', it, e_delta - print*, 'it, gradient= ',grad_non_hermit_right + + + print *, ' ***' + print *, ' iteration = ', it + + print *, ' TC HF total energy = ', TC_HF_energy + print *, ' TC HF 1 e energy = ', TC_HF_one_e_energy + print *, ' TC HF 2 non hermit = ', TC_HF_two_e_energy + if(three_body_h_tc) then + print *, ' TC HF 3 body = ', diag_three_elem_hf + endif + e_delta = dabs(TC_HF_energy - e_save) + + print *, ' delta E = ', e_delta + print *, ' gradient = ', grad_non_hermit + !print *, ' gradient= ', grad_non_hermit_right + + !rho_new = TCSCF_bi_ort_dm_ao + !!print*, rho_new + !rho_delta = 0.d0 + !do i = 1, ao_num + ! do j = 1, ao_num + ! rho_delta += dabs(rho_new(j,i) - rho_old(j,i)) + ! enddo + !enddo + !print *, ' rho_delta =', rho_delta + !rho_old = rho_new + e_save = TC_HF_energy mo_l_coef = fock_tc_leigvec_ao mo_r_coef = fock_tc_reigvec_ao - - rho_new = TCSCF_bi_ort_dm_ao - !print*, rho_new - rho_delta = 0.d0 - do i = 1, ao_num - do j = 1, ao_num - rho_delta += dabs(rho_new(j,i) - rho_old(j,i)) - enddo - enddo - print*, ' rho_delta =', rho_delta - rho_old = rho_new - call ezfio_set_bi_ortho_mos_mo_l_coef(mo_l_coef) call ezfio_set_bi_ortho_mos_mo_r_coef(mo_r_coef) TOUCH mo_l_coef mo_r_coef - call ezfio_set_tc_scf_bitc_energy(TC_HF_energy) + print *, ' ***' + print *, '' + enddo else do while( (grad_hermit.gt.dsqrt(thresh_tcscf)) .and. it .lt. n_it_tcscf_max ) print*,'grad_hermit = ',grad_hermit it += 1 - print*,'iteration = ', it - print*,'***' - print*,'TC HF total energy = ', TC_HF_energy - print*,'TC HF 1 e energy = ', TC_HF_one_e_energy - print*,'TC HF 2 e energy = ', TC_HF_two_e_energy - print*,'TC HF 3 body = ', diag_three_elem_hf - print*,'***' + print *, 'iteration = ', it + print *, '***' + print *, 'TC HF total energy = ', TC_HF_energy + print *, 'TC HF 1 e energy = ', TC_HF_one_e_energy + print *, 'TC HF 2 e energy = ', TC_HF_two_e_energy + print *, 'TC HF 3 body = ', diag_three_elem_hf + print *, '***' + print *, '' call save_good_hermit_tc_eigvectors TOUCH mo_coef call save_mos