Add GFN-FF charge initialization options (#315)

CMakeLists.txt

@@ -27,7 +27,7 @@ endif()
 project(
   crest
   LANGUAGES "C" "Fortran"
-  VERSION  3.0.1 
+  VERSION  3.0.2
   DESCRIPTION "A tool for the exploration of low-energy chemical space"
 )
 

meson.build

@@ -17,7 +17,7 @@
 project(
   'crest',
   'fortran', 'c',
-  version: '3.0.1',
+  version: '3.0.2',
   license: 'LGPL-3.0-or-later',
   meson_version: '>=0.63',
   default_options: [

src/calculator/api_helpers.F90

@@ -325,6 +325,9 @@ subroutine gfnff_init(calc,loadnew)
       if (allocated(calc%parametrisation)) then
         calc%ff_dat%refgeo = calc%parametrisation
       end if
+      if (allocated(calc%refcharges)) then
+        calc%ff_dat%refcharges = calc%refcharges
+      end if
 
     end if
     if (allocated(calc%solvent)) then

src/calculator/calc_type.f90

@@ -112,6 +112,7 @@ module calc_type
 
     !> atomic charges
     logical :: rdqat = .false.
+    logical :: dumpq = .false.
     real(wp),allocatable :: qat(:)
 
     !> dipole and dipole gradient
@@ -137,6 +138,7 @@ module calc_type
     logical  :: restart = .false.  !> restart option (some potentials can do this)
     character(len=:),allocatable :: restartfile
     character(len=:),allocatable :: refgeo
+    character(len=:),allocatable :: refcharges
 
 !>--- tblite data
     type(tblite_data),allocatable :: tblite

src/calculator/printouts.F90

@@ -106,6 +106,9 @@ subroutine calculation_summary(calc,mol,energy,grad,molnew,iounit,print)
             & calc%calcs(k)%qat(j), cn(j)
           enddo  
           write (iunit,*)
+          if(calc%calcs(k)%dumpq)then
+            call dumpq(k,mol%nat,calc%calcs(k)%qat)
+          endif
         end if
       end do
       write (iunit,'(a)') repeat('-',80)
@@ -214,6 +217,26 @@ subroutine calculation_summary(calc,mol,energy,grad,molnew,iounit,print)
 
   end subroutine calculation_summary
 
+!========================================================================================!
+
+  subroutine dumpq(id,nat,q)
+!********************************
+!* write atomic charges to file
+!********************************
+    implicit none
+    integer,intent(in) :: id,nat
+    real(wp),intent(in) :: q(nat)
+    integer :: i,ich
+    character(len=50) :: atmp
+
+    write(atmp,'("charges.",i0)') id
+    open(newunit=ich,file=trim(atmp))
+    do i=1,nat
+       write(ich,'(F20.10)') q(i)
+    enddo
+    close(ich)
+  end subroutine dumpq 
+
 !========================================================================================!
 !========================================================================================!
 end module crest_calculator_printout

src/parsing/parse_calcdata.f90

@@ -273,6 +273,8 @@ subroutine parse_setting_c(env,job,key,val)
         job%id = jobtype%tblite
         job%tblitelvl = xtblvl%ceh
         job%rdgrad = .false.
+        job%rdqat = .true.
+        job%rddip = .true.
       case ('gfn0','gfn0-xtb')
         job%id = jobtype%gfn0
       case ('gfn0*','gfn0*-xtb')
@@ -396,6 +398,15 @@ subroutine parse_setting_c(env,job,key,val)
         write(stderr,'(a,a,a)') 'specified parametrisation file ',val,' does not exist'
         error stop
       endif
+    case('refchrg','refcharges')
+      inquire(file=val,exist=ex)
+      if(ex)then
+        job%refcharges = val
+      else
+        write(stderr,'(a,a,a)') 'specified reference charge file ',val,' does not exist'
+        error stop
+      endif
+
 
     case ('print')
       select case (val)
@@ -427,6 +438,9 @@ subroutine parse_setting_bool(job,key,val)
       job%rddip = val
     case ('rdqat','rdchrg')
       job%rdqat = val
+    case ('dumpq','dumpchrg')
+      job%rdqat = val
+      job%dumpq = val
     case ('dipgrad')
       job%rddipgrad = val
     case ('rdgrad')

src/parsing/parse_xtbinput.f90

@@ -719,13 +719,15 @@ subroutine get_xtb_keyvalue(kv,str,io)
     character(len=*) :: str
     integer,intent(out) :: io
     character(len=:),allocatable :: tmpstr
+    character(len=:),allocatable :: tmpstr_rc
     character(len=:),allocatable :: ktmp
     character(len=:),allocatable :: vtmp
     integer :: i,j,k,na,plast
     integer :: l(3)
     call kv%deallocate()
     io = 0
     tmpstr = adjustl(lowercase(str))
+    tmpstr_rc=adjustl(trim(str))
 
     !> key-value conditions
     l(1) = index(tmpstr,'=')
@@ -747,7 +749,7 @@ subroutine get_xtb_keyvalue(kv,str,io)
     end if
 
     ktmp = trim(adjustl(tmpstr(:k-1)))
-    vtmp = trim(adjustl(tmpstr(k+1:)))
+    vtmp = trim(adjustl(tmpstr_rc(k+1:)))
     kv%key = ktmp !> the key as string
     kv%rawvalue = vtmp !> value as unformatted string