Skip to content

Commit

Permalink
Implementation of a generic docking algorithm (#677)
Browse files Browse the repository at this point in the history
* Implementation of docking algorithm

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

* Update Docking Algorithm

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

* Update docking algorithm for QCG use

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

Signed-off-by: cplett <[email protected]>
cplett authored Oct 10, 2022
1 parent f079c74 commit ee5b291
Showing 34 changed files with 7,297 additions and 103 deletions.
2 changes: 2 additions & 0 deletions src/CMakeLists.txt
Original file line number Diff line number Diff line change
@@ -18,10 +18,12 @@
add_subdirectory("api")
add_subdirectory("coulomb")
add_subdirectory("disp")
add_subdirectory("docking")
add_subdirectory("extern")
add_subdirectory("freq")
add_subdirectory("gfnff")
add_subdirectory("io")
add_subdirectory("iff")
add_subdirectory("main")
add_subdirectory("mctc")
add_subdirectory("param")
26 changes: 26 additions & 0 deletions src/docking/CMakeLists.txt
Original file line number Diff line number Diff line change
@@ -0,0 +1,26 @@
# This file is part of xtb.
#
# Copyright (C) 2022 Sebastian Ehlert, Christoph Plett
#
# xtb is free software: you can redistribute it and/or modify it under
# the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# xtb is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with xtb. If not, see <https://www.gnu.org/licenses/>.

set(dir "${CMAKE_CURRENT_SOURCE_DIR}")

list(APPEND srcs
"${dir}/param.f90"
"${dir}/search_nci.f90"
"${dir}/set_module.f90"
)

set(srcs ${srcs} PARENT_SCOPE)
22 changes: 22 additions & 0 deletions src/docking/meson.build
Original file line number Diff line number Diff line change
@@ -0,0 +1,22 @@
# This file is part of xtb.
#
# Copyright (C) 2022 Sebastian Ehlert, Christoph Plett
#
# xtb is free software: you can redistribute it and/or modify it under
# the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# xtb is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with xtb. If not, see <https://www.gnu.org/licenses/>.

srcs += files(
'search_nci.f90',
'param.f90',
'set_module.f90',
)
1,231 changes: 1,231 additions & 0 deletions src/docking/param.f90

Large diffs are not rendered by default.

1,168 changes: 1,168 additions & 0 deletions src/docking/search_nci.f90

Large diffs are not rendered by default.

417 changes: 417 additions & 0 deletions src/docking/set_module.f90

Large diffs are not rendered by default.

29 changes: 29 additions & 0 deletions src/iff/CMakeLists.txt
Original file line number Diff line number Diff line change
@@ -0,0 +1,29 @@
# This file is part of xtb.
#
# Copyright (C) 2022 Sebastian Ehlert, Christoph Plett
#
# xtb is free software: you can redistribute it and/or modify it under
# the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# xtb is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with xtb. If not, see <https://www.gnu.org/licenses/>.

set(dir "${CMAKE_CURRENT_SOURCE_DIR}")

list(APPEND srcs
"${dir}/calculator.f90"
"${dir}/data.f90"
"${dir}/iff_energy.f90"
"${dir}/iff_ini.f90"
"${dir}/iff_lmo.f90"
"${dir}/iff_prepare.f90"
)

set(srcs ${srcs} PARENT_SCOPE)
174 changes: 174 additions & 0 deletions src/iff/calculator.f90
Original file line number Diff line number Diff line change
@@ -0,0 +1,174 @@
! This file is part of xtb.
!
! Copyright (C) 2022 Stefan Grimme
!
! xtb is free software: you can redistribute it and/or modify it under
! the terms of the GNU Lesser General Public License as published by
! the Free Software Foundation, either version 3 of the License, or
! (at your option) any later version.
!
! xtb is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU Lesser General Public License for more details.
!
! You should have received a copy of the GNU Lesser General Public License
! along with xtb. If not, see <https://www.gnu.org/licenses/>.

!> Intermolecular force field calculator
module xtb_iff_calculator
use xtb_mctc_accuracy, only: wp
use xtb_type_calculator, only: TCalculator
use xtb_type_environment, only: TEnvironment
use xtb_type_molecule, only: TMolecule, init
use xtb_type_restart
use xtb_iff_data, only: TIFFData
use xtb_type_data, only: scc_results
use xtb_iff_iffenergy, only: iff_e
use xtb_docking_param
use xtb_iff_iffini, only : init_iff

implicit none

private

public :: TIFFCalculator, newIFFCalculator

!> Calculator interface for xTB based methods
type, extends(TCalculator) :: TIFFCalculator

type(TIFFData) :: dat

contains

!> Perform xTB single point calculationV
procedure :: singlepoint

!> Write informative printout
procedure :: writeInfo

end type TIFFCalculator

character(len=*), private, parameter :: outfmt = &
'(9x,"::",1x,a,f23.12,1x,a,1x,"::")'

contains

subroutine newIFFCalculator(env, comb, iff_data, calc)

character(len=*), parameter :: source = 'main_setup_newIFFCalculator'

type(TEnvironment), intent(inout) :: env

!> Combined structure of molA and molB (molA has to be first)
type(TMolecule), intent(in) :: comb

type(TIFFData), intent(in) :: iff_data

type(TIFFCalculator), intent(out) :: calc

integer, allocatable :: at(:)
real(wp), allocatable :: xyz(:,:)
logical :: exitRun
real(wp) :: molA_e,molB_e
character(len=:), allocatable :: fnam
real(wp) :: icoord(6), icoord0(6)
real(wp) :: rlmo2(4,(comb%n-natom_molA)*10)

call set_iff_param
call calc%dat%allocateIFFData(natom_molA,comb%n-natom_molA)
calc%dat = iff_data

xyz=calc%dat%xyz2
rlmo2=calc%dat%rlmo2

call init_iff(env,calc%dat%n1,calc%dat%n2,calc%dat%at1,calc%dat%at2,&
& calc%dat%neigh,calc%dat%xyz1,calc%dat%xyz2,calc%dat%q1,&
& calc%dat%q2,calc%dat%c6ab,calc%dat%z1,calc%dat%z2,&
& calc%dat%cprob,calc%dat%nlmo1,calc%dat%nlmo2,calc%dat%lmo1,calc%dat%lmo2,&
& calc%dat%qdr1,calc%dat%qdr2,calc%dat%rlmo1,calc%dat%rlmo2,&
& calc%dat%cn1,calc%dat%cn2,calc%dat%alp1,calc%dat%alp2,calc%dat%alpab,&
& calc%dat%den1,calc%dat%den2,calc%dat%gab1,calc%dat%gab2,calc%dat%qcm1,&
& calc%dat%qcm2,calc%dat%n,calc%dat%at,calc%dat%xyz,calc%dat%q,icoord,icoord0,&
& .false.)
calc%dat%xyz2=xyz
calc%dat%rlmo2=rlmo2

call env%check(exitRun)
if (exitRun) then
call env%error("Could not create IFF calculator", source)
return
end if

end subroutine newIFFCalculator

subroutine singlepoint(self, env, mol, chk, printlevel, restart, &
& energy, gradient, sigma, hlgap, results)

!> Source of the generated errors
character(len=*), parameter :: source = 'iff_calculator_singlepoint'

!> Calculator instance
class(TIFFCalculator), intent(inout) :: self

!> Computational environment
type(TEnvironment), intent(inout) :: env

!> Molecular structure data
type(TMolecule), intent(inout) :: mol

!> Wavefunction data
type(TRestart), intent(inout) :: chk

!> Print level for IO
integer, intent(in) :: printlevel

!> Restart from previous results
logical, intent(in) :: restart

!> Total energy
real(wp), intent(out) :: energy

!> Molecular gradient
real(wp), intent(out) :: gradient(:, :)

!> Strain derivatives
real(wp), intent(out) :: sigma(:, :)

!> HOMO-LUMO gab
real(wp), intent(out) :: hlgap

!> Detailed results
type(scc_results), intent(out) :: results

call iff_e(env, self%dat%n, self%dat%n1, self%dat%n2, self%dat%at1,&
& self%dat%at2, self%dat%neigh, self%dat%xyz1, self%dat%xyz2,&
& self%dat%q1, self%dat%q2, self%dat%c6ab, self%dat%z1, self%dat%z2,&
& self%dat%nlmo1, self%dat%nlmo2, self%dat%lmo1, self%dat%lmo2,&
& self%dat%rlmo1, self%dat%rlmo2,self%dat%qdr1, self%dat%qdr2,&
& self%dat%cn1, self%dat%cn2, self%dat%alp1,&
& self%dat%alp2, self%dat%alpab, self%dat%qct1, self%dat%qct2,&
& self%dat%den1, self%dat%den2, self%dat%gab1, self%dat%gab2,&
& .true., 0, energy)

results%e_total=energy
gradient = 0
sigma = 0
hlgap = 0

end subroutine singlepoint

subroutine writeInfo(self, unit, mol)

!> Calculator instance
class(TIFFCalculator), intent(in) :: self

!> Unit for I/O
integer, intent(in) :: unit

!> Molecular structure data
type(TMolecule), intent(in) :: mol

end subroutine writeInfo

end module xtb_iff_calculator
191 changes: 191 additions & 0 deletions src/iff/data.f90
Original file line number Diff line number Diff line change
@@ -0,0 +1,191 @@
! This file is part of xtb.
!
! Copyright (C) 2022 Christoph Plett
!
! xtb is free software: you can redistribute it and/or modify it under
! the terms of the GNU Lesser General Public License as published by
! the Free Software Foundation, either version 3 of the License, or
! (at your option) any later version.
!
! xtb is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU Lesser General Public License for more details.
!
! You should have received a copy of the GNU Lesser General Public License
! along with xtb. If not, see <https://www.gnu.org/licenses/>.

!> Topological data for force field type calculations
module xtb_iff_data
use xtb_mctc_accuracy, only: wp
use xtb_type_molecule, only: TMolecule
implicit none
private

public :: TIFFData

type :: TIFFData

!> Total number of atoms
integer :: n

!> Number of atoms molA and B
integer :: n1
integer :: n2

!> Number of LMOs molA and B
integer :: nlmo1
integer :: nlmo2

!> Coordinates molA and B and whole system
real(wp), allocatable :: xyz1(:, :)
real(wp), allocatable :: xyz2(:, :)
real(wp), allocatable :: xyz(:, :)

!> LMO positions molA and B
real(wp), allocatable :: rlmo1(:, :)
real(wp), allocatable :: rlmo2(:, :)

!> Charges
real(wp), allocatable :: q1(:)
real(wp), allocatable :: q2(:)
real(wp), allocatable :: q(:)

!> Deviation from atomic charge
real(wp), allocatable :: qdr1(:)
real(wp), allocatable :: qdr2(:)

!> Deviation from atomic position
real(wp), allocatable ::xyzdr1(:, :)
real(wp), allocatable ::xyzdr2(:, :)

!> Coordination numbers
real(wp), allocatable :: cn1(:)
real(wp), allocatable :: cn2(:)
real(wp), allocatable :: cn(:)

!> Electron count per element
real(wp), allocatable :: z1(:)
real(wp), allocatable :: z2(:)

!> Polarizabilities
real(wp), allocatable :: alp1(:)
real(wp), allocatable :: alp2(:)
real(wp), allocatable :: alpab(:, :)
real(wp), allocatable :: alp0(:)

!> Ordinary number
integer, allocatable :: at1(:)
integer, allocatable :: at2(:)
integer, allocatable :: at(:)

!> LMO values
integer, allocatable :: lmo1(:)
integer, allocatable :: lmo2(:)

!> C6 coeffs for the whole system
real(wp), allocatable :: c6ab(:, :)

!> C6 coeffs for molA with rare gas
real(wp), allocatable :: cprob(:)

real(wp), allocatable :: gab1(:, :)
real(wp), allocatable :: gab2(:, :)

!> Frontier orbital densities
real(wp), allocatable :: den1(:, :, :)
real(wp), allocatable :: den2(:, :, :)

!> Charge related stuff (from GFN2 calc.)
real(wp), allocatable :: qcm1(:)
real(wp), allocatable :: qcm2(:)
real(wp), allocatable :: qct1(:, :)
real(wp), allocatable :: qct2(:, :)

!> Neighbour list
integer, allocatable :: neigh(:, :)

contains

procedure :: delete
procedure :: allocateIFFData

end type TIFFData

contains

subroutine delete(self)
class(TIFFData), intent(out) :: self

if (allocated(self%xyz1)) deallocate (self%xyz1)
if (allocated(self%rlmo1)) deallocate (self%rlmo1)
if (allocated(self%q1)) deallocate (self%q1)
if (allocated(self%qdr1)) deallocate (self%qdr1)
if (allocated(self%xyzdr1)) deallocate (self%xyzdr1)
if (allocated(self%cn1)) deallocate (self%cn1)
if (allocated(self%z1)) deallocate (self%z1)
if (allocated(self%alp1)) deallocate (self%alp1)
if (allocated(self%qct1)) deallocate (self%qct1)
if (allocated(self%at1)) deallocate (self%at1)
if (allocated(self%lmo1)) deallocate (self%lmo1)
if (allocated(self%xyz2)) deallocate (self%xyz2)
if (allocated(self%rlmo2)) deallocate (self%rlmo2)
if (allocated(self%q2)) deallocate (self%q2)
if (allocated(self%qdr2)) deallocate (self%qdr2)
if (allocated(self%xyzdr2)) deallocate (self%xyzdr2)
if (allocated(self%cn2)) deallocate (self%cn2)
if (allocated(self%z2)) deallocate (self%z2)
if (allocated(self%alp2)) deallocate (self%alp2)
if (allocated(self%qct2)) deallocate (self%qct2)
if (allocated(self%at2)) deallocate (self%at2)
if (allocated(self%lmo2)) deallocate (self%lmo2)
if (allocated(self%c6ab)) deallocate (self%c6ab)
if (allocated(self%alpab)) deallocate (self%alpab)
if (allocated(self%cprob)) deallocate (self%cprob)
if (allocated(self%xyz)) deallocate (self%xyz)
if (allocated(self%q)) deallocate (self%q)
if (allocated(self%cn)) deallocate (self%cn)
if (allocated(self%alp0)) deallocate (self%alp0)
if (allocated(self%gab1)) deallocate (self%gab1)
if (allocated(self%gab2)) deallocate (self%gab2)
if (allocated(self%den1)) deallocate (self%den1)
if (allocated(self%den2)) deallocate (self%den2)
if (allocated(self%qcm1)) deallocate (self%qcm1)
if (allocated(self%qcm2)) deallocate (self%qcm2)
if (allocated(self%at)) deallocate (self%at)
if (allocated(self%neigh)) deallocate (self%neigh)

end subroutine delete

subroutine allocateIFFData(self, n1, n2)

class(TIFFData), intent(out) :: self

integer, intent(in) :: n1, n2

integer :: n

self%n1 = n1
self%n2 = n2

!Allocate Infos of molA
allocate(self%at1(n1),self%xyz1(3,n1),self%rlmo1(4,10*n1),self%q1(n1),self%lmo1(10*n1),&
&self%cn1(n1),self%alp1(n1),self%qct1(n1,2),self%qdr1(n1),self%xyzdr1(3,n1),&
&self%z1(n1), self%den1(2, 4, n1), self%gab1(n1, n1), self%qcm1(n1))
self%rlmo1 = 0.0_wp

!Allocate Infos of molB
allocate(self%at2(n2),self%xyz2(3,n2),self%rlmo2(4,10*n2),self%q2(n2),self%lmo2(10*n2),&
& self%cn2(n2),self%alp2(n2),self%qct2(n2,2),self%qdr2(n2),self%xyzdr2(3,n2),&
& self%z2(n2), self%den2(2, 4, n2), self%gab2(n2, n2), self%qcm2(n2))
self%rlmo2 = 0.0_wp

!Allocate combined Infos
self%n = n1 + n2
n = self%n
allocate(self%at(n),self%xyz(3,n),self%q(n),self%c6ab(n,n),self%alp0(n),self%cn(n),self%neigh(0:n,n),&
& self%alpab(n2, n1), self%cprob(n1))

end subroutine allocateIFFData

end module xtb_iff_data
972 changes: 972 additions & 0 deletions src/iff/iff_energy.f90

Large diffs are not rendered by default.

598 changes: 598 additions & 0 deletions src/iff/iff_ini.f90

Large diffs are not rendered by default.

591 changes: 591 additions & 0 deletions src/iff/iff_lmo.f90

Large diffs are not rendered by default.

241 changes: 241 additions & 0 deletions src/iff/iff_prepare.f90
Original file line number Diff line number Diff line change
@@ -0,0 +1,241 @@
! This file is part of xtb.
!
! Copyright (C) 2022 Christoph Plett
!
! xtb is free software: you can redistribute it and/or modify it under
! the terms of the GNU Lesser General Public License as published by
! the Free Software Foundation, either version 3 of the License, or
! (at your option) any later version.
!
! xtb is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU Lesser General Public License for more details.
!
! You should have received a copy of the GNU Lesser General Public License
! along with xtb. If not, see <https://www.gnu.org/licenses/>.

module xtb_iff_iffprepare
use xtb_mctc_accuracy, only: wp
use xtb_type_environment, only: TEnvironment
use xtb_splitparam
use xtb_iff_calculator, only: TIFFCalculator
use xtb_iff_data, only: TIFFData
use xtb_type_atomlist
use xtb_type_molecule, only: TMolecule
use xtb_setparam
use xtb_type_restart, only: TRestart
use xtb_type_calculator, only: TCalculator
use xtb_xtb_calculator, only: TxTBCalculator
use xtb_setmod, only: set_gfn
use xtb_readin, only: xfind
use xtb_solv_state, only: solutionState
use xtb_type_data, only: scc_results
use xtb_main_defaults, only: initDefaults
use xtb_disp_ncoord, only: ncoord_gfn, ncoord_erf
use xtb_scc_core, only: iniqshell
use xtb_eeq, only: goedecker_chrgeq
use xtb_main_setup, only: newCalculator
use xtb_single, only: singlepoint
use xtb_docking_param, only: chrg, uhf, gsolvstate_iff, pre_e_A, &
& pre_e_B, optlvl, ehomo, elumo, dipol, &
& natom_molA, natom_arg, split_mol

implicit none

private
public :: prepare_IFF, precomp

contains

subroutine prepare_IFF(env, comb, iff_data)

type(TEnvironment), intent(inout) :: env
!> Combined structure of molA and molB (molA has to be first)
type(TMolecule), intent(in) :: comb
!> IFF data
type(TIFFData) :: iff_data

character(len=*), parameter :: source = 'preperation_IFFCalculator'
type(TMolecule) :: molA, molB
integer, allocatable :: at(:)
real(wp), allocatable :: xyz(:, :)
real(wp) :: molA_e, molB_e
integer, allocatable :: list(:)
type(TAtomList) :: atl

call atl%resize(comb%n)

!> First make the argument natom_arg to a list of number of atoms
call atl%new(natom_arg)
if (atl%get_error()) then
call env%warning('something is wrong in the fixing list',source)
return
endif
call atl%to_list(list)

molA%n = size(list)
natom_molA = molA%n
molB%n = comb%n - molA%n

if (natom_molA == 0) call env%error('No atoms of Molecule A given')
if(natom_molA > comb%n) &
& call env%error('More atoms of Molecule A than contained in structure')

call split_mol(molA, molB, size(list), list, comb)
call iff_data%allocateIFFData(molA%n, molB%n)

call precomp(env, iff_data, molA, molA_e, 1)
call precomp(env, iff_data, molB, molB_e, 2)

end subroutine prepare_IFF

subroutine precomp(env, iff_data, mol, etot, mol_num)

!> Molecular structure data
type(TMolecule), intent(inout) :: mol
!> IFF data
type(TIFFData), intent(inout) :: iff_data
!> Calculation environment
type(TEnvironment), intent(inout) :: env
integer, intent(in) :: mol_num

class(TCalculator), allocatable :: calc
type(TRestart) :: chk
character(len=:), allocatable :: fnv !parameter file
logical :: restart = .false.
real(wp) :: acc = 1.0_wp !SCF accuracy
real(wp) :: egap
integer :: gsolvstate !Should be gas
real(wp), allocatable :: cn(:)
real(wp), allocatable :: dcn(:, :, :), dq(:, :, :), g(:, :)
real(wp) :: er
integer, external :: ncore
logical :: lgrad
logical :: exist
real(wp), intent(out) :: etot
real(wp) :: sigma(3, 3)
type(scc_results) :: res
integer :: i
integer :: extrun_tmp
integer, allocatable :: tmp_unit
logical :: newdisp_tmp
integer :: itemp = 48

allocate (cn(mol%n), g(3, mol%n))

set%pr_lmo = .true.
set%silent = .true.

!Setting up z
do i = 1, mol%n
mol%z(i) = mol%at(i) - ncore(mol%at(i))
! lanthanides without f are treated as La
if (mol%at(i) .gt. 57 .and. mol%at(i) .lt. 72) mol%z(i) = 3
end do

!> Set GFN2 settings
if (optlvl /= 'gfn1') then
set%gfn_method = 2
fnv = xfind('param_gfn2-xtb.txt')
extrun_tmp = set%mode_extrun
set%mode_extrun = p_ext_xtb
newdisp_tmp = set%newdisp
set%newdisp = .true.
else
fnv = xfind('param_gfn1-xtb.txt') !Other stuff already set
end if
if (mol_num .eq. 1) then
if (chrg(1) /= 0.0_wp) mol%chrg = chrg(1)
if (uhf(1) /= 0.0_wp) mol%uhf = uhf(1)
elseif (mol_num .eq. 2) then
if (chrg(2) /= 0.0_wp) mol%chrg = chrg(2)
if (uhf(2) /= 0.0_wp) mol%uhf = uhf(2)
end if

call open_file(itemp, 'tmp_lmo', 'w')
tmp_unit = env%unit
env%unit = itemp

!> New calculator
call newCalculator(env, mol, calc, fnv, restart, acc)
call env%checkpoint("Could not setup parameterisation")

! gsolvstate = solutionState%gsolv
call initDefaults(env, calc, mol, gsolvstate_iff)
!> initial guess, setup wavefunction
select type (calc)
type is (TxTBCalculator)
calc%etemp = set%etemp
calc%maxiter = set%maxscciter
call chk%wfn%allocate(mol%n, calc%basis%nshell, calc%basis%nao)
!> EN charges and CN
call ncoord_gfn(mol%n, mol%at, mol%xyz, cn)
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)
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)
end select

