diff --git a/src/docking/param.f90 b/src/docking/param.f90 index 5c872c705..8aa964f29 100644 --- a/src/docking/param.f90 +++ b/src/docking/param.f90 @@ -359,7 +359,7 @@ end subroutine diptot ! axis ! molw is the weigth, sum3 the CMA (all in a.u.) - subroutine axis(numat, nat, coord, sum3, sumw, eig, evec) + subroutine axis_docking(numat, nat, coord, sum3, sumw, eig, evec) integer, intent(in) ::numat, nat(numat) real(wp), intent(in) :: coord(3, *) @@ -368,7 +368,7 @@ subroutine axis(numat, nat, coord, sum3, sumw, eig, evec) real(wp) :: t(6) real(wp) :: x(numat), y(numat), z(numat) real(wp) :: sumwx, sumwy, sumwz - real(wp) :: atmass + real(wp) :: atmass_dock integer :: i real(wp) :: ams(107) data ams/1.00790d0, 4.00260d0, 6.94000d0, 9.01218d0,& @@ -396,11 +396,11 @@ subroutine axis(numat, nat, coord, sum3, sumw, eig, evec) sumwz = 0.d0 do i = 1, numat - atmass = ams(nat(i)) - sumw = sumw + atmass - sumwx = sumwx + atmass*coord(1, i) - sumwy = sumwy + atmass*coord(2, i) - sumwz = sumwz + atmass*coord(3, i) + atmass_dock = ams(nat(i)) + sumw = sumw + atmass_dock + sumwx = sumwx + atmass_dock*coord(1, i) + sumwy = sumwy + atmass_dock*coord(2, i) + sumwz = sumwz + atmass_dock*coord(3, i) end do sum3(1) = sumwx/sumw @@ -418,19 +418,19 @@ subroutine axis(numat, nat, coord, sum3, sumw, eig, evec) end do do i = 1, numat - atmass = ams(nat(i)) - t(1) = t(1) + atmass*(y(i)**2 + z(i)**2) - t(2) = t(2) - atmass*x(i)*y(i) - t(3) = t(3) + atmass*(z(i)**2 + x(i)**2) - t(4) = t(4) - atmass*z(i)*x(i) - t(5) = t(5) - atmass*y(i)*z(i) - t(6) = t(6) + atmass*(x(i)**2 + y(i)**2) + atmass_dock = ams(nat(i)) + t(1) = t(1) + atmass_dock*(y(i)**2 + z(i)**2) + t(2) = t(2) - atmass_dock*x(i)*y(i) + t(3) = t(3) + atmass_dock*(z(i)**2 + x(i)**2) + t(4) = t(4) - atmass_dock*z(i)*x(i) + t(5) = t(5) - atmass_dock*y(i)*z(i) + t(6) = t(6) + atmass_dock*(x(i)**2 + y(i)**2) end do call rsp(t, 3, 3, eig, evec) eig = eig/sumw - end subroutine axis + end subroutine axis_docking subroutine cmadock(n, numat, nat, coord, sum3) @@ -1225,6 +1225,7 @@ subroutine combine_mol(comb, molA, molB) call init(comb, at, xyz) comb%chrg = molA%chrg + molB%chrg comb%uhf = molA%uhf + molB%uhf + atmass=comb%atmass !Setting global atmass array required in axis module deallocate (at) deallocate (xyz) diff --git a/src/docking/search_nci.f90 b/src/docking/search_nci.f90 index e1d846e10..92b454f5b 100644 --- a/src/docking/search_nci.f90 +++ b/src/docking/search_nci.f90 @@ -48,13 +48,14 @@ module xtb_docking_search_nci & gff_print, gfnff_param_dealloc use xtb_constrain_param, only: read_userdata use xtb_fixparam - use xtb_disp_ncoord, only: ncoord_gfn, ncoord_erf + use xtb_disp_ncoord, only: ncoord_gfn, ncoord_erf, ncoord_d3 use xtb_scc_core, only: iniqshell use xtb_eeq, only: goedecker_chrgeq use xtb_basis, only: newBasisset use xtb_gfnff_neighbor, only: TNeigh use xtb_io_writer, only : writeMolecule use xtb_mctc_filetypes, only : generateFileName + use xtb_iniq, only: iniqcn implicit none private @@ -203,9 +204,9 @@ END SUBROUTINE Quicksort if (.not. fulle) write (env%unit, *) ! just output - call axis(n1, at1, xyz1, dum, pmass, dum2, dum3) + call axis_docking(n1, at1, xyz1, dum, pmass, dum2, dum3) r1 = sqrt(dum2(1) + dum2(2) + dum2(3)) - call axis(n2, at2, xyz2, dum, pmass, dum2, dum3) + call axis_docking(n2, at2, xyz2, dum, pmass, dum2, dum3) r2 = sqrt(dum2(1) + dum2(2) + dum2(3)) call rcma(n1, xyz1, at1, n2, xyz2, at2, r, rmin) write (env%unit, '('' Method for final opts. :'',1x,a )') optlvl @@ -500,7 +501,7 @@ END SUBROUTINE Quicksort ! xyztmp contains now each gridpoint belonging to fragment i ! same for attmp (anyway always probe_atom_type) ! determine size (=R) of gridpoint cluster belonging to fragment i - call axis(k, attmp, xyztmp, dum, pmass, dum2, dum3) + call axis_docking(k, attmp, xyztmp, dum, pmass, dum2, dum3) r = sqrt(dum2(1) + dum2(2) + dum2(3)) !> check if moleculeB is small enough to fit in gridpoint cluster of fragmen i @@ -1065,11 +1066,17 @@ subroutine restart_gff(env, mol, calc) call gfnff_param_dealloc(calc%topo) call newD3Model(calc%topo%dispm, mol%n, mol%at) call gfnff_set_param(mol%n, calc%gen, calc%param) - if (.not. allocated(calc%neigh%nb)) allocate (calc%neigh%nb(calc%neigh%numnb, mol%n, calc%neigh%numctr), source=0) - if (.not.allocated(calc%topo%qfrag)) & - & allocate( calc%topo%qfrag(mol%n), source = 0.0d0 ) - if (.not.allocated(calc%topo%fraglist)) & - & allocate( calc%topo%fraglist(mol%n), source = 0 ) + if (allocated(calc%neigh%nb)) deallocate(calc%neigh%nb) + allocate (calc%neigh%nb(calc%neigh%numnb, mol%n, calc%neigh%numctr), source=0) + if (allocated(calc%topo%qfrag)) deallocate(calc%topo%qfrag) + allocate( calc%topo%qfrag(mol%n), source = 0.0d0 ) + if (allocated(calc%topo%fraglist)) deallocate(calc%topo%fraglist) + allocate( calc%topo%fraglist(mol%n), source = 0 ) + if (allocated(calc%neigh%iTrUsed)) deallocate(calc%neigh%iTrUsed) + if (allocated(calc%neigh%bpair)) deallocate(calc%neigh%bpair) + if (allocated(calc%neigh%blist)) deallocate(calc%neigh%blist) + if (allocated(calc%neigh%vbond)) deallocate(calc%neigh%vbond) + if (allocated(calc%neigh%nr_hb)) deallocate(calc%neigh%nr_hb) calc%topo%qfrag(1) = set%ichrg calc%topo%qfrag(2:mol%n) = 0.0_wp call gfnff_ini(env, .false., ini, mol, calc%gen,& @@ -1116,15 +1123,30 @@ subroutine restart_xTB(env, mol, chk, calc, basisset) calc%maxiter = set%maxscciter call chk%wfn%allocate(mol%n, calc%basis%nshell, calc%basis%nao) + ! Make sure number of electrons is initialized an multiplicity is consistent + chk%wfn%nel = nint(sum(mol%z) - mol%chrg) + chk%wfn%nopen = mol%uhf + if (chk%wfn%nopen == 0 .and. mod(chk%wfn%nel, 2) /= 0) chk%wfn%nopen = 1 + !> EN charges and CN - call ncoord_gfn(mol%n, mol%at, mol%xyz, cn) + if (set%gfn_method < 2) then + call ncoord_d3(mol%n, mol%at, mol%xyz, cn) + else + call ncoord_gfn(mol%n, mol%at, mol%xyz, cn) + end if if (mol%npbc > 0) then chk%wfn%q = real(set%ichrg, wp)/real(mol%n, wp) else - call ncoord_erf(mol%n, mol%at, mol%xyz, cn) - call goedecker_chrgeq(mol%n,mol%at,mol%xyz,real(set%ichrg,wp),cn,dcn,chk%wfn%q,dq,er,g,& - & .false., .false., .false.) - chk%wfn%q = real(set%ichrg, wp)/real(mol%n, wp) + if (set%guess_charges == p_guess_gasteiger) then + call iniqcn(mol%n, mol%at, mol%z, mol%xyz, set%ichrg, 1.0_wp, chk%wfn%q, cn, set%gfn_method, .true.) + else if (set%guess_charges == p_guess_goedecker) then + call ncoord_erf(mol%n, mol%at, mol%xyz, cn) + call goedecker_chrgeq(mol%n, mol%at, mol%xyz, real(set%ichrg, wp), cn, dcn, chk%wfn%q, dq, er, g, & + .false., .false., .false.) + else + call ncoord_gfn(mol%n, mol%at, mol%xyz, cn) + chk%wfn%q = real(set%ichrg, wp) / real(mol%n, wp) + end if end if !> initialize shell charges from gasteiger charges call iniqshell(calc%xtbData,mol%n,mol%at,mol%z,calc%basis%nshell,chk%wfn%q,chk%wfn%qsh,set%gfn_method) diff --git a/src/iff/iff_prepare.f90 b/src/iff/iff_prepare.f90 index 15632a2f4..813be1b32 100644 --- a/src/iff/iff_prepare.f90 +++ b/src/iff/iff_prepare.f90 @@ -193,7 +193,7 @@ subroutine precomp(env, iff_data, mol, etot, mol_num) !> the SP call singlepoint & & (env, mol, chk, calc, egap, set%etemp, set%maxscciter, 2,& - & exist, lgrad, acc, etot, g, sigma, res) + & .false., lgrad, acc, etot, g, sigma, res) set%pr_lmo = .false. diff --git a/src/prog/dock.f90 b/src/prog/dock.f90 index a82c3dbc8..f68976267 100644 --- a/src/prog/dock.f90 +++ b/src/prog/dock.f90 @@ -45,7 +45,7 @@ module xtb_prog_dock use xtb_iff_iffprepare, only: precomp use xtb_iff_iffenergy, only : iff_e use xtb_docking_search_nci, only: docking_search - use xtb_sphereparam, only: sphere, rabc, boxr, init_walls, wpot, maxwalls + use xtb_sphereparam, only: init_walls, maxwalls use xtb_constrain_param, only: read_userdata use xtb_fixparam, only: init_fix use xtb_scanparam, only: init_constr, init_scan, maxconstr, maxscan @@ -232,7 +232,6 @@ subroutine xtbDock(env, argParser) & iff_data%qcm2, iff_data%n, iff_data%at, iff_data%xyz, iff_data%q, icoord, icoord0,& & .false.) - !> CONSTRAINTS & SCANS call init_fix(iff_data%n) call init_split(iff_data%n) diff --git a/test/unit/test_docking.f90 b/test/unit/test_docking.f90 index 3802d0116..a92a5b3b5 100644 --- a/test/unit/test_docking.f90 +++ b/test/unit/test_docking.f90 @@ -29,14 +29,16 @@ subroutine collect_docking(testsuite) type(unittest_type), allocatable, intent(out) :: testsuite(:) testsuite = [ & - new_unittest("eth_wat", test_dock_eth_wat)]!, & -! new_unittest("ellips", test_iff_ellips), & -! ] + new_unittest("dock_gfn2_eth_wat", test_dock_eth_wat_gfn2), & + new_unittest("dock_gfn2_wat_wat_wall", test_dock_wat_wat_gfn2_wall), & + new_unittest("dock_gfn2_wat_wat_attpot", test_dock_wat_wat_gfn2_attpot), & + new_unittest("dock_gfnff_wat_wat", test_dock_wat_wat_gfnff) & + ] end subroutine collect_docking -subroutine test_dock_eth_wat(error) +subroutine test_dock_eth_wat_gfn2(error) use xtb_type_environment, only: TEnvironment, init use xtb_mctc_accuracy, only: wp use xtb_type_environment, only: TEnvironment @@ -112,8 +114,9 @@ subroutine test_dock_eth_wat(error) set%pr_local = .false. call precomp(env, iff_data, molA, molA_e, 1) - call check_(error, molA_e,-11.3943358674_wp, thr=thr) + call check_(error, molA_e,-11.39433586739173_wp, thr=thr) call precomp(env, iff_data, molB, molB_e, 2) + call check_(error, molB_e,-5.070201841808753_wp, thr=thr) call env%checkpoint("LMO computation failed") @@ -150,27 +153,492 @@ subroutine test_dock_eth_wat(error) & iff_data%den1, iff_data%den2, iff_data%gab1, iff_data%gab2, & & set%verbose, 0, e, icoord) - call check_(error, e,0.07745330953243718_wp, thr=thr) + call check_(error, e,.7745330953199690E-01_wp, thr=thr) call env%checkpoint("xtb-IFF sp computation failed") + optlvl="gfn2" call set_optlvl(env) !Sets the optimization level for search in global parameters -! call docking_search(env, molA, molB, n, n1, n2, at1, at2, neigh, xyz1,& -! & xyz2, q1, q2, c6ab, z1, z2,& -! &nlmo1, nlmo2, lmo1, lmo2, rlmo1, rlmo2,& -! &qdr1, qdr2, cn1, cn2, alp1, alp2, alpab, qct1, qct2,& -! &den1, den2, gab1, gab2, molA_e, molB_e,& -! &cprob, e, icoord, comb) -! call env%checkpoint("Docking algorithm failed") + !Settings for fast + maxparent = 30 + maxgen = 3 + mxcma = 250 + stepr = 4.0 + stepa = 60 + n_opt = 2 -! call check_(error, calc%topo%nbond,6) -! -! call check_(error, res_gff%e_total,-0.76480130317838_wp, thr=thr) + 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%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,& + & iff_data%lmo2, iff_data%rlmo1, iff_data%rlmo2,& + & iff_data%qdr1, iff_data%qdr2, iff_data%cn1, iff_data%cn2, iff_data%alp1,& + & iff_data%alp2, iff_data%alpab, iff_data%qct1, iff_data%qct2,& + & iff_data%den1, iff_data%den2, iff_data%gab1, iff_data%gab2, molA_e, molB_e,& + & iff_data%cprob, e, icoord, comb) + + call env%checkpoint("Docking algorithm failed") + + call molA%deallocate + call molB%deallocate + +end subroutine test_dock_eth_wat_gfn2 + + +subroutine test_dock_wat_wat_gfn2_wall(error) + use xtb_type_environment, only: TEnvironment, init + use xtb_mctc_accuracy, only: wp + use xtb_type_environment, only: TEnvironment + use xtb_type_molecule + use xtb_setparam, only: initrand + use xtb_mctc_systools + use xtb_setmod + use xtb_docking_set_module + use xtb_docking_param + use xtb_iff_iffini, only: init_iff + use xtb_iff_iffprepare, only: precomp + use xtb_iff_iffenergy + use xtb_docking_search_nci, only: docking_search + use xtb_mctc_systools + use xtb_iff_data, only: TIFFData + use xtb_sphereparam, only: init_walls, maxwalls, set_sphere_radius,& + & spherepot_type, p_type_polynomial, clear_walls + + type(error_type), allocatable, intent(out) :: error + !> All important variables stored here + type(TIFFData) :: iff_data + real(wp),parameter :: thr = 1.0e-6_wp + real(wp) :: icoord0(6),icoord(6) + integer :: n1 + integer :: at1(3) + real(wp) :: xyz1(3,3) + integer :: n2 + integer :: at2(3) + real(wp) :: xyz2(3,3) + integer,parameter :: n = 6 + integer :: at(n) + real(wp) :: xyz(3,n) + real(wp) :: molA_e, molB_e + + logical, parameter :: restart = .false. + + type(TMolecule) :: molA,molB,comb + real(wp) :: r(3),e + integer :: i, j, k + type(TEnvironment) :: env + character(len=:), allocatable :: fnam + + logical :: exist + + n2 = 3 + at2 = [8,1,1] + xyz2 = reshape(& + &[-14.55824225787638_wp, 0.85763330814882_wp, 0.00000000000000_wp, & + & -12.72520790897730_wp, 0.85763330814882_wp, 0.00000000000000_wp, & + & -15.16924740842229_wp,-0.86379381534203_wp,-0.15285994688912_wp],& + & shape(xyz2)) + + n1=n2 + at1=at2 + xyz1=xyz2 + + qcg=.true. + call set_gbsa(env, 'solvent', "h2o") + call init(env) + + !Molecular stuff: + call init(molA, at1, xyz1) + call init(molB, at2, xyz2) + call iff_data%allocateIFFData(molA%n, molB%n) + + call initrand + + call set_iff_param + fnam = 'xtblmoinfo' + set%pr_local = .false. + + call precomp(env, iff_data, molA, molA_e, 1) +! call check_(error, molA_e,-5.084810260314862_wp, thr=thr) + call precomp(env, iff_data, molB, molB_e, 2) + call check_(error, molB_e,-5.08481026031486_wp, thr=thr) + + call env%checkpoint("LMO computation failed") + + call cmadock(molA%n, molA%n, molA%at, molA%xyz, r) + do i = 1, 3 + molA%xyz(i, 1:molA%n) = molA%xyz(i, 1:molA%n) - r(i) + end do + call cmadock(molB%n, molB%n, molB%at, molB%xyz, r) + do i = 1, 3 + molB%xyz(i, 1:molB%n) = molB%xyz(i, 1:molB%n) - r(i) + end do + + call init_iff(env, iff_data%n1, iff_data%n2, 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%cprob, iff_data%nlmo1, iff_data%nlmo2, iff_data%lmo1, iff_data%lmo2,& + & iff_data%qdr1, iff_data%qdr2, iff_data%rlmo1, iff_data%rlmo2,& + & iff_data%cn1, iff_data%cn2, iff_data%alp1, iff_data%alp2, iff_data%alpab,& + & iff_data%den1, iff_data%den2, iff_data%gab1, iff_data%gab2, iff_data%qcm1,& + & iff_data%qcm2, iff_data%n, iff_data%at, iff_data%xyz, iff_data%q, icoord, icoord0,& + & .false.) + + call check_(error, icoord(1) ,0.0000_wp, thr=1.0e-3_wp) + + !Wall potential stuff: + maxwalls=3 + call init_walls + spherepot_type = p_type_polynomial + call set_sphere_radius([6.7,4.7,3.8],[0.0,0.0,0.0],3,[1,2,3]) !"solute" + call set_sphere_radius([6.9,5.7,5.2],[0.0,0.0,0.0]) !"all" + + call env%checkpoint("Initializing xtb-IFF failed.") + + call iff_e(env, iff_data%n, iff_data%n1, iff_data%n2, 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, iff_data%lmo2, & + & iff_data%rlmo1, iff_data%rlmo2,& + & iff_data%qdr1, iff_data%qdr2, iff_data%cn1, iff_data%cn2, iff_data%alp1, & + & iff_data%alp2, iff_data%alpab, iff_data%qct1, iff_data%qct2, & + & iff_data%den1, iff_data%den2, iff_data%gab1, iff_data%gab2, & + & set%verbose, 0, e, icoord) + + call env%checkpoint("xtb-IFF sp computation failed") + + optlvl="gfn2" + call set_optlvl(env) !Sets the optimization level for search in global parameters + + !Settings for fast + maxparent = 30 + maxgen = 3 + mxcma = 250 + stepr = 4.0 + stepa = 60 + n_opt = 2 + + 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%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,& + & iff_data%lmo2, iff_data%rlmo1, iff_data%rlmo2,& + & iff_data%qdr1, iff_data%qdr2, iff_data%cn1, iff_data%cn2, iff_data%alp1,& + & iff_data%alp2, iff_data%alpab, iff_data%qct1, iff_data%qct2,& + & iff_data%den1, iff_data%den2, iff_data%gab1, iff_data%gab2, molA_e, molB_e,& + & iff_data%cprob, e, icoord, comb) + call env%checkpoint("Docking algorithm failed") + + call molA%deallocate + call molB%deallocate + call clear_walls + maxwalls=0 + +end subroutine test_dock_wat_wat_gfn2_wall + +subroutine test_dock_wat_wat_gfn2_attpot(error) + use xtb_type_environment, only: TEnvironment, init + use xtb_mctc_accuracy, only: wp + use xtb_type_environment, only: TEnvironment + use xtb_type_molecule + use xtb_setparam, only: initrand + use xtb_mctc_systools + use xtb_setmod + use xtb_docking_set_module + use xtb_docking_param + use xtb_iff_iffini, only: init_iff + use xtb_iff_iffprepare, only: precomp + use xtb_iff_iffenergy + use xtb_docking_search_nci, only: docking_search + use xtb_mctc_systools + use xtb_iff_data, only: TIFFData + use xtb_sphereparam, only: init_walls, maxwalls, set_sphere_radius,& + & spherepot_type, p_type_polynomial, clear_walls + + type(error_type), allocatable, intent(out) :: error + !> All important variables stored here + type(TIFFData) :: iff_data + real(wp),parameter :: thr = 1.0e-6_wp + real(wp) :: icoord0(6),icoord(6) + integer :: n1 + integer :: at1(3) + real(wp) :: xyz1(3,3) + integer :: n2 + integer :: at2(3) + real(wp) :: xyz2(3,3) + integer,parameter :: n = 6 + integer :: at(n) + real(wp) :: xyz(3,n) + real(wp) :: molA_e, molB_e + + logical, parameter :: restart = .false. + + type(TMolecule) :: molA,molB,comb + real(wp) :: r(3),e + integer :: i, j, k + type(TEnvironment) :: env + character(len=:), allocatable :: fnam + + logical :: exist + + n2 = 3 + at2 = [8,1,1] + xyz2 = reshape(& + &[-14.55824225787638_wp, 0.85763330814882_wp, 0.00000000000000_wp, & + & -12.72520790897730_wp, 0.85763330814882_wp, 0.00000000000000_wp, & + & -15.16924740842229_wp,-0.86379381534203_wp,-0.15285994688912_wp],& + & shape(xyz2)) + + n1=n2 + at1=at2 + xyz1=xyz2 + + qcg=.true. + call set_gbsa(env, 'solvent', "h2o") + call init(env) + + !Molecular stuff: + call init(molA, at1, xyz1) + call init(molB, at2, xyz2) + call iff_data%allocateIFFData(molA%n, molB%n) + + call initrand + + call set_iff_param + fnam = 'xtblmoinfo' + set%pr_local = .false. + + call precomp(env, iff_data, molA, molA_e, 1) +! call check_(error, molA_e,-5.084810260314862_wp, thr=thr) + call precomp(env, iff_data, molB, molB_e, 2) + call check_(error, molB_e,-5.08481026031486_wp, thr=thr) + + call env%checkpoint("LMO computation failed") + + call cmadock(molA%n, molA%n, molA%at, molA%xyz, r) + do i = 1, 3 + molA%xyz(i, 1:molA%n) = molA%xyz(i, 1:molA%n) - r(i) + end do + call cmadock(molB%n, molB%n, molB%at, molB%xyz, r) + do i = 1, 3 + molB%xyz(i, 1:molB%n) = molB%xyz(i, 1:molB%n) - r(i) + end do + + call init_iff(env, iff_data%n1, iff_data%n2, 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%cprob, iff_data%nlmo1, iff_data%nlmo2, iff_data%lmo1, iff_data%lmo2,& + & iff_data%qdr1, iff_data%qdr2, iff_data%rlmo1, iff_data%rlmo2,& + & iff_data%cn1, iff_data%cn2, iff_data%alp1, iff_data%alp2, iff_data%alpab,& + & iff_data%den1, iff_data%den2, iff_data%gab1, iff_data%gab2, iff_data%qcm1,& + & iff_data%qcm2, iff_data%n, iff_data%at, iff_data%xyz, iff_data%q, icoord, icoord0,& + & .false.) + + call check_(error, icoord(1) ,0.0000_wp, thr=1.0e-3_wp) + + directed_type = p_atom_qcg + allocate(directedset%expo(2), source=0.0_wp) + allocate(directedset%val(1), source=0.0_wp) + directedset%val(1)=0.011928 + directedset%expo(1)=0.01 + directedset%expo(2)=2.000 + directedset%n=3 + directedset%atoms=[1,2,3] + + call env%checkpoint("Initializing xtb-IFF failed.") + + call iff_e(env, iff_data%n, iff_data%n1, iff_data%n2, 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, iff_data%lmo2, & + & iff_data%rlmo1, iff_data%rlmo2,& + & iff_data%qdr1, iff_data%qdr2, iff_data%cn1, iff_data%cn2, iff_data%alp1, & + & iff_data%alp2, iff_data%alpab, iff_data%qct1, iff_data%qct2, & + & iff_data%den1, iff_data%den2, iff_data%gab1, iff_data%gab2, & + & set%verbose, 0, e, icoord) + + call env%checkpoint("xtb-IFF sp computation failed") + + optlvl="gfn2" + call set_optlvl(env) !Sets the optimization level for search in global parameters + + !Settings for fast + maxparent = 30 + maxgen = 3 + mxcma = 250 + stepr = 4.0 + stepa = 60 + n_opt = 2 + + 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%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,& + & iff_data%lmo2, iff_data%rlmo1, iff_data%rlmo2,& + & iff_data%qdr1, iff_data%qdr2, iff_data%cn1, iff_data%cn2, iff_data%alp1,& + & iff_data%alp2, iff_data%alpab, iff_data%qct1, iff_data%qct2,& + & iff_data%den1, iff_data%den2, iff_data%gab1, iff_data%gab2, molA_e, molB_e,& + & iff_data%cprob, e, icoord, comb) + call env%checkpoint("Docking algorithm failed") + + call molA%deallocate + call molB%deallocate + call directedset%deallocate + +end subroutine test_dock_wat_wat_gfn2_attpot + + +subroutine test_dock_wat_wat_gfnff(error) + use xtb_type_environment, only: TEnvironment, init + use xtb_mctc_accuracy, only: wp + use xtb_type_environment, only: TEnvironment + use xtb_type_molecule + use xtb_setparam, only: initrand + use xtb_mctc_systools + use xtb_setmod + use xtb_docking_set_module + use xtb_docking_param + use xtb_iff_iffini, only: init_iff + use xtb_iff_iffprepare, only: precomp + use xtb_iff_iffenergy + use xtb_docking_search_nci, only: docking_search + use xtb_mctc_systools + use xtb_iff_data, only: TIFFData + use xtb_mctc_global, only : persistentEnv + + type(error_type), allocatable, intent(out) :: error + !> All important variables stored here + type(TIFFData) :: iff_data + real(wp),parameter :: thr = 1.0e-6_wp + real(wp) :: icoord0(6),icoord(6) + integer :: n1 + integer :: at1(3) + real(wp) :: xyz1(3,3) + integer :: n2 + integer :: at2(3) + real(wp) :: xyz2(3,3) + integer,parameter :: n = 6 + integer :: at(n) + real(wp) :: xyz(3,n) + real(wp) :: molA_e, molB_e + + logical, parameter :: restart = .false. + + type(TMolecule) :: molA,molB,comb + real(wp) :: r(3),e + integer :: i, j, k + type(TEnvironment) :: env + character(len=:), allocatable :: fnam + + logical :: exist + + n2 = 3 + at2 = [8,1,1] + xyz2 = reshape(& + &[-14.55824225787638_wp, 0.85763330814882_wp, 0.00000000000000_wp, & + & -12.72520790897730_wp, 0.85763330814882_wp, 0.00000000000000_wp, & + & -15.16924740842229_wp,-0.86379381534203_wp,-0.15285994688912_wp],& + & shape(xyz2)) + + n1=n2 + at1=at2 + xyz1=xyz2 + + call set_gbsa(env, 'solvent', "h2o") + !GFN-FF + call set_exttyp('ff') + +! deallocate(persistentEnv%whoami) +! deallocate(persistentEnv%home) +! deallocate(persistentEnv%path) +! deallocate(persistentEnv%xtbpath) +! deallocate(persistentEnv%xtbhome) +! deallocate(persistentEnv%io%log) +! persistentEnv%io%count=0 + + + call init(env) + call init(molA, at1, xyz1) + call init(molB, at2, xyz2) + call iff_data%allocateIFFData(molA%n, molB%n) + + call initrand + + call set_iff_param + fnam = 'xtblmoinfo' + set%pr_local = .false. + + call precomp(env, iff_data, molA, molA_e, 1) +! call check_(error, molA_e,-5.084810260314862_wp, thr=thr) + call precomp(env, iff_data, molB, molB_e, 2) + call check_(error, molB_e,-5.08481026031486_wp, thr=thr) + + call env%checkpoint("LMO computation failed") + + call cmadock(molA%n, molA%n, molA%at, molA%xyz, r) + do i = 1, 3 + molA%xyz(i, 1:molA%n) = molA%xyz(i, 1:molA%n) - r(i) + end do + call cmadock(molB%n, molB%n, molB%at, molB%xyz, r) + do i = 1, 3 + molB%xyz(i, 1:molB%n) = molB%xyz(i, 1:molB%n) - r(i) + end do + + call init_iff(env, iff_data%n1, iff_data%n2, 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%cprob, iff_data%nlmo1, iff_data%nlmo2, iff_data%lmo1, iff_data%lmo2,& + & iff_data%qdr1, iff_data%qdr2, iff_data%rlmo1, iff_data%rlmo2,& + & iff_data%cn1, iff_data%cn2, iff_data%alp1, iff_data%alp2, iff_data%alpab,& + & iff_data%den1, iff_data%den2, iff_data%gab1, iff_data%gab2, iff_data%qcm1,& + & iff_data%qcm2, iff_data%n, iff_data%at, iff_data%xyz, iff_data%q, icoord, icoord0,& + & .false.) + + call check_(error, icoord(1) ,0.0000_wp, thr=1.0e-3_wp) + + call env%checkpoint("Initializing xtb-IFF failed.") + + call iff_e(env, iff_data%n, iff_data%n1, iff_data%n2, 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, iff_data%lmo2, & + & iff_data%rlmo1, iff_data%rlmo2,& + & iff_data%qdr1, iff_data%qdr2, iff_data%cn1, iff_data%cn2, iff_data%alp1, & + & iff_data%alp2, iff_data%alpab, iff_data%qct1, iff_data%qct2, & + & iff_data%den1, iff_data%den2, iff_data%gab1, iff_data%gab2, & + & set%verbose, 0, e, icoord) + + call env%checkpoint("xtb-IFF sp computation failed") + + optlvl = "gfnff" + call set_optlvl(env) !Sets the optimization level for search in global parameters + + !Settings for fast + maxparent = 30 + maxgen = 3 + mxcma = 250 + stepr = 4.0 + stepa = 60 + n_opt = 2 + + 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%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,& + & iff_data%lmo2, iff_data%rlmo1, iff_data%rlmo2,& + & iff_data%qdr1, iff_data%qdr2, iff_data%cn1, iff_data%cn2, iff_data%alp1,& + & iff_data%alp2, iff_data%alpab, iff_data%qct1, iff_data%qct2,& + & iff_data%den1, iff_data%den2, iff_data%gab1, iff_data%gab2, molA_e, molB_e,& + & iff_data%cprob, e, icoord, comb) + call env%checkpoint("Docking algorithm failed") call molA%deallocate call molB%deallocate -end subroutine test_dock_eth_wat +end subroutine test_dock_wat_wat_gfnff end module test_docking