Introducing a scaling factor for directed docking (#868)

* Introducing a scaling factor for directed docking

Signed-off-by: cplett <82893466+cplett@users.noreply.github.com>

* Typo

---------

Signed-off-by: cplett <82893466+cplett@users.noreply.github.com>
Co-authored-by: Marcel Stahn <70513124+MtoLStoN@users.noreply.github.com>

src/docking/param.f90

@@ -135,7 +135,7 @@ module xtb_docking_param
    integer, parameter :: p_atom_pot = 2
    !> Attractive atom-centered potential
    integer, parameter :: p_atom_att = 3 
-   !Wall pot for directed docking
+   !> Wall pot for directed docking (Not used)
    integer, parameter :: p_wall_pot = 1 
    integer :: place_wall_pot
    !QCG mode (special treatment of wall potentials)

src/docking/set_module.f90

@@ -129,7 +129,7 @@ subroutine rdblock_docking2(env,handler,line,id,nat,at,idMap,xyz,err)
       character(len=:),allocatable,intent(out) :: line
       character(len=:),allocatable :: key
       character(len=:),allocatable :: val
-      integer :: ie
+      integer :: ie, ia
       logical :: exitRun
       do
          call getline(id,line,err)
@@ -139,12 +139,18 @@ subroutine rdblock_docking2(env,handler,line,id,nat,at,idMap,xyz,err)
 
          ! find the first colon
          ie = index(line,colon)
-         if ((line.eq.'')) cycle
-         if (ie .eq. 0) then! cycle
+         ia = index(line,equal)
+         if ((line == '')) cycle
+         if (ie == 0 .AND. ia == 0) then !Only if no : or = is in string
             call set_logicals(env, line)
-         else
-            key = trim(line(:ie-1))
-            val = trim(adjustl(line(ie+1:)))
+         else 
+            if (ia /= 0) then
+               key = trim(line(:ia-1))
+               val = trim(adjustl(line(ia+1:)))
+            else
+               key = trim(line(:ie-1))
+               val = trim(adjustl(line(ie+1:)))
+            end if
             call handler(env,key,val,nat,at,idMap,xyz)
             call env%check(exitRun)
             if (exitRun) then
@@ -242,6 +248,8 @@ subroutine set_directed(env,key,val,nat,at,idMap,xyz)
 
       integer  :: narg
       character(len=p_str_length),dimension(p_arg_length) :: argv
+
+      logical, save :: set1 = .true.
 
       call atl%resize(nat)
 
@@ -254,6 +262,9 @@ subroutine set_directed(env,key,val,nat,at,idMap,xyz)
       endif
       select case(key)
       case default ! ignore, don't even think about raising them
+      case('scaling factor')
+         if (getValue(env,val,ddum).and.set1) directedset%fc = ddum
+         set1 = .false.
       case('elements')
          call atl%new
          do idum = 1, narg

src/prog/dock.f90

@@ -983,6 +983,7 @@ subroutine read_userdata_iff(fname,env,mol)
             select case(line(2:))
             case('directed'      )
                if (set%verbose) write(env%unit,'(">",1x,a)') line(2:)
+               directedset%fc = 1.0_wp !Default scaling. Can be changed by user
                call rdblock_docking2(env,set_directed,line,id,mol%n,mol%at,idMap,mol%xyz,err)
             case default ! unknown keyword -> ignore, we don't raise them
                call getline(id,line,err)
@@ -1031,14 +1032,15 @@ subroutine get_repulsive_pot(env,xyz,comb)
       end do
       !> Changing the distance to a repulsive potential sitting on every atom other then
       !  the defined docking atoms. This potentail is a damped exponential increase.
-      !  It is later in the energy calculation and RG screening added in sitance depdence to
+      !  It is later in the energy calculation and RG screening added in distance depdence to
       !  docked molecule via 1/r²
       do i=1, comb%n
          if(any(i == directedset%atoms)) cycle !Potential zero for atoms in defined docking region
          dist = directedset%val(i)
          rep_pot = 0.1*erf(0.07 * dist - 0.28) !Potential starts at distance of 4
          if(rep_pot < 0.0_wp) rep_pot = 0.0_wp
          directedset%val(i) = rep_pot !Overwrite distance with repulsive Potential
+         directedset%val(i) = directedset%val(i) * directedset%fc !Scaling if requested
       end do
    end subroutine get_repulsive_pot
 
@@ -1057,6 +1059,7 @@ subroutine get_attractive_pot(env,comb)
       do i = 1, comb%n
          if(any(i == directedset%atoms)) then
            directedset%val(i) = attractive_pot !attractive pot is negative
+           directedset%val(i) = directedset%val(i) * directedset%fc !Scaling if requested
          else 
            directedset%val(i) = 0.0_wp
          end if