Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add QChem z-matrix reader #73

Merged
merged 2 commits into from
Jan 30, 2025
Merged

Add QChem z-matrix reader #73

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
31f1b37
Add QChem z-matrix reader
awvwgk Jan 30, 2025
ddca908
Add extra test cases
awvwgk Jan 30, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
152 changes: 128 additions & 24 deletions src/mctc/io/read/qchem.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,7 @@
module mctc_io_read_qchem
use mctc_env_accuracy, only : wp
use mctc_env_error, only : error_type
use mctc_io_constants, only : pi
use mctc_io_convert, only : aatoau
use mctc_io_resize, only : resize
use mctc_io_symbols, only : symbol_length, to_number, to_symbol
Expand All @@ -41,14 +42,15 @@
!> Error handling
type(error_type), allocatable, intent(out) :: error

integer :: stat, pos, lnum, izp, iat
integer :: charge, multiplicity
integer :: stat, pos, lnum, izp, iat, iz, ij(3)
integer :: charge, multiplicity, zrepeat
type(token_type) :: token
character(len=:), allocatable :: line
real(wp) :: x, y, z
real(wp) :: x, y, z, zm(3)
character(len=symbol_length), allocatable :: sym(:)
real(wp), allocatable :: xyz(:, :), abc(:, :), lattice(:, :)
logical :: is_frac, periodic(3)
real(wp), parameter :: deg_to_rad = pi / 180.0_wp

iat = 0
lnum = 0
Expand Down Expand Up @@ -83,30 +85,95 @@
allocate(sym(initial_size), source=repeat(' ', symbol_length))
allocate(xyz(3, initial_size), source=0.0_wp)

do while(stat == 0)
call next_line(unit, line, pos, lnum, stat)
if (stat /= 0) exit
call next_line(unit, line, pos, lnum, stat)
if (stat /= 0) then
call io_error(error, "Failed to read molecule block", &
& line, token_type(0, 0), filename(unit), lnum, "unexpected end of input")
return

Check warning on line 92 in src/mctc/io/read/qchem.f90

View check run for this annotation

Codecov / codecov/patch

src/mctc/io/read/qchem.f90#L92 

Added line #L92 was not covered by tests
end if

call next_token(line, pos, token)
if (to_lower(line(token%first:token%last)) == '$end') exit
call next_token(line, pos, token)

if (iat >= size(sym)) call resize(sym)
if (iat >= size(xyz, 2)) call resize(xyz)
iat = iat + 1
iat = iat + 1

token%last = min(token%last, token%first + symbol_length - 1)
sym(iat) = line(token%first:token%last)
if (to_number(sym(iat)) == 0) then
call read_token(line, token, izp, stat)
sym(iat) = to_symbol(izp)
end if
if (stat /= 0) then
call io_error(error, "Cannot map symbol to atomic number", &
& line, token, filename(unit), lnum, "unknown element")
return
end if
token%last = min(token%last, token%first + symbol_length - 1)
sym(iat) = line(token%first:token%last)
if (to_number(sym(iat)) == 0) then
call read_token(line, token, izp, stat)
sym(iat) = to_symbol(izp)
end if
if (stat /= 0) then
call io_error(error, "Cannot map symbol to atomic number", &
& line, token, filename(unit), lnum, "unknown element")
return

Check warning on line 108 in src/mctc/io/read/qchem.f90

View check run for this annotation

Codecov / codecov/patch

src/mctc/io/read/qchem.f90#L108 

Added line #L108 was not covered by tests
end if

call read_next_token(line, pos, token, x, stat)
if (stat /= 0) then
stat = 0
xyz(:, iat) = [0, 0, 0] * aatoau
zrepeat = 1

do while(stat == 0)
call next_line(unit, line, pos, lnum, stat)
if (stat /= 0) exit

call next_token(line, pos, token)
if (to_lower(line(token%first:token%last)) == '$end') exit

if (iat >= size(sym)) call resize(sym)
if (iat >= size(xyz, 2)) call resize(xyz)
iat = iat + 1

