Skip to content

Commit

Permalink
outer region saturation (#848)
Browse files Browse the repository at this point in the history 
Signed-off-by: albert <92109627+Albkat@users.noreply.github.com>
  • Loading branch information
@Albkat
Albkat authored Aug 16, 2023
1 parent 31b6daf commit 3f25dfd
Show file tree
Hide file tree
Showing 4 changed files with 125 additions and 50 deletions.
139 changes: 98 additions & 41 deletions src/oniom.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -156,7 +156,7 @@ subroutine newOniomCalculator(self, env, mol, input)

endif

! write user-defined inner region list into array !
! write user-defined inner region list as raw string into array !
self%list = TAtomList(list=input%second_arg)
call self%list%to_list(self%idx)

Expand Down Expand Up @@ -582,7 +582,7 @@ subroutine writeInfo(self, unit, mol)
end subroutine writeInfo

!> create inner region
subroutine cutbond(self, env, mol, chk, topo, inner_mol,jacobian,idx2)
subroutine cutbond(self, env, mol, chk, topo, inner_mol, jacobian, idx2)

use xtb_type_molecule, only: init
use xtb_topology, only: makeBondTopology, topologyToNeighbourList
Expand Down Expand Up @@ -611,54 +611,86 @@ subroutine cutbond(self, env, mol, chk, topo, inner_mol,jacobian,idx2)
type(TTopology), allocatable, intent(inout) :: topo

!> jacobian matrix
real(wp),allocatable, intent(inout) :: jacobian(:,:)
real(wp), allocatable, intent(inout) :: jacobian(:,:)

!> list of inner region atom indices + host atoms
integer,allocatable, intent(out) :: idx2(:)
integer, allocatable, intent(out) :: idx2(:)


!> outer region geometry
type(TMolecule), allocatable :: outer_mol

!> neighbour list
type(TNeighbourList) :: neighList

!> pair-wise indices for cutbound
integer, allocatable :: brokenBondPairs(:,:)

integer, allocatable :: at(:), at2(:)
integer, allocatable :: at(:), at_out(:)
integer, allocatable :: bonded(:, :)
real(wp), allocatable :: xyz(:, :), xyz2(:, :)
real(wp), allocatable :: xyz(:, :), xyz_out(:, :)
character(len=:),allocatable :: fname_inner

!> if the both bonded atoms are inside the inner region
logical :: inside

integer :: i, j, n, k, pre_last,iterator
integer :: i, j, k, pre_last, pre_last_out, iterator
integer :: io

!> initial number of atoms in inner and outer regions
integer :: in, out

!> Loop indices
integer :: in_itr, out_itr


!-------!
! SETUP !
!-------!

inside = .FALSE.

! initial number of atoms in inner region without LAs !
n = len(self%list)

in = len(self%list)
out = mol%n - in
in_itr = 1
out_itr = 1

! allocate accordingly the basic molecular data !
allocate (at2(size(self%idx)))
allocate (xyz2(3, size(self%idx)))
allocate (at(n))
allocate (xyz(3, n))
allocate (idx2(size(self%idx)))
allocate (at(in))
allocate (xyz(3, in))
allocate (at_out(out))
allocate (xyz_out(3, out))

! save inner region list in the matching array !
idx2=self%idx

! assign initial inner region !
do i = 1, size(self%idx)
at(i) = mol%at(self%idx(i))
at2(i) = mol%at(self%idx(i))
xyz(:, i) = mol%xyz(:, self%idx(i))
xyz2(:, i) = mol%xyz(:, self%idx(i))
end do
! divide initial mol into inner and outer regions !
do i = 1, mol%n
if (i==self%idx(in_itr)) then

if (in_itr > in) then
call env%error("The internal error, inconsistent molecular dimensionality", source=source)
return
endif

at(in_itr) = mol%at(i)
xyz(:, in_itr) = mol%xyz(:, i)
in_itr = in_itr + 1

else

if (out_itr > out) then
call env%error("The internal error, inconsistent molecular dimensionality", source=source)
return
endif

at_out(out_itr) = mol%at(i)
xyz_out(:, out_itr) = mol%xyz(:, i)
out_itr = out_itr + 1

endif
enddo

! initialiaze jacobian as identity !
call create_jacobian(jacobian,at)
Expand Down Expand Up @@ -737,40 +769,44 @@ subroutine cutbond(self, env, mol, chk, topo, inner_mol,jacobian,idx2)
end select

endif

! increment array size by 1 !
pre_last = size(at)

! adust ordinal numbers !
call resize(at)
call resize(idx2)
call resize(at_out)

! save host atom index !
idx2(pre_last+1) = bonded(2,j)
! save index of the host atom -> jacobian !
idx2(size(idx2)) = bonded(2,j)

! assign new atom as H !
at(pre_last + 1) = 1

at(size(at)) = 1
at_out(size(at_out)) = 1

! adjust coordinate matrix !
call resize(xyz)
call resize(xyz_out)

! increment Jacobian !
call resize_jacobian(jacobian)

! determine new position of added H atom !
call coord(env,mol,xyz,self%idx(i),bonded(2,j),jacobian,i)
call coord(env,mol,xyz,xyz_out,self%idx(i),bonded(2,j),jacobian,i)

end if

inside = .FALSE.


! if atom in the list is bonded !
else if (bonded(2, j) == self%idx(i)) then

! iterate again through the list !
do k = 1, size(self%idx)

! inside inner region !
if (self%idx(k) == bonded(1, j)) then
! inside inner region !
inside = .TRUE.
end if

end do

! bond is broken !
Expand All @@ -788,26 +824,29 @@ subroutine cutbond(self, env, mol, chk, topo, inner_mol,jacobian,idx2)

endif

! increment array size by 1 !
pre_last = size(at)
call resize(at)
call resize(at_out)
call resize(idx2)

! save host atom index !
idx2(pre_last+1) = bonded(1,j)
! save index of the host atom -> jacobian !
idx2(size(idx2)) = bonded(1,j)

! assign new atom as H !
at(pre_last + 1) = 1
at(size(at)) = 1
at_out(size(at_out)) = 1

! adjust coordinate matrix !
call resize(xyz)
call resize(xyz_out)

! increment Jacobian matrix !
call resize_jacobian(jacobian)

! determine new position of added H atom !
call coord(env,mol, xyz, self%idx(i), bonded(1, j),jacobian,i)
call coord(env,mol, xyz, xyz_out, self%idx(i), bonded(1, j),jacobian,i)

end if

inside = .FALSE.

end if
Expand All @@ -817,8 +856,21 @@ subroutine cutbond(self, env, mol, chk, topo, inner_mol,jacobian,idx2)
!----------------!
! POSTPROCESSING !
!----------------!

! initialize new mol object !

! outer region saturation !
if (set%oniom_settings%outer) then

! initialize !
allocate(outer_mol)
call init(outer_mol, at_out, xyz_out)

! create Xmol file !
call open_file(io, "outer_region.xyz", "w")
call writeMolecule(outer_mol, io, filetype%xyz)
call close_file(io)
end if

! initialize inner region mol !
call init(inner_mol, at, xyz)

! create Xmol file for inner region !
Expand All @@ -832,6 +884,7 @@ subroutine cutbond(self, env, mol, chk, topo, inner_mol,jacobian,idx2)
call writeMolecule(inner_mol, io, filetype%xyz)
call close_file(io)


end subroutine cutbond

!> increase atomic number array by 1
Expand Down Expand Up @@ -908,7 +961,7 @@ subroutine resize_jacobian(matrix)
end subroutine resize_jacobian

!> calculate new postion for LA and corresponding J
subroutine newcoord(env,mol,xyz,idx1,idx2,jacobian,connectorPosition)
subroutine newcoord(env,mol,xyz,xyz_out,idx1,idx2,jacobian,connectorPosition)

implicit none

Expand All @@ -921,9 +974,12 @@ subroutine newcoord(env,mol,xyz,idx1,idx2,jacobian,connectorPosition)
!> molecular structure data
type(TMolecule), intent(in) :: mol

!> coordinate matrix
!> inreg coordinate matrix
real(wp), intent(inout) :: xyz(:, :)

!> outreg coordinate matrix
real(wp), intent(inout) :: xyz_out(:, :)

!> connector (one that stays in the inner region)
integer, intent(in) :: idx1

Expand Down Expand Up @@ -1194,6 +1250,7 @@ subroutine newcoord(env,mol,xyz,idx1,idx2,jacobian,connectorPosition)

! LA coordinates !
xyz(:, size(xyz, 2)) = xyz1 + (xyz2 - xyz1) * prefactor
xyz_out(:,size(xyz_out,2)) = xyz2 + (xyz1 - xyz2) * prefactor

! determine the difference between LA and LAH cooordinates !
! (yes -> change from derived to fixed) !
Expand Down
11 changes: 10 additions & 1 deletion src/set_module.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -781,7 +781,7 @@ subroutine rdcontrol(fname,env,copy_file)
case('write' ); call rdblock(env,set_write, line,id,copy,err,ncount)
case('gfn' ); call rdblock(env,set_gfn, line,id,copy,err,ncount)
case('scc' ); call rdblock(env,set_scc, line,id,copy,err,ncount)
case('oniom' ); call rdblock(env,set_oniom, line,id,copy,err,ncount)
case('oniom' ); call rdblock(env,set_oniom, line,id,copy,err,ncount)
case('opt' ); call rdblock(env,set_opt, line,id,copy,err,ncount)
case('hess' ); call rdblock(env,set_hess, line,id,copy,err,ncount)
case('md' ); call rdblock(env,set_md, line,id,copy,err,ncount)
Expand Down Expand Up @@ -1450,13 +1450,16 @@ end subroutine set_gfn

!> set ONIOM functionality
subroutine set_oniom(env,key,val)

implicit none

!> pointer to the error routine
character(len=*), parameter :: source = 'set_oniom'

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

!> parsed key
character(len=*), intent(in) :: key

!> key=val
Expand All @@ -1467,6 +1470,7 @@ subroutine set_oniom(env,key,val)
logical, save :: set2 = .true.
logical, save :: set3 = .true.
logical, save :: set4 = .true.
logical, save :: set5 = .true.

select case(key)
case default
Expand All @@ -1486,6 +1490,11 @@ subroutine set_oniom(env,key,val)
case('ignore topo')
if (getValue(env,val,ldum).and.set4) set%oniom_settings%ignore_topo = .true.
set4=.false.

case('outer')
if (getValue(env,val,ldum).and.set5) set%oniom_settings%outer = .true.
set5 = .false.

end select

end subroutine set_oniom
Expand Down
15 changes: 12 additions & 3 deletions src/setparam.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -102,6 +102,7 @@ module xtb_setparam
integer,parameter :: p_pcem_orca = 2

type oniom_settings

!> inner region charge
integer :: innerchrg

Expand All @@ -122,20 +123,28 @@ module xtb_setparam

!> print optimization logs for inner region calculations
logical :: logs = .false.

!> if saturate outer region
logical :: outer = .false.

!> log units
integer:: ilog1, ilog2
end type oniom_settings

end type oniom_settings

type qm_external

character(len=:),allocatable :: path
character(len=:),allocatable :: executable
character(len=:),allocatable :: input_file
character(len=:),allocatable :: input_string

!> if input_file exist
logical :: exist
!! if input_file exist

!> special case of the oniom embedding
logical :: oniom=.false.
!! special case of the oniom embedding

end type qm_external

integer, parameter :: p_ext_vtb = -1
Expand Down
10 changes: 5 additions & 5 deletions test/unit/test_oniom.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -179,8 +179,8 @@ subroutine test_oniom_cutbond(error)
type(TOniomCalculator) :: calc
type(oniom_input) :: oniom_input3
type(scc_results) :: results
real(wp),allocatable :: jacobian(:,:)
integer,allocatable :: idx_orig(:)
real(wp), allocatable :: jacobian(:,:)
integer, allocatable :: idx2(:)
real(wp) :: energy, hlgap, sigma(3,3)
real(wp), allocatable :: gradient(:,:)

Expand Down Expand Up @@ -220,13 +220,13 @@ subroutine test_oniom_cutbond(error)
call check_(error,norm2(gradient), 0.000327194200_wp, thr=thr)

!> cutbond
call calc%cutbond(env,mol,chk,calc%topo,inner_mol,jacobian,idx_orig)
call calc%cutbond(env,mol,chk,calc%topo,inner_mol,jacobian,idx2)

call check_(error,idx_orig(5),2)
call check_(error,idx2(5),2)
call check_(error,inner_mol%xyz(1,5),-4.38055107024201_wp,thr=thr)
call check_(error,inner_mol%xyz(2,5),4.42017847630110_wp,thr=thr)
call check_(error,inner_mol%xyz(3,5),0.00000120013337_wp,thr=thr)

end subroutine test_oniom_cutbond

!---------------------------------------------
Expand Down

0 comments on commit 3f25dfd

Please sign in to comment.