Add support for periodic boundary conditions in QCSchema IO (grimme-lab#93)

See MolSSI/QCSchema#86,
MolSSI/QCSchema#66

Author: awvwgk

.gitignore

@@ -35,3 +35,5 @@
 /build*/
 /_*/
 /docs*/
+/subprojects/*/
+/subprojects/test-drive.wrap

doc/format-qcschema.md

@@ -85,6 +85,24 @@ Caffeine molecule in ``qcschema_molecule`` format.
 }
 ```
 
+## Extensions
+
+The reader supports the following extensions:
+
+- Periodic boundary conditions are specified by providing the lattice vectors in Bohr
+  as extras to the molecule in periodic.lattice as flattened array.
+
+```json
+"extras": {
+  "periodic": {
+    "lattice": [
+       5.5900366437622173, 0.0000000000000000, 0.0000000000000000,
+       0.0000000000000000, 8.6808915904526547, 0.0000000000000000,
+       0.0000000000000000, 0.0000000000000000, 8.6808915904526547
+    ]
+  }
+}
+```
 
 ## Missing features
 

src/mctc/io/read/qcschema.F90

@@ -47,7 +47,7 @@ subroutine read_qcschema(self, unit, error)
 
 #if WITH_JSON
    class(json_value), allocatable :: root
-   type(json_object), pointer :: object, child
+   type(json_object), pointer :: object, child, extras, child2
    type(json_array), pointer :: array, child_array
    type(json_error), allocatable :: parse_error
    type(json_keyval), pointer :: val
@@ -58,8 +58,8 @@ subroutine read_qcschema(self, unit, error)
    character(len=:), allocatable :: symbol, message, schema_name, comment
    character(len=symbol_length), allocatable :: sym(:)
    integer, allocatable :: bond(:, :), list(:)
-   real(wp), allocatable, target :: geo(:)
-   real(wp), pointer :: xyz(:, :)
+   real(wp), allocatable, target :: geo(:), lat(:)
+   real(wp), pointer :: xyz(:, :), lattice(:, :)
 
    call json_load(root, unit, config=json_parser_config(context_detail=1), &
       & context=ctx, error=parse_error)
@@ -222,8 +222,45 @@ subroutine read_qcschema(self, unit, error)
       end do
    end if
 
+   nullify(lattice)
+   nullify(array)
+   call get_value(child, "extras", extras, stat=stat, origin=origin)
+   if (stat /= json_stat%success) then
+      call fatal_error(error, ctx%report("Could not read extras", origin=origin, &
+        & label="Expected object"))
+      return
+   end if
+   call get_value(extras, "periodic", child2, requested=.false., stat=stat, origin=origin)
+   if (stat /= json_stat%success) then
+      call fatal_error(error, ctx%report("Could not read periodic extras", origin=origin, &
+         & label="Expected object with 'lattice' key"))
+      return
+   end if
+   if (associated(child2)) then
+      call get_value(child2, "lattice", array, requested=.false., stat=stat, origin=origin)
+      if (stat /= json_stat%success) then
+         call fatal_error(error, ctx%report("Could not read lattice from extras", origin=origin, &
+            & label="Expected array with 9 elements"))
+         return
+      end if
+   end if
+   if (associated(array)) then
+      call get_value(array, lat, stat=stat, origin=origin)
+      if (stat /= json_stat%success) then
+         call fatal_error(error, ctx%report("Could not read lattice from extras", origin=origin, &
+            & label="Expected array with 9 elements"))
+         return
+      end if
+      if (size(lat) /= 9) then
+         call fatal_error(error, ctx%report("Lattice must have 9 elements", origin=array%origin, &
+            & label="Got "//to_string(size(lat))//" elements"))
+         return
+      end if
+      lattice(1:3, 1:3) => lat(1:9)
+   end if
+
    xyz(1:3, 1:size(geo)/3) => geo
-   call new(self, sym, xyz, charge=real(charge, wp), uhf=multiplicity-1)
+   call new(self, sym, xyz, charge=real(charge, wp), uhf=multiplicity-1, lattice=lattice)
    if (len(comment) > 0) self%comment = comment
    if (allocated(bond)) then
       self%nbd = size(bond, 2)

src/mctc/io/write/qcschema.f90

@@ -95,7 +95,7 @@ pure function json_string(mol, indent) result(string)
 
    string = string // ","
    if (present(indent)) string = string // nl // indent
-   string = string // json_key("geometry", indent) // json_array([mol%xyz], indent)
+   string = string // json_key("geometry", indent) // json_array([mol%xyz], indent, group=3)
 
    string = string // ","
    if (present(indent)) string = string // nl // indent
@@ -129,6 +129,28 @@ pure function json_string(mol, indent) result(string)
       string = string // "]"
    end if
 
+   if (allocated(mol%lattice) .and. any(mol%periodic)) then
+      ! place in "extras": {"periodic": {"lattice": [flattened list...]}
+      string = string // ","
+      if (present(indent)) string = string // nl // indent
+      string = string // json_key("extras", indent) // "{"
+      if (present(indent)) string = string // nl // indent // indent
+      string = string // json_key("periodic", indent) // "{"
+      if (present(indent)) string = string // nl // indent // indent // indent
+      block
+         character(len=:), allocatable :: indent3, indent4
+         if (present(indent)) then
+            indent3 = indent // indent // indent
+            indent4 = indent // indent // indent // indent
+         endif
+         string = string // json_key("lattice", indent) // json_array([mol%lattice], indent3, indent4, group=3)
+      end block
+      if (present(indent)) string = string // nl // indent // indent
+      string = string // "}"
+      if (present(indent)) string = string // nl // indent
+      string = string // "}"
+   end if
+
    if (present(indent)) string = string // nl
    string = string // "}"
 end function json_string
@@ -150,18 +172,29 @@ pure function json_array_int_1(array, indent) result(string)
    string = string // "]"
 end function json_array_int_1
 
-pure function json_array_real_1(array, indent) result(string)
+pure function json_array_real_1(array, indent, indent2, group) result(string)
    real(wp), intent(in) :: array(:)
    character(len=*), intent(in), optional :: indent
+   character(len=*), intent(in), optional :: indent2
+   integer, intent(in), optional :: group
    character(len=:), allocatable :: string
 
-   integer :: i
+   integer :: i, j, step
+
+   step = 1
+   if (present(group)) step = group
 
    string = "["
-   do i = 1, size(array)
-      if (present(indent)) string = string // nl // indent // indent
-      string = string // json_value(array(i), '(es23.16)')
-      if (i /= size(array)) string = string // ","
+   do i = 1, size(array), step
+      if (present(indent2)) then
+         string = string // nl // indent2
+      else if (present(indent)) then
+         string = string // nl // indent // indent
+      end if
+      do j = 1, step, 1
+         string = string // json_value(array(i + j - 1), '(es23.16)')
+         if (i + j - 1 /= size(array)) string = string // ","
+      end do
    end do
    if (present(indent)) string = string // nl // indent
    string = string // "]"

test/CMakeLists.txt

@@ -40,6 +40,7 @@ set(
   "write-genformat"
   "write-pdb"
   "write-qchem"
+  "write-qcschema"
   "write-turbomole"
   "write-vasp"
   "write-xyz"

test/main.f90

@@ -43,6 +43,7 @@ program tester
    use test_write_genformat, only : collect_write_genformat
    use test_write_pdb, only : collect_write_pdb
    use test_write_qchem, only : collect_write_qchem
+   use test_write_qcschema, only : collect_write_qcschema
    use test_write_turbomole, only : collect_write_turbomole
    use test_write_vasp, only : collect_write_vasp
    use test_write_xyz, only : collect_write_xyz
@@ -80,6 +81,7 @@ program tester
       & new_testsuite("write-genformat", collect_write_genformat), &
       & new_testsuite("write-pdb", collect_write_pdb), &
       & new_testsuite("write-qchem", collect_write_qchem), &
+      & new_testsuite("write-qcschema", collect_write_qcschema), &
       & new_testsuite("write-turbomole", collect_write_turbomole), &
       & new_testsuite("write-vasp", collect_write_vasp), &
       & new_testsuite("write-xyz", collect_write_xyz) &

test/meson.build

@@ -38,6 +38,7 @@ tests = [
   'write-genformat',
   'write-pdb',
   'write-qchem',
+  'write-qcschema',
   'write-turbomole',
   'write-vasp',
   'write-xyz',