token%last = min(token%last, token%first + symbol_length - 1)
sym(iat) = line(token%first:token%last)
if (to_number(sym(iat)) == 0) then
call read_token(line, token, izp, stat)
sym(iat) = to_symbol(izp)
end if
if (stat /= 0) then
call io_error(error, "Cannot map symbol to atomic number", &
& line, token, filename(unit), lnum, "unknown element")
return
end if

do iz = 1, zrepeat
if (stat == 0) &
call read_next_token(line, pos, token, ij(iz), stat)
if (stat == 0) &
call read_next_token(line, pos, token, zm(iz), stat)
end do
if (stat /= 0) then
call io_error(error, "Cannot read coordinates", &
& line, token, filename(unit), lnum, "expected value")
return
end if

select case(zrepeat)
case(1)
x = 0.0_wp
y = 0.0_wp
z = zm(1)
zrepeat = zrepeat + 1
case(2)
x = 0.0_wp
y = zm(1) * sin(zm(2) * deg_to_rad)
z = zm(1) * cos(zm(2) * deg_to_rad)
zrepeat = zrepeat + 1
case default
x = zm(1) * sin(zm(2) * deg_to_rad) * cos(zm(3) * deg_to_rad)
y = zm(1) * sin(zm(2) * deg_to_rad) * sin(zm(3) * deg_to_rad)
z = zm(1) * cos(zm(2) * deg_to_rad)
end select
if (ij(1) >= iat) then
call io_error(error, "Cannot read coordinates", &
& line, token, filename(unit), lnum, "invalid atom index")
return
end if
xyz(:, iat) = xyz(:, ij(1)) + [x, y, z] * aatoau

call read_next_token(line, pos, token, x, stat)
end do
else
if (stat == 0) &
call read_next_token(line, pos, token, y, stat)
if (stat == 0) &
Expand All @@ -118,7 +185,44 @@
end if

xyz(:, iat) = [x, y, z] * aatoau
end do

do while(stat == 0)
call next_line(unit, line, pos, lnum, stat)
if (stat /= 0) exit

call next_token(line, pos, token)
if (to_lower(line(token%first:token%last)) == '$end') exit

if (iat >= size(sym)) call resize(sym)
if (iat >= size(xyz, 2)) call resize(xyz)
iat = iat + 1

token%last = min(token%last, token%first + symbol_length - 1)
sym(iat) = line(token%first:token%last)
if (to_number(sym(iat)) == 0) then
call read_token(line, token, izp, stat)
sym(iat) = to_symbol(izp)
end if
if (stat /= 0) then
call io_error(error, "Cannot map symbol to atomic number", &
& line, token, filename(unit), lnum, "unknown element")
return
end if

call read_next_token(line, pos, token, x, stat)
if (stat == 0) &
call read_next_token(line, pos, token, y, stat)
if (stat == 0) &
call read_next_token(line, pos, token, z, stat)
if (stat /= 0) then
call io_error(error, "Cannot read coordinates", &
& line, token, filename(unit), lnum, "expected real value")
return

Check warning on line 220 in src/mctc/io/read/qchem.f90

View check run for this annotation

Codecov / codecov/patch

src/mctc/io/read/qchem.f90#L220 

Added line #L220 was not covered by tests
end if

xyz(:, iat) = [x, y, z] * aatoau
end do
end if

if (stat /= 0) then
call io_error(error, "Failed to read molecule block", &
Expand Down
137 changes: 136 additions & 1 deletion test/test_read_qchem.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -36,11 +36,16 @@ subroutine collect_read_qchem(testsuite)
& new_unittest("valid1-qchem", test_valid1_qchem), &
& new_unittest("valid2-qchem", test_valid2_qchem), &
& new_unittest("valid3-qchem", test_valid3_qchem), &
& new_unittest("valid4-qchem", test_valid4_qchem), &
& new_unittest("invalid1-qchem", test_invalid1_qchem, should_fail=.true.), &
& new_unittest("invalid2-qchem", test_invalid2_qchem, should_fail=.true.), &
& new_unittest("invalid3-qchem", test_invalid3_qchem, should_fail=.true.), &
& new_unittest("invalid4-qchem", test_invalid4_qchem, should_fail=.true.), &
& new_unittest("invalid5-qchem", test_invalid5_qchem, should_fail=.true.) &
& new_unittest("invalid5-qchem", test_invalid5_qchem, should_fail=.true.), &
& new_unittest("invalid6-qchem", test_invalid6_qchem, should_fail=.true.), &
& new_unittest("invalid7-qchem", test_invalid7_qchem, should_fail=.true.), &
& new_unittest("invalid8-qchem", test_invalid8_qchem, should_fail=.true.), &
& new_unittest("invalid9-qchem", test_invalid9_qchem, should_fail=.true.) &
& ]

