diff --git a/src/axis_trafo.f90 b/src/axis_trafo.f90 index 32dd280ac..4d084e939 100644 --- a/src/axis_trafo.f90 +++ b/src/axis_trafo.f90 @@ -112,96 +112,128 @@ subroutine axis(numat,nat,xyz,aa,bb,cc) return end subroutine axis - subroutine axis2(numat,nat,xyz,aa,bb,cc,avmom,sumw) - use xtb_splitparam - implicit double precision (a-h,o-z) - dimension xyz(3,numat) - integer nat(numat) - PARAMETER (BOHR=0.52917726) - dimension t(6), rot(3), xyzmom(3), eig(3), evec(3,3) - dimension x(numat),y(numat),z(numat),coord(3,numat) - data t /6*0.d0/ - !************************************************************************ - !* const1 = 10**40/(n*a*a) - !* n = avergadro's number - !* a = cm in an angstrom - !* 10**40 is to allow units to be 10**(-40)gram-cm**2 - !* - !************************************************************************ - const1 = 1.66053d0 - !************************************************************************ - !* - !* const2 = conversion factor from angstrom-amu to cm**(-1) - !* - !* = (planck's constant*n*10**16)/(8*pi*pi*c) - !* = 6.62618*10**(-27)[erg-sec]*6.02205*10**23*10**16/ - !* (8*(3.1415926535)**2*2.997925*10**10[cm/sec]) - !* - !************************************************************************ - const2=16.8576522d0 - - sumw=1.d-20 - sumwx=0.d0 - sumwy=0.d0 - sumwz=0.d0 - coord(1:3,1:numat)=xyz(1:3,1:numat)*bohr +subroutine axis2(n,xyz,aa,bb,cc,avmom,sumw) + + use xtb_splitparam, only : atmass + use xtb_mctc_convert, only : autoaa + + implicit double precision (a-h,o-z) + + !> number of atoms + integer, intent(in) :: n + + !> cartesian coordinates + real(wp), intent(in) :: xyz(:,:) + + !> rotational constants + real(wp), intent(out) :: aa,bb,cc + + !> average moment of inertia + real(wp), intent(out) :: avmom + + !> sum of atomic masses + real(wp), intent(out) :: sumw + + !> const1 = 10**40/(n*a*a) + !> n = avergadro's number + !> a = cm in an angstrom + !> 10**40 is to allow units to be 10**(-40)gram-cm**2 + real(wp), parameter :: const1 = 1.66053_wp + + !> const2 = conversion factor from angstrom-amu to cm**(-1) + !> = (planck's constant*n*10**16)/(8*pi*pi*c) + !> = 6.62618*10**(-27)[erg-sec]*6.02205*10**23*10**16/ + !> (8*(3.1415926535)**2*2.997925*10**10[cm/sec]) + real(wp), parameter :: const2 = 16.8576522_wp + + !> temporary variables + real(wp) :: xsum,eps,ams + + !> weighted sum of coordinates + real(wp) :: sumwx,sumwy,sumwz + + !> loop variables + integer :: i + + !> temporary arrays + real(wp) :: t(6), rot(3), xyzmom(3), eig(3), evec(3,3) + real(wp) :: x(n),y(n),z(n) + + !> Cartesian coordinates in Bohr + real(wp) :: coord(3,n) + + t(:) = 0.0_wp + sumw = 0.0_wp + sumwx = 0.0_wp + sumwy = 0.0_wp + sumwz = 0.0_wp + + ! convert Bohr to Angstrom ! + coord(:,:) = xyz(:,:) * autoaa + + ! cma ! + do i=1,n + sumw=sumw+atmass(i) + sumwx=sumwx+atmass(i)*coord(1,i) + sumwy=sumwy+atmass(i)*coord(2,i) + sumwz=sumwz+atmass(i)*coord(3,i) + enddo + + sumwx=sumwx/sumw + sumwy=sumwy/sumw + sumwz=sumwz/sumw + + do i=1,n + x(i)=coord(1,i)-sumwx + y(i)=coord(2,i)-sumwy + z(i)=coord(3,i)-sumwz + enddo + + !************************************************************************ + !* matrix for moments of inertia is of form + !* + !* | y**2+z**2 | + !* | -y*x z**2+x**2 | -i =0 + !* | -z*x -z*y x**2+y**2 | + !* + !************************************************************************ + do i=1,6 + t(i)=dble(i)*1.0d-10 + enddo + + do i=1,n + t(1)=t(1)+atmass(i)*(y(i)**2+z(i)**2) + t(2)=t(2)-atmass(i)*x(i)*y(i) + t(3)=t(3)+atmass(i)*(z(i)**2+x(i)**2) + t(4)=t(4)-atmass(i)*z(i)*x(i) + t(5)=t(5)-atmass(i)*y(i)*z(i) + t(6)=t(6)+atmass(i)*(x(i)**2+y(i)**2) + enddo + + call rsp(t,3,3,eig,evec) + + do i=1,3 + + if(eig(i).lt.3.d-4) then + eig(i)=0.d0 + rot(i)=0.d0 + else + rot(i)=2.9979245d+4*const2/eig(i) + endif - do i=1,numat - sumw=sumw+atmass(i) - sumwx=sumwx+atmass(i)*coord(1,i) - sumwy=sumwy+atmass(i)*coord(2,i) - sumwz=sumwz+atmass(i)*coord(3,i) - enddo + xyzmom(i)=eig(i)*const1 - sumwx=sumwx/sumw - sumwy=sumwy/sumw - sumwz=sumwz/sumw - f=1.0d0/bohr - do i=1,numat - x(i)=coord(1,i)-sumwx - y(i)=coord(2,i)-sumwy - z(i)=coord(3,i)-sumwz - enddo + enddo - !************************************************************************ - !* matrix for moments of inertia is of form - !* - !* | y**2+z**2 | - !* | -y*x z**2+x**2 | -i =0 - !* | -z*x -z*y x**2+y**2 | - !* - !************************************************************************ - do i=1,6 - t(i)=dble(i)*1.0d-10 - enddo - do i=1,numat - t(1)=t(1)+atmass(i)*(y(i)**2+z(i)**2) - t(2)=t(2)-atmass(i)*x(i)*y(i) - t(3)=t(3)+atmass(i)*(z(i)**2+x(i)**2) - t(4)=t(4)-atmass(i)*z(i)*x(i) - t(5)=t(5)-atmass(i)*y(i)*z(i) - t(6)=t(6)+atmass(i)*(x(i)**2+y(i)**2) - enddo - call rsp(t,3,3,eig,evec) - do i=1,3 - if(eig(i).lt.3.d-4) then - eig(i)=0.d0 - rot(i)=0.d0 - else - rot(i)=2.9979245d+4*const2/eig(i) - endif - xyzmom(i)=eig(i)*const1 - enddo + aa=rot(3)/2.9979245d+4 + bb=rot(2)/2.9979245d+4 + cc=rot(1)/2.9979245d+4 - aa=rot(3)/2.9979245d+4 - bb=rot(2)/2.9979245d+4 - cc=rot(1)/2.9979245d+4 - avmom=1.d-47*(xyzmom(1)+xyzmom(2)+xyzmom(3))/3. + avmom=1.d-47*(xyzmom(1)+xyzmom(2)+xyzmom(3))/3. - return - end subroutine axis2 +end subroutine axis2 subroutine axisvec(numat,nat,xyz,aa,bb,cc,evec) use xtb_splitparam diff --git a/src/detrotra.f90 b/src/detrotra.f90 index 0b697ba30..7812d733a 100644 --- a/src/detrotra.f90 +++ b/src/detrotra.f90 @@ -77,62 +77,71 @@ subroutine detrotra4(linear,mol,h,eig) enddo end subroutine detrotra4 - - subroutine detrotra8(linear,n,xyz,h,eig) - use xtb_mctc_accuracy, only : wp - use xtb_type_molecule - implicit none - integer, intent(in) :: n - real(wp), intent(in) :: xyz(3,n) ! values from projected Lindh diag - real(wp), intent(in) :: h(3*n,3*n) ! values from projected Lindh diag - real(wp), intent(inout) :: eig(3*n) ! eigenvectors from projected Lindh diag - logical, intent(in) :: linear - - integer :: i,j,k,kk,ii,nn,n3,nend - integer,allocatable :: ind(:) - real(wp), allocatable :: tmp(:,:) - real(wp), allocatable :: e(:) - real(wp) :: a0,b0,c0 - - n3 = 3*n - - allocate(tmp(3,n),e(n3),ind(n3)) - - nn = 0 - do ii=1, n3 - if(eig(ii).gt.0.05) cycle ! only lowest checked - - kk=0 - do j=1,n - do k=1,3 - kk=kk+1 - tmp(k,j)=xyz(k,j)+h(kk,ii) ! distort along mode ii - enddo - enddo - - c0=0 ! compared all interatomic distances of original and distortet geom. - do i=2,n - do j=1,i-1 - a0 = sqrt((xyz(1,i)-xyz(1,j))**2+(xyz(2,i)-xyz(2,j))**2+(xyz(3,i)-xyz(3,j))**2) - b0 = sqrt((tmp(1,i)-tmp(1,j))**2+(tmp(2,i)-tmp(2,j))**2+(tmp(3,i)-tmp(3,j))**2) - c0 = c0+(a0-b0)**2 - enddo - enddo - nn = nn + 1 - e(nn)=sqrt(c0/n)*abs(eig(ii)) ! weight by Lindh eigenvalue - ind(nn)=nn - - enddo - - call qsort(e, 1, nn, ind) ! sort - - nend = 6 - if (linear) nend = 5 - - do i=1,nend - eig(ind(i)) = 0.0 ! identifier for rot/tra - enddo - - end subroutine detrotra8 + +!> determine rotational and translational modes +subroutine detrotra8(linear,n,xyz,h,eig) + + use xtb_mctc_accuracy, only : wp + use xtb_type_molecule + implicit none + + integer, intent(in) :: n + real(wp), intent(in) :: xyz(3,n) ! values from projected Lindh diag + real(wp), intent(in) :: h(3*n,3*n) ! values from projected Lindh diag + real(wp), intent(inout) :: eig(3*n) ! eigenvalues from projected Lindh diag + logical, intent(in) :: linear + + integer :: i,j,k,kk,ii,nn,n3,nend + integer,allocatable :: ind(:) + real(wp), allocatable :: tmpxyz(:,:) + real(wp), allocatable :: e(:) + real(wp) :: a0,b0,c0 + + n3 = 3*n + + allocate(tmpxyz(3,n),e(n3),ind(n3)) + + nn = 0 + do ii=1, n3 + + if(eig(ii).gt.0.05) cycle ! check only low-lying modes + + ! distort initial geometry along ii-th mode ! + kk=0 + do j=1,n + do k=1,3 + kk=kk+1 + tmpxyz(k,j)=xyz(k,j)+h(kk,ii) + enddo + enddo + + ! compare all interatomic distances of original and distortet geom. ! + c0=0 + do i=2,n + do j=1,i-1 + a0 = sqrt((xyz(1,i)-xyz(1,j))**2+(xyz(2,i)-xyz(2,j))**2+(xyz(3,i)-xyz(3,j))**2) + b0 = sqrt((tmpxyz(1,i)-tmpxyz(1,j))**2+(tmpxyz(2,i)-tmpxyz(2,j))**2+(tmpxyz(3,i)-tmpxyz(3,j))**2) + c0 = c0+(a0-b0)**2 ! sum of squared differences + enddo + enddo + + nn = nn + 1 + e(nn)=sqrt(c0/n)*abs(eig(ii)) ! weight by Lindh eigenvalue + ind(nn)=nn + + enddo + + ! sort in ascending order ! + call qsort(e,1,nn,ind) + + nend = 6 + if (linear) nend = 5 + + ! set lowest eigenvalues to zero ! + do i=1,nend + eig(ind(i)) = 0.0 + enddo + +end subroutine detrotra8 end module xtb_detrotra diff --git a/src/main/property.F90 b/src/main/property.F90 index b4a42abe3..ab22c3717 100644 --- a/src/main/property.F90 +++ b/src/main/property.F90 @@ -1239,7 +1239,7 @@ subroutine print_thermo(iunit, nat, nvib_in, at, xyz, freq, etot, htot, gtot, ni nvib = 0 nimag = 0 - call axis2(nat, at, xyz, aa, bb, cc, avmom, wt) + call axis2(nat,xyz,aa,bb,cc,avmom,wt) nvib_theo = 3 * nat - 6 if (cc < 1.d-10) linear = .true. diff --git a/src/optimizer.f90 b/src/optimizer.f90 index 6c36d7ca5..e6fdc9a4d 100644 --- a/src/optimizer.f90 +++ b/src/optimizer.f90 @@ -22,9 +22,11 @@ module xtb_optimizer use xtb_type_environment, only : TEnvironment use xtb_extern_turbomole, only : TTMCalculator use xtb_bfgs + use xtb_hessian, only : numhess use xtb_david2 implicit none + !> time profiling logical,private,parameter :: profile = .true. type :: convergence_log @@ -42,6 +44,7 @@ module xtb_optimizer contains +!> construct optimizer settings from optimization level subroutine get_optthr(n,olev,ethr,gthr,maxcycle,acc) use xtb_setparam implicit none @@ -51,66 +54,74 @@ subroutine get_optthr(n,olev,ethr,gthr,maxcycle,acc) real(wp),intent(out) :: gthr integer, intent(out) :: maxcycle real(wp),intent(out) :: acc + select case(olev) -! very approximate = crude + ! very approximate = crude ! case(p_olev_crude) ethr = 5.d-4 gthr = 1.d-2 maxcycle=n acc=3.00d0 -! approximate = sloopy + + ! approximate = sloopy ! case(p_olev_sloppy) ethr = 1.d-4 gthr = 6.d-3 maxcycle=n acc=3.00d0 -! loose + + ! loose ! case(p_olev_loose) ethr = 5.d-5 gthr = 4.d-3 maxcycle=n*2 acc=2.00d0 -! for DCOSMO-RS opts with TM i.e. between loose and normal, keyword "lax" + + ! for DCOSMO-RS opts with TM i.e. between loose and normal ! case(p_olev_lax) ethr = 2.d-5 gthr = 2.5d-3 maxcycle=n*2 acc=2.00d0 -! normal + + ! normal ! case default ethr = 5.d-6 gthr = 1.d-3 maxcycle=n*3 acc=1.0d0 -! tight + + ! tight ! case(p_olev_tight) ethr = 1.d-6 gthr = 8.d-4 maxcycle=n*5 acc=0.20d0 -! very tight + + ! very tight ! case(p_olev_vtight) ethr = 1.d-7 gthr = 2.d-4 maxcycle=n*20 acc=0.05d0 -! extreme + + ! extreme ! case(p_olev_extreme) ethr = 5.d-8 gthr = 5.d-5 maxcycle=n*20 acc=0.01d0 end select + maxcycle=min(maxcycle,10000) maxcycle=max(maxcycle,200) end subroutine get_optthr -!---------------------------------------- -! Approximate Normal Coordinate Optimizer -!---------------------------------------- +!> Approximate Normal Coordinate rational function optimizer subroutine ancopt(env,ilog,mol,chk,calc, & & egap,et,maxiter,maxcycle_in,etot,g,sigma,tight,pr,fail) + use xtb_mctc_convert use xtb_mctc_la @@ -132,67 +143,150 @@ subroutine ancopt(env,ilog,mol,chk,calc, & use xtb_lsrmsd implicit none - + + !> traceback for error handling character(len=*), parameter :: source = "optimizer_ancopt" - !! source of errors in the main program unit + + !> calculation environment type(TEnvironment), intent(inout) :: env - !! calculation environment + + !> molecular structure data type(TMolecule), intent(inout) :: mol - !! molecular structure data + + !> optimization level integer, intent(in) :: tight - !! optimization level + + !> max number of SCC cycles integer, intent(in) :: maxiter - !! max SCC cycles + + !> max number of optimization cycles integer, intent(in) :: maxcycle_in - !! max GEOOPT cycles + + !> wrapper for changing info during SCC type(TRestart),intent(inout) :: chk - !! wrapper for changing info during SCC + + !> polymorphic calculator instance class(TCalculator), intent(inout) :: calc - !! calculator instance + + !> electronic energy real(wp) :: eel - !! electronic energy + + !> total energy real(wp),intent(inout) :: etot - !! total energy + + !> electronic temperature real(wp),intent(in) :: et - !! electronic temperature + + !> HOMO-LUMO gap real(wp),intent(inout) :: egap - !! HOMO-LUMO gap + + !> gradients real(wp),intent(inout) :: g(3,mol%n) - !! gradients + + !> strain derivatives real(wp),intent(inout) :: sigma(3,3) - !! strain derivatives + + !> printlevel logical, intent(in) :: pr - !! if printed + + !> optimization failure logical, intent(out) :: fail - !! if failed + +!-----------------! +! local variables ! +!-----------------! - !> local variables type(TMolecule) :: molopt type(scc_results) :: res type(tb_anc) :: anc type(tb_timer) :: timer - real(wp) :: step,amu2au,au2cm,dumi,dumj,damp,hlow,edum,s6,thr,aaa,bbb - real(wp) :: maxdispl,gthr,ethr,hmax,energy,acc,rij(3),t1,t0,w1,w0,ccc - integer :: n3,i,j,k,l,jjj,ic,jc,ia,ja,ii,jj,info,lwork,nat3,liwork - integer :: nvar,iter,nread,maxcycle,maxmicro,itry,maxopt,iupdat,iii + + !> maximum coordinate displacement + real(wp) :: maxdispl + + !> lowest value for force constant + real(wp) :: hlow + + !> highest value for force constant + real(wp) :: hmax + + !> dispersion scaling + real(wp) :: s6 + + !> total energy from initial single point callculation + real(wp) :: estart + + !> energy & gradient convergence thresholds + real(wp) :: ethr, gthr + + + !> number of modes followed + integer :: modef + + !> max number of opt cycles + integer :: maxopt + + !> max number of opt cycles before generation of new ANC + integer :: maxmicro + + !> algorithm for updating Hessian + integer :: iupdat ! 0 = BFGS, 1 = Powell + + !> current iteration + integer :: iter + + !> number of atomic coordinates + integer :: nat3 + + !> deegrees of freedom + integer :: nvar + + !> work array size + integer :: lwork, liwork + + !> step size for numerical Hessian + real(wp) :: step_hess + + real(wp) :: step,amu2au,au2cm,dumi,dumj,damp,edum,thr,aaa,bbb + real(wp) :: energy,acc,rij(3),t1,t0,w1,w0,ccc + + integer :: n3,i,j,k,l,jjj,ic,jc,ia,ja,ii,jj,info + integer :: nread,itry,iii integer :: id,ihess,error integer, intent(in) :: ilog integer, external :: lin - real(wp),allocatable :: h (:,:) + + !> Hessian matrix + real(wp),allocatable :: h (:,:), hess_tmp(:,:) + real(wp),allocatable :: b (:,:) + real(wp),allocatable :: pmode(:,:) + real(wp),allocatable :: grmsd(:,:) + + !> force constants real(wp),allocatable :: fc(:) + + !> eigenvalues real(wp),allocatable :: eig(:) + real(wp),allocatable :: aux(:) real(wp),allocatable :: hess(:) integer, allocatable :: iwork(:) integer, allocatable :: totsym(:) - real(wp),allocatable :: pmode(:,:) - real(wp),allocatable :: grmsd(:,:) type(convergence_log), allocatable :: avconv + real(wp) :: U(3,3), x_center(3), y_center(3), rmsdval - integer :: modef - logical :: restart,ex,converged,linear - real(wp) :: estart,esave + real(wp) :: esave + logical :: restart,ex,converged + + !> if linear molecule + logical :: linear + + !> dipole + real(wp),allocatable :: dipgrad(:,:) + integer, allocatable :: list(:) + + !> formatting strings for output character(len=*),parameter :: scifmt = & '(10x,":",3x,a,e22.7,1x,a,1x,":")' character(len=*),parameter :: dblfmt = & @@ -201,43 +295,51 @@ subroutine ancopt(env,ilog,mol,chk,calc, & '(10x,":",3x,a,i18, 10x,":")' character(len=*),parameter :: chrfmt = & '(10x,":",3x,a,a18, 10x,":")' + + !> error string + character(len=128) :: errStr - ! Print ANCopt header ! + !> debug mode + logical, parameter :: debug(2) = [.false.,.false.] + character(len=9):: hessfmt + + ! print ANCopt header ! call ancopt_header(env%unit,set%veryverbose) - if(mol%n.eq.1) return - !! do not optimize for 1 molecule + if(mol%n.eq.1) return ! skip optimization for 1 atom + + ! performance profiling timer ! if (profile) call timer%new(8,.false.) if (profile) call timer%measure(1,'optimizer setup') - -! defaults - fail =.false. - modef=0 - hmax = 5.0_wp - maxdispl=set%optset%maxdispl_opt - hlow = set%optset%hlow_opt!0.01 in ancopt, 0.002 too small - s6 = set%mhset%s6 !slightly better than 30 for various proteins -! initial number of steps before new ANC are made by model Hessian -! increased during opt. - maxmicro=set%optset%micro_opt + + ! defaults ! + iter = 0 + fail = .false. + modef = 0 + iupdat= 0 + hmax = 5.0_wp + nat3 = 3 * mol%n + maxdispl = set%optset%maxdispl_opt ! def: 1.0 + hlow = set%optset%hlow_opt ! def: 0.01 + s6 = set%mhset%s6 ! def: 20.0 + maxmicro=set%optset%micro_opt ! def: 20 estart = etot - - iupdat=0 !0=BFGS, 1=Powell - + call get_optthr(mol%n,tight,ethr,gthr,maxopt,acc) + + ! if provided by user ! + if(maxcycle_in.gt.0)then + maxopt=maxcycle_in + endif + if(maxopt.lt.maxmicro) maxmicro=maxopt + + ! use Powell update if TS optimization ! if(set%tsopt)then hlow=max(hlow,0.250d0) iupdat=1 endif - - call get_optthr(mol%n,tight,ethr,gthr,maxcycle,acc) - - if(maxcycle_in.le.0)then - maxopt=maxcycle - else - maxopt=maxcycle_in - endif - if(maxopt.lt.maxmicro) maxmicro=maxopt - if (set%optset%average_conv) then + + ! energy and gradient averaging ! + if (set%optset%average_conv) then ! default: .false. select type(calc) class is(TTMCalculator) avconv = load_turbomole_log(maxopt) @@ -250,22 +352,23 @@ subroutine ancopt(env,ilog,mol,chk,calc, & end select end if - call axis2(mol%n,mol%at,mol%xyz,aaa,bbb,ccc,dumi,dumj) - - !call open_file(ilog,'xtbopt.log','w') - iter = 0 - nat3 = 3 * mol%n - nvar = nat3 - 6 - linear = .false. - if(ccc.lt.1.d-10) then + ! determine if linear molecule via rotational constants ! + call axis2(mol%n,mol%xyz,aaa,bbb,ccc,dumi,dumj) + if (ccc.lt.1.d-10) then linear = .true. nvar = nat3 - 5 + else + linear = .false. + nvar = nat3 - 6 endif - if(fixset%n.gt.0) then ! exact fixing - nvar=nat3-3*fixset%n-3 + + ! exact fixing case ! + if(fixset%n.gt.0) then + nvar = nat3 - 3*fixset%n - 3 if(nvar.le.0) nvar=1 endif + ! adjust if restrated ! call open_binary(id,'.xtbtmpmode','r') if(id.ne.-1)then read (id) modef @@ -277,43 +380,52 @@ subroutine ancopt(env,ilog,mol,chk,calc, & allocate(pmode(nat3,1)) ! dummy allocated endif + ! print ANCopt settings ! if(pr)then write(env%unit,'(/,10x,51("."))') write(env%unit,'(10x,":",22x,a,22x,":")') "SETUP" write(env%unit,'(10x,":",49("."),":")') - write(env%unit,chrfmt) "optimization level",int2optlevel(tight) - write(env%unit,intfmt) "max. optcycles ",maxopt - write(env%unit,intfmt) "ANC micro-cycles ",maxmicro - write(env%unit,intfmt) "degrees of freedom",nvar + write(env%unit,chrfmt) "optimization level", int2optlevel(tight) + write(env%unit,intfmt) "max. optcycles ", maxopt + write(env%unit,intfmt) "ANC micro-cycles ", maxmicro + write(env%unit,intfmt) "degrees of freedom", nvar + if (modef>0) then - write(env%unit,intfmt) "# mode follow ",modef + write(env%unit,intfmt) "# mode follow ", modef endif + write(env%unit,'(10x,":",49("."),":")') + if (set%optset%exact_rf) then - write(env%unit,chrfmt) "RF solver ","spevx" + write(env%unit,chrfmt) "RF solver ", "spevx" else - write(env%unit,chrfmt) "RF solver ","davidson" + write(env%unit,chrfmt) "RF solver ", "davidson" endif - write(env%unit,chrfmt) "write xtbopt.log ",bool2string(ilog.ne.-1) + + write(env%unit,chrfmt) "write xtbopt.log ", bool2string(ilog.ne.-1) + if (linear) then - write(env%unit,chrfmt) "linear (good luck)",bool2string(linear) + write(env%unit,chrfmt) "linear (good luck)", bool2string(linear) else - write(env%unit,chrfmt) "linear? ",bool2string(linear) + write(env%unit,chrfmt) "linear? ", bool2string(linear) endif - write(env%unit,scifmt) "energy convergence",ethr, "Eh " - write(env%unit,scifmt) "grad. convergence ",gthr, "Eh/α" - 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, " " + + write(env%unit,scifmt) "energy convergence", ethr, "Eh " + write(env%unit,scifmt) "grad. convergence ", gthr, "Eh/α" + 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, " " write(env%unit,'(10x,51("."))') endif + ! work arrays ! lwork = 1 + 6*nat3 + 2*nat3**2 liwork = 8 * nat3 allocate(h(nat3,nat3),fc(nat3*(nat3+1)/2),eig(nat3)) + ! read in Hessian from file ! if (set%mhset%model == p_modh_read) then call open_file(ihess, 'hessian', 'r') if (ihess == -1) then @@ -332,39 +444,71 @@ subroutine ancopt(env,ilog,mol,chk,calc, & ex = .false. endif - call anc%allocate(mol%n,nvar,hlow,hmax) - - molopt = mol - - if (profile) call timer%measure(1) + call anc%allocate(mol%n,nvar,hlow,hmax) ! allocate ANC + molopt = mol ! copy molecular information + if (profile) call timer%measure(1) ! start opt timer ! ====================================================================== ANC_microiter: do ! ====================================================================== - if (profile) call timer%measure(2,'model hessian') +!----------------------------------------------------------------! +!--------------------- Hessian generation -----------------------! +!----------------------------------------------------------------! + + if (profile) call timer%measure(2,'model hessian') ! start timer for Hessian + + ! initial guess for Hessian ! if (.not.ex)then ! normal case - if(pr)write(env%unit,'(/,''generating ANC from model Hessian ...'')') - call modhes(env, calc, set%mhset, molopt%n, molopt%xyz, molopt%at, fc, pr) ! WBO (array wb) not used in present version - call env%check(fail) - if (fail) then - call env%error("Calculation of model hessian failed", source) - return - end if - !call qpothess(molopt%n,fc,molopt%xyz) - thr=1.d-11 - else - if(pr)write(env%unit,'(/,''generating ANC from read Hessian ...'')') - k=0 - do i=1,nat3 - do j=1,i - k=k+1 - fc(k)=h(j,i) - enddo - enddo - thr=1.d-10 + + if (pr) & + & write(env%unit,'(/,''generating ANC from exact Hessian ...'')') + + call modhes(env, calc, set%mhset, molopt%n, molopt%xyz, molopt%at, fc, pr) ! def: Lindh model (1995) + call env%check(fail) + if (fail) then + call env%error("Calculation of model hessian failed", source) + return + end if + + ! blow up Hessian ! + k=0 + do i=1,nat3 + do j=1,i + k=k+1 + h(i,j)=fc(k) + h(j,i)=fc(k) + enddo + enddo + + thr=1.d-11 + + ! read in Hessian from file ! + else + + if(pr) & + & write(env%unit,'(/,''generating ANC from read Hessian ...'')') + + ! pack Hessian ! + k=0 + do i=1,nat3 + do j=1,i + k=k+1 + fc(k)=h(j,i) + enddo + enddo + + thr=1.d-10 + + endif + + if (debug(2) .and. nat3 <= 30) then !######## DEBUG ######## + write(env%unit,'(/,''Hessian matrix'')') + write(hessfmt,'(a,i0,a)') '(', nat3, 'F10.6)' + write(env%unit,hessfmt) (h(:,i), i=1,nat3) endif + ! project out translational and rotational modes ! if(modef.eq.0)then if(fixset%n.gt.0)then call trproj(molopt%n,nat3,molopt%xyz,fc,.false., -1 ,pmode,1) ! exact fixing @@ -375,42 +519,41 @@ subroutine ancopt(env,ilog,mol,chk,calc, & else call trproj(molopt%n,nat3,molopt%xyz,fc,.false.,modef,pmode,modef) ! NMF endif - if (profile) call timer%measure(2) - - if (profile) call timer%measure(3,'ANC generation') - ! this is completely useless, we blow up the Hessian, just to pack it again... - k=0 - do i=1,nat3 - do j=1,i - k=k+1 - h(i,j)=fc(k) - h(j,i)=fc(k) - enddo - enddo + + if (profile) call timer%measure(2) ! stop timer for model Hessian -! initialize hessian for opt. - call anc%new(env%unit,molopt%xyz,h,pr,linear) +!----------------------------------------------------------------! +!---------------------- ANC generation --------------------------! +!----------------------------------------------------------------! - if (profile) call timer%measure(3) + if (profile) call timer%measure(3,'ANC generation') ! start timer for ANC generation - esave = etot + ! diagonalize Hessian and sort eigenvalues ! + call anc%new(env%unit,molopt%xyz,h,pr,linear) + if (profile) call timer%measure(3) ! stop timer for ANC generation + esave = etot ! save energy -! now everything is prepared for the optimization +!----------------------------------------------------------------! +!---------------------- Optimization ----------------------------! +!----------------------------------------------------------------! call relax(env,iter,molopt,anc,restart,maxmicro,maxdispl,ethr,gthr, & & iii,chk,calc,egap,acc,et,maxiter,iupdat,etot,g,sigma,ilog,pr,fail, & & converged,timer,set%optset%exact_rf,avconv) + ! check for errors ! call env%check(fail) if (fail) then - call env%error("Could not relax structure", source) + call env%error("Could not relax/optimize structure", source) return endif + ! dynamically adjust number of micro iterations ! maxmicro=min(int(maxmicro*1.1),2*set%optset%micro_opt) + ! assess the optimization by RMSD change ! call rmsd(molopt%n,anc%xyz,molopt%xyz,1,U,x_center,y_center,rmsdval,.false.,grmsd) - ! this comes close to a goto, but it's not a goto ... it's even worse + ! this comes close to a goto, but it's not a goto ... it's even worse ! if (restart.and.iter.lt.maxopt) then if (pr) then write(env%unit,'(" * RMSD in coord.:",f14.7,1x,"α")',advance='no') rmsdval @@ -423,6 +566,7 @@ subroutine ancopt(env,ilog,mol,chk,calc, & enddo ANC_microiter ! ====================================================================== + ! convergence report ! if (converged) then if(pr) then call rmsd(mol%n,mol%xyz,molopt%xyz,1,U,x_center,y_center,rmsdval,.false.,grmsd) @@ -440,20 +584,19 @@ subroutine ancopt(env,ilog,mol,chk,calc, & write(env%unit,'(72("-"))') endif else -! not converging in the given cycles is a FAILURE, we should make this clearer -! This is still no ERROR, since we want the geometry written afterwards + ! not converging in the given cycles is a FAILURE, we should make this clearer ! + ! This is still no ERROR, since we want the geometry written afterwards ! if(pr) then write(env%unit,'(/,3x,"***",1x,a,1x,i0,1x,a,1x,"***",/)') & "FAILED TO CONVERGE GEOMETRY OPTIMIZATION IN",iter,"ITERATIONS" endif endif - mol = molopt + mol = molopt ! copy optimized geometry back to molecule - !call close_file(ilog) - - if (pr.and.profile) call timer%write(env%unit,'ANCopt') + if (pr.and.profile) call timer%write(env%unit,'ANCopt') ! write timer report + ! cleanup ! if (profile) call timer%deallocate if (allocated(pmode)) deallocate(pmode) if (allocated(h)) deallocate(h) @@ -485,36 +628,75 @@ subroutine relax(env,iter,mol,anc,restart,maxcycle,maxdispl,ethr,gthr, & implicit none - !> Source of errors in the main program unit + !> source of errors in the main program unit character(len=*), parameter :: source = "optimizer_relax" - !> Calculation environment + !> calculation environment type(TEnvironment), intent(inout) :: env + !> molecular structure data type(TMolecule), intent(inout) :: mol + + !> timer instance type(tb_timer), intent(inout) :: timer + + !> ANC instance type(tb_anc), intent(inout) :: anc + + !> WFN type(TRestart),intent(inout) :: chk + + !> polymorphic calculator instance class(TCalculator), intent(inout) :: calc + + !> number of micro iterations integer, intent(in) :: maxiter + + !> Hessian update algorithm integer, intent(in) :: iupdat + + !> iunit for log file integer, intent(in) :: ilog + + !> temperature real(wp),intent(in) :: et + + !> total energy real(wp),intent(inout) :: etot + + !> SCC accuracy real(wp),intent(in) :: acc_in + + !> gradients real(wp),intent(inout) :: g(3,mol%n) + + !> strain derivatives real(wp),intent(inout) :: sigma(3,3) + + !> HOMO-LUMO gap real(wp),intent(inout) :: egap + + !> conv check logical, intent(out) :: converged + + !> solver type logical, intent(in) :: exact + + !> averaging type(convergence_log), intent(inout), optional :: avconv + !> SCC results storage type(scc_results) :: res + + !> dimensions for RFO integer :: nvar1,npvar,npvar1 + + !> local booleans to control optimization logical :: restart, first, pr, fail logical :: econverged logical :: gconverged logical :: lowered + integer :: maxcycle, iter, prlevel integer :: i, j, ii, jj, k, lwork, info, m, idum, imax(3) real(wp) :: energy,ethr,gthr,dsnrm,maxdispl,t0,w0,t1,w1 @@ -522,56 +704,81 @@ subroutine relax(env,iter,mol,anc,restart,maxcycle,maxdispl,ethr,gthr, & real(wp) :: depred,echng,dummy,maxd,alp,gchng,smallreal,gnold real(wp),allocatable :: gold(:) real(wp),allocatable :: displ(:), gint(:) + + !> RFO eigenvalues real(sp),allocatable :: eaug(:) + + !> RFO eigenvectors real(sp),allocatable :: Uaug(:,:) + + !> Augmented Hessian real(sp),allocatable :: Aaug(:) + real(sp) :: r4dum,sdot parameter (r4dum=1.e-8) parameter (smallreal=1.d-14) - allocate( gold(anc%nvar), displ(anc%nvar), gint(anc%nvar), source = 0.0_wp ) +!----------------------------------------------------------------! +!--------------------- Initialization ---------------------------! +!----------------------------------------------------------------! + - prlevel=0 - if(pr)prlevel=1 - gnorm =0.0_wp - depred =0.0_wp - echng =0.0_wp - maxd =maxdispl - first =.true. - acc =acc_in + ! set printlevel ! + if (pr) then + prlevel = 1 + else + prlevel=0 + endif + + ! initialize variables ! + gnorm = 0.0_wp + depred = 0.0_wp + echng = 0.0_wp + alp = 1.0_wp + maxd = maxdispl + acc = acc_in energy = etot e_in = etot - alp =1.0_wp + first = .true. converged = .false. nvar1 = anc%nvar+1 ! dimension of RF calculation npvar = anc%nvar*(nvar1)/2 ! packed size of Hessian (note the abuse of nvar1!) npvar1 = nvar1*(nvar1+1)/2 ! packed size of augmented Hessian - allocate(Uaug(nvar1,1),eaug(nvar1),Aaug(npvar1)) + + allocate( gold(anc%nvar), displ(anc%nvar), gint(anc%nvar), source = 0.0_wp ) + allocate( Uaug(nvar1,1),eaug(nvar1),Aaug(npvar1) ) !! ======================================================================== main_loop: do ii=1,maxcycle !! ======================================================================== - iter=iter+1 + + iter=iter+1 ! iteration counter if(pr) & - write(env%unit,'(/,72("."),/,30(".")," CYCLE",i5,1x,30("."),/,72("."))')iter + write(env%unit,'(/,72("."),/,30(".")," CYCLE",i5,1x,30("."),/,72("."))')iter + ! save values from previous iteration ! gold = gint gnold= gnorm eold = energy -! calc predicted energy change based on E = E0 + delta * G + delta^2 * H + + ! dE via 2-nd order Taylor expansion ! + ! E = E0 + delta * G + delta^2 * H ! if (ii > 1) & - call prdechng(anc%nvar,gold,displ,anc%hess,depred) -! get gradient + call prdechng(anc%nvar,gold,displ,anc%hess,depred) + + ! displace cartesian coordinates ! if (profile) call timer%measure(4,'coordinate transformation') call anc%get_cartesian(mol%xyz) if (profile) call timer%measure(4) + + ! single point + analytical gradients ! if (profile) call timer%measure(5,'single point calculation') g = 0.0_wp call calc%singlepoint(env,mol,chk,prlevel,iter.eq.1,energy,g,sigma,egap,res) if (profile) call timer%measure(5) - - ! something went wrong in SCC or diag + + ! something went wrong in SCC or diag ! call env%check(fail) if (fail) then call env%error('SCF not converged, aborting...', source) @@ -582,24 +789,28 @@ subroutine relax(env,iter,mol,anc,restart,maxcycle,maxdispl,ethr,gthr, & fail=.true. return endif - if (profile) call timer%measure(6,'optimization log') + ! write geometry to log file ! + 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 + + + ! transform xyz (g) to internal gradient (gint) ! if (profile) call timer%measure(4) call dgemv('t',anc%n3,anc%nvar,1.0_wp,anc%B,anc%n3,g,1,0.0_wp,gint,1) if (profile) call timer%measure(4) - gnorm = norm2(gint) - if(gnorm.gt.500.) then + ! check Euclidean norm of the normal gradient ! + gnorm = norm2(gint) + if(gnorm.gt.500.) then call env%error('|grad| > 500, something is totally wrong!', source) fail=.true. return endif + ! average energy and gradient ! if (present(avconv)) then call avconv%set_eg_log(energy, gnorm) energy = avconv%get_averaged_energy() @@ -612,7 +823,7 @@ subroutine relax(env,iter,mol,anc,restart,maxcycle,maxdispl,ethr,gthr, & end if end if -! adapt SCC acuracy + ! adjust SCC accuracy dynamically ! if(gnorm .lt.0.004)then acc=acc_in elseif(gnorm.lt.0.02)then @@ -621,33 +832,36 @@ subroutine relax(env,iter,mol,anc,restart,maxcycle,maxdispl,ethr,gthr, & acc=6.0d0*acc_in endif - first =.false. + first =.false. ! distinguish between first and subsequent iterations + ! calculate the change in energy and gradient norm ! gchng = gnorm -gnold echng = energy-eold - ! check for convergence + ! check for convergence ! econverged = abs(echng).lt.ethr gconverged = gnorm.lt.gthr lowered = echng.lt.0.0_wp + ! print energy and gradient norm of a current cycle ! if(pr) then - !write(env%unit,'(" E :",F16.8,2x,"G :",F10.6,4x,"pred/act E change:",2D11.3)')& - !energy,gnorm,depred,echng write(env%unit,'(" * total energy :",f14.7,1x,"Eh")',advance='no') energy write(env%unit,'(5x,"change ",e18.7,1x,"Eh")') echng write(env%unit,'(3x,"gradient norm :",f14.7,1x,"Eh/α")',advance='no') gnorm write(env%unit,'(3x,"predicted",e18.7)',advance='no') depred write(env%unit,'(1x,"("f7.2"%)")') (depred-echng)/echng*100 endif + + ! check 0 energy case ! if ( energy .eq. 0 ) then call env%error('external program error', source) return end if - if(ii.eq.1) estart=energy + if(ii.eq.1) estart=energy ! save energy of first iteration - if(gnorm.lt.0.002)then ! 0.002 + ! adjust step size dynamically ! + if(gnorm.lt.0.002)then alp = 1.5d0 ! 1.5 elseif(gnorm.lt.0.0006)then alp = 2.0d0 ! 2 @@ -663,11 +877,10 @@ subroutine relax(env,iter,mol,anc,restart,maxcycle,maxdispl,ethr,gthr, & if (profile) call timer%measure(7,'hessian update') if(ii.gt.1)then -! hessian update if(iupdat.eq.0)then - call bfgs (anc%nvar,gnorm,gint,gold,displ,anc%hess) + call bfgs (anc%nvar,gnorm,gint,gold,displ,anc%hess) else - call powell(anc%nvar,gnorm,gint,gold,displ,anc%hess) + call powell(anc%nvar,gnorm,gint,gold,displ,anc%hess) endif endif if (profile) call timer%measure(7) @@ -682,14 +895,15 @@ subroutine relax(env,iter,mol,anc,restart,maxcycle,maxdispl,ethr,gthr, & ! ⎛ H g ⎞ ⎛ dx ⎞ ⎛ dx ⎞ ! ⎝ g 0 ⎠ ⎝ 1 ⎠ = λ ⎝ 1 ⎠ ! first augment hessian by gradient, we keep everything nicely packed, no blowup - Aaug(1:npvar) = anc%hess - Aaug(npvar+1:npvar1-1) = gint - Aaug(npvar1) = 0.0_sp -! chose your favourite solver + Aaug(1:npvar) = anc%hess ! pack hessian + Aaug(npvar+1:npvar1-1) = gint ! augment with gradient + Aaug(npvar1) = 0.0_sp ! add zero + + ! chose your favourite solver ! if (exact .or. nvar1.lt.50) then call solver_sspevx(nvar1,r4dum,Aaug,Uaug,eaug,fail) else - ! steepest decent guess for displacement + ! steepest decent guess for displacement ! if (ii.eq.1) then Uaug(:,1)=[-real(gint(1:anc%nvar),sp),1.0_sp] dsnrm = sqrt(sdot(nvar1,Uaug,1,Uaug,1)) @@ -699,65 +913,74 @@ subroutine relax(env,iter,mol,anc,restart,maxcycle,maxdispl,ethr,gthr, & if (fail) & ! retry with better solver call solver_sspevx(nvar1,r4dum,Aaug,Uaug,eaug,fail) endif -! divide by last element to get the displacement vector + + ! error handling after RF calculation ! if (fail .or. abs(Uaug(nvar1,1)).lt.1.e-10) then call env%error("internal rational function error", source) return end if + + ! divide by last element to get the displacement vector ! displ(1:anc%nvar) = Uaug(1:anc%nvar,1)/Uaug(nvar1,1) -! check if step is too large, just cut off everything thats to large + + ! check if step is too large, cut off everything that is too large ! do j=1,anc%nvar if(abs(displ(j)).gt.maxd) then if(displ(j) < 0) displ(j)=-maxd if(displ(j) > 0) displ(j)= maxd endif enddo -! now some output - dsnrm=sqrt(ddot(anc%nvar,displ,1,displ,1)) + + ! now some output ! + dsnrm=sqrt(ddot(anc%nvar,displ,1,displ,1)) ! displacement norm if(pr)then - ! this array is currently not used and will be overwritten in next step + + ! find the largest displacement ! gold = abs(displ) imax(1) = maxloc(gold,1); gold(imax(1)) = 0.0_wp imax(2) = maxloc(gold,1); gold(imax(2)) = 0.0_wp imax(3) = maxloc(gold,1) + write(env%unit,'(3x,"displ. norm :",f14.7,1x,"α")',advance='no') & dsnrm*alp - write(env%unit,'(6x,"lambda ",e18.7)') eaug(1) + write(env%unit,'(6x,"lambda ",e18.7)') eaug(1) ! eigenvalue of the RF write(env%unit,'(3x,"maximum displ.:",f14.7,1x,"α")',advance='no') & abs(displ(imax(1)))*alp write(env%unit,'(6x,"in ANC''s ",3("#",i0,", "),"...")') imax - !call prdispl(anc%nvar,displ) endif + if (profile) call timer%measure(8) -! ------------------------------------------------------------------------ -! 2nd: exit and redo hessian (internal restart) + ! 2nd: exit and redo hessian (internal restart) ! if(ii.gt.2.and.dsnrm.gt.2.0) then if (pr) write(*,*) 'exit because of too large step' exit main_loop endif -! new coordinates + ! new coordinates ! anc%coord = anc%coord + displ * alp -! conv ? + ! conv check ! if(abs(echng).lt.ethr.and.gnorm.lt.gthr.and.echng.lt.1.0e-10_wp) then restart=.false. converged = .true. etot=energy return endif - !! ======================================================================== enddo main_loop !! ======================================================================== - if (allocated(Uaug)) deallocate(Uaug) - if (allocated(eaug)) deallocate(eaug) - if (allocated(Aaug)) deallocate(Aaug) + + ! cleanup ! + if (allocated(Uaug)) deallocate(Uaug) + if (allocated(eaug)) deallocate(eaug) + if (allocated(Aaug)) deallocate(Aaug) + ! post micro cycle processing ! restart=.true. etot=energy call anc%get_cartesian(mol%xyz) + end subroutine relax pure subroutine solver_ssyevx(n,thr,A,U,e,fail) @@ -872,6 +1095,7 @@ subroutine sort(nat3,nvar,hess,b) end subroutine sort +!> calculate predicted energy change subroutine prdechng(nat3,grad,displ,hess,depred) !--------------------------------------------------------------------- ! Purpose: @@ -908,11 +1132,11 @@ subroutine prdechng(nat3,grad,displ,hess,depred) !--------------------------------------------------------------------- allocate( hdx(nat3), source = 0.0_wp ) - call blas_spmv('u',nat3,0.5d0,hess,displ,1,0.0d0,hdx,1) + call blas_spmv('u',nat3,0.5d0,hess,displ,1,0.0d0,hdx,1) ! hdx = 0.5*H*displ - gtmp = ddot(nat3,displ,1,grad,1) + gtmp = ddot(nat3,displ,1,grad,1) ! gtmp = grad*displ - htmp = ddot(nat3,displ,1,hdx,1) + htmp = ddot(nat3,displ,1,hdx,1) ! htmp = hdx*displ depred = htmp + gtmp @@ -991,6 +1215,7 @@ subroutine prdispl(nvar,displ) end subroutine prdispl +!> generate model Hessian subroutine modhes(env, calc, modh, natoms, xyz, chg, Hess, pr) use xtb_type_setvar use xtb_modelhessian @@ -1016,21 +1241,22 @@ subroutine modhes(env, calc, modh, natoms, xyz, chg, Hess, pr) type(modhess_setvar),intent(in) :: modh logical, intent(in) :: pr -! Other variables + !> other variables integer :: i integer :: nhess integer, intent(in) :: natoms real(wp),intent(in) :: xyz(3,natoms) - real(wp),intent(out) :: hess((natoms*3)*((natoms*3)+1)/2) + + !> model hessian + real(wp),intent(out) :: Hess((natoms*3)*((natoms*3)+1)/2) integer, intent(in) :: chg(natoms) -! initialize - nhess=3*natoms - Hess=0.d0 + ! initialize ! + Hess=0.0_wp select type(calc) - class default - select case(modh%model) + class default ! all calculator cases except GFN-FF + select case(modh%model) case default call env%error("internal error in model hessian!", source) return @@ -1047,7 +1273,7 @@ subroutine modhes(env, calc, modh, natoms, xyz, chg, Hess, pr) if (pr) write(env%unit,'(a)') "Using Swart-Hessian" call mh_swart(xyz, natoms, Hess, chg, modh) end select - type is(TGFFCalculator) + type is(TGFFCalculator) ! GFN-FF case select case(modh%model) case default call env%error("internal error in model hessian!", source) diff --git a/src/qsort.f90 b/src/qsort.f90 index f6ddccced..5d14a2f0f 100644 --- a/src/qsort.f90 +++ b/src/qsort.f90 @@ -15,32 +15,50 @@ ! You should have received a copy of the GNU Lesser General Public License ! along with xtb. If not, see . +!> QuickSort algorithm recursive subroutine qsort(a, first, last, ind) - use xtb_mctc_accuracy, only : wp - implicit none - real(wp) :: a(*), x, t - integer :: ind(*) - integer :: first, last - integer :: i, j, ii + use xtb_mctc_accuracy, only : wp + implicit none + + !> input array + real(wp) :: a(*) + + !> range to sort + integer :: first, last - x = a( (first+last) / 2 ) - i = first - j = last - do - do while (a(i) < x) + !> index array + integer :: ind(*) + + !> pivot postion + real(wp) :: x + + !> temporary variables + integer :: i, j, ii + real(wp) :: t + + x = a( (first+last) / 2 ) ! choose a middle element as pivot + i = first + j = last + + ! parallel swapping of elements ! + do + do while (a(i) < x) i=i+1 - end do - do while (x < a(j)) + end do + do while (x < a(j)) j=j-1 - end do - if (i >= j) exit - t = a(i); a(i) = a(j); a(j) = t - ii=ind(i); ind(i) = ind(j); ind(j) = ii - i=i+1 - j=j-1 - end do - if (first < i-1) call qsort(a, first, i-1, ind) - if (j+1 < last) call qsort(a, j+1, last, ind) + end do + if (i >= j) exit + t = a(i); a(i) = a(j); a(j) = t + ii=ind(i); ind(i) = ind(j); ind(j) = ii + i=i+1 + j=j-1 + end do + + ! recursive calls ! + if (first < i-1) call qsort(a, first, i-1, ind) + if (j+1 < last) call qsort(a, j+1, last, ind) + end subroutine qsort recursive subroutine qqsort(a, first, last) diff --git a/src/setparam.f90 b/src/setparam.f90 index 06d33fb1c..58ecc0a89 100644 --- a/src/setparam.f90 +++ b/src/setparam.f90 @@ -245,17 +245,16 @@ module xtb_setparam character(len=:), allocatable :: opt_outfile character(len=:), allocatable :: opt_logfile integer, allocatable :: opt_engine + + !> ANCopt settings type(ancopt_setvar) :: optset = ancopt_setvar (& optlev = p_olev_normal, & -! number of opt. cycles before new ANC are made - micro_opt = 20, & -! total number of opt. cycles, 0 means automatically determined + micro_opt = 20, & ! increased during opt. maxoptcycle = 0, & ! det. in ancopt routine if not read in -! maximum coordinate displacement in ancopt maxdispl_opt = 1.000_wp, & -! lowest force constant in ANC generation (should be > 0.005) - hlow_opt = 0.010_wp, & + hlow_opt = 0.010_wp, & ! 0.002 is too small average_conv = .false.) + type(modhess_setvar) :: mhset = modhess_setvar (& model = p_modh_old, & ! force constants for stretch, bend and torsion diff --git a/src/type/anc.f90 b/src/type/anc.f90 index e50652a85..836d70d80 100644 --- a/src/type/anc.f90 +++ b/src/type/anc.f90 @@ -25,22 +25,39 @@ module xtb_type_anc private type tb_anc - integer :: n !< number of atoms - integer :: n3 !< dimension of hessian - integer :: nvar !< actual dimension + integer :: n ! number of atoms + integer :: n3 ! dimension of hessian + integer :: nvar ! actual dimension + + !> lower bound for eigenvalues real(wp) :: hlow + + !> upper bound for eigenvalues real(wp) :: hmax + + !> hessian matrix real(wp),allocatable :: hess(:) + + !> transformation matrix real(wp),allocatable :: B(:,:) + + !> eigenvalues of hessian real(wp),allocatable :: eigv(:) + + !> internal coordinates (approximate normal coordinates) real(wp),allocatable :: coord(:) + + !> cartesian coordinates real(wp),allocatable :: xyz(:,:) contains procedure :: allocate => allocate_anc procedure :: deallocate => deallocate_anc - procedure :: write => write_anc + procedure :: write_anc + procedure :: write_anc_2 + generic :: write => write_anc, write_anc_2 procedure :: new => generate_anc_blowup procedure :: get_cartesian + procedure :: get_normal end type tb_anc contains @@ -81,16 +98,15 @@ subroutine deallocate_anc(self) end subroutine deallocate_anc !> @brief print information about current approximate normal coordinates to unit -subroutine write_anc(self,iunit,comment) +subroutine write_anc(self,iunit) + implicit none class(tb_anc), intent(in) :: self !< approximate normal coordinates integer, intent(in) :: iunit !< file handle - character(len=*),intent(in) :: comment !< name of the variable character(len=*),parameter :: dfmt = '(1x,a,1x,"=",1x,g0)' write(iunit,'(72(">"))') write(iunit,'(1x,"*",1x,a)') "Writing 'tb_anc' class" - write(iunit,'( "->",1x,a)') comment write(iunit,'(72("-"))') write(iunit,'(1x,"*",1x,a)') "status of the fields" write(iunit,dfmt) "integer :: n ",self%n @@ -125,19 +141,69 @@ subroutine write_anc(self,iunit,comment) write(iunit,dfmt) "size(2) :: xyz(:,*) ",size(self%xyz,2) endif write(iunit,'(72("<"))') + end subroutine write_anc +!> print information about current approximate normal coordinates to unit +subroutine write_anc_2(self,iunit,nvar) + + use xtb_mctc_accuracy, only : wp + + implicit none + + !> optimization state holder + class(tb_anc), intent(in) :: self + + !> output unit + integer, intent(in) :: iunit + + !> number of variables to print + integer, intent(in) :: nvar + + integer :: i,j + + write(iunit,*) 'Transformation matrix B' + do i = 1, self%n3 + do j = 1, nvar + write(iunit,'(f8.3,1x)',advance='no') self%B(i,j) + enddo + write(iunit,*) + enddo + + ! print internal coordinates ! + write(iunit,*) 'Internal coordinates' + do i = 1, self%nvar + write(iunit,'(f8.3,1x)',advance='no') self%coord(i) + enddo + write(iunit,*) + +end subroutine write_anc_2 + +!> initialize subroutine generate_anc_blowup(self,iunit,xyz,hess,pr,linear) + use xtb_mctc_accuracy, only : wp use xtb_mctc_la use xtb_detrotra, only : detrotra8 implicit none - class(tb_anc),intent(inout) :: self - integer, intent(in) :: iunit - real(wp), intent(in) :: xyz(3,self%n) - real(wp), intent(inout) :: hess(self%n3,self%n3) - logical, intent(in) :: pr - logical, intent(in) :: linear + + !> ANC object + class(tb_anc),intent(inout) :: self ! anc + + !> output unit + integer, intent(in) :: iunit ! env%unit + + !> cartesian coordinates + real(wp), intent(in) :: xyz(3,self%n) ! molopt%xyz + + !> Hessian matrix + real(wp), intent(inout) :: hess(self%n3,self%n3) ! h + + !> printlevel + logical, intent(in) :: pr ! pr + + !> if linear structure + logical, intent(in) :: linear ! linear real(wp),parameter :: thr1 = 1.0e-10_wp real(wp),parameter :: thr2 = 1.0e-11_wp @@ -145,39 +211,46 @@ subroutine generate_anc_blowup(self,iunit,xyz,hess,pr,linear) integer :: i,itry integer :: nvar integer :: info + + !> LAPACK, length of auxilary fp workspace integer :: lwork + + !> LAPACK, length of auxilary integer workspace integer :: liwork - logical :: fail + + !> LAPACK, auxilary workspace for integer operations integer, allocatable :: iwork(:) + real(wp) :: elow,damp,thr - real(wp),allocatable :: aux(:) + + !> LAPACK, auxilary workspace for floating point operations + real(wp),allocatable :: aux(:) ! = work + logical :: fail + + ! Intialize ! self%xyz = xyz - thr = thr2 lwork = 1 + 6*self%n3 + 2*self%n3**2 liwork = 8 * self%n3 + allocate(iwork(liwork), source = 0) + allocate(aux(lwork), source = 0.0_wp) - allocate(iwork(liwork), source = 0 ) - allocate(aux(lwork), source = 0.0_wp ) - + ! Diagonalize Hessian ! call lapack_syevd('V','U',self%n3,hess,self%n3,self%eigv, & & aux,lwork,iwork,liwork,info) + ! determine, sort, and nullify rot/trans modes ! call detrotra8(linear,self%n,self%xyz,hess,self%eigv) - !elow = 1.0e+99_wp - elow = minval(self%eigv,mask=(abs(self%eigv) > thr1)) -! do i = 1, self%n3 -! if (abs(self%eigv(i)) > thr1 ) elow = min(elow,self%eigv(i)) -! enddo + ! find lowest eigenvalue (ignore nullified ones) ! + elow = minval(self%eigv,mask=(abs(self%eigv) > thr1)) + ! shift eigenvalues to hlow ! 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 -! enddo + ! print eigenvalue spectrum ! if(pr)then write(iunit,*) 'Shifting diagonal of input Hessian by ', damp write(iunit,*) 'Lowest eigenvalues of input Hessian' @@ -188,11 +261,13 @@ subroutine generate_anc_blowup(self,iunit,xyz,hess,pr,linear) endif fail = .true. - get_anc: do itry = 1, maxtry + get_anc: do itry = 1, maxtry ! 4 times + self%B = 0.0_wp self%hess = 0.0_wp nvar = 0 - ! take largest (positive) first + + ! take largest (positive) first ! do i = self%n3, 1, -1 if (abs(self%eigv(i)) > thr .and. nvar < self%nvar) then nvar = nvar+1 @@ -202,6 +277,7 @@ subroutine generate_anc_blowup(self,iunit,xyz,hess,pr,linear) endif enddo + ! reduce thr if not enough modes found ! if (nvar.ne.self%nvar) then thr = thr * 0.1_wp cycle get_anc @@ -212,14 +288,15 @@ subroutine generate_anc_blowup(self,iunit,xyz,hess,pr,linear) enddo get_anc - if (fail) then - write(*,*) 'nvar, selv%nvar',nvar,self%nvar - call raise('E',"ANC generation failed!") - end if + ! nvar .ne. self%nvar ! + if (fail) & + & call raise('E',"ANC generation failed, not enough modes!") + ! sort eigenvalues and eigenvectors in ascending order ! call sort(self%n3,self%nvar,self%hess,self%B) self%coord = 0.0_wp + end subroutine generate_anc_blowup subroutine generate_anc_packed(self,xyz,hess,pr) @@ -285,7 +362,7 @@ subroutine generate_anc_packed(self,xyz,hess,pr) self%B = 0.0_wp self%hess = 0.0_wp nvar = 0 - ! take largest (positive) first + ! take largest (positive) first ! do i = self%n3, 1, -1 if (abs(self%eigv(i)) > thr .and. nvar < self%nvar) then nvar = nvar+1 @@ -311,39 +388,52 @@ subroutine generate_anc_packed(self,xyz,hess,pr) call sort(self%n3,self%nvar,self%hess,self%B) self%coord = 0.0_wp + end subroutine generate_anc_packed +!> sort eigenvalues and eigenvectors pure subroutine sort(nat3,nvar,hess,b) + implicit none - integer :: ii,k,j,m,i - integer, intent(in) :: nat3,nvar + integer, intent(in) :: nat3, nvar real(wp),intent(inout) :: hess(nvar*(nvar+1)/2) real(wp),intent(inout) :: b(nat3,nat3) + real(wp) :: pp,sc1 real(wp),allocatable :: edum(:) + integer :: ii,k,j,m,i + allocate( edum(nvar), source = 0.0_wp ) + ! copy diagonal elements ! do k=1,nvar edum(k)=hess(k+k*(k-1)/2) enddo -! sort + + ! sort eigenvalues and eigenvectors ! do ii = 2, nvar + i = ii - 1 k = i pp= edum(i) + do j = ii, nvar if (edum(j) .gt. pp) cycle k = j pp= edum(j) enddo + if (k .eq. i) cycle + edum(k) = edum(i) edum(i) = pp + do m=1,nat3 sc1=b(m,i) b(m,i)=b(m,k) b(m,k)=sc1 enddo + enddo do k=1,nvar @@ -352,18 +442,63 @@ pure subroutine sort(nat3,nvar,hess,b) end subroutine sort +!> transform and add displacement vector to Cartesian coordinates subroutine get_cartesian(self,xyz) + + use xtb_mctc_io, only : stdout use xtb_mctc_accuracy, only : wp + implicit none - class(tb_anc),intent(in) :: self - integer :: m,i,j,k - real(wp),intent(out) :: xyz (3,self%n) - real(wp) :: dum -! generate cartesian displacement vector + !> optimization state holder + class(tb_anc),intent(inout) :: self + + !> cartesian coordinates to be transformed + real(wp),intent(out) :: xyz(3,self%n) + + !> temporary storage + real(wp), allocatable :: displ_cartesian(:) + integer :: i + + !> debug mode + logical, parameter :: debug = .false. + + allocate(displ_cartesian(3*self%n), source = 0.0_wp) xyz = self%xyz - call dgemv('n',self%n3,self%nvar,1.0_wp,self%B,self%n3,self%coord,1,1.0_wp,xyz,1) + + if (debug) & !####### DEBUG ####### + call self%write(stdout,self%nvar) + + call dgemv('n',self%n3,self%nvar,1.0_wp,self%B,self%n3,self%coord,1,0.0_wp,displ_cartesian,1) ! B * coord + xyz = xyz + reshape(displ_cartesian,(/3,self%n/)) ! xyz + displ_cartesian end subroutine get_cartesian -end module xtb_type_anc +!> transform gradients from Cartesian to normal coordinates +subroutine get_normal(self,g_cartesian, g_normal) + + use xtb_mctc_io, only : stdout + use xtb_mctc_accuracy, only : wp + + implicit none + + !> optimization state holder + class(tb_anc),intent(inout) :: self + + !> cartesian coordinates to be transformed + real(wp),intent(in) :: g_cartesian(3,self%n) + + !> cartesian coordinates to be transformed + real(wp),intent(out) :: g_normal(:) + + !> temporary storage + real(wp), allocatable :: array_form(:) + + allocate(array_form(3*self%n), source = 0.0_wp) + array_form = reshape(g_cartesian,(/3*self%n/)) + + call dgemv('t',self%n3,self%nvar,1.0_wp,self%B,self%n3,array_form,1,0.0_wp,g_normal,1) + +end subroutine get_normal + +end module xtb_type_anc \ No newline at end of file diff --git a/src/type/setvar.f90 b/src/type/setvar.f90 index 4b8501847..385192cef 100644 --- a/src/type/setvar.f90 +++ b/src/type/setvar.f90 @@ -70,39 +70,51 @@ module xtb_type_setvar real(wp) :: broydamp = 0.40_wp end type scc_setvar -!! ------------------------------------------------------------------------ -! approximate normal coordinate rational function optimizer -!! ------------------------------------------------------------------------ + !> approximate normal coordinate rational function optimizer type :: ancopt_setvar -! default optimization level -! crude = -3, sloppy = -2, loose = -1, normal = 0, -! tight = 1, verytight = 2, extreme = 3 + + !> default optimization level + !> crude = -3, sloppy = -2, loose = -1, normal = 0, + !> tight = 1, verytight = 2, extreme = 3 integer :: optlev = 0 -! number of opt. cycles before new ANC are made + + !> number of opt. cycles before new ANC are made by model Hessian integer :: micro_opt = 0 -! total number of opt. cycles, 0 means automatically determined - integer :: maxoptcycle = 0 ! det. in ancopt routine if not read in -! maximum coordinate displacement in ancopt + + !> total number of opt. cycles, 0 means automatically determined + integer :: maxoptcycle = 0 + + !> maximum coordinate displacement in ancopt real(wp) :: maxdispl_opt = 0.0_wp -! lowest force constant in ANC generation (should be > 0.005) + + !> lowest force constant in ANC generation (should be > 0.005) real(wp) :: hlow_opt = 0.0_wp - logical :: exact_rf = .false. + + logical :: exact_rf = .false. + + !> average energy and gradient before checking for convergence + !> to accelerate numerically noisy potential energy surfaces logical :: average_conv = .false. + end type ancopt_setvar type modhess_setvar integer :: model = 0 -! force constants for stretch, bend and torsion + + !> cutoff for constructing internal coordinates + real(wp) :: rcut = 0.0_wp + + !> dispersion scaling in ANC generation + real(wp) :: s6 = 0.0_wp + + !> force constants for stretch, bend and torsion real(wp) :: kr = 0.0_wp real(wp) :: kf = 0.0_wp real(wp) :: kt = 0.0_wp real(wp) :: ko = 0.0_wp real(wp) :: kd = 0.0_wp real(wp) :: kq = 0.0_wp -! cutoff for constructing internal coordinates - real(wp) :: rcut = 0.0_wp -! dispersion scaling in ANC generation - real(wp) :: s6 = 0.0_wp + end type modhess_setvar !! ------------------------------------------------------------------------ diff --git a/test/unit/test_thermo.f90 b/test/unit/test_thermo.f90 index f1b676f5f..8da8499af 100644 --- a/test/unit/test_thermo.f90 +++ b/test/unit/test_thermo.f90 @@ -94,7 +94,7 @@ subroutine test_axis(error) call check(error, rot1(2), .226266337664493_wp, thr=thr2) call check(error, rot1(3), 7.01216792608729_wp, thr=thr2) - call axis2(mol%n, mol%at, mol%xyz, rot2(1), rot2(2), rot2(3), avmom2, mass2) + call axis2(mol%n, mol%xyz, rot2(1), rot2(2), rot2(3), avmom2, mass2) call check(error, rot2(1), .226251131473004_wp, thr=thr2) call check(error, rot2(2), .226266337664493_wp, thr=thr2) @@ -128,7 +128,7 @@ subroutine test_axis(error) call check(error, rot1(2), .19263614502269_wp, thr=thr2) call check(error, rot1(3), 4.7581644454539_wp, thr=thr2) - call axis2(mol%n, mol%at, mol%xyz, rot2(1), rot2(2), rot2(3), avmom2, mass2) + call axis2(mol%n, mol%xyz, rot2(1), rot2(2), rot2(3), avmom2, mass2) call check(error, rot2(1), .19017218374861_wp, thr=thr2) call check(error, rot2(2), .19263614502269_wp, thr=thr2)