E type 2DM is successfully created from S type DM

src/CMakeLists.txt

@@ -35,6 +35,8 @@ add_executable(r4dcascicoexe
 	tramo_ty.f90
 	uramda_s_half.f90
 	DMRG_interface.f90
+	module_error.f90
+	module_index_utils.f90
 )
 
 add_executable(r4dcaspt2ocoexe
@@ -67,6 +69,8 @@ add_executable(r4dcaspt2ocoexe
 	solvF_dmrg.f90
 	solvG_dmrg.f90
 	DMRG_interface.f90
+	module_error.f90
+	module_index_utils.f90
 
 
 

src/DM4_interface.f90

@@ -0,0 +1,206 @@
+module DM4_interface
+
+    implicit none
+
+    interface create_permutation_list
+
+        module procedure create_permutation_list_DM2, create_permutation_list_DM3, create_permutation_list_DM4
+
+    end interface create_permutation_list
+
+contains
+
+    subroutine create_permutation_list_DM2(i, j, list)
+        implicit none
+
+        integer, intent(inout) :: i, j
+
+        integer, intent(out) :: list(3, 2)
+
+        list(:, 1) = (/i, j, 1/)
+
+        list(:, 2) = (/j, i, -1/)
+
+    end subroutine create_permutation_list_DM2
+
+    subroutine create_permutation_list_DM3(i, j, k, list)
+        implicit none
+
+        integer, intent(inout) :: i, j, k
+
+        integer, intent(out) :: list(4, 6)
+
+        list(:, 1) = (/i, j, k, 1/)
+
+        list(:, 2) = (/i, k, j, -1/)
+
+        list(:, 3) = (/j, i, k, -1/)
+
+        list(:, 4) = (/j, k, i, 1/)
+
+        list(:, 5) = (/k, i, j, 1/)
+
+        list(:, 6) = (/k, j, i, -1/)
+
+    end subroutine create_permutation_list_DM3
+
+    subroutine create_permutation_list_DM4(i, j, k, l, list)
+        implicit none
+
+        integer, intent(inout) :: i, j, k, l
+
+        integer, intent(out) :: list(5, 24)
+        list(:, 1) = (/i, j, k, l, 1/)
+        list(:, 2) = (/j, i, k, l, -1/)
+        list(:, 3) = (/k, i, j, l, 1/)
+        list(:, 4) = (/i, k, j, l, -1/)
+        list(:, 5) = (/j, k, i, l, 1/)
+        list(:, 6) = (/k, j, i, l, -1/)
+        list(:, 7) = (/k, j, l, i, 1/)
+        list(:, 8) = (/j, k, l, i, -1/)
+        list(:, 9) = (/l, k, j, i, 1/)
+        list(:, 10) = (/k, l, j, i, -1/)
+        list(:, 11) = (/j, l, k, i, 1/)
+        list(:, 12) = (/l, j, k, i, -1/)
+        list(:, 13) = (/l, i, k, j, 1/)
+        list(:, 14) = (/i, l, k, j, -1/)
+        list(:, 15) = (/k, l, i, j, 1/)
+        list(:, 16) = (/l, k, i, j, -1/)
+        list(:, 17) = (/i, k, l, j, 1/)
+        list(:, 18) = (/k, i, l, j, -1/)
+        list(:, 19) = (/j, i, l, k, 1/)
+        list(:, 20) = (/i, j, l, k, -1/)
+        list(:, 21) = (/l, j, i, k, 1/)
+        list(:, 22) = (/j, l, i, k, -1/)
+        list(:, 23) = (/i, l, j, k, 1/)
+        list(:, 24) = (/l, i, j, k, -1/)
+
+    end subroutine create_permutation_list_DM4
+
+! +=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=
+! +=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=
+
+    SUBROUTINE bAmat_interface(sc, bc) ! Assume C1 molecule, overlap matrix B in space A
+        !
+        !  B(xyz,tuv) = Siguma_w eps(w){-<0|EzyEtxEuvEww|0> + d(tx)<0|EzyEuvEww|0>}
+        !
+        !               + S(xyz,tuv)(eps(u)+eps(t)-eps(v))
+        !
+        ! +=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=
+        ! +=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=
+
+        use four_caspt2_module
+        use module_file_manager
+
+        Implicit NONE
+#ifdef HAVE_MPI
+        include 'mpif.h'
+#endif
+        integer :: iu, iv
+        integer :: jt, ju, jv, jw
+        integer :: iw1, iw2, iostat
+
+        complex*16, intent(in)  :: sc(dimn, dimn)
+        complex*16, intent(out) :: bc(dimn, dimn)
+        integer, intent(in) ::isym
+        integer :: counter_use_3dm_Bamat
+
+        real*8               :: e
+        complex*16           :: dm3, dm4
+        character(50)        :: chr_isym, chr_rank, dmat4A_filename, dmat3A_B_filename
+
+        integer :: DM3_unit_num = 300, DM4_unit_num = 400
+        bc(:, :) = 0.0d+00
+
+        ! Get a 4DM file name and open it
+
+        dmat4A_filename = "required_S_type_dmat4A"
+        call open_unformatted_file(unit=DM4_unit_num, file=trim(dmat4A_filename), status="old", optional_action="read")
+
+        ! Do!iostatを使ってファイルのすべてを読む
+        !     read ijklmnop, S型4DM1個
+
+        !     Do 複製type = 1,N
+        !  !(N個複製する必要があるとする、4次のDMならこの前の計算があっていれば24個)
+        !  !(事前に配列で、type=1はindexとしてijklmnopみたいなのを作っておく。)
+        !  !(ここから処理部)
+
+        !                              !(dirac_orderからenergy_orderに変換)
+        !                              call convert dirac_order to energy_order(intent(in)D_i..,intent(out)E_i..)
+
+        !                              !（インデックスをxyztuv...に変換）
+        !                            call convert to A_subspace_index&
+        !              (intent(in)ijklmn... intent(out)xyztuv...)
+
+        !                               !（部分空間AorCの判定条件を満たすか確かめる）
+        !                             call sym_check_ABmat(intent (out)sym_check,intent(in)xyz...)
+
+        !                              if (sym_check==.true.)then
+
+        if (rank == 0) print *, "solv_A_dmrg_bAmat_read 4DM file"
+        do
+            read (DM4_unit_num, iostat=iostat) i, j, k, l, m, n, o, p, dm4
+            if (iostat < 0) then
+                print *, "IOSTAT < 0, End read a file. filename: ", trim(dmat4A_filename)
+                exit
+            else if (iostat > 0) then
+                print *, "ERROR: error occured while reading a file. filename: ", trim(dmat4A_filename)
+                print *, "Exit PROGRAM"
+                stop
+
+            end if
+
+            !(dirac_orderからenergy_orderに変換)
+            call convert_index_dirac_order_to_energy_order(i, j, k, l, m, n, o, p)
+
+            call create_permutation_list(i, j, k, l, left_list)
+            call create_permutation_list(m, n, o, p, right_list)
+
+            do left_idx = 1, 24 !左側のS型４DMのインデックスパターンを指定するもの
+                do right_idx = 1, 24 !右側のS型４DMのインデックスパターンを指定するもの
+
+                    !（インデックスをxyztuv...に変換）
+            call convert_to_A_subspace_index(left_list(1:4, left_idx), right_list(1:4, right_idx), iz, iy, it, ix, iu, iv, iw1, iw2)
+
+                    !（部分空間AorCの判定条件を満たすか確かめる）
+                    call sym_check_ABmat(intent(out) sym_check, intent(in) xyz...)
+
+                    if (sym_check) then
+
+                        B_i = indsym_reverse(ix, iy, iz)
+                        B_j = indsym_reverse(it, iu, iv)
+
+                        if (B_i > B_j) print *, "i is larger than B_j 4dm", B_i, B_j
+                        if (B_i == 0 .or. B_j == 0 .or. B_i > B_j) cycle
+
+                        bc(B_i, B_j) = bc(B_i, B_j) - dm4*eps(iw2 + ninact)
+
+                    end if
+
+                end do
+
+            end do
+        end do
+
+    end SUBROUTINE bAmat_interface
+
+    subroutine convert_to_A_subspace_index(left_list, right_list, iz, iy, it, ix, iu, iv, iw1, iw2)
+
+        implicit NONE
+
+        integer, intent(out) ::iz, iy, it, ix, iu, iv, iw1, iw2
+        integer, intent(in) :: left_list(4), right_list(4)
+
+        iz = left_list(1)
+        iy = left_list(2)
+        it = left_list(3)
+        ix = left_list(4)
+        iu = right_list(1)
+        iv = right_list(2)
+        iw1 = right_list(3)
+        iw2 = right_list(4)
+
+    end subroutine convert_to_A_subspace_index
+
+end module DM4_interface
+

src/DMRG_interface.f90

@@ -2,6 +2,7 @@ module DMRG_interface
 
     use four_caspt2_module
     use module_file_manager
+    use module_index_utils
 
     implicit none
 
@@ -721,6 +722,17 @@ SUBROUTINE bAmat_interface_for_debug(sc, bc, isym, indsym) ! Assume C1 molecule,
             ! 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)
 
+            call create_permutation_list(i1, i2, left_list)
+            call create_permutation_list(i3, i4, right_list)
+
+            do left_index = 1, 2
+                do right_index = 1, 2
+                    sign = left_list(3, left_index)*right_list(3, right_index)
+                 S_2DM(left_list(1, left_index), left_list(2, left_index), right_list(1, right_index), right_list(2, right_index)) &
+                        = S_2DM(i1, i2, i3, i4)*sign
+                end do
+            end do
+
             if (iostat < 0) then
                 print *, "IOSTAT < 0, End read a file. filename: ", trim(S_type_dm2_filename)
                 exit
@@ -738,6 +750,8 @@ SUBROUTINE bAmat_interface_for_debug(sc, bc, isym, indsym) ! Assume C1 molecule,
             ! read (DM4_unit_num, iostat=iostat) i, j, dm1
             read (dmrg_dm1_unit_num, *, iostat=iostat) i1, i2, S_1DM(i1, i2)
 
+            S_1DM(i2, i1) = -S_1DM(i1, i2)
+
             if (iostat < 0) then
                 print *, "IOSTAT < 0, End read a file. filename: ", trim(S_type_dm1_filename)
                 exit
@@ -768,15 +782,16 @@ SUBROUTINE bAmat_interface_for_debug(sc, bc, isym, indsym) ! Assume C1 molecule,
                     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)
