File: main.f90

package info (click to toggle)
mctc-lib 0.3.2-1
  • links: PTS, VCS
  • area: main
  • in suites: forky, sid, trixie
  • size: 1,088 kB
  • sloc: f90: 12,894; python: 23; ansic: 15; makefile: 2
file content (346 lines) | stat: -rw-r--r-- 11,057 bytes parent folder | download | duplicates (2)
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
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
! This file is part of mctc-lib.
!
! Licensed under the Apache License, Version 2.0 (the "License");
! you may not use this file except in compliance with the License.
! You may obtain a copy of the License at
!
!     http://www.apache.org/licenses/LICENSE-2.0
!
! Unless required by applicable law or agreed to in writing, software
! distributed under the License is distributed on an "AS IS" BASIS,
! WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
! See the License for the specific language governing permissions and
! limitations under the License.

!> Example application using tool chain library.
!>
!> This program uses the [[read_structure]] and [[write_structure]] procedures
!> to implement a structure converter.
!> Usually, the input structure can be inferred by the name of the input file.
!> To allow formats with non-standard extensions (because most geometry formats
!> are not really standardized) additional hints can be passed by the command
!> line to determine the read/write formats.
!>
!> To add support for piping standard input and standard output reading and
!> writing from units is combined with the additional format hints.
!>
!> Additional filters or modifications can also be implemented in an intermediary
!> step, this program implements an element symbol normalization. Other filters
!> like folding back to central cells or removing lattice vector could be added
!> in a similar manner.
program main
   use, intrinsic :: iso_fortran_env, only : output_unit, error_unit, input_unit
   use mctc_env, only : error_type, fatal_error, get_argument, wp
   use mctc_io, only : structure_type, read_structure, write_structure, &
      & filetype, get_filetype, to_symbol
   use mctc_version, only : get_mctc_version
   implicit none
   character(len=*), parameter :: prog_name = "mctc-convert"

   character(len=:), allocatable :: input, output, template, filename
   integer, allocatable :: input_format, output_format, template_format
   type(structure_type) :: mol
   type(structure_type), allocatable :: mol_template
   type(error_type), allocatable :: error
   logical :: normalize, read_dot_files
   integer :: charge, unpaired

   call get_arguments(input, input_format, output, output_format, normalize, &
      & template, template_format, read_dot_files, error)
   if (allocated(error)) then
      write(error_unit, '(a)') error%message
      error stop
   end if

   if (allocated(template)) then
      allocate(mol_template)
      if (template == "-") then
         if (.not.allocated(template_format)) then
            template_format = merge(output_format, filetype%xyz, allocated(output_format))
         end if
         call read_structure(mol_template, input_unit, template_format, error)
      else
         call read_structure(mol_template, template, error, template_format)
      end if
      if (allocated(error)) then
         write(error_unit, '(a)') error%message
         error stop
      end if
   end if

   if (input == "-") then
      if (.not.allocated(input_format)) input_format = filetype%xyz
      call read_structure(mol, input_unit, input_format, error)
   else
      call read_structure(mol, input, error, input_format)

      if (read_dot_files) then
         charge = nint(mol%charge)
         if (.not.allocated(error)) then
            filename = join(dirname(input), ".CHRG")
            if (exists(filename)) call read_file(filename, charge, error)
         end if
         mol%charge = charge

         unpaired = mol%uhf
         if (.not.allocated(error)) then
            filename = join(dirname(input), ".UHF")
            if (exists(filename)) call read_file(filename, unpaired, error)
         end if
         mol%uhf = unpaired
      end if
   end if
   if (allocated(error)) then
      write(error_unit, '(a)') error%message
      error stop
   end if

   if (allocated(mol_template)) then
      if (mol%nat /= mol_template%nat) then
         write(error_unit, '(*(a, 1x))') &
            "Number of atoms missmatch in", template, "and", input
         error stop
      end if

      ! move_alloc can also move non-allocated objects
      call move_alloc(mol_template%lattice, mol%lattice)
      call move_alloc(mol_template%periodic, mol%periodic)
      call move_alloc(mol_template%bond, mol%bond)
      call move_alloc(mol_template%comment, mol%comment)
      call move_alloc(mol_template%pdb, mol%pdb)
      call move_alloc(mol_template%sdf, mol%sdf)
   end if

   if (normalize) then
      mol%sym = to_symbol(mol%num)
   end if

   if (output == "-") then
      if (.not.allocated(output_format)) output_format = filetype%xyz
      call write_structure(mol, output_unit, output_format, error)
   else
      call write_structure(mol, output, error, output_format)
   end if
   if (allocated(error)) then
      write(error_unit, '(a)') error%message
      error stop
   end if


contains


subroutine help(unit)
   integer, intent(in) :: unit

   write(unit, '(a, *(1x, a))') &
      "Usage: "//prog_name//" [options] <input> <output>"

   write(unit, '(a)') &
      "", &
      "Read structure from input file and writes it to output file.", &
      "The format is determined by the file extension or the format hint", &
      ""

   write(unit, '(2x, a, t25, a)') &
      "-i, --input <format>", "Hint for the format of the input file", &
      "-o, --output <format>", "Hint for the format of the output file", &
      "--normalize", "Normalize all element symbols to capitalized format", &
      "--template <file>", "File to use as template to fill in meta data", &
      "", "(useful to add back SDF or PDB annotions)", &
      "--template-format <format>", "", "", "Hint for the format of the template file", &
      "--ignore-dot-files", "Do not read charge and spin from .CHRG and .UHF files", &
      "--version", "Print program version and exit", &
      "--help", "Show this help message"

   write(unit, '(a)')

