Switch to V3K connection table format for 1000+ atoms or bonds (grimme-lab#107)

Co-authored-by: awvwgk <28669218+awvwgk@users.noreply.github.com>

Author: Copilot

src/mctc/io/write/ctfile.f90

@@ -21,6 +21,9 @@ module mctc_io_write_ctfile
 
    public :: write_molfile, write_sdf
 
+   !> Maximum number of atoms/bonds supported by V2000 format
+   integer, parameter :: v2k_max = 999
+
 
 contains
 
@@ -57,7 +60,27 @@ subroutine write_sdf(self, unit, energy, gnorm)
 end subroutine write_sdf
 
 
+!> Write molecular structure data to a molfile
 subroutine write_molfile(self, unit, comment_line)
+   !> Instance of the molecular structure data
+   class(structure_type), intent(in) :: self
+   !> File handle
+   integer, intent(in) :: unit
+   !> Optional comment line
+   character(len=*), intent(in), optional :: comment_line
+
+   ! Switch to V3000 format if we exceed V2000 limits
+   if (self%nat > v2k_max .or. self%nbd > v2k_max) then
+      call write_molfile_v3k(self, unit, comment_line)
+   else
+      call write_molfile_v2k(self, unit, comment_line)
+   end if
+
+end subroutine write_molfile
+
+
+!> Write molecular structure data to a V2000 molfile
+subroutine write_molfile_v2k(self, unit, comment_line)
    class(structure_type), intent(in) :: self
    integer, intent(in) :: unit
    character(len=*), intent(in), optional :: comment_line
@@ -129,7 +152,88 @@ subroutine write_molfile(self, unit, comment_line)
 
    write(unit, '(a)') "M  END"
 
-end subroutine write_molfile
+end subroutine write_molfile_v2k
+
+
+!> Write molecular structure data to a V3000 molfile
+subroutine write_molfile_v3k(self, unit, comment_line)
+   class(structure_type), intent(in) :: self
+   integer, intent(in) :: unit
+   character(len=*), intent(in), optional :: comment_line
+   integer :: iatom, ibond, btype
+   logical :: has_sdf_data
+   character(len=8)  :: date
+   character(len=10) :: time
+   character(len=256) :: line
+
+   call date_and_time(date, time)
+
+   ! Header block (3 lines)
+   if (present(comment_line)) then
+      write(unit, '(a)') comment_line
+   else
+      if (allocated(self%comment)) then
+         write(unit, '(a)') self%comment
+      else
+         write(unit, '(a)')
+      end if
+   end if
+   write(unit, '(2x, 3x, 5x, 3a2, a4, "3D")') &
+      &  date(5:6), date(7:8), date(3:4), time(:4)
+   write(unit, '(a)')
+
+   ! Counts line for V3000 (atoms and bonds are set to 0 in header)
+   write(unit, '(a)') "  0  0  0     0  0            999 V3000"
+
+   ! V3000 block
+   write(unit, '(a)') "M  V30 BEGIN CTAB"
+   write(unit, '(a, i0, a, i0, a)') "M  V30 COUNTS ", self%nat, " ", self%nbd, " 0 0 0"
+
+   has_sdf_data = allocated(self%sdf)
+
+   ! Atom block
+   write(unit, '(a)') "M  V30 BEGIN ATOM"
+   do iatom = 1, self%nat
+      ! Build basic atom line: index, symbol, x, y, z, aamap (0)
+      write(line, '(a, i0, a, a, 3(a, f0.6), a)') &
+         & "M  V30 ", iatom, " ", trim(self%sym(self%id(iatom))), &
+         & " ", self%xyz(1, iatom)*autoaa, &
+         & " ", self%xyz(2, iatom)*autoaa, &
+         & " ", self%xyz(3, iatom)*autoaa, " 0"
+
+      ! Add optional properties
+      if (has_sdf_data) then
+         if (self%sdf(iatom)%charge /= 0) then
+            write(line, '(a, a, i0)') trim(line), " CHG=", self%sdf(iatom)%charge
+         end if
+         if (self%sdf(iatom)%valence /= 0) then
+            write(line, '(a, a, i0)') trim(line), " VAL=", self%sdf(iatom)%valence
+         end if
+      end if
+
+      write(unit, '(a)') trim(line)
+   end do
+   write(unit, '(a)') "M  V30 END ATOM"
+
+   ! Bond block
+   if (self%nbd > 0) then
+      write(unit, '(a)') "M  V30 BEGIN BOND"
+      do ibond = 1, self%nbd
+         if (size(self%bond, 1) > 2) then
+            btype = self%bond(3, ibond)
+         else
+            btype = 1
+         end if
+         write(unit, '(a, i0, a, i0, a, i0, a, i0)') &
+            & "M  V30 ", ibond, " ", btype, " ", self%bond(1, ibond), " ", self%bond(2, ibond)
+      end do
+      write(unit, '(a)') "M  V30 END BOND"
+   end if
+
+   write(unit, '(a)') "M  V30 END CTAB"
+   write(unit, '(a)') "M  END"
+
+end subroutine write_molfile_v3k
 
 
 end module mctc_io_write_ctfile

subprojects/.wraplock

@@ -13,6 +13,7 @@
 ! limitations under the License.
 
 module test_write_ctfile
+   use mctc_env_accuracy, only : wp
    use mctc_env_testing, only : new_unittest, unittest_type, error_type, check
    use testsuite_structure, only : get_structure
    use mctc_io_write_ctfile
@@ -35,7 +36,10 @@ subroutine collect_write_ctfile(testsuite)
 
    testsuite = [ &
       & new_unittest("valid1-mol", test_valid1_mol), &
-      & new_unittest("valid1-sdf", test_valid1_sdf) &
+      & new_unittest("valid1-sdf", test_valid1_sdf), &
+      & new_unittest("v3k-large-mol", test_v3k_large_mol), &
+      & new_unittest("v3k-large-sdf", test_v3k_large_sdf), &
+      & new_unittest("v3k-with-bonds", test_v3k_with_bonds) &
       & ]
 
 end subroutine collect_write_ctfile
@@ -97,4 +101,147 @@ subroutine test_valid1_sdf(error)
 end subroutine test_valid1_sdf
 
 
+subroutine test_v3k_large_mol(error)
+
+   !> Error handling
+   type(error_type), allocatable, intent(out) :: error
+
+   type(structure_type) :: struc, struc_read
+   integer :: unit, nat, i
+   character(len=2), allocatable :: sym(:)
+   real(wp), allocatable :: xyz(:, :)
+
+   ! Create a structure with 1001 atoms to trigger V3000 format
+   nat = 1001
+   allocate(sym(nat), xyz(3, nat))
+
+   ! Create alternating H and C atoms on a line
+   do i = 1, nat
+      if (mod(i, 2) == 0) then
+         sym(i) = "C"
+      else
+         sym(i) = "H"
+      end if
+      xyz(1, i) = real(i, wp)
+      xyz(2, i) = 0.0_wp
+      xyz(3, i) = 0.0_wp
+   end do
+
+   call new(struc, sym, xyz)
+
+   open(status='scratch', newunit=unit)
+   call write_molfile(struc, unit)
+   rewind(unit)
+
+   call read_molfile(struc_read, unit, error)
+   close(unit)
+   if (allocated(error)) return
+
+   call check(error, struc_read%nat, nat, "Number of atoms does not match")
+   if (allocated(error)) return
+   call check(error, struc_read%nid, 2, "Number of species does not match")
+   if (allocated(error)) return
+
+end subroutine test_v3k_large_mol
+
+
+subroutine test_v3k_large_sdf(error)
+
+   !> Error handling
+   type(error_type), allocatable, intent(out) :: error
+
+   type(structure_type) :: struc, struc_read
+   integer :: unit, nat, i
+   character(len=2), allocatable :: sym(:)
+   real(wp), allocatable :: xyz(:, :)
+
+   ! Create a structure with 1001 atoms to trigger V3000 format
+   nat = 1001
+   allocate(sym(nat), xyz(3, nat))
+
+   ! Create alternating H and C atoms on a line
+   do i = 1, nat
+      if (mod(i, 2) == 0) then
+         sym(i) = "C"
+      else
+         sym(i) = "H"
+      end if
+      xyz(1, i) = real(i, wp)
+      xyz(2, i) = 0.0_wp
+      xyz(3, i) = 0.0_wp
+   end do
+
+   call new(struc, sym, xyz)
+
+   open(status='scratch', newunit=unit)
+   call write_sdf(struc, unit)
+   rewind(unit)
+
+   call read_sdf(struc_read, unit, error)
+   close(unit)
+   if (allocated(error)) return
+
+   call check(error, struc_read%nat, nat, "Number of atoms does not match")
+   if (allocated(error)) return
+   call check(error, struc_read%nid, 2, "Number of species does not match")
+   if (allocated(error)) return
+
+end subroutine test_v3k_large_sdf
+
+
+subroutine test_v3k_with_bonds(error)
+
+   !> Error handling
+   type(error_type), allocatable, intent(out) :: error
+
+   type(structure_type) :: struc, struc_read
+   integer :: unit, nat, nbd, i
+   character(len=2), allocatable :: sym(:)
+   real(wp), allocatable :: xyz(:, :)
+   integer, allocatable :: bond(:, :)
+
+   ! Create a structure with 1001 atoms and 1000 bonds to trigger V3000 format
+   nat = 1001
+   nbd = 1000
+   allocate(sym(nat), xyz(3, nat), bond(3, nbd))
+
+   ! Create a linear chain of alternating H and C atoms
+   do i = 1, nat
+      if (mod(i, 2) == 0) then
+         sym(i) = "C"
+      else
+         sym(i) = "H"
+      end if
+      xyz(1, i) = real(i, wp)
+      xyz(2, i) = 0.0_wp
+      xyz(3, i) = 0.0_wp
+   end do
+
+   ! Create bonds connecting each consecutive atom pair
+   do i = 1, nbd
+      bond(1, i) = i
+      bond(2, i) = i + 1
+      bond(3, i) = 1  ! single bond
+   end do
+
+   call new(struc, sym, xyz, bond=bond)
+
+   open(status='scratch', newunit=unit)
+   call write_molfile(struc, unit)
+   rewind(unit)
+
+   call read_molfile(struc_read, unit, error)
+   close(unit)
+   if (allocated(error)) return
+
+   call check(error, struc_read%nat, nat, "Number of atoms does not match")
+   if (allocated(error)) return
+   call check(error, struc_read%nid, 2, "Number of species does not match")
+   if (allocated(error)) return
+   call check(error, struc_read%nbd, nbd, "Number of bonds does not match")
+   if (allocated(error)) return
+
+end subroutine test_v3k_with_bonds
+
+
 end module test_write_ctfile