call delete_file('.sccnotconverged')
call delete_file('charges')

call env%checkpoint("Setup for calculation failed")

allocate(res%iff_results)

!> the SP
call singlepoint &
& (env, mol, chk, calc, egap, set%etemp, set%maxscciter, 2,&
& exist, lgrad, acc, etot, g, sigma, res)

set%pr_lmo = .false.

!> Save the results
if (mol_num .eq. 1) then
pre_e_A = etot
iff_data%n1 = res%iff_results%n
iff_data%at1 = res%iff_results%at
iff_data%xyz1 = res%iff_results%xyz
iff_data%q1 = res%iff_results%q
iff_data%nlmo1 = res%iff_results%nlmo
iff_data%rlmo1 = res%iff_results%rlmo
iff_data%lmo1 = res%iff_results%lmo
iff_data%qct1 = res%iff_results%qct
ehomo(1) = res%iff_results%ehomo
elumo(1) = res%iff_results%elumo
dipol(1) = res%iff_results%dipol
elseif (mol_num .eq. 2) then
pre_e_B = etot
iff_data%n2 = res%iff_results%n
iff_data%at2 = res%iff_results%at
iff_data%xyz2 = res%iff_results%xyz
iff_data%q2 = res%iff_results%q
iff_data%nlmo2 = res%iff_results%nlmo
iff_data%rlmo2 = res%iff_results%rlmo
iff_data%lmo2 = res%iff_results%lmo
iff_data%qct2 = res%iff_results%qct
ehomo(2) = res%iff_results%ehomo
elumo(2) = res%iff_results%elumo
dipol(2) = res%iff_results%dipol
end if

