Support initial charge and moment for FHI-aims I/O

Author: Copilot

doc/format-aims.md

@@ -28,6 +28,10 @@ The format is identified by:
 | `atom` | Cartesian coordinates in Ångström |
 | `atom_frac` | Fractional coordinates (periodic systems) |
 | `lattice_vector` | Lattice vector (3 reals, one per line) |
+| `initial_charge` | Global initial total molecular charge |
+| `initial_moment` | Global initial spin moment / unpaired-electron count |
+
+The reader and writer support the optional global `initial_charge` and `initial_moment` keywords for preserving molecular charge and spin information in `geometry.in` files.
 
 ### Atom Specification
 

src/mctc/io/read/aims.f90

@@ -20,7 +20,7 @@ module mctc_io_read_aims
    use mctc_io_symbols, only : symbol_length, to_number
    use mctc_io_structure, only : structure_type, new
    use mctc_io_utils, only : next_line, token_type, next_token, io_error, filename, &
-      read_next_token, to_string
+      read_next_token, to_lower, to_string
    implicit none
    private
 
@@ -44,7 +44,7 @@ subroutine read_aims(mol, unit, error)
    integer :: stat, pos, lnum, ilt, iat
    type(token_type) :: token
    character(len=:), allocatable :: line
-   real(wp) :: x, y, z
+   real(wp) :: x, y, z, charge, moment
    character(len=symbol_length), allocatable :: sym(:)
    real(wp), allocatable :: xyz(:, :), abc(:, :), lattice(:, :)
    logical :: is_frac, periodic(3)
@@ -58,6 +58,8 @@ subroutine read_aims(mol, unit, error)
    iat = 0
    ilt = 0
    periodic(:) = .false.
+   charge = 0.0_wp
+   moment = 0.0_wp
 
    lnum = 0
    stat = 0
@@ -68,7 +70,7 @@ subroutine read_aims(mol, unit, error)
       if (line(1:1) == "#") cycle
 
       call next_token(line, pos, token)
-      select case(line(token%first:token%last))
+      select case(to_lower(line(token%first:token%last)))
       case("atom", "atom_frac")
          is_frac = token%last - token%first + 1 > 4
          call read_next_token(line, pos, token, x, stat)
@@ -127,6 +129,22 @@ subroutine read_aims(mol, unit, error)
          if (.not.allocated(lattice)) allocate(lattice(3, 3), source=0.0_wp)
          lattice(:, ilt) = [x, y, z] * aatoau
 
+      case("initial_charge")
+         call read_next_token(line, pos, token, charge, stat)
+         if (stat /= 0) then
+            call io_error(error, "Cannot read initial charge", &
+               & line, token, filename(unit), lnum, "expected real value")
+            exit
+         end if
+
+      case("initial_moment")
+         call read_next_token(line, pos, token, moment, stat)
+         if (stat /= 0) then
+            call io_error(error, "Cannot read initial moment", &
+               & line, token, filename(unit), lnum, "expected real value")
+            exit
+         end if
+
       case default
          call io_error(error, "Unexpected keyword found", &
             & line, token, filename(unit), lnum, "invalid in this context")
@@ -149,7 +167,8 @@ subroutine read_aims(mol, unit, error)
       periodic(:ilt) = .true.
    end if
 
-   call new(mol, sym(:iat), xyz, lattice=lattice, periodic=periodic)
+   call new(mol, sym(:iat), xyz, charge=charge, uhf=nint(moment), &
+      & lattice=lattice, periodic=periodic)
 
 end subroutine read_aims
 

src/mctc/io/write/aims.f90

@@ -34,6 +34,13 @@ subroutine write_aims(self, unit)
    integer :: iat, ilt
    logical :: expo
 
+   if (self%charge /= 0.0_wp) then
+      write(unit, '(a, 1x, es24.14)') "initial_charge", self%charge
+   end if
+   if (self%uhf /= 0) then
+      write(unit, '(a, 1x, es24.14)') "initial_moment", real(self%uhf, wp)
+   end if
+
    expo = maxval(self%xyz) > 1.0e+5 .or. minval(self%xyz) < -1.0e+5
    if (expo) then
       do iat = 1, self%nat

test/test_read_aims.f90

