Update Turbomole reader

- do not backspace on unit to allow reading from stdin
- support eht group in coord format for reading of charge/uhf info

doc/format-tmol.md

@@ -136,12 +136,14 @@ The original format does only allow for the ``periodic`` or ``eht`` group to
 appear in the ``control`` file, to make the format self-contained, all groups
 must appear in the same file.
 
+The ``coord`` group only supports the ``frac`` modifier in Turbomole, but this
+reader also allows ``angs`` and ``bohr``.
+
 ## Missing Features
 
 The following features are currently not supported:
 
 - Preserving information about frozen atoms from ``coord`` data group
-- Reading charge and spin from the ``eht`` data group
 
 @Note Feel free to contribute support for missing features
       or bring missing features to our attention by opening an issue.

src/mctc/io/read/turbomole.f90

@@ -49,11 +49,12 @@ subroutine read_coord(mol, unit, error)
    integer, parameter :: p_initial_size = 100
    integer, parameter :: p_nlv(3) = [1, 4, 9], p_ncp(3) = [1, 3, 6]
 
-   logical :: has_coord, has_periodic, has_lattice, has_cell
+   logical :: has_coord, has_periodic, has_lattice, has_cell, has_eht
    logical :: cartesian, coord_in_bohr, lattice_in_bohr, pbc(3)
-   integer :: stat, iatom, i, j, natoms, periodic, cell_vectors
+   integer :: stat, iatom, i, j, natoms, periodic, cell_vectors, icharge
+   integer, allocatable :: unpaired
    real(wp) :: latvec(9), conv, cellpar(6), lattice(3, 3)
-   real(wp), allocatable :: coord(:, :), xyz(:, :)
+   real(wp), allocatable :: coord(:, :), xyz(:, :), charge
    character(len=:), allocatable :: line, cell_string, lattice_string
    character(len=symbol_length), allocatable :: sym(:)
    type(structure_info) :: info
@@ -64,6 +65,7 @@ subroutine read_coord(mol, unit, error)
    iatom = 0
    periodic = 0
    cell_vectors = 0
+   has_eht = .false.
    has_coord = .false.
    has_periodic = .false.
    has_lattice = .false.
@@ -75,27 +77,38 @@ subroutine read_coord(mol, unit, error)
    pbc = .false.
 
    stat = 0
+   call getline(unit, line, stat)
    do while(stat == 0)
-      call getline(unit, line, stat)
       if (index(line, flag) == 1) then
          if (index(line, 'end') == 2) then
             exit
 
+         else if (.not.has_eht .and. index(line, 'eht') == 2) then
+            has_eht = .true.
+            i = index(line, 'charge=')
+            if (i > 0) then
+               read(line(i+7:), *, iostat=stat) icharge
+               charge = real(icharge, wp)
+            end if
+            j = index(line, 'unpaired=')
+            if (j > 0) then
+               allocate(unpaired)
+               read(line(j+9:), *, iostat=stat) unpaired
+            end if
+
          else if (.not.has_coord .and. index(line, 'coord') == 2) then
             has_coord = .true.
             ! $coord angs / $coord bohr / $coord frac
             call select_unit(line, coord_in_bohr, cartesian)
             coord_group: do while(stat == 0)
                call getline(unit, line, stat)
-               if (index(line, flag) == 1) then
-                  backspace(unit)
-                  exit coord_group
-               end if
+               if (index(line, flag) == 1) exit coord_group
                if (iatom >= size(coord, 2)) call resize(coord)
                if (iatom >= size(sym)) call resize(sym)
                iatom = iatom + 1
                read(line, *, iostat=stat) coord(:, iatom), sym(iatom)
             end do coord_group
+            cycle
 
          else if (.not.has_periodic .and. index(line, 'periodic') == 2) then
             has_periodic = .true.
@@ -110,13 +123,11 @@ subroutine read_coord(mol, unit, error)
             lattice_string = ''
             lattice_group: do while(stat == 0)
                call getline(unit, line, stat)
-               if (index(line, flag) == 1) then
-                  backspace(unit)
-                  exit lattice_group
-               end if
+               if (index(line, flag) == 1) exit lattice_group
                cell_vectors = cell_vectors + 1
                lattice_string = lattice_string // ' ' // line
             end do lattice_group
+            cycle
 
          else if (.not.has_cell .and. index(line, 'cell') == 2) then
             has_cell = .true.
@@ -127,6 +138,7 @@ subroutine read_coord(mol, unit, error)
 
          end if
       end if
+      call getline(unit, line, stat)
    end do
 
    if (.not.has_coord .or. iatom == 0) then
@@ -217,7 +229,8 @@ subroutine read_coord(mol, unit, error)
    ! save data on input format
    info = structure_info(cartesian=cartesian, lattice=has_lattice, &
       & angs_lattice=.not.lattice_in_bohr, angs_coord=.not.coord_in_bohr)
-   call new(mol, sym(:natoms), xyz, lattice=lattice, periodic=pbc, info=info)
+   call new(mol, sym(:natoms), xyz, charge=charge, uhf=unpaired, &
+      & lattice=lattice, periodic=pbc, info=info)
 
 contains
 

src/mctc/io/write/turbomole.f90

@@ -33,6 +33,8 @@ subroutine write_coord(mol, unit)
    do iat = 1, mol%nat
       write(unit, '(3es24.14, 6x, a)') mol%xyz(:, iat), trim(mol%sym(mol%id(iat)))
    enddo
+   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"

test/test_read_turbomole.f90

@@ -13,6 +13,7 @@
 ! limitations under the License.
 
 module test_read_turbomole