if (optlvl /= 'gfn1') then
set%mode_extrun = extrun_tmp
set%newdisp = newdisp_tmp
end if

set%silent = .false.
env%unit = tmp_unit
call remove_file(itemp)
deallocate (cn, g)

end subroutine precomp

end module xtb_iff_iffprepare
25 changes: 25 additions & 0 deletions src/iff/meson.build
Original file line number Diff line number Diff line change
@@ -0,0 +1,25 @@
# This file is part of xtb.
#
# Copyright (C) 2022 Sebastian Ehlert, Christoph Plett
#
# xtb is free software: you can redistribute it and/or modify it under
# the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# xtb is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with xtb. If not, see <https://www.gnu.org/licenses/>.

srcs += files(
'calculator.f90',
'data.f90',
'iff_energy.f90',
'iff_ini.f90',
'iff_lmo.f90',
'iff_prepare.f90',
)
223 changes: 131 additions & 92 deletions src/local.f90
Original file line number Diff line number Diff line change
@@ -17,7 +17,7 @@
module xtb_local
contains

subroutine local(nat,at,nbf,nao,ihomo,xyz,z,focc,s,p,cmo,eig,q,etot,gbsa,basis)
subroutine local(nat,at,nbf,nao,ihomo,xyz,z,focc,s,p,cmo,eig,q,etot,gbsa,basis,results)
use xtb_mctc_accuracy, only : wp, sp
use xtb_mctc_constants, only : pi
use xtb_mctc_convert, only : autoev,autoaa
@@ -29,13 +29,17 @@ subroutine local(nat,at,nbf,nao,ihomo,xyz,z,focc,s,p,cmo,eig,q,etot,gbsa,basis)
use xtb_dtrafo
use xtb_onetri
use xtb_dipole
use xtb_type_data, only : scc_results
use xtb_docking_param, only : dipol, ehomo, elumo

