8000 Update Turbomole reader by awvwgk · Pull Request #18 · grimme-lab/mctc-lib · GitHub
[go: up one dir, main page]
More Web Proxy on the site http://driver.im/
Skip to content

Update Turbomole reader #18

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

Merged
merged 1 commit into from
Aug 5, 2021
Merged

Update Turbomole reader #18

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
4 changes: 3 additions & 1 deletion doc/format-tmol.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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.
39 changes: 26 additions & 13 deletions src/mctc/io/read/turbomole.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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.
Expand All @@ -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.
Expand All @@ -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.
Expand All @@ -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
Expand Down Expand Up @@ -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

Expand Down
2 changes: 2 additions & 0 deletions src/mctc/io/write/turbomole.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand Down
95 changes: 95 additions & 0 deletions test/test_read_turbomole.f90
View file
Open in desktop
6D47
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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.), &
Expand Down Expand Up @@ -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"
Expand Down Expand Up @@ -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", &
Expand Down Expand Up @@ -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
Expand Down
0