Allows writing of Turbomole format with angstrom.

src/mctc/io/write/turbomole.f90

@@ -15,30 +15,42 @@
 module mctc_io_write_turbomole
    use mctc_env_accuracy, only : wp
    use mctc_io_structure, only : structure_type
+   use mctc_io_convert, only : autoaa
    implicit none
    private
 
    public :: write_coord
 
-
 contains
 
-
 subroutine write_coord(mol, unit)
    class(structure_type), intent(in) :: mol
    integer, intent(in) :: unit
    integer :: iat, ilt, npbc
    logical :: expo
 
-   write(unit, '(a)') "$coord"
+   real(wp) :: conv_fac
+   logical :: angs
+
+   angs = mol%info%angs_coord
+   conv_fac = 1.0_wp
+   if (angs) conv_fac = autoaa
+
+   if (angs) then
+      write(unit, '(a)') "$coord angs"
+   else
+      write(unit, '(a)') "$coord"
+   end if
    expo = maxval(mol%xyz) > 1.0e+5 .or. minval(mol%xyz) < -1.0e+5
    if (expo) then
       do iat = 1, mol%nat
-         write(unit, '(3es24.14, 6x, a)') mol%xyz(:, iat), trim(mol%sym(mol%id(iat)))
+         write(unit, '(3es24.14, 6x, a)') mol%xyz(:, iat) * conv_fac, &
+            trim(mol%sym(mol%id(iat)))
       end do
    else
       do iat = 1, mol%nat
-         write(unit, '(3f24.14, 6x, a)') mol%xyz(:, iat), trim(mol%sym(mol%id(iat)))
+         write(unit, '(3f24.14, 6x, a)') mol%xyz(:, iat) * conv_fac, &
+            trim(mol%sym(mol%id(iat)))
       end do
    end if
    if (any([nint(mol%charge), mol%uhf] /= 0)) then
@@ -49,15 +61,18 @@ subroutine write_coord(mol, unit)
       write(unit, '(a, 1x, i0)') "$periodic", count(mol%periodic)
       npbc = count(mol%periodic)
       if (size(mol%lattice, 2) == 3) then
-         write(unit, '(a)') "$lattice bohr"
+         if (angs) then
+            write(unit, '(a)') "$lattice angs"
+         else
+            write(unit, '(a)') "$lattice bohr"
+         end if
          do ilt = 1, npbc
-            write(unit, '(3f20.14)') mol%lattice(:npbc, ilt)
+            write(unit, '(3f20.14)') mol%lattice(:npbc, ilt) * conv_fac
          end do
       end if
    end if
    write(unit, '(a)') "$end"
 
 end subroutine write_coord
 
-
-end module mctc_io_write_turbomole
+end module mctc_io_write_turbomole

test/test_write_turbomole.f90

@@ -18,6 +18,7 @@ module test_write_turbomole
    use mctc_io_write_turbomole
    use mctc_io_read_turbomole
    use mctc_io_structure
+   use mctc_io_convert, only : autoaa
    implicit none
    private
 
@@ -35,7 +36,9 @@ subroutine collect_write_turbomole(testsuite)
 
    testsuite = [ &
       & new_unittest("valid1-coord", test_valid1_coord), &
-      & new_unittest("valid2-coord", test_valid2_coord) &
+      & new_unittest("valid2-coord", test_valid2_coord), &
+      & new_unittest("valid1-coord-angs", test_valid1_coord_angs), &
+      & new_unittest("valid2-coord-angs", test_valid2_coord_angs) &
       & ]
 
 end subroutine collect_write_turbomole
@@ -96,5 +99,69 @@ subroutine test_valid2_coord(error)
 
 end subroutine test_valid2_coord
 
+subroutine test_valid1_coord_angs(error)
+
+   !> Error handling
+   type(error_type), allocatable, intent(out) :: error
+
+   type(structure_type) :: struc
+   integer :: unit, nat, nid
+
+   call get_structure(struc, "mindless01")
+   nat = struc%nat
+   nid = struc%nid
+
+   struc%xyz = struc%xyz * autoaa
+   struc%info%angs_coord = .true.
+
+   open(status='scratch', newunit=unit)
+   call write_coord(struc, unit)
+   rewind(unit)
+
+   call read_coord(struc, unit, error)
+   close(unit)
+   if (allocated(error)) return
+
+   call check(error, struc%nat, nat, "Number of atoms does not match")
+   if (allocated(error)) return
+   call check(error, struc%nid, nid, "Number of species does not match")
+   if (allocated(error)) return
+   call check(error, struc%info%angs_coord, .true., "Coordinates are not written in angstrom.")
+   if (allocated(error)) return
+
+end subroutine test_valid1_coord_angs
+
+subroutine test_valid2_coord_angs(error)
+
+   !> Error handling
+   type(error_type), allocatable, intent(out) :: error
+
+   type(structure_type) :: struc
+   integer :: unit, nat, nid
+
+   call get_structure(struc, "x01")
+   nat = struc%nat
+   nid = struc%nid
+
+   struc%xyz = struc%xyz * autoaa
+   struc%lattice = struc%lattice * autoaa
+   struc%info%angs_coord = .true.
+
+   open(status='scratch', newunit=unit)
+   call write_coord(struc, unit)
+   rewind(unit)
+
+   call read_coord(struc, unit, error)
+   close(unit)
+   if (allocated(error)) return
+
+   call check(error, struc%nat, nat, "Number of atoms does not match")
+   if (allocated(error)) return
+   call check(error, struc%nid, nid, "Number of species does not match")
+   if (allocated(error)) return
+   call check(error, struc%info%angs_coord, .true., "Coordinates are not written in angstrom.")
+   if (allocated(error)) return
+
+end subroutine test_valid2_coord_angs
 
 end module test_write_turbomole