implicit none
type(TBasisset), intent(in) :: basis
integer, intent(in) :: nao,ihomo,nat,at(nat),nbf
real(wp),intent(in) :: cmo(nao,nao),eig(nao),focc(nao)
real(wp),intent(in) :: s(nao,nao),xyz(3,nat),z(nat),q(nat)
real(wp),intent(in) :: p(nao,nao),etot
logical, intent(in) :: gbsa
type(scc_results), intent(inout) :: results

real(wp),allocatable :: op(:,:)
real(wp),allocatable :: oc(:,:,:)
@@ -64,7 +68,7 @@ subroutine local(nat,at,nbf,nao,ihomo,xyz,z,focc,s,p,cmo,eig,q,etot,gbsa,basis)
integer i1,i2,i3,is3,ipi,piset(nat),ncyc,j1,j2,smo,pmo,nl,nm,nr
integer imo,new
real(wp) :: dd,dum,det,pp,dtot(3),diptot,thr,t0,w0,t1,w1,vec1(3),v(6)
real(wp) :: elumo,ehomo,qhl(nat,2),ps,efh,efl,r1,r2,pithr,vec2(3)
real(wp) :: enlumo,enhomo,qhl(nat,2),ps,efh,efl,r1,r2,pithr,vec2(3)
real(wp) :: praxis(3,3),aa,bb,cc
character(len=80) :: atmp
character(len=5) :: lmostring(4)
@@ -74,34 +78,36 @@ subroutine local(nat,at,nbf,nao,ihomo,xyz,z,focc,s,p,cmo,eig,q,etot,gbsa,basis)
integer :: iscreen,icoord,ilmoi,icent ! file handles

call timing(t0,w0)
write(*,*)
write(*,*)'localization/xTB-IFF output generation'
if(set%pr_local) then
write(*,*)
write(*,*)'localization/xTB-IFF output generation'
end if
n=ihomo
ilumo=n+1

if(ilumo.gt.nao)then
elumo=1000.
enlumo=1000.
else
elumo=eig(ilumo)
enlumo=eig(ilumo)
endif
efh=eig(n)
efl=elumo
efl=enlumo

! HOMO/LUMO populations for xTB FF (CT terms)
qhl = 0
thr = 0.01

if(ihomo.eq.0) then
! the HOMO is non-existing, place at unrealistically low energy
ehomo=-999.999999990d0
enhomo=-999.999999990d0
qhl(1:nat,1)=0.0d0
else
klev=0
ehomo=0
enhomo=0
do nlev=ihomo,1,-1 ! loop over highest levels
if(efh-eig(nlev).gt.thr) exit
klev=klev+1
ehomo=ehomo+eig(nlev)
enhomo=enhomo+eig(nlev)
do i=1,nao
ii=basis%aoat2(i)
do j=1,i-1
@@ -115,17 +121,17 @@ subroutine local(nat,at,nbf,nao,ihomo,xyz,z,focc,s,p,cmo,eig,q,etot,gbsa,basis)
enddo
enddo
dum=1./float(klev)
ehomo=ehomo*dum
enhomo=enhomo*dum
qhl(1:nat,1)=qhl(1:nat,1)*dum
write(*,*)'averaging CT terms over ',klev,' occ. levels'
if(set%pr_local) write(*,*)'averaging CT terms over ',klev,' occ. levels'
endif

klev=0
elumo=0
enlumo=0
do nlev=ilumo,nao ! loop over highest levels
if(eig(nlev)-efl.gt.thr) exit
klev=klev+1
elumo=elumo+eig(nlev)
enlumo=enlumo+eig(nlev)
do i=1,nao
ii=basis%aoat2(i)
do j=1,i-1
@@ -139,11 +145,11 @@ subroutine local(nat,at,nbf,nao,ihomo,xyz,z,focc,s,p,cmo,eig,q,etot,gbsa,basis)
enddo
enddo
dum=1./float(klev)
elumo=elumo*dum
enlumo=enlumo*dum
qhl(1:nat,2)=qhl(1:nat,2)*dum
write(*,*)'averaging CT terms over ',klev,' virt. levels'
if(set%pr_local) write(*,*)'averaging CT terms over ',klev,' virt. levels'
if(ilumo.gt.nao)then
elumo=1000.
enlumo=1000.
qhl(1:nat,2)=0
endif

@@ -201,18 +207,20 @@ subroutine local(nat,at,nbf,nao,ihomo,xyz,z,focc,s,p,cmo,eig,q,etot,gbsa,basis)
dtot(3)=dtot(3)+op(k,3)*focc(i)
enddo
diptot=sqrt(dtot(1)**2+dtot(2)**2+dtot(3)**2)
write(*,*)'dipole moment from electron density (au)'
write(*,*)' X Y Z '
write(*,'(3f9.4,'' total (Debye): '',f8.3)') &
& dtot(1),dtot(2),dtot(3),diptot*2.5418
if(set%pr_local) then
write(*,*)'dipole moment from electron density (au)'
write(*,*)' X Y Z '
write(*,'(3f9.4,'' total (Debye): '',f8.3)') &
& dtot(1),dtot(2),dtot(3),diptot*2.5418
end if

call timing(t1,w1)
call prtime(6,t1-t0,w1-w0,'init local')
if(set%pr_local) call prtime(6,t1-t0,w1-w0,'init local')
! do optimization
write(*,*) 'doing rotations ...'
if(set%pr_local) write(*,*) 'doing rotations ...'
call lopt(.true.,n,3,1.d-6,op,d)

write(*,*) 'doing transformations ...'
if(set%pr_local) write(*,*) 'doing transformations ...'
! MO charge center
k=0
do i=1,n
@@ -265,9 +273,9 @@ subroutine local(nat,at,nbf,nao,ihomo,xyz,z,focc,s,p,cmo,eig,q,etot,gbsa,basis)
call mocent(nat,nao,n,cmo,s,qmo,xcen,basis%aoat2)

allocate(rklmo(5,2*n))
write(*,*) 'lmo centers(Z=2) and atoms on file <lmocent.coord>'
write(*,*) 'LMO Fii/eV ncent charge center contributions...'
call open_file(iscreen,'xtbscreen.xyz','w')
if(set%pr_local) write(*,*) 'lmo centers(Z=2) and atoms on file <lmocent.coord>'
if(set%pr_local) write(*,*) 'LMO Fii/eV ncent charge center contributions...'
if(set%pr_local) call open_file(iscreen,'xtbscreen.xyz','w')

allocate(tmpq(nat,n))
tmpq(1:nat,1:n)=qmo(1:nat,1:n)
@@ -313,32 +321,34 @@ subroutine local(nat,at,nbf,nao,ihomo,xyz,z,focc,s,p,cmo,eig,q,etot,gbsa,basis)
call lmotype(nat,at,xyz,ecent(i,1),ecent(i,2),ecent(i,3), &
& imem(1),imem(2),xcen(i),.true.,pithr,jdum)
endif
write(*,'(i5,1x,a5,2f7.2,3f10.5,12(i5,a2,'':'',f6.2))') &
if(set%pr_local) write(*,'(i5,1x,a5,2f7.2,3f10.5,12(i5,a2,'':'',f6.2))') &
& i,lmostring(jdum),autoev*f(i),xcen(i),ecent(i,1:3), &
& (imem(j),toSymbol(at(imem(j))),qmo(j,i),j=1,idum)