end subroutine collect_read_qchem
Expand Down Expand Up @@ -157,6 +162,36 @@ subroutine test_valid3_qchem(error)

end subroutine test_valid3_qchem

subroutine test_valid4_qchem(error)

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

type(structure_type) :: struc
integer :: unit

open(status='scratch', newunit=unit)
write(unit, '(a)') &
"$molecule", &
"0 1", &
"P", &
"H 1 1.407461", &
"H 1 1.407521 2 100.786448", &
"H 1 1.407521 2 100.786448 3 103.310033 0", &
"O 1 1.487800 2 117.174945 3 -128.344983 0", &
"$end"
rewind(unit)

call read_qchem(struc, unit, error)
close(unit)
if (allocated(error)) return

call check(error, struc%nat, 5, "Number of atoms does not match")
if (allocated(error)) return
call check(error, struc%nid, 3, "Number of species does not match")
if (allocated(error)) return

end subroutine test_valid4_qchem


subroutine test_invalid1_qchem(error)
Expand Down Expand Up @@ -323,5 +358,105 @@ subroutine test_invalid5_qchem(error)

end subroutine test_invalid5_qchem

subroutine test_invalid6_qchem(error)

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

type(structure_type) :: struc
integer :: unit

open(status='scratch', newunit=unit)
write(unit, '(a)') &
"$molecule", &
"0 1", &
"P", &
"H", &
"H 1 1.407521 2 100.786448", &
"H 1 1.407521 2 100.786448 3 103.310033 0", &
"O 1 1.487800 2 117.174945 3 -128.344983 0", &
"$end"
rewind(unit)

call read_qchem(struc, unit, error)
close(unit)

end subroutine test_invalid6_qchem

subroutine test_invalid7_qchem(error)

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

type(structure_type) :: struc
integer :: unit

open(status='scratch', newunit=unit)
write(unit, '(a)') &
"$molecule", &
"0 1", &
"P", &
"H 1", &
"H 1 1.407521 2 100.786448", &
"H 1 1.407521 2 100.786448 3 103.310033 0", &
"O 1 1.487800 2 117.174945 3 -128.344983 0", &
"$end"
rewind(unit)

call read_qchem(struc, unit, error)
close(unit)

end subroutine test_invalid7_qchem

subroutine test_invalid8_qchem(error)

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

type(structure_type) :: struc
integer :: unit

open(status='scratch', newunit=unit)
write(unit, '(a)') &
"$molecule", &
"0 1", &
"P", &
"H 2 1.407461", &
"H 1 1.407521 2 100.786448", &
"H 1 1.407521 2 100.786448 3 103.310033 0", &
"O 1 1.487800 2 117.174945 3 -128.344983 0", &
"$end"
rewind(unit)

call read_qchem(struc, unit, error)
close(unit)

end subroutine test_invalid8_qchem

subroutine test_invalid9_qchem(error)

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

type(structure_type) :: struc
integer :: unit

open(status='scratch', newunit=unit)
write(unit, '(a)') &
"$molecule", &
"0 1", &
"P", &
"*", &
"H 1 1.407521 2 100.786448", &
"H 1 1.407521 2 100.786448 3 103.310033 0", &
"O 1 1.487800 2 117.174945 3 -128.344983 0", &
"$end"
rewind(unit)

call read_qchem(struc, unit, error)
close(unit)

end subroutine test_invalid9_qchem


end module test_read_qchem