@@ -40,6 +40,8 @@ subroutine collect_read_aims(testsuite)
       & new_unittest("valid5-aims", test_valid5_aims), &
       & new_unittest("valid6-aims", test_valid6_aims), &
       & new_unittest("valid7-aims", test_valid7_aims), &
+      & new_unittest("valid8-aims", test_valid8_aims), &
+      & new_unittest("valid9-aims", test_valid9_aims), &
       & new_unittest("invalid1-aims", test_invalid1_aims, should_fail=.true.), &
       & new_unittest("invalid2-aims", test_invalid2_aims, should_fail=.true.), &
       & new_unittest("invalid3-aims", test_invalid3_aims, should_fail=.true.), &
@@ -355,6 +357,61 @@ end subroutine test_valid7_aims
 
 
 
+subroutine test_valid8_aims(error)
+
+   !> Error handling
+   type(error_type), allocatable, intent(out) :: error
+
+   type(structure_type) :: struc
+   integer :: unit
+
+   open(status='scratch', newunit=unit)
+   write(unit, '(a)') &
+      "initial_charge 1.0", &
+      "initial_moment 2.0", &
+      "atom 0.0 0.0 0.0 H", &
+      "atom 0.0 0.0 0.0 He"
+   rewind(unit)
+
+   call read_aims(struc, unit, error)
+   close(unit)
+   if (allocated(error)) return
+
+   call check(error, struc%charge, 1.0_wp, "Charge does not match")
+   if (allocated(error)) return
+   call check(error, struc%uhf, 2, "Spin information does not match")
+   if (allocated(error)) return
+
+end subroutine test_valid8_aims
+
+
+subroutine test_valid9_aims(error)
+
+   !> Error handling
+   type(error_type), allocatable, intent(out) :: error
+
+   type(structure_type) :: struc
+   integer :: unit
+
+   open(status='scratch', newunit=unit)
+   write(unit, '(a)') &
+      "initial_charge -1.0", &
+      "initial_moment -2.0", &
+      "atom 0.0 0.0 0.0 H"
+   rewind(unit)
+
+   call read_aims(struc, unit, error)
+   close(unit)
+   if (allocated(error)) return
+
+   call check(error, struc%charge, -1.0_wp, "Charge does not match")
+   if (allocated(error)) return
+   call check(error, struc%uhf, -2, "Spin information does not match")
+   if (allocated(error)) return
+
+end subroutine test_valid9_aims
+
+
 subroutine test_invalid1_aims(error)
 
    !> Error handling

test/test_write_aims.f90

@@ -13,6 +13,7 @@
 ! limitations under the License.
 
 module test_write_aims
+   use mctc_env, only : wp
    use mctc_env_testing, only : new_unittest, unittest_type, error_type, check
    use testsuite_structure, only : get_structure
    use mctc_io_write_aims
@@ -35,7 +36,8 @@ subroutine collect_write_aims(testsuite)
 
    testsuite = [ &
       & new_unittest("valid1-aims", test_valid1_aims), &
-      & new_unittest("valid2-aims", test_valid2_aims) &
+      & new_unittest("valid2-aims", test_valid2_aims), &
+      & new_unittest("valid3-aims", test_valid3_aims) &
       & ]
 
 end subroutine collect_write_aims
@@ -97,4 +99,31 @@ subroutine test_valid2_aims(error)
 end subroutine test_valid2_aims
 
 
+subroutine test_valid3_aims(error)
+
+   !> Error handling
+   type(error_type), allocatable, intent(out) :: error
+
+   type(structure_type) :: struc
+   integer :: unit
+
+   call new(struc, [6, 1], ["C", "H"], reshape([0.0_wp, 0.0_wp, 0.0_wp, 0.5_wp, 0.0_wp, 0.0_wp], [3, 2]), &
+      & charge=1.0_wp, uhf=2)
+
+   open(status='scratch', newunit=unit)
+   call write_aims(struc, unit)
+   rewind(unit)
+
+   call read_aims(struc, unit, error)
+   close(unit)
+   if (allocated(error)) return
+
+   call check(error, struc%charge, 1.0_wp, "Charge does not match")
+   if (allocated(error)) return
+   call check(error, struc%uhf, 2, "Spin information does not match")
+   if (allocated(error)) return
+
+end subroutine test_valid3_aims
+
+
 end module test_write_aims