! write + LP/pi as H for protonation search
if(jdum.gt.1) then
if( i.eq.maxlp .or. (i.eq.maxpi.and.maxlp.eq.0) )then
call open_file(icoord,'coordprot.0','w')
write(icoord,'(''$coord'')')
do ii=1,nat
write(icoord,'(3F24.10,5x,a2)') xyz(1:3,ii),toSymbol(at(ii))
enddo
write(icoord,'(3F24.10,5x,a2)') ecent(i,1:3),toSymbol(1)
write(icoord,'(''$end'')')
write(icoord,'(''$set'')')
write(icoord,'('' chrg '',i2)')set%ichrg+1
write(icoord,'('' ewin_conf 50.0 '')')
write(icoord,'(''$end'')')
call close_file(icoord)
else
write(iscreen,*) nat+1
write(iscreen,*)
do ii=1,nat
write(iscreen,'(a2,3F24.10)')toSymbol(at(ii)),xyz(1:3,ii)*autoaa
enddo
write(iscreen,'(a2,3F24.10)') toSymbol(1),ecent(i,1:3)*autoaa
if(set%pr_local) then
if(jdum.gt.1) then
if( i.eq.maxlp .or. (i.eq.maxpi.and.maxlp.eq.0) )then
call open_file(icoord,'coordprot.0','w')
write(icoord,'(''$coord'')')
do ii=1,nat
write(icoord,'(3F24.10,5x,a2)') xyz(1:3,ii),toSymbol(at(ii))
enddo
write(icoord,'(3F24.10,5x,a2)') ecent(i,1:3),toSymbol(1)
write(icoord,'(''$end'')')
write(icoord,'(''$set'')')
write(icoord,'('' chrg '',i2)')set%ichrg+1
write(icoord,'('' ewin_conf 50.0 '')')
write(icoord,'(''$end'')')
call close_file(icoord)
else
write(iscreen,*) nat+1
write(iscreen,*)
do ii=1,nat
write(iscreen,'(a2,3F24.10)')toSymbol(at(ii)),xyz(1:3,ii)*autoaa
enddo
write(iscreen,'(a2,3F24.10)') toSymbol(1),ecent(i,1:3)*autoaa
endif
endif
endif

@@ -348,7 +358,7 @@ subroutine local(nat,at,nbf,nao,ihomo,xyz,z,focc,s,p,cmo,eig,q,etot,gbsa,basis)

new=0
if(maxval(rklmo(5,1:n)).ge.3.)then
write(*,*) 'starting deloc pi regularization ...'
if(set%pr_local) write(*,*) 'starting deloc pi regularization ...'

call atomneigh(n,nat,xyz,ecent,aneigh) ! the nearest atom neighbors for each LMO
call lmoneigh (n,rklmo,ecent,aneigh,lneigh) ! the nearest LMO neighbor but not LP
@@ -417,7 +427,7 @@ subroutine local(nat,at,nbf,nao,ihomo,xyz,z,focc,s,p,cmo,eig,q,etot,gbsa,basis)
enddo
enddo

write(*,*) 'thr ',pithr, '# pi deloc LMO',npi
if(set%pr_local) write(*,*) 'thr ',pithr, '# pi deloc LMO',npi


allocate(wbo(nat,nat))
@@ -454,12 +464,14 @@ subroutine local(nat,at,nbf,nao,ihomo,xyz,z,focc,s,p,cmo,eig,q,etot,gbsa,basis)
rklmo(1:3,n+k)= dtot(1:3)
rklmo(4:5,n+k)=rklmo(4:5,pmo)
! add to screen file, protomer search
write(iscreen,*) nat+1
write(iscreen,*)
do ii=1,nat
write(iscreen,'(a2,3F24.10)')toSymbol(at(ii)),xyz(1:3,ii)*autoaa
enddo
write(iscreen,'(a2,3F24.10)') toSymbol(1),dtot(1:3)*autoaa
if(set%pr_local) then
write(iscreen,*) nat+1
write(iscreen,*)
do ii=1,nat
write(iscreen,'(a2,3F24.10)')toSymbol(at(ii)),xyz(1:3,ii)*autoaa
enddo
write(iscreen,'(a2,3F24.10)') toSymbol(1),dtot(1:3)*autoaa
endif
imo=lneigh(1,pmo)
vec1(1:3)=xyz(1:3,nm)
vec2(1:3)=(xyz(1:3,nl)+xyz(1:3,nr))*0.5
@@ -470,59 +482,86 @@ subroutine local(nat,at,nbf,nao,ihomo,xyz,z,focc,s,p,cmo,eig,q,etot,gbsa,basis)
rklmo(4:5,n+k)=rklmo(4:5,imo)
rklmo( 5,smo)=0 ! remove sigma
! add to screen file, protomer search
write(iscreen,*) nat+1
write(iscreen,*)
do ii=1,nat
write(iscreen,'(a2,3F24.10)')toSymbol(at(ii)),xyz(1:3,ii)*autoaa
enddo
write(iscreen,'(a2,3F24.10)') toSymbol(1),dtot(1:3)*autoaa
if(set%pr_local) then
write(iscreen,*) nat+1
write(iscreen,*)
do ii=1,nat
write(iscreen,'(a2,3F24.10)')toSymbol(at(ii)),xyz(1:3,ii)*autoaa
enddo
write(iscreen,'(a2,3F24.10)') toSymbol(1),dtot(1:3)*autoaa
end if
m=m+1
enddo
enddo
new=k

call close_file(iscreen)
if(set%pr_local) call close_file(iscreen)

deallocate(wbo)


endif
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccc

! all LMO on file
call open_file(ilmoi,'xtblmoinfo','w')
call open_file(icent,'lmocent.coord','w')
write(icent,'(''$coord'')')
write(ilmoi,*) nat
do i=1,nat
write(ilmoi,'(i2,3F20.10,E18.10)') at(i),xyz(1:3,i),q(i)
write(icent,'(3F24.10,5x,a2)') xyz(1:3,i),toSymbol(at(i))
enddo
!> If the normal xtb mode is used with --lmo, set%pr_local is true and
! a file with all information is written.
! If the docking mode is used, information are stored intrenally

if(set%pr_local) then
call open_file(ilmoi,set%lmoinfo_fname,'w')
write(ilmoi,*) nat
if(set%pr_local) call open_file(icent,'lmocent.coord','w')
if(set%pr_local) write(icent,'(''$coord'')')
do i=1,nat
write(ilmoi,'(i2,3F20.10,E18.10)') at(i),xyz(1:3,i),q(i)
if(set%pr_local) write(icent,'(3F24.10,5x,a2)') xyz(1:3,i),toSymbol(at(i))
enddo
k=0
do i=1,n+new
if(int(rklmo(5,i)).gt.0)k=k+1
enddo
write(ilmoi,'(i5,5f14.8)')k,diptot,enlumo,enhomo,etot
do i=1,n+new
if(int(rklmo(5,i)).gt.0)then
write(ilmoi,'(i2,3F20.10,f14.8)')int(rklmo(5,i)),rklmo(1:3,i)
if(set%pr_local) write(icent,'(3F24.10,5x,a2)') rklmo(1:3,i),toSymbol(2)
endif
enddo
write(ilmoi,'(10(F10.6))') (qhl(i,1),i=1,nat) ! HOMO atom pop
write(ilmoi,'(10(F10.6))') (qhl(i,2),i=1,nat) ! LUMO atom pop
if(set%pr_local) write(icent,'(''$end'')')
call close_file(ilmoi)
if(set%pr_local) call close_file(icent)
end if

