From 3069ad496235ba810960832e57a788be27bf8043 Mon Sep 17 00:00:00 2001
From: Sebastian Ehlert <28669218+awvwgk@users.noreply.github.com>
Date: Wed, 26 Jan 2022 21:45:33 +0100
Subject: [PATCH] Implement support for helical boundary conditions in gen
 format

- gen 1D periodicity is along the z-axis stored in the first lattice vector
- the helical angle and the order of the screw axis are stored in the second vector
---
 doc/format-gen.md               |   4 +-
 src/mctc/io/read/aims.f90       |  15 +--
 src/mctc/io/read/genformat.f90  |  65 ++++++++++--
 src/mctc/io/read/turbomole.f90  |  26 +++--
 src/mctc/io/utils.f90           |   5 +
 src/mctc/io/write/aims.f90      |   9 +-
 src/mctc/io/write/genformat.f90 |  24 +++--
 src/mctc/io/write/turbomole.f90 |  15 ++-
 src/mctc/io/write/vasp.f90      |  22 ++--
 test/test_read_aims.f90         | 105 +++++++++++++++++++
 test/test_read_genformat.f90    | 179 ++++++++++++++++++++++++++++++--
 test/test_read_turbomole.f90    |  98 +++++++++++++++++
 12 files changed, 506 insertions(+), 61 deletions(-)

diff --git a/doc/format-gen.md b/doc/format-gen.md
index fba36215..699bb9ff 100644
--- a/doc/format-gen.md
+++ b/doc/format-gen.md
@@ -98,9 +98,7 @@ No extension implemented to the original format.
 
 ## Missing Features
 
-The following features are currently not supported:
-
-- Helical boundary conditions
+The implementation of this format is (to our knowledge) feature-complete.
 
 @Note Feel free to contribute support for missing features
       or bring missing features to our attention by opening an issue.
diff --git a/src/mctc/io/read/aims.f90 b/src/mctc/io/read/aims.f90
index f61bb78d..98ea5b05 100644
--- a/src/mctc/io/read/aims.f90
+++ b/src/mctc/io/read/aims.f90
@@ -47,7 +47,7 @@ subroutine read_aims(mol, unit, error)
    real(wp) :: x, y, z
    character(len=symbol_length), allocatable :: sym(:)
    real(wp), allocatable :: xyz(:, :), abc(:, :), lattice(:, :)
-   logical :: is_frac
+   logical :: is_frac, periodic(3)
    logical, allocatable :: frac(:)
 
    allocate(sym(initial_size), source=repeat(' ', symbol_length))
@@ -57,6 +57,7 @@ subroutine read_aims(mol, unit, error)
 
    iat = 0
    ilt = 0
+   periodic(:) = .false.
 
    lnum = 0
    stat = 0
@@ -99,7 +100,9 @@ subroutine read_aims(mol, unit, error)
          frac(iat) = is_frac
          if (frac(iat)) then
             abc(:, iat) = [x, y, z]
+            xyz(:, iat) = 0.0_wp
          else
+            abc(:, iat) = 0.0_wp
             xyz(:, iat) = [x, y, z] * aatoau
          end if
 
@@ -141,14 +144,12 @@ subroutine read_aims(mol, unit, error)
    end if
 
    if (allocated(lattice)) then
-      if (ilt /= 3) then
-         call fatal_error(error, "Not enough lattice vectors provided")
-         return
-      end if
-      xyz(:, :iat) = xyz(:, :iat) + matmul(lattice, abc(:, :iat))
+      xyz(ilt+1:3, :iat) = xyz(ilt+1:3, :iat) + abc(ilt+1:3, :iat) * aatoau
+      xyz(:ilt, :iat) = xyz(:ilt, :iat) + matmul(lattice(:ilt, :ilt), abc(:ilt, :iat))
+      periodic(:ilt) = .true.
    end if
 
-   call new(mol, sym(:iat), xyz, lattice=lattice)
+   call new(mol, sym(:iat), xyz, lattice=lattice, periodic=periodic)
 
 end subroutine read_aims
 
diff --git a/src/mctc/io/read/genformat.f90 b/src/mctc/io/read/genformat.f90
index e9066638..ac2d58f9 100644
--- a/src/mctc/io/read/genformat.f90
+++ b/src/mctc/io/read/genformat.f90
@@ -15,6 +15,7 @@
 module mctc_io_read_genformat
    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_structure, only : structure_type, new
    use mctc_io_structure_info, only : structure_info
@@ -43,12 +44,12 @@ subroutine read_genformat(mol, unit, error)
 
    character(len=:), allocatable :: line
    integer :: natoms, nspecies, iatom, dummy, isp, ilat, stat, istart, iend
-   logical :: cartesian, periodic
-   real(wp) :: coord(3), lattice(3, 3)
+   logical :: cartesian, periodic(3)
+   real(wp) :: coord(3), origin(3)
    character(len=1) :: variant
    type(token_type) :: token
    character(len=symbol_length), allocatable :: species(:), sym(:)
-   real(wp), allocatable :: xyz(:, :), abc(:, :)
+   real(wp), allocatable :: xyz(:, :), abc(:, :), lattice(:, :)
    type(structure_info) :: info
    integer :: pos, lnum
 
@@ -74,14 +75,20 @@ subroutine read_genformat(mol, unit, error)
    case('s', 'S')
       cartesian = .true.
       periodic = .true.
