-
Notifications
You must be signed in to change notification settings - Fork 0
/
convert_mdcint.f90
240 lines (222 loc) · 8.97 KB
/
convert_mdcint.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Program convert_mdcint
! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Implicit None
Character :: datex*10, timex*8
character(50) :: mdcint_filename, mdcint_output_filename, argv
integer :: nkr, nmo, mdcint_unit = 10, file_output_unit = 11
integer :: i, i0, inz, nz, ikr, jkr, stat, mdcint_idx, split_num
integer, allocatable :: indk(:), indl(:), kr(:)
double precision, allocatable :: rklr(:), rkli(:)
logical :: realonly, exists, combine = .false., split = .false.
call get_command_argument(1, argv)
if (trim(argv) == "--combine") then
print *, "Combine all MDCINT files into one file (--combine option is detected)"
combine = .true.
else if (trim(argv) == "--split") then
print *, "Split MDCINT files into multiple files (--split option is detected)"
split = .true.
end if
if (split) then
call get_command_argument(2, argv)
read (argv, *) split_num
if (split_num <= 0) then
print *, "Error: split_num <= 0"
call exit(1)
end if
end if
inquire (file="MDCINT", exist=exists)
if (.not. exists) then
print *, "Error: MDCINT file does not exist"
call exit(1)
end if
open (mdcint_unit, file="MDCINT", form="unformatted", status="old")
read (mdcint_unit) datex, timex, nkr
rewind (mdcint_unit)
nmo = 2*nkr
allocate (kr(-nkr:nkr))
kr(:) = 0
read (mdcint_unit) datex, timex, nkr, (kr(i0), kr(-1*i0), i0=1, nkr)
allocate (indk(nmo**2))
allocate (indl(nmo**2))
allocate (rklr(nmo**2))
allocate (rkli(nmo**2))
realonly = is_realonly(mdcint_unit)
close(mdcint_unit)
if (split) then
call split_mdcint
else
call dump_formatted_mdcint
endif
deallocate (kr)
deallocate (indk)
deallocate (indl)
deallocate (rklr)
deallocate (rkli)
print *, "End program normally"
contains
logical function is_realonly(unit_num)
implicit none
integer, intent(in) :: unit_num
is_realonly = .false.
rewind (unit_num)
read (unit_num) ! Skip the first line
read (unit_num, iostat=stat) ikr, jkr, nz, (indk(inz), indl(inz), inz=1, nz), (rklr(inz), rkli(inz), inz=1, nz)
if (stat /= 0) then
is_realonly = .true.
else
is_realonly = .false.
end if
print *, "realonly = ", realonly
end function is_realonly
subroutine dump_formatted_mdcint
implicit none
if (combine) then
open (file_output_unit, file="formatted_MDCINT", form="formatted", status="replace")
write (file_output_unit, *) datex, timex, nkr, (kr(i0), kr(-1*i0), i0=1, nkr)
end if
mdcint_idx = 0
datex = ""
timex = ""
nkr = 0
do
call get_filename(mdcint_idx, mdcint_filename, mdcint_output_filename)
inquire (file=mdcint_filename, exist=exists)
if (.not. exists) then
exit
end if
open (mdcint_unit, file=trim(mdcint_filename), form="unformatted", status="old")
rewind (mdcint_unit)
read (mdcint_unit) datex, timex, nkr, (kr(i0), kr(-1*i0), i0=1, nkr)
if (.not. combine) then
open (file_output_unit, file=trim(mdcint_output_filename), form="formatted", status="replace")
write (file_output_unit, *) datex, timex, nkr, (kr(i0), kr(-1*i0), i0=1, nkr)
end if
do
if (realonly) then
read (mdcint_unit, iostat=stat) ikr, jkr, nz, (indk(inz), indl(inz), inz=1, nz), &
(rklr(inz), inz=1, nz)
else
read (mdcint_unit, iostat=stat) ikr, jkr, nz, (indk(inz), indl(inz), inz=1, nz), &
(rklr(inz), rkli(inz), inz=1, nz)
end if
if (ikr == 0) then
exit
end if
if (stat < 0) then
print *, "End of file: ", trim(mdcint_filename)
exit
else if (stat > 0) then
print *, "Error while reading file. iostat = ", stat
call exit(stat)
end if
if (realonly) then
do i = 1, nz
write (file_output_unit, '(4I6, E20.7)') ikr, jkr, indk(i), indl(i), rklr(i)
end do
else
do i = 1, nz
write (file_output_unit, '(4I6, 2E20.7)') ikr, jkr, indk(i), indl(i), rklr(i), rkli(i)
end do
end if
end do
close (mdcint_unit)
if (.not. combine) close (file_output_unit)
mdcint_idx = mdcint_idx + 1
end do
if (combine) close (file_output_unit)
end subroutine dump_formatted_mdcint
subroutine split_mdcint
implicit none
integer, allocatable :: mdcint_files_unit(:)
integer, parameter :: start_unit = 21
integer :: write_file_idx, cnt
mdcint_idx = 0
cnt = 0
write_file_idx = 0
datex = ""
timex = ""
nkr = 0
allocate (mdcint_files_unit(split_num))
! mkdir
call system("mkdir -p split")
do i = 0, split_num-1
mdcint_files_unit(i+1) = start_unit + i
call get_filename(i, mdcint_filename, mdcint_output_filename)
! Open split files into split directory
open (mdcint_files_unit(i+1), file="split/"//trim(mdcint_filename), form="unformatted", status="replace")
end do
do
call get_filename(mdcint_idx, mdcint_filename, mdcint_output_filename)
inquire (file=mdcint_filename, exist=exists)
if (.not. exists) then
exit
end if
open (mdcint_unit, file=trim(mdcint_filename), form="unformatted", status="old")
rewind (mdcint_unit)
read (mdcint_unit) datex, timex, nkr, (kr(i0), kr(-1*i0), i0=1, nkr)
do i = 1, split_num
write (mdcint_files_unit(i)) datex, timex, nkr, (kr(i0), kr(-1*i0), i0=1, nkr)
end do
do
if (realonly) then
read (mdcint_unit, iostat=stat) ikr, jkr, nz, (indk(inz), indl(inz), inz=1, nz), &
(rklr(inz), inz=1, nz)
else
read (mdcint_unit, iostat=stat) ikr, jkr, nz, (indk(inz), indl(inz), inz=1, nz), &
(rklr(inz), rkli(inz), inz=1, nz)
end if
if (ikr == 0) then
print *, "End of file: ", trim(mdcint_filename)
exit
end if
write_file_idx = mod(cnt, split_num)
if (realonly) then
write(mdcint_files_unit(write_file_idx+1)) ikr, jkr, nz, (indk(inz), indl(inz), inz=1, nz), &
(rklr(inz), inz=1, nz)
else
write(mdcint_files_unit(write_file_idx+1)) ikr, jkr, nz, (indk(inz), indl(inz), inz=1, nz), &
(rklr(inz), rkli(inz), inz=1, nz)
end if
cnt = cnt + 1
end do
close (mdcint_unit)
mdcint_idx = mdcint_idx + 1
end do
do i = 1, split_num
close (mdcint_files_unit(i))
end do
end subroutine split_mdcint
subroutine get_filename(idx, filename, output_filename)
! Get the filename of the MDCINT file
! idx = 0: filename = "MDCINT", output_filename = "formatted_MDCINT"
! idx > 0: filename = "MDCIN"//trim(adjustl(padding_x))//trim(adjustl(chr_idx)), output_filename = "formatted_MDCINT"//trim(adjustl(chr_idx))
! (e.g.) idx = 1: filename = "MDCINXXXX1", output_filename = "formatted_MDCINT1"
! idx = 100: filename = "MDCINXX100", output_filename = "formatted_MDCINT100"
integer, intent(in) :: idx
character(len=*), intent(out) :: filename, output_filename
character(50) :: chr_idx, padding_x
if (idx < 0) then
print *, "Error: idx < 0"
call exit(1)
else if (idx == 0) then
filename = "MDCINT"
output_filename = "formatted_MDCINT"
else
write (chr_idx, *) idx
if (idx < 10) then
padding_x = "XXXX"
else if (idx < 100) then
padding_x = "XXX"
else if (idx < 1000) then
padding_x = "XX"
else if (idx < 10000) then
padding_x = "X"
else if (idx < 100000) then
padding_x = ""
end if
filename = "MDCIN"//trim(adjustl(padding_x))//trim(adjustl(chr_idx))
output_filename = "formatted_MDCINT"//trim(adjustl(chr_idx))
end if
end subroutine get_filename
end Program convert_mdcint