Skip to content

Commit

Permalink
calc E_type2DM from S_type1-2DM
Browse files Browse the repository at this point in the history
  • Loading branch information
yasuto-masuda committed Jul 9, 2023
1 parent f2b57e0 commit 2bbbfaa
Show file tree
Hide file tree
Showing 2 changed files with 182 additions and 7 deletions.
180 changes: 175 additions & 5 deletions src/DMRG_interface.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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>
! ==========================================================================
Expand Down Expand Up @@ -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
Expand Down
9 changes: 7 additions & 2 deletions src/r4dcasci_co.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down

0 comments on commit 2bbbfaa

Please sign in to comment.