+      allocate(lattice(3, 3), source=0.0_wp)
    case('f', 'F')
       cartesian = .false.
       periodic = .true.
+      allocate(lattice(3, 3), source=0.0_wp)
+   case('h', 'H')
+      cartesian = .true.
+      periodic = [.false., .false., .true.]
+      allocate(lattice(3, 1), source=0.0_wp)
    case default
       call io_error(error, "Invalid input version found", &
          & line, token, filename(unit), lnum, "unknown identifier")
       return
-   endselect
+   end select
 
    call advance_line(unit, line, pos, lnum, stat)
    isp = 0
@@ -124,13 +131,27 @@ subroutine read_genformat(mol, unit, error)
       end if
    end do
 
-   if (periodic) then
+   if (any(periodic)) then
       call advance_line(unit, line, pos, lnum, stat)
       if (stat /= 0) then
          call io_error(error, "Unexpected end of file", &
             & line, token_type(0, 0), filename(unit), lnum, "missing lattice information")
          return
       end if
+      if (stat == 0) &
+         call read_next_token(line, pos, token, origin(1), stat)
+      if (stat == 0) &
+         call read_next_token(line, pos, token, origin(2), stat)
+      if (stat == 0) &
+         call read_next_token(line, pos, token, origin(3), stat)
+      if (stat /= 0) then
+         call io_error(error, "Cannot read origin", &
+            & line, token, filename(unit), lnum, "expected real value")
+         return
+         end if
+   end if
+
+   if (all(periodic)) then
       do ilat = 1, 3
          call advance_line(unit, line, pos, lnum, stat)
          if (stat == 0) &
@@ -149,12 +170,38 @@ subroutine read_genformat(mol, unit, error)
       if (.not.cartesian) then
          xyz = matmul(lattice, abc)
       end if
-      info = structure_info(cartesian=cartesian)
-      call new(mol, sym, xyz, lattice=lattice, info=info)
-   else
-      call new(mol, sym, xyz)
    end if
 
+   if (count(periodic) == 1) then
+      call advance_line(unit, line, pos, lnum, stat)
+      if (stat == 0) &
+         call read_next_token(line, pos, token, coord(1), stat)
+      if (stat == 0) &
+         call read_next_token(line, pos, token, coord(2), stat)
+      if (stat == 0) &
+         call read_next_token(line, pos, token, coord(3), stat)
+      if (stat /= 0) then
+         call io_error(error, "Cannot read lattice vector", &
+            & line, token, filename(unit), lnum, "expected real value")
+         return
+      end if
+      if (coord(3) < 1) then
+         call io_error(error, "Invalid helical axis rotation order", &
+            & line, token, filename(unit), lnum, "expected positive value")
+         return
+      end if
+
+      ! Store helical axis in *first* lattice vector, however it is not an actual
+      ! lattice vector as on would expect but a screw axis
+      lattice(:, 1) = [coord(1) * aatoau, coord(2) * pi / 180.0_wp, coord(3)]
+   end if
+
+   if (any(periodic)) then
+      xyz(:, :) = xyz - spread(origin, 2, natoms)
+   end if
+
+   info = structure_info(cartesian=cartesian)
+   call new(mol, sym, xyz, lattice=lattice, periodic=periodic, info=info)
 
 contains
 
diff --git a/src/mctc/io/read/turbomole.f90 b/src/mctc/io/read/turbomole.f90
index 64168323..df438115 100644
--- a/src/mctc/io/read/turbomole.f90
+++ b/src/mctc/io/read/turbomole.f90
@@ -78,18 +78,19 @@ subroutine read_coord(mol, unit, error)
    cartesian = .true.
    coord_in_bohr = .true.
    lattice_in_bohr = .true.
-   lattice = 0.0_wp
-   pbc = .false.
+   lattice(:, :) = 0.0_wp
+   pbc(:) = .false.
 
    stat = 0
    call next_line(unit, line, pos, lnum, stat)
    do while(stat == 0)
       if (index(line, flag) == 1) then
          call next_token(line, pos, token)
-         if (index(line, 'end') == 2) then
+         select case(line(token%first:token%last))
+         case('$end')
             exit
 
-         else if (index(line, 'eht') == 2) then
+         case('$eht')
             if (has_eht) then
                pos = 0
                call next_token(line_eht, pos, token2)
@@ -120,7 +121,7 @@ subroutine read_coord(mol, unit, error)
                return
             end if
 
-         else if (index(line, 'coord') == 2) then
+         case('$coord')
             if (has_coord) then
                pos = 0
                call next_token(line_coord, pos, token2)
@@ -164,7 +165,7 @@ subroutine read_coord(mol, unit, error)
             end do coord_group
             cycle
 
-         else if (index(line, 'periodic') == 2) then
+         case('$periodic')
             if (has_periodic) then
                pos = 0
                call next_token(line_periodic, pos, token2)
@@ -185,7 +186,7 @@ subroutine read_coord(mol, unit, error)
                return
             end if
 
-         else if (index(line, 'lattice') == 2) then
+         case('$lattice')
             if (has_lattice) then
                pos = 0
                call next_token(line_lattice, pos, token2)
@@ -210,7 +211,7 @@ subroutine read_coord(mol, unit, error)
             end do lattice_group
             cycle
 
