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

@@ -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.

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
 

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
 

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
 

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))

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
 

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

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