Skip to content

Commit

Permalink
aISS docking update for next QCG version (#956)
Browse files Browse the repository at this point in the history 
* Update aISS docking for next generation QCG

Signed-off-by: cplett <[email protected]>

* Cleanup

Signed-off-by: cplett <[email protected]>

---------

Signed-off-by: cplett <[email protected]>
  • Loading branch information
@cplett
cplett authored Feb 1, 2024
1 parent 25d87c3 commit a523bf1
Show file tree
Hide file tree
Showing 5 changed files with 103 additions and 11 deletions.
2 changes: 2 additions & 0 deletions src/docking/param.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -135,6 +135,8 @@ module xtb_docking_param
integer, parameter :: p_atom_pot = 2
!> Attractive atom-centered potential
integer, parameter :: p_atom_att = 3
!> Attractive atom-centered potential for QCG mode
integer, parameter :: p_atom_qcg = 4
!> Wall pot for directed docking (Not used)
integer, parameter :: p_wall_pot = 1
integer :: place_wall_pot
Expand Down
2 changes: 2 additions & 0 deletions src/docking/search_nci.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -880,11 +880,13 @@ END SUBROUTINE Quicksort
do j = 1, molA%n
comb%xyz(1:3, j) = molA%xyz(1:3, j) !comb overwritten with A, as it is changed upon geo_opt
end do

counter = 0
do j = molA%n + 1, molA%n + molB%n
counter = counter + 1
comb%xyz(1:3, j) = xyz_tmp(1:3, counter) !combined molA and shifted molB
end do

select type (calc)
type is (TGFFCalculator)
call restart_gff(env, comb, calc)
Expand Down
24 changes: 21 additions & 3 deletions src/docking/set_module.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,6 +28,7 @@ module xtb_docking_set_module
use xtb_type_atomlist
use xtb_mctc_strings, only : parse
use xtb_setmod
use xtb_mctc_convert, only : autokcal

implicit none

Expand Down Expand Up @@ -250,6 +251,9 @@ subroutine set_directed(env,key,val,nat,at,idMap,xyz)
character(len=p_str_length),dimension(p_arg_length) :: argv

logical, save :: set1 = .true.
logical, save :: set2 = .true.
logical, save :: set3 = .true.
logical, save :: set4 = .true.

call atl%resize(nat)

Expand Down Expand Up @@ -303,6 +307,16 @@ subroutine set_directed(env,key,val,nat,at,idMap,xyz)
call atl%to_list(list)
directedset%atoms = list
directedset%n = size(list)
case('expo')
if (getValue(env,val,ddum).and.set2) directedset%expo(1) = ddum
set2 = .false.
case('prefac')
!Prefactor is given in kcal, but xtb energies are in hartree
if (getValue(env,val,ddum).and.set3) directedset%val(1) = ddum / autokcal
set3 = .false.
case('midpoint')
if (getValue(env,val,ddum).and.set4) directedset%expo(2) = ddum
set4 = .false.
end select
call write_set_directed(env%unit)

Expand Down Expand Up @@ -331,6 +345,7 @@ subroutine set_logicals(env, key)
logical, save :: set13 = .true.
logical, save :: set14 = .true.
logical, save :: set15 = .true.
logical, save :: set16 = .true.
select case (key)
case default ! do nothing
call env%warning("the key '"//key//"' is not recognized by scc", source)
Expand Down Expand Up @@ -367,11 +382,9 @@ subroutine set_logicals(env, key)
set8 = .false.
case('qcg')
if (set8) then
!The following setups will become obsolete with next Crest version
maxparent = 50 ! # of parents in gene pool 100
maxgen = 7 ! # of generations 10
! mxcma = 250 ! R points in CMA search 1000
! stepr = 4.0 ! R grid step in Bohr 2.5
! stepa = 60 ! angular grid size in deg. 45
n_opt = 5 ! # of final grad opts 15
qcg = .true.
end if
Expand All @@ -397,6 +410,11 @@ subroutine set_logicals(env, key)
case('repulsive')
if (set15) directed_type = p_atom_pot
set15 = .false.
case('solv_tool')
if (set16) directed_type = p_atom_qcg
allocate(directedset%expo(2), source=0.0_wp)
allocate(directedset%val(1), source=0.0_wp)
set16 = .false.
end select
end subroutine set_logicals

Expand Down
84 changes: 77 additions & 7 deletions src/iff/iff_energy.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,7 +19,7 @@ module xtb_iff_iffenergy
use xtb_mctc_accuracy, only: wp
use xtb_type_environment, only: TEnvironment
use xtb_docking_param
!use xtb_sphereparam
use xtb_sphereparam, only: number_walls
use xtb_sphereparam, only : sphere_alpha, wpot, polynomial_cavity_list, sphere, cavity_egrad,&
& cavitye, maxwalls, rabc, sphere
implicit none
Expand Down Expand Up @@ -78,6 +78,10 @@ subroutine iff_e(env, n, n1, n2, at1, at2, neigh, A1, c02, q01, q02, c6ab,&
real(wp) :: c2(3, n2), cl2(4, 10*n2), dip(3), h, zz1, zz2, nn1, nn2, ees
real(wp) :: rab(n2, n1), rotm(3, 3), q1(n1), q2(n2), AL2(4, n2*10), A2(3, n2)
real(wp) :: r2ab(n2, n1), r0tmp(n2, n1), r0tmp2(n2, n1), r, oner, sab(n2, n1)

!> QCG stuff
real(wp),allocatable :: rab_solu_solv(:,:)

!> Thermo stuff
real(wp) :: aa, bb, cc, avmom, wt, h298, g298, ts298, zp, symn, freq(3*n)
real(wp) :: xyz(3, n)
Expand All @@ -86,7 +90,7 @@ subroutine iff_e(env, n, n1, n2, at1, at2, neigh, A1, c02, q01, q02, c6ab,&
real(wp), parameter ::au = 627.509541d0
integer :: pair(2, n1*n2)
integer :: npair
integer :: i, j, k, kk, i1, j2, iat, jat, n3, metal(86)
integer :: i, j, k, kk, i1, j2, iat, jat, n3, metal(86), counter
integer :: at(n)
logical :: linear, ATM, ijh, ijx
character(len=4) :: pgroup
Expand Down Expand Up @@ -125,13 +129,34 @@ subroutine iff_e(env, n, n1, n2, at1, at2, neigh, A1, c02, q01, q02, c6ab,&
dgrrho = 0
npair = 0
pair = 0
counter = 0

!> Distances of solute and screened solvent positions for QCG attractive potential
if(directedset%n > 0 .and. directed_type == p_atom_qcg) then
allocate(rab_solu_solv(directedset%n, n2), source=0.0_wp)
end if

!> Distance between all atoms of molA and molB
do i = 1, n1
!> Check if qcg mode and i is a solute atom
if(directedset%n > 0 .and. directed_type == p_atom_qcg) then
if(any(i == directedset%atoms)) then
counter = counter + 1
end if
end if
do j = 1, n2
r2 = (A1(1, i) - A2(1, j))**2&
&+ (A1(2, i) - A2(2, j))**2&
&+ (A1(3, i) - A2(3, j))**2
r2ab(j, i) = r2 + 1.d-12
rab(j, i) = sqrt(r2) + 1.d-12
if(directedset%n > 0 .and. directed_type == p_atom_qcg) then
if(any(i == directedset%atoms)) then
!Safe distance in Angström only, if solute atom
! (important if molA is solute with already added solvents)
rab_solu_solv (counter,j) = sqrt(r2)*autoang
end if
end if
if (r2 .gt. 12000.0d0 .or. r2 .lt. 5.0d0) cycle ! induction cut-off
if (metal(at1(i)) .eq. 1 .or. metal(at2(j)) .eq. 1) cycle
npair = npair + 1
Expand All @@ -140,7 +165,6 @@ subroutine iff_e(env, n, n1, n2, at1, at2, neigh, A1, c02, q01, q02, c6ab,&
end do
end do


! atomic density overlap
call ovlp(n1, n2, r2ab, den1, den2, sab)

Expand Down Expand Up @@ -187,6 +211,7 @@ subroutine iff_e(env, n, n1, n2, at1, at2, neigh, A1, c02, q01, q02, c6ab,&
ed = ed - directedset%val(i) * (exp(-0.01*((rab(j,i)-2)**2)))
!- Because ed is substracted at the end
end if

if (rab(j, i) .gt. 25) cycle
nn12 = nn1*(z2(j) - q2(j))
ep = ep + nn12*sab(j, i)/rab(j, i) ! ovlp S(S+1) is slightly worse, S(S-1) very slightly better
Expand Down Expand Up @@ -280,7 +305,7 @@ subroutine iff_e(env, n, n1, n2, at1, at2, neigh, A1, c02, q01, q02, c6ab,&

! Cavity Energies (For Wall Potentials)
esph = 0.0_wp
if(qcg) then
if(qcg .and. (number_walls > 0)) then
do i=1, maxwalls
if(.not. allocated(wpot(i)%list)) then
rabc(1:3) = wpot(i)%radius(1:3)
Expand All @@ -303,6 +328,20 @@ subroutine iff_e(env, n, n1, n2, at1, at2, neigh, A1, c02, q01, q02, c6ab,&

e = -ed + es + ep + esl + ei + gsolv + eabc + ect + esph

!For QCG v2, the closer the solvent to the solute, the more attractive the energy is
if(directedset%n > 0 .and. directed_type == p_atom_qcg) then
!Energie in Hartree
!Potential of type (a/2)*erf(-(b*r-c*b))+(a/2)
!a=E_int solv-solv and is the maximum value of potential
!b=slope of error function
!c=mid point of error function
e = e - &
& ((directedset%val(1)/2) * &
& erf(-(directedset%expo(1) * minval(rab_solu_solv)) + (directedset%expo(2) * directedset%expo(1))) &
& + (directedset%val(1)/2))
! &directedset%val(1) * (exp(-directedset%expo(1) * (minval(rab_solu_solv))**2))
end if

if(isnan(e)) e = 1.0_wp**42

! all done now output
Expand All @@ -327,7 +366,6 @@ subroutine iff_e(env, n, n1, n2, at1, at2, neigh, A1, c02, q01, q02, c6ab,&
write (env%unit, '('' '',F14.8,'' <== Gint total'')') e*au
end if


end subroutine iff_e

subroutine drudescf(n1,n2,nl1,nl2,at1,at2,npair,pair,xyz1,qat1,qdr1,xyz2,&
Expand Down Expand Up @@ -848,14 +886,26 @@ subroutine intermole_probe(n1, nl1, at1, q1, c1, cl1, c2, l1, cn1, c6xx, e, eqp)
real(wp), intent(in) :: cn1(n1)
real(wp), intent(out) :: e, eqp

integer :: i, j, iat
!> QCG stuff
real(wp),allocatable :: rab_solu_solv(:)

integer :: i, j, iat, counter
real(wp) :: ed, ep
real(wp) :: r, r0, r2, rr, av, c6, r6, r8, r42, tmpr
real(wp) :: R06, R08, dc6, t6, t8
real(wp), parameter ::au = 627.509541d0
real(wp), parameter ::autoang = 0.52917726d0

ed = 0
ep = 0
eqp = 0
counter = 0

!> Distances of solute and screened solvent positions for QCG attractive potential
if(directedset%n > 0 .and. directed_type == p_atom_qcg) then
allocate(rab_solu_solv(directedset%n), source=0.0_wp)
end if

do i = 1, n1 !cycle over every atom of molA
iat = at1(i)
r2 = (c1(1, i) - c2(1))**2 + (c1(2, i) - c2(2))**2 + (c1(3, i) - c2(3))**2
Expand All @@ -866,6 +916,15 @@ subroutine intermole_probe(n1, nl1, at1, q1, c1, cl1, c2, l1, cn1, c6xx, e, eqp)
r42 = rrab(iat, probe_atom_type)
R06 = r0ab6(iat, probe_atom_type)
R08 = r0ab8(iat, probe_atom_type)
if(directedset%n > 0 .and. directed_type == p_atom_qcg) then
if(any(i == directedset%atoms)) then
counter=counter + 1
!Safe distance in Angström only, if solute atom
! (important if molA is solute with already added solvents)
rab_solu_solv (counter) = r *autoang
end if
end if

dc6 = 1./(r6 + R06) + r42/(r6*r2 + R08)
ed = ed + dc6*c6xx(i)
if(directedset%n > 0 .and. directed_type == p_atom_att) then
Expand All @@ -889,8 +948,19 @@ subroutine intermole_probe(n1, nl1, at1, q1, c1, cl1, c2, l1, cn1, c6xx, e, eqp)
end do

e = 8.0d0*ep - ed*0.70d0
if(directedset%n > 0 .and. directed_type == p_atom_qcg) then
!Energie in Hartree
!Potential of type (a/2)*erf(-(b*r-c*b))+(a/2)
!a=E_int solv-solv and is the maximum value of potential
!b=slope of error function
!c=mid point of error function
e = e - &
& (directedset%val(1)/2) * &
& erf(-((directedset%expo(1) * minval(rab_solu_solv)) - (directedset%expo(2) * directedset%expo(1)))) &
& + (directedset%val(1)/2)
end if

eqp = eqp*0.1 ! 0.1=probe charge
eqp = eqp*0.1 ! 0.1=probe charge, is added to e in search_nci.f90

end subroutine intermole_probe

Expand Down
2 changes: 1 addition & 1 deletion src/prog/dock.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -282,7 +282,7 @@ subroutine xtbDock(env, argParser)
!> First set optimization level in global parameters
call set_optlvl(env)
call docking_search(env, molA, molB, iff_data%n, iff_data%n1, iff_data%n2,&
iff_data%at1, iff_data%at2, iff_data%neigh, iff_data%xyz1,&
& iff_data%at1, iff_data%at2, iff_data%neigh, iff_data%xyz1,&
& iff_data%xyz2, iff_data%q1, iff_data%q2, iff_data%c6ab,&
& iff_data%z1, iff_data%z2,&
& iff_data%nlmo1, iff_data%nlmo2, iff_data%lmo1,&
Expand Down

0 comments on commit a523bf1

Please sign in to comment.