-         else if (index(line, 'cell') == 2) then
+         case('$cell')
             if (has_cell) then
                pos = 0
                call next_token(line_cell, pos, token2)
@@ -228,7 +229,7 @@ subroutine read_coord(mol, unit, error)
             call next_line(unit, cell_string, pos, lnum, stat)
             if (debug) print*, cell_string
 
-         end if
+         end select
       end if
       token = token_type(0, 0)
       call next_line(unit, line, pos, lnum, stat)
@@ -358,7 +359,10 @@ subroutine read_coord(mol, unit, error)
       end if
       xyz(:, :) = coord(:, :natoms) * conv
    else
-      xyz = matmul(lattice, coord)
+      ! Non-periodic coordinates are in Bohr
+      xyz(periodic+1:3, :) = coord(periodic+1:3, :natoms)
+      ! Periodic coordinates must still be transformed with lattice
+      xyz(:periodic, :) = matmul(lattice(:periodic, :periodic), coord(:periodic, :natoms))
    end if
 
    ! save data on input format
@@ -387,7 +391,7 @@ pure subroutine cell_to_dlat(cellpar, lattice)
    real(wp), intent(in)  :: cellpar(6)
 
    !> Direct lattice
-   real(wp), intent(out) :: lattice(3, 3)
+   real(wp), intent(out) :: lattice(:, :)
 
    real(wp) :: dvol
 
diff --git a/src/mctc/io/utils.f90 b/src/mctc/io/utils.f90
index 2a1da780..005966b8 100644
--- a/src/mctc/io/utils.f90
+++ b/src/mctc/io/utils.f90
@@ -130,6 +130,11 @@ subroutine next_token(string, pos, token)
 
    integer :: start
 
+   if (pos >= len(string)) then
+      token = token_type(len(string)+1, len(string)+1)
+      return
+   end if
+
    do while(pos < len(string))
       pos = pos + 1
       select case(string(pos:pos))
diff --git a/src/mctc/io/write/aims.f90 b/src/mctc/io/write/aims.f90
index 3820f687..ded0a50f 100644
--- a/src/mctc/io/write/aims.f90
+++ b/src/mctc/io/write/aims.f90
@@ -48,9 +48,12 @@ subroutine write_aims(self, unit)
    end if
 
    if (any(self%periodic)) then
-      do ilt = 1, size(self%lattice, 2)
-         write(unit, '(a, 1x, 3f24.14)') &
-            "lattice_vector", self%lattice(:, ilt) * autoaa
+      if (size(self%lattice, 2) /= 3) return
+      do ilt = 1, 3
+         if (self%periodic(ilt)) then
+            write(unit, '(a, 1x, 3f24.14)') &
+               "lattice_vector", self%lattice(:, ilt) * autoaa
+         end if
       end do
    end if
 
diff --git a/src/mctc/io/write/genformat.f90 b/src/mctc/io/write/genformat.f90
index 15b5cac0..277f6d29 100644
--- a/src/mctc/io/write/genformat.f90
+++ b/src/mctc/io/write/genformat.f90
@@ -14,6 +14,7 @@
 
 module mctc_io_write_genformat
    use mctc_env_accuracy, only : wp
+   use mctc_io_constants, only : pi
    use mctc_io_convert, only : autoaa
    use mctc_io_math, only : matinv_3x3
    use mctc_io_symbols, only : to_symbol
@@ -34,16 +35,23 @@ subroutine write_genformat(mol, unit)
    real(wp), parameter :: zero3(3) = 0.0_wp
    real(wp), allocatable :: inv_lat(:, :)
    real(wp), allocatable :: abc(:, :)
+   logical :: helical
 
+   helical = .false.
    write(unit, '(i0, 1x)', advance='no') mol%nat
    if (.not.any(mol%periodic)) then
       write(unit, '("C")') ! cluster
    else
-      if (mol%info%cartesian) then
-         write(unit, '("S")') ! supercell
+      helical = count(mol%periodic) == 1 .and. mol%periodic(3) .and. size(mol%lattice, 2) == 1
+      if (helical) then
+         write(unit, '("H")') ! helical
       else
-         write(unit, '("F")') ! fractional
-      endif
+         if (mol%info%cartesian) then
+            write(unit, '("S")') ! supercell
+         else
+            write(unit, '("F")') ! fractional
+         endif
+      end if
    endif
 
    do izp = 1, mol%nid
@@ -66,10 +74,14 @@ subroutine write_genformat(mol, unit)
    endif
 
    if (any(mol%periodic)) then
-      ! scaling factor for lattice parameters is always one
       write(unit, '(3f20.14)') zero3
       ! write the lattice parameters
-      write(unit, '(3f20.14)') mol%lattice(:, :)*autoaa
+      if (helical) then
+         write(unit, '(2f20.14,1x,i0)') &
+            & mol%lattice(3, 1)*autoaa, mol%lattice(2, 2)*180.0_wp/pi, nint(mol%lattice(3, 2))
+      else
+         write(unit, '(3f20.14)') mol%lattice(:, :)*autoaa
+      end if
    endif
 
 end subroutine write_genformat
