From 2bbbfaaae62ac859de63895463350ff6dfc38c48 Mon Sep 17 00:00:00 2001 From: yasuto-masuda Date: Sun, 9 Jul 2023 22:23:10 +0900 Subject: [PATCH] calc E_type2DM from S_type1-2DM --- src/DMRG_interface.f90 | 180 +++++++++++++++++++++++++++++++++++++++-- src/r4dcasci_co.f90 | 9 ++- 2 files changed, 182 insertions(+), 7 deletions(-) diff --git a/src/DMRG_interface.f90 b/src/DMRG_interface.f90 index 09aeb840..2bd753b5 100644 --- a/src/DMRG_interface.f90 +++ b/src/DMRG_interface.f90 @@ -664,17 +664,19 @@ SUBROUTINE bAmat_interface_for_debug(sc, bc, isym, indsym) ! Assume C1 molecule, integer :: B_i, B_j integer::dimn, sign - complex*16, allocatable :: E_type3DM(:, :, :, :, :, :), E_type4DM(:, :, :, :, :, :, :, :) + complex*16, allocatable :: E_type2DM(:, :, :, :), E_type3DM(:, :, :, :, :, :), E_type4DM(:, :, :, :, :, :, :, :) + complex*16, allocatable :: S_1DM(:, :), S_2DM(:, :, :, :) complex*16 :: S_type_4DM, S_type_3DM, S_type_2DM, S_type_1DM character(50) :: S_type_dm4_filename, S_type_dm3_filename, S_type_dm2_filename, S_type_dm1_filename - character(100) :: E_type_4dm_filename + character(100) :: E_type_4dm_filename, E_type_2dm_filename integer :: dmrg_dm4_unit_num = 500, dmrg_dm3_unit_num = 550, dmrg_dm2_unit_num = 600, dmrg_dm1_unit_num = 650 - integer :: E_type_4dm_unit_num = 700 + integer :: E_type_2dm_unit_num = 700, E_type_3dm_unit_num = 750, E_type_4dm_unit_num = 800 - integer ::i1, i2, i3, i4, i5, i6, i7, i8 - integer ::t_i1, t_i2, t_i3, t_i4, t_i5, t_i6, t_i7, t_i8 + integer :: i1, i2, i3, i4, i5, i6, i7, i8 + integer :: t_i1, t_i2, t_i3, t_i4, t_i5, t_i6, t_i7, t_i8 + integer :: j1, j2, j3, j4 logical :: sym_check character(30) ::chr_isym @@ -683,6 +685,173 @@ SUBROUTINE bAmat_interface_for_debug(sc, bc, isym, indsym) ! Assume C1 molecule, bc(:, :) = 0.0d+00 + ! ========================================================================== + ! E型2DM <0|EijEkl|0> + ! ========================================================================== + + allocate (left_list(3, 2), right_list(3, 2)) + + allocate (E_type2DM(nact, nact, nact, nact)) + allocate (S_2DM(nact, nact, nact, nact)) + allocate (S_1DM(nact, nact)) + + ! initialization + + E_type2DM(:, :, :, :) = 0.0d+00 + S_2DM(:, :, :, :) = 0.0d+00 + S_1DM(:, :) = 0.0d+00 + + ! dmrgで出したDMのopen + + S_type_dm2_filename = "dmrg_2dm" + S_type_dm1_filename = "dmrg_1dm" + + call open_formatted_file(unit=dmrg_dm2_unit_num, file=trim(S_type_dm2_filename), status="old", optional_action="read") + call open_formatted_file(unit=dmrg_dm1_unit_num, file=trim(S_type_dm1_filename), status="old", optional_action="read") + + ! =!=!=!=!=!=!=!=!=!=!=!=!=!=!=!=!=!=!=!==!=!=!=!=!=!=!=!=!==!=!=!=!=!=!=!=!=!==!=!=!=!=!=!=!=!=!= + + if (rank == 0) print *, "start_add_S_type2dm_to_E_type2dm" + +! open & readS型2DMの繰り返しの中 + ! S型2DMのreadとE型2DMへの足しこみ + + if (rank == 0) print *, "read S_type2DM" + do + ! read (DM4_unit_num, iostat=iostat) i, j, k, l, dm2 + read (dmrg_dm2_unit_num, *, iostat=iostat) i1, i2, i3, i4, S_2DM(i1, i2, i3, i4) + + if (iostat < 0) then + print *, "IOSTAT < 0, End read a file. filename: ", trim(S_type_dm2_filename) + exit + else if (iostat > 0) then + print *, "ERROR: error occured while reading a file. filename: ", trim(S_type_dm2_filename) + print *, "Exit PROGRAM" + stop + + end if + + end do! read S_type_2DM + + if (rank == 0) print *, "read S_type1DM" + do + ! read (DM4_unit_num, iostat=iostat) i, j, dm1 + read (dmrg_dm1_unit_num, *, iostat=iostat) i1, i2, S_1DM(i1, i2) + + if (iostat < 0) then + print *, "IOSTAT < 0, End read a file. filename: ", trim(S_type_dm1_filename) + exit + else if (iostat > 0) then + print *, "ERROR: error occured while reading a file. filename: ", trim(S_type_dm1_filename) + print *, "Exit PROGRAM" + stop + + end if + + end do! read S_type_1DM + + ! =)=)=)=)=)=)=)=)=)=)=)=)=)=)=)=)=)=)=)=)=)=)=)=)=)=)=)=)=)=)=)=)=)=)=)=)=)=) + ! + ! E型2DMをS型1-2DMで計算する式 + ! + ! E_type2DM(i,j,k,l) = d_kj*S_type1DM(i,l) - S_type2DM(i,k,j,l) + ! + ! =)=)=)=)=)=)=)=)=)=)=)=)=)=)=)=)=)=)=)=)=)=)=)=)=)=)=)=)=)=)=)=)=)=)=)=)=)=) + + print *, 'from_caspt2_idx_to_dmrg' + print *, from_caspt2_idx_to_dmrg + + ! j:caspt2での番号(active内の番号) + do j1 = 1, nact + do j2 = 1, nact + do j3 = 1, nact + do j4 = 1, nact + + ! i:S型DMにアクセスするためにdmrgの番号(active内の番号)へ変換 + i1 = from_caspt2_idx_to_dmrg(j1) + i2 = from_caspt2_idx_to_dmrg(j2) + i3 = from_caspt2_idx_to_dmrg(j3) + i4 = from_caspt2_idx_to_dmrg(j4) + + if (i1 == 0 .or. i2 == 0 .or. i3 == 0 .or. i4 == 0) cycle + + ! S型1DMの足しこみ d_kj*S_type1DM(i,l) + if (i2 == i3) then + + E_type2DM(j1, j2, j3, j4) = E_type2DM(j1, j2, j3, j4) + S_1DM(i1, i4) + + print *, 'j1, j2, j3, j4,E_type2DM(j1, j2, j3, j4)' + print *, j1, j2, j3, j4, E_type2DM(j1, j2, j3, j4) + + end if + + ! S型2DMの足しこみ - S_type2DM(i,k,j,l) + + ! S型2DMの左の2つのindexの置換 + call create_permutation_list(i1, i2, left_list) + ! S型2DMの右の2つのindexの置換 + call create_permutation_list(i3, i4, right_list) + + ! 複製したすべてのパターンのS型2DMを考え、それぞれE型2DMに必要であれば足しこむ + do left_index = 1, 2 + do right_index = 1, 2 + + sign = left_list(3, left_index)*right_list(3, right_index) + + E_type2DM(j1, j2, j3, j4) = E_type2DM(j1, j2, j3, j4) & + - S_2DM(i1, i2, i3, i4)*sign + + print *, 'after - S_2DM j1, j2, j3, j4,E_type2DM(j1, j2, j3, j4)' + print *, j1, j2, j3, j4, E_type2DM(j1, j2, j3, j4) + + end do!right_index + end do!left_index + + end do! j4 + end do! j3 + end do! j2 + end do! j1 + + deallocate (left_list, right_list) + + ! close S_type_1-2DM + close (dmrg_dm2_unit_num) + close (dmrg_dm1_unit_num) + + ! Write E_type2DM created from S_type1-2DM + + write (chr_isym, '(I2)') isym + + E_type_2dm_filename = "E_type_2dm_created_from_S_type1-2DM_"//trim(adjustl(chr_isym)) + + call open_formatted_file(unit=E_type_2dm_unit_num, file=trim(E_type_2dm_filename), status="replace", optional_action="write") + + if (rank == 0) print *, "Write_E_type2DM created_from_S_type1-2DM" + + ! caspt2のindex + do j1 = 1, nact + do j2 = 1, nact + do j3 = 1, nact + do j4 = 1, nact + + if (E_type2DM(j1, j2, j3, j4) /= 0.0d+00) then + write (E_type_2dm_unit_num, '(2e20.10,4I4)') E_type2DM(j1, j2, j3, j4), j1, j2, j3, j4 + ! E_type_2dm created from CASCI + ! write (dmat_unit_nonzero, '(e20.10,4I4)') & + ! dr_t, i_t, j_t, k_t, l_t + end if + + End do + End do + End do + End do + + ! E型2DMのclose + close (E_type_2dm_unit_num) + if (rank == 0) print *, "Close_E_type2DM_file_created_from_S_type1-2DM" + + ! =!=!=!=!=!=!=!=!=!=!=!=!=!=!=!=!=!=!=!==!=!=!=!=!=!=!=!=!==!=!=!=!=!=!=!=!=!==!=!=!=!=!=!=!=!=!= + ! ========================================================================== ! E型4DM <0|EzyEtxEuvEww|0> ! ========================================================================== @@ -731,6 +900,7 @@ SUBROUTINE bAmat_interface_for_debug(sc, bc, isym, indsym) ! Assume C1 molecule, end if + ! caspt2のactiveの番号へ変換 t_i1 = from_dmrg_idx_to_caspt2(i1) - ninact t_i2 = from_dmrg_idx_to_caspt2(i2) - ninact t_i3 = from_dmrg_idx_to_caspt2(i3) - ninact diff --git a/src/r4dcasci_co.f90 b/src/r4dcasci_co.f90 index 38f51b45..b8914a3d 100644 --- a/src/r4dcasci_co.f90 +++ b/src/r4dcasci_co.f90 @@ -318,10 +318,15 @@ PROGRAM r4dcasci_co ! DO CASCI CALC IN THIS PROGRAM! write (dmat_unit) i_t, j_t, k_t, l_t, DM2 ! end if - if (dr_t /= 0.0d+00) write (dmat_unit_nonzero, '(e20.10,4I4)') & - dr_t, i_t, j_t, k_t, l_t + ! if (dr_t /= 0.0d+00) write (dmat_unit_nonzero, '(e20.10,4I4)') & + ! dr_t, i_t, j_t, k_t, l_t ! write (dmat_unit) dr_t 12_09オリジナル + + if (DM2 /= 0.0d+00) then + write (dmat_unit_nonzero, '(2e20.10,4I4)') DM2, i_t, j_t, k_t, l_t + end if + End do End do End do