!> Saving results
k=0
do i=1,n+new
if(int(rklmo(5,i)).gt.0)k=k+1
enddo
! dum=dble(m)/dble(m+npi)
! if(dum.lt.0.01) dum=1.0d0
write(ilmoi,'(i5,5f14.8)')k,diptot,elumo,ehomo,etot
results%iff_results%at = at
results%iff_results%xyz = xyz
results%iff_results%q = q
results%iff_results%nlmo = k
results%iff_results%dipol = diptot
results%iff_results%elumo = enlumo
results%iff_results%ehomo = enhomo
results%iff_results%qct(1:nat,1) = qhl(1:nat,1)
results%iff_results%qct(1:nat,2) = qhl(1:nat,2)
do i=1,n+new
if(int(rklmo(5,i)).gt.0)then
write(ilmoi,'(i2,3F20.10,f14.8)')int(rklmo(5,i)),rklmo(1:3,i)
write(icent,'(3F24.10,5x,a2)') rklmo(1:3,i),toSymbol(2)
if(int(rklmo(5,i)).gt.0) then
results%iff_results%lmo(i) = int(rklmo(5,i))
results%iff_results%rlmo(1:3,i) = rklmo(1:3,i)
endif
enddo
write(ilmoi,'(10(F10.6))') (qhl(i,1),i=1,nat) ! HOMO atom pop
write(ilmoi,'(10(F10.6))') (qhl(i,2),i=1,nat) ! LUMO atom pop
write(icent,'(''$end'')')
call close_file(ilmoi)
call close_file(icent)
! call camm(nat,at,nbf,nao,xyz,s,p)

write(*,*)'files:'
write(*,*)'coordprot.0/xtbscreen.xyz/xtblmoinfo/lmocent.coord'
write(*,*)'with protonation site input, xtbdock and'
write(*,*)'LMO center info written'
write(*,*)
end do

if(set%pr_local) then
write(*,*)'files:'
write(*,*)'coordprot.0/xtbscreen.xyz/xtblmoinfo/lmocent.coord'
write(*,*)'with protonation site input, xtbdock and'
write(*,*)'LMO center info written'
write(*,*)
end if

deallocate(xcen,cca,d,f,qmo,ecent,rklmo)

12 changes: 7 additions & 5 deletions src/lopt.f90
Original file line number Diff line number Diff line change
@@ -31,6 +31,8 @@ subroutine lopt(init, n, no, accr, op, d)
! sum(k) op(ij,k)*(op(ii,k)-op(jj,k)) = 0 all i.ne.j
!
use xtb_mctc_accuracy, only : wp
use xtb_setparam

implicit real(wp)(a-h,o-z)

logical init
@@ -75,7 +77,7 @@ subroutine lopt(init, n, no, accr, op, d)
! set d to unit matrix
!
if(init)then
write(*,*) 'initialization of trafo matrix to unity'
if(set%pr_local) write(*,*) 'initialization of trafo matrix to unity'
do i=1,n
do j=1,n
10 d(j,i)=a0-accr*i+accr*i*j
@@ -228,16 +230,16 @@ subroutine lopt(init, n, no, accr, op, d)
if (thsw.le.th) go to 210
ths=dmin1(thsw**2,ths*aqu)
ths=dmax1(th,ths)
if(mod(isw,50).eq.0)&
if(mod(isw,50).eq.0 .AND. set%pr_local)&
&write(*,'(''iteration'',i7,2x,''convergence'',d16.8)')isw,thsw
200 continue
end do
write (*,300) isw,thsw
if(set%pr_local) write (*,300) isw,thsw
return
210 write (*,310) isw,thsw
210 if(set%pr_local) write (*,310) isw,thsw

return
300 format (/' not converged in',i7,' iterations, threshold :',d16.8)
310 format (/' converged in',i7,' iterations, threshold : ',d16.8)
end subroutine
end subroutine lopt

26 changes: 24 additions & 2 deletions src/main/setup.f90
Original file line number Diff line number Diff line change
@@ -31,6 +31,8 @@ module xtb_main_setup
use xtb_type_wavefunction, only : TWavefunction
use xtb_xtb_calculator, only : TxTBCalculator, newXTBcalculator, newWavefunction
use xtb_gfnff_calculator, only : TGFFCalculator, newGFFCalculator
use xtb_iff_calculator, only : TIFFCalculator, newIFFCalculator
use xtb_iff_data, only : TIFFData
use xtb_oniom, only : TOniomCalculator, newOniomCalculator, oniom_input
use xtb_setparam
implicit none
@@ -43,7 +45,7 @@ module xtb_main_setup
contains


subroutine newCalculator(env, mol, calc, fname, restart, accuracy, input)
subroutine newCalculator(env, mol, calc, fname, restart, accuracy, input, iff_data)

character(len=*), parameter :: source = 'main_setup_newCalculator'

@@ -59,10 +61,13 @@ subroutine newCalculator(env, mol, calc, fname, restart, accuracy, input)

real(wp), intent(in) :: accuracy

type(oniom_input), intent(in) :: input
type(oniom_input), intent(in), optional :: input

type(TIFFData), intent(in), optional, allocatable :: iff_data

type(TxTBCalculator), allocatable :: xtb
type(TGFFCalculator), allocatable :: gfnff
type(TIFFCalculator), allocatable :: iff
type(TOrcaCalculator), allocatable :: orca
type(TMopacCalculator), allocatable :: mopac
type(TTMCalculator), allocatable :: turbo
@@ -98,6 +103,23 @@ subroutine newCalculator(env, mol, calc, fname, restart, accuracy, input)
end if

call move_alloc(gfnff, calc)
case(p_ext_iff)
allocate(iff)

if (.not. allocated(iff_data)) then
call env%error("IFF Data not present for Calculator", source)
end if

call newIFFCalculator(env, mol, iff_data, iff)

call env%check(exitRun)
if (exitRun) then
call env%error("Could not construct new calculator", source)
return
end if

call move_alloc(iff, calc)

case(p_ext_orca)
allocate(orca)
call newOrcaCalculator(orca, env, set%ext_orca)
2 changes: 2 additions & 0 deletions src/meson.build
Original file line number Diff line number Diff line change
@@ -18,10 +18,12 @@
subdir('api')
subdir('coulomb')
subdir('disp')
subdir('docking')
subdir('extern')
subdir('freq')
subdir('gfnff')
subdir('io')
subdir('iff')
subdir('main')
subdir('mctc')
subdir('param')
1 change: 1 addition & 0 deletions src/prog/CMakeLists.txt
Original file line number Diff line number Diff line change
@@ -19,6 +19,7 @@ set(dir "${CMAKE_CURRENT_SOURCE_DIR}")

