diff --git a/src/calculator/calc_type.f90 b/src/calculator/calc_type.f90
index 77de8a97..3b61bc49 100644
--- a/src/calculator/calc_type.f90
+++ b/src/calculator/calc_type.f90
@@ -97,6 +97,8 @@ module calc_type
character(len=:),allocatable :: shortflag !> shorter job description
!>--- gradient format specifications
+ logical :: numgrad = .false. !> run numerical gradient (expensive!)
+ real(wp) :: gradstep = 0.0005_wp !> displacement for numerical gradient
logical :: rdgrad = .true.
integer :: gradtype = 0
integer :: gradfmt = 0
@@ -1091,6 +1093,9 @@ subroutine calculation_settings_info(self,iunit)
character(len=*),parameter :: fmt3 = '(" :",2x,a20," : ",a)'
character(len=*),parameter :: fmt4 = '(" :",1x,a)'
character(len=20) :: atmp
+ logical :: gxtbwarn
+
+ gxtbwarn=.false.
if (allocated(self%description)) then
write (iunit,'(" :",1x,a)') trim(self%description)
@@ -1111,7 +1116,12 @@ subroutine calculation_settings_info(self,iunit)
end if
if (any((/jobtype%orca,jobtype%xtbsys,jobtype%turbomole, &
& jobtype%generic,jobtype%terachem/) == self%id)) then
- write (iunit,'(" :",3x,a,a)') 'selected binary : ',trim(self%binary)
+ if(index(self%binary,'gxtb').ne.0)then
+ write(iunit,fmt4) 'g-xTB (development version)'
+ gxtbwarn = .true.
+ else
+ write (iunit,'(" :",3x,a,a)') 'selected binary : ',trim(self%binary)
+ endif
end if
if (self%refine_lvl > 0) then
write (atmp,*) 'refinement stage'
@@ -1169,6 +1179,11 @@ subroutine calculation_settings_info(self,iunit)
endif
end if
+ if(gxtbwarn)then
+ write(iunit,fmt4) 'WARNING: This currently is the development version of g-xTB.'
+ write(iunit,fmt4) 'WARNING: Gradients are NUMERICAL (i.e., expensive and noisy!)'
+ endif
+
end subroutine calculation_settings_info
!=========================================================================================!
@@ -1199,10 +1214,15 @@ subroutine create_calclevel_shortcut(self,levelstring)
self%id = jobtype%turbomole
self%rdgrad = .false.
self%binary = 'gp3'
- case ('gxtb')
+ case ('gxtb','gxtb_dev')
self%id = jobtype%turbomole
self%rdgrad = .false.
- self%binary = 'gxtb'
+ self%binary = 'gxtb'
+ self%rdwbo = .false.
+ if(index(levelstring,'_dev').ne.0)then
+ self%other = '-grad'
+ self%rdgrad=.true.
+ endif
case ('orca')
self%id = jobtype%orca
diff --git a/src/calculator/calculator.F90 b/src/calculator/calculator.F90
index 9cd47931..85e02162 100644
--- a/src/calculator/calculator.F90
+++ b/src/calculator/calculator.F90
@@ -179,6 +179,9 @@ subroutine engrad_mol(mol,calc,energy,gradient,iostatus)
!==========================================!
call potential_core(molptr,calc,i,iostatus)
!==========================================!
+ !> and numerical gradient, if selected
+ !==========================================!
+ call numgrad_core(molptr,calc,i,iostatus)
!==========================================!
if (iostatus /= 0) then
@@ -382,6 +385,68 @@ subroutine potential_core(molptr,calc,id,iostatus)
end subroutine potential_core
+ subroutine numgrad_core(molptr,calc,id,iostatus)
+!*******************************************************
+!* subroutine numgrad
+!* routine to perform a numerical gradient calculation
+!*******************************************************
+ implicit none
+ type(coord),intent(in) :: molptr
+ type(calcdata),intent(inout) :: calc
+ integer,intent(in) :: id
+ integer,intent(out) :: iostatus
+
+ integer :: i,j,k,l,ich,och,io,pnat
+ type(coord),allocatable :: moltmp
+ real(wp) :: energy,el,er, step,step2
+ real(wp),allocatable :: ngrd(:,:)
+ !real(wp),parameter :: step = 0.0005_wp
+ !real(wp),parameter :: step2 = 0.5_wp/step
+
+ if (id > calc%ncalculations) return
+ if (.not.calc%calcs(id)%numgrad) return
+
+ pnat = molptr%nat
+ step = calc%calcs(id)%gradstep
+ step2 = 0.5_wp/step
+
+ !> back up energy
+ energy = calc%etmp(id)
+
+ !> allocate temprorary gradient space
+ !$omp critical
+ allocate(ngrd(3,pnat), source=0.0_wp)
+ allocate(moltmp, source=molptr)
+ !$omp end critical
+
+ do i = 1,molptr%nat
+ do j = 1,3
+ moltmp%xyz(j,i) = moltmp%xyz(j,i)+step
+ call potential_core(moltmp,calc,id,iostatus)
+ er = calc%etmp(id)
+
+ moltmp%xyz(j,i) = moltmp%xyz(j,i)-2*step
+ call potential_core(moltmp,calc,id,iostatus)
+ el = calc%etmp(id)
+
+ moltmp%xyz(j,i) = moltmp%xyz(j,i)+step
+ ngrd(j,i) = step2*(er-el)
+ end do
+ end do
+
+ !> transfer tmp gradient to the calc object
+ calc%grdtmp(:,1:pnat,id) = ngrd(:,1:pnat)
+ !$omp critical
+ deallocate(moltmp)
+ deallocate(ngrd)
+ !$omp end critical
+
+ !> restore the energy
+ calc%etmp(id) = energy
+
+ return
+ end subroutine numgrad_core
+
!========================================================================================!
!========================================================================================!
!========================================================================================!
diff --git a/src/calculator/gradreader.f90 b/src/calculator/gradreader.f90
index aa90514d..424c7a04 100644
--- a/src/calculator/gradreader.f90
+++ b/src/calculator/gradreader.f90
@@ -176,7 +176,7 @@ subroutine rd_grad_tm(iunit,nat,energy,grad,iostatus)
integer,intent(out) :: iostatus
integer :: c,io,n,i,j
character(len=128) :: atmp
- character(len=20) :: btmp(8)
+ character(len=20) :: btmp(10)
real(wp) :: dum
logical :: readblock
@@ -184,7 +184,7 @@ subroutine rd_grad_tm(iunit,nat,energy,grad,iostatus)
energy = 0.0_wp
grad(:,:) = 0.0_wp
- c = 0
+ c = 1
readblock = .false.
do
read (iunit,'(a)',iostat=io) atmp
@@ -192,11 +192,14 @@ subroutine rd_grad_tm(iunit,nat,energy,grad,iostatus)
atmp = adjustl(atmp)
if (atmp(1:4) == '$end') readblock = .false.
if( readblock ) then
+
if(index(atmp,'cycle').ne.0)then
- read(atmp,*) btmp(1:2),j,btmp(3:6),energy,btmp(7:8),dum
+ read(atmp,*) btmp(1:8)
+ read(btmp(7),*) energy
elseif(c < nat)then !> skip coords
c = c + 1
else !> read grad
+ !backspace(iunit)
call rd_grad_n3(iunit,nat,grad,iostatus)
exit
endif
@@ -270,7 +273,7 @@ subroutine rd_grad_3n(iunit,nat,grad,iostatus)
grad(:,:) = 0.0_wp
c = 0
- do i = 1,n
+ do i = 1,nat
do j = 1,3
read (iunit,*,iostat=io) dum
if (io < 0) then
@@ -301,7 +304,7 @@ subroutine rd_grad_n3(iunit,nat,grad,iostatus)
grad(:,:) = 0.0_wp
c = 0
- do i = 1,n
+ do i = 1,nat
read (iunit,*,iostat=io) dum(1:3)
if (io < 0) then
iostatus = 3
diff --git a/src/calculator/printouts.F90 b/src/calculator/printouts.F90
index b2658be6..e426c39a 100644
--- a/src/calculator/printouts.F90
+++ b/src/calculator/printouts.F90
@@ -131,7 +131,8 @@ subroutine calculation_summary(calc,mol,energy,grad,molnew,iounit,print)
end if
!>--- gradients
- if (all(calc%calcs(:)%rdgrad.eqv..false.)) then
+ if (all(calc%calcs(:)%rdgrad.eqv..false.) .and. &
+ & all(calc%calcs(:)%numgrad.eqv..false.) ) then
write (iunit,*)
write (iunit,'(a)') '> No gradients calculated'
else if (present(grad)) then
diff --git a/src/confparse.f90 b/src/confparse.f90
index ddb05e29..f06d639d 100644
--- a/src/confparse.f90
+++ b/src/confparse.f90
@@ -94,12 +94,12 @@ subroutine parseflags(env,arg,nra)
!>--- check if help is requested or citations shall be diplayed
do i = 1,nra
if (any((/character(6)::'-h','-H','--h','--H','--help'/) == trim(arg(i)))) then
- if(nra > i)then
- ctmp=trim(arg(i+1))
- if(ctmp(1:1).ne.'-')then
+ if (nra > i) then
+ ctmp = trim(arg(i+1))
+ if (ctmp(1:1) .ne. '-') then
call confscript_morehelp(ctmp)
- endif
- endif
+ end if
+ end if
call confscript_help()
end if
if (any((/character(10)::'-cite','--cite','--citation'/) == trim(arg(i)))) then
@@ -262,7 +262,6 @@ subroutine parseflags(env,arg,nra)
error stop
end if
-
!>--- options for constrained conformer sampling
env%fixfile = 'none selected'
@@ -423,7 +422,6 @@ subroutine parseflags(env,arg,nra)
env%inputcoords = env%ensemblename !> just for a printout
exit
-
case ('-pka','-pKa') !> pKa calculation script
env%crestver = crest_pka
env%runver = 33
@@ -529,7 +527,7 @@ subroutine parseflags(env,arg,nra)
case ('-solvtool','-qcg')
!> Set solute file if present
- if(i == 2) env%solu_file = trim(arg(i-1))
+ if (i == 2) env%solu_file = trim(arg(i-1))
!> Set solvent file if prensent
!> If it is another argument, it doesent matter as solvent file is checke in solvtool
if (nra >= i+1) env%solv_file = trim(arg(i+1))
@@ -555,7 +553,7 @@ subroutine parseflags(env,arg,nra)
env%autozsort = .false.
exit
- case ('-msreact')
+ case ('-msreact')
env%crestver = crest_msreac
env%preopt = .false.
env%presp = .true.
@@ -643,7 +641,7 @@ subroutine parseflags(env,arg,nra)
case ('-rmsd','-rmsdheavy','-hrmsd')
ctmp = trim(arg(i+1))
dtmp = trim(arg(i+2))
- if ((argument == '-rmsdheavy').or.(argument=='-hrmsd')) then
+ if ((argument == '-rmsdheavy').or.(argument == '-hrmsd')) then
call quick_rmsd_tool(ctmp,dtmp,.true.)
else
call quick_rmsd_tool(ctmp,dtmp,.false.)
@@ -715,7 +713,7 @@ subroutine parseflags(env,arg,nra)
env%preopt = .false.
env%crestver = crest_optimize
env%legacy = .false.
- if(argument.eq.'-ohess') env%crest_ohess=.true.
+ if (argument .eq. '-ohess') env%crest_ohess = .true.
exit
case ('-hess','-numhess') !> Numerical hessian
@@ -1039,40 +1037,40 @@ subroutine parseflags(env,arg,nra)
!========================================================================================!
if (env%crestver == crest_msreac) then
select case (argument) !> msreact
- case('-msei')
- env%msei=.true.
- case('-mscid')
- env%mscid=.true.
- env%msei=.false.
- case('-msnoiso') !> filter out non fragmentated structures in msreact
- env%msnoiso=.true.
- case('-msiso') !> filter out fragmentated structures in msreact
- env%msiso=.true.
- case('-msnbonds') ! give number of bonds up to which bias potential is added between atoms default 3
- call readl(arg(i + 1),xx,j)
+ case ('-msei')
+ env%msei = .true.
+ case ('-mscid')
+ env%mscid = .true.
+ env%msei = .false.
+ case ('-msnoiso') !> filter out non fragmentated structures in msreact
+ env%msnoiso = .true.
+ case ('-msiso') !> filter out fragmentated structures in msreact
+ env%msiso = .true.
+ case ('-msnbonds') ! give number of bonds up to which bias potential is added between atoms default 3
+ call readl(arg(i+1),xx,j)
env%msnbonds = xx(1)
- case('-msnshifts') ! give number of times atoms are randomly shifted before optimization
- call readl(arg(i + 1),xx,j)
+ case ('-msnshifts') ! give number of times atoms are randomly shifted before optimization
+ call readl(arg(i+1),xx,j)
env%msnshifts = xx(1)
- case('-msnshifts2') ! give number of times atoms are randomly shifted before applying the constrained optimization default 0
- call readl(arg(i + 1),xx,j)
+ case ('-msnshifts2') ! give number of times atoms are randomly shifted before applying the constrained optimization default 0
+ call readl(arg(i+1),xx,j)
env%msnshifts2 = xx(1)
- case('-msnfrag') ! give number of structures that should be generated
- call readl(arg(i + 1),xx,j)
+ case ('-msnfrag') ! give number of structures that should be generated
+ call readl(arg(i+1),xx,j)
env%msnfrag = xx(1)
- case('-msmolbar') !> filter out structures with same molbar code in msreact
- env%msmolbar=.true.
- case('-msinchi') !> filter out structures with same inchi code in msreact
- env%msinchi=.true.
- case('-msnoattrh') !> add attractive potential for H-atoms
- env%msattrh=.false.
- case('-mslargeprint') !> additional printouts and keep MSDIR
- env%mslargeprint=.true.
- case('-msinput') ! give number of times atoms are randomly shifted before applying the constrained optimization default 0
- ctmp = trim(arg(i+1))
- if (ctmp(1:1) .ne. '-') then
- env%msinput = trim(ctmp)
- end if
+ case ('-msmolbar') !> filter out structures with same molbar code in msreact
+ env%msmolbar = .true.
+ case ('-msinchi') !> filter out structures with same inchi code in msreact
+ env%msinchi = .true.
+ case ('-msnoattrh') !> add attractive potential for H-atoms
+ env%msattrh = .false.
+ case ('-mslargeprint') !> additional printouts and keep MSDIR
+ env%mslargeprint = .true.
+ case ('-msinput') ! give number of times atoms are randomly shifted before applying the constrained optimization default 0
+ ctmp = trim(arg(i+1))
+ if (ctmp(1:1) .ne. '-') then
+ env%msinput = trim(ctmp)
+ end if
end select !> msreact
end if
!========================================================================================!
@@ -1092,7 +1090,7 @@ subroutine parseflags(env,arg,nra)
env%performCross = .true. !> do the genetic crossing
env%autozsort = .true.
case ('-keepdir','-keeptmp') !> Do not delete temporary directories at the end
- env%keepModef = .true.
+ env%keepModef = .true.
case ('-opt','-optlev') !> settings for optimization level of GFN-xTB
env%optlev = optlevnum(arg(i+1))
write (*,'(2x,a,1x,a)') trim(arg(i)),optlevflag(env%optlev)
@@ -1111,7 +1109,7 @@ subroutine parseflags(env,arg,nra)
write (*,'(2x,a,'' : Use of GFN1-xTB requested.'')') env%gfnver
case ('-gfn2')
env%gfnver = '--gfn2'
- write (*,'(2x,a,'' : Use of GFN2-xTB requested.'')') env%gfnver
+ write (*,'(2x,a,'' : Use of GFN2-xTB requested.'')') env%gfnver
case ('-gfn0')
env%gfnver = '--gfn0'
write (*,'(2x,a,'' : Use of GFN0-xTB requested.'')') env%gfnver
@@ -1123,7 +1121,7 @@ subroutine parseflags(env,arg,nra)
ctype = 5 !> bond constraint activated
if (any((/crest_imtd,crest_imtd2/) == env%crestver)) then
bondconst = .true.
- endif
+ end if
env%cts%cbonds_md = .true.
env%checkiso = .true.
case ('stereoisomers')
@@ -1131,6 +1129,12 @@ subroutine parseflags(env,arg,nra)
case default
env%gfnver = '--gfn2'
end select !> GFN
+
+ case ('-gxtb')
+ call gxtb_dev_warning()
+ case ('-gxtb_dev')
+ env%gfnver = 'gxtb_dev'
+
case ('-gfn2@gfn0','-gfn2@gfn1','-gfn2@gff','-gfn2@ff','-gfn2@gfnff')
if (.not.env%legacy) then !TODO
write (*,'("> ",a,1x,a)') argument,'option not yet available with new calculator'
@@ -1173,14 +1177,14 @@ subroutine parseflags(env,arg,nra)
write (*,'(2x,a,a)') argument,' : energy reweighting'
end if
- case('-refine','-rsp','-ropt') !> add one refinement step (via cmd only one is possible)
+ case ('-refine','-rsp','-ropt') !> add one refinement step (via cmd only one is possible)
env%legacy = .false. !> new calculators only!
- if(nra >= i+1)then
+ if (nra >= i+1) then
env%gfnver2 = trim(arg(i+1))
write (*,'(2x,a,1x,a,a)') argument,trim(env%gfnver2), &
& ' : adding refinement step (singlepoint on optimized structures)'
- endif
-
+ end if
+
case ('-charges') !> read charges from file for GFN-FF calcs.
ctmp = trim(arg(i+1))
if ((len_trim(ctmp) < 1).or.(ctmp(1:1) == '-')) then
@@ -1209,8 +1213,8 @@ subroutine parseflags(env,arg,nra)
if (io .eq. 0) env%cts%dscal = rdum
end if
case ('-mtd_kscal','-mtdkscal')
- call readl(arg(i+1),xx,j)
- env%mtd_kscal = xx(1)
+ call readl(arg(i+1),xx,j)
+ env%mtd_kscal = xx(1)
case ('-norestart')
env%allowrestart = .false.
case ('-readbias')
@@ -1436,10 +1440,10 @@ subroutine parseflags(env,arg,nra)
env%potpad = xx(1)
case ('-watoms','-wat')
ctmp = arg(i+1)
- if(ctmp(1:1) .ne. '-')then
- env%potatlist = trim(ctmp)
- write(*,*) env%potatlist
- endif
+ if (ctmp(1:1) .ne. '-') then
+ env%potatlist = trim(ctmp)
+ write (*,*) env%potatlist
+ end if
case ('-wall')
env%wallsetup = .true.
write (*,'(2x,a,1x,a)') '--wall:','requesting setup of wall potential'
@@ -1599,8 +1603,8 @@ subroutine parseflags(env,arg,nra)
env%protb%threshsort = .true.
ctmp = trim(arg(i+1))
if (ctmp(1:1) .ne. '-') then
- read(ctmp,*,iostat=io) idum
- if(io.eq.0) env%protb%amount = idum
+ read (ctmp,*,iostat=io) idum
+ if (io .eq. 0) env%protb%amount = idum
end if
case ('-swel') !> switch out H+ to something else in protonation script
if (env%properties .eq. -3) then
@@ -1612,8 +1616,8 @@ subroutine parseflags(env,arg,nra)
env%protb%threshsort = .true.
ctmp = trim(arg(i+1))
if (ctmp(1:1) .ne. '-') then
- read(ctmp,*,iostat=io) idum
- if(io.eq.0) env%protb%amount = idum
+ read (ctmp,*,iostat=io) idum
+ if (io .eq. 0) env%protb%amount = idum
end if
case ('-tautomerize') !> tautomerization tool
env%properties = p_tautomerize
@@ -1811,7 +1815,7 @@ subroutine parseflags(env,arg,nra)
env%final_gfn2_opt = .false.
case ('-directed') !> specify the directed list
env%qcg_flag = .true.
- ctmp = trim(arg(i + 1))
+ ctmp = trim(arg(i+1))
if (ctmp(1:1) .ne. '-') then
env%directed_file = trim(ctmp)
write (*,'(2x,a,1x,a)') trim(argument)//' :',trim(ctmp)
@@ -2067,7 +2071,7 @@ subroutine parseflags(env,arg,nra)
end if
!>--- automatic wall potential for the LEGACY version
- if (env%NCI.or.env%wallsetup .and. env%legacy) then
+ if (env%NCI.or.env%wallsetup.and.env%legacy) then
call wallpot(env)
if (env%wallsetup) then
write (*,'(2x,a)') 'Automatically generated ellipsoide potential:'
@@ -2139,21 +2143,21 @@ subroutine parseflags(env,arg,nra)
env%lmover = env%gfnver
end if
end if
- if (env%ensemble_opt == '--gfn2' .or. env%gfnver == '--gfn2') &
+ if (env%ensemble_opt == '--gfn2'.or.env%gfnver == '--gfn2') &
& env%final_gfn2_opt = .false. !Prevent additional opt.
if (env%useqmdff) then
env%autozsort = .false.
end if
- if (.not.env%preopt .and. env%crestver.ne.crest_trialopt) then
+ if (.not.env%preopt.and.env%crestver .ne. crest_trialopt) then
if (allocated(env%ref%topo)) deallocate (env%ref%topo)
end if
!>-- turn off niceprint if we are not writing to terminal
- if(env%niceprint)then
+ if (env%niceprint) then
env%niceprint = myisatty(output_unit)
- endif
+ end if
!>-- driver for optimization along trajectory, additional settings
if (.not.any((/crest_mfmdgc,crest_imtd,crest_imtd2,crest_compr/) == env%crestver) &
@@ -2181,22 +2185,22 @@ subroutine parseflags(env,arg,nra)
if (env%sdfformat) then
env%autozsort = .false.
end if
-
+
!>--- 2023/08/19 moved zsort to a standalone property tool
- if(env%autozsort)then
+ if (env%autozsort) then
env%properties = p_zsort
- endif
+ end if
!>--- convert ProgName to absolute path (to make legacy routines more stable)
ctmp = absolute_filepath(trim(env%ProgName))
env%ProgName = ctmp
!>--- for legacy runtypes, check if xtb is present
- if(env%legacy.or.env%QCG)then
+ if (env%legacy.or.env%QCG) then
call checkprog_silent(env%ProgName,.true.,iostat=io)
- if(io /= 0 ) error stop
- write(stdout,'(/,a,a)') 'Selected path to xtb binary: ',trim(env%Progname)
- endif
+ if (io /= 0) error stop
+ write (stdout,'(/,a,a)') 'Selected path to xtb binary: ',trim(env%Progname)
+ end if
!========================================================================================!
!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>!
@@ -2207,25 +2211,25 @@ subroutine parseflags(env,arg,nra)
write (stdout,'(/,a)',advance='no') '> Setting up backup calculator ...'
flush (stdout)
call env2calc_setup(env)
- write(stdout,*) 'done.'
+ write (stdout,*) 'done.'
call env%calc%info(stdout)
end if
!>--- pass on opt-level to new calculator
- if(.not.env%legacy)then
- env%calc%optlev = nint(env%optlev)
- endif
+ if (.not.env%legacy) then
+ env%calc%optlev = nint(env%optlev)
+ end if
!>--- ONIOM setup from toml file
- if (allocated(env%ONIOM_toml))then
- allocate(env%calc%ONIOM)
- call ONIOM_read_toml(env%ONIOM_toml,env%nat,env%ref%at,env%ref%xyz,env%calc%ONIOM)
+ if (allocated(env%ONIOM_toml)) then
+ allocate (env%calc%ONIOM)
+ call ONIOM_read_toml(env%ONIOM_toml,env%nat,env%ref%at,env%ref%xyz,env%calc%ONIOM)
call env%calc%ONIOMexpand()
- endif
+ end if
!>--- important printouts
- if( .not.env%legacy)then
+ if (.not.env%legacy) then
call print_frozen(env)
- endif
+ end if
return
end subroutine parseflags
@@ -2283,7 +2287,7 @@ subroutine parseRC2(env,bondconst)
else
env%cts%used = .false.
return
- end if
+ end if
!>--- read the data
call read_constrainbuffer(env%constraints,env%cts)
@@ -2298,9 +2302,9 @@ subroutine parseRC2(env,bondconst)
end if
end do
end if
- if(.not.env%legacy)then
+ if (.not.env%legacy) then
call parse_xtbinputfile(env,env%constraints)
- endif
+ end if
!>--- some settings
create = .false.
@@ -2450,7 +2454,7 @@ subroutine inputcoords(env,arg)
else
inputfile = 'coord'
end if
- if(.not.allocated(env%inputcoords)) env%inputcoords = inputfile
+ if (.not.allocated(env%inputcoords)) env%inputcoords = inputfile
!>-- if the input was a SDF file, special handling
env%sdfformat = .false.
@@ -2464,16 +2468,16 @@ subroutine inputcoords(env,arg)
if (.not.allocated(env%inputcoords)) env%inputcoords = 'coord'
call mol%open('coord')
!>-- shift to CMA and/or align according to rot.const. We have to be careful about this.
- if (any((/ crest_sp, crest_optimize, crest_numhessian, crest_trialopt /) == env%crestver))then
+ if (any((/crest_sp,crest_optimize,crest_numhessian,crest_trialopt/) == env%crestver)) then
!> some runtypes should only do a CMA translation, but no rotation
call CMAtrf(mol%nat,mol%nat,mol%at,mol%xyz)
- else if (env%crestver == crest_solv)then
+ else if (env%crestver == crest_solv) then
!> runtypes like qcg must not modify input coordinates!
continue
else
!> all other can align with rot. axis
call axis(mol%nat,mol%at,mol%xyz)
- endif
+ end if
!>-- overwrite coord
call mol%write('coord')
diff --git a/src/legacy_wrappers.f90 b/src/legacy_wrappers.f90
index cc4d7e3f..e32cd2cd 100644
--- a/src/legacy_wrappers.f90
+++ b/src/legacy_wrappers.f90
@@ -52,7 +52,8 @@ subroutine env2calc(env,calc,molin)
cal%rdwbo = .false.
cal%rddip = .false.
!> except for SP runtype (from command line!)
- if (env%crestver == crest_sp) then
+ if (env%crestver == crest_sp.and. &
+ & cal%id .ne. jobtype%turbomole) then
cal%rdwbo = .true.
cal%rddip = .true.
cal%rdqat = .true.
@@ -406,13 +407,12 @@ subroutine tautomerize(env,tim)
end if
end subroutine tautomerize
-
!========================================================================================!
subroutine catchdiatomic(env)
!****************************************
!* subroutine catchdiatomic
-!* if we only have one or two atoms just
+!* if we only have one or two atoms just
!* write the "optimized" structure
!****************************************
use crest_data
diff --git a/src/parsing/parse_calcdata.f90 b/src/parsing/parse_calcdata.f90
index a317e5d4..00768671 100644
--- a/src/parsing/parse_calcdata.f90
+++ b/src/parsing/parse_calcdata.f90
@@ -239,8 +239,13 @@ subroutine parse_setting_auto(env,job,kv,rd)
job%id = jobtype%gfn0occ
case ('gfnff','gff','gfn-ff')
job%id = jobtype%gfnff
- case ('pvol','libpvol', 'pv')
+ case ('pvol','libpvol','pv')
job%id = jobtype%libpvol
+ case ('gxtb_dev')
+ job%id = jobtype%turbomole
+ job%rdgrad = .true.
+ job%binary = 'gxtb'
+ job%other ='-grad'
case ('none')
job%id = jobtype%unknown
case ('lj','lennard-jones')
@@ -290,6 +295,11 @@ subroutine parse_setting_auto(env,job,kv,rd)
case ('gradmt')
job%gradfmt = conv2gradfmt(kv%value_c)
+ case ('numgrad')
+ job%numgrad = kv%value_b
+ case ('gradstep')
+ job%gradstep = kv%value_f
+
case ('efile')
job%efile = kv%value_c
@@ -822,7 +832,7 @@ subroutine parse_constraint_auto(env,calc,constr,kv,success,rd)
dum4 = kv%value_fa(6)
call constr%bondrangeconstraint(atm1,atm2,dum1,dum2,beta=dum3,T=dum4)
case default
- write(stdout,'(a)') '**ERROR** wrong number of arguments in bondrange constraint'
+ write (stdout,'(a)') '**ERROR** wrong number of arguments in bondrange constraint'
call creststop(status_config)
end select
success = .true.
diff --git a/src/printouts.f90 b/src/printouts.f90
index 00592cc7..f34ceada 100644
--- a/src/printouts.f90
+++ b/src/printouts.f90
@@ -80,30 +80,45 @@ subroutine box3(version,date,commit,author)
character(len=*) :: date
character(len=*) :: commit
character(len=*) :: author
- character(len=200) :: logo(10)
+ character(len=200) :: logo(13)
character(len=200) :: info(2)
integer,parameter :: pad_left = 7
integer :: i,lcount
write (*,*)
- write (logo(1),'(''╔════════════════════════════════════════════╗'')')
- write (logo(2),'(''║ ___ ___ ___ ___ _____ ║'')')
- write (logo(3),'(''║ / __| _ \ __/ __|_ _| ║'')')
- write (logo(4),'(''║ | (__| / _|\__ \ | | ║'')')
- write (logo(5),'(''║ \___|_|_\___|___/ |_| ║'')')
- write (logo(6),'(''║ ║'')')
- write (logo(7),'(''║ Conformer-Rotamer Ensemble Sampling Tool ║'')')
- write (logo(8),'(''║ based on the xTB methods ║'')')
- write (logo(9),'(''║ ║'')')
- write (logo(10),'("╚════════════════════════════════════════════╝")')
- do i = 1,10
+ !write (logo(1),'(''╔════════════════════════════════════════════╗'')')
+ !write (logo(2),'(''║ ___ ___ ___ ___ _____ ║'')')
+ !write (logo(3),'(''║ / __| _ \ __/ __|_ _| ║'')')
+ !write (logo(4),'(''║ | (__| / _|\__ \ | | ║'')')
+ !write (logo(5),'(''║ \___|_|_\___|___/ |_| ║'')')
+ !write (logo(6),'(''║ ║'')')
+ !write (logo(7),'(''║ Conformer-Rotamer Ensemble Sampling Tool ║'')')
+ !write (logo(8),'(''║ based on the xTB methods ║'')')
+ !write (logo(9),'(''║ ║'')')
+ !write (logo(10),'("╚════════════════════════════════════════════╝")')
+
+ write (logo(1),'(''╔════════════════════════════════════════════════╗'')')
+ write (logo(2),'(''║ ║'')')
+ write (logo(3),'(''║ ██████╗██████╗ ███████╗███████╗████████╗ ║'')')
+ write (logo(4),'(''║ ██╔════╝██╔══██╗██╔════╝██╔════╝╚══██╔══╝ ║'')')
+ write (logo(5),'(''║ ██║ ██████╔╝█████╗ ███████╗ ██║ ║'')')
+ write (logo(6),'(''║ ██║ ██╔══██╗██╔══╝ ╚════██║ ██║ ║'')')
+ write (logo(7),'(''║ ╚██████╗██║ ██║███████╗███████║ ██║ ║'')')
+ write (logo(8),'(''║ ╚═════╝╚═╝ ╚═╝╚══════╝╚══════╝ ╚═╝ ║'')')
+ write (logo(9),'(''║ ║'')')
+ write (logo(10),'(''║ Conformer-Rotamer Ensemble Sampling Tool ║'')')
+ write (logo(11),'(''║ based on the xTB methods ║'')')
+ write (logo(12),'(''║ ║'')')
+ write (logo(13),'(''╚════════════════════════════════════════════════╝'')')
+
+ do i = 1,13
write (*,'(a,a)') repeat(" ",pad_left),trim(logo(i))
end do
- write (*,'(a,''Version '',a,'', '',a)') repeat(" ",pad_left),trim(version),trim(date)
- if(author(1:2).eq."'@")then
- write (*,'(a,"commit (",a,") compiled by ",a)') repeat(" ",pad_left),commit,"'usr"//author(2:)
+ write (*,'(a,'' Version '',a,'', '',a)') repeat(" ",pad_left),trim(version),trim(date)
+ if (author(1:2) .eq. "'@") then
+ write (*,'(a," commit (",a,") compiled by ",a)') repeat(" ",pad_left),commit,"'usr"//author(2:)
else
- write (*,'(a,"commit (",a,") compiled by ",a)') repeat(" ",pad_left),commit,author
- endif
+ write (*,'(a," commit (",a,") compiled by ",a)') repeat(" ",pad_left),commit,author
+ end if
end subroutine box3
subroutine disclaimer
@@ -145,7 +160,7 @@ subroutine confscript_morehelp(flag)
character(len=*),intent(in) :: flag
write (*,'(80("-"))')
- write(*,*)
+ write (*,*)
select case (flag)
case default
write (*,'(/,1x,''General, technical, and calculation options:'')')
@@ -307,22 +322,22 @@ subroutine confscript_morehelp(flag)
write (*,'(5x,''-freqscal : defines frequency scale factor. Only for outprint'')')
write (*,'(5x,''-freqlvl [method] : define a method for frequency computation. All gfn versions are supported'')')
write (*,*)
-
- case('msreact')
+
+ case ('msreact')
write (*,'(1x,'' mass spectral fragment generator (msreact)'')')
write (*,'(1x,''General usage :'')')
write (*,'(5x,'' -msreact [options]'')')
write (*,'(1x,''options:'')')
- write(*,'(5x,''-msnoattrh : deactivate attractive potential between hydrogen and LMO centers)'')')
- write(*,'(5x,''-msnshifts [int] : perform n optimizations with randomly shifted atom postions (default 0) '')')
- write(*,'(5x,''-msnshifts2 [int] : perform n optimizations with randomly shifted atom postions and repulsive potential applied to bonds (default 0) '')')
- write(*,'(5x ''-msnbonds [int] : maximum number of bonds between atoms pairs for applying repulsive potential (default 3)'')')
- write(*,'(5x,''-msmolbar : sort out topological duplicates by molbar codes (requires sourced "molbar")'')')
- write(*,'(5x,''-msinchi : sort out topological duplicates by inchi codes (requires sourced "obabel")'')')
- write(*,'(5x ''-msnfrag [int] : number of fragments that are printed by msreact (random selection)'')')
- write(*,'(5x,''-msiso : print only non-dissociated structures (isomers)'')')
- write(*,'(5x,''-msnoiso : print only dissociated structures'')')
- write(*,'(5x,''-mslargeprint : do not remove temporary files and MSDIR do not remove temporary files and MSDIR with constrained optimizations'')')
+ write (*,'(5x,''-msnoattrh : deactivate attractive potential between hydrogen and LMO centers)'')')
+ write (*,'(5x,''-msnshifts [int] : perform n optimizations with randomly shifted atom postions (default 0) '')')
+ write (*,'(5x,''-msnshifts2 [int] : perform n optimizations with randomly shifted atom postions and repulsive potential applied to bonds (default 0) '')')
+ write (*,'(5x ''-msnbonds [int] : maximum number of bonds between atoms pairs for applying repulsive potential (default 3)'')')
+ write (*,'(5x,''-msmolbar : sort out topological duplicates by molbar codes (requires sourced "molbar")'')')
+ write (*,'(5x,''-msinchi : sort out topological duplicates by inchi codes (requires sourced "obabel")'')')
+ write (*,'(5x ''-msnfrag [int] : number of fragments that are printed by msreact (random selection)'')')
+ write (*,'(5x,''-msiso : print only non-dissociated structures (isomers)'')')
+ write (*,'(5x,''-msnoiso : print only dissociated structures'')')
+ write (*,'(5x,''-mslargeprint : do not remove temporary files and MSDIR do not remove temporary files and MSDIR with constrained optimizations'')')
write (*,'(5x,''-chrg : set the molecules´ charge'')')
write (*,'(5x,''-ewin : set energy window in for sorting out fragments kcal/mol,'')')
write (*,'(5x,'' [default: 200.0 kcal/mol] '')')
@@ -334,7 +349,7 @@ subroutine confscript_morehelp(flag)
write (*,'(5x,'' fc_rep : force constant for repulsive potential between atom pairs (default 0.5) '')')
write (*,'(5x,'' fc_attr : force constant for attractive potential between hydrogen and LMO centers (default -0.5) '')')
write (*,'(5x,'' etemp : electronic temperature in xTB optimizations'')')
-
+
case ('other')
write (*,'(1x,''Other tools for standalone use:'')')
write (*,'(5x,''-zsort : use only the zsort subroutine'')')
@@ -406,8 +421,8 @@ subroutine crestcite
write (*,'( 5x,'' JCTC, 2022, 18 (5), 3174-3189.'')')
write (*,'(/5x,''• P.Pracht, C.Bannwarth, JCTC, 2022, 18 (10), 6370-6385.'')')
write (*,'(/3x,''• P.Pracht, S.Grimme, C.Bannwarth, F.Bohle, S.Ehlert,'')')
- write (*,'( 3x,'' G.Feldmann, J.Gorges, M.Müller, T.Neudecker, C.Plett,'')')
- write (*,'( 3x,'' S.Spicher, P.Steinbach, P.Wesołowski, F.Zeller,'')')
+ write (*,'( 3x,'' G.Feldmann, J.Gorges, M.Müller, T.Neudecker, C.Plett,'')')
+ write (*,'( 3x,'' S.Spicher, P.Steinbach, P.Wesołowski, F.Zeller,'')')
write (*,'( 3x,'' J. Chem. Phys., 2024, 160, 114110.'')')
write (*,'(/,/)')
@@ -545,23 +560,23 @@ end subroutine qcg_head
!========================================================================================!
subroutine msreact_head()
- implicit none
- write (*,*)
- write (*,'(2x,''========================================'')')
- write (*,'(2x,''| |'')')
- write (*,'(2x,''| MSREACT |'')')
- write (*,'(2x,''| automated MS fragment generator |'')')
- write (*,'(2x,''| |'')')
- write (*,'(2x,''| University of Bonn, MCTC |'')')
- write (*,'(2x,''========================================'')')
- write (*,'(2x,'' S. Grimme, P. Pracht, J. Gorges.'')')
- write (*,*)
- write (*,'(3x,''Cite work conducted with this code as'')')
- write (*,'(/,3x,''Philipp Pracht, Stefan Grimme, Christoph Bannwarth, Fabian Bohle, Sebastian Ehlert, Gereon Feldmann,'')')
- write (*,'(3x,''Johannes Gorges, Marcel Müller, Tim Neudecker, Christoph Plett, Sebastian Spicher, Pit Steinbach,'')')
- write (*,'(3x,''Patryk A. Wesolowski, and Felix Zeller J. Chem. Phys., 2024, submitted.'')')
- write (*,*)
- end subroutine msreact_head
+ implicit none
+ write (*,*)
+ write (*,'(2x,''========================================'')')
+ write (*,'(2x,''| |'')')
+ write (*,'(2x,''| MSREACT |'')')
+ write (*,'(2x,''| automated MS fragment generator |'')')
+ write (*,'(2x,''| |'')')
+ write (*,'(2x,''| University of Bonn, MCTC |'')')
+ write (*,'(2x,''========================================'')')
+ write (*,'(2x,'' S. Grimme, P. Pracht, J. Gorges.'')')
+ write (*,*)
+ write (*,'(3x,''Cite work conducted with this code as'')')
+ write (*,'(/,3x,''Philipp Pracht, Stefan Grimme, Christoph Bannwarth, Fabian Bohle, Sebastian Ehlert, Gereon Feldmann,'')')
+ write (*,'(3x,''Johannes Gorges, Marcel Müller, Tim Neudecker, Christoph Plett, Sebastian Spicher, Pit Steinbach,'')')
+ write (*,'(3x,''Patryk A. Wesolowski, and Felix Zeller J. Chem. Phys., 2024, submitted.'')')
+ write (*,*)
+end subroutine msreact_head
!========================================================================================!
@@ -651,17 +666,17 @@ end subroutine mtdwarning
subroutine printiter
implicit none
write (*,*)
- write (*,'(90("*"))')
- write (*,'("**",25x,"N E W I T E R A T I O N C Y C L E",25x,"**")')
- write (*,'(90("*"))')
+ write (*,'(80("*"))')
+ write (*,'("**",20x,"N E W I T E R A T I O N C Y C L E",20x,"**")')
+ write (*,'(80("*"))')
end subroutine printiter
subroutine printiter2(i)
implicit none
integer :: i
write (*,*)
- write (*,'(90("*"))')
- write (*,'("**",26x,"I T E R A T I O N C Y C L E ",i3,23x,"**")') i
- write (*,'(90("*"))')
+ write (*,'(80("*"))')
+ write (*,'("**",21x,"I T E R A T I O N C Y C L E ",i3,18x,"**")') i
+ write (*,'(80("*"))')
end subroutine printiter2
!========================================================================================!
@@ -723,12 +738,12 @@ subroutine print_crest_metadata()
write (*,'(2x,a,1x,a)') 'CREST version :',version
write (*,'(2x,a,1x,a)') 'timestamp :',date
write (*,'(2x,a,1x,a)') 'commit :',commit
- if(author(1:2).eq."'@")then
- l = len_trim(author)
- write (*,'(2x,a,1x,a)') 'compiled by :',"'usr"//author(2:l)
+ if (author(1:2) .eq. "'@") then
+ l = len_trim(author)
+ write (*,'(2x,a,1x,a)') 'compiled by :',"'usr"//author(2:l)
else
- write (*,'(2x,a,1x,a)') 'compiled by :',author
- endif
+ write (*,'(2x,a,1x,a)') 'compiled by :',author
+ end if
write (*,'(2x,a,1x,a)') 'Fortran compiler :',fcompiler
write (*,'(2x,a,1x,a)') 'C compiler :',ccompiler
write (*,'(2x,a,1x,a)') 'build system :',bsystem
@@ -1214,46 +1229,64 @@ end subroutine print_frozen
!========================================================================================!
subroutine progbar(percent,bar)
- use crest_parameters
- implicit none
- real(wp),intent(in) :: percent
- character(len=52),intent(inout) :: bar
- integer :: i
- integer :: done,notdone
-
- bar='['
+ use crest_parameters
+ implicit none
+ real(wp),intent(in) :: percent
+ character(len=52),intent(inout) :: bar
+ integer :: i
+ integer :: done,notdone
- done=nint(percent/2)
- notdone=50-done
+ bar = '['
- do i=1,done
- bar=trim(bar)//'#'
- enddo
+ done = nint(percent/2)
+ notdone = 50-done
+ do i = 1,done
+ bar = trim(bar)//'#'
+ end do
- do i=1,notdone
- bar=trim(bar)//'-'
- enddo
+ do i = 1,notdone
+ bar = trim(bar)//'-'
+ end do
- bar=trim(bar)//']'
+ bar = trim(bar)//']'
end subroutine progbar
subroutine printprogbar(percent)
- use crest_parameters
- implicit none
- real(wp),intent(in) :: percent
- character(len=52) :: bar
+ use crest_parameters
+ implicit none
+ real(wp),intent(in) :: percent
+ character(len=52) :: bar
- if(percent>0.0_wp)then
- call progbar(percent,bar)
- else
- call progbar(0.0_wp,bar)
- endif
- write(0,FMT="(A1,A52,2x,F6.2,A)",ADVANCE="NO") achar(13), &
- & bar, percent, '% finished.'
-
- flush(0)
+ if (percent > 0.0_wp) then
+ call progbar(percent,bar)
+ else
+ call progbar(0.0_wp,bar)
+ end if
+ write (0,FMT="(A1,A52,2x,F6.2,A)",ADVANCE="NO") achar(13), &
+ & bar,percent,'% finished.'
+
+ flush (0)
end subroutine printprogbar
!========================================================================================!
!========================================================================================!
+
+subroutine gxtb_dev_warning
+ use crest_parameters
+ use crest_data, only: status_ioerr
+ write (stdout,*)
+ write (stdout,'(a)') "!!! WARNING !!!"
+ write (stdout,'(a)') "You have selected g-xTB for your calculations, but currently only the"
+ write (stdout,'(a)') "preliminary binary version is available."
+ write (stdout,'(a)') "This version does NOT HAVE ANALYTICAL GRADIENTS available and uses"
+ write (stdout,'(a)') "NUMERICAL gradients which are SLOW and NOISY."
+ write (stdout,*)
+ write (stdout,'(a)') 'The cmd argument "--gxtb" will be disabled until an implementation'
+ write(stdout,'(a)') 'with analytical gradients is available'
+ write(stdout,*)
+ write (stdout,'(a)') 'Please use "--gxtb_dev" in the mean time.'
+ write (stdout,'(a)') "Make sure you have the dev version gxtb installed (https://github.com/grimme-lab/g-xtb)"
+ write(stdout,*)
+ call creststop(status_ioerr)
+end subroutine gxtb_dev_warning