|
| 1 | +! This file is part of mctc-lib. |
| 2 | +! |
| 3 | +! Licensed under the Apache License, Version 2.0 (the "License"); |
| 4 | +! you may not use this file except in compliance with the License. |
| 5 | +! You may obtain a copy of the License at |
| 6 | +! |
| 7 | +! http://www.apache.org/licenses/LICENSE-2.0 |
| 8 | +! |
| 9 | +! Unless required by applicable law or agreed to in writing, software |
| 10 | +! distributed under the License is distributed on an "AS IS" BASIS, |
| 11 | +! WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. |
| 12 | +! See the License for the specific language governing permissions and |
| 13 | +! limitations under the License. |
| 14 | + |
| 15 | +module mctc_io_read_qchem |
| 16 | + use mctc_env_accuracy, only : wp |
| 17 | + use mctc_env_error, only : error_type |
| 18 | + use mctc_io_convert, only : aatoau |
| 19 | + use mctc_io_resize, only : resize |
| 20 | + use mctc_io_symbols, only : symbol_length, to_number, to_symbol |
| 21 | + use mctc_io_structure, only : structure_type, new |
| 22 | + use mctc_io_utils, only : next_line, token_type, next_token, io_error, filename, & |
| 23 | + read_next_token, read_token |
| 24 | + implicit none |
| 25 | + private |
| 26 | + |
| 27 | + public :: read_qchem |
| 28 | + |
| 29 | + integer, parameter :: initial_size = 64 |
| 30 | + |
| 31 | +contains |
| 32 | + |
| 33 | +subroutine read_qchem(mol, unit, error) |
| 34 | + |
| 35 | + !> Instance of the molecular structure data |
| 36 | + type(structure_type), intent(out) :: mol |
| 37 | + |
| 38 | + !> File handle |
| 39 | + integer, intent(in) :: unit |
| 40 | + |
| 41 | + !> Error handling |
| 42 | + type(error_type), allocatable, intent(out) :: error |
| 43 | + |
| 44 | + integer :: stat, pos, lnum, izp, iat |
| 45 | + integer :: charge, multiplicity |
| 46 | + type(token_type) :: token |
| 47 | + character(len=:), allocatable :: line |
| 48 | + real(wp) :: x, y, z |
| 49 | + character(len=symbol_length), allocatable :: sym(:) |
| 50 | + real(wp), allocatable :: xyz(:, :), abc(:, :), lattice(:, :) |
| 51 | + logical :: is_frac, periodic(3) |
| 52 | + |
| 53 | + iat = 0 |
| 54 | + lnum = 0 |
| 55 | + stat = 0 |
| 56 | + |
| 57 | + do while(stat == 0) |
| 58 | + call next_line(unit, line, pos, lnum, stat) |
| 59 | + if (stat /= 0) exit |
| 60 | + |
| 61 | + call next_token(line, pos, token) |
| 62 | + if (token%first > len(line)) cycle |
| 63 | + if (to_lower(line(token%first:token%last)) == '$molecule') exit |
| 64 | + end do |
| 65 | + |
| 66 | + if (stat /= 0) then |
| 67 | + call io_error(error, "No atoms found", & |
| 68 | + & line, token_type(0, 0), filename(unit), lnum+1, "expected molecule block") |
| 69 | + return |
| 70 | + end if |
| 71 | + |
| 72 | + call next_line(unit, line, pos, lnum, stat) |
| 73 | + if (stat == 0) & |
| 74 | + call read_next_token(line, pos, token, charge, stat) |
| 75 | + if (stat == 0) & |
| 76 | + call read_next_token(line, pos, token, multiplicity, stat) |
| 77 | + if (stat /= 0) then |
| 78 | + call io_error(error, "Failed to read charge and multiplicity", & |
| 79 | + & line, token, filename(unit), lnum, "expected integer value") |
| 80 | + return |
| 81 | + end if |
| 82 | + |
| 83 | + allocate(sym(initial_size), source=repeat(' ', symbol_length)) |
| 84 | + allocate(xyz(3, initial_size), source=0.0_wp) |
| 85 | + |
| 86 | + do while(stat == 0) |
| 87 | + call next_line(unit, line, pos, lnum, stat) |
| 88 | + if (stat /= 0) exit |
| 89 | + |
| 90 | + call next_token(line, pos, token) |
| 91 | + if (to_lower(line(token%first:token%last)) == '$end') exit |
| 92 | + |
| 93 | + if (iat >= size(sym)) call resize(sym) |
| 94 | + if (iat >= size(xyz, 2)) call resize(xyz) |
| 95 | + iat = iat + 1 |
| 96 | + |
| 97 | + token%last = min(token%last, token%first + symbol_length - 1) |
| 98 | + sym(iat) = line(token%first:token%last) |
| 99 | + if (to_number(sym(iat)) == 0) then |
| 100 | + call read_token(line, token, izp, stat) |
| 101 | + sym(iat) = to_symbol(izp) |
| 102 | + end if |
| 103 | + if (stat /= 0) then |
| 104 | + call io_error(error, "Cannot map symbol to atomic number", & |
| 105 | + & line, token, filename(unit), lnum, "unknown element") |
| 106 | + return |
| 107 | + end if |
| 108 | + |
| 109 | + call read_next_token(line, pos, token, x, stat) |
| 110 | + if (stat == 0) & |
| 111 | + call read_next_token(line, pos, token, y, stat) |
| 112 | + if (stat == 0) & |
| 113 | + call read_next_token(line, pos, token, z, stat) |
| 114 | + if (stat /= 0) then |
| 115 | + call io_error(error, "Cannot read coordinates", & |
| 116 | + & line, token, filename(unit), lnum, "expected real value") |
| 117 | + return |
| 118 | + end if |
| 119 | + |
| 120 | + xyz(:, iat) = [x, y, z] * aatoau |
| 121 | + end do |
| 122 | + |
| 123 | + if (stat /= 0) then |
| 124 | + call io_error(error, "Failed to read molecule block", & |
| 125 | + & line, token_type(0, 0), filename(unit), lnum, "unexpected end of input") |
| 126 | + return |
| 127 | + end if |
| 128 | + |
| 129 | + call new(mol, sym(:iat), xyz, charge=real(charge, wp), uhf=multiplicity-1) |
| 130 | +end subroutine read_qchem |
| 131 | + |
| 132 | +!> Convert input string to lowercase |
| 133 | +elemental function to_lower(str) result(lcstr) |
| 134 | + |
| 135 | + !> Input string |
| 136 | + character(len=*), intent(in) :: str |
| 137 | + |
| 138 | + !> Lowercase version of string |
| 139 | + character(len=len(str)):: lcstr |
| 140 | + |
| 141 | + integer :: ilen, iquote, i, iav, iqc |
| 142 | + integer, parameter :: offset = iachar('A') - iachar('a') |
| 143 | + |
| 144 | + ilen = len(str) |
| 145 | + iquote = 0 |
| 146 | + lcstr = str |
| 147 | + |
| 148 | + do i = 1, ilen |
| 149 | + iav = iachar(str(i:i)) |
| 150 | + if (iquote == 0 .and. (iav == 34 .or.iav == 39)) then |
| 151 | + iquote = 1 |
| 152 | + iqc = iav |
| 153 | + cycle |
| 154 | + end if |
| 155 | + if (iquote == 1 .and. iav==iqc) then |
| 156 | + iquote=0 |
| 157 | + cycle |
| 158 | + end if |
| 159 | + if (iquote == 1) cycle |
| 160 | + if (iav >= iachar('A') .and. iav <= iachar('Z')) then |
| 161 | + lcstr(i:i) = achar(iav - offset) |
| 162 | + else |
| 163 | + lcstr(i:i) = str(i:i) |
| 164 | + end if |
| 165 | + end do |
| 166 | + |
| 167 | +end function to_lower |
| 168 | + |
| 169 | +end module mctc_io_read_qchem |
0 commit comments