+                        i1 = from_caspt2_idx_to_dmrg(convert_active_to_global_idx(j1))
+                        i2 = from_caspt2_idx_to_dmrg(convert_active_to_global_idx(j2))
+                        i3 = from_caspt2_idx_to_dmrg(convert_active_to_global_idx(j3))
+                        i4 = from_caspt2_idx_to_dmrg(convert_active_to_global_idx(j4))
 
-                        if (i1 == 0 .or. i2 == 0 .or. i3 == 0 .or. i4 == 0) cycle
+                        ! 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
+                            ! if (j2 == j3) thenとするべきかも？
 
                             E_type2DM(j1, j2, j3, j4) = E_type2DM(j1, j2, j3, j4) + S_1DM(i1, i4)
 
@@ -2697,11 +2712,11 @@ subroutine fcidump_integral_transform
 
         implicit none
 ! FCIDUMP is 1,2-integral file to calculate DMRG(Qcmaquis)
-! This subrutine read FCIDUM data and create the list which convert DMRG spinor indices into CASPT2 spinor indices.
+! This subrutine read FCIDUMP data and create the list which convert DMRG spinor indices into CASPT2 spinor indices.
 ! This list is named from_dmrg_idx_to_caspt2(dmrg_idx)
 
 ! FCIDUMP is 1,2-integral file to calculate DMRG(Qcmaquis)
