From ce88d0309c5a1e146146281bc3cdc54bb9308245 Mon Sep 17 00:00:00 2001 From: cplett <82893466+cplett@users.noreply.github.com> Date: Fri, 26 Jan 2024 14:23:46 +0100 Subject: [PATCH 1/2] Update aISS docking for next generation QCG Signed-off-by: cplett <82893466+cplett@users.noreply.github.com> --- src/docking/param.f90 | 2 + src/docking/search_nci.f90 | 2 + src/docking/set_module.f90 | 24 +++++++++-- src/iff/iff_energy.f90 | 85 ++++++++++++++++++++++++++++++++++---- src/prog/dock.f90 | 2 +- 5 files changed, 104 insertions(+), 11 deletions(-) diff --git a/src/docking/param.f90 b/src/docking/param.f90 index ce1053988..5c872c705 100644 --- a/src/docking/param.f90 +++ b/src/docking/param.f90 @@ -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 diff --git a/src/docking/search_nci.f90 b/src/docking/search_nci.f90 index 82db2814a..736b4e595 100644 --- a/src/docking/search_nci.f90 +++ b/src/docking/search_nci.f90 @@ -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) diff --git a/src/docking/set_module.f90 b/src/docking/set_module.f90 index 879cd0349..c1305218a 100644 --- a/src/docking/set_module.f90 +++ b/src/docking/set_module.f90 @@ -245,11 +245,15 @@ subroutine set_directed(env,key,val,nat,at,idMap,xyz) integer, allocatable :: list(:) real(wp) :: ddum logical :: ldum + real(wp), parameter ::au = 627.509541d0 integer :: narg 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) @@ -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 / au + 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) @@ -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) @@ -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 @@ -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 diff --git a/src/iff/iff_energy.f90 b/src/iff/iff_energy.f90 index 756e9a9d4..75f39194f 100644 --- a/src/iff/iff_energy.f90 +++ b/src/iff/iff_energy.f90 @@ -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 @@ -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) @@ -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 @@ -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 @@ -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) @@ -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 @@ -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) @@ -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 @@ -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,& @@ -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 @@ -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 @@ -889,8 +948,20 @@ 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) +! &directedset%val(1) * (exp(-directedset%expo(1) * (minval(rab_solu_solv))**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 diff --git a/src/prog/dock.f90 b/src/prog/dock.f90 index a26bab7f1..a82c3dbc8 100644 --- a/src/prog/dock.f90 +++ b/src/prog/dock.f90 @@ -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,& From 51eeca214fdbb2f2ee91b85de50e42e1bf473b45 Mon Sep 17 00:00:00 2001 From: cplett <82893466+cplett@users.noreply.github.com> Date: Wed, 31 Jan 2024 10:02:11 +0100 Subject: [PATCH 2/2] Cleanup Signed-off-by: cplett <82893466+cplett@users.noreply.github.com> --- src/docking/set_module.f90 | 4 ++-- src/iff/iff_energy.f90 | 1 - 2 files changed, 2 insertions(+), 3 deletions(-) diff --git a/src/docking/set_module.f90 b/src/docking/set_module.f90 index c1305218a..e412926a3 100644 --- a/src/docking/set_module.f90 +++ b/src/docking/set_module.f90 @@ -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 @@ -245,7 +246,6 @@ subroutine set_directed(env,key,val,nat,at,idMap,xyz) integer, allocatable :: list(:) real(wp) :: ddum logical :: ldum - real(wp), parameter ::au = 627.509541d0 integer :: narg character(len=p_str_length),dimension(p_arg_length) :: argv @@ -312,7 +312,7 @@ subroutine set_directed(env,key,val,nat,at,idMap,xyz) 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 / au + 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 diff --git a/src/iff/iff_energy.f90 b/src/iff/iff_energy.f90 index 75f39194f..dcecd3506 100644 --- a/src/iff/iff_energy.f90 +++ b/src/iff/iff_energy.f90 @@ -958,7 +958,6 @@ subroutine intermole_probe(n1, nl1, at1, q1, c1, cl1, c2, l1, cn1, c6xx, e, eqp) & (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 eqp = eqp*0.1 ! 0.1=probe charge, is added to e in search_nci.f90