diff --git a/src/mctc/io/write/turbomole.f90 b/src/mctc/io/write/turbomole.f90
index 564d20c8..b666d156 100644
--- a/src/mctc/io/write/turbomole.f90
+++ b/src/mctc/io/write/turbomole.f90
@@ -27,19 +27,24 @@ module mctc_io_write_turbomole
 subroutine write_coord(mol, unit)
    class(structure_type), intent(in) :: mol
    integer, intent(in) :: unit
-   integer :: iat
+   integer :: iat, ilt, npbc
 
    write(unit, '(a)') "$coord"
    do iat = 1, mol%nat
       write(unit, '(3es24.14, 6x, a)') mol%xyz(:, iat), trim(mol%sym(mol%id(iat)))
-   enddo
+   end do
    write(unit, '(a, *(1x, a, "=", i0))') &
       "$eht", "charge", nint(mol%charge), "unpaired", mol%uhf
    write(unit, '(a, 1x, i0)') "$periodic", count(mol%periodic)
    if (any(mol%periodic)) then
-      write(unit, '(a)') "$lattice bohr"
-      write(unit, '(3f20.14)') mol%lattice
-   endif
+      npbc = count(mol%periodic)
+      if (size(mol%lattice, 2) == 3) then
+         write(unit, '(a)') "$lattice bohr"
+         do ilt = 1, npbc
+            write(unit, '(3f20.14)') mol%lattice(:npbc, ilt)
+         end do
+      end if
+   end if
    write(unit, '(a)') "$end"
 
 end subroutine write_coord
diff --git a/src/mctc/io/write/vasp.f90 b/src/mctc/io/write/vasp.f90
index 98b7d3b9..6e09481d 100644
--- a/src/mctc/io/write/vasp.f90
+++ b/src/mctc/io/write/vasp.f90
@@ -47,8 +47,8 @@ subroutine write_vasp(self, unit, comment_line)
          j = j+1
          izp = self%id(i)
          species(j) = self%id(i)
-      endif
-   enddo
+      end if
+   end do
 
    ! use vasp 5.x format
    if (present(comment_line)) then
@@ -65,22 +65,24 @@ subroutine write_vasp(self, unit, comment_line)
    write(unit, '(f20.14)') self%info%scale
    ! write the lattice parameters
    if (any(self%periodic)) then
-      do i = 1, size(self%lattice, 2)
-         write(unit, '(3f20.14)') self%lattice(:, i)*autoaa/self%info%scale
-      enddo
+      if (size(self%lattice, 2) == 3) then
+         write(unit, '(3f20.14)') self%lattice
+      else
+         write(unit, '(3f20.14)') spread(0.0_wp, 1, 9)
+      end if
    else
       write(unit, '(3f20.14)') spread(0.0_wp, 1, 9)
    end if
 
    do i = 1, j
       write(unit, '(1x, a)', advance='no') self%sym(species(i))
-   enddo
+   end do
    write(unit, '(a)')
 
    ! write the count of the consequtive atom types
    do i = 1, j
       write(unit, '(1x, i0)', advance='no') kinds(i)
-   enddo
+   end do
    write(unit, '(a)')
    deallocate(kinds, species)
 
@@ -93,7 +95,7 @@ subroutine write_vasp(self, unit, comment_line)
       ! now write the cartesian coordinates
       do i = 1, self%nat
          write(unit, '(3f20.14)') self%xyz(:, i)*autoaa/self%info%scale
-      enddo
+      end do
    else
       write(unit, '("Direct")')
       inv_lat = matinv_3x3(self%lattice)
@@ -102,8 +104,8 @@ subroutine write_vasp(self, unit, comment_line)
       ! now write the fractional coordinates
       do i = 1, self%nat
          write(unit, '(3f20.14)') abc(:, i)
-      enddo
-   endif
+      end do
+   end if
 
 end subroutine write_vasp
 
diff --git a/test/test_read_aims.f90 b/test/test_read_aims.f90
index 63cb5e63..e5d2d806 100644
--- a/test/test_read_aims.f90
+++ b/test/test_read_aims.f90
@@ -13,6 +13,7 @@
 ! limitations under the License.
 
 module test_read_aims
+   use mctc_env, only : wp
    use mctc_env_testing, only : new_unittest, unittest_type, error_type, check
    use mctc_io_read_aims
    use mctc_io_structure
@@ -37,6 +38,8 @@ subroutine collect_read_aims(testsuite)
       & new_unittest("valid3-aims", test_valid3_aims), &
       & new_unittest("valid4-aims", test_valid4_aims), &
       & new_unittest("valid5-aims", test_valid5_aims), &
+      & new_unittest("valid6-aims", test_valid6_aims), &
+      & new_unittest("valid7-aims", test_valid7_aims), &
       & new_unittest("invalid1-aims", test_invalid1_aims, should_fail=.true.), &
       & new_unittest("invalid2-aims", test_invalid2_aims, should_fail=.true.), &
       & new_unittest("invalid3-aims", test_invalid3_aims, should_fail=.true.), &
@@ -250,6 +253,108 @@ subroutine test_valid5_aims(error)
 end subroutine test_valid5_aims
 
 
