Allow QCG site directed solvation (#226)

* Make final GFN2 optimizations in QCG optional

* Small QCG Bugfix

* Allow solvation at defined position via directed docking

* Updating the help printout

---------

Signed-off-by: cplett <[email protected]>
Co-authored-by: Philipp Pracht <[email protected]>

src/classes.f90

@@ -375,6 +375,9 @@ module crest_data
       integer           :: nsolv = 0            !Number of solventmolecules
       integer           :: nqcgclust = 0        !Number of cluster to be taken
       integer           :: max_solv = 0         !Maximal number of solvents added, if none is given
+      character(len=:), allocatable :: directed_file !name of the directed list
+      character(len=64), allocatable :: directed_list(:,:) !How many solvents at which atom to add
+      integer, allocatable :: directed_number(:) !Numbers of solvents added per defined atom
       integer           :: ensemble_method = -1 !Default -1 for qcgmtd, 0= crest, 1= standard MD, 2= MTD
       character(len=20) :: ensemble_opt         !Method for ensemble optimization in qcg mode
       character(len=20) :: freqver              !Method for frequency computation in qcg mode

src/confparse.f90

@@ -1646,6 +1646,13 @@ subroutine parseflags(env,arg,nra)
         env%docking_qcg_flag = ''
       case ('-fin_opt_gfn2')
         env%final_gfn2_opt = .true.
+      case ('-directed') !specify the directed list
+        env%qcg_flag = .true.
+        ctmp = trim(arg(i + 1))
+        if (ctmp(1:1) .ne. '-') then
+          env%directed_file = trim(ctmp)
+          write (*,'(2x,a,1x,a)') trim(argument)//' :',trim(ctmp)
+        end if
       case ('-nclus')
         env%qcg_flag = .true.
         call readl(arg(i + 1),xx,j)

src/printouts.f90

@@ -234,6 +234,7 @@ subroutine confscript_help()
       write(*,'(5x,''-wscal <FLOAT>     : Scaling factor for outer wall potential'')')
       write(*,'(5x,''-samerand          : use same random number for every xtbiff run'')')
       write(*,'(5x,''-fin_opt_gfn2      : perform GFN2-xTB optimizations for final grow and ensemble structures'')')
+      write(*,'(5x,''-directed <FILE>   : Perform directed solvation at positions defined in <FILE>'')')
       write(*,'(5x,''-ensemble          : ensemble generation'')')
       write(*,'(5x,''-qcgmtd            : NCI-MTD CREST ensemble generation (Default)'')')
       write(*,'(5x,''-ncimtd            : NCI-MTD CREST ensemble generation'')')

src/qcg/solvtool.f90

@@ -170,7 +170,9 @@ subroutine qcg_setup(env, solu, solv)
    logical :: e_there, tmp, used_tmp
    character(len=512) :: thispath, tmp_grow
    character(len=40)  :: solv_tmp
-   character(len=80) :: atmp
+   character(len=80)  :: atmp
+   character(len=20)  :: gfnver_tmp
+
 
    call getcwd(thispath)
 
@@ -224,6 +226,11 @@ subroutine qcg_setup(env, solu, solv)
    env%ref%xyz = solu%xyz
 
 !---- Geometry preoptimization solute
+   if (env%final_gfn2_opt) then !If GFN2 final opt, solute also GFN2 optimized
+      gfnver_tmp = env%gfnver
+      env%gfnver = '--gfn2'
+   end if
+
    if ((.not. env%nopreopt) .and. (solu%nat /= 1)) then
       call wrc0('coord', solu%nat, solu%at, solu%xyz) !write coord for xtbopt routine
       if (env%cts%used) then
@@ -263,6 +270,10 @@ subroutine qcg_setup(env, solu, solv)
       call xtbsp3(env, 'solute')
    end if
 
+   if (env%final_gfn2_opt) then !If GFN2 final opt, solute also GFN2 optimized
+      env%gfnver = gfnver_tmp
+   end if
+
    call grepval('xtb.out', '| TOTAL ENERGY', e_there, solu%energy)
    if (.not. e_there) then
       write (*, *) 'Total Energy of solute not found'
@@ -358,9 +369,9 @@ subroutine read_qcg_input(env, solu, solv)
    type(systemdata)               :: env
    type(zmolecule), intent(inout) :: solu, solv
    logical                        :: pr
-   real(wp), parameter             :: amutokg = 1.66053886E-27
-   real(wp), parameter             :: third = 1.0d0/3.0d0
-   integer                        :: i
+   real(wp), parameter            :: amutokg = 1.66053886E-27
+   real(wp), parameter            :: third = 1.0d0/3.0d0
+   integer                        :: i, ich
    real(wp)                       :: r_solu, r_solv
 
    pr = .true.
@@ -403,8 +414,89 @@ subroutine read_qcg_input(env, solu, solv)
    solu%mass = solu%mass*amutokg
    solv%mass = solv%mass*amutokg
 
+!--- If directed docking is requested, it is read in here:
+   if(allocated(env%directed_file)) then
+      if (env%use_xtbiff) error stop 'xTB-IFF does not support directed docking. &
+              &Please use the aISS algorithm of xtb.'
+      call read_directed_input(env)
+   end if
+
+
 end subroutine read_qcg_input
 
+!> Read input for directed docking
+subroutine read_directed_input(env)
+   use iso_fortran_env, wp => real64
+   use crest_data
+   implicit none
+
+   type(systemdata)           :: env
+
+   integer                    :: nlines
+   integer                    :: io, ich, i, i_check
+   integer                    :: index
+   character(len=512)         :: dum
+   character(len=1), parameter :: delim_space = ' ', delim_tab = achar(9)
+
+     open (newunit=ich, file=env%directed_file) 
+        !First check number of lines
+        nlines = 0
+        do
+           read(ich,*,iostat=io)
+           if (io /= 0) exit
+           nlines = nlines + 1
+        end do
+        !Allocate directed list 
+        !First entry is the atom number, Second how many solvents to add to this atom
+        allocate(env%directed_list(nlines,2))
+        allocate(env%directed_number(nlines), source = 0)
+        !Now read lines into directed_list
+        rewind(ich)
+        do i=1, nlines
+           read(ich,'(A)') dum
+           !> Remove leading tab and spaces first
+           dum = adjustl(dum) !Leading spaces are removed
+           index = SCAN(trim(dum), delim_tab) 
+           if (index == 1) then !Leading tab -> remove it
+              dum = dum(2:)
+           end if
+           index = SCAN(trim(dum), delim_space)
+           if (index == 0) then !No space = check for tab
+              index = SCAN(trim(dum), delim_tab)
+           end if
+           if (index == 0) then  !Second value is missing
+              write(*,'(a,1x,i0)') "No second value found in directed list on line", i
+              error stop
+           end if
+           env%directed_list(i, 1) = dum(1:index-1)
+           env%directed_list(i, 2) = dum(index+1:)
+           !Remove multiple spaces
+           env%directed_list(i, 2) = adjustl(env%directed_list(i, 2))
+           !Check, if spaces are still in second argument (e.g. a third number is giveb)
+           index = SCAN(trim(env%directed_list(i, 2)), delim_space)
+           if (index == 0) index = SCAN(trim(dum), delim_tab)
+           if (index /= 0) then
+              write(*,'(a,1x,i0)') "Too many values at line", i
+              error stop
+           end if
+           !> Make array with which solvent molecule at which atom to add
+           read(env%directed_list(i,2), *, iostat=io) env%directed_number(i)
+           env%directed_number(i) = sum(env%directed_number)
+           if (io/= 0) then
+              write(*,'(a,1x,i0)') "Second value is no number in line", i
+              error stop
+           end if
+        end do
+     close(ich)
+     write(*,*) 'Performing directed docking'
+     do i=1, nlines
+        write(*,'(a,1x,a,1x,a,1x,a)') 'Docking', trim(env%directed_list(i,2)),&
+               & 'solvent molecules at', trim(env%directed_list(i,1))
+     end do
+
+end subroutine read_directed_input
+
+
 subroutine qcg_grow(env, solu, solv, clus, tim)
    use iso_fortran_env, wp => real64
    use crest_data

src/qcg/solvtool_misc.f90

@@ -150,9 +150,10 @@ subroutine xtb_dock(env, fnameA, fnameB, solu, clus)
 
    type(systemdata)                :: env
    type(zmolecule), intent(in)     :: solu, clus
-   character(len=*), intent(in)     :: fnameA, fnameB
+   character(len=*), intent(in)    :: fnameA, fnameB
    character(len=80)               :: pipe
    character(len=512)              :: jobcall
+   integer                         :: i, ich
 
    call remove('xtb_dock.out')
    call remove('xcontrol')
@@ -162,6 +163,21 @@ subroutine xtb_dock(env, fnameA, fnameB, solu, clus)
 !---- writing wall pot in xcontrol
    call write_wall(env, solu%nat, solu%ell_abc, clus%ell_abc, 'xcontrol')
 
+!---- Write directed stuff, if requested
+   if (allocated(env%directed_file)) then
+      do i=1, size(env%directed_number)
+         if &
+         & ((i==1 .and. env%directed_number(i) >= clus%nmol) .OR. &
+         & (env%directed_number(i) >= clus%nmol .and. env%directed_number(i-1) < clus%nmol)) &
+         & then
+            open(newunit=ich, file='xcontrol', status='old', position='append', action='write')
+                 write(ich,'("$directed")')
+                 write(ich,'(a,1x,a)') 'atoms:', trim(env%directed_list(i,1))
+                 write(ich,'("$end")')
+         end if
+      end do
+   end if
+
 !--- Setting threads
 !   if(env%autothreads)then
    call ompautoset(env%threads, 7, env%omp, env%MAXRUN, 1) !set the global OMP/MKL variables for the xtb jobs