Skip to content

Commit

Permalink
aISS update (#1089)
Browse files Browse the repository at this point in the history 
* aISS bug fixes

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

* Refactoring aISS docking code

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

* More unit tests for aISS docking algorithm

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

* aISS bug fix for GFN-FF optimizations

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

* Docking code clean up

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

---------

Signed-off-by: cplett <[email protected]>
  • Loading branch information
@cplett
cplett authored Sep 2, 2024
1 parent 043d587 commit 7662ad3
Show file tree
Hide file tree
Showing 5 changed files with 539 additions and 49 deletions.
31 changes: 16 additions & 15 deletions src/docking/param.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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, *)
Expand All @@ -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,&
Expand Down Expand Up @@ -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
Expand All @@ -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)

Expand Down Expand Up @@ -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)

Expand Down
50 changes: 36 additions & 14 deletions src/docking/search_nci.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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,&
Expand Down Expand Up @@ -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)
Expand Down
2 changes: 1 addition & 1 deletion src/iff/iff_prepare.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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.

Expand Down
3 changes: 1 addition & 2 deletions src/prog/dock.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand Down
Loading

0 comments on commit 7662ad3

Please sign in to comment.