+subroutine test_valid6_aims(error)
+
+   !> Error handling
+   type(error_type), allocatable, intent(out) :: error
+
+   type(structure_type) :: struc1, struc2
+   integer :: unit
+
+   open(status='scratch', newunit=unit)
+   write(unit, '(a)') &
+      "atom            0.00000000000000   0.00000000000000   0.00000000000000 Mg", &
+      "atom            1.48881365396205  -1.48881365396205   0.00000000000000 O", &
+      "atom            1.48881365396205   1.48881365396205   0.00000000000000 O", &
+      "atom            0.00000000000000   0.00000000000000   2.10550046127949 O", &
+      "atom            2.97762730792410   0.00000000000000   0.00000000000000 Mg", &
+      "atom            1.48881365396205  -1.48881365396205   2.10550046127949 Mg", &
+      "atom            1.48881365396205   1.48881365396205   2.10550046127949 Mg", &
+      "atom            2.97762730792410   0.00000000000000   2.10550046127949 O", &
+      "lattice_vector  2.97762730792410  -2.97762730792410   0.00000000000000", &
+      "lattice_vector  2.97762730792410   2.97762730792410   0.00000000000000"
+   rewind(unit)
+
+   call read_aims(struc1, unit, error)
+   close(unit)
+   if (allocated(error)) return
+
+   open(status='scratch', newunit=unit)
+   write(unit, '(a)') &
+      "atom_frac       0.00000000000000   0.00000000000000   0.00000000000000 Mg", &
+      "atom_frac       0.50000000000000   0.00000000000000   0.00000000000000 O", &
+      "atom_frac       0.00000000000000   0.50000000000000   0.00000000000000 O", &
+      "atom_frac       0.00000000000000   0.00000000000000   2.10550046127949 O", &
+      "atom_frac       0.50000000000000   0.50000000000000   0.00000000000000 Mg", &
+      "atom_frac       0.50000000000000   0.00000000000000   2.10550046127949 Mg", &
+      "atom_frac       0.00000000000000   0.50000000000000   2.10550046127949 Mg", &
+      "atom_frac       0.50000000000000   0.50000000000000   2.10550046127949 O", &
+      "lattice_vector  2.97762730792410  -2.97762730792410   0.00000000000000", &
+      "lattice_vector  2.97762730792410   2.97762730792410   0.00000000000000"
+   rewind(unit)
+
+   call read_aims(struc2, unit, error)
+   close(unit)
+   if (allocated(error)) return
+
+   call check(error, norm2(struc1%xyz - struc2%xyz), 0.0_wp, "Coordinates do not match")
+   if (allocated(error)) return
+
+end subroutine test_valid6_aims
+
+
+subroutine test_valid7_aims(error)
+
+   !> Error handling
+   type(error_type), allocatable, intent(out) :: error
+
+   type(structure_type) :: struc
+   integer :: unit
+
+   open(status='scratch', newunit=unit)
+   write(unit, '(a)') &
+      "atom             -1.05835465887935   1.85522662363901   0.00000000000000 B", &
+      "atom             -1.05835465887935   1.57910813351869   1.38575958673374 N", &
+      "atom             -1.05835465887935   0.79318285794365  -0.93200541748473 N", &
+      "atom             -1.05835465887935   2.94621239911295  -0.36993970607209 H", &
+      "atom             -1.05835465887935   0.24094589146163   1.83951376127452 B", &
+      "atom             -1.05835465887935  -0.54497939258025  -0.47825125458585 B", &
+      "atom             -1.05835465887935  -0.82109787740880   0.90750834908157 N", &
+      "atom             -1.05835465887935   0.01583015330186   2.96930500972370 H", &
+      "atom             -1.05835465887935  -1.41084943428660  -1.23810284035547 H", &
+      "atom             -1.05835465887935   2.39762594336943   2.10405675666812 H", &
+      "atom             -1.05835465887935  -1.85242037510793   1.25721697940438 H", &
+      "atom             -1.05835465887935   1.00598756166737  -2.00001121245014 H", &
+      "atom              1.05835465887935   0.82109787740880  -0.90750833320625 B", &
+      "atom              1.05835465887935   0.54497938728848   0.47825125193996 N", &
+      "atom              1.05835465887935  -0.24094588775739  -1.83951375598275 N", &
+      "atom              1.05835465887935   1.91208365288273  -1.27744804245341 H", &
+      "atom              1.05835465887935  -0.79318285794365   0.93200542277650 B", &
+      "atom              1.05835465887935  -1.57910813881047  -1.38575959202551 B", &
+      "atom              1.05835465887935  -1.85522662893078   0.00000001308910 N", &
+      "atom              1.05835465887935  -1.01829859700301   2.06179667651746 H", &
+      "atom              1.05835465887935  -2.44497818051681  -2.14561117356172 H", &
+      "atom              1.05835465887935   1.36349719184744   1.19654842346187 H", &
+      "atom              1.05835465887935  -2.88654912133814   0.34970864461060 H", &
+      "atom              1.05835465887935  -0.02814118763207  -2.90751955094817 H", &
+      "lattice_vector    4.23341864610095   0.00000000000000   0.00000000000000"
+   rewind(unit)
+
+   call read_aims(struc, unit, error)
+   close(unit)
+   if (allocated(error)) return
+
+   call check(error, struc%nat, 24, "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
+   call check(error, count(struc%periodic), 1, "Periodicity does not match")
+   if (allocated(error)) return
+
+end subroutine test_valid7_aims
+
+
+
 subroutine test_invalid1_aims(error)
 
    !> Error handling
diff --git a/test/test_read_genformat.f90 b/test/test_read_genformat.f90
index d38b8cf7..b1952056 100644
--- a/test/test_read_genformat.f90
+++ b/test/test_read_genformat.f90
@@ -36,12 +36,17 @@ subroutine collect_read_genformat(testsuite)
       & new_unittest("valid2-gen", test_valid2_gen), &
       & new_unittest("valid3-gen", test_valid3_gen), &
       & new_unittest("valid4-gen", test_valid4_gen), &
+      & new_unittest("valid5-gen", test_valid5_gen), &
+      & new_unittest("valid6-gen", test_valid6_gen), &
       & new_unittest("invalid1-gen", test_invalid1_gen, should_fail=.true.), &
       & new_unittest("invalid2-gen", test_invalid2_gen, should_fail=.true.), &
       & new_unittest("invalid3-gen", test_invalid3_gen, should_fail=.true.), &
       & new_unittest("invalid4-gen", test_invalid4_gen, should_fail=.true.), &
       & new_unittest("invalid5-gen", test_invalid5_gen, should_fail=.true.), &
-      & new_unittest("invalid6-gen", test_invalid6_gen, should_fail=.true.) &
+      & new_unittest("invalid6-gen", test_invalid6_gen, should_fail=.true.), &
+      & new_unittest("invalid7-gen", test_invalid7_gen, should_fail=.true.), &
+      & new_unittest("invalid8-gen", test_invalid8_gen, should_fail=.true.), &
+      & new_unittest("invalid9-gen", test_invalid9_gen, should_fail=.true.) &
       & ]
 
 end subroutine collect_read_genformat
