diff --git a/man/xcontrol.7.adoc b/man/xcontrol.7.adoc
index efc08e9eb..ba4efef5a 100644
--- a/man/xcontrol.7.adoc
+++ b/man/xcontrol.7.adoc
@@ -326,6 +326,19 @@ $modef
*mode*='int'::
can set by --modef via cmdline
+$oniom
+~~~~~~
+*inner logs*='bool'::
+ to print optimization log files for model region geometry (`high.inner_region.log` and `low.inner_region.log`)
+
+*derived k*='bool'::
+ to calculate prefactor *k* and create jacobian dynamically (see more )
+
+*silent*='bool'::
+ to hide the execution runs of external software
+
+
+
$opt
~~~~
*engine*='method'::
diff --git a/man/xtb.1.adoc b/man/xtb.1.adoc
index c867bed5f..5b32bb116 100644
--- a/man/xtb.1.adoc
+++ b/man/xtb.1.adoc
@@ -24,8 +24,8 @@ geometries, frequencies and non-covalent interactions (GFN)
as well as with the ionisation potential and
electron affinity (IPEA) parametrisation of the GFN1 Hamiltonian.
The generalized born (GB) model with solvent accessable surface area (SASA)
-is also available available in this version.
-Ground state calculations for the simplified Tamm-Danceoff approximation (sTDA)
+is also available in this version.
+Ground state calculations for the simplified Tamm-Dancoff approximation (sTDA)
with the vTB model are currently not implemented.
GEOMETRY INPUT
@@ -66,6 +66,9 @@ OPTIONS
*-c, --chrg* 'INT'::
specify molecular charge as 'INT', overrides `.CHRG` file and `xcontrol` option
+*-c, --chrg* 'INT:INT'::
+ specify charges for 'inner region:outer region' for `oniom` calculation, overrides `.CHRG` file and `xcontrol` option
+
*-u, --uhf* 'INT'::
specify number of unpaired electrons as 'INT', overrides `.UHF` file and `xcontrol` option
@@ -79,10 +82,10 @@ OPTIONS
use tblite library as implementation of for xTB
*--oniom* 'METHOD' 'LIST'::
- use subtractive embedding via ONIOM method. 'METHOD' is given as `inner:outer`
- where `inner` can be 'orca', 'turbomole', 'gfn2', 'gfn1', or 'gfnff' and
- `outer` can be 'gfn2', 'gfn1', or 'gfnff'.
- The inner region is given as a comma separated indices directly in the commandline
+ use subtractive embedding via ONIOM method. 'METHOD' is given as `high:low`
+ where `high` can be 'orca', 'turbomole', 'gfn2', 'gfn1', or 'gfnff' and
+ `low` can be 'gfn2', 'gfn1', or 'gfnff'.
+ The inner region is given as comma-separated indices directly in the commandline
or in a file with each index on a separate line.
*--etemp* 'REAL'::
@@ -263,6 +266,9 @@ GENERAL
*-h, --help*::
show help page
+*--cut*::
+ create inner region for `oniom` calculation without performing any calcultion
+
ENVIRONMENT VARIABLES
---------------------
`xtb(1)` accesses a path-like variable to determine the location of its
@@ -421,7 +427,7 @@ After `xtb(1)` has evaluated the all input sources it immediately enters
the production mode. Severe errors will lead to an abnormal termination
which is signalled by the printout to STDERR and a non-zero return value
(usually 128). All non-fatal errors are summerized in the end of the calculation
-in one block, right bevor the timing analysis.
+in one block, right before the timing analysis.
To aid the user to fix the problems generating these warnings a brief
summary of each warning with its respective string representation in the
diff --git a/src/bias_path.f90 b/src/bias_path.f90
index ea27fc3f4..9db3821a7 100644
--- a/src/bias_path.f90
+++ b/src/bias_path.f90
@@ -282,7 +282,7 @@ subroutine bias_path(env, mol, chk, calc, egap, et, maxiter, epot, grd, sigma)
enddo
factor2= factor2 + alp_change
enddo bias_loop
-
+
!! ------------------------------------------------------------------------
! output and find path yielding product
!! ------------------------------------------------------------------------
diff --git a/src/extern/driver.f90 b/src/extern/driver.f90
index 914a18326..471aa2bbd 100644
--- a/src/extern/driver.f90
+++ b/src/extern/driver.f90
@@ -117,7 +117,7 @@ subroutine singlepoint(self, env, mol, chk, printlevel, restart, &
!$omp critical (turbo_lock)
inquire (file='gradient', exist=exist)
if (exist) then
- call rdtm(mol%n, .true., energy, gradient, xyz_cached)
+ call rdtm(env,mol%n, .true., energy, gradient, xyz_cached)
cache = all(abs(xyz_cached - mol%xyz) < 1.e-10_wp)
end if
if (.not. cache) then
@@ -135,7 +135,7 @@ subroutine singlepoint(self, env, mol, chk, printlevel, restart, &
end if
write (env%unit, '(72("="))')
- call rdtm(mol%n, .true., energy, gradient, xyz_cached)
+ call rdtm(env,mol%n, .true., energy, gradient, xyz_cached)
end if
!$omp end critical (turbo_lock)
diff --git a/src/extern/orca.f90 b/src/extern/orca.f90
index 4238abd21..5cb2b5d32 100644
--- a/src/extern/orca.f90
+++ b/src/extern/orca.f90
@@ -57,27 +57,38 @@ module xtb_extern_orca
contains
-
+!-----------------------------------------------------
+! Check the ORCA execution settings
+!-----------------------------------------------------
subroutine checkOrca(env, ext)
character(len=*), parameter :: source = 'extern_orca_checkOrca'
+ !! the name of the error producer routine
+
type(TEnvironment), intent(inout) :: env
+ !! Calculation environment to handle I/O stream and error log
type(qm_external), intent(inout) :: ext
+ !! Settings for the Orca calculation
+
character(len=:),allocatable :: homedir,syspath
character(len=:),allocatable :: line
character(len=5) :: chdum
logical :: exist
logical :: chk_xyzfile
logical :: chk_engrad
- integer :: iorca ! file handle
+ integer :: iorca
+ !! file handle
integer :: err
integer :: idx_bang,idx_star,idx_engr,idx_file,idx_xmol,idx_spce
!$ integer,external :: omp_get_num_threads
- ! check for ORCA executable
- if (allocated(ext%executable)) then ! user input
+ !> check if the ORCA executable exist
+ if (allocated(ext%executable)) then
+ !! user input
+
if (ext%executable(1:1).eq.'~') then
- ! this is relative to the users home, expand it
- call rdvar('HOME',homedir)
+ !! this is relative to the users home, expand it
+ call rdvar('HOME',homedir)
+ !! read environment variable
ext%executable = homedir // ext%executable(2:)
endif
inquire(file=ext%executable,exist=exist)
@@ -85,15 +96,19 @@ subroutine checkOrca(env, ext)
call env%error("'"//ext%executable//"' was not found, please check",source)
return
endif
- else ! no executable provided, lets find it
+
+ else
+ !! no executable provided, lets find it
call rdvar('PATH',syspath)
call rdpath(syspath,'orca',ext%executable,exist)
+ !! find ORCA excutable in PATH variable
if (.not.exist) then
call env%error('Could not locate orca executable',source)
return
endif
endif
-
+
+ !> check for wrong executable with the same name
if (index(ext%executable,'/usr') == 1) then
inquire(file=ext%executable//'_scf',exist=exist)
if (.not. exist) then
@@ -103,25 +118,31 @@ subroutine checkOrca(env, ext)
endif
endif
- ! see if there is a preference for an input file
+ !> see if there is a preference for an input file
if (allocated(ext%input_file)) then
inquire(file=ext%input_file,exist=exist)
ext%exist = exist
else
ext%input_file = "orca-"//get_random_name()//'.inp'
- inquire(file=ext%input_file,exist=exist)
+ !! Create random file for the ORCA input file
+ inquire(file=ext%input_file,exist=exist)
endif
- ! sanity check
+
+ !> sanity check
if (ext%exist) then
+
call open_file(iorca,ext%input_file,'r')
+ !> if exists, but empty
if (iorca.eq.-1) then
call env%error("ORCA input file '"//ext%input_file//"' just vanished!",source)
return
endif
+
chk_engrad = .false.
chk_xyzfile = .false.
do
call strip_line(iorca,line,err)
+ !! to read line from iorca unit
if (err.ne.0) exit
idx_bang = index(line,'!')
idx_star = index(line,'*')
@@ -131,10 +152,12 @@ subroutine checkOrca(env, ext)
idx_spce = index(trim(line),' ',back=.true.)
if (idx_bang.gt.0 .and. idx_engr.gt.0 .and. idx_engr.gt.idx_bang) &
chk_engrad = .true.
+
if (idx_star.gt.0 .and. idx_file.gt.0 .and. idx_xmol.gt.0 .and. &
&idx_file.gt.idx_star .and. idx_xmol.gt.idx_file) then
chk_xyzfile = .true.
ext%input_string = trim(adjustl(line(idx_spce:)))
+ !! calculation method
endif
enddo
call close_file(iorca)
@@ -145,7 +168,7 @@ subroutine checkOrca(env, ext)
endif
else
- ! check for the input line
+ !! if not exist, check for the input line
if (allocated(ext%input_string)) then
if (index(ext%input_string,'!') == 1) then
ext%input_string = ext%input_string(2:)
@@ -158,28 +181,39 @@ subroutine checkOrca(env, ext)
end subroutine checkOrca
-
-!> Create a new calculator for driving the Orca program
-subroutine newOrcaCalculator(self, env, ext)
- !> Instance of the Orca calculator
+!-----------------------------------------------------
+! Create a new calculator for driving the Orca program
+!-----------------------------------------------------
+subroutine newOrcaCalculator(self, env, ext,oniom)
+
+ implicit none
type(TOrcaCalculator), intent(out) :: self
- !> Calculation environment
+ !! Instance of the Orca calculator
type(TEnvironment), intent(inout) :: env
- !> Settings for the Orca calculator
+ !! Calculation environment
type(qm_external), intent(in) :: ext
+ !! Settings for the Orca calculator
+ logical, intent(in),optional :: oniom
self%ext = ext
+ !! to save external settings from TSet%ext_orca to the new caluclator
self%threadsafe = .false.
+ !! not to call orca in parallel
+ self%ext%oniom=oniom
+ !! if oniom calc
call checkOrca(env, self%ext)
-end subroutine newOrcaCalculator
-subroutine writeOrcaInp(io, mol, input, mode)
+end subroutine newOrcaCalculator
+
+!-----------------------------------------------------
+! Create *.inp for the ORCA run
+!-----------------------------------------------------
+subroutine writeOrcaInp(io,mol,input,mode)
integer, intent(in) :: io
type(TMolecule), intent(in) :: mol
character(len=*), intent(in) :: input
character(len=*), intent(in) :: mode
-
integer :: i, nel, mult
write(io,'("#",1x,a)') &
@@ -187,7 +221,7 @@ subroutine writeOrcaInp(io, mol, input, mode)
!$omp parallel
!$omp master
!$ write(io,'("%",a,/,3x,a,1x,i0,/,a)') &
- !$ "pal","nprocs",omp_get_num_threads(),"end"
+ !$ "pal","nprocs",omp_get_num_threads()-1,"end"
!$omp end master
!$omp end parallel
write(io,'("%",a,1x,i0)') &
@@ -215,6 +249,7 @@ subroutine writeOrcaInp(io, mol, input, mode)
write(io,'(3x,a2,3(2x,F20.14))') toSymbol(mol%at(i)), mol%xyz(:,i)*autoaa
end do
write(io,'("*",/)')
+
end subroutine writeOrcaInp
@@ -274,7 +309,7 @@ subroutine readOrcaHess(io, hess, dipgrad)
end if
if (line == "$dipole_derivatives") then
- call getline(io, line, stat)
+ call getline(io, line, stat)
read(line, *) ndim
do ii = 1, ndim
read(io, *) dipgrad(:, ii)
@@ -284,38 +319,56 @@ subroutine readOrcaHess(io, hess, dipgrad)
end do
end subroutine readOrcaHess
-
+!-----------------------------------------------------
+! To create or modify input file and execute the ORCA
+!-----------------------------------------------------
subroutine runOrca(env,ext,mol,energy,gradient)
+
character(len=*), parameter :: source = 'extern_orca_runOrca'
+ !! the name of the error producer routine
type(TEnvironment), intent(inout) :: env
+ !! Calculation environment to handle I/O stream and error log
type(qm_external), intent(in) :: ext
+ !! Settings for the Orca calculation
type(TMolecule), intent(in) :: mol
+ !! Molecular structure data
real(wp),intent(out) :: energy
real(wp),intent(out) :: gradient(:, :)
integer :: i,j,err
- integer :: iorca ! file handle
+ integer :: iorca
+ !! file IO unit
logical :: exist
character(len=:),allocatable :: line
character(len=:),allocatable :: outfile
+ character(len=:),allocatable :: tmpfile
!$omp critical (orca_lock)
+
+ !> To decide whether to create or modify the input_file
if (ext%exist) then
- ! we dump the name of the external xyz file to input_string... not cool
+ !! we dump the name of the external xyz file to input_string... not cool
call open_file(iorca,ext%input_string,'w')
call writeMolecule(mol, iorca, format=fileType%xyz)
call close_file(iorca)
else
+ !! create new.inp file
call open_file(iorca,ext%input_file,'w')
call writeOrcaInp(iorca,mol,ext%input_string, "engrad")
call close_file(iorca)
endif
-
+ !> Actual orca cml run
write(env%unit,'(72("="))')
write(env%unit,'(1x,"*",1x,a)') &
"letting orca take over the control..."
- call execute_command_line('exec 2>&1 '//ext%executable//' '// &
+ if (set%oniom_settings%silent) then
+ call execute_command_line('exec 2>&1 '//ext%executable//' '// &
+ ext%input_file//'>orca.out',exitstat=err)
+ else
+ call execute_command_line('exec 2>&1 '//ext%executable//' '// &
ext%input_file,exitstat=err)
+ endif
+
if (err.ne.0) then
call env%error('orca returned with non-zero exit status, doing the same',source)
else
@@ -323,12 +376,16 @@ subroutine runOrca(env,ext,mol,energy,gradient)
"successful orca run, taking over control again..."
endif
write(env%unit,'(72("="))')
-
+
+ !> find output .engrad file
i = index(ext%input_file,'.inp')
+
if (i > 0) then
outfile = ext%input_file(:i-1)//'.engrad'
+ tmpfile = ext%input_file(:i-1)
else
outfile = ext%input_file//'.engrad'
+ tmpfile = ext%input_file
endif
inquire(file=outfile,exist=exist)
if (.not.exist) then
@@ -336,8 +393,18 @@ subroutine runOrca(env,ext,mol,energy,gradient)
else
call open_file(iorca,outfile,'r')
call readOrcaEngrad(iorca, energy, gradient)
- call close_file(iorca)
+ if (set%ceasefiles) then
+ call remove_file(iorca)
+ call env%io%deleteFile(tmpfile//".gbw")
+ call env%io%deleteFile(tmpfile//".densities")
+ call env%io%deleteFile(tmpfile//".rr")
+ call env%io%deleteFile(tmpfile//"_property.txt")
+ call env%io%deleteFile("orca.out")
+ else
+ call close_file(iorca)
+ endif
endif
+
!$omp end critical (orca_lock)
end subroutine runOrca
diff --git a/src/extern/turbomole.f90 b/src/extern/turbomole.f90
index b3966beaa..271250ef9 100644
--- a/src/extern/turbomole.f90
+++ b/src/extern/turbomole.f90
@@ -47,6 +47,7 @@ module xtb_extern_turbomole
public :: TTMCalculator, newTMCalculator
type, extends(TCalculator) :: TTMCalculator
+ type(qm_external) :: ext
integer :: extcode
integer :: extmode
contains
@@ -60,17 +61,27 @@ module xtb_extern_turbomole
contains
-!> Create a new calculator for driving the Turbomole program
-subroutine newTMCalculator(self, extcode, extmode)
- !> Instance of the Turbomole calculator
+!----------------------------------------------------------
+! Create a new calculator for driving the Turbomole program
+!----------------------------------------------------------
+subroutine newTMCalculator(self,extcode,extmode)
+
+ implicit none
type(TTMCalculator), intent(out) :: self
- integer, intent(in) :: extcode, extmode
+ !! Instance of the Turbomole calculator
+ integer, intent(in) :: extcode
+ !! RI, RI+d3+gCP,NoRI
+ integer, intent(in) :: extmode
+ !! ridft+rdgrad, ridft+rdgrad+dftd3+gcp
self%threadsafe = .false.
+ !! No parallelization
self%extcode = extcode
self%extmode = extmode
+
end subroutine newTMCalculator
+
subroutine singlepoint(self, env, mol, chk, printlevel, restart, &
& energy, gradient, sigma, hlgap, results)
!> Source of the generated errors
@@ -118,7 +129,7 @@ subroutine singlepoint(self, env, mol, chk, printlevel, restart, &
dipole(:) = 0.0_wp
call external_turbomole(env,mol%n,mol%at,mol%xyz,chk%wfn%nel,chk%wfn%nopen, &
- & self%extcode,self%extmode,.true.,energy,gradient,results%dipole,self%lSolv)
+ & self%extcode,self%extmode,.true.,energy,gradient,results%dipole,self%lSolv,mol%chrg,mol%uhf)
call env%check(exitRun)
if (exitRun) then
@@ -163,58 +174,128 @@ subroutine singlepoint(self, env, mol, chk, printlevel, restart, &
endif
end subroutine singlepoint
-subroutine external_turbomole(env,n,at,xyz,nel,nopen,extcode,extmode,grd,eel,g,dip,lsolv)
+!---------------------------------------------------------
+! To run TURBOMOLE, and cefine if needed
+!---------------------------------------------------------
+subroutine external_turbomole(env,n,at,xyz,nel,nopen,extcode,extmode,grd,eel,g,dip,lsolv,chrg,uhf)
+
use xtb_mctc_accuracy, only : wp
use xtb_setparam
+
implicit none
- !> Computational environment
+ character(len=*),parameter :: source = 'external_turbomole'
type(TEnvironment), intent(inout) :: env
- integer n, at(n), nel, nopen
- logical grd,lsolv
+ !! Calculation environment to handle I/O stream and error log
+ integer :: n, at(n), nel, nopen
+ !! structural information from mol(Tmolecule) and chk(TRestart)
+ logical :: grd
+ !! if gradient calcultion is needed
+ logical :: lsolv
+ !! if solvated
integer, intent(in) :: extcode, extmode
- real(wp) xyz(3,n)
- real(wp) xyz_cached(3,n)
- real(wp) g (3,n)
- real(wp) eel
- real(wp) dip(3)
+ !! turbomole settings
+ real(wp) :: xyz(3,n)
+ !! coordinates
+ real(wp) :: xyz_cached(3,n)
+ !! to save xyz values
+ real(wp) :: g (3,n)
+ !! gradients
+ real(wp) :: eel
+ !! electronic energy
+ real(wp) :: dip(3)
+ !! dipole moments
+ real(wp), intent(in) :: chrg
+ !! charge
+ integer, intent(in) :: uhf
+ !! multiplicity, number of unpaired electrons
+
character(len=255) atmp
character(len=:), allocatable :: syspath, cefine
- logical :: cache, exist
+ logical :: cache
+ !! to check if molecule moves b/n optimization runs
+ logical :: exist
+ integer :: err
cache = .false.
dip=0
-
+
+ !> TM input
inquire(file="control", exist=exist)
+ !! if control file present in the calc dir
if (.not.exist) then
- call rdvar("PATH", syspath)
- call rdpath(syspath, "cefine", cefine, exist)
- if (exist) then
- call wrtm(n,at,xyz)
- call execute_command_line("exec "//cefine//" --func tpss --def2/SVP --cosmo 2.38 --d4 -sym c1 --noopt")
- end if
+
+ !> cefine
+ if(len(set%ext_turbo%input_string).ne.0) then
+ call rdvar("PATH", syspath)
+ call rdpath(syspath, "cefine", cefine, exist)
+ if (exist) then
+ call wrtm(n,at,xyz)
+ if (allocated(set%ext_turbo%input_string)) then
+ call execute_command_line("exec "//cefine//" "//set%ext_turbo%input_string)
+ else
+ call execute_command_line("exec "//cefine//" --func b97-3c")
+ endif
+ else
+ call env%error("No cefine binary is found", source)
+ return
+ end if
+
+ else
+ call writeControl(chrg,uhf)
+ !! default control file
+ endif
+
end if
+ inquire(file="control", exist=exist)
if (.not.exist) then
call env%error("No 'control' file in current directory")
return
end if
-
+
+
! TM (RI)
if(extcode.eq.1)then
!$omp critical (turbo_lock)
inquire(file='gradient', exist=exist)
if (exist .and. grd) then
- call rdtm(n,grd,eel,g,xyz_cached)
- cache = all(abs(xyz_cached - xyz) < 1.e-10_wp)
+ call rdtm(env,n,grd,eel,g,xyz_cached)
+ cache = all(abs(xyz_cached - xyz) < 1.e-7_wp)
end if
if (.not.cache) then
call wrtm(n,at,xyz)
if(extmode.eq.1)then
- call execute_command_line('exec ridft > job.last 2>> /dev/null')
- if(grd)call execute_command_line('exec rdgrad >> job.last 2>> /dev/null ')
- endif
+
+ write(env%unit,'(72("="))')
+ write(env%unit,'(1x,"*",1x,a)') &
+ "letting TURBOMOLE take over the control..."
+
+ if (set%oniom_settings%silent) then
+ call execute_command_line('exec ridft > job.last 2>> /dev/null',exitstat=err)
+ if(grd)call execute_command_line('exec rdgrad >> job.last 2>> /dev/null ',exitstat=err)
+ else
+ call execute_command_line('exec ridft | tee job.last 2>> /dev/null',exitstat=err)
+ if(grd)call execute_command_line('exec rdgrad |tee -a job.last 2>> /dev/null ',exitstat=err)
+ endif
+
+ if (err.ne.0) then
+ call env%error('TURBOMOLE returned with non-zero exit status, doing the same',source)
+ else
+ write(env%unit,'(1x,"*",1x,a)') &
+ "successful TURBOMOLE run, taking over control again..."
+ endif
+ write(env%unit,'(72("="))')
+ endif
call extcodeok(extcode)
- call rdtm(n,grd,eel,g,xyz_cached)
+ call rdtm(env,n,grd,eel,g,xyz_cached)
+ if (set%ceasefiles) then
+ call env%io%deleteFile('job.last')
+ call env%io%deleteFile('mos')
+ call env%io%deleteFile('statistics')
+ call env%io%deleteFile('basis')
+ call env%io%deleteFile('auxbasis')
+ call env%io%deleteFile('coord')
+ endif
end if
!$omp end critical (turbo_lock)
return
@@ -225,8 +306,8 @@ subroutine external_turbomole(env,n,at,xyz,nel,nopen,extcode,extmode,grd,eel,g,d
!$omp critical (turbo_lock)
inquire(file='gradient', exist=exist)
if (exist .and. grd) then
- call rdtm(n,grd,eel,g,xyz_cached)
- cache = all(abs(xyz_cached - xyz) < 1.e-10_wp)
+ call rdtm(env,n,grd,eel,g,xyz_cached)
+ cache = all(abs(xyz_cached - xyz) < 1.e-7_wp)
end if
if (.not.cache) then
call wrtm(n,at,xyz)
@@ -237,7 +318,16 @@ subroutine external_turbomole(env,n,at,xyz,nel,nopen,extcode,extmode,grd,eel,g,d
call execute_command_line('exec gcp coord -file -grad >>job.last 2>>/dev/null')
endif
call extcodeok(extcode)
- call rdtm(n,.true.,eel,g,xyz_cached)
+ call rdtm(env,n,.true.,eel,g,xyz_cached)
+
+ if (set%ceasefiles) then
+ call env%io%deleteFile('job.last')
+ call env%io%deleteFile('mos')
+ call env%io%deleteFile('statistics')
+ call env%io%deleteFile('basis')
+ call env%io%deleteFile('auxbasis')
+ call env%io%deleteFile('coord')
+ endif
end if
!$omp end critical (turbo_lock)
return
@@ -248,8 +338,8 @@ subroutine external_turbomole(env,n,at,xyz,nel,nopen,extcode,extmode,grd,eel,g,d
!$omp critical (turbo_lock)
inquire(file='gradient', exist=exist)
if (exist .and. grd) then
- call rdtm(n,grd,eel,g,xyz_cached)
- cache = all(abs(xyz_cached - xyz) < 1.e-10_wp)
+ call rdtm(env,n,grd,eel,g,xyz_cached)
+ cache = all(abs(xyz_cached - xyz) < 1.e-7_wp)
end if
if (.not.cache) then
call wrtm(n,at,xyz)
@@ -258,7 +348,15 @@ subroutine external_turbomole(env,n,at,xyz,nel,nopen,extcode,extmode,grd,eel,g,d
if(grd)call execute_command_line('exec grad >> job.last 2>> /dev/null')
endif
call extcodeok(extcode)
- call rdtm(n,grd,eel,g,xyz_cached)
+ call rdtm(env,n,grd,eel,g,xyz_cached)
+ if (set%ceasefiles) then
+ call env%io%deleteFile('job.last')
+ call env%io%deleteFile('mos')
+ call env%io%deleteFile('statistics')
+ call env%io%deleteFile('basis')
+ call env%io%deleteFile('auxbasis')
+ call env%io%deleteFile('coord')
+ endif
end if
!$omp end critical (turbo_lock)
return
@@ -269,9 +367,41 @@ subroutine external_turbomole(env,n,at,xyz,nel,nopen,extcode,extmode,grd,eel,g,d
end subroutine external_turbomole
-!ccccccccccccccccccccccccccccccccc
-! TM
-!ccccccccccccccccccccccccccccccccc
+!-------------------------------------------------------
+! create default TM control -> future TMprep
+!-------------------------------------------------------
+subroutine writeControl(chrg,uhf)
+
+ use xtb_mctc_accuracy, only : wp
+ use xtb_mctc_symbols, only : toSymbol
+ implicit none
+ real(wp), intent(in) :: chrg
+ !! charge
+ integer, intent(in) :: uhf
+ !! number of unpaired electrons
+ integer :: iunit
+ !! number of electrins
+
+ iunit=39
+ open(newunit=iunit,file='control')
+
+ write(iunit,'(a)')'$coord file=coord'
+ write(iunit,'(a,i0,a,i0)')'$eht charge=',nint(chrg),' unpaired=',uhf
+ write(iunit,'(a)')'$symmetry c1'
+ write(iunit,'(a)')'$atoms'
+ write(iunit,'(3x,a)')'basis =def2-mTZVP'
+ write(iunit,'(a)')'$dft'
+ write(iunit,'(3x,a)')'functional b97-3c'
+ write(iunit,'(3x,a)')'gridsize m4'
+ write(iunit,'(a)')'$rij'
+ write(iunit,'(a)')'$energy file=energy'
+ write(iunit,'(a)')'$grad file=gradient'
+ write(iunit,'(a)')'$disp3 -bj -abc'
+ write(iunit,'(a)')'$end'
+
+ close(iunit)
+
+end subroutine writeControl
subroutine wrtm(n,at,xyz)
use xtb_mctc_accuracy, only : wp
@@ -294,9 +424,46 @@ subroutine wrtm(n,at,xyz)
end subroutine wrtm
-subroutine rdtm(n,grd,e,g,xyz)
+!----------------------------------------
+! read TM energy and gradient files
+!----------------------------------------
+!subroutine readTM(n,ifgrd,energy,gradient,xyz)
+
+ ! use xtb_mctc_accuracy, only : wp
+ !use xtb_filetools, only : open_file, close_file, remove_file
+ ! use xtb_readin , only : strip_line,getValue
+
+! implicit none
+! integer, intent(in) :: n
+ !! number of atoms
+! logical, intent(in) :: ifgrd
+ !! if gradient fileshould be read
+! real(wp), intent(out) :: energy, gradient(3,n), xyz(3,n)
+ !! TM output
+
+! integer :: grd, e
+ !! file units for 'gradient' and 'energy' files
+! character(len=:), allocatable :: line
+ !! tmp line
+! logical :: exist
+
+
+! if (.not.ifgrd) then
+! call open_file(e,'energy','r')
+! do
+! call strip_line(e,line,exist)
+! if (exist)
+! call readline()
+! enddo
+! endif
+!
+!
+!end subroutine readTM
+subroutine rdtm(env,n,grd,e,g,xyz)
+
use xtb_mctc_accuracy, only : wp
implicit none
+ type(TEnvironment), intent(inout) :: env
integer n, iunit, i, nl, j, nn
logical grd
real(wp) g(3,n), e, xx(10), x, y, z, xyz(3,n)
@@ -310,8 +477,14 @@ subroutine rdtm(n,grd,e,g,xyz)
101 read(iunit,'(a)',end=102)a1
call readl(a1,xx,nn)
if(nn.ge.4) e=xx(2)
+ !! to assign the last entry
goto 101
- 102 close(iunit)
+ 102 continue
+ if (set%ceasefiles) then
+ call env%io%deleteFile('energy')
+ else
+ close(iunit)
+ endif
return
endif
@@ -344,8 +517,12 @@ subroutine rdtm(n,grd,e,g,xyz)
do i=1,n
read(iunit,*)g(1,i),g(2,i),g(3,i)
enddo
-
- close(iunit)
+ if (set%ceasefiles) then
+ call env%io%deleteFile('energy')
+ call env%io%deleteFile('gradient')
+ else
+ close(iunit)
+ endif
end subroutine rdtm
diff --git a/src/geoopt_driver.f90 b/src/geoopt_driver.f90
index 69cfd7fc1..55985c4e0 100644
--- a/src/geoopt_driver.f90
+++ b/src/geoopt_driver.f90
@@ -35,34 +35,51 @@ subroutine geometry_optimization &
use xtb_setparam
+ ! Declaration section
implicit none
-
+
character(len=*), parameter :: source = 'xtb_geoopt'
-
- !> Calculation environment
+
+ !> Dummy-argument list
+ class(TCalculator), intent(inout) :: calc
+ !! instance of calculatior
type(TEnvironment), intent(inout) :: env
-
+ !! calculation environment
type(TMolecule), intent(inout) :: mol
+ !! molecular information
+ type(TRestart),intent(inout) :: wfn
+ !! wavefunction
integer, intent(in) :: tight
+ !! optimization level
integer, intent(in) :: maxiter
+ !! max number of SCC cycles
integer, intent(in) :: maxcycle_in
- type(TRestart),intent(inout) :: wfn
- class(TCalculator), intent(inout) :: calc
+ !! max number of optimization cycles
real(wp),intent(inout) :: etot
+ !! total energy
real(wp),intent(in) :: et
+ !! electronic temprature
real(wp),intent(inout) :: egap
+ !! HOMO-LUMO gap
real(wp),intent(inout) :: g(3,mol%n)
+ !! gradients
real(wp),intent(inout) :: sigma(3,3)
+ !! strain derivatives
logical, intent(in) :: pr
+ !! printlevel
logical, intent(out) :: fail
+ !! true-> optimization not converged
logical, intent(in) :: initial_sp
+ !! true-> perform initial single point
+ !> Local variables
type(scc_results) :: res
logical :: final_sp, exitRun
integer :: printlevel
integer :: ilog
final_sp = pr
+
if (pr) then
call delete_file('NOT_CONVERGED')
call delete_file('.xtboptok')
@@ -70,9 +87,11 @@ subroutine geometry_optimization &
else
printlevel = 0
endif
-
+
+ !> create optimization log file
if (.not.allocated(set%opt_logfile)) then
call open_file(ilog,'xtbopt.log','w')
+ !! default
else
if (set%opt_logfile == '-') then
ilog = stdout
@@ -82,22 +101,27 @@ subroutine geometry_optimization &
"Optimization log is written to '"//set%opt_logfile//"'"
endif
endif
+
+ if (set%oniom_settings%logs) then
+ call open_file(set%oniom_settings%ilog1,"low.inner_region.log",'w')
+ call open_file(set%oniom_settings%ilog2,"high.inner_region.log",'w')
+ endif
call mol%update
+ !! update interatomic and minimum image(periodic) distances
- if (initial_sp) then
- call singlepoint &
+ if (initial_sp) call singlepoint &
&(env,mol,wfn,calc, &
& egap,et,maxiter,printlevel-1,.false.,.false.,1.0_wp,etot,g,sigma,res)
- endif
+ !> actual calculation
select case(set%opt_engine)
case(p_engine_rf)
call ancopt &
&(env,ilog,mol,wfn,calc, &
& egap,et,maxiter,maxcycle_in,etot,g,sigma,tight,pr,fail)
- ! required since ANCopt might perform an untracked displacement
final_sp = .true.
+ !! required since ANCopt might perform an untracked displacement
case(p_engine_lbfgs)
call l_ancopt &
&(env,ilog,mol,wfn,calc, &
@@ -107,7 +131,7 @@ subroutine geometry_optimization &
&(env,ilog,mol,wfn,calc, &
& tight,maxcycle_in,etot,egap,g,sigma,printlevel,fail)
end select
-
+
if (pr) then
if (fail) then
call touch_file('NOT_CONVERGED')
@@ -118,11 +142,13 @@ subroutine geometry_optimization &
end if
call env%check(exitRun)
+ !! check if any errors occur
+
if (exitRun) then
call env%rescue("Trying to recover from failed geometry optimization", source)
fail = .true.
end if
-
+
if (pr.and.set%pr_finalstruct) then
write(env%unit,'(''================'')')
write(env%unit,*) 'final structure:'
@@ -130,15 +156,22 @@ subroutine geometry_optimization &
call writeMolecule(mol, env%unit)
endif
if (pr.and.set%pr_geosum) call geosum(mol%n,mol%at,mol%xyz)
-
+
+ !> usually final SP is performed for all three engines
if (final_sp) then
if (pr) call generic_header(env%unit,'Final Singlepoint',49,10)
call singlepoint &
&(env,mol,wfn,calc, &
& egap,et,maxiter,printlevel,.false.,.false.,1.0_wp,etot,g,sigma,res)
endif
-
+
if (ilog .ne.stdout) call close_file(ilog)
+ !! close log file if it is not printed out on stdout
+
+ if (set%oniom_settings%logs) then
+ call close_file(set%oniom_settings%ilog1)
+ call close_file(set%oniom_settings%ilog2)
+ endif
end subroutine geometry_optimization
diff --git a/src/io/writer.f90 b/src/io/writer.f90
index d86f842e7..0dafbc504 100644
--- a/src/io/writer.f90
+++ b/src/io/writer.f90
@@ -24,28 +24,37 @@ module xtb_io_writer
use xtb_mctc_filetypes, only : fileType
use xtb_mctc_version, only : version
use xtb_type_molecule, only : TMolecule, assignment(=)
+
implicit none
private
-
public :: writeMolecule
-
contains
-
-
+!------------------------------------------------------------------
+! To write molecular structure data into a file
+!------------------------------------------------------------------
subroutine writeMolecule(self, unit, format, energy, gnorm, number)
+
class(TMolecule), intent(in) :: self
+ !! dynamic type of TMolecule
integer, intent(in) :: unit
+ !! IO unit for file
integer, intent(in), optional :: format
+ !! the format: xyz, pdb, sdf, etc.
real(wp), intent(in), optional :: energy
+ !! energy output for comment line
real(wp), intent(in), optional :: gnorm
+ !! gradient norm output for comment line
integer, intent(in), optional :: number
+ !!
+
character(len=:), allocatable :: comment_line
character(len=20) :: energy_line
character(len=20) :: gnorm_line
type(structure_type) :: struc
type(error_type), allocatable :: error
integer :: ftype
+
if (present(format)) then
ftype = format
else
diff --git a/src/main/setup.f90 b/src/main/setup.f90
index 75cba40a2..b47438673 100644
--- a/src/main/setup.f90
+++ b/src/main/setup.f90
@@ -80,10 +80,20 @@ subroutine newCalculator(env, mol, calc, fname, restart, accuracy, input, iff_da
type(TDriverCalculator), allocatable :: driver
logical :: exitRun
-
+
select case(set%mode_extrun)
case default
call env%error("Unknown calculator type", source)
+
+ case(p_ext_oniom)
+ if (.not.present(input)) then
+ call env%error("ONIOM calculator requires input", source)
+ return
+ end if
+ allocate(oniom)
+ call newOniomCalculator(oniom, env, mol, input)
+ call move_alloc(oniom, calc)
+
case(p_ext_eht, p_ext_xtb)
allocate(xtb)
@@ -149,26 +159,21 @@ subroutine newCalculator(env, mol, calc, fname, restart, accuracy, input, iff_da
allocate(orca)
call newOrcaCalculator(orca, env, set%ext_orca)
call move_alloc(orca, calc)
+
case(p_ext_mopac)
allocate(mopac)
call newMopacCalculator(mopac, env, set%ext_mopac)
call move_alloc(mopac, calc)
+
case(p_ext_turbomole)
allocate(turbo)
call newTMCalculator(turbo, set%extcode, set%extmode)
call move_alloc(turbo, calc)
+
case(p_ext_driver)
allocate(driver)
call newDriverCalculator(driver, env, set%ext_driver)
call move_alloc(driver, calc)
- case(p_ext_oniom)
- if (.not.present(input)) then
- call env%error("ONIOM calculator requires input", source)
- return
- end if
- allocate(oniom)
- call newOniomCalculator(oniom, env, mol, input)
- call move_alloc(oniom, calc)
end select
end subroutine newCalculator
diff --git a/src/mctc/systools.F90 b/src/mctc/systools.F90
index 5394f1f6e..8c3555a4b 100644
--- a/src/mctc/systools.F90
+++ b/src/mctc/systools.F90
@@ -29,8 +29,13 @@ module xtb_mctc_systools
contains
+!> To read line
subroutine getline(unit,line,iostat)
+
use, intrinsic :: iso_fortran_env, only : iostat_eor
+ !! iostat_eor - assigned to the iostat variable if eof occurs
+ implicit none
+
integer,intent(in) :: unit
character(len=:),allocatable,intent(out) :: line
integer,intent(out),optional :: iostat
diff --git a/src/oniom.f90 b/src/oniom.f90
index 68f3df89f..7a0423859 100644
--- a/src/oniom.f90
+++ b/src/oniom.f90
@@ -14,7 +14,11 @@
! You should have received a copy of the GNU Lesser General Public Licen
! along with xtb. If not, see .
+!-----------------------------------------------------------------------
+!> ONIOM implementaion
+!-----------------------------------------------------------------------
module xtb_oniom
+
use xtb_mctc_accuracy, only: wp
use xtb_mctc_convert, only: aatoau
use xtb_type_atomlist, only: TAtomlist, len
@@ -34,40 +38,46 @@ module xtb_oniom
implicit none
private
- public :: TOniomCalculator, newOniomCalculator, oniom_input
-
+ public :: TOniomCalculator, newOniomCalculator, oniom_input, calculateCharge
+
+ !> To handle CML arguments
type :: oniom_input
character(len=:), allocatable :: first_arg
character(len=:), allocatable :: second_arg
- character(len=:), allocatable :: chrg
+ !character(len=:), allocatable :: chrg
+ !logical :: g
end type oniom_input
+
+ !> ONIOM calculator
+ type, extends(TCalculator) :: TOniomCalculator
- type, extends(TCalculator) :: TOniomCalculator !> new calculator type
-
- !> The method specified by the usernan
integer :: method_low, method_high
+ !! methods specified by the user, or default: gfnff, gfn2
- !> Charge, if present
- integer :: chrg_real, chrg_model
+ integer :: chrg_model
+ !! charge for inner region
+ !! rewrite
- type(TAtomList) :: list !> list of atoms in inner region
+ type(TAtomList) :: list
+ !! list of atoms in inner region
- !> QM 1, low level, whole region
class(TCalculator), allocatable :: real_low
+ !! whole region low level
- !> QM 1, low level, inner region
class(TCalculator), allocatable :: model_low
+ !! inner region low level
- !> QM 2, high level, inner region
class(TCalculator), allocatable :: model_high
+ !! inner region high level
- !> Index list of the inner region
integer, allocatable :: idx(:)
-
- !> Bond topology
+ !! index list of the inner region
+
type(TTopology), allocatable :: topo
+ !! topology used for cutting bonds
type(TRestart) :: chk_low, chk_high
+ !! wavefunctions for inner region calculations
contains
@@ -87,53 +97,63 @@ module xtb_oniom
module procedure :: new_coordinates
end interface
+
contains
-!> Create new calculator
+!--------------------------------------------------
+! Create ONIOM Calcultor
+!--------------------------------------------------
subroutine newOniomCalculator(self, env, mol, input)
- !> Calculator instance
+
+ implicit none
type(TOniomCalculator), intent(out) :: self
- !> Computation environment
+ !! new calculator created in this routine
type(TEnvironment), intent(inout) :: env
- !> Molecular structure data
+ !! calculation environment
type(TMolecule), intent(in) :: mol
- !> Restart checkpoint
- !type(TRestart), intent(inout) :: chk
- !> Input for Oniom calculation
+ !! molecular structure data
type(oniom_input), intent(in) :: input
-
+ !! CML input
+
+ !> Local variables
type(TxTBCalculator), allocatable :: xtb
type(TOrcaCalculator), allocatable :: orca
type(TGFFCalculator), allocatable :: gff
integer :: icol
integer :: i
-
!> Method identification
if (allocated(input%first_arg)) then
!! if methods specified
+
icol = index(input%first_arg, ':')
if (icol == 0) then
call env%error("Invalid method '"//input%first_arg//"' provided")
return
end if
+
self%method_high = string_to_id(input%first_arg(:icol - 1))
self%method_low = string_to_id(input%first_arg(icol + 1:))
+
if (self%method_high < 0 .or. self%method_high > 5) then
- call env%error("Invalid inner method")
+ call env%error("Invalid high-level method")
return
end if
+
if (self%method_low < 0 .or. self%method_low > 3) then
- call env%error("Invalid outer method")
+ call env%error("Invalid low-level method")
return
end if
+
else
!! default, gfn2:gfnff
+
self%method_high = 2
self%method_low = 3
- end if
-
-
+
+ endif
+
+ !> Write user-defined inner region list into array
self%list = TAtomList(list=input%second_arg)
call self%list%to_list(self%idx)
@@ -142,20 +162,22 @@ subroutine newOniomCalculator(self, env, mol, input)
return
end if
- ! Check if the user-defined inner region is valid
+ !> Check if the user-defined inner region is valid
if (any(self%idx < 1) .or. any(self%idx > mol%n)) then
- call env%error('The inner region specification is not correct')
+ call env%error('The specification of inner region is not correct')
return
end if
- ! The allocation of new calculator
+ !> Whole system calculator allocation
select case (self%method_low)
case default
+ !! GFN1/2
allocate (xtb)
call newXTBCalculator(env, mol, xtb, method=self%method_low)
call move_alloc(xtb, self%real_low)
case (3)
+ !! GFN-FF
allocate (gff)
call newGFFCalculator(env, mol, gff, "", .false.)
call move_alloc(gff, self%real_low)
@@ -164,179 +186,209 @@ subroutine newOniomCalculator(self, env, mol, input)
end subroutine newOniomCalculator
+!---------------------------------------------------------
+! 3 singlepoint energy calculations
+!---------------------------------------------------------
subroutine singlepoint(self, env, mol, chk, printlevel, restart, energy, gradient, sigma, hlgap, results)
- !? Is it right? What is inside init function? What interface/module procedure do?
+ use xtb_io_writer
+ use xtb_mctc_filetypes, only : fileType
use xtb_xtb_calculator, only: TxTBCalculator
- use xtb_setmod, only: set_chrg
+ implicit none
+ !> Dummy-argument list
class(TOniomCalculator), intent(inout) :: self
-
+ !! instance of TOniomCalculator
type(TEnvironment), intent(inout) :: env
-
+ !! calculation environment
type(TMolecule), intent(inout) :: mol
-
+ !! molecular structure data
type(TRestart), intent(inout) :: chk
-
- integer, intent(in) :: printlevel !> Print level for IO
-
+ !! restart data wrapper
+ integer, intent(in) :: printlevel
+ !! print level for IO
logical, intent(in) :: restart
-
+ !! if restarted
real(wp), intent(out) :: energy
-
+ !! ONIOM total energy
real(wp), intent(out) :: gradient(:, :)
-
- real(wp), intent(out) :: sigma(:, :) !> Strain derivatives
-
+ !! ONIOM gradients
+ real(wp), intent(out) :: sigma(:, :)
+ !! strain derivatives
real(wp), intent(out) :: hlgap
-
+ !! HOMO-LUMO gap
type(scc_results), intent(out) :: results
-
- type(TxTBCalculator), allocatable :: tmp !> Creating dummy object
- type(scc_results) :: results_low, results_high
+ !! final output stream
+
+ !> Local data
+ !> Temporary storages
+ type(TxTBCalculator), allocatable :: tmp
type(TGFFCalculator), allocatable :: gff
type(TOrcaCalculator), allocatable :: orca
type(TTMCalculator), allocatable :: turbo
+
+ !> For inner region high- and low-level calculation
+ type(scc_results) :: results_low, results_high
type(TMolecule) :: inner_mol
- logical :: exitRun
real(wp), allocatable :: gradient_low(:, :), gradient_high(:, :)
real(wp) :: energy_model_low, energy_model_high, hlgap_low, hlgap_high
real(wp) :: sigma_low(3, 3), sigma_high(3, 3)
- integer :: i
- !type(TRestart), allocatable :: chk0
- ! Check whether the calculator is initialized
+ real(wp),allocatable :: arr_gh(:), arr_gl(:)
+ !! 1-dim arrays for matmul of gradient matrix and Jacobian
+ real(wp), allocatable :: jacobian(:,:)
+ !! Jacobian matrix
+ integer,allocatable :: idx2(:)
+ integer :: i, coord_unit
+ logical :: exitRun
+ !! if any errors occur
+
+ !> Check whether the calculator is initialized
if (.not. allocated(self%real_low)) then
call env%error("Outer region calculator not provided")
return
end if
- ! Forward solvation to outer region
+ !> Forward solvation to outer region
if (allocated(self%solvation)) then
call move_alloc(self%solvation, self%real_low%solvation)
end if
+
+ !> First singlepoint of whole system with low-level
+ if (.not.set%oniom_settings%cut_inner) then
+
+ if (printlevel > 0) then
+ write(env%unit,'(a)')
+ write(env%unit,'(2x,72("-"))')
+ write (env%unit, '(/6x,a/)') "Singlepoint calculation of whole system with low-level method"
+ write(env%unit,'(2x,72("-"))')
+ write(env%unit,'(a)')
+ endif
+
+ call self%real_low%singlepoint(env, mol, chk, printlevel, restart, &
+ & energy, gradient, sigma, hlgap, results)
+
+ endif
+
+ call self%cutbond(env, mol, chk, self%topo, inner_mol,jacobian,idx2)
+ !! creating Linked Atoms
+ inner_mol%chrg = real(set%oniom_settings%innerchrg)
+ !! define inner region charge
- ! Check whether the low-level calculator needs a wavefunction
- if (.not.allocated(chk%wfn%qsh)) then
- ! print *, 'overwritten' - for restart
- select type(xtb => self%real_low)
- type is(TxTBCalculator)
- call newWavefunction(env, mol, xtb, chk)
- end select
- end if
-
-
- ! Perform calculation on outer region
- call self%real_low%singlepoint(env, mol, chk, printlevel, restart, &
- & energy, gradient, sigma, hlgap, results)
-
- !> check for charges
- call calculateCharge(self, env, mol, chk)
-
-
- !> Creating Linked atoms
- call self%cutbond(env, mol, chk, self%topo, inner_mol)
+ !> --cut flag termination
+ if (set%oniom_settings%cut_inner) then
+
+ write(env%unit,'(a)')
+ write(env%unit,'(2x,72("-"))')
+ write(env%unit,'(2x,"|",24x,a,1x,i0,22x,"|")') "INNER REGION CHARGE = ", nint(inner_mol%chrg)
+ write(env%unit,'(2x,72("-"))')
+ write(env%unit,'(a)')
+ call terminate(0)
+
+ endif
+
call env%check(exitRun)
if (exitRun) then
- call env%error("Could not create linked atoms")
+ call env%error("Could not create the linked atoms")
return
end if
- inner_mol%chrg = float(self%chrg_model)
- ! Inner region, low method
+ !> Inner region, low level method
if (.not. allocated(self%model_low)) then
- if (printlevel > 0) &
- write (env%unit, '(/a/)') "Initializing low level method"
-
select case (self%method_low)
case default
call env%error("Invalid low-level inner method")
+ return
case (1, 2)
-
+ !! GFN1/2
allocate (tmp)
call newXTBCalculator(env, inner_mol, tmp, method=self%method_low)
call move_alloc(tmp, self%model_low)
case (3)
-
+ !! GFN-FF
allocate (gff)
call newGFFCalculator(env, inner_mol, gff, "", .false.)
call move_alloc(gff, self%model_low)
-
+
end select
- if (allocated(self%real_low%solvation)) then
- allocate(self%model_low%solvation)
- self%model_low%solvation = self%real_low%solvation
- end if
call env%check(exitRun)
if (exitRun) then
- call env%error("Could not setup low level method")
+ call env%error("Could not setup low-level method")
return
end if
+
end if
-
- select type (calc => self%model_low)
- type is (TxTBCalculator)
- call newWavefunction(env, inner_mol, calc, self%chk_low)
- end select
-
- ! Inner region, high method
+
+ !> Inner region, high-level method
if (.not. allocated(self%model_high)) then
- if (printlevel > 0) &
- write (env%unit, '(/a/)') "Initializing high level method"
-
select case (self%method_high)
+ case default
+ call env%error("Invalid high-level inner method")
+ return
+
case (1, 2)
+ !! GFN1/2
allocate (tmp)
call newXTBCalculator(env, inner_mol, tmp, method=self%method_high)
call move_alloc(tmp, self%model_high)
case (3)
-
+ !! GFN-FF
allocate (gff)
call newGFFCalculator(env, inner_mol, gff, "", .false.)
call move_alloc(gff, self%model_high)
case (4)
-
+ !! ORCA
allocate (orca)
- call newOrcaCalculator(orca, env, set%ext_orca)
+ call newOrcaCalculator(orca, env, set%ext_orca,oniom=.true.)
call move_alloc(orca, self%model_high)
- case (5)
-
+ case (5)
+ !! TURBOMOLE
allocate (turbo)
call newTMCalculator(turbo, 1, 1)
call move_alloc(turbo, self%model_high)
-
+ call protectCoord(env)
+ !! to copy coord file into origin.coord
end select
- if (allocated(self%real_low%solvation)) then
-
- allocate(self%model_high%solvation)
- self%model_high%solvation = self%real_low%solvation
- end if
call env%check(exitRun)
if (exitRun) then
- call env%error("Could not setup high level method")
+ call env%error("Could not setup high-level method")
return
end if
end if
+
+ !> To setup partial and shell charges for GFN1/2
+ if (.not.allocated(self%chk_low%wfn%qsh)) then
+ select type (calc => self%model_low)
+ type is (TxTBCalculator)
+ call newWavefunction(env, inner_mol, calc, self%chk_low)
+ end select
+ end if
+ if (.not.allocated(self%chk_high%wfn%qsh)) then
+ select type (calc => self%model_high)
+ type is (TxTBCalculator)
+ call newWavefunction(env, inner_mol, calc, self%chk_high)
+ end select
+ end if
+
- select type (calc => self%model_high)
- type is (TxTBCalculator)
- call newWavefunction(env, inner_mol, calc, self%chk_high)
- end select
-
- !> Inner region low method, actual calculatoin
-
- if (printlevel > 0) &
- write (env%unit, '(/a/)') "SP of inner region with low method"
+ !> Inner region low-level method, singlepoint calculation
+ if (printlevel > 0) then
+ write(env%unit,'(a)')
+ write(env%unit,'(2x,72("-"))')
+ write (env%unit, '(/6x,a/)') "Singlepoint calculation of inner region with low-level method"
+ write(env%unit,'(2x,72("-"))')
+ write(env%unit,'(a)')
+ endif
allocate (gradient_low(3, inner_mol%n))
energy_model_low = 0.0_wp
@@ -345,10 +397,15 @@ subroutine singlepoint(self, env, mol, chk, printlevel, restart, energy, gradien
call self%model_low%singlepoint(env, inner_mol, self%chk_low, printlevel, restart, &
& energy_model_low, gradient_low, sigma_low, hlgap_low, results_low)
- !> Inner region high method, actual calculatoin
- if (printlevel > 0) &
- write (env%unit, '(/a/)') "SP of inner region with high method"
+ !> Inner region high-level method, singlepoint calculation
+ if (printlevel > 0) then
+ write(env%unit,'(a)')
+ write(env%unit,'(2x,72("-"))')
+ write (env%unit, '(/6x,a/)') "Singlepoint calculation of inner region with high-level method"
+ write(env%unit,'(2x,72("-"))')
+ write(env%unit,'(a)')
+ endif
allocate (gradient_high(3, inner_mol%n))
energy_model_high = 0.0_wp
@@ -356,14 +413,42 @@ subroutine singlepoint(self, env, mol, chk, printlevel, restart, energy, gradien
call self%model_high%singlepoint(env, inner_mol, self%chk_high, printlevel, restart, &
& energy_model_high, gradient_high, sigma_high, hlgap_high, results_high)
+
- results%dipole = results%dipole + results_high%dipole - results_low%dipole
- energy = energy + energy_model_high - energy_model_low !> The actual Oniom energy
- do i = 1, size(self%idx)
- gradient(:, self%idx(i)) = gradient(:, self%idx(i)) + gradient_high(:, i) - gradient_low(:, i)
- end do
+ !> Write opt logs for inner region
+ if (set%oniom_settings%logs)then
+ call writeMolecule(inner_mol, set%oniom_settings%ilog1, format=filetype%xyz,energy=energy_model_low)
+ call writeMolecule(inner_mol, set%oniom_settings%ilog2, format=filetype%xyz,energy=energy_model_high)
+ endif
+
+ results%dipole = results%dipole + results_high%dipole - results_low%dipole
+ sigma = sigma - sigma_low + sigma_high
+ results%hl_gap = hlgap - hlgap_low + hlgap_high
+
+ !> ONIOM energy
+ energy = energy + energy_model_high - energy_model_low
results%e_total = energy
+
+ !> [gradient*Jacobian] with forward and backward transformation
+ !! rewrite for hessian
+ call matrix_to_array(gradient_high,arr_gh)
+ call matrix_to_array(gradient_low,arr_gl)
+
+ arr_gh=matmul(jacobian,arr_gh)
+ arr_gl=matmul(jacobian,arr_gl)
+
+ call array_to_matrix(arr_gh,gradient_high)
+ call array_to_matrix(arr_gl,gradient_low)
+
+
+ !> ONIOM gradients
+ do i = 1, size(idx2)
+ gradient(:, idx2(i)) = gradient(:, idx2(i)) + gradient_high(:, i) - gradient_low(:, i)
+ end do
+ results%gnorm=norm2(gradient)
+ deallocate(gradient_high)
+ deallocate(gradient_low)
end subroutine singlepoint
@@ -387,6 +472,9 @@ subroutine hessian(self, env, mol0, chk0, list, step, hess, dipgrad)
!> Array to add dipole gradient to
real(wp), intent(inout) :: dipgrad(:, :)
+ real(wp), allocatable :: jacobian(:,:)
+ integer,allocatable :: idx2(:)
+
integer :: ii, jj, ic, jc, im, jm
real(wp), allocatable :: hess_model(:, :), dipgrad_model(:, :)
type(TMolecule) :: mol_model
@@ -397,9 +485,8 @@ subroutine hessian(self, env, mol0, chk0, list, step, hess, dipgrad)
& hess, dipgrad)
!call chk0%wfn%allocate(mol0%n,chk0%basis%nshell,chk0%real_low%basis%nao)
- print*,"cutbond"
! Creating Linked atoms
- call self%cutbond(env, mol0, chk0, self%topo, mol_model)
+ call self%cutbond(env, mol0, chk0, self%topo, mol_model,jacobian,idx2)
mol_model%chrg = float(self%chrg_model)
list_model = [(ii, ii = 1, size(self%idx))]
@@ -453,144 +540,217 @@ subroutine writeInfo(self, unit, mol)
end subroutine writeInfo
-subroutine cutbond(self, env, mol, chk, topo, inner_mol)
+!---------------------------------------------------------------------
+! Create inner region
+!---------------------------------------------------------------------
+subroutine cutbond(self, env, mol, chk, topo, inner_mol,jacobian,idx2)
- !> to initialize mol object
use xtb_type_molecule, only: init
-
- !> To get topology info from wfn
use xtb_topology, only: makeBondTopology, topologyToNeighbourList
+ use xtb_io_writer
+ use xtb_mctc_filetypes, only : fileType
+
+ implicit none
+ character(len=*), parameter :: source = "xtb_oniom_cutbond"
- !> actual calculator
+ !> Dummy-argument list
class(TOniomCalculator), intent(in) :: self
-
- !> Environment
+ !! polymorphic calculator
type(TEnvironment), intent(inout) :: env
-
- !> Wavefunction
+ !! calculation environment
type(TRestart), intent(in) :: chk
-
- !> Molecular data storage
+ !! wavefunction wrapper
type(TMolecule), intent(in) :: mol
-
- !> Inner region geometry
+ !! molecular structure data
type(TMolecule), intent(inout) :: inner_mol
-
- !> Some buffer geometry
- type(TMolecule) :: cash
-
- !> topology info
+ !! inner region
type(TTopology), allocatable, intent(inout) :: topo
-
- !> neighbour lists
+ !! topology info
+ real(wp),allocatable, intent(inout) :: jacobian(:,:)
+ !! jacobian matrix
+ integer,allocatable, intent(out) :: idx2(:)
+ !! list of inner region atom indices + host atoms
+
+ !> Local variables
type(TNeighbourList) :: neighList
-
- integer :: i, j, n, k, pre_last,iterator
-
- !> Arrays for inner region atoms and their atomic numbers
+ !! neighbour lists
+ integer, allocatable :: brokenBondPairs(:,:)
+ !! pair-wise indices for cutbound
integer, allocatable :: at(:), at2(:)
+ integer, allocatable :: bonded(:, :)
real(wp), allocatable :: xyz(:, :), xyz2(:, :)
+ character(len=:),allocatable :: fname_inner
logical :: inside
- integer, allocatable :: bonded(:, :)
+ !! if the both bonded atoms are inside the inner region
+ integer :: i, j, n, k, pre_last,iterator
+ integer :: io
+
+
inside = .FALSE.
n = len(self%list)
+ !! initial number of atoms in inner region without LAs
+
+ !> Allocate accordingly the basic molecular data
allocate (at2(size(self%idx)))
allocate (xyz2(3, size(self%idx)))
allocate (at(n))
allocate (xyz(3, n))
-
- !> Assignment of initial inner region
+ allocate (idx2(size(self%idx)))
+
+ idx2=self%idx
+ !! to save inner atom list in the matching array
+
+ !> 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
-
- !> To identify bonded atoms and save them into an array + iterator
+
+ call create_jacobian(jacobian,at)
+ !! initialiaze jacobian as identity matrix
+
+ !> To identify bonded atoms and save them into an array + assign iterator
select type (calc => self%real_low)
+ class default
+ call env%error("Topology information could not be derived from the given calculator",source)
+ return
+
type is (TGFFCalculator)
- bonded = calc%topo%blist !> automatic allocation
+ !! GFN-FF
+ bonded = calc%topo%blist
+ !! bonded atom list
iterator = size(bonded,2)
+ !! iterator = number of bonds
type is (TxTBCalculator)
+ !! GFN1/2
if (.not. allocated(topo)) then
allocate (topo)
- call makeBondTopology(topo, mol, chk%wfn%wbo) !> return assigned topo%list
-
- call topologyToNeighbourList(topo, neighList, mol) !> return neighList
+ call makeBondTopology(topo, mol, chk%wfn%wbo)
+ !! return assigned topo%list
+ call topologyToNeighbourList(topo, neighList, mol)
+ !! return neighList
end if
+
allocate (bonded(2, len(topo)))
do i = 1, len(topo)
bonded(:, i) = topo%list(1:2, i)
end do
+
iterator = size(bonded,2)
+ !! iterator = number of bonds
end select
- !> Actual bond cutting and creating linked atom
+
+ !> Algorithm to identify broken covalent bonds due to the ONIOM boundary
do i = 1, size(self%idx)
+ !! iterate for all atoms in the user-provided list
do j = 1, iterator
+ !! iterate through pair of atoms that are bonded
if (bonded(1, j) == self%idx(i)) then
+ !! if atom in the list is bonded
do k = 1, size(self%idx)
+ !! iterate again through the list
if (self%idx(k) == bonded(2, j)) then
+ !! if it is found -> inside inner region
inside = .TRUE.
end if
end do
+
if (.not. inside) then
+ !! if bond is broken
+
+ !> Check if single bond is broken
select type (calc => self%real_low)
class default
call checkfororder(env, mol, self%idx(i), bonded(2, j), bond=topo%list(3, j))
type is (TGFFCalculator)
call checkfororder(env, mol, self%idx(i), bonded(2, j), hybrid=calc%topo%hyb)
end select
+
+ !> Increase no of corresponding arrays by 1
pre_last = size(at)
call resize(at)
+ call resize(idx2)
+ idx2(pre_last+1) = bonded(2,j)
+ !! to save index of the host atom
at(pre_last + 1) = 1
+ !! transform the newly created atom to hydrogen atom(Link Atom)
call resize(xyz)
- call coord(mol, xyz, at(pre_last), self%idx(i), bonded(2, j))
+ !! sdjust accordingly the coordinate matrix
+ call resize_jacobian(jacobian)
+ !! increase matrix size in both directions
+ call coord(env,mol,xyz,self%idx(i),bonded(2,j),jacobian,i)
+ !! determine new position of added H atom
+
end if
inside = .FALSE.
else if (bonded(2, j) == self%idx(i)) then
+ !! if atom in the list is bonded
do k = 1, size(self%idx)
+ !! iterate again through the list
if (self%idx(k) == bonded(1, j)) then
inside = .TRUE.
+ !! if it is found -> inside inner region
end if
end do
if (.not. inside) then
+ !! if bond is broken
select type (calc => self%real_low)
+
+ !> Check if single bond is broken
class default
call checkfororder(env, mol, self%idx(i), bonded(1, j), bond=topo%list(3, j))
type is (TGFFCalculator)
call checkfororder(env, mol, self%idx(i), bonded(1, j), hybrid=calc%topo%hyb)
end select
+
+ !> Increase no of corresponding arrays by 1
pre_last = size(at)
call resize(at)
+ call resize(idx2)
+ idx2(pre_last+1) = bonded(1,j)
+ !! to save index of the host atom
at(pre_last + 1) = 1
+ !! transform the newly created atom to hydrogen atom(Link Atom)
call resize(xyz)
- call coord(mol, xyz, at(pre_last), self%idx(i), bonded(1, j))
+ !! adjust accordingly the coordinate matrix
+ call resize_jacobian(jacobian)
+ !! increase matrix size in both directions
+ call coord(env,mol, xyz, self%idx(i), bonded(1, j),jacobian,i)
+ !! determine new position of added H atom
end if
inside = .FALSE.
end if
end do
end do
- !call init(cash,at2,xyz2)
+
call init(inner_mol, at, xyz)
- !block
- !use xtb_io_writer
- !use xtb_mctc_filetypes, only : fileType
- !integer :: io
- !call open_file(io, "inner-region_with_h.xyz", "w")
- !call writeMolecule(inner_mol, io, filetype%xyz)
- !call close_file(io)
- !stop
- !end block
+ !! initialize mol object
+
+ !> Create Xmol file for inner region
+ if (set%oniom_settings%cut_inner) then
+ fname_inner = "inner_region_without_h.xyz"
+ else
+ fname_inner = "inner_region.xyz"
+ endif
+ call open_file(io, fname_inner, "w")
+ call writeMolecule(inner_mol, io, filetype%xyz)
+ call close_file(io)
end subroutine cutbond
+!---------------------------------------
+! increase atomic number array by 1
+!---------------------------------------
subroutine new_atom(at)
-
+
+ implicit none
integer, allocatable, intent(inout) :: at(:)
integer, allocatable :: tmp2(:)
+
allocate (tmp2(size(at) + 1))
tmp2(:size(at)) = at(:size(at))
deallocate (at)
@@ -598,10 +758,14 @@ subroutine new_atom(at)
end subroutine new_atom
+!--------------------------------------
+! increase coordinate matrix by 1
+!--------------------------------------
subroutine new_coordinates(xyz)
-
- real, allocatable :: xyz(:, :)
- real, allocatable :: tmp1(:, :)
+
+ implicit none
+ real, allocatable :: xyz(:,:)
+ real, allocatable :: tmp1(:,:)
integer :: atom_num
atom_num = size(xyz, 2)
@@ -612,51 +776,403 @@ subroutine new_coordinates(xyz)
end subroutine new_coordinates
-!call coord(mol, xyz, pre_last, at(self%idx(i)), self%idx(i), bonded(1, j))
-subroutine newcoord(mol, xyz, at, idx1, idx2)
+!------------------------------------
+! create identity matrix
+!------------------------------------
+subroutine create_jacobian(matrix,at)
+
+ implicit none
+ real(wp),allocatable :: matrix(:,:)
+ integer,intent(in) :: at(:)
+ integer :: i, j
+
+ allocate (matrix(size(at)*3,size(at)*3))
+ do i=1,size(at)*3
+ do j=1,size(at)*3
+ if (i==j) then
+ matrix(i,j)=1.0_wp
+ else
+ matrix(i,j)=0.0_wp
+ endif
+ enddo
+ enddo
+
+end subroutine create_jacobian
- type(TMolecule), intent(in) :: mol
+!----------------------------------------------------
+! To increase matrix dimensionality
+! (3 new entries in diagonal and subsequent increase)
+!----------------------------------------------------
+subroutine resize_jacobian(matrix)
+
+ implicit none
+ real(wp), allocatable :: matrix(:,:)
+ real(wp), allocatable :: tmp(:,:)
+ !! To store old values while reallocation
+ integer :: coord_num
+ !! The current number of coordinates
- real(wp), intent(inout) :: xyz(:, :)
+ coord_num = size(matrix,1)
- !integer, intent(in) :: pre_last
+ allocate(tmp(coord_num+3,coord_num+3))
+ tmp(:coord_num,:coord_num) = matrix(:coord_num,:coord_num)
+ deallocate(matrix)
+ call move_alloc(tmp, matrix)
- integer, intent(in) :: at, idx1, idx2
+end subroutine resize_jacobian
- real(wp) :: dist, prefactor
+!----------------------------------------------------
+! calculate new postion for LA and corresponding J
+!----------------------------------------------------
+subroutine newcoord(env,mol,xyz,idx1,idx2,jacobian,connectorPosition)
+
+ implicit none
+ !> Dummy-argument list
+ character(len=*), parameter :: source = "oniom_newcoord"
+ !! name of error producer routine
+ type(TEnvironment), intent(inout) :: env
+ !! calculation environment
+ type(TMolecule), intent(in) :: mol
+ !! molecular structure data
+ real(wp), intent(inout) :: xyz(:, :)
+ !! coordinate matrix
+ integer, intent(in) :: idx1
+ !! connector, one that stays in the inner region
+ integer, intent(in) :: idx2
+ !! host atom, the one being substitued
+ real(wp), intent(inout) :: jacobian(:,:)
+ !! jacobian maxtrix
+ integer, intent(in) :: connectorPosition
+ !! ordinal number of connector atom in mol%at
+
+ !> Local data
+ real(wp) :: dist
+ !! standard bond length between LAC(linked atom connector)-H
+ real(wp) :: dist2
+ !! standard bond length between LAC(linked atom connector)-LAH(linked atom host)
+ real(wp) :: prefactor
+ !! scaling parameter
+ real(wp) :: dist_12
+ !! vector length squared
real(wp) :: xyz1(3), xyz2(3)
-
- select case (at)
+ !! coordinates of the cleaved atom pair
+ character(len=:), allocatable :: warning
+ !! message showed in case of default parameter usage
+ character(len=:), allocatable :: warning2
+ !! message showed in case of the switcing between derived and fixed prefactor values
+ logical :: def
+ !! to check if default is used
+ logical, save :: rep=.false.
+ !! to control warnings
+ character(len=3) :: dummy1, dummy2
+ !! to write ordinal number as character
+ real(wp) :: xyz_difference(3)
+ !! the difference between new and old coordinates
+ integer :: i,j,k
+
+ write (dummy1, '(I3)') mol%at(idx1)
+ write (dummy2, '(I3)') mol%at(idx2)
+
+ !> initiale initialization
+ def = .false.
+ warning = "Atoms "//dummy1//" and "//dummy2//" are not accounted in the parameter suite(S,P,N,C,O), the default distance values will be used."
+ warning2 = "The distance between atoms "//dummy1//" and "//dummy2//" is almost the same, switching to fixed regime"
+
+ xyz1 = mol%xyz(:, idx1)
+ !! coordinates of connector atom
+ xyz2 = mol%xyz(:, idx2)
+ !! coordinates of host atom
+ dist_12=sum((xyz1 - xyz2)**2)
+ !! distance squared
+
+ !> To identify average bond distances
+ !> for connector-H and connector-host
+ select case (mol%at(idx1))
case default
- dist = 1.09*aatoau
- !> Carbon
+ def = .true.
+ dist = 1.084*aatoau
+ !! general case: C-H
+ dist2 = 1.528*aatoau
+ !! general case: C-C
+
+ case (1)
+ !! H
+ dist = 0.740*aatoau
+ !! H-H
+
+ select case (mol%at(idx2))
+ case default
+ dist2 = 1.084*aatoau
+ def = .true.
+ case(1)
+ dist2 = 0.740*aatoau
+ !! H-H
+ case (6)
+ dist2 = 1.084*aatoau
+ !! H-C
+ case (8)
+ dist2 = 0.964*aatoau
+ !! H-O
+ case (7)
+ dist2 = 1.024*aatoau
+ !! H-N
+ case (15)
+ dist2 = 1.414*aatoau
+ !! H-P
+ case (16)
+ dist2 = 1.389*aatoau
+ !! H-S
+ end select
+
case (6)
- dist = 1.09*aatoau
- !> Oxygen
+ !! C
+ dist = 1.084*aatoau
+ !! C-H
+
+ select case (mol%at(idx2))
+ case default
+ dist2 = 1.528*aatoau
+ def = .true.
+ case(1)
+ dist2 = 1.084*aatoau
+ !! C-H
+ case (6)
+ dist2 = 1.528*aatoau
+ !! C-C
+ case (8)
+ dist2 = 1.430*aatoau
+ !! C-O
+ case (7)
+ dist2 = 1.475*aatoau
+ !! C-N
+ case (15)
+ dist2 = 1.860*aatoau
+ !! C-P
+ case (16)
+ dist2 = 1.750*aatoau
+ !! C-S
+ end select
+
+ case (7)
+ !! N
+ dist = 1.024*aatoau
+ !! N-H
+
+ select case (mol%at(idx2))
+ case default
+ dist2 = 1.470*aatoau
+ def = .true.
+ case(1)
+ dist2 = 1.024*aatoau
+ !! N-H
+ case (6)
+ dist2 = 1.475*aatoau
+ !! N-C
+ case (8)
+ dist2 = 1.360*aatoau
+ !! N-O
+ case (7)
+ dist2 = 1.470*aatoau
+ !! N-N
+ case (15)
+ dist2 = 1.770*aatoau
+ !! N-P
+ case (16)
+ dist2 = 1.650*aatoau
+ !! N-S
+ end select
+
+
case (8)
+ !! O
dist = 0.964*aatoau
- !> Nitrogen
- case (7)
- dist = 1.024*aatoau
- !> Phosphorus
+ !! O-H
+
+ select case (mol%at(idx2))
+ case default
+ dist2 = 1.450*aatoau
+ def = .true.
+ case(1)
+ dist2 = 0.964*aatoau
+ !! O-H
+ case (6)
+ dist2 = 1.430*aatoau
+ !! O-C
+ case (8)
+ dist2 = 1.450*aatoau
+ !! O-O
+ case (7)
+ dist2 = 1.360*aatoau
+ !! O-N
+ case (15)
+ dist2 = 1.750*aatoau
+ !! O-P
+ case (16)
+ dist2 = 1.500*aatoau
+ !! O-S
+ end select
+
case (15)
+ !! P
dist = 1.414*aatoau
- !> Sulfur
+ !! P-H
+
+ select case (mol%at(idx2))
+ case default
+ dist2 = 1.770*aatoau
+ def = .true.
+ case(1)
+ dist2 = 1.414*aatoau
+ !! P-H
+ case (6)
+ dist2 = 1.860*aatoau
+ !! P-C
+ case (8)
+ dist2 = 1.750*aatoau
+ !! P-O
+ case (7)
+ dist2 = 1.770*aatoau
+ !! P-N
+ endselect
+
case (16)
+ !! S
dist = 1.389*aatoau
- end select
-
- xyz1 = mol%xyz(:, idx1)
- xyz2 = mol%xyz(:, idx2)
- prefactor = dist/sqrt(sum((xyz1 - xyz2)**2))
+ !! S-H
+
+ select case (mol%at(idx2))
+ case default
+ dist2 = 1.650*aatoau
+ def = .true.
+ case(1)
+ dist2 = 1.389*aatoau
+ !! S-H
+ case (6)
+ dist2 = 1.750*aatoau
+ !! S-C
+ case (8)
+ dist2 = 1.500*aatoau
+ !! S-O
+ case (7)
+ dist2 = 1.650*aatoau
+ !! S-N
+ case (16)
+ dist2 = 2.040*aatoau
+ !! S-S
+ endselect
+
+ end select
+
+ !> different ways of computing scaling parameter prefactor(k)
+ if (set%oniom_settings%derived) then
+ !! prefactor can change
+ prefactor = dist/sqrt(dist_12)
+ else
+ !! prefactor is fixed
+ prefactor = dist/dist2
+ !> if default values are used
+ if(def.and. .not.rep) then
+ rep=.true.
+ call env%warning(warning,source)
+ endif
+ endif
- !> new H coordinates
xyz(:, size(xyz, 2)) = xyz1 + (xyz2 - xyz1) * prefactor
+ !! LA (linked atom) coordinates
+
+ !> To determine if the difference between the coordinates of LA and LAH is small
+ !> if yes -> change from derived to fixed
+ xyz_difference=xyz2-xyz(:,size(xyz,2))
+ if (all(xyz_difference<1.0E-5).and.set%oniom_settings%derived) then
+ set%oniom_settings%derived=.false.
+ call env%warning(warning2,source)
+ prefactor = dist/dist2
+ endif
+
+ call derivative(jacobian,connectorPosition,size(xyz,2),prefactor,mol%xyz,idx1,idx2,dist_12,set%oniom_settings%derived)
+ !! take derivatives of model system coordinates wrt whole molecule
end subroutine newcoord
+!----------------------------------------------------------------------------
+! redefine Jacobian matrix for newly added atoms
+!----------------------------------------------------------------------------
+subroutine derivative(jacobian,con,link,prefactor,xyz,idx1,idx2,dist_12,derived)
+
+ implicit none
+ !> Dummy-argument list
+ real(wp), intent(inout) :: jacobian(:,:)
+ !! jacobian matrix
+ integer,intent(in) :: con
+ !! position of connector atom in the model system atom list
+ integer, intent(in) :: link
+ !! position of linked atom in the model system atom list
+ real(wp),intent(in) :: prefactor
+ !! scaling factor k
+ real(wp),intent(in) :: dist_12
+ !! square of distance between idx1 and idx2
+ real(wp),intent(in) :: xyz(:,:)
+ !! coordinates of whole molecule
+ integer, intent(in) :: idx1
+ !! connector, one that stays in the inner region
+ integer, intent(in) :: idx2
+ !! host atom, the one being substitued
+ logical, intent(in) :: derived
+ !! To fix value of prefactor variable
+ integer :: con3, link3
+ !! to account for all 3 coordinates
+ integer :: counter1(3), counter2(3)
+ !! to save the positions of changed matrix elements
+ integer :: i, j
+
+
+ con3=con*3
+ link3=link*3
+
+ do i=1,3
+ counter1(i)=con3-i+1
+ counter2(i)=link3-i+1
+ enddo
+
+ !> Assign all new matrix elements to 0
+ jacobian(counter2(3):,:) = 0.0_wp
+ jacobian(:,counter2(3):) = 0.0_wp
+
+ !> Algorithm to find non-zero elements of the matrix
+ do i=1,3
+
+ if(.not.derived) then
+ !> Fixed Jacobian
+ jacobian(counter1(i),counter2(i))=1-prefactor
+ jacobian(counter2(i),counter2(i))=prefactor
+
+ else
+ !> Derived Jacobian
+ if (i==1) then
+ !> x coordinate
+ jacobian(counter2(i),counter2(i)) =(prefactor*( (xyz(2,idx1)**2) - 2*xyz(2,idx1)*xyz(2,idx2) + (xyz(2,idx2)**2) + ((xyz(3,idx1)-xyz(3,idx2))**2) ))/(dist_12**(3.0_wp/2.0_wp))
+ else if (i==2) then
+ !> y coordinate
+ jacobian(counter2(i),counter2(i)) =(prefactor*( (xyz(1,idx1)**2) - 2*xyz(1,idx1)*xyz(1,idx2) + (xyz(1,idx2)**2) + ((xyz(3,idx1)-xyz(3,idx2))**2) ))/(dist_12**(3.0_wp/2.0_wp))
+ else
+ !> z coordinate
+ jacobian(counter2(i),counter2(i)) =(prefactor*( (xyz(1,idx1)**2) - 2*xyz(1,idx1)*xyz(1,idx2) + (xyz(1,idx2)**2) + ((xyz(2,idx1)-xyz(2,idx2))**2) ))/(dist_12**(3.0_wp/2.0_wp))
+ endif
+
+ jacobian(counter1(i),counter2(i)) = (prefactor*((xyz(i,idx2)-xyz(i,idx1))**2))/(dist_12**(3.0_wp/2.0_wp)) - (prefactor/(sqrt(dist_12))) + 1.0_wp
+
+ endif
+
+ enddo
+
+end subroutine derivative
+
+!--------------------------------------
+! To allocate method
+!--------------------------------------
function string_to_id(string) result(id)
+ implicit none
character(len=*), intent(in) :: string
integer :: id
@@ -668,10 +1184,8 @@ function string_to_id(string) result(id)
id = 2
case ('gfn1')
id = 1
- !> FF
case ('gfnff')
id = 3
- !> how it is specified in program??
case ('orca')
id = 4
case ('turbomole')
@@ -680,20 +1194,29 @@ function string_to_id(string) result(id)
end function string_to_id
+!----------------------------------------------------------
+! To check if (bond order > 1)
+!----------------------------------------------------------
subroutine checkfororder(env, mol, idx1, idx2, bond, hybrid)
-
- character(len=*), parameter :: source = 'oniom_cutbond'
-
+
+ implicit none
+ !> Dummy-argument list
+ character(len=*), parameter :: source = 'xtb_oniom_checkorder'
+ !! name of error producer routine
integer, intent(in), optional :: hybrid(:)
-
+ !! hybridization from GFN-FF; topo%hyb
integer, intent(in), optional :: bond
-
+ !! wiberg bond order
type(TEnvironment), intent(inout) :: env
-
+ !! calculation environment
type(TMolecule), intent(in) :: mol
-
- integer, intent(in) :: idx1, idx2
-
+ !! molecular staructure data
+ integer, intent(in) :: idx1
+ !! connector, one that stays in the inner region
+ integer, intent(in) :: idx2
+ !! host atom, the one being substitued
+
+ !> Local data
character(len=:), allocatable :: warning
integer :: b
character(len=5) :: dummy1, dummy2
@@ -701,6 +1224,7 @@ subroutine checkfororder(env, mol, idx1, idx2, bond, hybrid)
write (dummy1, '(I5)') idx1
write (dummy2, '(I5)') idx2
warning = "You are cutting a bond with the order higher than 1 between "//dummy1//" and "//dummy2
+
if (present(bond)) then
if (bond > 1) then
call env%error(warning, source)
@@ -717,43 +1241,125 @@ subroutine checkfororder(env, mol, idx1, idx2, bond, hybrid)
end subroutine checkfororder
-subroutine calculateCharge(self, env, mol, chk)
-
- character(len=*), parameter :: source = 'charge for inner region'
-
+! -------------------------------------------------
+! Automatic ONIOM inner region charge determination
+!--------------------------------------------------
+function calculateCharge(self, env, mol, chk) result(chrg_model)
+
+ implicit none
+ !> Dummy-argument list
+ character(len=*), parameter :: source = 'xtb_oniom_calculateCharge'
+ !! name of error producer routine
class(TOniomCalculator), intent(inout) :: self
-
+ !! polymorhic calculator
type(TEnvironment), intent(inout) :: env
-
+ !! calculation environment
type(TMolecule), intent(in) :: mol
-
+ !! molecular structure data
type(TRestart), intent(in) :: chk
-
+ !! wavefuntion wrapper
+ real(wp) :: charge
+ !! inner region charge
+
+ !> Local data
integer :: i, j, n, k, pre_last
+ integer :: chrg_model
integer, allocatable :: at(:)
- real(wp) :: charge
charge = 0.0_wp
select type (calc => self%real_low)
type is (TGFFCalculator)
-
+ !! GFN-FF
do i = 1, size(self%idx)
charge = charge + calc%topo%qa(self%idx(i))
-
end do
- type is (TxTBCalculator)
- do i = 1, size(self%idx)
+
+ type is (TxTBCalculator)
+ !! GFN1/2
+ do i = 1, size(self%idx)
charge = charge + chk%wfn%q(self%idx(i))
- end do
+
+ end do
class default
call env%error("Not possible to calculate with external methods for real region", source)
return
end select
+
+ chrg_model = nint(charge)
+
+end function calculateCharge
- self%chrg_model = nint(charge)
+!---------------------------------------------
+! To transform matrix into 1 dimensional array
+!---------------------------------------------
+subroutine matrix_to_array(mtrx,arr)
+
+ implicit none
+ real(wp), intent(in) :: mtrx (:,:)
+ real(wp), allocatable, intent(out) :: arr(:)
+ integer :: i, j, k
-end subroutine calculateCharge
+ allocate(arr(size(mtrx,1)*size(mtrx,2)))
+
+ k=1
+ do i=1, size(mtrx,2)
+ do j=1, size(mtrx,1)
+ arr(k)=mtrx(j,i)
+ k=k+1
+ enddo
+ enddo
+
+end subroutine matrix_to_array
+
+!--------------------------------------------
+! To transform 1 dimensional array to matrix
+!--------------------------------------------
+subroutine array_to_matrix(arr,mtrx)
+
+ implicit none
+ real(wp), intent(out) :: mtrx (:,:)
+ real(wp), allocatable, intent(inout) :: arr(:)
+ integer :: i, j, k
+
+ k=1
+ do i=1, size(mtrx,2)
+ do j=1, size(mtrx,1)
+ mtrx(j,i)=arr(k)
+ k=k+1
+ enddo
+ enddo
+
+ deallocate(arr)
+
+end subroutine array_to_matrix
+
+!----------------------------------------------------
+! check if the coord is present -> create origin.coord
+!----------------------------------------------------
+subroutine protectCoord(env)
+
+ use xtb_readin, only : mirror_line
+ implicit none
+ character(len=*),parameter :: source = "xtb_oniom_protectCoord"
+ type(TEnvironment),intent(inout) :: env
+ integer :: cunit, new_cunit, err
+ character(len=:),allocatable :: line
+ logical :: exist
+
+ inquire(file='coord',exist=exist)
+ if(exist) then
+ call open_file(cunit,'coord','r')
+ call open_file(new_cunit,'origin.coord','w')
+ do
+ call mirror_line(cunit,new_cunit,line,err)
+ if(is_iostat_end(err)) exit
+ enddo
+ call env%warning("coord file will be renamed as unopt.coord",source)
+
+ endif
+
+end subroutine protectCoord
end module xtb_oniom
diff --git a/src/optimizer.f90 b/src/optimizer.f90
index 3ebac2ae7..300e0f004 100644
--- a/src/optimizer.f90
+++ b/src/optimizer.f90
@@ -106,6 +106,9 @@ subroutine get_optthr(n,olev,ethr,gthr,maxcycle,acc)
end subroutine get_optthr
+!----------------------------------------
+! Approximate Normal Coordinate Optimizer
+!----------------------------------------
subroutine ancopt(env,ilog,mol,chk,calc, &
& egap,et,maxiter,maxcycle_in,etot,g,sigma,tight,pr,fail)
use xtb_mctc_convert
@@ -130,27 +133,40 @@ subroutine ancopt(env,ilog,mol,chk,calc, &
implicit none
- !> Source of errors in the main program unit
character(len=*), parameter :: source = "optimizer_ancopt"
-
- !> Calculation environment
+ !! source of errors in the main program unit
type(TEnvironment), intent(inout) :: env
-
+ !! calculation environment
type(TMolecule), intent(inout) :: mol
+ !! molecular structure data
integer, intent(in) :: tight
+ !! optimization level
integer, intent(in) :: maxiter
+ !! max SCC cycles
integer, intent(in) :: maxcycle_in
+ !! max GEOOPT cycles
type(TRestart),intent(inout) :: chk
+ !! wrapper for changing info during SCC
class(TCalculator), intent(inout) :: calc
+ !! calculator instance
real(wp) :: eel
+ !! electronic energy
real(wp),intent(inout) :: etot
+ !! total energy
real(wp),intent(in) :: et
+ !! electronic temperature
real(wp),intent(inout) :: egap
+ !! HOMO-LUMO gap
real(wp),intent(inout) :: g(3,mol%n)
+ !! gradients
real(wp),intent(inout) :: sigma(3,3)
+ !! strain derivatives
logical, intent(in) :: pr
+ !! if printed
logical, intent(out) :: fail
-
+ !! if failed
+
+ !> local variables
type(TMolecule) :: molopt
type(scc_results) :: res
type(tb_anc) :: anc
@@ -185,8 +201,10 @@ subroutine ancopt(env,ilog,mol,chk,calc, &
'(10x,":",3x,a,i18, 10x,":")'
character(len=*),parameter :: chrfmt = &
'(10x,":",3x,a,a18, 10x,":")'
+
if(mol%n.eq.1) return
+ !! do not optimize for 1 molecule
if (profile) call timer%new(8,.false.)
if (profile) call timer%measure(1,'optimizer setup')
@@ -282,7 +300,7 @@ subroutine ancopt(env,ilog,mol,chk,calc, &
endif
write(env%unit,scifmt) "energy convergence",ethr, "Eh "
write(env%unit,scifmt) "grad. convergence ",gthr, "Eh/α"
- write(env%unit,dblfmt) "maximium RF displ.",maxdispl," "
+ write(env%unit,dblfmt) "maxmium RF displ. ",maxdispl," "
write(env%unit,scifmt) "Hlow (freq-cutoff)",hlow, " "
write(env%unit,dblfmt) "Hmax (freq-cutoff)",hmax, " "
write(env%unit,dblfmt) "S6 in model hess. ",s6, " "
@@ -563,8 +581,10 @@ subroutine relax(env,iter,mol,anc,restart,maxcycle,maxdispl,ethr,gthr, &
return
endif
if (profile) call timer%measure(6,'optimization log')
+
call writeMolecule(mol, ilog, format=fileType%xyz, energy=res%e_total, &
& gnorm=res%gnorm)
+ !! to write log file
if (profile) call timer%measure(6)
! transform xyz to internal gradient
if (profile) call timer%measure(4)
diff --git a/src/prog/main.F90 b/src/prog/main.F90
index bba163fc2..200e2dd37 100644
--- a/src/prog/main.F90
+++ b/src/prog/main.F90
@@ -90,7 +90,7 @@ module xtb_prog_main
use xtb_kopt
use xtb_iff_iffprepare, only : prepare_IFF
use xtb_iff_data, only : TIFFData
- use xtb_oniom, only : oniom_input, TOniomCalculator
+ use xtb_oniom, only : oniom_input, TOniomCalculator, calculateCharge
use xtb_tblite_calculator, only : TTBLiteCalculator, TTBLiteInput, newTBLiteWavefunction
implicit none
private
@@ -392,6 +392,7 @@ subroutine xtbMain(env, argParser)
endif
mol%chrg = real(set%ichrg, wp)
+ !! To assign charge
mol%uhf = set%nalphabeta
call initrand
@@ -478,6 +479,7 @@ subroutine xtbMain(env, argParser)
! ------------------------------------------------------------------------
!> Print the method header and select the parameter file
+
if (.not.allocated(fnv)) then
select case(set%runtyp)
case default
@@ -590,10 +592,17 @@ subroutine xtbMain(env, argParser)
select type(xtb => calc%real_low)
type is(TxTBCalculator)
call chk%wfn%allocate(mol%n,xtb%basis%nshell,xtb%basis%nao)
+ call newWavefunction(env,mol,xtb,chk)
+ !! assigns only partial charges q and shell charges
if (restart) then ! only in first run
call readRestart(env,chk%wfn,'xtbrestart',mol%n,mol%at,set%gfn_method,exist,.true.)
endif
- end select
+ end select
+ if (.not.set%oniom_settings%fixed_chrgs) then
+ set%oniom_settings%innerchrg = calculateCharge(calc,env,mol,chk)
+ endif
+
+
end select
! ========================================================================
@@ -655,21 +664,32 @@ subroutine xtbMain(env, argParser)
! ------------------------------------------------------------------------
- ! ANCopt
+ !> Geometry optimization(ANCopt,L_ANCopt,FIRE)
+ ! ------------------------------------------------------------------------
if ((set%runtyp.eq.p_run_opt).or.(set%runtyp.eq.p_run_ohess).or. &
& (set%runtyp.eq.p_run_omd).or.(set%runtyp.eq.p_run_screen).or. &
& (set%runtyp.eq.p_run_metaopt)) then
+
if (set%opt_engine.eq.p_engine_rf) &
call ancopt_header(env%unit,set%veryverbose)
+ !! Print ANCopt header
+
call start_timing(3)
+ !! the system_clock and cpu_time calls for the optimization start
+
call geometry_optimization &
& (env, mol,chk,calc, &
& egap,set%etemp,set%maxscciter,set%optset%maxoptcycle,etot,g,sigma,set%optset%optlev,.true.,.false.,murks)
+ !! Optimization
+
+ !> write results
res%e_total = etot
res%gnorm = norm2(g)
if (nscan.gt.0) then
call relaxed_scan(env,mol,chk,calc)
endif
+
+ !> if geo opt fails -> xtblast file
if (murks) then
call generateFileName(tmpname, 'xtblast', extension, mol%ftype)
write(env%unit,'(/,a,1x,a,/)') &
@@ -679,8 +699,10 @@ subroutine xtbMain(env, argParser)
call close_file(ich)
call env%terminate("Geometry optimization failed")
end if
- call stop_timing(3)
- endif
+
+ call stop_timing(3)
+ !! the system_clock and cpu_time calls for the optimization end
+ endif
! ------------------------------------------------------------------------
@@ -1164,6 +1186,7 @@ subroutine parseArguments(env, args, inputFile, paramFile, accuracy, lgrad, &
real(wp) :: ddum
character(len=:), allocatable :: flag, sec
logical :: exist
+
set%gfn_method = 2
coffee = .false.
@@ -1345,10 +1368,9 @@ subroutine parseArguments(env, args, inputFile, paramFile, accuracy, lgrad, &
else
call env%error("Compiled without support for tblite library", source)
cycle
- end if
+ endif
case('--color')
- call args%nextArg(sec)
if (allocated(sec)) then
select case(sec)
case('auto')
@@ -1363,14 +1385,14 @@ subroutine parseArguments(env, args, inputFile, paramFile, accuracy, lgrad, &
else
call env%error("No color scheme provided for --color option", source)
end if
-
+
case('--oniom')
call set_exttyp('oniom')
call args%nextArg(sec)
!> To handle no argument case
if (.not.allocated(sec)) then
- call env%error("No inner region provided for ONIOM", source)
+ call env%error("No inner region is provided for ONIOM", source)
return
end if
call move_alloc(sec, oniom%first_arg)
@@ -1379,14 +1401,17 @@ subroutine parseArguments(env, args, inputFile, paramFile, accuracy, lgrad, &
if (.not.allocated(sec)) then
call env%warning("No method is specified for the ONIOM calculation, default gfn2:gfnff combination will be used", source)
call move_alloc(oniom%first_arg, sec)
- !return
end if
+
inquire(file=sec, exist=exist)
if (exist) then
sec = read_whole_file(sec)
end if
call move_alloc(sec, oniom%second_arg)
+ case('--cut')
+ call set_cut
+
case('--etemp')
call args%nextArg(sec)
if (allocated(sec)) then
diff --git a/src/readin.f90 b/src/readin.f90
index 5d675faab..c4e34f83f 100644
--- a/src/readin.f90
+++ b/src/readin.f90
@@ -16,6 +16,7 @@
! along with xtb. If not, see .
module xtb_readin
+
use xtb_mctc_accuracy, only : wp
use xtb_mctc_strings, only : value
use xtb_type_environment, only : TEnvironment
@@ -34,11 +35,11 @@ module xtb_readin
! this function returns a logical and is always evaluated for its
! side effect (parsing the given string for its real/int/bool value)
interface getValue
- module procedure getIntValue
- module procedure getIntArray
- module procedure getRealValue
- module procedure getRealArray
- module procedure getBoolValue
+ module procedure getIntValue
+ module procedure getIntArray
+ module procedure getRealValue
+ module procedure getRealArray
+ module procedure getBoolValue
end interface getValue
contains
@@ -68,19 +69,23 @@ function xfind(name) result(fname)
end function xfind
!! ------------------------------------------------------------------[SAW]-
-! wrapper around getline from the MCTC lib that strips comments
-! automatically und removes all leading and trailing whitespace
+!> wrapper around getline from the MCTC lib that strips comments
+!> automatically und removes all leading and trailing whitespace
subroutine strip_line(in,line,err)
+
use xtb_mctc_systools, only : getline
implicit none
integer,intent(in) :: in
+ !! input file unit
character(len=:),allocatable,intent(out) :: line
+ !! output line
integer,intent(out) :: err
integer :: ic
call getline(in,line,iostat=err)
if (err.ne.0) return
-! check for comment characters
+
+ !> check for comment characters
ic = index(line,hash)
if (ic.eq.1) then
line = ''
@@ -88,6 +93,7 @@ subroutine strip_line(in,line,err)
else if (ic.gt.1) then
line = line(:ic-1)
endif
+ !> to remove all leading and trailing whitespaces
line = trim(adjustl(line))
end subroutine strip_line
@@ -98,19 +104,25 @@ end subroutine strip_line
! files you plan to replace in the next step of you program.
! Funnily this subroutine exist way before strip_line...
subroutine mirror_line(in,out,line,err)
+
use xtb_mctc_systools, only : getline
implicit none
+
integer,intent(in) :: in
integer,intent(in) :: out
character(len=:),allocatable,intent(out) :: line
integer,intent(out) :: err
+
integer :: ic
call getline(in,line,iostat=err)
+ !! returns the line from input file
if (err.ne.0) return
+ !! iostat=-1 will return err=0
! now write the line to the copy, we write what we read, not what we see
! if (out.ne.-1) write(out,'(a)') trim(line)
-! check for comment characters
+
+ !> Check for comment characters
ic = index(line,hash)
if (ic.eq.1) then
line = ''
diff --git a/src/set_module.f90 b/src/set_module.f90
index 7d47f5405..4e76d3b4b 100644
--- a/src/set_module.f90
+++ b/src/set_module.f90
@@ -708,20 +708,24 @@ subroutine close_set(ictrl)
call close_file(ictrl)
end subroutine close_set
+!> Read xcontrol file
subroutine rdcontrol(fname,env,copy_file)
+
use xtb_readin, only : find_new_name
use xtb_splitparam, only : maxfrag
use xtb_scanparam, only : maxconstr,maxscan
use xtb_sphereparam, only : maxwalls
implicit none
+
character(len=*), parameter :: source = 'set_rdcontrol'
character(len=*),intent(in) :: fname
type(TEnvironment), intent(inout) :: env
+ logical,intent(in),optional :: copy_file
+
character(len=:),allocatable :: line
character(len=:),allocatable :: key
character(len=:),allocatable :: val
character(len=:),allocatable :: newname
- logical,intent(in),optional :: copy_file
integer :: i
integer :: id
integer :: ic
@@ -734,6 +738,7 @@ subroutine rdcontrol(fname,env,copy_file)
logical :: exitRun
if (present(copy_file)) then
+ !! in case of "--copy" command line argument
do_copy=copy_file
else
do_copy=.false.
@@ -747,31 +752,36 @@ subroutine rdcontrol(fname,env,copy_file)
if (do_copy) then
newname = find_new_name(fname)
+ !! newname = #1.fname
! newunit would be nice to have here, but it would always open to a
! negative number, currently I'm checking for a negative number in
! in mirror_line to avoid copying to another unit, so I go with unit 41
call open_file(copy,newname,'w')
else
- copy = -1 ! deactivate copy in mirror_line
+ copy = -1
+ !! deactivate copy in mirror_line
endif
! read first line before the readloop starts, I have to do this
! to avoid using backspace on id (dammit Turbomole format)
call mirror_line(id,copy,line,err)
+ !! to read a line id-unit, and copy it to copy-unit
readflags: do
- ! check if there is a $ in the *first* column
+ !! check if there is a $ in the *first* column
if (index(line,flag).eq.1) then
select case(line(2:))
- ! logical
+ !> logical
case('fit' ); call set_fit; call mirror_line(id,copy,line,err)
+ case('derived' ); call set_derived; call mirror_line(id,copy,line,err)
case('samerand'); call set_samerand; call mirror_line(id,copy,line,err)
case('cma' ); call set_cma; call mirror_line(id,copy,line,err)
- ! data
+ !> data
case('cube' ); call rdblock(env,set_cube, line,id,copy,err,ncount)
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('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)
@@ -787,7 +797,7 @@ subroutine rdcontrol(fname,env,copy_file)
case('path' ); call rdblock(env,set_path, line,id,copy,err,ncount)
case('reactor' ); call rdblock(env,set_reactor, line,id,copy,err,ncount)
case('stm' ); call rdblock(env,set_stm, line,id,copy,err,ncount)
- ! data + user data which is read later, but we start counting here
+ !> data + user data which is read later, but we start counting here
case('fix' ); call rdblock(env,set_fix, line,id,copy,err,ncount)
case('wall' ); call rdblock(env,set_wall, line,id,copy,err,ncount)
maxwalls = maxwalls + ncount
@@ -798,31 +808,35 @@ subroutine rdcontrol(fname,env,copy_file)
maxconstr = maxconstr + ncount
case('split' ); call rdblock(env,set_split, line,id,copy,err,ncount)
maxfrag = maxfrag + ncount
- ! legacy
+ !> legacy
case('set' ); call rdsetbl(env,set_legacy,line,id,copy,err)
- case default ! unknown keyword -> ignore, we don't raise them
- ! except for chrg and spin which you actually can set here
- ! read them here because select case might not catch them that easy
+ case default
+ !! unknown keyword -> ignore, we don't raise them
+ !! except for chrg and spin which you actually can set here
+ !! read them here because select case might not catch them that easy
if (index(line(2:),'chrg').eq.1) call set_chrg(env,line(7:))
if (index(line(2:),'spin').eq.1) call set_spin(env,line(7:))
-! get a new line
call mirror_line(id,copy,line,err)
+ !! get a new line
end select
- else ! not a keyword -> ignore
+ else
+ !! not a keyword -> ignore
call mirror_line(id,copy,line,err)
endif
- ! check for end of file, which I will tolerate as alternative to $end
if (is_iostat_end(err)) exit readflags
+ !! check for end of file, which I will tolerate as alternative to $end
! if (index(line,flag_end).ne.0) exit readflags ! compatibility reasons
call env%check(exitRun)
if (exitRun) then
call env%error("processing of data group failed", source)
exit
end if
+
enddo readflags
if (do_copy) call close_file(copy)
call close_file(id)
+
end subroutine rdcontrol
subroutine rdsetbl(env,handler,line,id,copy,err)
@@ -851,17 +865,21 @@ subroutine rdsetbl(env,handler,line,id,copy,err)
val = trim(adjustl(line(ie+1:)))
call handler(env,key,val)
call env%check(exitRun)
+
if (exitRun) then
call env%error("handler could not process input", source)
return
end if
+
enddo
end subroutine rdsetbl
+!> to read blocks of certain instruction ($)
+!> in other words options of certain flag
subroutine rdblock(env,handler,line,id,copy,err,ncount)
+
implicit none
- character(len=*), parameter :: source = 'set_rdblock'
type(TEnvironment), intent(inout) :: env
procedure(handlerInterface) :: handler
integer,intent(in) :: id
@@ -869,20 +887,26 @@ subroutine rdblock(env,handler,line,id,copy,err,ncount)
integer,intent(out) :: err
integer,intent(out) :: ncount
character(len=:),allocatable,intent(out) :: line
+
+ character(len=*), parameter :: source = 'set_rdblock'
character(len=:),allocatable :: key
character(len=:),allocatable :: val
integer :: ie
logical :: exitRun
+
ncount = 0
do
call mirror_line(id,copy,line,err)
if (is_iostat_end(err)) return
+ !! to check if EOF
if (index(line,flag).ne.0) return
-
- ! find the equal sign
+ !! to check if new flag
ie = index(line,equal)
- if (line.eq.'') cycle ! skip empty lines
- ncount = ncount + 1 ! but count non-empty lines first
+ !! find the equal sign
+ if (line.eq.'') cycle
+ !! skip empty lines
+ ncount = ncount + 1
+ !! but count non-empty lines first
if (ie.eq.0) cycle
key = trim(line(:ie-1))
val = trim(adjustl(line(ie+1:)))
@@ -931,6 +955,7 @@ subroutine set_exttyp(typ)
set%mode_extrun = p_ext_iff
end select
set1 = .false.
+
end subroutine set_exttyp
subroutine set_geopref(typ)
@@ -1018,6 +1043,11 @@ subroutine set_runtyp(typ)
set1 = .false.
end subroutine set_runtyp
+subroutine set_derived
+ implicit none
+ set%oniom_settings%derived = .true.
+end subroutine set_derived
+
subroutine set_fit
implicit none
set%fit = .true.
@@ -1043,22 +1073,63 @@ subroutine set_define
set%define = .true.
end subroutine set_define
+!-----------------------------------
+! Just cut molecule in ONIOM rotuine
+!-----------------------------------
+subroutine set_cut
+ implicit none
+ set%oniom_settings%cut_inner = .true.
+end subroutine set_cut
+
+
+!-----------------------------------
+! Specify charge
+!-----------------------------------
subroutine set_chrg(env,val)
+
implicit none
character(len=*), parameter :: source = 'set_chrg'
+ !! Name of error producer routine
type(TEnvironment), intent(inout) :: env
+ !! Calculation environment to handle I/O stream and error log
character(len=*),intent(in) :: val
+ !! Charge as character
integer :: err
integer :: idum
+ integer :: ind, idum1, idum2
logical,save :: set1 = .true.
+
+
if (set1) then
- if (getValue(env,val,idum)) then
- set%ichrg = idum
+ ind = index(val,":")
+
+
+ if (ind.ne.0) then
+ !! inner:outer
+ if (getValue(env,val(:ind-1),idum1) .and. &
+ & getValue(env,val(ind+1:),idum2)) then
+ set%oniom_settings%fixed_chrgs = .true.
+ set%oniom_settings%innerchrg = idum1
+ set%ichrg = idum2
+ else
+ call env%error('Charge could not be read from your argument',source)
+ endif
+
else
- call env%error('Charge could not be read from your argument',source)
+ !! usual case
+ if (getValue(env,val,idum)) then
+ !! to transform character into int
+ set%ichrg = idum
+ else
+ call env%error('Charge could not be read from your argument',source)
+ endif
+
+
endif
endif
+
set1 = .false.
+
end subroutine set_chrg
subroutine set_spin(env,val)
@@ -1366,6 +1437,42 @@ subroutine set_gfn(env,key,val)
end select
end subroutine set_gfn
+!-----------------------------------
+! set ONIOM functionality
+!-----------------------------------
+subroutine set_oniom(env,key,val)
+ implicit none
+ character(len=*), parameter :: source = 'set_oniom'
+ !! pointer to the error routine
+ type(TEnvironment), intent(inout) :: env
+ !! Calculation environment to handle IO stream and error log
+ character(len=*), intent(in) :: key
+ character(len=*), intent(in) :: val
+ !! key=val
+ logical :: ldum
+ logical, save :: set1 = .true.
+ logical, save :: set2 = .true.
+ logical, save :: set3 = .true.
+
+ select case(key)
+ case default
+ call env%warning("the key '"//key//"' is not recognized by oniom",source)
+ case('inner logs')
+ if (getValue(env,val,ldum).and.set1) set%oniom_settings%logs = .true.
+ set1=.false.
+
+ case('derived k')
+ if (getValue(env,val,ldum).and.set2) set%oniom_settings%derived = .true.
+ set2=.false.
+
+ case('silent')
+ if (getValue(env,val,ldum).and.set2) set%oniom_settings%silent = .true.
+ set3=.false.
+
+ end select
+
+end subroutine set_oniom
+
subroutine set_scc(env,key,val)
implicit none
character(len=*), parameter :: source = 'set_scc'
@@ -2075,12 +2182,15 @@ subroutine set_symmetry(env,key,val)
end select
end subroutine set_symmetry
+!> Options for $external
subroutine set_external(env,key,val)
+
implicit none
- character(len=*), parameter :: source = 'set_external'
type(TEnvironment), intent(inout) :: env
character(len=*),intent(in) :: key
character(len=*),intent(in) :: val
+
+ character(len=*), parameter :: source = 'set_external'
integer :: err
integer :: idum
real(wp) :: ddum
@@ -2093,6 +2203,7 @@ subroutine set_external(env,key,val)
logical,save :: set6 = .true.
logical,save :: set7 = .true.
logical,save :: set8 = .true.
+
select case(key)
case default
call env%warning("the key '"//key//"' is not recognized by external",source)
@@ -2117,7 +2228,11 @@ subroutine set_external(env,key,val)
case('turbodir')
if (set7) set%ext_turbo%path = val
set7 = .false.
+ case ('cefine')
+ if (set8) set%ext_turbo%input_string = val
+ set8 = .false.
end select
+
end subroutine set_external
subroutine set_fix(env,key,val)
diff --git a/src/setparam.f90 b/src/setparam.f90
index 7f6114fd7..6d65f817c 100644
--- a/src/setparam.f90
+++ b/src/setparam.f90
@@ -100,6 +100,23 @@ module xtb_setparam
! interface mode
integer,parameter :: p_pcem_legacy = 1
integer,parameter :: p_pcem_orca = 2
+
+ type oniom_settings
+ integer :: innerchrg
+ !! inner region charge
+ logical :: derived = .false.
+ !! set ONIOM optimization parameter g to derived value
+ logical :: cut_inner = .false.
+ !! to execute xtb just for checking inner region cut
+ logical :: fixed_chrgs= .false.
+ !! if charges for oniom explicitely given
+ logical :: silent = .false.
+ !! zo mute external output
+ logical :: logs = .false.
+ !! if optimization logs of inner regions are needed
+ integer:: ilog1, ilog2
+ !! log units
+ end type oniom_settings
type qm_external
character(len=:),allocatable :: path
@@ -107,6 +124,9 @@ module xtb_setparam
character(len=:),allocatable :: input_file
character(len=:),allocatable :: input_string
logical :: exist
+ !! if input_file exist
+ logical :: oniom=.false.
+ !! special case of the oniom embedding
end type qm_external
integer, parameter :: p_ext_vtb = -1
@@ -391,7 +411,13 @@ module xtb_setparam
real(wp) :: ex_open ! set to 0.5/-0.5 in .xtbrc, respectively
!! ------------------------------------------------------------------------
+! ONIOM
+!! ------------------------------------------------------------------------
+ type(oniom_settings) :: oniom_settings
+!! ------------------------------------------------------------------------
+! External settings
+!! ------------------------------------------------------------------------
type(qm_external) :: ext_driver
type(qm_external) :: ext_orca
type(qm_external) :: ext_turbo
@@ -432,8 +458,11 @@ module xtb_setparam
! character(len=80) :: inputname = ''
character(len= 4) :: pgroup = 'C1 '
-
+!! ------------------------------------------------------------------------
+
end type
+
+
type(TSet) :: set
type(env_setvar) :: xenv
@@ -443,6 +472,7 @@ module xtb_setparam
contains
+
subroutine initrand
implicit none
integer :: i,j
diff --git a/src/type/anc.f90 b/src/type/anc.f90
index 71355e2f6..e50652a85 100644
--- a/src/type/anc.f90
+++ b/src/type/anc.f90
@@ -172,7 +172,7 @@ subroutine generate_anc_blowup(self,iunit,xyz,hess,pr,linear)
! if (abs(self%eigv(i)) > thr1 ) elow = min(elow,self%eigv(i))
! enddo
- damp = max(self%hlow - elow,0.0_wp)
+ damp = max(self%hlow - elow,0.0_wp)
where(abs(self%eigv) > thr2) self%eigv = self%eigv + damp
! do i = 1, self%n3
! if (abs(self%eigv(i)) > thr2 ) self%eigv(i) = self%eigv(i) + damp
diff --git a/src/type/iohandler.f90 b/src/type/iohandler.f90
index a06018f4a..dbf28d2df 100644
--- a/src/type/iohandler.f90
+++ b/src/type/iohandler.f90
@@ -436,6 +436,7 @@ end subroutine closeFile
subroutine deleteFile(self, file, iostat)
+ ! Declaration section
class(TIOHandler), intent(inout) :: self
character(len=*), intent(in) :: file
integer, intent(out), optional :: iostat
@@ -443,18 +444,23 @@ subroutine deleteFile(self, file, iostat)
integer :: unit
logical :: exist
integer :: error
-
+
+ ! Execution section
unit = -1
error = 0
call self%getName(file, name)
+
!$omp critical(io)
+ !! to restrict execution to 1 thread at a time
inquire(file=name, exist=exist)
if (exist) then
open(file=name, newunit=unit, status='old', iostat=error)
+ !! open existing file withh automatically chosen logical unit number
if (error == 0) then
call self%pushBack(TFileHandle(name, fileStatus%deleted, unit, .false.))
+ !! to record deletion in log
close(unit, status='delete', iostat=error)
end if
else
diff --git a/src/type/molecule.f90 b/src/type/molecule.f90
index c3c77902e..b02d09d3d 100644
--- a/src/type/molecule.f90
+++ b/src/type/molecule.f90
@@ -488,9 +488,12 @@ end subroutine deallocate_molecule
subroutine update(self)
use xtb_mctc_accuracy, only : wp
use xtb_pbc_tools
- implicit none
- class(TMolecule),intent(inout) :: self !< molecular structure information
+ implicit none
+ class(TMolecule),intent(inout) :: self
+ !! molecular structure information
+
+ !> For periodic calculations
if (self%npbc > 0) then
call dlat_to_cell(self%lattice,self%cellpar)
call dlat_to_rlat(self%lattice,self%rec_lat)
@@ -498,19 +501,22 @@ subroutine update(self)
call self%wrap_back
endif
-
+
call self%calculate_distances
end subroutine update
!> calculates all distances for molecular structures and minimum
-! image distances for peridic structures
+!> image distances for periodic structures
subroutine mol_calculate_distances(self)
use xtb_mctc_accuracy, only : wp
use xtb_pbc_tools
+
implicit none
- class(TMolecule),intent(inout) :: self !< molecular structure information
+ class(TMolecule),intent(inout) :: self
+ !! molecular structure information
integer :: i,j
+
if (self%npbc > 0) then
do i = 1, self%n
do j = 1, i-1
@@ -530,6 +536,7 @@ subroutine mol_calculate_distances(self)
self%dist(i,i) = 0.0_wp
enddo
endif
+
end subroutine mol_calculate_distances
!> get all nuclear charges
diff --git a/src/type/wavefunction.f90 b/src/type/wavefunction.f90
index 391563817..71712d450 100644
--- a/src/type/wavefunction.f90
+++ b/src/type/wavefunction.f90
@@ -25,23 +25,41 @@ module xtb_type_wavefunction
type :: TWavefunction
integer :: n = 0
+ !! Number of atoms
integer :: nel = 0
+ !! Number of elctrons
integer :: nopen = 0
+ !! Number of unpaired electrons
integer :: nao = 0
+ !! Number of atomic orbitals
integer :: nshell = 0
- real(wp),allocatable :: P(:,:) ! density matrix
- real(wp),allocatable :: q(:) ! partial charges
- real(wp),allocatable :: qsh(:) ! shell charges
- real(wp),allocatable :: dipm(:,:) ! dipole moments
- real(wp),allocatable :: qp(:,:) ! quadrupole moments
- real(wp),allocatable :: wbo(:,:) ! wiberg bond orders
- integer :: ihomo = 0,ihomoa = 0,ihomob = 0 ! HOMO position
+ !! Number of shells
+ real(wp),allocatable :: P(:,:)
+ !! Density matrix
+ real(wp),allocatable :: q(:)
+ !! Partial charges
+ real(wp),allocatable :: qsh(:)
+ !! Shell charges
+ real(wp),allocatable :: dipm(:,:)
+ !! Dipole moments
+ real(wp),allocatable :: qp(:,:)
+ !! Quadrupole moments
+ real(wp),allocatable :: wbo(:,:)
+ !! Wiberg bond orders
+ integer :: ihomo = 0,ihomoa = 0,ihomob = 0
+ !! HOMO position
real(wp) :: efa = 0.0_wp, efb = 0.0_wp
- real(wp),allocatable :: focc(:) ! fractional occupation
- real(wp),allocatable :: focca(:) ! for alpha space
- real(wp),allocatable :: foccb(:) ! for beta space
- real(wp),allocatable :: emo(:) ! orbital energies
- real(wp),allocatable :: C(:,:) ! molecular orbitals
+ real(wp),allocatable :: focc(:)
+ !! Fractional occupation
+ real(wp),allocatable :: focca(:)
+ !! For alpha space
+ real(wp),allocatable :: foccb(:)
+ !! For beta space
+ real(wp),allocatable :: emo(:)
+ !! Orbital energies
+ real(wp),allocatable :: C(:,:)
+ !! Molecular orbitals
+
contains
procedure :: allocate => allocate_wavefunction
procedure :: deallocate => deallocate_wavefunction
diff --git a/src/xtb/calculator.f90 b/src/xtb/calculator.f90
index 799cd9ac5..135c6a489 100644
--- a/src/xtb/calculator.f90
+++ b/src/xtb/calculator.f90
@@ -427,33 +427,58 @@ subroutine writeInfo(self, unit, mol)
end subroutine writeInfo
-
+!---------------------------------------------
+! Initialize new wavefunction
+!---------------------------------------------
subroutine newWavefunction(env, mol, calc, chk)
+
+ implicit none
character(len=*), parameter :: source = 'xtb_calculator_newWavefunction'
+ !! Name of error producer routine
type(TEnvironment), intent(inout) :: env
+ !! Calculation environment to handle I/O stream and error log
type(TRestart), intent(inout) :: chk
+ !! Restart data wrapper for wfn and nlist
type(TxTBCalculator), intent(in) :: calc
+ !! Instance of xTB Calculator
type(TMolecule), intent(in) :: mol
+ !! Molecular structure data
real(wp), allocatable :: cn(:)
+ !! Coordination number
type(chrg_parameter) :: chrgeq
+ !! guess charges(gasteiger/goedecker/sad)
logical :: exitRun
+ !! if it is recommended to terminate the run
associate(wfn => chk%wfn)
allocate(cn(mol%n))
call wfn%allocate(mol%n,calc%basis%nshell,calc%basis%nao)
-
+
+ !> find partial charges
if (mol%npbc > 0) then
+ !! if periodic
wfn%q = mol%chrg/real(mol%n,wp)
+ !! evenly distribute charge with the equal partial charges
else
if (set%guess_charges.eq.p_guess_gasteiger) then
call iniqcn(mol%n,mol%at,mol%z,mol%xyz,nint(mol%chrg),1.0_wp, &
& wfn%q,cn,calc%xtbData%level,.true.)
else if (set%guess_charges.eq.p_guess_goedecker) then
+ !! default
+
call new_charge_model_2019(chrgeq,mol%n,mol%at)
+ !! to get parametrized values for q (en,gam,kappa,alpha)
+
call ncoord_erf(mol%n,mol%at,mol%xyz,cn)
+ !! to obtain CN
+ !! (49) Extended Tight-Binding Quantum Chemistry Mehods 2020
+
call eeq_chrgeq(mol,env,chrgeq,cn,wfn%q)
+ !! to obtain partial charges q
+ !! (47) Extended Tight-Binding Quantum Chemistry Mehods 2020
call env%check(exitRun)
+ !! to check status of environment
if (exitRun) then
call env%rescue("EEQ guess failed, falling back to SAD guess", source)
wfn%q = mol%chrg/real(mol%n,wp)
@@ -462,10 +487,13 @@ subroutine newWavefunction(env, mol, calc, chk)
wfn%q = mol%chrg/real(mol%n,wp)
end if
end if
-
+
+ !> find shell charges
call iniqshell(calc%xtbData,mol%n,mol%at,mol%z,calc%basis%nshell, &
& wfn%q,wfn%qsh,calc%xtbData%level)
+
end associate
+
end subroutine newWavefunction
end module xtb_xtb_calculator
diff --git a/test/unit/CMakeLists.txt b/test/unit/CMakeLists.txt
index 52246c262..a24416929 100644
--- a/test/unit/CMakeLists.txt
+++ b/test/unit/CMakeLists.txt
@@ -34,6 +34,7 @@ set(
"iff"
"latticepoint"
"molecule"
+ "oniom"
"pbc-tools"
"peeq"
"repulsion"
diff --git a/test/unit/main.f90 b/test/unit/main.f90
index 4af45c975..9aaae663c 100644
--- a/test/unit/main.f90
+++ b/test/unit/main.f90
@@ -34,6 +34,7 @@ program tester
use test_iff, only : collect_iff
use test_latticepoint, only : collect_latticepoint
use test_molecule, only : collect_molecule
+ use test_oniom, only : collect_oniom
use test_pbc_tools, only : collect_pbc_tools
use test_peeq, only : collect_peeq
use test_repulsion, only : collect_repulsion
@@ -67,6 +68,7 @@ program tester
new_testsuite("iff", collect_iff), &
new_testsuite("latticepoint", collect_latticepoint), &
new_testsuite("molecule", collect_molecule), &
+ new_testsuite("oniom", collect_oniom), &
new_testsuite("pbc-tools", collect_pbc_tools), &
new_testsuite("peeq", collect_peeq), &
new_testsuite("repulsion", collect_repulsion), &
diff --git a/test/unit/meson.build b/test/unit/meson.build
index 6d9226135..b7244b2b8 100644
--- a/test/unit/meson.build
+++ b/test/unit/meson.build
@@ -47,6 +47,7 @@ tests = [
'iff',
'latticepoint',
'molecule',
+ 'oniom',
'pbc-tools',
'peeq',
'repulsion',
diff --git a/test/unit/test_oniom.f90 b/test/unit/test_oniom.f90
new file mode 100644
index 000000000..db2cdadb0
--- /dev/null
+++ b/test/unit/test_oniom.f90
@@ -0,0 +1,331 @@
+! 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 .
+
+module test_oniom
+ use testdrive, only : new_unittest, unittest_type, error_type, check_ => check, test_failed
+ implicit none
+ private
+
+ public :: collect_oniom
+
+contains
+
+!> Collect all exported unit tests
+subroutine collect_oniom(testsuite)
+ type(unittest_type),allocatable,intent(out) :: testsuite(:)
+ !! array of objects(unit tests)
+
+ testsuite = [ &
+ new_unittest("calculateCharge", test_oniom_calculateCharge), &
+ !! function that registers a new unit test -> result unittest_type object
+ new_unittest("cutbond", test_oniom_cutbond), &
+ new_unittest("singlepoint", test_oniom_singlepoint) &
+ ]
+
+end subroutine collect_oniom
+
+!---------------------------------------------
+! Unit test for automatic charge determination
+!---------------------------------------------
+subroutine test_oniom_calculateCharge(error)
+ use xtb_mctc_accuracy, only : wp
+ use xtb_mctc_io, only : stdout
+ use xtb_mctc_convert, only : aatoau
+
+ use xtb_type_environment
+ use xtb_type_molecule
+ use xtb_type_restart
+ use xtb_type_basisset
+ use xtb_type_environment
+ use xtb_setparam
+
+ use xtb_oniom, only : TOniomCalculator, calculateCharge, newOniomCalculator, oniom_input
+ use xtb_xtb_calculator, only : TxTBCalculator,newWavefunction
+ use xtb_io_writer
+ use xtb_mctc_filetypes, only : filetype
+ type(error_type), allocatable, intent(out) :: error
+ !! error message type, with stat amd message
+ real(wp), parameter :: thr = 1.0e-7_wp
+ !> Molecular structure data
+ integer, parameter :: nat = 17
+ !! 1,3 pentadiene-1-cation + h2o
+ integer, parameter :: at(nat) = [6,6,6,6,6,1,1,1,1,1,1,1,1,1, 8,1,1]
+ real(wp), parameter :: xyz(3,nat) = reshape(&
+ [ 7.11179039031184_wp, -0.65533236978132_wp, 0.64998029924754_wp, &
+ & 9.54574522537207_wp, 0.30133566481849_wp, 0.45943898636196_wp, &
+ & 5.13954887436106_wp, 0.33928833942508_wp, -0.70036877110532_wp, &
+ &11.69126318297630_wp, -0.71972100226778_wp, 1.75001053205592_wp, &
+ & 2.55237710171780_wp, -0.58098675992530_wp, -0.59290041363421_wp, &
+ & 6.76811222294413_wp, -2.23038959587740_wp, 1.90836157103276_wp, &
+ & 9.86172384803294_wp, 1.88061446715386_wp, -0.80779818687663_wp, &
+ & 5.50044858469288_wp, 1.92224580974752_wp, -1.95079140534105_wp, &
+ &11.22134447610320_wp, -1.94366600251515_wp, 3.32992090742500_wp, &
+ &12.66603076756040_wp, -1.91198818602751_wp, 0.30482336978588_wp, &
+ &13.10010470005330_wp, 0.69862396283326_wp, 2.23170047213733_wp, &
+ & 1.27494163620898_wp, 0.97527170197973_wp, -0.13473721665480_wp, &
+ & 1.96940212226880_wp, -1.18150778329094_wp, -2.48288874184834_wp, &
+ & 2.29849092954603_wp, -2.12388855284911_wp, 0.73684528757996_wp, &
+ &14.20884549186960_wp, -3.85859246461132_wp, -2.12437520158188_wp, &
+ &13.93153897438870_wp, -5.65462186833272_wp, -2.15399113818471_wp, &
+ &15.98310331300670_wp, -3.59190426260830_wp, -2.41348990450454_wp &
+ & ], shape(xyz))
+
+ type(TMolecule) :: mol
+ type(TRestart) :: chk1,chk2
+ type(TEnvironment) :: env
+ type(TOniomCalculator) :: calc1,calc2
+ type(oniom_input) :: oniom_input1, oniom_input2
+
+
+ integer :: innchrg, innchrg2
+
+ call init(env)
+ !! construct calculation environment
+
+ call init(mol, at, xyz)
+ !! construct molecular structure type
+ !! interface to initMoleculeNumbers
+ mol%chrg = 1.0_wp
+ !> allocate calculation settings
+ oniom_input1%first_arg = 'turbomole:gfn2'
+ oniom_input2%first_arg = 'orca:gfn2'
+ oniom_input1%second_arg = '1-14'
+ oniom_input2%second_arg = '15-17'
+
+ call newOniomCalculator(calc1,env,mol,oniom_input1)
+ select type(xtb => calc1%real_low)
+ type is(TxTBCalculator)
+ call chk1%wfn%allocate(mol%n,xtb%basis%nshell,xtb%basis%nao)
+ call newWavefunction(env,mol,xtb,chk1)
+ end select
+ innchrg=calculateCharge(calc1,env,mol,chk1)
+ call check_(error,calc1%method_low,2)
+ call check_(error,calc1%method_high,5)
+ call check_(error,chk1%wfn%q(6),0.159091489555113_wp, thr=thr)
+ call check_(error,chk1%wfn%q(17),0.340373819431617_wp, thr=thr)
+ call check_(error,innchrg,1)
+
+ call newOniomCalculator(calc2,env,mol,oniom_input2)
+ select type(xtb => calc2%real_low)
+ type is(TxTBCalculator)
+ call chk2%wfn%allocate(mol%n,xtb%basis%nshell,xtb%basis%nao)
+ call newWavefunction(env,mol,xtb,chk2)
+ end select
+
+ innchrg2=calculateCharge(calc2,env,mol,chk2)
+ call check_(error,calc2%method_low,2)
+ call check_(error,calc2%method_high,4)
+ call check_(error,chk2%wfn%q(6),0.159091489555113_wp, thr=thr)
+ call check_(error,chk2%wfn%q(17),0.340373819431617_wp, thr=thr)
+ call check_(error,innchrg2,0)
+
+endsubroutine test_oniom_calculateCharge
+
+!---------------------------------------------
+! Unit test for cutbond
+!---------------------------------------------
+subroutine test_oniom_cutbond(error)
+ use xtb_mctc_accuracy, only : wp
+ use xtb_mctc_io, only : stdout
+ use xtb_mctc_convert, only : aatoau
+
+ use xtb_type_environment
+ use xtb_type_molecule
+ use xtb_type_restart
+ use xtb_type_basisset
+ use xtb_type_environment
+ use xtb_setparam
+
+ use xtb_oniom, only : TOniomCalculator, calculateCharge, newOniomCalculator, oniom_input
+ use xtb_xtb_calculator, only : TxTBCalculator,newWavefunction
+ use xtb_io_writer
+ use xtb_mctc_filetypes, only : filetype
+ use xtb_type_data, only : scc_results
+
+ type(error_type), allocatable, intent(out) :: error
+ !! error message type, with stat amd message
+ real(wp), parameter :: thr = 1.0e-7_wp
+ !! molecular structure data
+ integer, parameter :: nat = 8
+ !! ethane
+ integer, parameter :: at(nat) = [6,6,1,1,1,1,1,1]
+ real(wp), parameter :: xyz(3,nat) = reshape(&
+ [-6.42056194812604_wp, 4.42017706707601_wp, 0.00000514267503_wp, &
+ &-3.54497465898797_wp, 4.42017905349380_wp, -0.00000041470841_wp, &
+ &-7.14339394244295_wp, 4.96696898911297_wp, 1.84484217131586_wp, &
+ &-7.14339651464075_wp, 2.54910578073934_wp, -0.44887528991923_wp, &
+ &-7.14340046326272_wp, 5.74445490934900_wp, -1.39594753088679_wp, &
+ &-2.82213644603571_wp, 3.09590294303751_wp, 1.39595405551176_wp, &
+ &-2.82214050550548_wp, 6.29125104633740_wp, 0.44887759570039_wp, &
+ &-2.82214295551184_wp, 3.87338472886870_wp, -1.84483683242736_wp &
+ & ], shape(xyz))
+
+ type(TMolecule) :: mol, inner_mol
+ type(TRestart) :: chk
+ type(TEnvironment) :: env
+ type(TOniomCalculator) :: calc
+ type(oniom_input) :: oniom_input3
+ type(scc_results) :: results
+ real(wp),allocatable :: jacobian(:,:)
+ integer,allocatable :: idx_orig(:)
+ real(wp) :: energy, hlgap, sigma(3,3)
+ real(wp), allocatable :: gradient(:,:)
+
+
+ call init(env)
+ !! To initialize environment
+ call init(mol,at,xyz)
+ !! construct molecular structure type
+ !! interface to initMoleculeNumbers
+ mol%chrg = 0.0_wp
+ allocate(gradient(3,mol%n))
+
+ oniom_input3%first_arg = "gfn2:gfn2"
+ oniom_input3%second_arg = '1,3-5'
+
+ call newOniomCalculator(calc,env,mol,oniom_input3)
+
+ select type(xtb => calc%real_low)
+ type is(TxTBCalculator)
+ call chk%wfn%allocate(mol%n,xtb%basis%nshell,xtb%basis%nao)
+ call newWavefunction(env,mol,xtb,chk)
+ end select
+
+ !> Check newOniomCalculator
+ call check_(error,calc%method_low,2)
+ call check_(error,calc%method_high,2)
+ call check_(error,calc%idx(4),5)
+
+ !> check newWavefunction
+ call check_(error,chk%wfn%q(1),-0.258519131601850_wp, thr=thr)
+
+ !> SP for whole molecule
+ call calc%real_low%singlepoint(env,mol,chk,0,.false.,energy,gradient,sigma,hlgap,results)
+
+ call check_(error,energy, -7.336370550119_wp, thr=thr)
+ call check_(error,hlgap, 15.501740190914_wp, thr=thr)
+ 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 check_(error,idx_orig(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
+
+!---------------------------------------------
+! Unit test ONIOM sp
+!---------------------------------------------
+subroutine test_oniom_singlepoint(error)
+ use xtb_mctc_accuracy, only : wp
+ use xtb_mctc_io, only : stdout
+ use xtb_mctc_convert, only : autoaa
+
+ use xtb_type_environment
+ use xtb_type_molecule
+ use xtb_type_restart
+ use xtb_type_basisset
+ use xtb_type_environment
+ use xtb_setparam
+
+ use xtb_oniom, only : TOniomCalculator, calculateCharge, newOniomCalculator, oniom_input
+ use xtb_xtb_calculator, only : TxTBCalculator,newWavefunction
+ use xtb_io_writer
+ use xtb_mctc_filetypes, only : filetype
+ use xtb_type_data, only : scc_results
+ use xtb_gfnff_calculator, only :TGFFCalculator
+
+
+ type(error_type), allocatable, intent(out) :: error
+ !! error message type, with stat amd message
+ real(wp), parameter :: thr = 1.0e-7_wp
+ integer, parameter :: nat = 12
+ !! water cluster
+ integer, parameter :: at(nat) = [8,1,1, 8,1,1, 8,1,1, 8,1,1]
+ real(wp),parameter :: xyz(3,nat) = reshape([&
+ &-2.75237178376284_wp, 2.43247309226225_wp,-0.01392519847964_wp, &
+ &-0.93157260886974_wp, 2.79621404458590_wp,-0.01863384029005_wp, &
+ &-3.43820531288547_wp, 3.30583608421060_wp, 1.42134539425148_wp, &
+ &-2.43247309226225_wp,-2.75237178376284_wp, 0.01392519847964_wp, &
+ &-2.79621404458590_wp,-0.93157260886974_wp, 0.01863384029005_wp, &
+ &-3.30583608421060_wp,-3.43820531288547_wp,-1.42134539425148_wp, &
+ & 2.75237178376284_wp,-2.43247309226225_wp,-0.01392519847964_wp, &
+ & 0.93157260886974_wp,-2.79621404458590_wp,-0.01863384029005_wp, &
+ & 3.43820531288547_wp,-3.30583608421060_wp, 1.42134539425148_wp, &
+ & 2.43247309226225_wp, 2.75237178376284_wp, 0.01392519847964_wp, &
+ & 2.79621404458590_wp, 0.93157260886974_wp, 0.01863384029005_wp, &
+ & 3.30583608421060_wp, 3.43820531288547_wp,-1.42134539425148_wp], shape(xyz))
+
+ type(TMolecule) :: mol, inner_mol
+ type(TRestart) :: chk
+ type(TEnvironment) :: env
+ type(TOniomCalculator) :: calc
+ type(oniom_input) :: oniom_input4
+ type(scc_results) :: results
+ real(wp),allocatable :: jacobian(:,:)
+ integer,allocatable :: idx_orig(:)
+ real(wp) :: energy, hlgap, sigma(3,3)
+ real(wp), allocatable :: gradient(:,:)
+ integer :: io
+
+ call init(env)
+ !! To initialize environment
+ call init(mol,at,xyz)
+ !! construct molecular structure type
+ !! interface to initMoleculeNumbers
+ mol%chrg = 0.0_wp
+ allocate(gradient(3,mol%n))
+ call open_file(io,"w.coord","w")
+ call writeMolecule(mol,io,filetype%tmol)
+ call close_file(io)
+ oniom_input4%second_arg = '1-6'
+
+ call newOniomCalculator(calc,env,mol,oniom_input4)
+
+ select type(xtb => calc%real_low)
+ type is(TGFFCalculator)
+ call check_(error,xtb%topo%qa(1),-0.570632781827558_wp, thr=thr)
+ call check_(error,xtb%topo%qa(2),0.285316390913779_wp, thr=thr)
+ end select
+
+ !> Check newOniomCalculator
+ call check_(error,calc%method_low,3)
+ call check_(error,calc%method_high,2)
+ call check_(error,mol%at(calc%idx(4)),8)
+
+
+ !> SP for whole molecule
+ call calc%singlepoint(env,mol,chk,0,.false.,energy,gradient,sigma,hlgap,results)
+
+ call check_(error,energy, -10.831044042094_wp, thr=thr)
+ call check_(error,norm2(gradient), 0.034482361663_wp, thr=thr)
+ !call check_(error,hlgap, 15.501740190914_wp, thr=thr)
+
+ !> cutbond
+ !call calc%cutbond(env,mol,chk,calc%topo,inner_mol,jacobian,idx_orig)
+
+ !call check_(error,idx_orig(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_singlepoint
+
+
+end module test_oniom