end subroutine help


subroutine version(unit)
   integer, intent(in) :: unit
   character(len=:), allocatable :: version_string

   call get_mctc_version(string=version_string)
   write(unit, '(a, *(1x, a))') &
      & prog_name, "version", version_string

end subroutine version


subroutine get_arguments(input, input_format, output, output_format, normalize, &
      & template, template_format, read_dot_files, error)

   !> Input file name
   character(len=:), allocatable :: input

   !> Input file format
   integer, allocatable, intent(out) :: input_format

   !> Output file name
   character(len=:), allocatable :: output

   !> Output file format
   integer, allocatable, intent(out) :: output_format

   !> Template file name
   character(len=:), allocatable :: template

   !> Template file format
   integer, allocatable, intent(out) :: template_format

   !> Normalize element symbols
   logical, intent(out) :: normalize

   !> Read information from .CHRG and .UHF files
   logical, intent(out) :: read_dot_files

   !> Error handling
   type(error_type), allocatable, intent(out) :: error

   integer :: iarg, narg
   character(len=:), allocatable :: arg

   normalize = .false.
   read_dot_files = .true.
   iarg = 0
   narg = command_argument_count()
   do while(iarg < narg)
      iarg = iarg + 1
      call get_argument(iarg, arg)
      select case(arg)
      case("--help")
         call help(output_unit)
         stop
      case("--version")
         call version(output_unit)
         stop
      case default
         if (.not.allocated(input)) then
            call move_alloc(arg, input)
            cycle
         end if
         if (.not.allocated(output)) then
            call move_alloc(arg, output)
            cycle
         end if
         call fatal_error(error, "Too many positional arguments present")
         exit
      case("-i", "--input")
         iarg = iarg + 1
         call get_argument(iarg, arg)
         if (.not.allocated(arg)) then
            call fatal_error(error, "Missing argument for input format")
            exit
         end if
         if (index(arg, ".") == 0) arg = "."//arg
         input_format = get_filetype(arg)
      case("-o", "--output")
         iarg = iarg + 1
         call get_argument(iarg, arg)
         if (.not.allocated(arg)) then
            call fatal_error(error, "Missing argument for output format")
            exit
         end if
         output_format = get_filetype("."//arg)
      case("--normalize")
         normalize = .true.
      case("--template")
         iarg = iarg + 1
         call get_argument(iarg, template)
         if (.not.allocated(template)) then
            call fatal_error(error, "Missing argument for template file")
            exit
         end if
      case("--template-format")
         iarg = iarg + 1
         call get_argument(iarg, arg)
         if (.not.allocated(arg)) then
            call fatal_error(error, "Missing argument for template format")
            exit
         end if
         template_format = get_filetype("."//arg)
      case("--ignore-dot-files")
         read_dot_files = .false.
      end select
   end do

   if (.not.(allocated(input).and.(allocated(output)))) then
      if (.not.allocated(error)) then
         call help(output_unit)
         error stop
      end if
   end if

end subroutine get_arguments


!> Extract dirname from path
function dirname(filename)
   character(len=*), intent(in) :: filename
   character(len=:), allocatable :: dirname

   dirname = filename(1:scan(filename, "/\", back=.true.))
   if (len_trim(dirname) == 0) dirname = "."
end function dirname


!> Construct path by joining strings with os file separator
function join(a1, a2) result(path)
   use mctc_env_system, only : is_windows
   character(len=*), intent(in) :: a1, a2
   character(len=:), allocatable :: path
   character :: filesep

   if (is_windows()) then
      filesep = '\'
   else
      filesep = '/'
   end if

   path = a1 // filesep // a2
end function join


!> test if pathname already exists
function exists(filename)
    character(len=*), intent(in) :: filename
    logical :: exists
    inquire(file=filename, exist=exists)
end function exists


subroutine read_file(filename, val, error)
   use mctc_io_utils, only : next_line, read_next_token, io_error, token_type
   character(len=*), intent(in) :: filename
   integer, intent(out) :: val
   type(error_type), allocatable, intent(out) :: error

   integer :: io, stat, lnum, pos
   type(token_type) :: token
   character(len=:), allocatable :: line

   lnum = 0

   open(file=filename, newunit=io, status='old', iostat=stat)
   if (stat /= 0) then
      call fatal_error(error, "Error: Could not open file '"//filename//"'")
      return
   end if

   call next_line(io, line, pos, lnum, stat)
   if (stat == 0) &
      call read_next_token(line, pos, token, val, stat)
   if (stat /= 0) then
      call io_error(error, "Cannot read value from file", line, token, &
         filename, lnum, "expected integer value")
      return
   end if

   close(io, iostat=stat)

end subroutine read_file


end program main