@@ -182,6 +187,88 @@ subroutine test_valid4_gen(error)
 end subroutine test_valid4_gen
 
 
+subroutine test_valid5_gen(error)
+
+   !> Error handling
+   type(error_type), allocatable, intent(out) :: error
+
+   type(structure_type) :: struc
+   integer :: unit
+
+   open(status='scratch', newunit=unit)
+   write(unit, '(a)') &
+      "   20  H", &
+      "  C", &
+      "    1 1    0.2756230044E+01    0.2849950460E+01    0.1794011798E+01", &
+      "    2 1    0.2656226397E+01    0.2949964389E+01    0.3569110265E+00", &
+      "    3 1    0.4149823216E+00    0.3947943175E+01    0.1774023191E+01", &
+      "    4 1   -0.1984731085E+01    0.3437783679E+01    0.1784008240E+01", &
+      "    5 1   -0.3626350732E+01    0.1614594006E+01    0.1784022394E+01", &
+      "    6 1   -0.3882893767E+01   -0.8252790149E+00    0.1784006601E+01", &
+      "    7 1   -0.2656230041E+01   -0.2949950471E+01    0.1784011798E+01", &
+      "    8 1   -0.4149823216E+00   -0.3947943175E+01    0.1784023191E+01", &
+      "    9 1    0.1984731085E+01   -0.3437783679E+01    0.1784008240E+01", &
+      "   10 1    0.3626350732E+01   -0.1614594006E+01    0.1784022394E+01", &
+      "   11 1    0.3882893767E+01    0.8252790258E+00    0.1784006601E+01", &
+      "   12 1    0.4149905833E+00    0.3947943870E+01    0.3569255177E+00", &
+      "   13 1   -0.1984725150E+01    0.3437762712E+01    0.3569151866E+00", &
+      "   14 1   -0.3626358050E+01    0.1614595957E+01    0.3569260541E+00", &
+      "   15 1   -0.3882900023E+01   -0.8252696970E+00    0.3569133218E+00", &
+      "   16 1   -0.2656226396E+01   -0.2949964400E+01    0.3569110265E+00", &
+      "   17 1   -0.4149905833E+00   -0.3947943870E+01    0.3569255177E+00", &
+      "   18 1    0.1984725150E+01   -0.3437762712E+01    0.3569151866E+00", &
+      "   19 1    0.3626358050E+01   -0.1614595957E+01    0.3569260541E+00", &
+      "   20 1    0.3882900026E+01    0.8252697074E+00    0.3569133218E+00", &
+      "    0 0 0", &
+      "    0.2140932670E+01   18.0 1"
+   rewind(unit)
+
+   call read_genformat(struc, unit, error)
+   close(unit)
+   if (allocated(error)) return
+
+   call check(error, count(struc%periodic), 1, "Incorrect periodicity")
+   if (allocated(error)) return
+   call check(error, struc%nat, 20, "Number of atoms does not match")
+   if (allocated(error)) return
+   call check(error, struc%nid, 1, "Number of species does not match")
+   if (allocated(error)) return
+
+end subroutine test_valid5_gen
+
+
+subroutine test_valid6_gen(error)
+
+   !> Error handling
+   type(error_type), allocatable, intent(out) :: error
+
+   type(structure_type) :: struc
+   integer :: unit
+
+   open(status='scratch', newunit=unit)
+   write(unit, '(a)') &
+      "    2  H", &
+      "  C", &
+      "    1 1    0.0 0.0 1.4271041431", &
+      "    2 1    0.0 0.0 0.0", &
+      "   -0.2703556133E+01  -0.2906666140E+01 -0.3618948259E+00", &
+      "    0.2140932670E+01   18.00000000 10"
+   rewind(unit)
+
+   call read_genformat(struc, unit, error)
+   close(unit)
+   if (allocated(error)) return
+
+   call check(error, count(struc%periodic), 1, "Incorrect periodicity")
+   if (allocated(error)) return
+   call check(error, struc%nat, 2, "Number of atoms does not match")
+   if (allocated(error)) return
+   call check(error, struc%nid, 1, "Number of species does not match")
+   if (allocated(error)) return
+
+end subroutine test_valid6_gen
+
+
 subroutine test_invalid1_gen(error)
 
    !> Error handling