+   use mctc_env, only : wp
    use mctc_env_testing, only : new_unittest, unittest_type, error_type, check
    use mctc_io_read_turbomole
    use mctc_io_structure
@@ -39,6 +40,8 @@ subroutine collect_read_turbomole(testsuite)
       & new_unittest("valid5-coord", test_valid5_coord), &
       & new_unittest("valid6-coord", test_valid6_coord), &
       & new_unittest("valid7-coord", test_valid7_coord), &
+      & new_unittest("valid8-coord", test_valid8_coord), &
+      & new_unittest("valid9-coord", test_valid9_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.), &
@@ -130,6 +133,7 @@ subroutine test_valid2_coord(error)
       "    8.33964228963694     10.16660428027860      3.43902155668011      c", &
       "    8.33965118613331     12.33211762632282      5.02051902430387      c", &
       "$periodic 1", &
+      "$eht charge=0 unpaired=0", &
       "$cell", &
       " 9.29556285275863798006", &
       "$end"
@@ -162,6 +166,7 @@ subroutine test_valid3_coord(error)
       "    0.12856915667443     -0.07403227791901      4.02358027265954  c", &
       "   -0.12317720857511      2.75170732207802     -2.13345350602279  c", &
       "    2.44816466162280      1.28612566399214      4.02317048854901  c", &
+      "$eht unpaired=0 charge=0", &
       "$periodic 2", &
       "$cell  angs", &
       "    2.4809835980     2.4811430162   120.2612191150", &
@@ -318,6 +323,96 @@ subroutine test_valid7_coord(error)
 end subroutine test_valid7_coord
 
 
+subroutine test_valid8_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", &
+      "   -1.79537625851198     -3.77866422935275     -1.07883558363403      h", &
+      "   -2.68278833302782      0.38892666265890      1.66214865238427      s", &
+      "    0.11484649791305      1.48857933226955      3.65660396510375      b", &
+      "   -1.07998879593946     -0.16259121615748     -4.55703065871422      o", &
+      "    0.60302832999383      4.08816149622342     -0.02589373148029      mg", &
+      "   -1.22534089315880     -1.79981382478068     -3.70773173318592      h", &
+      "   -1.33460982049866     -4.24819082475503      2.72791902701083      h", &
+      "   -0.16278082578516      2.41267994179303      5.69030695190570      h", &
+      "    2.87802444057103     -0.33120525058830      1.88311373530297      si", &
+      "    0.68489327931487      0.32790204044961     -4.20547693710673      h", &
+      "   -1.20919773588330     -2.87253762561437      0.94064204223101      b", &
+      "   -3.25572604597922      2.21241092990940     -2.86715549314771      li", &
+      "   -1.83147468262373      5.20527293771933     -2.26976270603341      f", &
+      "    4.90885865772880     -1.92576561961811      2.99069919443735      h", &
+      "    1.26806242248758     -2.60409341782411      0.55162805282247      h", &
+      "    4.11956976339902      1.59892866766766     -1.39117477789609      s", &
+      "$eht unpaired=1", &
+      "$end"
+   rewind(unit)
+
+   call read_coord(struc, unit, error)
+   close(unit)
+   if (allocated(error)) return
+
+   call check(error, struc%nat, 16, "Number of atoms does not match")
+   if (allocated(error)) return
+   call check(error, struc%nid, 8, "Number of species does not match")
+   if (allocated(error)) return
+   call check(error, struc%uhf, 1, "Number of unpaired electrons does not match")
+   if (allocated(error)) return
+
+end subroutine test_valid8_coord
+
+
+subroutine test_valid9_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", &
+      "    4.82824919102333E-02    5.71831000079710E-02    1.73514614763116E-01      C", &
+      "    4.82824919102333E-02    5.71831000079710E-02    2.78568246476372E+00      N", &
+      "    2.46093310136750E+00    5.71831000079710E-02    3.59954953387915E+00      C", &
+      "    3.99138416000780E+00   -2.21116805417472E-01    1.58364683739854E+00      N", &
+      "    2.54075511539052E+00   -1.18599185608072E-01   -5.86344093538442E-01      C", &
+      "   -2.06104824371096E+00    8.28021114689117E-01    4.40357113204146E+00      C", &
+      "    6.72173545596011E+00    2.10496546922931E-01    1.72565972456309E+00      C", &
+      "    3.05878562448454E+00    7.09403031823937E-02    5.55721088395376E+00      H", &
+      "    3.36822820962351E+00   -2.07680855613880E-01   -2.46191575873710E+00      H", &
+      "   -1.68465267663933E+00    1.48551338123814E-01   -9.21486948343917E-01      H", &
+      "   -3.83682349412373E+00    3.78984491295393E-01    3.43261116458953E+00      H", &
+      "   -1.96215889726624E+00   -2.17412943024358E-01    6.19219651728748E+00      H", &
+      "   -1.85966017471395E+00    2.87036107386343E+00    4.74746341688781E+00      H", &
+      "    7.49947096948557E+00   -8.77758695396645E-01    3.31081834253025E+00      H", &
+      "    7.58490546886959E+00   -4.29156708916399E-01   -4.73754235690626E-02      H", &
+      "    7.00829346274163E+00    2.24769645216395E+00    2.03795579552532E+00      H", &
+      "$eht charge=1", &
+      "$end"
+   rewind(unit)
+
+   call read_coord(struc, unit, error)
+   close(unit)
+   if (allocated(error)) return
+
+   call check(error, struc%nat, 16, "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, struc%charge, 1.0_wp, "Total charge does not match")
+   if (allocated(error)) return
+
+end subroutine test_valid9_coord
+
+
 subroutine test_invalid1_coord(error)
 
    !> Error handling