LCPQ/qp2
9
1
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-12-21 11:03:29 +01:00
Code Releases Activity

Merge branch 'dev-stable-tc-scf' of https://github.com/AbdAmmar/qp2 into dev-stable-tc-scf

Browse Source
This commit is contained in:
AbdAmmar 2023-06-22 22:07:26 +02:00
commit 49066b2125
11 changed files with 628 additions and 174 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

84
src/ao_basis/aos_in_r.irp.f
View File

@ -65,46 +65,60 @@ double precision function primitive_value(i,j,r)
end end
! ---
subroutine give_all_aos_at_r(r,aos_array) subroutine give_all_aos_at_r(r, tmp_array)
implicit none
BEGIN_dOC
! input : r == r(1) = x and so on
!
! output : aos_array(i) = aos(i) evaluated in $\textbf{r}$
END_DOC
double precision, intent(in) :: r(3)
double precision, intent(out):: aos_array(ao_num)
integer :: power_ao(3) BEGIN_dOC
integer :: i,j,k,l,m !
double precision :: dx,dy,dz,r2 ! input : r == r(1) = x and so on
double precision :: dx2,dy2,dz2 !
double precision :: center_ao(3) ! output : tmp_array(i) = aos(i) evaluated in $\textbf{r}$
double precision :: beta !
do i = 1, nucl_num END_DOC
center_ao(1:3) = nucl_coord(i,1:3)
dx = (r(1) - center_ao(1)) implicit none
dy = (r(2) - center_ao(2)) double precision, intent(in) :: r(3)
dz = (r(3) - center_ao(3)) double precision, intent(out) :: tmp_array(ao_num)
r2 = dx*dx + dy*dy + dz*dz integer :: p_ao(3)
do j = 1,Nucl_N_Aos(i) integer :: i, j, k, l, m
k = Nucl_Aos_transposed(j,i) ! index of the ao in the ordered format double precision :: dx, dy, dz, r2
aos_array(k) = 0.d0 double precision :: dx2, dy2, dz2
power_ao(1:3)= ao_power_ordered_transp_per_nucl(1:3,j,i) double precision :: c_ao(3)
dx2 = dx**power_ao(1) double precision :: beta
dy2 = dy**power_ao(2)
dz2 = dz**power_ao(3) do i = 1, nucl_num
do l = 1,ao_prim_num(k)
beta = ao_expo_ordered_transp_per_nucl(l,j,i) c_ao(1:3) = nucl_coord(i,1:3)
if(dabs(beta*r2).gt.40.d0)cycle dx = r(1) - c_ao(1)
aos_array(k)+= ao_coef_normalized_ordered_transp_per_nucl(l,j,i) * dexp(-beta*r2) dy = r(2) - c_ao(2)
enddo dz = r(3) - c_ao(3)
aos_array(k) = aos_array(k) * dx2 * dy2 * dz2 r2 = dx*dx + dy*dy + dz*dz
do j = 1, Nucl_N_Aos(i)
k = Nucl_Aos_transposed(j,i) ! index of the ao in the ordered format
p_ao(1:3) = ao_power_ordered_transp_per_nucl(1:3,j,i)
dx2 = dx**p_ao(1)
dy2 = dy**p_ao(2)
dz2 = dz**p_ao(3)
tmp_array(k) = 0.d0
do l = 1,ao_prim_num(k)
beta = ao_expo_ordered_transp_per_nucl(l,j,i)
if(dabs(beta*r2).gt.40.d0) cycle
tmp_array(k) += ao_coef_normalized_ordered_transp_per_nucl(l,j,i) * dexp(-beta*r2)
enddo
tmp_array(k) = tmp_array(k) * dx2 * dy2 * dz2
enddo
enddo enddo
enddo
return
end end
! ---
subroutine give_all_aos_and_grad_at_r(r,aos_array,aos_grad_array) subroutine give_all_aos_and_grad_at_r(r,aos_array,aos_grad_array)
implicit none implicit none

38
src/ao_basis/aos_transp.irp.f
View File