list(APPEND prog
"${dir}/argparser.f90"
"${dir}/dock.f90"
"${dir}/info.f90"
"${dir}/main.F90"
"${dir}/primary.f90"
1,047 changes: 1,047 additions & 0 deletions src/prog/dock.f90

Large diffs are not rendered by default.

28 changes: 27 additions & 1 deletion src/prog/main.F90
Original file line number Diff line number Diff line change
@@ -67,6 +67,7 @@ module xtb_prog_main
use xtb_screening, only : screen
use xtb_xtb_calculator
use xtb_gfnff_calculator
use xtb_iff_calculator, only : TIFFCalculator
use xtb_paramset
use xtb_xtb_gfn0
use xtb_xtb_gfn1
@@ -85,6 +86,8 @@ module xtb_prog_main
use xtb_gfnff_convert, only : struc_convert
use xtb_scan
use xtb_kopt
use xtb_iff_iffprepare, only : prepare_IFF
use xtb_iff_data, only : TIFFData
use xtb_oniom, only : oniom_input, TOniomCalculator
implicit none
private
@@ -112,6 +115,7 @@ subroutine xtbMain(env, argParser)
type(freq_results) :: fres
type(TRestart) :: chk
type(chrg_parameter) :: chrgeq
type(TIFFData), allocatable :: iff_data
type(oniom_input) :: oniom
! store important names and stuff like that in FORTRAN strings
character(len=:),allocatable :: fname ! geometry input file
@@ -503,10 +507,17 @@ subroutine xtbMain(env, argParser)
end select
endif

!-------------------------------------------------------------------------
!> Perform a precomputation of electronic properties for xTB-IFF
if(set%mode_extrun == p_ext_iff) then
allocate(iff_data)
call prepare_IFF(env, mol, iff_data)
call env%checkpoint("Could not generate electronic properties")
end if

! ------------------------------------------------------------------------
!> Obtain the parameter data
call newCalculator(env, mol, calc, fnv, restart, acc, oniom)
call newCalculator(env, mol, calc, fnv, restart, acc, oniom, iff_data)
call env%checkpoint("Could not setup single-point calculator")

call initDefaults(env, calc, mol, gsolvstate)
@@ -1126,6 +1137,12 @@ subroutine parseArguments(env, args, inputFile, paramFile, accuracy, lgrad, &
!> Input for ONIOM model
type(oniom_input), intent(out) :: oniom

!> Stuff for second argument parser
! integer :: narg
! character(len=p_str_length), dimension(p_arg_length) :: argv
! type(TAtomList) :: atl
! integer, allocatable :: list(:)

!$ integer :: omp_get_num_threads, nproc
integer :: nFlags
integer :: idum, ndum
@@ -1297,6 +1314,9 @@ subroutine parseArguments(env, args, inputFile, paramFile, accuracy, lgrad, &
case('--gff')
call set_exttyp('ff')

case('--iff')
call set_exttyp('iff')

case('--oniom')
call set_exttyp('oniom')
call args%nextArg(sec) ! 'gfn2:gfnff'
@@ -1571,6 +1591,12 @@ subroutine parseArguments(env, args, inputFile, paramFile, accuracy, lgrad, &
call set_opt(env,'optlevel',sec)
end if

case('--nat')
call args%nextArg(sec)
if (allocated(sec)) then
call set_natom(env,sec)
end if

case('--bias-input', '--gesc')
call args%nextArg(sec)
if (allocated(sec)) then
1 change: 1 addition & 0 deletions src/prog/meson.build
Original file line number Diff line number Diff line change
@@ -17,6 +17,7 @@

prog += files(
'argparser.f90',
'dock.f90',
'info.f90',
'main.F90',
'primary.f90',
5 changes: 5 additions & 0 deletions src/prog/primary.f90
Original file line number Diff line number Diff line change
@@ -18,6 +18,7 @@
!> The primary driver for the xtb program package
program xtb_prog_primary
use xtb_prog_argparser
use xtb_prog_dock, only : xtbDock
use xtb_prog_info, only : xtbInfo
use xtb_prog_main, only : xtbMain
use xtb_prog_thermo, only : xtbThermo
@@ -66,6 +67,10 @@ program xtb_prog_primary
!> Run the thermo submodule
call xtbTopology(env, argParser)

case(xtbSubmodule%dock)
!> Run the docking submodule
call xtbDock(env, argParser)

end select

contains
6 changes: 6 additions & 0 deletions src/prog/submodules.f90
Original file line number Diff line number Diff line change
@@ -42,6 +42,9 @@ module xtb_prog_submodules
!> Force field topology generator
integer :: topo = 4

!> Fragment docking
integer :: dock = 5

end type TSubmoduleEnum

!> Actual enumerator for the submodules
@@ -73,6 +76,9 @@ function getSubmodule(argument) result(submod)
case('topo')
submod = xtbSubmodule%topo

case('dock')
submod = xtbSubmodule%dock

end select

end function getSubmodule
8 changes: 6 additions & 2 deletions src/scf_module.F90
Original file line number Diff line number Diff line change
@@ -838,9 +838,13 @@ subroutine scf(env, mol, wfn, basis, pcem, xtbData, solvation, &
! LMO /xTB-IFF
if (set%pr_lmo) then
tmp=wfn%emo*evtoau
allocate(res%iff_results)
! if(.not.set%pr_local) then
call res%iff_results%allocateIFFResults(mol%n)
! endif
call local(mol%n, mol%at, basis%nbf, basis%nao, wfn%ihomoa, mol%xyz, &
& mol%z, wfn%focc, S, wfn%P, wfn%C, tmp, wfn%q, eel, &
& allocated(solvation), basis)
& mol%z, wfn%focc, S, wfn%P, wfn%C, tmp, wfn%q, eel, &
& allocated(solvation), basis, res)
endif

endif printing
12 changes: 12 additions & 0 deletions src/set_module.f90
Original file line number Diff line number Diff line change
@@ -924,6 +924,8 @@ subroutine set_exttyp(typ)
set%mode_extrun = p_ext_mopac
case('ff')
set%mode_extrun = p_ext_gfnff
case('iff')
set%mode_extrun = p_ext_iff
end select
set1 = .false.
end subroutine set_exttyp
@@ -2386,6 +2388,16 @@ subroutine set_wall(env,key,val)

end subroutine set_wall

subroutine set_natom(env,arg)
use xtb_docking_param
implicit none
type(TEnvironment), intent(inout) :: env
character(len=*),intent(in) :: arg
integer :: idum
character(len=*), parameter :: source = 'set_natom'
natom_arg = arg
end subroutine set_natom

! this is a dummy routine
subroutine set_split(env,key,val)
use xtb_splitparam
3 changes: 3 additions & 0 deletions src/setparam.f90
Original file line number Diff line number Diff line change
@@ -118,6 +118,7 @@ module xtb_setparam
integer, parameter :: p_ext_mopac = 12
integer, parameter :: p_ext_gfnff = 13
integer, parameter :: p_ext_oniom = 14
integer, parameter :: p_ext_iff = 15

integer, parameter :: p_run_scc = 2
integer, parameter :: p_run_grad = 3
@@ -307,8 +308,10 @@ module xtb_setparam
character(len=:),allocatable :: property_file
logical :: pr_esp = .false.
character(len=:),allocatable :: esp_gridfile
character(len=10) :: lmoinfo_fname='xtblmoinfo'
logical :: pr_molden_input = .false.
logical :: pr_lmo = .false.
logical :: pr_local = .true.
logical :: pr_density = .false.
logical :: pr_spin_population = .false.
logical :: pr_spin_density = .false.
2 changes: 1 addition & 1 deletion src/sphereparam.f90
Original file line number Diff line number Diff line change
@@ -37,7 +37,7 @@ module xtb_sphereparam

integer :: maxwalls = 0

! old, not needed anymore
! old, only needed for QCG mode of docking
real(wp) :: boxr = -1.0_wp
real(wp) :: rabc(3) = (/-1.0_wp,-1.0_wp,-1.0_wp/)
integer :: sphere = -1
52 changes: 52 additions & 0 deletions src/type/data.f90
Original file line number Diff line number Diff line change
@@ -18,6 +18,7 @@
module xtb_type_data
use xtb_mctc_accuracy, only : wp
use xtb_type_pcem
use xtb_iff_data, only : TIFFData

implicit none

@@ -26,6 +27,35 @@ module xtb_type_data

private

type :: TIFFResults
!> Number of atoms
integer :: n
!> Ordinary numbers
integer, allocatable :: at(:)
!> Coordinates
real(wp), allocatable :: xyz(:, :)
!> Charges
real(wp), allocatable :: q(:)
!> Number of LMOs
integer :: nlmo
!> LMO positions
real(wp), allocatable :: rlmo(:, :)
!> LMO values
integer, allocatable :: lmo(:)
!> Charge related stuff
real(wp), allocatable :: qct(:, :)
!> HOMO, LUMO, Dipol
real(wp) :: elumo
real(wp) :: ehomo
real(wp) :: dipol

contains

procedure :: delete
procedure :: allocateIFFResults

end type TIFFResults

type :: scc_results
real(wp) :: e_atom = 0.0_wp
real(wp) :: e_elec = 0.0_wp
@@ -54,6 +84,7 @@ module xtb_type_data
real(wp) :: e_hb = 0.0_wp
real(wp) :: e_batm = 0.0_wp
real(wp) :: e_ext = 0.0_wp
type(TIFFResults), allocatable :: iff_results
end type scc_results

type freq_results
@@ -83,6 +114,27 @@ module xtb_type_data
end type freq_results

contains

subroutine delete(self)
class(TIFFResults), intent(out) :: self
if (allocated(self%at)) deallocate (self%at)
if (allocated(self%xyz)) deallocate (self%xyz)
if (allocated(self%q)) deallocate (self%q)
if (allocated(self%lmo)) deallocate (self%lmo)
if (allocated(self%rlmo)) deallocate (self%rlmo)
if (allocated(self%qct)) deallocate (self%qct)
end subroutine delete

subroutine allocateIFFResults(self, n)
use xtb_type_molecule, only: TMolecule
class(TIFFResults), intent(out) :: self
integer, intent(in) :: n
self%n = n
allocate(self%at(n), self%lmo(10*n), source = 0)

allocate(self%xyz(3,n), self%q(n), self%rlmo(4,10*n),&
& self%qct(n,2), source= 0.0_wp)
end subroutine allocateIFFResults

subroutine allocate_freq_results(self,n)
implicit none
2 changes: 2 additions & 0 deletions test/unit/CMakeLists.txt
Original file line number Diff line number Diff line change
@@ -28,12 +28,14 @@ set(
"coulomb"
"dftd3"
"dftd4"
"docking"
"eeq"
"gfn0"
"gfn1"
"gfn2"
"gfnff"
"hessian"
"iff"
"latticepoint"
"molecule"
"pbc-tools"
4 changes: 4 additions & 0 deletions test/unit/main.f90
Original file line number Diff line number Diff line change
@@ -24,12 +24,14 @@ program tester
use test_coulomb, only : collect_coulomb
use test_dftd3, only : collect_dftd3
use test_dftd4, only : collect_dftd4
use test_docking, only : collect_docking
use test_eeq, only : collect_eeq
use test_gfn0, only : collect_gfn0
use test_gfn1, only : collect_gfn1
use test_gfn2, only : collect_gfn2
use test_gfnff, only : collect_gfnff
use test_hessian, only : collect_hessian
use test_iff, only : collect_iff
use test_latticepoint, only : collect_latticepoint
use test_molecule, only : collect_molecule
use test_pbc_tools, only : collect_pbc_tools
@@ -54,12 +56,14 @@ program tester
new_testsuite("coulomb", collect_coulomb), &
new_testsuite("dftd3", collect_dftd3), &
new_testsuite("dftd4", collect_dftd4), &
new_testsuite("docking", collect_docking), &
new_testsuite("eeq", collect_eeq), &
new_testsuite("gfn0", collect_gfn0), &
new_testsuite("gfn1", collect_gfn1), &
new_testsuite("gfn2", collect_gfn2), &
new_testsuite("gfnff", collect_gfnff), &
new_testsuite("hessian", collect_hessian), &
new_testsuite("iff", collect_iff), &
new_testsuite("latticepoint", collect_latticepoint), &
new_testsuite("molecule", collect_molecule), &
new_testsuite("pbc-tools", collect_pbc_tools), &
2 changes: 2 additions & 0 deletions test/unit/meson.build
Original file line number Diff line number Diff line change
@@ -37,12 +37,14 @@ tests = [
'coulomb',
'dftd3',
'dftd4',
'docking',
'eeq',
'gfn0',
'gfn1',
'gfn2',
'gfnff',
'hessian',
'iff',
'latticepoint',
'molecule',
'pbc-tools',
176 changes: 176 additions & 0 deletions test/unit/test_docking.f90
Original file line number Diff line number Diff line change
@@ -0,0 +1,176 @@
! This file is part of xtb.
! SPDX-Identifier: LGPL-3.0-or-later
!
! xtb is free software: you can redistribute it and/or modify it under
! the terms of the GNU Lesser General Public License as published by
! the Free Software Foundation, either version 3 of the License, or
! (at your option) any later version.
!
! xtb is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU Lesser General Public License for more details.
!
! You should have received a copy of the GNU Lesser General Public License
! along with xtb. If not, see <https://www.gnu.org/licenses/>.

module test_docking
use testdrive, only : new_unittest, unittest_type, error_type, check_ => check, test_failed
implicit none
private

public :: collect_docking

contains

!> Collect all exported unit tests
subroutine collect_docking(testsuite)
!> Collection of tests
type(unittest_type), allocatable, intent(out) :: testsuite(:)

testsuite = [ &
new_unittest("eth_wat", test_dock_eth_wat)]!, &
! new_unittest("ellips", test_iff_ellips), &
! ]

end subroutine collect_docking


subroutine test_dock_eth_wat(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
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(9)
real(wp) :: xyz1(3,9)
integer :: n2
integer :: at2(3)
real(wp) :: xyz2(3,3)
integer,parameter :: n = 12
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

n1 = 9
at1 = [6,6,1,1,1,1,8,1,1]
xyz1 = reshape(&
&[-10.23731657240310_wp, 4.49526140160057_wp,-0.01444418438295_wp, &
& -8.92346521409264_wp, 1.93070071439448_wp,-0.06395255870166_wp, &
& -11.68223026601517_wp, 4.59275319594179_wp,-1.47773402464476_wp, &
& -8.87186386250594_wp, 5.99323351202269_wp,-0.34218865473603_wp, &
& -11.13173905821587_wp, 4.79659470137971_wp, 1.81087156124027_wp, &
& -8.00237319881732_wp, 1.64827398274591_wp,-1.90766384850460_wp, &
& -10.57950417459641_wp,-0.08678859323640_wp, 0.45375344239209_wp, &
& -7.47634978432660_wp, 1.83563324786661_wp, 1.40507877109963_wp, &
& -11.96440476139715_wp,-0.02296384600667_wp,-0.72783557254661_wp],&
& shape(xyz1))

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))

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,-11.3943358674_wp, thr=thr)
call precomp(env, iff_data, molB, molB_e, 2)

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) ,-4.5265_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 check_(error, e,0.07745330953243718_wp, thr=thr)

call env%checkpoint("xtb-IFF sp computation failed")

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")

! call check_(error, calc%topo%nbond,6)
!
! call check_(error, res_gff%e_total,-0.76480130317838_wp, thr=thr)

call molA%deallocate
call molB%deallocate

end subroutine test_dock_eth_wat

end module test_docking
101 changes: 101 additions & 0 deletions test/unit/test_iff.f90
Original file line number Diff line number Diff line change
@@ -0,0 +1,101 @@
! This file is part of xtb.
! SPDX-Identifier: LGPL-3.0-or-later
!
! xtb is free software: you can redistribute it and/or modify it under
! the terms of the GNU Lesser General Public License as published by
! the Free Software Foundation, either version 3 of the License, or
! (at your option) any later version.
!
! xtb is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU Lesser General Public License for more details.
!
! You should have received a copy of the GNU Lesser General Public License
! along with xtb. If not, see <https://www.gnu.org/licenses/>.

module test_iff
use testdrive, only : new_unittest, unittest_type, error_type, check_ => check, test_failed
implicit none
private

public :: collect_iff

contains

!> Collect all exported unit tests
subroutine collect_iff(testsuite)
!> Collection of tests
type(unittest_type), allocatable, intent(out) :: testsuite(:)

testsuite = [ &
new_unittest("iff_sp", test_iff_sp) &
]

end subroutine collect_iff


subroutine test_iff_sp(error)
use xtb_mctc_accuracy, only : wp
use xtb_mctc_systools
use xtb_type_environment
use xtb_type_options
use xtb_type_molecule
use xtb_type_data
use xtb_setparam
use xtb_setmod
use xtb_docking_param
use xtb_type_restart, only: TRestart
use xtb_iff_calculator, only : TIFFCalculator, newIFFCalculator
use xtb_iff_data, only : TIFFData
use xtb_iff_iffprepare, only : prepare_IFF
type(error_type), allocatable, intent(out) :: error
real(wp),parameter :: thr = 1.0e-10_wp
integer, parameter :: nat = 6
integer, parameter :: at(nat) = [8,1,1,8,1,1]
real(wp),parameter :: xyz(3,nat) = reshape(&
&[-0.04394313126682_wp,-0.15162589575685_wp,-0.12104386666899_wp, &
& 1.57021196791551_wp, 0.39852223676357_wp,-0.74238980424452_wp, &
& -0.52096572705612_wp,-1.61453119575215_wp,-1.08142949421996_wp, &
& -3.35527313086215_wp, 3.63311912367937_wp,-1.97688072772302_wp, &
& -4.33218397586988_wp, 4.42399249914350_wp,-0.67195282333972_wp, &
& -2.30284007732108_wp, 2.38679302539156_wp,-1.14422985849389_wp],&
& shape(xyz))

type(TMolecule) :: mol
type(TEnvironment) :: env
type(TIFFCalculator) :: calc
type(TIFFData) :: iff_data

real(wp) :: etot,egap
real(wp), allocatable :: g(:,:)
type(TRestart) :: chk
real(wp) :: sigma(3,3)
type(scc_results) :: res


logical :: exist

call init(env)
call init(mol,at,xyz)

natom_arg = '1-3'
call prepare_IFF(env, mol, iff_data)
call env%checkpoint("Could not generate electronic properties for IFF")

call newIFFCalculator(env, mol, iff_data, calc)

call env%checkpoint("xtb-IFF parameter setup failed")

allocate( g(3,mol%n), source = 0.0_wp )

call calc%singlepoint(env, mol, chk, 2, exist, etot, g, sigma, egap, res)
write(*,*) 'WW', etot

call check_(error, etot,-0.002628051951328639_wp, thr=thr)

call mol%deallocate

end subroutine test_iff_sp

end module test_iff

0 comments on commit ee5b291

Please sign in to comment.