From 5a5071f248ca55516f37b9f09f01f279c7c6faea Mon Sep 17 00:00:00 2001 From: Abdallah Ammar Date: Thu, 22 Jun 2023 18:26:52 +0200 Subject: [PATCH] fixed bug in nucl_aos --- src/ao_basis/aos_in_r.irp.f | 84 ++++--- src/ao_basis/aos_transp.irp.f | 38 +-- src/bi_ort_ints/one_e_bi_ort.irp.f | 2 +- src/bi_ortho_mos/mos_rl.irp.f | 2 + src/dft_utils_in_r/ao_in_r.irp.f | 81 ++++--- src/non_h_ints_mu/total_tc_int.irp.f | 3 + src/tc_bi_ortho/tc_utils.irp.f | 13 +- src/tc_scf/fock_three_hermit.irp.f | 52 +++-- src/tc_scf/molden_lr_mos.irp.f | 332 ++++++++++++++++++++++++++- 9 files changed, 490 insertions(+), 117 deletions(-) diff --git a/src/ao_basis/aos_in_r.irp.f b/src/ao_basis/aos_in_r.irp.f index 7fcb980a..1b1595a3 100644 --- a/src/ao_basis/aos_in_r.irp.f +++ b/src/ao_basis/aos_in_r.irp.f @@ -65,46 +65,60 @@ double precision function primitive_value(i,j,r) end +! --- -subroutine give_all_aos_at_r(r,aos_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) +subroutine give_all_aos_at_r(r, tmp_array) - integer :: power_ao(3) - integer :: i,j,k,l,m - double precision :: dx,dy,dz,r2 - double precision :: dx2,dy2,dz2 - double precision :: center_ao(3) - double precision :: beta - do i = 1, nucl_num - center_ao(1:3) = nucl_coord(i,1:3) - dx = (r(1) - center_ao(1)) - dy = (r(2) - center_ao(2)) - dz = (r(3) - center_ao(3)) - 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 - aos_array(k) = 0.d0 - power_ao(1:3)= ao_power_ordered_transp_per_nucl(1:3,j,i) - dx2 = dx**power_ao(1) - dy2 = dy**power_ao(2) - dz2 = dz**power_ao(3) - 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 - aos_array(k)+= ao_coef_normalized_ordered_transp_per_nucl(l,j,i) * dexp(-beta*r2) - enddo - aos_array(k) = aos_array(k) * dx2 * dy2 * dz2 + BEGIN_dOC + ! + ! input : r == r(1) = x and so on + ! + ! output : tmp_array(i) = aos(i) evaluated in $\textbf{r}$ + ! + END_DOC + + implicit none + double precision, intent(in) :: r(3) + double precision, intent(out) :: tmp_array(ao_num) + integer :: p_ao(3) + integer :: i, j, k, l, m + double precision :: dx, dy, dz, r2 + double precision :: dx2, dy2, dz2 + double precision :: c_ao(3) + double precision :: beta + + do i = 1, nucl_num + + c_ao(1:3) = nucl_coord(i,1:3) + dx = r(1) - c_ao(1) + dy = r(2) - c_ao(2) + dz = r(3) - c_ao(3) + 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 + + return end +! --- subroutine give_all_aos_and_grad_at_r(r,aos_array,aos_grad_array) implicit none diff --git a/src/ao_basis/aos_transp.irp.f b/src/ao_basis/aos_transp.irp.f index ae6193bf..4e44a9f6 100644 --- a/src/ao_basis/aos_transp.irp.f +++ b/src/ao_basis/aos_transp.irp.f @@ -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 - END_DOC - integer :: i - integer, allocatable :: nucl_tmp(:) - allocate(nucl_tmp(nucl_num)) - nucl_tmp = 0 - Nucl_Aos = 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) + +! --- + +BEGIN_PROVIDER [ integer, Nucl_Aos_transposed, (N_AOs_max,nucl_num)] + + BEGIN_DOC + ! List of AOs attached on each atom + END_DOC + + implicit none + integer :: i + integer, allocatable :: nucl_tmp(:) + + allocate(nucl_tmp(nucl_num)) + 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 +! --- + BEGIN_PROVIDER [double precision, ao_expo_ordered_transp_per_nucl, (ao_prim_num_max,N_AOs_max,nucl_num) ] implicit none integer :: i,j,k,l diff --git a/src/bi_ort_ints/one_e_bi_ort.irp.f b/src/bi_ort_ints/one_e_bi_ort.irp.f index 5f2795f1..49181182 100644 --- a/src/bi_ort_ints/one_e_bi_ort.irp.f +++ b/src/bi_ort_ints/one_e_bi_ort.irp.f @@ -6,7 +6,7 @@ BEGIN_PROVIDER [double precision, ao_one_e_integrals_tc_tot, (ao_num,ao_num)] implicit none 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 diff --git a/src/bi_ortho_mos/mos_rl.irp.f b/src/bi_ortho_mos/mos_rl.irp.f index c69309d1..13eedfb7 100644 --- a/src/bi_ortho_mos/mos_rl.irp.f +++ b/src/bi_ortho_mos/mos_rl.irp.f @@ -136,6 +136,7 @@ BEGIN_PROVIDER [ double precision, mo_r_coef, (ao_num, mo_num) ] mo_r_coef(j,i) = mo_coef(j,i) enddo enddo + call ezfio_set_bi_ortho_mos_mo_r_coef(mo_r_coef) endif 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) enddo enddo + call ezfio_set_bi_ortho_mos_mo_l_coef(mo_l_coef) endif END_PROVIDER diff --git a/src/dft_utils_in_r/ao_in_r.irp.f b/src/dft_utils_in_r/ao_in_r.irp.f index b8beea76..16414f39 100644 --- a/src/dft_utils_in_r/ao_in_r.irp.f +++ b/src/dft_utils_in_r/ao_in_r.irp.f @@ -1,53 +1,64 @@ - BEGIN_PROVIDER[double precision, aos_in_r_array, (ao_num,n_points_final_grid)] - implicit none - BEGIN_DOC - ! aos_in_r_array(i,j) = value of the ith ao on the jth grid point - END_DOC - integer :: i,j - double precision :: aos_array(ao_num), r(3) - !$OMP PARALLEL DO & - !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (i,r,aos_array,j) & - !$OMP SHARED(aos_in_r_array,n_points_final_grid,ao_num,final_grid_points) - do i = 1, n_points_final_grid - 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,aos_array) - do j = 1, ao_num - aos_in_r_array(j,i) = aos_array(j) +! --- + +BEGIN_PROVIDER[double precision, aos_in_r_array, (ao_num,n_points_final_grid)] + + BEGIN_DOC + ! aos_in_r_array(i,j) = value of the ith ao on the jth grid point + END_DOC + + implicit none + integer :: i, j + double precision :: tmp_array(ao_num), r(3) + + !$OMP PARALLEL DO & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (i,r,tmp_array,j) & + !$OMP SHARED(aos_in_r_array,n_points_final_grid,ao_num,final_grid_points) + do i = 1, n_points_final_grid + 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 - !$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)] - implicit none - BEGIN_DOC - ! aos_in_r_array_transp(i,j) = value of the jth ao on the ith grid point - END_DOC - integer :: i,j - double precision :: aos_array(ao_num), r(3) - 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) +BEGIN_PROVIDER[double precision, aos_in_r_array_transp, (n_points_final_grid,ao_num)] + + BEGIN_DOC + ! aos_in_r_array_transp(i,j) = value of the jth ao on the ith grid point + END_DOC + + implicit none + integer :: i, j + double precision :: aos_array(ao_num), r(3) + + 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 - 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 ! 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 END_DOC + + implicit none integer :: i,j,m double precision :: aos_array(ao_num), r(3) double precision :: aos_grad_array(3,ao_num) diff --git a/src/non_h_ints_mu/total_tc_int.irp.f b/src/non_h_ints_mu/total_tc_int.irp.f index afa10305..2bdf39f0 100644 --- a/src/non_h_ints_mu/total_tc_int.irp.f +++ b/src/non_h_ints_mu/total_tc_int.irp.f @@ -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)] diff --git a/src/tc_bi_ortho/tc_utils.irp.f b/src/tc_bi_ortho/tc_utils.irp.f index 9023e2f0..53fe5884 100644 --- a/src/tc_bi_ortho/tc_utils.irp.f +++ b/src/tc_bi_ortho/tc_utils.irp.f @@ -5,27 +5,34 @@ subroutine write_tc_energy() integer :: i, j, k double precision :: hmono, htwoe, hthree, htot double precision :: E_TC, O_TC + double precision :: E_1e, E_2e, E_3e do k = 1, n_states E_TC = 0.d0 + E_1e = 0.d0 + E_2e = 0.d0 + E_3e = 0.d0 do i = 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) 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 O_TC = 0.d0 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) enddo print *, ' state :', k 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 enddo diff --git a/src/tc_scf/fock_three_hermit.irp.f b/src/tc_scf/fock_three_hermit.irp.f index fe8fbfd7..89e6f620 100644 --- a/src/tc_scf/fock_three_hermit.irp.f +++ b/src/tc_scf/fock_three_hermit.irp.f @@ -1,30 +1,36 @@ + +! --- + BEGIN_PROVIDER [ double precision, fock_3_mat, (mo_num, mo_num)] - implicit none + + implicit none integer :: i,j double precision :: contrib + fock_3_mat = 0.d0 - if(.not.bi_ortho.and.three_body_h_tc)then - call give_fock_ia_three_e_total(1,1,contrib) -!! !$OMP PARALLEL & -!! !$OMP DEFAULT (NONE) & -!! !$OMP PRIVATE (i,j,m,integral) & -!! !$OMP SHARED (mo_num,three_body_3_index) -!! !$OMP DO SCHEDULE (guided) COLLAPSE(3) - do i = 1, mo_num - do j = 1, mo_num - call give_fock_ia_three_e_total(j,i,contrib) - fock_3_mat(j,i) = -contrib - enddo - enddo - else if(bi_ortho.and.three_body_h_tc)then -!! !$OMP END DO -!! !$OMP END PARALLEL -!! do i = 1, mo_num -!! do j = 1, i-1 -!! mat_three(j,i) = mat_three(i,j) -!! enddo -!! enddo - endif + if(.not.bi_ortho .and. three_body_h_tc) then + + call give_fock_ia_three_e_total(1, 1, contrib) + !! !$OMP PARALLEL & + !! !$OMP DEFAULT (NONE) & + !! !$OMP PRIVATE (i,j,m,integral) & + !! !$OMP SHARED (mo_num,three_body_3_index) + !! !$OMP DO SCHEDULE (guided) COLLAPSE(3) + do i = 1, mo_num + do j = 1, mo_num + call give_fock_ia_three_e_total(j,i,contrib) + fock_3_mat(j,i) = -contrib + enddo + enddo + !else if(bi_ortho.and.three_body_h_tc) then + !! !$OMP END DO + !! !$OMP END PARALLEL + !! do i = 1, mo_num + !! do j = 1, i-1 + !! mat_three(j,i) = mat_three(i,j) + !! enddo + !! enddo + endif END_PROVIDER diff --git a/src/tc_scf/molden_lr_mos.irp.f b/src/tc_scf/molden_lr_mos.irp.f index 735349ba..e12fcd1c 100644 --- a/src/tc_scf/molden_lr_mos.irp.f +++ b/src/tc_scf/molden_lr_mos.irp.f @@ -1,6 +1,9 @@ -program molden +! --- + +program molden_lr_mos + BEGIN_DOC -! TODO : Put the documentation of the program here + ! TODO : Put the documentation of the program here END_DOC implicit none @@ -14,13 +17,21 @@ program molden ! 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 molden_lr + !call molden_lr + call molden_l() + call molden_r() + end + +! --- + subroutine molden_lr - implicit none + BEGIN_DOC ! Produces a Molden file END_DOC + + implicit none character*(128) :: output integer :: i_unit_output,getUnitAndOpen integer :: i,j,k,l @@ -37,7 +48,7 @@ subroutine molden_lr 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))') & + 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)), & @@ -174,3 +185,314 @@ subroutine molden_lr close(i_unit_output) 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 +