mirror of
https://github.com/QuantumPackage/qp2.git
synced 2024-09-01 21:53:40 +02:00
minor modifs the minimize_tc_angles routines
This commit is contained in:
parent
a402d509a9
commit
c2a8a17572
@ -1281,10 +1281,10 @@ subroutine impose_orthog_svd_overlap(n, m, C,overlap)
|
|||||||
, C, size(C, 1), Stmp, size(Stmp, 1) &
|
, C, size(C, 1), Stmp, size(Stmp, 1) &
|
||||||
, 0.d0, S, size(S, 1) )
|
, 0.d0, S, size(S, 1) )
|
||||||
|
|
||||||
print *, ' eigenvec overlap bef SVD: '
|
! print *, ' eigenvec overlap bef SVD: '
|
||||||
do i = 1, m
|
! do i = 1, m
|
||||||
write(*, '(1000(F16.10,X))') S(i,:)
|
! write(*, '(1000(F16.10,X))') S(i,:)
|
||||||
enddo
|
! enddo
|
||||||
|
|
||||||
! ---
|
! ---
|
||||||
|
|
||||||
@ -1340,10 +1340,10 @@ subroutine impose_orthog_svd_overlap(n, m, C,overlap)
|
|||||||
, C, size(C, 1), C, size(C, 1) &
|
, C, size(C, 1), C, size(C, 1) &
|
||||||
, 0.d0, S, size(S, 1) )
|
, 0.d0, S, size(S, 1) )
|
||||||
|
|
||||||
print *, ' eigenvec overlap aft SVD: '
|
! print *, ' eigenvec overlap aft SVD: '
|
||||||
do i = 1, m
|
! do i = 1, m
|
||||||
write(*, '(1000(F16.10,X))') S(i,:)
|
! write(*, '(1000(F16.10,X))') S(i,:)
|
||||||
enddo
|
! enddo
|
||||||
|
|
||||||
deallocate(S)
|
deallocate(S)
|
||||||
|
|
||||||
@ -2516,7 +2516,7 @@ subroutine impose_biorthog_svd_overlap(n, m, overlap, L, R)
|
|||||||
|
|
||||||
print *, ' overlap bef SVD: '
|
print *, ' overlap bef SVD: '
|
||||||
do i = 1, m
|
do i = 1, m
|
||||||
write(*, '(1000(F16.10,X))') S(i,:)
|
write(*, '(1000(F25.16,X))') S(i,:)
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
! ---
|
! ---
|
||||||
@ -2530,6 +2530,7 @@ subroutine impose_biorthog_svd_overlap(n, m, overlap, L, R)
|
|||||||
threshold = 1.d-6
|
threshold = 1.d-6
|
||||||
num_linear_dependencies = 0
|
num_linear_dependencies = 0
|
||||||
do i = 1, m
|
do i = 1, m
|
||||||
|
print*,'D(i) = ',D(i)
|
||||||
if(abs(D(i)) <= threshold) then
|
if(abs(D(i)) <= threshold) then
|
||||||
D(i) = 0.d0
|
D(i) = 0.d0
|
||||||
num_linear_dependencies += 1
|
num_linear_dependencies += 1
|
||||||
@ -2585,11 +2586,18 @@ subroutine impose_biorthog_svd_overlap(n, m, overlap, L, R)
|
|||||||
! ---
|
! ---
|
||||||
|
|
||||||
allocate(S(m,m))
|
allocate(S(m,m))
|
||||||
|
! call dgemm( 'T', 'N', m, m, n, 1.d0 &
|
||||||
|
! , L, size(L, 1), R, size(R, 1) &
|
||||||
|
! , 0.d0, S, size(S, 1) )
|
||||||
|
! S = C.T x overlap x C
|
||||||
|
call dgemm( 'N', 'N', n, m, n, 1.d0 &
|
||||||
|
, overlap, size(overlap, 1), R, size(R, 1) &
|
||||||
|
, 0.d0, Stmp, size(Stmp, 1) )
|
||||||
call dgemm( 'T', 'N', m, m, n, 1.d0 &
|
call dgemm( 'T', 'N', m, m, n, 1.d0 &
|
||||||
, L, size(L, 1), R, size(R, 1) &
|
, L, size(L, 1), Stmp, size(Stmp, 1) &
|
||||||
, 0.d0, S, size(S, 1) )
|
, 0.d0, S, size(S, 1) )
|
||||||
|
|
||||||
print *, ' overlap aft SVD: '
|
print *, ' overlap aft SVD with overlap: '
|
||||||
do i = 1, m
|
do i = 1, m
|
||||||
write(*, '(1000(F16.10,X))') S(i,:)
|
write(*, '(1000(F16.10,X))') S(i,:)
|
||||||
enddo
|
enddo
|
||||||
|
@ -131,3 +131,30 @@ END_PROVIDER
|
|||||||
|
|
||||||
! ---
|
! ---
|
||||||
|
|
||||||
|
BEGIN_PROVIDER [ double precision, grad_non_hermit_left]
|
||||||
|
&BEGIN_PROVIDER [ double precision, grad_non_hermit_right]
|
||||||
|
&BEGIN_PROVIDER [ double precision, grad_non_hermit]
|
||||||
|
implicit none
|
||||||
|
integer :: i, k
|
||||||
|
grad_non_hermit_left = 0.d0
|
||||||
|
grad_non_hermit_right = 0.d0
|
||||||
|
do i = 1, elec_beta_num ! doc --> SOMO
|
||||||
|
do k = elec_beta_num+1, elec_alpha_num
|
||||||
|
grad_non_hermit_left+= dabs(Fock_matrix_tc_mo_tot(k,i))
|
||||||
|
grad_non_hermit_right+= dabs(Fock_matrix_tc_mo_tot(i,k))
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
do i = 1, elec_beta_num ! doc --> virt
|
||||||
|
do k = elec_alpha_num+1, mo_num
|
||||||
|
grad_non_hermit_left+= dabs(Fock_matrix_tc_mo_tot(k,i))
|
||||||
|
grad_non_hermit_right+= dabs(Fock_matrix_tc_mo_tot(i,k))
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
do i = elec_beta_num+1, elec_alpha_num ! SOMO --> virt
|
||||||
|
do k = elec_alpha_num+1, mo_num
|
||||||
|
grad_non_hermit_left+= dabs(Fock_matrix_tc_mo_tot(k,i))
|
||||||
|
grad_non_hermit_right+= dabs(Fock_matrix_tc_mo_tot(i,k))
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
grad_non_hermit = grad_non_hermit_left + grad_non_hermit_right
|
||||||
|
END_PROVIDER
|
||||||
|
10
src/tc_scf/minimize_tc_angles.irp.f
Normal file
10
src/tc_scf/minimize_tc_angles.irp.f
Normal file
@ -0,0 +1,10 @@
|
|||||||
|
program print_angles
|
||||||
|
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 = 14 ! small grid for quick debug
|
||||||
|
call minimize_tc_orb_angles
|
||||||
|
end
|
||||||
|
|
@ -5,16 +5,5 @@ program print_angles
|
|||||||
! my_n_pt_a_grid = 50
|
! my_n_pt_a_grid = 50
|
||||||
my_n_pt_r_grid = 10 ! small grid for quick debug
|
my_n_pt_r_grid = 10 ! small grid for quick debug
|
||||||
my_n_pt_a_grid = 14 ! small grid for quick debug
|
my_n_pt_a_grid = 14 ! small grid for quick debug
|
||||||
call routine
|
call print_angles_tc
|
||||||
end
|
|
||||||
subroutine routine
|
|
||||||
implicit none
|
|
||||||
integer :: i,j
|
|
||||||
double precision :: left,right
|
|
||||||
print*,'energy,product of norms, angle between vectors'
|
|
||||||
do i = 1, mo_num
|
|
||||||
left = overlap_mo_l(i,i)
|
|
||||||
right = overlap_mo_r(i,i)
|
|
||||||
print*,Fock_matrix_tc_mo_tot(i,i),left*right,angle_left_right(i)
|
|
||||||
enddo
|
|
||||||
end
|
end
|
||||||
|
@ -1,10 +1,44 @@
|
|||||||
subroutine routine_save_rotated_mos
|
subroutine minimize_tc_orb_angles
|
||||||
implicit none
|
implicit none
|
||||||
|
double precision :: thr_deg
|
||||||
|
logical :: good_angles
|
||||||
|
integer :: i
|
||||||
|
good_angles = .False.
|
||||||
|
thr_deg = thr_degen_tc
|
||||||
|
call print_energy_and_mos
|
||||||
|
i = 1
|
||||||
|
do while (.not. good_angles)
|
||||||
|
print*,'iteration = ',i
|
||||||
|
call routine_save_rotated_mos(thr_deg,good_angles)
|
||||||
|
thr_deg *= 10.d0
|
||||||
|
i+=1
|
||||||
|
if(i.gt.100)then
|
||||||
|
print*,'minimize_tc_orb_angles does not seem to converge ..'
|
||||||
|
print*,'Something is weird in the tc orbitals ...'
|
||||||
|
print*,'STOPPING'
|
||||||
|
endif
|
||||||
|
enddo
|
||||||
|
print*,'Converged ANGLES MINIMIZATION !!'
|
||||||
|
call print_angles_tc
|
||||||
|
call print_energy_and_mos
|
||||||
|
end
|
||||||
|
|
||||||
|
subroutine routine_save_rotated_mos(thr_deg,good_angles)
|
||||||
|
implicit none
|
||||||
|
double precision, intent(in) :: thr_deg
|
||||||
|
logical, intent(out) :: good_angles
|
||||||
|
good_angles = .False.
|
||||||
integer :: i,j,k,n_degen_list,m,n,n_degen,ilast,ifirst
|
integer :: i,j,k,n_degen_list,m,n,n_degen,ilast,ifirst
|
||||||
double precision, allocatable :: mo_r_coef_good(:,:),mo_l_coef_good(:,:)
|
double precision, allocatable :: mo_r_coef_good(:,:),mo_l_coef_good(:,:)
|
||||||
allocate(mo_l_coef_good(ao_num, mo_num), mo_r_coef_good(ao_num,mo_num))
|
allocate(mo_l_coef_good(ao_num, mo_num), mo_r_coef_good(ao_num,mo_num))
|
||||||
double precision, allocatable :: mo_r_coef_new(:,:)
|
double precision, allocatable :: mo_r_coef_new(:,:)
|
||||||
double precision :: norm
|
double precision :: norm
|
||||||
|
print*,'***************************************'
|
||||||
|
print*,'***************************************'
|
||||||
|
print*,'THRESHOLD FOR DEGENERACIES ::: ',thr_deg
|
||||||
|
print*,'***************************************'
|
||||||
|
print*,'***************************************'
|
||||||
|
print*,'Starting with the following TC energy gradient :',grad_non_hermit
|
||||||
mo_r_coef_good = mo_r_coef
|
mo_r_coef_good = mo_r_coef
|
||||||
mo_l_coef_good = mo_l_coef
|
mo_l_coef_good = mo_l_coef
|
||||||
allocate(mo_r_coef_new(ao_num, mo_num))
|
allocate(mo_r_coef_new(ao_num, mo_num))
|
||||||
@ -19,22 +53,19 @@ subroutine routine_save_rotated_mos
|
|||||||
integer, allocatable :: list_degen(:,:)
|
integer, allocatable :: list_degen(:,:)
|
||||||
allocate(list_degen(2,mo_num),s_mat(mo_num,mo_num),fock_diag(mo_num))
|
allocate(list_degen(2,mo_num),s_mat(mo_num,mo_num),fock_diag(mo_num))
|
||||||
do i = 1, mo_num
|
do i = 1, mo_num
|
||||||
fock_diag(i) = fock_matrix_mo(i,i)
|
fock_diag(i) = Fock_matrix_tc_mo_tot(i,i)
|
||||||
enddo
|
enddo
|
||||||
! compute the overlap between the left and rescaled right
|
! compute the overlap between the left and rescaled right
|
||||||
call build_s_matrix(ao_num,mo_num,mo_r_coef_new,mo_r_coef_new,ao_overlap,s_mat)
|
call build_s_matrix(ao_num,mo_num,mo_r_coef_new,mo_r_coef_new,ao_overlap,s_mat)
|
||||||
call give_degen(fock_diag,mo_num,thr_degen_tc,list_degen,n_degen_list)
|
call give_degen(fock_diag,mo_num,thr_deg,list_degen,n_degen_list)
|
||||||
print*,'fock_matrix_mo'
|
print*,'fock_matrix_mo'
|
||||||
do i = 1, mo_num
|
do i = 1, mo_num
|
||||||
print*,i,fock_diag(i),angle_left_right(i)
|
print*,i,fock_diag(i),angle_left_right(i)
|
||||||
enddo
|
enddo
|
||||||
print*,'Overlap '
|
|
||||||
do i = 1, mo_num
|
|
||||||
write(*,'(I2,X,100(F8.4,X))')i,s_mat(:,i)
|
|
||||||
enddo
|
|
||||||
|
|
||||||
do i = 1, n_degen_list
|
do i = 1, n_degen_list
|
||||||
ifirst = list_degen(1,i)
|
ifirst = list_degen(1,i)
|
||||||
|
! if(ifirst.ne.12)cycle
|
||||||
ilast = list_degen(2,i)
|
ilast = list_degen(2,i)
|
||||||
n_degen = ilast - ifirst +1
|
n_degen = ilast - ifirst +1
|
||||||
print*,'ifirst,n_degen = ',ifirst,n_degen
|
print*,'ifirst,n_degen = ',ifirst,n_degen
|
||||||
@ -48,16 +79,29 @@ subroutine routine_save_rotated_mos
|
|||||||
mo_l_coef_tmp(1:ao_num,j) = mo_l_coef(1:ao_num,j+ifirst-1)
|
mo_l_coef_tmp(1:ao_num,j) = mo_l_coef(1:ao_num,j+ifirst-1)
|
||||||
enddo
|
enddo
|
||||||
! Orthogonalization of right functions
|
! Orthogonalization of right functions
|
||||||
print*,'Orthogonalization of right functions'
|
print*,'Orthogonalization of RIGHT functions'
|
||||||
|
print*,'------------------------------------'
|
||||||
call orthog_functions(ao_num,n_degen,mo_r_coef_tmp,ao_overlap)
|
call orthog_functions(ao_num,n_degen,mo_r_coef_tmp,ao_overlap)
|
||||||
|
|
||||||
! Orthogonalization of left functions
|
! Orthogonalization of left functions
|
||||||
print*,'Orthogonalization of left functions'
|
print*,'Orthogonalization of LEFT functions'
|
||||||
call orthog_functions(ao_num,n_degen,mo_r_coef_tmp,ao_overlap)
|
print*,'------------------------------------'
|
||||||
|
call orthog_functions(ao_num,n_degen,mo_l_coef_tmp,ao_overlap)
|
||||||
print*,'Overlap lef-right '
|
print*,'Overlap lef-right '
|
||||||
call build_s_matrix(ao_num,n_degen,mo_r_coef_tmp,mo_l_coef_tmp,ao_overlap,stmp)
|
call build_s_matrix(ao_num,n_degen,mo_r_coef_tmp,mo_l_coef_tmp,ao_overlap,stmp)
|
||||||
do j = 1, n_degen
|
do j = 1, n_degen
|
||||||
write(*,'(100(F8.4,X))')stmp(:,j)
|
write(*,'(100(F8.4,X))')stmp(:,j)
|
||||||
enddo
|
enddo
|
||||||
|
call build_s_matrix(ao_num,n_degen,mo_l_coef_tmp,mo_l_coef_tmp,ao_overlap,stmp)
|
||||||
|
print*,'LEFT/LEFT OVERLAP '
|
||||||
|
do j = 1, n_degen
|
||||||
|
write(*,'(100(F16.10,X))')stmp(:,j)
|
||||||
|
enddo
|
||||||
|
call build_s_matrix(ao_num,n_degen,mo_r_coef_tmp,mo_r_coef_tmp,ao_overlap,stmp)
|
||||||
|
print*,'RIGHT/RIGHT OVERLAP '
|
||||||
|
do j = 1, n_degen
|
||||||
|
write(*,'(100(F16.10,X))')stmp(:,j)
|
||||||
|
enddo
|
||||||
if(maxovl_tc)then
|
if(maxovl_tc)then
|
||||||
T = 0.d0
|
T = 0.d0
|
||||||
Snew = 0.d0
|
Snew = 0.d0
|
||||||
@ -77,6 +121,16 @@ subroutine routine_save_rotated_mos
|
|||||||
else
|
else
|
||||||
mo_l_coef_new = mo_l_coef_tmp
|
mo_l_coef_new = mo_l_coef_tmp
|
||||||
endif
|
endif
|
||||||
|
call build_s_matrix(ao_num,n_degen,mo_l_coef_new,mo_l_coef_new,ao_overlap,stmp)
|
||||||
|
print*,'LEFT/LEFT OVERLAP '
|
||||||
|
do j = 1, n_degen
|
||||||
|
write(*,'(100(F16.10,X))')stmp(:,j)
|
||||||
|
enddo
|
||||||
|
call build_s_matrix(ao_num,n_degen,mo_r_coef_tmp,mo_r_coef_tmp,ao_overlap,stmp)
|
||||||
|
print*,'RIGHT/RIGHT OVERLAP '
|
||||||
|
do j = 1, n_degen
|
||||||
|
write(*,'(100(F16.10,X))')stmp(:,j)
|
||||||
|
enddo
|
||||||
call impose_biorthog_svd_overlap(ao_num, n_degen, ao_overlap, mo_l_coef_new, mo_r_coef_tmp)
|
call impose_biorthog_svd_overlap(ao_num, n_degen, ao_overlap, mo_l_coef_new, mo_r_coef_tmp)
|
||||||
call build_s_matrix(ao_num,n_degen,mo_l_coef_new,mo_r_coef_tmp,ao_overlap,stmp)
|
call build_s_matrix(ao_num,n_degen,mo_l_coef_new,mo_r_coef_tmp,ao_overlap,stmp)
|
||||||
print*,'LAST OVERLAP '
|
print*,'LAST OVERLAP '
|
||||||
@ -131,6 +185,13 @@ subroutine routine_save_rotated_mos
|
|||||||
call ezfio_set_bi_ortho_mos_mo_l_coef(mo_l_coef)
|
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_bi_ortho_mos_mo_r_coef(mo_r_coef)
|
||||||
TOUCH mo_l_coef mo_r_coef
|
TOUCH mo_l_coef mo_r_coef
|
||||||
|
double precision, allocatable :: new_angles(:)
|
||||||
|
allocate(new_angles(mo_num))
|
||||||
|
new_angles(1:mo_num) = dabs(angle_left_right(1:mo_num))
|
||||||
|
double precision :: max_angle
|
||||||
|
max_angle = maxval(new_angles)
|
||||||
|
good_angles = max_angle.lt.45.d0
|
||||||
|
|
||||||
end
|
end
|
||||||
|
|
||||||
subroutine build_s_matrix(m,n,C1,C2,overlap,smat)
|
subroutine build_s_matrix(m,n,C1,C2,overlap,smat)
|
||||||
@ -175,3 +236,28 @@ subroutine orthog_functions(m,n,coef,overlap)
|
|||||||
write(*,'(100(F16.10,X))')stmp(:,j)
|
write(*,'(100(F16.10,X))')stmp(:,j)
|
||||||
enddo
|
enddo
|
||||||
end
|
end
|
||||||
|
|
||||||
|
subroutine print_angles_tc
|
||||||
|
implicit none
|
||||||
|
integer :: i,j
|
||||||
|
double precision :: left,right
|
||||||
|
print*,'product of norms, angle between vectors'
|
||||||
|
do i = 1, mo_num
|
||||||
|
left = overlap_mo_l(i,i)
|
||||||
|
right = overlap_mo_r(i,i)
|
||||||
|
! print*,Fock_matrix_tc_mo_tot(i,i),left*right,angle_left_right(i)
|
||||||
|
print*,left*right,angle_left_right(i)
|
||||||
|
enddo
|
||||||
|
end
|
||||||
|
|
||||||
|
subroutine print_energy_and_mos
|
||||||
|
implicit none
|
||||||
|
integer :: i
|
||||||
|
print*,''
|
||||||
|
print*,'TC energy = ', TC_HF_energy
|
||||||
|
print*,'TC SCF energy gradient = ',grad_non_hermit
|
||||||
|
print*,'Diag Fock elem, product of left/right norm, angle left/right '
|
||||||
|
do i = 1, mo_num
|
||||||
|
write(*,'(I3,X,100(F16.10,X))')i,Fock_matrix_tc_mo_tot(i,i),overlap_mo_l(i,i)*overlap_mo_r(i,i),angle_left_right(i)
|
||||||
|
enddo
|
||||||
|
end
|
||||||
|
@ -19,7 +19,7 @@ program tc_scf
|
|||||||
!call orthonormalize_mos
|
!call orthonormalize_mos
|
||||||
|
|
||||||
call routine_scf()
|
call routine_scf()
|
||||||
call routine_save_rotated_mos
|
call minimize_tc_orb_angles
|
||||||
call print_energy_and_mos
|
call print_energy_and_mos
|
||||||
|
|
||||||
|
|
||||||
@ -126,6 +126,7 @@ subroutine routine_scf()
|
|||||||
print*,'***'
|
print*,'***'
|
||||||
e_delta = dabs( TC_HF_energy - e_save )
|
e_delta = dabs( TC_HF_energy - e_save )
|
||||||
print*, 'it, delta E = ', it, e_delta
|
print*, 'it, delta E = ', it, e_delta
|
||||||
|
print*, 'it, gradient= ',grad_non_hermit_right
|
||||||
e_save = TC_HF_energy
|
e_save = TC_HF_energy
|
||||||
mo_l_coef = fock_tc_leigvec_ao
|
mo_l_coef = fock_tc_leigvec_ao
|
||||||
mo_r_coef = fock_tc_reigvec_ao
|
mo_r_coef = fock_tc_reigvec_ao
|
||||||
@ -181,13 +182,3 @@ end subroutine routine_scf
|
|||||||
|
|
||||||
! ---
|
! ---
|
||||||
|
|
||||||
subroutine print_energy_and_mos
|
|
||||||
implicit none
|
|
||||||
integer :: i
|
|
||||||
print*,'Energy converged !'
|
|
||||||
print*,'Final TC energy = ', TC_HF_energy
|
|
||||||
print*,'Diag Fock elem, product of left/right norm, angle left/right '
|
|
||||||
do i = 1, mo_num
|
|
||||||
write(*,'(I3,X,100(F16.10,X))')i,Fock_matrix_tc_mo_tot(i,i),overlap_mo_l(i,i)*overlap_mo_r(i,i),angle_left_right(i)
|
|
||||||
enddo
|
|
||||||
end
|
|
||||||
|
@ -1,187 +0,0 @@
|
|||||||
program print_angles
|
|
||||||
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 = 14 ! small grid for quick debug
|
|
||||||
call routine
|
|
||||||
end
|
|
||||||
subroutine routine
|
|
||||||
implicit none
|
|
||||||
integer :: i,j,k,n_degen_list,m,n,n_degen,ilast,ifirst
|
|
||||||
double precision, allocatable :: mo_r_coef_good(:,:),mo_l_coef_good(:,:)
|
|
||||||
allocate(mo_l_coef_good(ao_num, mo_num), mo_r_coef_good(ao_num,mo_num))
|
|
||||||
double precision, allocatable :: mo_r_coef_new(:,:)
|
|
||||||
double precision :: norm
|
|
||||||
mo_r_coef_good = mo_r_coef
|
|
||||||
mo_l_coef_good = mo_l_coef
|
|
||||||
allocate(mo_r_coef_new(ao_num, mo_num))
|
|
||||||
mo_r_coef_new = mo_r_coef
|
|
||||||
do i = 1, mo_num
|
|
||||||
norm = 1.d0/dsqrt(overlap_mo_r(i,i))
|
|
||||||
do j = 1, ao_num
|
|
||||||
mo_r_coef_new(j,i) *= norm
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
double precision, allocatable :: fock_diag(:),s_mat(:,:)
|
|
||||||
integer, allocatable :: list_degen(:,:)
|
|
||||||
allocate(list_degen(2,mo_num),s_mat(mo_num,mo_num),fock_diag(mo_num))
|
|
||||||
do i = 1, mo_num
|
|
||||||
fock_diag(i) = fock_matrix_mo(i,i)
|
|
||||||
enddo
|
|
||||||
! compute the overlap between the left and rescaled right
|
|
||||||
call build_s_matrix(ao_num,mo_num,mo_r_coef_new,mo_r_coef_new,ao_overlap,s_mat)
|
|
||||||
call give_degen(fock_diag,mo_num,thr_degen_tc,list_degen,n_degen_list)
|
|
||||||
print*,'fock_matrix_mo'
|
|
||||||
do i = 1, mo_num
|
|
||||||
print*,i,fock_diag(i),angle_left_right(i)
|
|
||||||
enddo
|
|
||||||
print*,'Overlap '
|
|
||||||
do i = 1, mo_num
|
|
||||||
write(*,'(I2,X,100(F8.4,X))')i,s_mat(:,i)
|
|
||||||
enddo
|
|
||||||
|
|
||||||
do i = 1, n_degen_list
|
|
||||||
ifirst = list_degen(1,i)
|
|
||||||
ilast = list_degen(2,i)
|
|
||||||
n_degen = ilast - ifirst +1
|
|
||||||
print*,'ifirst,n_degen = ',ifirst,n_degen
|
|
||||||
double precision, allocatable :: stmp(:,:),T(:,:),Snew(:,:),smat2(:,:)
|
|
||||||
double precision, allocatable :: mo_l_coef_tmp(:,:),mo_r_coef_tmp(:,:),mo_l_coef_new(:,:)
|
|
||||||
allocate(stmp(n_degen,n_degen),smat2(n_degen,n_degen))
|
|
||||||
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,j+ifirst-1)
|
|
||||||
mo_l_coef_tmp(1:ao_num,j) = mo_l_coef(1:ao_num,j+ifirst-1)
|
|
||||||
enddo
|
|
||||||
! Orthogonalization of right functions
|
|
||||||
print*,'Orthogonalization of right functions'
|
|
||||||
call orthog_functions(ao_num,n_degen,mo_r_coef_tmp,ao_overlap)
|
|
||||||
! Orthogonalization of left functions
|
|
||||||
print*,'Orthogonalization of left functions'
|
|
||||||
call orthog_functions(ao_num,n_degen,mo_r_coef_tmp,ao_overlap)
|
|
||||||
print*,'Overlap lef-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)
|
|
||||||
enddo
|
|
||||||
T = 0.d0
|
|
||||||
Snew = 0.d0
|
|
||||||
call maxovl(n_degen, n_degen, stmp, T, Snew)
|
|
||||||
print*,'overlap after'
|
|
||||||
do j = 1, n_degen
|
|
||||||
write(*,'(100(F16.10,X))')Snew(:,j)
|
|
||||||
enddo
|
|
||||||
! mo_l_coef_new = 0.D0
|
|
||||||
! do j = 1, n_degen
|
|
||||||
! do k = 1, n_degen
|
|
||||||
! do m = 1, ao_num
|
|
||||||
! mo_l_coef_new(m,j) += T(k,j) * mo_l_coef_tmp(m,k)
|
|
||||||
! enddo
|
|
||||||
! enddo
|
|
||||||
! enddo
|
|
||||||
call dgemm( 'N', 'N', ao_num, n_degen, n_degen, 1.d0 &
|
|
||||||
, mo_l_coef_tmp, size(mo_l_coef_tmp, 1), T(1,1), size(T, 1) &
|
|
||||||
, 0.d0, mo_l_coef_new, size(mo_l_coef_new, 1) )
|
|
||||||
call build_s_matrix(ao_num,n_degen,mo_l_coef_new,mo_r_coef_tmp,ao_overlap,stmp)
|
|
||||||
print*,'Overlap test'
|
|
||||||
do j = 1, n_degen
|
|
||||||
write(*,'(100(F16.10,X))')stmp(:,j)
|
|
||||||
enddo
|
|
||||||
call impose_biorthog_svd_overlap(ao_num, n_degen, ao_overlap, mo_l_coef_new, mo_r_coef_tmp)
|
|
||||||
call build_s_matrix(ao_num,n_degen,mo_l_coef_new,mo_r_coef_tmp,ao_overlap,stmp)
|
|
||||||
print*,'LAST OVERLAP '
|
|
||||||
do j = 1, n_degen
|
|
||||||
write(*,'(100(F16.10,X))')stmp(:,j)
|
|
||||||
enddo
|
|
||||||
call build_s_matrix(ao_num,n_degen,mo_l_coef_new,mo_l_coef_new,ao_overlap,stmp)
|
|
||||||
print*,'LEFT OVERLAP '
|
|
||||||
do j = 1, n_degen
|
|
||||||
write(*,'(100(F16.10,X))')stmp(:,j)
|
|
||||||
enddo
|
|
||||||
call build_s_matrix(ao_num,n_degen,mo_r_coef_tmp,mo_r_coef_tmp,ao_overlap,stmp)
|
|
||||||
print*,'RIGHT OVERLAP '
|
|
||||||
do j = 1, n_degen
|
|
||||||
write(*,'(100(F16.10,X))')stmp(:,j)
|
|
||||||
enddo
|
|
||||||
do j = 1, n_degen
|
|
||||||
mo_l_coef_good(1:ao_num,j+ifirst-1) = mo_l_coef_new(1:ao_num,j)
|
|
||||||
mo_r_coef_good(1:ao_num,j+ifirst-1) = mo_r_coef_tmp(1:ao_num,j)
|
|
||||||
enddo
|
|
||||||
deallocate(stmp,smat2)
|
|
||||||
deallocate(mo_r_coef_tmp,mo_l_coef_tmp,mo_l_coef_new)
|
|
||||||
deallocate(T,Snew)
|
|
||||||
enddo
|
|
||||||
|
|
||||||
allocate(stmp(mo_num, mo_num))
|
|
||||||
print*,'l coef'
|
|
||||||
do i = 1, mo_num
|
|
||||||
write(*,'(100(F8.4,X))')mo_l_coef_good(:,i)
|
|
||||||
enddo
|
|
||||||
print*,'r coef'
|
|
||||||
do i = 1, mo_num
|
|
||||||
write(*,'(100(F8.4,X))')mo_r_coef_good(:,i)
|
|
||||||
enddo
|
|
||||||
call build_s_matrix(ao_num,mo_num,mo_l_coef_good,mo_r_coef_good,ao_overlap,stmp)
|
|
||||||
print*,'LEFT/RIGHT OVERLAP '
|
|
||||||
do j = 1, mo_num
|
|
||||||
write(*,'(100(F16.10,X))')stmp(:,j)
|
|
||||||
enddo
|
|
||||||
call build_s_matrix(ao_num,mo_num,mo_l_coef_good,mo_l_coef_good,ao_overlap,stmp)
|
|
||||||
print*,'LEFT/LEFT OVERLAP '
|
|
||||||
do j = 1, mo_num
|
|
||||||
write(*,'(100(F16.10,X))')stmp(:,j)
|
|
||||||
enddo
|
|
||||||
call build_s_matrix(ao_num,mo_num,mo_r_coef_good,mo_r_coef_good,ao_overlap,stmp)
|
|
||||||
print*,'RIGHT/RIGHT OVERLAP '
|
|
||||||
do j = 1, mo_num
|
|
||||||
write(*,'(100(F16.10,X))')stmp(:,j)
|
|
||||||
enddo
|
|
||||||
call ezfio_set_bi_ortho_mos_mo_l_coef(mo_l_coef_good)
|
|
||||||
call ezfio_set_bi_ortho_mos_mo_r_coef(mo_r_coef_good)
|
|
||||||
end
|
|
||||||
|
|
||||||
subroutine build_s_matrix(m,n,C1,C2,overlap,smat)
|
|
||||||
implicit none
|
|
||||||
integer, intent(in) :: m,n
|
|
||||||
double precision, intent(in) :: C1(m,n),C2(m,n),overlap(m,m)
|
|
||||||
double precision, intent(out):: smat(n,n)
|
|
||||||
integer :: i,j,k,l
|
|
||||||
smat = 0.D0
|
|
||||||
do i = 1, n
|
|
||||||
do j = 1, n
|
|
||||||
do k = 1, m
|
|
||||||
do l = 1, m
|
|
||||||
smat(i,j) += C1(k,i) * overlap(l,k) * C2(l,j)
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
end
|
|
||||||
|
|
||||||
subroutine orthog_functions(m,n,coef,overlap)
|
|
||||||
implicit none
|
|
||||||
integer, intent(in) :: m,n
|
|
||||||
double precision, intent(in) :: overlap(m,m)
|
|
||||||
double precision, intent(inout) :: coef(m,n)
|
|
||||||
double precision, allocatable :: stmp(:,:)
|
|
||||||
integer :: j
|
|
||||||
allocate(stmp(n,n))
|
|
||||||
call build_s_matrix(m,n,coef,coef,overlap,stmp)
|
|
||||||
print*,'overlap before'
|
|
||||||
do j = 1, n
|
|
||||||
write(*,'(100(F16.10,X))')stmp(:,j)
|
|
||||||
enddo
|
|
||||||
call impose_orthog_svd_overlap(m, n, coef,overlap)
|
|
||||||
call build_s_matrix(m,n,coef,coef,overlap,stmp)
|
|
||||||
do j = 1, n
|
|
||||||
coef(1,:m) *= 1.d0/dsqrt(stmp(j,j))
|
|
||||||
enddo
|
|
||||||
print*,'overlap after'
|
|
||||||
call build_s_matrix(m,n,coef,coef,overlap,stmp)
|
|
||||||
do j = 1, n
|
|
||||||
write(*,'(100(F16.10,X))')stmp(:,j)
|
|
||||||
enddo
|
|
||||||
end
|
|
Loading…
Reference in New Issue
Block a user