@ -1,20 +1,28 @@
BEGIN_PROVIDER [ integer, Nucl_Aos_transposed, (N_AOs_max,nucl_num)]
implicit none ! ---
BEGIN_DOC
! List of AOs attached on each atom BEGIN_PROVIDER [ integer, Nucl_Aos_transposed, (N_AOs_max,nucl_num)]
END_DOC
integer :: i BEGIN_DOC
integer, allocatable :: nucl_tmp(:) ! List of AOs attached on each atom
allocate(nucl_tmp(nucl_num)) END_DOC
nucl_tmp = 0
Nucl_Aos = 0 implicit none
do i = 1, ao_num integer :: i
nucl_tmp(ao_nucl(i))+=1 integer, allocatable :: nucl_tmp(:)
Nucl_Aos_transposed(nucl_tmp(ao_nucl(i)),ao_nucl(i)) = i
enddo allocate(nucl_tmp(nucl_num))
deallocate(nucl_tmp) nucl_tmp = 0
do i = 1, ao_num
nucl_tmp(ao_nucl(i)) += 1
Nucl_Aos_transposed(nucl_tmp(ao_nucl(i)),ao_nucl(i)) = i
enddo
deallocate(nucl_tmp)
END_PROVIDER END_PROVIDER
! ---
BEGIN_PROVIDER [double precision, ao_expo_ordered_transp_per_nucl, (ao_prim_num_max,N_AOs_max,nucl_num) ] BEGIN_PROVIDER [double precision, ao_expo_ordered_transp_per_nucl, (ao_prim_num_max,N_AOs_max,nucl_num) ]
implicit none implicit none
integer :: i,j,k,l integer :: i,j,k,l

2
src/bi_ort_ints/one_e_bi_ort.irp.f
View File

@ -6,7 +6,7 @@ BEGIN_PROVIDER [double precision, ao_one_e_integrals_tc_tot, (ao_num,ao_num)]
implicit none implicit none
integer :: i, j integer :: i, j
ao_one_e_integrals_tc_tot = ao_one_e_integrals ao_one_e_integrals_tc_tot = ao_one_e_integrals
!provide j1b_type !provide j1b_type

150
src/bi_ort_ints/total_twoe_pot.irp.f
View File