-! This subrutine read FCIDUM data and create the list which convert DMRG spinor indices into CASPT2 spinor indices.
+! This subrutine read FCIDUMP data and create the list which convert DMRG spinor indices into CASPT2 spinor indices.
 ! This list is named from_dmrg_idx_to_caspt2(dmrg_idx)
 
         integer :: file_unit, file_unit_trans
@@ -2723,14 +2738,17 @@ subroutine fcidump_integral_transform
 
         character(100) :: header(4), form1
         integer :: start_idx, end_idx
+        character(:), allocatable  :: FCIDUMP_filename
+
+        FCIDUMP_filename = "FCIDUMP.save"
 
 ! read header
 
-        call open_formatted_file(file_unit, 'FCIDUMP', 'old', 'read')
+        call open_formatted_file(file_unit, FCIDUMP_filename, 'old', 'read')
 
         do idx = 1, 4
             read (file_unit, '(A)', iostat=iostat) header(idx)
-            call check_iostat(iostat, 'FCIDUMP', end_of_file)
+            call check_iostat(iostat, FCIDUMP_filename, end_of_file)
 
         end do
 
@@ -2745,7 +2763,7 @@ subroutine fcidump_integral_transform
         do
             read (file_unit, *, iostat=iostat) integral_r, i, j, k, l
 
-            call check_iostat(iostat, 'FCIDUMP', end_of_file)
+            call check_iostat(iostat, FCIDUMP_filename, end_of_file)
 
             if (end_of_file) exit
 
@@ -2833,7 +2851,6 @@ subroutine fcidump_integral_transform
         close (file_unit_trans)
 
     contains
-
         subroutine create_from_dmrg_idx_to_caspt2
             implicit none
 
@@ -2915,6 +2932,7 @@ subroutine create_from_dmrg_idx_to_caspt2
             end do
 
         end subroutine create_from_dmrg_idx_to_caspt2
+
         subroutine transform_1_integral
             implicit none
 
@@ -2970,8 +2988,6 @@ subroutine transform_2_integral
 ! use module_file_manager
 ! use four_caspt2_module
 
-# 1260
-
             integer  :: spi, spj, spk, spl
             character(50)   :: fname, chr_rank
             logical :: is_end_of_file