diff --git a/src/approxrab.f90 b/src/approxrab.f90 index f0ab00b07..924d200c5 100644 --- a/src/approxrab.f90 +++ b/src/approxrab.f90 @@ -46,7 +46,7 @@ module xtb_approxrab 2.75372921, 2.62540906, 2.55860939, 3.32492356, 2.65140898,& ! 80 1.52014458, 2.54984804, 1.72021963, 2.69303422, 1.81031095,& 2.34224386,& - 2.52387890,& !@thomas TODO + 2.52387890,& 2.30204667, 1.60119300, 2.00000000, 2.00000000, 2.00000000,& ! 92 2.00000000, 2.00000000, 2.00000000, 2.00000000, 2.00000000,& 2.00000000, 2.00000000, 2.00000000, 2.00000000, 2.00000000,& ! 102 @@ -71,7 +71,7 @@ module xtb_approxrab 2.33473532, 2.19498900, 2.12678348, 2.34895048, 2.33422774,& 2.86560827, 2.62488837, 2.88376127, 2.75174124, 2.83054552,& 2.63264944,& - 4.24537037,& !@thomas TODO + 4.24537037,& 3.66542289, 4.20000000, 4.20000000, 4.20000000, 4.20000000,& ! 92 4.20000000, 4.20000000, 4.20000000, 4.20000000, 4.20000000,& 4.20000000, 4.20000000, 4.20000000, 4.20000000, 4.20000000,& ! 102 @@ -96,7 +96,7 @@ module xtb_approxrab 0.04671159, 0.06758819, 0.09488437, 0.07556405, 0.13384502,& ! 90 0.03203572, 0.04235009, 0.03153769,-0.00152488, 0.02714675,& 0.04800662,& - 0.04582912,& !@thomas TODO + 0.04582912,& 0.10557321, 0.02167468, 0.05463616, 0.05370913, 0.05985441,& ! 92 0.02793994, 0.02922983, 0.02220438, 0.03340460,-0.04110969,& -0.01987240, 0.07260201, 0.07700000, 0.07700000, 0.07700000,& ! 102 diff --git a/src/geoopt_driver.f90 b/src/geoopt_driver.f90 index b05bb942f..33c4359c6 100644 --- a/src/geoopt_driver.f90 +++ b/src/geoopt_driver.f90 @@ -162,7 +162,7 @@ subroutine geometry_optimization & ! create new Limited-memory BFGS optimizer call new_lbfgs_optimizer(lbfgs_opt, env, opt_input, filter) ! run optimization - call relax_pbc(lbfgs_opt, env, mol, wfn, calc, filter, printlevel) + call relax_pbc(lbfgs_opt, env, mol, wfn, calc, filter, printlevel, fail) case(p_engine_inertial) call fire & ! FIRE ! &(env,ilog,mol,wfn,calc, & diff --git a/src/gfnff/gdisp0.f90 b/src/gfnff/gdisp0.f90 index 48aead5e4..042572b21 100644 --- a/src/gfnff/gdisp0.f90 +++ b/src/gfnff/gdisp0.f90 @@ -130,6 +130,7 @@ subroutine weight_references_d4(dispm, nat, atoms, wf, cn, gwvec, gwdcn) norm = norm + gw dnorm = dnorm + 2*wf*(dispm%cn(iref, ati) - cn(iat)) * gw end do + if (norm.eq.0.0_wp) cycle norm = 1.0_wp / norm do iref = 1, dispm%nref(ati) expw = weight_cn(wf, cn(iat), dispm%cn(iref, ati)) diff --git a/src/gfnff/gfnff_ini.f90 b/src/gfnff/gfnff_ini.f90 index 0dcc2c80b..2335b00d6 100644 --- a/src/gfnff/gfnff_ini.f90 +++ b/src/gfnff/gfnff_ini.f90 @@ -2207,6 +2207,7 @@ subroutine adjust_NB_LnH_AnH(param, mol, topo, neigh) type(TNeigh), intent(inout) :: neigh ! main type for introducing PBC integer nb_tmp(neigh%numnb) integer :: i,iTr,iTrH,idx,inb,l,k,m,nnb,count_idx + real(wp), parameter :: qthr = -0.0281_wp ! charge threshold ! loop over all atoms do i=1, mol%n @@ -2224,7 +2225,7 @@ subroutine adjust_NB_LnH_AnH(param, mol, topo, neigh) ! check if neighbor is H if (mol%at(inb).eq.1) then ! remove hydrogen as neighbor if charge is gt threshold - if (topo%qa(inb).gt.-0.0282) then + if (topo%qa(inb).gt.qthr) then ! setup copy of neighbor list for this An or Ln nb_tmp = neigh%nb(:,i,iTr) nb_tmp(neigh%numnb) = neigh%nb(neigh%numnb,i,iTr) - 1 @@ -2452,10 +2453,10 @@ subroutine gfnff_topo_changes(env, neigh) ! check if hardcoded size of ffnb is still up to date if (size(set%ffnb, dim=1).ne.neigh%numnb) call env%error('The array set%ffnb has not been adjusted to changes in neigh%numnb.', source) ! only do something if there are changes stored in set%ffnb - if(set%ffnb(1,1).ne.0) then + if(set%ffnb(1,1).ne.-1) then d2=size(set%ffnb, dim=2) do i=1, d2 - if (set%ffnb(1,i).eq.0) exit + if (set%ffnb(1,i).eq.-1) exit idx=set%ffnb(1,i) int_tmp = set%ffnb(2:41,i) neigh%nb(1:40,idx,1) = int_tmp diff --git a/src/gfnff/gfnff_ini2.f90 b/src/gfnff/gfnff_ini2.f90 index 24d010fd8..09e16f2b5 100644 --- a/src/gfnff/gfnff_ini2.f90 +++ b/src/gfnff/gfnff_ini2.f90 @@ -43,7 +43,8 @@ subroutine gfnff_neigh(env,makeneighbor,natoms,at,xyz,rab,fq,f_in,f2_in,lintr, & type(TMolecule), intent(in) :: mol type(TNeigh), intent(inout) :: neigh ! contains nb, nbf and nbm logical, intent(in) :: makeneighbor, nb_call - integer at(natoms),natoms + integer, intent(in) :: natoms + integer at(natoms) integer hyb (natoms) integer itag(natoms) real*8 rab (natoms*(natoms+1)/2) @@ -180,7 +181,7 @@ subroutine gfnff_neigh(env,makeneighbor,natoms,at,xyz,rab,fq,f_in,f2_in,lintr, & elseif(nm.eq.1)then ! distinguish M-CR2-R i.e. not an eta coord. ncm=0 do iTr=1, numctr - do k=1,numnbf ! + do k=1,neigh%nbf(neigh%numnb,i,iTr) ! if(neigh%nbf(k,i,iTr).ne.im)then ! all neighbors that are not the metal im kk=neigh%nbf(k,i,iTr) do l=1,sum(neigh%nbf(neigh%numnb,kk,:)) diff --git a/src/gfnff/neighbor.f90 b/src/gfnff/neighbor.f90 index 843d80cf0..e07006431 100644 --- a/src/gfnff/neighbor.f90 +++ b/src/gfnff/neighbor.f90 @@ -115,10 +115,10 @@ subroutine get_nb(self, mol, rab, rtmp, mchar, icase, f_in, f2_in, param) !integer, intent(inout) :: nbf(20,mol%n) type(TGFFData), intent(in) :: param !real(wp) :: latThresh, maxNBr - real(wp) :: dist(mol%n,mol%n,self%numctr) + real(wp), allocatable :: dist(:,:,:) integer :: i,j,iTr - dist=0.0_wp + allocate(dist(mol%n,mol%n,self%numctr), source=0.0_wp) !$omp parallel do collapse(3) default(none) shared(dist,mol,self) & !$omp private(iTr,i,j) do iTr=1, self%numctr @@ -404,14 +404,15 @@ end function id2v subroutine fillnb(self,n,at,rad,dist,mchar,icase,f,f2,param) implicit none type(TGFFData), intent(in) :: param - class(TNeigh), intent(inout) :: self + type(TNeigh), intent(inout) :: self integer, intent(in) :: n,at(n) real(wp), intent(in) :: rad(n*(n+1)/2) real(wp), intent(in) :: dist(n,n,self%numctr) real(wp), intent(in) :: mchar(n),f,f2 integer i,j,k,iTr,nn,icase,hc_crit,nnfi,nnfj,lin - integer tag(n,n,self%numctr) + integer, allocatable :: tag(:,:,:) real(wp) rco,fm + allocate(tag(n,n,self%numctr), source=0) if(icase.eq.1) then if(.not. allocated(self%nbf)) then allocate(self%nbf(self%numnb,n,self%numctr),source=0) @@ -434,7 +435,6 @@ subroutine fillnb(self,n,at,rad,dist,mchar,icase,f,f2,param) endif endif - tag = 0 do iTr=1, self%numctr do i=1,n do j=1,n diff --git a/src/lbfgs_anc/driver.f90 b/src/lbfgs_anc/driver.f90 index 78b486424..9e28bd613 100644 --- a/src/lbfgs_anc/driver.f90 +++ b/src/lbfgs_anc/driver.f90 @@ -36,7 +36,7 @@ module xtb_pbc_optimizer_driver !> Driver for performing geometry optimization !subroutine relax(self, ctx, mol, calc, filter, accuracy, verbosity) -subroutine relax_pbc(self, env, mol, chk, calc, filter, printlevel) +subroutine relax_pbc(self, env, mol, chk, calc, filter, printlevel, optfail) use xtb_pbc, only: cross_product use xtb_gfnff_calculator, only : TGFFCalculator, newGFFCalculator !> Instance of the optimization driver @@ -64,6 +64,8 @@ subroutine relax_pbc(self, env, mol, chk, calc, filter, printlevel) type(scc_results) :: results !> optimization cycle failed logical :: fail + !> optimization failed + logical, intent(out) :: optfail !> singlepoint energy real(wp) :: energy !> gradient and stress tensor @@ -84,6 +86,7 @@ subroutine relax_pbc(self, env, mol, chk, calc, filter, printlevel) real(wp) :: emin_global, xyz_emin(3,mol%n), latt_emin(3,3) character(6), parameter :: fname = "noName" + optfail = .false. nvar = filter%get_dimension() allocate(gcurr(nvar), glast(nvar), displ(nvar)) allocate(gxyz(3, mol%n), sigma(3, 3), source=0.0_wp) @@ -203,6 +206,7 @@ subroutine relax_pbc(self, env, mol, chk, calc, filter, printlevel) ! Output message if not converged after maximum steps if (.not.converged) then + optfail = .true. write(*,*) '' write(*,*) 'Could not converge in',step,' steps.' call env%error("Could not converge geometry.",source) diff --git a/src/set_module.f90 b/src/set_module.f90 index 889a54ef9..c81703d56 100644 --- a/src/set_module.f90 +++ b/src/set_module.f90 @@ -1508,16 +1508,18 @@ subroutine set_ffnb(env,key,val) character(len=*),intent(in) :: val integer :: i,j,k,l integer :: i_start,i_end,i_ffnb + logical :: nonb k=1 ! start at 1 since first entry of ffnb is the atom index that the following NBs belong to i_start=1 i_end=1 i_ffnb=0 + nonb = .false. ! expect that the atom should have neighbors do i=1, len(val) ! get next empty row in ffnb and read atom index into first entry if (val(i:i).eq.":") then do j=1,size(set%ffnb, dim=2) - if (set%ffnb(1,j).eq.0 .and. i_ffnb.eq.0) then + if (set%ffnb(1,j).eq.-1 .and. i_ffnb.eq.0) then i_ffnb = j ! index of next empty row l=i-1 ! we take care of ":" and "," and trust read() to handle whitespaces @@ -1531,17 +1533,30 @@ subroutine set_ffnb(env,key,val) if (val(i:i).eq.",") then k = k + 1 i_end=i-1 - read(val(i_start:i_end), *) set%ffnb(k,i_ffnb) - i_start=i+1 + read(val(i_start:i_end), *) set%ffnb(k,i_ffnb) + ! if the first neighbor index (k=2) is zero the atom (k=1) has no neighbors + if (set%ffnb(k,i_ffnb) == 0 .and. k == 2) then + nonb = .true. + endif + i_start=i+1 endif ! read the last neighbor into ffnb if (i.eq.len(val)) then k = k + 1 read(val(i_start:), *) set%ffnb(k,i_ffnb) + ! if the first neighbor index (k=2) is zero the atom (k=1) has no neighbors + if (set%ffnb(k,i_ffnb) == 0 .and. k == 2) then + nonb = .true. + endif + ! set rest of ffnb to 0 + set%ffnb(k+1:41, i_ffnb) = 0 endif enddo - set%ffnb(42,i_ffnb) = k - 1 ! number of neighbors of atom set%ffnb(1,i_ffnb) - + if (nonb) then + set%ffnb(2:,i_ffnb) = 0 + else + set%ffnb(42,i_ffnb) = k - 1 ! number of neighbors of atom set%ffnb(1,i_ffnb) + endif end subroutine set_ffnb diff --git a/src/setparam.f90 b/src/setparam.f90 index c209e2edc..13047c724 100644 --- a/src/setparam.f90 +++ b/src/setparam.f90 @@ -519,7 +519,7 @@ module xtb_setparam !> GFN-FF manual setup of nb list via xcontrol ! allows a maximum of 164 atoms neighbors to be changed ! ffnb(42,i) stores the number of neighbors of atom i - integer :: ffnb(42,164) = 0 + integer :: ffnb(42,164) = -1 end type TSet type(TSet) :: set diff --git a/test/unit/molstock.f90 b/test/unit/molstock.f90 index ed4759afc..2123d48e4 100644 --- a/test/unit/molstock.f90 +++ b/test/unit/molstock.f90 @@ -57,6 +57,8 @@ subroutine getMolecule(mol, name) case ('mgh2'); call MgH2(mol) case ('x06_b'); call x06_benzene(mol) case ('mcv15'); call mcv15(mol) + case ('Th_15'); call Th_1519394(mol) + case ('Ce_0a7745'); call Ce_0a7745(mol) end select end subroutine getMolecule @@ -1264,6 +1266,350 @@ subroutine mcv15(mol) call init(mol, sym, xyz, charge, uhf, lattice, pbc) end subroutine mcv15 + subroutine Th_1519394(mol) + type(TMolecule), intent(out) :: mol + integer, parameter :: nat = 270 + character(len=*), parameter :: sym(nat) = [character(len=4) :: & + & "th","th","n", "n", "n", "n","si","si","si","si","si","si","si","si","o","o","o","o","c","c","h","h","h","h","h","h","c", & + & "c", "h", "h", "h", "h", "h", "h", "c", "c", "h", "h", "h", "h", "h","h","c","c","c","c","h","h","c","c","h","h","c","c", & + & "c", "c", "c", "c", "c", "c", "h", "h", "c", "c", "c", "c", "h", "h","c","c","c","c","c","c","h","h","c","c","h","h","c", & + & "c", "h", "h", "c", "c", "h", "h", "c", "c", "h", "h", "c", "c", "h","h","c","c","h","h","c","c","h","h","c","c","h","h", & + & "c", "c", "h", "h", "c", "c", "h", "h", "c", "c", "h", "h", "c", "c","h","h","c","c","h","h","c","c","h","h","h","h","h", & + & "h", "c", "c", "h", "h", "h", "h", "h", "h", "c", "c", "h", "h", "h","h","h","h","c","c","h","h","c","c","h","h","c","c", & + & "h", "h", "c", "c", "h", "h", "c", "c", "h", "h", "c", "c", "h", "h","c","c","h","h","c","c","h","h","c","c","c","c","h", & + & "h", "c", "c", "h", "h", "h", "h", "h", "h", "c", "c", "h", "h", "h","h","h","h","c","c","h","h","h","h","h","h","c","c", & + & "h", "h", "c", "c", "h", "h", "c", "c", "h", "h", "c", "c", "h", "h","h","h","h","h","c","c","h","h","h","h","h","h","c", & + & "c", "h", "h", "h", "h", "h", "h", "c", "c", "h", "h", "h", "h", "c","c","h","h","h","h","c","c","h","h","h","h","h","h" ] + real(wp), parameter :: xyz(3, nat) = reshape([& + & 16.4817133_wp, 19.3954377_wp, 22.6174173_wp, & + & 19.6280928_wp, 5.91957031_wp, 8.44137637_wp, & + & 19.8241200_wp, 17.0053068_wp, 20.8957352_wp, & + & 16.2856861_wp, 8.30970120_wp, 10.1630584_wp, & + & 13.6007319_wp, 21.8607354_wp, 20.2469170_wp, & + & 22.5090743_wp, 3.45427265_wp, 10.8118766_wp, & + & 20.4959967_wp, 13.8364585_wp, 21.3808736_wp, & + & 15.6138095_wp, 11.4785495_wp, 9.67792013_wp, & + & 8.27156201_wp, 0.77721684_wp, 20.8867282_wp, & + & 27.8382441_wp, 24.5377912_wp, 10.1720655_wp, & + & 21.6515396_wp, 18.8241709_wp, 18.8828148_wp, & + & 14.4582665_wp, 6.49083710_wp, 12.1759789_wp, & + & 12.3949139_wp, 20.2296757_wp, 17.7115877_wp, & + & 23.7148922_wp, 5.08533237_wp, 13.3472060_wp, & + & 18.5780928_wp, 21.5825152_wp, 25.2703663_wp, & + & 17.5317134_wp, 3.73249278_wp, 5.78842739_wp, & + & 13.9284701_wp, 17.2358969_wp, 24.8346114_wp, & + & 22.1813360_wp, 8.07911115_wp, 6.22418227_wp, & + & 8.64853421_wp, 1.41152752_wp, 24.3252472_wp, & + & 27.4612719_wp, 23.9034805_wp, 6.73354648_wp, & + & 12.1094242_wp, 24.4599706_wp, 25.2880698_wp, & + & 24.0003819_wp, 0.85503747_wp, 5.77072388_wp, & + & 8.17965451_wp, 3.16620524_wp, 24.6855292_wp, & + & 27.9301516_wp, 22.1488028_wp, 6.37326448_wp, & + & 10.4130731_wp, 1.13179002_wp, 24.8159762_wp, & + & 25.6967330_wp, 24.1832180_wp, 6.24281754_wp, & + & 24.8634746_wp, 19.6520502_wp, 20.1074630_wp, & + & 11.2463315_wp, 5.66295783_wp, 10.9513306_wp, & + & 24.6978153_wp, 20.5474271_wp, 21.7194144_wp, & + & 11.4119908_wp, 4.76758091_wp, 9.33937928_wp, & + & 25.7302282_wp, 20.7373232_wp, 18.8837466_wp, & + & 10.3795779_wp, 4.57768482_wp, 12.1750471_wp, & + & 25.8444155_wp, 18.1010833_wp, 20.3621451_wp, & + & 10.2653906_wp, 7.21392472_wp, 10.6966485_wp, & + & 20.0008905_wp, 21.9541389_wp, 18.4861940_wp, & + & 16.1089156_wp, 3.36086916_wp, 12.5725997_wp, & + & 18.2878845_wp, 21.6720095_wp, 17.8432770_wp, & + & 17.8219216_wp, 3.64299850_wp, 13.2155167_wp, & + & 20.9436483_wp, 23.0014974_wp, 17.2842187_wp, & + & 15.1661578_wp, 2.31351064_wp, 13.7745750_wp, & + & 19.9096982_wp, 22.8164391_wp, 20.1198865_wp, & + & 16.2001079_wp, 2.49856891_wp, 10.9389071_wp, & + & 12.5057544_wp, 13.0880009_wp, 25.9371986_wp, & + & 23.6040517_wp, 12.2270071_wp, 5.12159509_wp, & + & 15.0643565_wp, 19.7864888_wp, 0.12734105_wp, & + & 21.0454497_wp, 5.52851918_wp, 30.9314526_wp, & + & 14.6697944_wp, 21.3351769_wp, 0.94418733_wp, & + & 21.4400118_wp, 3.97983112_wp, 30.1146064_wp, & + & 14.5762345_wp, 2.37951694_wp, 29.8754537_wp, & + & 21.5335716_wp, 22.9354911_wp, 1.18334004_wp, & + & 16.0938634_wp, 3.14714139_wp, 29.2977601_wp, & + & 20.0159427_wp, 22.1678666_wp, 1.76103361_wp, & + & 18.4441588_wp, 24.2114237_wp, 28.9747486_wp, & + & 17.6656473_wp, 1.10358437_wp, 2.08404506_wp, & + & 20.1395262_wp, 22.6798111_wp, 27.1609151_wp, & + & 15.9702799_wp, 2.63519694_wp, 3.89787861_wp, & + & 12.2742914_wp, 15.9160074_wp, 26.5459510_wp, & + & 23.8355147_wp, 9.39900058_wp, 4.51284273_wp, & + & 4.51156592_wp, 16.3343038_wp, 0.04969407_wp, & + & 31.5982402_wp, 8.98070421_wp, 31.0090996_wp, & + & 34.3331051_wp, 17.7370508_wp, 30.7388881_wp, & + & 1.77670110_wp, 7.57795725_wp, 0.31990558_wp, & + & 21.4435420_wp, 20.5554255_wp, 28.6796901_wp, & + & 14.6662641_wp, 4.75958250_wp, 2.37910360_wp, & + & 23.2444153_wp, 16.3689110_wp, 28.8971017_wp, & + & 12.8653908_wp, 8.94609700_wp, 2.16169204_wp, & + & 23.5791781_wp, 14.7875966_wp, 28.1144201_wp, & + & 12.5306280_wp, 10.5274113_wp, 2.94437365_wp, & + & 22.1055709_wp, 24.3688537_wp, 25.8440222_wp, & + & 14.0042353_wp, 0.94615431_wp, 5.21477147_wp, & + & 13.1656210_wp, 16.5260927_wp, 29.2294307_wp, & + & 22.9441851_wp, 8.78891534_wp, 1.82936295_wp, & + & 16.2273261_wp, 23.1146667_wp, 29.8816654_wp, & + & 19.8824800_wp, 2.20034132_wp, 1.17712828_wp, & + & 15.7453867_wp, 21.4873418_wp, 29.2977601_wp, & + & 20.3644194_wp, 3.82766619_wp, 1.76103361_wp, & + & 24.6741664_wp, 24.2097241_wp, 26.3254335_wp, & + & 11.4356397_wp, 1.10528392_wp, 4.73336017_wp, & + & 25.2815040_wp, 23.0044607_wp, 27.5087736_wp, & + & 10.8283021_wp, 2.31054732_wp, 3.55002012_wp, & + & 9.55551602_wp, 15.5674812_wp, 1.18334004_wp, & + & 26.5542901_wp, 9.74752679_wp, 29.8754537_wp, & + & 11.2913450_wp, 15.3054608_wp, 1.56225733_wp, & + & 24.8184611_wp, 10.0095472_wp, 29.4965364_wp, & + & 15.7218572_wp, 16.1690660_wp, 29.8568184_wp, & + & 20.3879489_wp, 9.14594205_wp, 1.20197532_wp, & + & 16.8520406_wp, 15.5528297_wp, 28.6051490_wp, & + & 19.2577656_wp, 9.76217835_wp, 2.45364471_wp, & + & 16.3244938_wp, 17.9186718_wp, 1.45044567_wp, & + & 19.7853123_wp, 7.39633621_wp, 29.6083480_wp, & + & 16.8144037_wp, 18.2078602_wp, 3.15557344_wp, & + & 19.2954024_wp, 7.10714781_wp, 27.9032203_wp, & + & 5.97973456_wp, 2.51484087_wp, 0.57769356_wp, & + & 30.1300716_wp, 22.8001672_wp, 30.4811001_wp, & + & 6.46260731_wp, 4.13097188_wp, 1.18644592_wp, & + & 29.6471988_wp, 21.1840361_wp, 29.8723478_wp, & + & 16.8776846_wp, 15.6375511_wp, 0.31369382_wp, & + & 19.2321215_wp, 9.67745694_wp, 30.7450999_wp, & + & 17.7363311_wp, 14.3547630_wp, 1.22682235_wp, & + & 18.3734750_wp, 10.9602450_wp, 29.8319713_wp, & + & 16.7075220_wp, 1.95425792_wp, 24.0829886_wp, & + & 19.4022841_wp, 23.3607501_wp, 6.97580507_wp, & + & 14.9526263_wp, 2.08380107_wp, 23.7196007_wp, & + & 21.1571799_wp, 23.2312070_wp, 7.33919296_wp, & + & 7.93576576_wp, 16.4761400_wp, 3.02823239_wp, & + & 28.1740404_wp, 8.83886798_wp, 28.0305613_wp, & + & 8.55300001_wp, 16.8310314_wp, 4.67745434_wp, & + & 27.5568061_wp, 8.48397658_wp, 26.3813394_wp, & + & 22.0427437_wp, 18.2789174_wp, 27.5367265_wp, & + & 14.0670624_wp, 7.03609058_wp, 3.52206721_wp, & + & 21.6312331_wp, 18.0155490_wp, 25.8098576_wp, & + & 14.4785730_wp, 7.29945905_wp, 5.24893614_wp, & + & 12.1205797_wp, 12.3397613_wp, 23.4462834_wp, & + & 23.9892264_wp, 12.9752467_wp, 7.61251035_wp, & + & 11.6922860_wp, 13.5697187_wp, 22.2101434_wp, & + & 24.4175201_wp, 11.7452892_wp, 8.84865034_wp, & + & 7.64276815_wp, 23.2198378_wp, 0.55284653_wp, & + & 28.4670380_wp, 2.09517022_wp, 30.5059472_wp, & + & 6.13581692_wp, 22.4480511_wp, 1.14917537_wp, & + & 29.9739892_wp, 2.86695698_wp, 29.9096183_wp, & + & 13.0133989_wp, 11.2293649_wp, 27.7137616_wp, & + & 23.0964072_wp, 14.0856431_wp, 3.34503209_wp, & + & 13.2276122_wp, 11.6837149_wp, 29.4375247_wp, & + & 22.8821939_wp, 13.6312931_wp, 1.62126903_wp, & + & 5.40598312_wp, 16.8690581_wp, 2.45675059_wp, & + & 30.7038230_wp, 8.44594992_wp, 28.6020431_wp, & + & 4.28878075_wp, 17.5022268_wp, 3.70841997_wp, & + & 31.8210254_wp, 7.81278121_wp, 27.3503737_wp, & + & 21.9917352_wp, 17.4112350_wp, 15.6784790_wp, & + & 14.1180709_wp, 7.90377307_wp, 15.3803146_wp, & + & 22.8998354_wp, 15.8048548_wp, 15.8058201_wp, & + & 13.2099707_wp, 9.51015322_wp, 15.2529736_wp, & + & 22.9368436_wp, 18.5739166_wp, 14.5883154_wp, & + & 13.1729625_wp, 6.74109139_wp, 16.4704783_wp, & + & 20.3116677_wp, 17.1163901_wp, 14.9517033_wp, & + & 15.7981384_wp, 8.19861792_wp, 16.1070904_wp, & + & 19.1079269_wp, 12.6791087_wp, 24.3997883_wp, & + & 17.0018793_wp, 12.6358993_wp, 6.65900538_wp, & + & 20.0235779_wp, 13.4250145_wp, 25.8253870_wp, & + & 16.0862282_wp, 11.8899935_wp, 5.23340675_wp, & + & 19.2445079_wp, 10.8337611_wp, 24.4743294_wp, & + & 16.8652983_wp, 14.4812469_wp, 6.58446427_wp, & + & 17.3248810_wp, 13.1708817_wp, 24.4867529_wp, & + & 18.7849251_wp, 12.1441263_wp, 6.57204076_wp, & + & 19.1095423_wp, 11.7881322_wp, 18.9054877_wp, & + & 17.0002638_wp, 13.5268758_wp, 12.1533059_wp, & + & 17.2818644_wp, 12.0706367_wp, 18.8247348_wp, & + & 18.8279418_wp, 13.2443713_wp, 12.2340588_wp, & + & 19.4436158_wp, 10.0129979_wp, 19.3154638_wp, & + & 16.6661903_wp, 15.3020101_wp, 11.7433299_wp, & + & 19.8729327_wp, 12.1980677_wp, 17.2686893_wp, & + & 16.2368734_wp, 13.1169402_wp, 13.7901044_wp, & + & 20.9675058_wp, 3.36716584_wp, 23.3562128_wp, & + & 15.1423003_wp, 21.9478422_wp, 7.70258085_wp, & + & 22.1234677_wp, 4.45166106_wp, 22.5114137_wp, & + & 13.9863384_wp, 20.8633470_wp, 8.54738004_wp, & + & 21.8375205_wp, 1.62272370_wp, 25.0799759_wp, & + & 14.2722856_wp, 23.6922843_wp, 5.97881780_wp, & + & 23.5926940_wp, 1.49987436_wp, 25.4278344_wp, & + & 12.5171121_wp, 23.8151337_wp, 5.63095931_wp, & + & 18.4091842_wp, 3.53771449_wp, 22.8561663_wp, & + & 17.7006219_wp, 21.7772935_wp, 8.20262743_wp, & + & 17.8102244_wp, 4.74275272_wp, 21.6666145_wp, & + & 18.2995817_wp, 20.5722553_wp, 9.39217923_wp, & + & 3.78539572_wp, 1.40604394_wp, 1.41938687_wp, & + & 32.3244104_wp, 23.9089641_wp, 29.6394068_wp, & + & 2.74145793_wp, 2.26872621_wp, 2.59651516_wp, & + & 33.3683482_wp, 23.0462818_wp, 28.4622785_wp, & + & 12.9077485_wp, 8.01038440_wp, 24.5271294_wp, & + & 23.2020576_wp, 17.3046236_wp, 6.53166432_wp, & + & 13.0855448_wp, 6.28782465_wp, 24.0519298_wp, & + & 23.0242613_wp, 19.0271834_wp, 7.00686387_wp, & + & 13.2105660_wp, 8.70715130_wp, 26.9994094_wp, & + & 22.8992401_wp, 16.6078567_wp, 4.05938434_wp, & + & 13.5600309_wp, 7.45795512_wp, 28.2386552_wp, & + & 22.5497752_wp, 17.8570529_wp, 2.82013847_wp, & + & 12.3490360_wp, 9.80534053_wp, 22.7288252_wp, & + & 23.7607701_wp, 15.5096675_wp, 8.32996848_wp, & + & 12.1196937_wp, 9.32472022_wp, 21.0143798_wp, & + & 23.9901124_wp, 15.9902878_wp, 10.0444139_wp, & + & 33.3494870_wp, 18.8789064_wp, 24.8035526_wp, & + & 2.76031916_wp, 6.43610167_wp, 6.25524106_wp, & + & 10.1576736_wp, 19.9436614_wp, 24.1482121_wp, & + & 25.9521325_wp, 5.37134666_wp, 6.91058161_wp, & + & 9.52332313_wp, 16.6786849_wp, 26.0987043_wp, & + & 26.5864830_wp, 8.63632311_wp, 4.96008936_wp, & + & 32.0655798_wp, 15.1969284_wp, 27.0615269_wp, & + & 4.04422638_wp, 10.1180796_wp, 3.99726675_wp, & + & 32.4737602_wp, 13.6982941_wp, 27.9591261_wp, & + & 3.63604592_wp, 11.6167139_wp, 3.09966762_wp, & + & 13.9428210_wp, 17.0114357_wp, 17.7283594_wp, & + & 22.1669851_wp, 8.30357230_wp, 13.3304342_wp, & + & 15.7714787_wp, 17.2061804_wp, 17.4954185_wp, & + & 20.3383274_wp, 8.10882765_wp, 13.5633752_wp, & + & 13.2493431_wp, 15.9866248_wp, 16.3493490_wp, & + & 22.8604630_wp, 9.32838321_wp, 14.7094447_wp, & + & 13.6115786_wp, 16.1750959_wp, 19.3465226_wp, & + & 22.4982275_wp, 9.13991211_wp, 11.7122711_wp, & + & 8.92566846_wp, 19.5807174_wp, 17.8060064_wp, & + & 27.1841377_wp, 5.73429062_wp, 13.2527872_wp, & + & 8.52499886_wp, 18.5905770_wp, 19.3185697_wp, & + & 27.5848073_wp, 6.72443103_wp, 11.7402240_wp, & + & 8.43494146_wp, 18.6285830_wp, 16.2965490_wp, & + & 27.6748647_wp, 6.68642499_wp, 14.7622446_wp, & + & 7.99975850_wp, 21.1838394_wp, 17.8557005_wp, & + & 28.1100477_wp, 4.13116863_wp, 13.2030932_wp, & + & 23.9660621_wp, 13.1652927_wp, 21.5299558_wp, & + & 12.1437440_wp, 12.1497153_wp, 9.52883792_wp, & + & 24.7395053_wp, 13.5412509_wp, 19.8900515_wp, & + & 11.3703008_wp, 11.7737571_wp, 11.1687422_wp, & + & 24.2336286_wp, 11.3807079_wp, 21.9461436_wp, & + & 11.8761776_wp, 13.9343000_wp, 9.11265008_wp, & + & 24.7455805_wp, 14.2237978_wp, 22.8375310_wp, & + & 11.3642256_wp, 11.0912102_wp, 8.22126270_wp, & + & 28.9369548_wp, 18.0418670_wp, 25.4682108_wp, & + & 7.17285137_wp, 7.27314107_wp, 5.59058287_wp, & + & 27.2156057_wp, 18.5188944_wp, 25.2818581_wp, & + & 8.89420048_wp, 6.79611366_wp, 5.77693564_wp, & + & 29.5768601_wp, 15.8435015_wp, 26.7478331_wp, & + & 6.53294606_wp, 9.47150653_wp, 4.31096057_wp, & + & 28.2871598_wp, 14.7838095_wp, 27.4093854_wp, & + & 7.82264633_wp, 10.5311985_wp, 3.64940826_wp, & + & 30.8359788_wp, 19.5269976_wp, 24.4650118_wp, & + & 5.27382733_wp, 5.78801038_wp, 6.59378191_wp, & + & 30.4153401_wp, 21.0052674_wp, 23.5394597_wp, & + & 5.69446602_wp, 4.30974067_wp, 7.51933397_wp, & + & 10.3694938_wp, 3.00326359_wp, 19.1228993_wp, & + & 25.7403123_wp, 22.3117444_wp, 11.9358944_wp, & + & 12.1331631_wp, 2.70383701_wp, 19.6043106_wp, & + & 23.9766430_wp, 22.6111710_wp, 11.4544831_wp, & + & 9.90546554_wp, 4.74776072_wp, 19.5359812_wp, & + & 26.2043406_wp, 20.5672473_wp, 11.5228124_wp, & + & 10.1692546_wp, 2.72672925_wp, 17.2997481_wp, & + & 25.9405515_wp, 22.5882788_wp, 13.7590456_wp, & + & 4.93529212_wp, 1.63075053_wp, 20.0298160_wp, & + & 31.1745140_wp, 23.6842575_wp, 11.0289776_wp, & + & 4.65881696_wp, 1.31388278_wp, 18.2284060_wp, & + & 31.4509892_wp, 24.0011252_wp, 12.8303876_wp, & + & 29.1271684_wp, 3.42424344_wp, 20.3932039_wp, & + & 6.98263771_wp, 21.8907646_wp, 10.6655897_wp, & + & 32.7961403_wp, 24.7917020_wp, 21.0268033_wp, & + & 3.31366589_wp, 0.52330599_wp, 10.0319903_wp, & + & 13.1301069_wp, 21.7530442_wp, 14.6069507_wp, & + & 22.9796993_wp, 3.56196381_wp, 16.4518430_wp, & + & 12.3833913_wp, 23.4493199_wp, 14.5572566_wp, & + & 23.7264148_wp, 1.86568809_wp, 16.5015371_wp, & + & 12.4262585_wp, 20.7308774_wp, 13.2341520_wp, & + & 23.6835476_wp, 4.58413064_wp, 17.8246417_wp, & + & 14.9668141_wp, 21.8833789_wp, 14.3988567_wp, & + & 21.1429920_wp, 3.43162912_wp, 16.6599369_wp, & + & 29.7702042_wp, 12.2918867_wp, 16.7872780_wp, & + & 6.33960200_wp, 13.0231213_wp, 14.2715157_wp, & + & 28.0720651_wp, 11.5494986_wp, 16.5263841_wp, & + & 8.03774107_wp, 13.7655094_wp, 14.5324095_wp, & + & 29.5538680_wp, 13.8408614_wp, 17.8122182_wp, & + & 6.55593812_wp, 11.4741466_wp, 13.2465755_wp, & + & 6.85725737_wp, 10.4490503_wp, 18.2656766_wp, & + & 29.2525488_wp, 14.8659576_wp, 12.7931171_wp, & + & 8.55037771_wp, 11.1965014_wp, 18.5327822_wp, & + & 27.5594284_wp, 14.1185066_wp, 12.5260115_wp, & + & 7.08043417_wp, 8.90468835_wp, 17.2345246_wp, & + & 29.0293720_wp, 16.4103197_wp, 13.8242690_wp, & + & 30.2602414_wp, 9.67663966_wp, 20.8466623_wp, & + & 5.84956479_wp, 15.6383684_wp, 10.2121313_wp, & + & 30.0239765_wp, 11.1852809_wp, 21.8933437_wp, & + & 6.08582963_wp, 14.1297271_wp, 9.16545003_wp, & + & 6.95557442_wp, 8.51935187_wp, 21.6852497_wp, & + & 29.1542317_wp, 16.7956562_wp, 9.37354395_wp, & + & 28.6282565_wp, 8.83102642_wp, 20.6075096_wp, & + & 7.48154961_wp, 16.4839816_wp, 10.4512840_wp], shape(xyz)) + real(wp), parameter :: lattice(3, 3) = reshape([& + & 24.48139122_wp, 0.00000000_wp, 0.00000000_wp, & + & 4.55484485_wp, 24.18913280_wp, 0.00000000_wp, & + & 7.07357013_wp, 1.12587527_wp, 31.05879374_wp & + ], shape(lattice)) + real(wp), parameter :: charge = 0.0_wp + integer, parameter :: uhf = 0 + logical, parameter :: pbc(3) = [.true., .true., .true.] + call init(mol, sym, xyz, charge, uhf, lattice, pbc) + end subroutine Th_1519394 + + subroutine Ce_0a7745(mol) + type(TMolecule), intent(out) :: mol + integer, parameter :: nat = 33 + character(len=*), parameter :: sym(nat) = [character(len=4) ::& + & "Ce", "N", "C", "C", "N", "C", "C", "N", "H", "H", & + & "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", & + & "H", "N", "O", "O", "O", "C", "N", "C", "O", "H", & + & "H", "H", "H" ] + real(wp), parameter :: xyz(3, nat) = reshape([& + & 0.293361019_wp, -0.373088547_wp, 1.367046476_wp, & + & 0.809123858_wp, -0.017139812_wp, -3.579971976_wp, & + & -1.007866311_wp, 1.702434996_wp, -4.786089410_wp, & + & -3.594201610_wp, 1.503843721_wp, -3.551059172_wp, & + & -3.581691626_wp, 1.904484469_wp, -0.794705250_wp, & + & -3.845327260_wp, 4.532847268_wp, 0.066234886_wp, & + & -1.387134261_wp, 6.009837986_wp, -0.039929904_wp, & + & 0.617713540_wp, 4.694267638_wp, 1.363985121_wp, & + & 2.563167262_wp, 0.242716370_wp, -4.339755095_wp, & + & 0.312371660_wp, -1.850891847_wp, -4.005424822_wp, & + & -1.233197204_wp, 1.292969229_wp, -6.818319337_wp, & + & -0.257342847_wp, 3.635246453_wp, -4.664750123_wp, & + & -4.335087472_wp, -0.411941308_wp, -3.883896562_wp, & + & -4.902761076_wp, 2.830998087_wp, -4.491575658_wp, & + & -5.077976444_wp, 0.912189498_wp, -0.087928937_wp, & + & -4.514195675_wp, 4.474511435_wp, 2.038484919_wp, & + & -5.303156159_wp, 5.549349625_wp, -1.030259560_wp, & + & -1.734484744_wp, 7.926530083_wp, 0.702788992_wp, & + & -0.760255551_wp, 6.262418725_wp, -2.003165945_wp, & + & 0.408445315_wp, 5.511706289_wp, 3.023655624_wp, & + & 2.390068387_wp, 3.861900154_wp, 0.699311897_wp, & + & 1.043714407_wp, 1.007374982_wp, 6.603419729_wp, & + & 1.353969575_wp, 1.450855812_wp, 8.769839079_wp, & + & -1.092393742_wp, 1.208328414_wp, 5.504411948_wp, & + & 2.834059443_wp, 0.293228739_wp, 5.140847618_wp, & + & -1.767725019_wp, -4.790964903_wp, -1.672785199_wp, & + & -2.152454276_wp, -3.966496472_wp, 0.367287089_wp, & + & 5.882980699_wp, -4.703376116_wp, -0.492065678_wp, & + & 3.976115176_wp, -2.948047698_wp, 0.283515548_wp, & + & 6.543704397_wp, -5.826024366_wp, 1.126257626_wp, & + & 7.482652411_wp, -3.703049811_wp, -1.363701662_wp, & + & 5.024856254_wp, -5.991696619_wp, -1.872529207_wp, & + & 4.542598253_wp, -1.981282921_wp, 1.732122586_wp ],& + & shape(xyz)) + real(wp), parameter :: charge = 1.0_wp + integer, parameter :: uhf = 1 + call init(mol, sym, xyz, chrg=charge, uhf=1) + end subroutine Ce_0a7745 + subroutine h2o(mol) type(TMolecule), intent(out) :: mol integer, parameter :: nat = 3 diff --git a/test/unit/test_cpx.f90 b/test/unit/test_cpx.f90 index ef5411ec0..532dc0556 100644 --- a/test/unit/test_cpx.f90 +++ b/test/unit/test_cpx.f90 @@ -70,8 +70,6 @@ subroutine test_cpx_solv(error) character(len=*), parameter :: cpxsolvent="water" - call skip_test(error, 'Excluding CPCM-X test.') - return if (.not.get_xtb_feature('cpcmx')) then call skip_test(error, 'CPCM-X libary not available.') return diff --git a/test/unit/test_gfnff.f90 b/test/unit/test_gfnff.f90 index 6f2c3a4a3..0afec9760 100644 --- a/test/unit/test_gfnff.f90 +++ b/test/unit/test_gfnff.f90 @@ -37,7 +37,8 @@ subroutine collect_gfnff(testsuite) new_unittest("scaleup", test_gfnff_scaleup), & new_unittest("pdb", test_gfnff_pdb), & new_unittest("sdf", test_gfnff_sdf), & - new_unittest("pbc", test_gfnff_pbc) & + new_unittest("pbc", test_gfnff_pbc), & + new_unittest("Ln_An", test_gfnff_LnAn_H) & ] end subroutine collect_gfnff @@ -702,14 +703,14 @@ subroutine test_gfnff_pbc(error) integer, allocatable :: attmp(:) character(len=symbolLength), allocatable :: symtmp(:) logical, parameter :: pbc(3) = [.true., .true., .true. ] - ! structure names from X23 and mcVOL22 benchmark - character(len=*), parameter :: pbc_strucs(2) = [& - & "x06_b", "mcv15"] + ! structures from X23, mcVOL22, and GFN-FF for Ln/An paper + character(len=*), parameter :: pbc_strucs(3) = [& + & "x06_b", "mcv15", "Th_15"] !> references for original GFN-FF - real(wp), parameter :: ref_energies(2) = & - & [-9.522300429916_wp, -11.059826732607_wp ] - real(wp), parameter :: ref_gnorms(2) = & - & [ 0.083513313043_wp, 0.083870455567_wp ] + real(wp), parameter :: ref_energies(3) = & + & [-9.522300429916_wp, -11.059826732607_wp, -46.265672914_wp ] + real(wp), parameter :: ref_gnorms(3) = & + & [ 0.083513313043_wp, 0.083870455567_wp, 1.6127503408_wp ] real(wp), parameter :: ref_snorm(2) = & & [ 1.410826672_wp, 0.419545156_wp ] !> references for mcGFN-FF @@ -774,6 +775,30 @@ subroutine test_gfnff_pbc(error) call check_(error, norm2(sigma), mcref_snorm(iMol), thr=thr) end do + ! check GFN-FF calculation on periodic system with actinide (Th) + ! load molecule "Th_15" from molstock + call getMolecule(mol, pbc_strucs(3)) + ! reset energy, gradient, sigma + energy=0.0_wp + if (allocated(gradient)) deallocate(gradient) + allocate(gradient(3, mol%n)) + sigma = 0.0_wp + ! setup new calculator + call delete_file('charges') + ! original angewChem2020_2 version is default + call newGFFCalculator(env, mol, calc, '.param_gfnff.xtb', .false.) + call env%check(exitRun) + call check_(error, .not.exitRun) + ! run single point calculation + call calc%singlepoint(env, mol, chk, 2, .false., energy, gradient, sigma, & + & hl_gap, res) + call env%check(exitRun) + call check_(error, .not.exitRun) + ! check energy and gradient versus reference calculation + call check_(error, energy, ref_energies(3), thr=thr) + call check_(error, norm2(gradient), ref_gnorms(3), thr=thr) + + ! check super cell scaling ! ! load molecule from molstock call getMolecule(mol, "mcv15") @@ -854,5 +879,106 @@ subroutine test_gfnff_pbc(error) end subroutine test_gfnff_pbc +subroutine test_gfnff_LnAn_H(error) + use xtb_mctc_accuracy, only : wp + use xtb_test_molstock, only : getMolecule + + use xtb_type_molecule + use xtb_type_param + use xtb_type_pcem + use xtb_type_data, only : scc_results + use xtb_type_environment, only : TEnvironment, init + use xtb_type_restart, only : TRestart + use xtb_mctc_symbols, only : symbolLength + + use xtb_gfnff_calculator, only : TGFFCalculator, newGFFCalculator + + type(error_type), allocatable, intent(out) :: error + + real(wp), parameter :: thr = 1.0e-8_wp ! threshold for energy + real(wp), parameter :: qthr = 1.0e-4_wp ! threshold for charge + + type(TEnvironment) :: env + type(TMolecule) :: mol + type(TRestart) :: chk + type(TGFFCalculator) :: calc + type(scc_results) :: res + + integer :: i + integer, parameter :: nat = 33 + logical :: exitRun + real(wp) :: energy, charges(nat), hl_gap, sigma(3,3) + real(wp), allocatable :: gradient(:, :) + real(wp), allocatable :: xyztmp(:,:), lattmp(:,:) + integer, allocatable :: attmp(:) + character(len=symbolLength), allocatable :: symtmp(:) + logical, parameter :: pbc(3) = [.true., .true., .true. ] + ! structures from X23, mcVOL22, and GFN-FF for Ln/An paper + character(len=*), parameter :: struc(1) = ["Ce_0a7745"] + !> references for original GFN-FF + real(wp), parameter :: ref_charges(nat) = [& + & 0.93714798_wp, -0.39685427_wp, -0.00104516_wp, -0.01815703_wp, -0.26122075_wp, & + & -0.01774628_wp, -0.00287629_wp, -0.36261674_wp, 0.19471344_wp, 0.20497869_wp, & + & 0.09972845_wp, 0.08148308_wp, 0.10674954_wp, 0.09486206_wp, 0.17232582_wp, & + & 0.11414769_wp, 0.09722964_wp, 0.10512679_wp, 0.08798849_wp, 0.22945410_wp, & + & 0.15310280_wp, 0.88629154_wp, -0.48033324_wp, -0.48478265_wp, -0.49522170_wp, & + & -0.14264511_wp, -0.12548232_wp, -0.01061975_wp, -0.40598965_wp, 0.13119126_wp, & + & 0.11822500_wp, 0.13303034_wp, 0.25781424_wp & + & ] + real(wp), parameter :: ref_energy = -3.550331926678 ! reference energy + real(wp), parameter :: ref_gnorm = 0.213662542141 ! reference gradient norm + ! number of neighbors and sum of neighbor indices for atom 1 (Ce) + integer, parameter :: ref_numnb_1 = 9 + integer, parameter :: ref_sumnb_1 = 177 ! 2+5+8+22+24+25+26+27+29 +9 + ! number of neighbors and sum of neighbor indices for atom 21 (H) + integer, parameter :: ref_numnb_21 = 1 + integer, parameter :: ref_sumnb_21 = 9 ! 8 +1 + ! number of neighbors and sum of neighbor indices for atom 33 (H) + integer, parameter :: ref_numnb_33 = 1 + integer, parameter :: ref_sumnb_33 = 30 ! 29 +1 + + + + call init(env) + + ! check GFN-FF calculation on periodic system with actinide (Th) + ! load molecule from molstock + call getMolecule(mol, struc(1)) + ! reset energy and gradient + energy=0.0_wp + if (allocated(gradient)) deallocate(gradient) + allocate(gradient(3, mol%n)) + ! setup new calculator + call delete_file('charges') + ! original angewChem2020_2 version is default + call newGFFCalculator(env, mol, calc, '.param_gfnff.xtb', .false.) + call env%check(exitRun) + call check_(error, .not.exitRun) + ! run single point calculation + call calc%singlepoint(env, mol, chk, 2, .false., energy, gradient, sigma, & + & hl_gap, res) + call env%check(exitRun) + call check_(error, .not.exitRun) + ! check energy and gradient versus reference calculation + call check_(error, energy, ref_energy, thr=thr) + call check_(error, norm2(gradient), ref_gnorm, thr=thr) + ! get charges from SP calculation + charges = chk%nlist%q ! these charges are in the gfnff_charges file + ! check charges of all atoms + do i=1, nat + ! compare to reference + call check_(error, charges(i), ref_charges(i), thr=qthr) + enddo + ! The Ln/An-H bonds in GFN-FF are charge dependent + ! But it is compared in the gfnff_ini where only topo%qa is available + ! Therefore we check the neighbor list here additionally + call check_(error, calc%neigh%nb(42,1,1), ref_numnb_1) ! check numnb of atom 1 + call check_(error, sum(calc%neigh%nb(:,1,1)), ref_sumnb_1) ! check sumnb of atom 1 + call check_(error, calc%neigh%nb(42,21,1), ref_numnb_21) ! check numnb of atom 21 + call check_(error, sum(calc%neigh%nb(:,21,1)), ref_sumnb_21) ! check sumnb of atom 21 + call check_(error, calc%neigh%nb(42,33,1), ref_numnb_33) ! check numnb of atom 33 + call check_(error, sum(calc%neigh%nb(:,33,1)), ref_sumnb_33) ! check sumnb of atom 33 + +end subroutine test_gfnff_LnAn_H end module test_gfnff