@@ -204,7 +291,6 @@ subroutine test_invalid1_gen(error)
 
    call read_genformat(struc, unit, error)
    close(unit)
-   if (allocated(error)) return
 
 end subroutine test_invalid1_gen
 
@@ -231,7 +317,6 @@ subroutine test_invalid2_gen(error)
 
    call read_genformat(struc, unit, error)
    close(unit)
-   if (allocated(error)) return
 
 end subroutine test_invalid2_gen
 
@@ -258,7 +343,6 @@ subroutine test_invalid3_gen(error)
 
    call read_genformat(struc, unit, error)
    close(unit)
-   if (allocated(error)) return
 
 end subroutine test_invalid3_gen
 
@@ -279,7 +363,6 @@ subroutine test_invalid4_gen(error)
 
    call read_genformat(struc, unit, error)
    close(unit)
-   if (allocated(error)) return
 
 end subroutine test_invalid4_gen
 
@@ -302,7 +385,6 @@ subroutine test_invalid5_gen(error)
 
    call read_genformat(struc, unit, error)
    close(unit)
-   if (allocated(error)) return
 
 end subroutine test_invalid5_gen
 
@@ -329,9 +411,92 @@ subroutine test_invalid6_gen(error)
 
    call read_genformat(struc, unit, error)
    close(unit)
-   if (allocated(error)) return
 
 end subroutine test_invalid6_gen
 
 
+subroutine test_invalid7_gen(error)
+
+   !> Error handling
+   type(error_type), allocatable, intent(out) :: error
+
+   type(structure_type) :: struc
+   integer :: unit
+
+   open(status='scratch', newunit=unit)
+   write(unit, '(a)') &
+      "12 H", &
+      " C H ", &
+      " 1 1   1.39792890   0.00000000  -0.00000000", &
+      " 2 2   2.49455487  -0.00000000   0.00000000", &
+      " 3 1   0.69896445   1.21064194  -0.00000000", &
+      " 4 2   1.24727743   2.16034789   0.00000000", &
+      " 5 1  -0.69896445   1.21064194  -0.00000000", &
+      " 6 2  -1.24727743   2.16034789   0.00000000", &
+      " 7 1  -1.39792890  -0.00000000  -0.00000000", &
+      " 8 2  -2.49455487   0.00000000   0.00000000", &
+      " 9 1  -0.69896445  -1.21064194  -0.00000000", &
+      "10 2  -1.24727743  -2.16034789   0.00000000", &
+      "11 1   0.69896445  -1.21064194  -0.00000000", &
+      "12 2   1.24727743  -2.16034789   0.00000000", &
+      "  0 0 0", &
+      "  3.0 ***** 1"
+   rewind(unit)
+
+   call read_genformat(struc, unit, error)
+   close(unit)
+
+end subroutine test_invalid7_gen
+
+
+subroutine test_invalid8_gen(error)
+
+   !> Error handling
+   type(error_type), allocatable, intent(out) :: error
+
+   type(structure_type) :: struc
+   integer :: unit
+
+   open(status='scratch', newunit=unit)
+   write(unit, '(a)') &
+      "   2  H", &
+      "  C", &
+      "    1 1    0.2756230044E+01    0.2849950460E+01    0.1794011798E+01", &
+      "    2 1    0.2656226397E+01    0.2949964389E+01    0.3569110265E+00", &
+      "    0 0 0", &
+      "    0.2140932670E+01 18.0 -10"
+   rewind(unit)
+
+   call read_genformat(struc, unit, error)
+   close(unit)
+
+end subroutine test_invalid8_gen
+
+
+subroutine test_invalid9_gen(error)
+
+   !> Error handling
+   type(error_type), allocatable, intent(out) :: error
+
+   type(structure_type) :: struc
+   integer :: unit
+
+   open(status='scratch', newunit=unit)
+   write(unit, '(a)') &
+      "    2  S", &
+      "  C", &
+      "    1 1    0.0000000000E+00    0.0000000000E+00    0.0000000000E+00", &
+      "    2 1    1.5000000000E+00    0.0000000000E+00    0.0000000000E+00", &
+      "    0.0000000000E+00    0.0000000000E+00", &
+      "    0.2000000000E+01    0.0000000000E+00    0.0000000000E+00", &
+      "    0.0000000000E+00    0.1000000000E+03    0.0000000000E+00", &
+      "    0.0000000000E+00    0.0000000000E+00    0.1000000000E+03"
+   rewind(unit)
+
+   call read_genformat(struc, unit, error)
+   close(unit)
+
+end subroutine test_invalid9_gen
+
+
 end module test_read_genformat