@ -80,6 +80,8 @@ BEGIN_PROVIDER [double precision, ao_two_e_tc_tot, (ao_num, ao_num, ao_num, ao_n
enddo enddo
enddo enddo
FREE ao_tc_int_chemist
endif endif
END_PROVIDER END_PROVIDER
@ -128,69 +130,99 @@ BEGIN_PROVIDER [double precision, mo_bi_ortho_tc_two_e_chemist, (mo_num, mo_num,
implicit none implicit none
integer :: i, j, k, l, m, n, p, q integer :: i, j, k, l, m, n, p, q
double precision, allocatable :: mo_tmp_1(:,:,:,:), mo_tmp_2(:,:,:,:) double precision, allocatable :: a1(:,:,:,:), a2(:,:,:,:)
allocate(mo_tmp_1(mo_num,ao_num,ao_num,ao_num)) PROVIDE mo_r_coef mo_l_coef
mo_tmp_1 = 0.d0
do m = 1, ao_num allocate(a2(ao_num,ao_num,ao_num,mo_num))
do p = 1, ao_num
do n = 1, ao_num
do q = 1, ao_num
do k = 1, mo_num
! (k n|p m) = sum_q c_qk * (q n|p m)
mo_tmp_1(k,n,p,m) += mo_l_coef_transp(k,q) * ao_two_e_tc_tot(q,n,p,m)
enddo
enddo
enddo
enddo
enddo
allocate(mo_tmp_2(mo_num,mo_num,ao_num,ao_num)) call dgemm( 'T', 'N', ao_num*ao_num*ao_num, mo_num, ao_num, 1.d0 &
mo_tmp_2 = 0.d0 , ao_two_e_tc_tot(1,1,1,1), ao_num, mo_l_coef(1,1), ao_num &
, 0.d0 , a2(1,1,1,1), ao_num*ao_num*ao_num)
do m = 1, ao_num allocate(a1(ao_num,ao_num,mo_num,mo_num))
do p = 1, ao_num
do n = 1, ao_num
do i = 1, mo_num
do k = 1, mo_num
! (k i|p m) = sum_n c_ni * (k n|p m)
mo_tmp_2(k,i,p,m) += mo_r_coef_transp(i,n) * mo_tmp_1(k,n,p,m)
enddo
enddo
enddo
enddo
enddo
deallocate(mo_tmp_1)
allocate(mo_tmp_1(mo_num,mo_num,mo_num,ao_num)) call dgemm( 'T', 'N', ao_num*ao_num*mo_num, mo_num, ao_num, 1.d0 &
mo_tmp_1 = 0.d0 , a2(1,1,1,1), ao_num, mo_r_coef(1,1), ao_num &
do m = 1, ao_num , 0.d0, a1(1,1,1,1), ao_num*ao_num*mo_num)
do p = 1, ao_num
do l = 1, mo_num
do i = 1, mo_num
do k = 1, mo_num
mo_tmp_1(k,i,l,m) += mo_l_coef_transp(l,p) * mo_tmp_2(k,i,p,m)
enddo
enddo
enddo
enddo
enddo
deallocate(mo_tmp_2)
mo_bi_ortho_tc_two_e_chemist = 0.d0 deallocate(a2)
do m = 1, ao_num allocate(a2(ao_num,mo_num,mo_num,mo_num))
do j = 1, mo_num
do l = 1, mo_num call dgemm( 'T', 'N', ao_num*mo_num*mo_num, mo_num, ao_num, 1.d0 &
do i = 1, mo_num , a1(1,1,1,1), ao_num, mo_l_coef(1,1), ao_num &
do k = 1, mo_num , 0.d0, a2(1,1,1,1), ao_num*mo_num*mo_num)
mo_bi_ortho_tc_two_e_chemist(k,i,l,j) += mo_r_coef_transp(j,m) * mo_tmp_1(k,i,l,m)
enddo deallocate(a1)
enddo
enddo call dgemm( 'T', 'N', mo_num*mo_num*mo_num, mo_num, ao_num, 1.d0 &
enddo , a2(1,1,1,1), ao_num, mo_r_coef(1,1), ao_num &
enddo , 0.d0, mo_bi_ortho_tc_two_e_chemist(1,1,1,1), mo_num*mo_num*mo_num)
deallocate(mo_tmp_1)
deallocate(a2)
!allocate(a1(mo_num,ao_num,ao_num,ao_num))
!a1 = 0.d0
!do m = 1, ao_num
! do p = 1, ao_num
! do n = 1, ao_num
! do q = 1, ao_num
! do k = 1, mo_num
! ! (k n|p m) = sum_q c_qk * (q n|p m)
! a1(k,n,p,m) += mo_l_coef_transp(k,q) * ao_two_e_tc_tot(q,n,p,m)
! enddo
! enddo
! enddo
! enddo
!enddo
!allocate(a2(mo_num,mo_num,ao_num,ao_num))
!a2 = 0.d0
!do m = 1, ao_num
! do p = 1, ao_num
! do n = 1, ao_num
! do i = 1, mo_num
! do k = 1, mo_num
! ! (k i|p m) = sum_n c_ni * (k n|p m)
! a2(k,i,p,m) += mo_r_coef_transp(i,n) * a1(k,n,p,m)
! enddo
! enddo
! enddo
! enddo
!enddo
!deallocate(a1)
!allocate(a1(mo_num,mo_num,mo_num,ao_num))
!a1 = 0.d0
!do m = 1, ao_num
! do p = 1, ao_num
! do l = 1, mo_num
! do i = 1, mo_num
! do k = 1, mo_num
! a1(k,i,l,m) += mo_l_coef_transp(l,p) * a2(k,i,p,m)
! enddo
! enddo
! enddo
! enddo
!enddo
!deallocate(a2)
!mo_bi_ortho_tc_two_e_chemist = 0.d0
!do m = 1, ao_num
! do j = 1, mo_num
! do l = 1, mo_num
! do i = 1, mo_num
! do k = 1, mo_num
! mo_bi_ortho_tc_two_e_chemist(k,i,l,j) += mo_r_coef_transp(j,m) * a1(k,i,l,m)
! enddo
! enddo
! enddo
! enddo
!enddo
!deallocate(a1)
END_PROVIDER END_PROVIDER
@ -209,6 +241,8 @@ BEGIN_PROVIDER [double precision, mo_bi_ortho_tc_two_e, (mo_num, mo_num, mo_num,
implicit none implicit none
integer :: i, j, k, l integer :: i, j, k, l
PROVIDE mo_bi_ortho_tc_two_e_chemist
do j = 1, mo_num do j = 1, mo_num
do i = 1, mo_num do i = 1, mo_num
do l = 1, mo_num do l = 1, mo_num
@ -220,6 +254,8 @@ BEGIN_PROVIDER [double precision, mo_bi_ortho_tc_two_e, (mo_num, mo_num, mo_num,
enddo enddo
enddo enddo
FREE mo_bi_ortho_tc_two_e_chemist
END_PROVIDER END_PROVIDER
! --- ! ---

2
src/bi_ortho_mos/mos_rl.irp.f
View File

@ -136,6 +136,7 @@ BEGIN_PROVIDER [ double precision, mo_r_coef, (ao_num, mo_num) ]
mo_r_coef(j,i) = mo_coef(j,i) mo_r_coef(j,i) = mo_coef(j,i)
enddo enddo
enddo enddo
call ezfio_set_bi_ortho_mos_mo_r_coef(mo_r_coef)
endif endif
END_PROVIDER END_PROVIDER
@ -191,6 +192,7 @@ BEGIN_PROVIDER [ double precision, mo_l_coef, (ao_num, mo_num) ]
mo_l_coef(j,i) = mo_coef(j,i) mo_l_coef(j,i) = mo_coef(j,i)
enddo enddo
enddo enddo
call ezfio_set_bi_ortho_mos_mo_l_coef(mo_l_coef)
endif endif
END_PROVIDER END_PROVIDER

81
src/dft_utils_in_r/ao_in_r.irp.f
View File

@ -1,53 +1,64 @@
BEGIN_PROVIDER[double precision, aos_in_r_array, (ao_num,n_points_final_grid)] ! ---
implicit none
BEGIN_DOC BEGIN_PROVIDER[double precision, aos_in_r_array, (ao_num,n_points_final_grid)]
! aos_in_r_array(i,j) = value of the ith ao on the jth grid point
END_DOC BEGIN_DOC
integer :: i,j ! aos_in_r_array(i,j) = value of the ith ao on the jth grid point
double precision :: aos_array(ao_num), r(3) END_DOC
!$OMP PARALLEL DO &
!$OMP DEFAULT (NONE) & implicit none
!$OMP PRIVATE (i,r,aos_array,j) & integer :: i, j
!$OMP SHARED(aos_in_r_array,n_points_final_grid,ao_num,final_grid_points) double precision :: tmp_array(ao_num), r(3)
do i = 1, n_points_final_grid
r(1) = final_grid_points(1,i) !$OMP PARALLEL DO &
r(2) = final_grid_points(2,i) !$OMP DEFAULT (NONE) &
r(3) = final_grid_points(3,i) !$OMP PRIVATE (i,r,tmp_array,j) &
call give_all_aos_at_r(r,aos_array) !$OMP SHARED(aos_in_r_array,n_points_final_grid,ao_num,final_grid_points)
do j = 1, ao_num do i = 1, n_points_final_grid
aos_in_r_array(j,i) = aos_array(j) r(1) = final_grid_points(1,i)
r(2) = final_grid_points(2,i)
r(3) = final_grid_points(3,i)
call give_all_aos_at_r(r, tmp_array)
do j = 1, ao_num
aos_in_r_array(j,i) = tmp_array(j)
enddo
enddo enddo
enddo !$OMP END PARALLEL DO
!$OMP END PARALLEL DO
END_PROVIDER END_PROVIDER
! ---
BEGIN_PROVIDER[double precision, aos_in_r_array_transp, (n_points_final_grid,ao_num)] BEGIN_PROVIDER[double precision, aos_in_r_array_transp, (n_points_final_grid,ao_num)]
implicit none
BEGIN_DOC BEGIN_DOC
! aos_in_r_array_transp(i,j) = value of the jth ao on the ith grid point ! aos_in_r_array_transp(i,j) = value of the jth ao on the ith grid point
END_DOC END_DOC
integer :: i,j
double precision :: aos_array(ao_num), r(3) implicit none
do i = 1, n_points_final_grid integer :: i, j
do j = 1, ao_num double precision :: aos_array(ao_num), r(3)
aos_in_r_array_transp(i,j) = aos_in_r_array(j,i)
do i = 1, n_points_final_grid
do j = 1, ao_num
aos_in_r_array_transp(i,j) = aos_in_r_array(j,i)
enddo
enddo enddo
enddo
END_PROVIDER END_PROVIDER
! ---
BEGIN_PROVIDER[double precision, aos_grad_in_r_array, (ao_num,n_points_final_grid,3)]
BEGIN_PROVIDER[double precision, aos_grad_in_r_array, (ao_num,n_points_final_grid,3)]
implicit none
BEGIN_DOC BEGIN_DOC
! aos_grad_in_r_array(i,j,k) = value of the kth component of the gradient of ith ao on the jth grid point ! aos_grad_in_r_array(i,j,k) = value of the kth component of the gradient of ith ao on the jth grid point
! !
! k = 1 : x, k= 2, y, k 3, z ! k = 1 : x, k= 2, y, k 3, z
END_DOC END_DOC
implicit none
integer :: i,j,m integer :: i,j,m
double precision :: aos_array(ao_num), r(3) double precision :: aos_array(ao_num), r(3)
double precision :: aos_grad_array(3,ao_num) double precision :: aos_grad_array(3,ao_num)

4
src/non_h_ints_mu/total_tc_int.irp.f
View File

@ -1,4 +1,7 @@
! TODO
! remove ao_two_e_coul and use map directly
! --- ! ---
BEGIN_PROVIDER [double precision, ao_vartc_int_chemist, (ao_num, ao_num, ao_num, ao_num)] BEGIN_PROVIDER [double precision, ao_vartc_int_chemist, (ao_num, ao_num, ao_num, ao_num)]
@ -116,6 +119,7 @@ BEGIN_PROVIDER [double precision, ao_tc_int_chemist_no_cycle, (ao_num, ao_num, a
call wall_time(wall1) call wall_time(wall1)
print *, ' wall time for ao_tc_int_chemist_no_cycle ', wall1 - wall0 print *, ' wall time for ao_tc_int_chemist_no_cycle ', wall1 - wall0
END_PROVIDER END_PROVIDER
! --- ! ---

44
src/tc_bi_ortho/print_tc_dump.irp.f Normal file
View File

@ -0,0 +1,44 @@
program tc_bi_ortho
BEGIN_DOC
! TODO
END_DOC
implicit none
my_grid_becke = .True.
my_n_pt_r_grid = 30
my_n_pt_a_grid = 50
touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
call KMat_tilde_dump()
end
! ---
subroutine KMat_tilde_dump()
implicit none
integer :: i, j, k, l
PROVIDE mo_bi_ortho_tc_two_e_chemist
print *, ' Kmat_tilde in chem notation'
open(33, file='Kmat_tilde.dat', action='write')
do l = 1, mo_num
do k = 1, mo_num
do j = 1, mo_num
do i = 1, mo_num
write(33, '(4(I4, 2X), 4X, E15.7)') i, j, k, l, mo_bi_ortho_tc_two_e_chemist(i,j,k,l)
! TCHint convention
!write(33, '(4(I4, 2X), 4X, E15.7)') i, j, k, l, mo_bi_ortho_tc_two_e_chemist(j,i,l,k)
enddo
enddo
enddo
enddo
close(33)
return
end subroutine KMat_tilde_dump
! ---

13
src/tc_bi_ortho/tc_utils.irp.f
View File

@ -5,27 +5,34 @@ subroutine write_tc_energy()
integer :: i, j, k integer :: i, j, k
double precision :: hmono, htwoe, hthree, htot double precision :: hmono, htwoe, hthree, htot
double precision :: E_TC, O_TC double precision :: E_TC, O_TC
double precision :: E_1e, E_2e, E_3e
do k = 1, n_states do k = 1, n_states
E_TC = 0.d0 E_TC = 0.d0
E_1e = 0.d0
E_2e = 0.d0
E_3e = 0.d0
do i = 1, N_det do i = 1, N_det
do j = 1, N_det do j = 1, N_det
!htot = htilde_matrix_elmt_bi_ortho(i,j)
call htilde_mu_mat_bi_ortho_slow(psi_det(1,1,i), psi_det(1,1,j), N_int, hmono, htwoe, hthree, htot) call htilde_mu_mat_bi_ortho_slow(psi_det(1,1,i), psi_det(1,1,j), N_int, hmono, htwoe, hthree, htot)
E_TC = E_TC + psi_l_coef_bi_ortho(i,k) * psi_r_coef_bi_ortho(j,k) * htot E_TC = E_TC + psi_l_coef_bi_ortho(i,k) * psi_r_coef_bi_ortho(j,k) * htot
!E_TC = E_TC + leigvec_tc_bi_orth(i,k) * reigvec_tc_bi_orth(j,k) * htot E_1e = E_1e + psi_l_coef_bi_ortho(i,k) * psi_r_coef_bi_ortho(j,k) * hmono
E_2e = E_2e + psi_l_coef_bi_ortho(i,k) * psi_r_coef_bi_ortho(j,k) * htwoe
E_3e = E_3e + psi_l_coef_bi_ortho(i,k) * psi_r_coef_bi_ortho(j,k) * hthree
enddo enddo
enddo enddo
O_TC = 0.d0 O_TC = 0.d0
do i = 1, N_det do i = 1, N_det
!O_TC = O_TC + leigvec_tc_bi_orth(i,k) * reigvec_tc_bi_orth(i,k)
O_TC = O_TC + psi_l_coef_bi_ortho(i,k) * psi_r_coef_bi_ortho(i,k) O_TC = O_TC + psi_l_coef_bi_ortho(i,k) * psi_r_coef_bi_ortho(i,k)
enddo enddo
print *, ' state :', k print *, ' state :', k
print *, " E_TC = ", E_TC / O_TC print *, " E_TC = ", E_TC / O_TC
print *, " E_1e = ", E_1e / O_TC
print *, " E_2e = ", E_2e / O_TC
print *, " E_3e = ", E_3e / O_TC
print *, " O_TC = ", O_TC print *, " O_TC = ", O_TC
enddo enddo

52
src/tc_scf/fock_three_hermit.irp.f
View File

@ -1,30 +1,36 @@
! ---
BEGIN_PROVIDER [ double precision, fock_3_mat, (mo_num, mo_num)] BEGIN_PROVIDER [ double precision, fock_3_mat, (mo_num, mo_num)]
implicit none
implicit none
integer :: i,j integer :: i,j
double precision :: contrib double precision :: contrib
fock_3_mat = 0.d0 fock_3_mat = 0.d0
if(.not.bi_ortho.and.three_body_h_tc)then if(.not.bi_ortho .and. three_body_h_tc) then
call give_fock_ia_three_e_total(1,1,contrib)
!! !$OMP PARALLEL & call give_fock_ia_three_e_total(1, 1, contrib)
!! !$OMP DEFAULT (NONE) & !! !$OMP PARALLEL &
!! !$OMP PRIVATE (i,j,m,integral) & !! !$OMP DEFAULT (NONE) &
!! !$OMP SHARED (mo_num,three_body_3_index) !! !$OMP PRIVATE (i,j,m,integral) &
!! !$OMP DO SCHEDULE (guided) COLLAPSE(3) !! !$OMP SHARED (mo_num,three_body_3_index)
do i = 1, mo_num !! !$OMP DO SCHEDULE (guided) COLLAPSE(3)
do j = 1, mo_num do i = 1, mo_num
call give_fock_ia_three_e_total(j,i,contrib) do j = 1, mo_num
fock_3_mat(j,i) = -contrib call give_fock_ia_three_e_total(j,i,contrib)
enddo fock_3_mat(j,i) = -contrib
enddo enddo
else if(bi_ortho.and.three_body_h_tc)then enddo
!! !$OMP END DO !else if(bi_ortho.and.three_body_h_tc) then
!! !$OMP END PARALLEL !! !$OMP END DO
!! do i = 1, mo_num !! !$OMP END PARALLEL
!! do j = 1, i-1 !! do i = 1, mo_num
!! mat_three(j,i) = mat_three(i,j) !! do j = 1, i-1
!! enddo !! mat_three(j,i) = mat_three(i,j)
!! enddo !! enddo
endif !! enddo
endif
END_PROVIDER END_PROVIDER

332
src/tc_scf/molden_lr_mos.irp.f
View File

@ -1,6 +1,9 @@
program molden ! ---
program molden_lr_mos
BEGIN_DOC BEGIN_DOC
! TODO : Put the documentation of the program here ! TODO : Put the documentation of the program here
END_DOC END_DOC
implicit none implicit none
@ -14,13 +17,21 @@ program molden
! my_n_pt_a_grid = 26 ! 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 touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
call molden_lr !call molden_lr
call molden_l()
call molden_r()
end end
! ---
subroutine molden_lr subroutine molden_lr
implicit none
BEGIN_DOC BEGIN_DOC
! Produces a Molden file ! Produces a Molden file
END_DOC END_DOC
implicit none
character*(128) :: output character*(128) :: output
integer :: i_unit_output,getUnitAndOpen integer :: i_unit_output,getUnitAndOpen
integer :: i,j,k,l integer :: i,j,k,l
@ -37,7 +48,7 @@ subroutine molden_lr
write(i_unit_output,'(A)') '[Atoms] Angs' write(i_unit_output,'(A)') '[Atoms] Angs'
do i = 1, nucl_num do i = 1, nucl_num
write(i_unit_output,'(A2,2X,I4,2X,I4,3(2X,F15.10))') & write(i_unit_output,'(A2,2X,I4,2X,I4,3(2X,F15.10))') &
trim(element_name(int(nucl_charge(i)))), & trim(element_name(int(nucl_charge(i)))), &
i, & i, &
int(nucl_charge(i)), & int(nucl_charge(i)), &
@ -174,3 +185,314 @@ subroutine molden_lr
close(i_unit_output) close(i_unit_output)
end end
! ---
subroutine molden_l()
BEGIN_DOC
! Produces a Molden file
END_DOC
implicit none
character*(128) :: output
integer :: i_unit_output, getUnitAndOpen
integer :: i, j, k, l
double precision, parameter :: a0 = 0.529177249d0
PROVIDE ezfio_filename
PROVIDE mo_l_coef
output=trim(ezfio_filename)//'_left.mol'
print*,'output = ',trim(output)
i_unit_output = getUnitAndOpen(output,'w')
write(i_unit_output,'(A)') '[Molden Format]'
write(i_unit_output,'(A)') '[Atoms] Angs'
do i = 1, nucl_num
write(i_unit_output,'(A2,2X,I4,2X,I4,3(2X,F15.10))') &
trim(element_name(int(nucl_charge(i)))), &
i, &
int(nucl_charge(i)), &
nucl_coord(i,1)*a0, nucl_coord(i,2)*a0, nucl_coord(i,3)*a0
enddo
write(i_unit_output,'(A)') '[GTO]'
character*(1) :: character_shell
integer :: i_shell,i_prim,i_ao
integer :: iorder(ao_num)
integer :: nsort(ao_num)
i_shell = 0
i_prim = 0
do i=1,nucl_num
write(i_unit_output,*) i, 0
do j=1,nucl_num_shell_aos(i)
i_shell +=1
i_ao = nucl_list_shell_aos(i,j)
character_shell = trim(ao_l_char(i_ao))
write(i_unit_output,*) character_shell, ao_prim_num(i_ao), '1.00'
do k = 1, ao_prim_num(i_ao)
i_prim +=1
write(i_unit_output,'(E20.10,2X,E20.10)') ao_expo(i_ao,k), ao_coef(i_ao,k)
enddo
l = i_ao
do while ( ao_l(l) == ao_l(i_ao) )
nsort(l) = i*10000 + j*100
l += 1
if (l > ao_num) exit
enddo
enddo
write(i_unit_output,*)''
enddo
do i=1,ao_num
iorder(i) = i
! p
if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 0 )) then
nsort(i) += 1
else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 0 )) then
nsort(i) += 2
else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 1 )) then
nsort(i) += 3
! d
else if ((ao_power(i,1) == 2 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 0 )) then
nsort(i) += 1
else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 2 ).and.(ao_power(i,3) == 0 )) then
nsort(i) += 2
else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 2 )) then
nsort(i) += 3
else if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 0 )) then
nsort(i) += 4
else if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 1 )) then
nsort(i) += 5
else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 1 )) then
nsort(i) += 6
! f
else if ((ao_power(i,1) == 3 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 0 )) then
nsort(i) += 1
else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 3 ).and.(ao_power(i,3) == 0 )) then
nsort(i) += 2
else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 3 )) then
nsort(i) += 3
else if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 2 ).and.(ao_power(i,3) == 0 )) then
nsort(i) += 4
else if ((ao_power(i,1) == 2 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 0 )) then
nsort(i) += 5
else if ((ao_power(i,1) == 2 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 1 )) then
nsort(i) += 6
else if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 2 )) then
nsort(i) += 7
else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 2 )) then
nsort(i) += 8
else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 2 ).and.(ao_power(i,3) == 1 )) then
nsort(i) += 9
else if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 1 )) then
nsort(i) += 10
! g
else if ((ao_power(i,1) == 4 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 0 )) then
nsort(i) += 1
else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 4 ).and.(ao_power(i,3) == 0 )) then
nsort(i) += 2
else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 4 )) then
nsort(i) += 3
else if ((ao_power(i,1) == 3 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 0 )) then
nsort(i) += 4
else if ((ao_power(i,1) == 3 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 1 )) then
nsort(i) += 5
else if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 3 ).and.(ao_power(i,3) == 0 )) then
nsort(i) += 6
else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 3 ).and.(ao_power(i,3) == 1 )) then
nsort(i) += 7
else if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 3 )) then
nsort(i) += 8
else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 3 )) then
nsort(i) += 9
else if ((ao_power(i,1) == 2 ).and.(ao_power(i,2) == 2 ).and.(ao_power(i,3) == 0 )) then
nsort(i) += 10
else if ((ao_power(i,1) == 2 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 2 )) then
nsort(i) += 11
else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 2 ).and.(ao_power(i,3) == 2 )) then
nsort(i) += 12
else if ((ao_power(i,1) == 2 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 1 )) then
nsort(i) += 13
else if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 2 ).and.(ao_power(i,3) == 1 )) then
nsort(i) += 14
else if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 2 )) then
nsort(i) += 15
endif
enddo
call isort(nsort,iorder,ao_num)
write(i_unit_output,'(A)') '[MO]'
do i=1,mo_num
write (i_unit_output,*) 'Sym= 1'
write (i_unit_output,*) 'Ene=', Fock_matrix_tc_mo_tot(i,i)
write (i_unit_output,*) 'Spin= Alpha'
write (i_unit_output,*) 'Occup=', mo_occ(i)
do j=1,ao_num
write(i_unit_output, '(I6,2X,E20.10)') j, mo_l_coef(iorder(j),i)
enddo
enddo
close(i_unit_output)
end
! ---
subroutine molden_r()
BEGIN_DOC
! Produces a Molden file
END_DOC
implicit none
character*(128) :: output
integer :: i_unit_output, getUnitAndOpen
integer :: i, j, k, l
double precision, parameter :: a0 = 0.529177249d0
PROVIDE ezfio_filename
output=trim(ezfio_filename)//'_right.mol'
print*,'output = ',trim(output)
i_unit_output = getUnitAndOpen(output,'w')
write(i_unit_output,'(A)') '[Molden Format]'
write(i_unit_output,'(A)') '[Atoms] Angs'
do i = 1, nucl_num
write(i_unit_output,'(A2,2X,I4,2X,I4,3(2X,F15.10))') &
trim(element_name(int(nucl_charge(i)))), &
i, &
int(nucl_charge(i)), &
nucl_coord(i,1)*a0, nucl_coord(i,2)*a0, nucl_coord(i,3)*a0
enddo
write(i_unit_output,'(A)') '[GTO]'
character*(1) :: character_shell
integer :: i_shell,i_prim,i_ao
integer :: iorder(ao_num)
integer :: nsort(ao_num)
i_shell = 0
i_prim = 0
do i=1,nucl_num
write(i_unit_output,*) i, 0
do j=1,nucl_num_shell_aos(i)
i_shell +=1
i_ao = nucl_list_shell_aos(i,j)
character_shell = trim(ao_l_char(i_ao))
write(i_unit_output,*) character_shell, ao_prim_num(i_ao), '1.00'
do k = 1, ao_prim_num(i_ao)
i_prim +=1
write(i_unit_output,'(E20.10,2X,E20.10)') ao_expo(i_ao,k), ao_coef(i_ao,k)
enddo
l = i_ao
do while ( ao_l(l) == ao_l(i_ao) )
nsort(l) = i*10000 + j*100
l += 1
if (l > ao_num) exit
enddo
enddo
write(i_unit_output,*)''
enddo
do i=1,ao_num
iorder(i) = i
! p
if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 0 )) then
nsort(i) += 1
else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 0 )) then
nsort(i) += 2
else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 1 )) then
nsort(i) += 3
! d
else if ((ao_power(i,1) == 2 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 0 )) then
nsort(i) += 1
else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 2 ).and.(ao_power(i,3) == 0 )) then
nsort(i) += 2
else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 2 )) then
nsort(i) += 3
else if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 0 )) then
nsort(i) += 4
else if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 1 )) then
nsort(i) += 5
else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 1 )) then
nsort(i) += 6
! f
else if ((ao_power(i,1) == 3 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 0 )) then
nsort(i) += 1
else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 3 ).and.(ao_power(i,3) == 0 )) then
nsort(i) += 2
else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 3 )) then
nsort(i) += 3
else if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 2 ).and.(ao_power(i,3) == 0 )) then
nsort(i) += 4
else if ((ao_power(i,1) == 2 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 0 )) then
nsort(i) += 5
else if ((ao_power(i,1) == 2 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 1 )) then
nsort(i) += 6
else if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 2 )) then
nsort(i) += 7
else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 2 )) then
nsort(i) += 8
else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 2 ).and.(ao_power(i,3) == 1 )) then
nsort(i) += 9
else if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 1 )) then
nsort(i) += 10
! g
else if ((ao_power(i,1) == 4 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 0 )) then
nsort(i) += 1
else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 4 ).and.(ao_power(i,3) == 0 )) then
nsort(i) += 2
else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 4 )) then
nsort(i) += 3
else if ((ao_power(i,1) == 3 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 0 )) then
nsort(i) += 4
else if ((ao_power(i,1) == 3 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 1 )) then
nsort(i) += 5
else if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 3 ).and.(ao_power(i,3) == 0 )) then
nsort(i) += 6
else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 3 ).and.(ao_power(i,3) == 1 )) then
nsort(i) += 7
else if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 3 )) then
nsort(i) += 8
else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 3 )) then
nsort(i) += 9
else if ((ao_power(i,1) == 2 ).and.(ao_power(i,2) == 2 ).and.(ao_power(i,3) == 0 )) then
nsort(i) += 10
else if ((ao_power(i,1) == 2 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 2 )) then
nsort(i) += 11
else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 2 ).and.(ao_power(i,3) == 2 )) then
nsort(i) += 12
else if ((ao_power(i,1) == 2 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 1 )) then
nsort(i) += 13
else if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 2 ).and.(ao_power(i,3) == 1 )) then
nsort(i) += 14
else if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 2 )) then
nsort(i) += 15
endif
enddo
call isort(nsort, iorder, ao_num)
write(i_unit_output,'(A)') '[MO]'
do i=1,mo_num
write (i_unit_output,*) 'Sym= 1'
write (i_unit_output,*) 'Ene=', Fock_matrix_tc_mo_tot(i,i)
write (i_unit_output,*) 'Spin= Alpha'
write (i_unit_output,*) 'Occup=', mo_occ(i)
do j=1,ao_num
write(i_unit_output, '(I6,2X,E20.10)') j, mo_r_coef(iorder(j),i)
enddo
enddo
close(i_unit_output)
end