diff --git a/test/test_read_turbomole.f90 b/test/test_read_turbomole.f90
index 07caf199..966b14e7 100644
--- a/test/test_read_turbomole.f90
+++ b/test/test_read_turbomole.f90
@@ -42,6 +42,8 @@ subroutine collect_read_turbomole(testsuite)
       & new_unittest("valid7-coord", test_valid7_coord), &
       & new_unittest("valid8-coord", test_valid8_coord), &
       & new_unittest("valid9-coord", test_valid9_coord), &
+      & new_unittest("valid10-coord", test_valid10_coord), &
+      & new_unittest("valid11-coord", test_valid11_coord), &
       & new_unittest("invalid1-coord", test_invalid1_coord, should_fail=.true.), &
       & new_unittest("invalid2-coord", test_invalid2_coord, should_fail=.true.), &
       & new_unittest("invalid3-coord", test_invalid3_coord, should_fail=.true.), &
@@ -422,6 +424,102 @@ subroutine test_valid9_coord(error)
 end subroutine test_valid9_coord
 
 
+subroutine test_valid10_coord(error)
+
+   !> Error handling
+   type(error_type), allocatable, intent(out) :: error
+
+   type(structure_type) :: struc
+   integer :: unit
+
+   open(status='scratch', newunit=unit)
+   write(unit, '(a)') &
+      "$cell", &
+      " 8.00000006 ", &
+      "$periodic 1", &
+      "$coord", &
+      "   -2.00000001000000      3.50586945000000      0.00000000000000      b", &
+      "   -2.00000001000000      2.98408124000000      2.61870552000000      n", &
+      "   -2.00000001000000      1.49889804000000     -1.76123460000000      n", &
+      "   -2.00000001000000      5.56753332000000     -0.69908457400000      h", &
+      "   -2.00000001000000      0.45532164600000      3.47617645000000      b", &
+      "   -2.00000001000000     -1.02986157000000     -0.90376369200000      b", &
+      "   -2.00000001000000     -1.55164977000000      1.71494186000000      n", &
+      "   -2.00000001000000      0.02991464770000      5.61117202000000      h", &
+      "   -2.00000001000000     -2.66611845000000     -2.33967477000000      h", &
+      "   -2.00000001000000      4.53085539000000      3.97609015000000      h", &
+      "   -2.00000001000000     -3.50056641000000      2.37579525000000      h", &
+      "   -2.00000001000000      1.90104056000000     -3.77947261000000      h", &
+      "    2.00000001000000      1.55164977000000     -1.71494183000000      b", &
+      "    2.00000001000000      1.02986156000000      0.90376368700000      n", &
+      "    2.00000001000000     -0.45532163900000     -3.47617644000000      n", &
+      "    2.00000001000000      3.61331364000000     -2.41402641000000      h", &
+      "    2.00000001000000     -1.49889804000000      1.76123461000000      b", &
+      "    2.00000001000000     -2.98408125000000     -2.61870553000000      b", &
+      "    2.00000001000000     -3.50586946000000      0.00000002473480      n", &
+      "    2.00000001000000     -1.92430504000000      3.89623019000000      h", &
+      "    2.00000001000000     -4.62033813000000     -4.05461660000000      h", &
+      "    2.00000001000000      2.57663570000000      2.26114832000000      h", &
+      "    2.00000001000000     -5.45478609000000      0.66085341700000      h", &
+      "    2.00000001000000     -0.05317912580000     -5.49441445000000      h", &
+      "$user-defined bonds", &
+      "$end"
+   rewind(unit)
+
+   call read_coord(struc, unit, error)
+   close(unit)
+   if (allocated(error)) return
+
+   call check(error, struc%nat, 24, "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
+   call check(error, count(struc%periodic), 1, "Periodic of system does not match")
+   if (allocated(error)) return
+
+end subroutine test_valid10_coord
+
+
+subroutine test_valid11_coord(error)
+
+   !> Error handling
+   type(error_type), allocatable, intent(out) :: error
+
+   type(structure_type) :: struc
+   integer :: unit
+
+   open(status='scratch', newunit=unit)
+   write(unit, '(a)') &
+      "$coord frac", &
+      "    0.00000000000000      0.00000000000000      0.00000000000000      mg", &
+      "    0.50000000000000      0.00000000000000      0.00000000000000      o", &
+      "    0.00000000000000      0.50000000000000      0.00000000000000      o", &
+      "    0.00000000000000      0.00000000000000      3.97881835572287      o", &
+      "    0.50000000000000      0.50000000000000      0.00000000000000      mg", &
+      "    0.50000000000000      0.00000000000000      3.97881835572287      mg", &
+      "    0.00000000000000      0.50000000000000      3.97881835572287      mg", &
+      "    0.50000000000000      0.50000000000000      3.97881835572287      o", &
+      "$periodic 2", &
+      "$lattice", &
+      " 5.626898880882 -5.626898880882", &
+      " 5.626898880882  5.626898880882", &
+      "$end"
+   rewind(unit)
+
+   call read_coord(struc, unit, error)
+   close(unit)
+   if (allocated(error)) return
+
+   call check(error, struc%nat, 8, "Number of atoms does not match")
+   if (allocated(error)) return
+   call check(error, struc%nid, 2, "Number of species does not match")
+   if (allocated(error)) return
+   call check(error, count(struc%periodic), 2, "Periodic of system does not match")
+   if (allocated(error)) return
+
+end subroutine test_valid11_coord
+
+
 subroutine test_invalid1_coord(error)
 
    !> Error handling