From 1779d8afafb75af831569c4d7ef54163d9bbd577 Mon Sep 17 00:00:00 2001 From: Thomas3R <39367840+Thomas3R@users.noreply.github.com> Date: Wed, 18 Sep 2024 11:15:03 +0200 Subject: [PATCH] GFN-FF reparameterization for lanthanides and extension to actinides (#1103) enables GFN-FF calculations for systems including elements up to lawrencium (103). performance has been improved for lanthanides. neighbor list can be manually adjusted via the detailed input file unit tests have been added to test Ln/An energies, gradients, and the Ln/An-H neighbor list setup. --- .param_gfnff.xtb | 17 ++ src/approxrab.f90 | 141 +++++++++------- src/geoopt_driver.f90 | 2 +- src/gfnff/gdisp0.f90 | 1 + src/gfnff/gfnff_ini.f90 | 146 ++++++++++++++-- src/gfnff/gfnff_ini2.f90 | 9 +- src/gfnff/gfnff_param.f90 | 213 ++++++++++++++--------- src/gfnff/gfnff_rab.f | 146 ++++++++++++---- src/gfnff/neighbor.f90 | 10 +- src/lbfgs_anc/driver.f90 | 6 +- src/model_hessian.f90 | 7 +- src/optimizer.f90 | 1 + src/pqn.f90 | 2 + src/set_module.f90 | 61 +++++++ src/setparam.f90 | 4 + src/type/anc.f90 | 2 +- src/type/molecule.f90 | 2 + test/unit/molstock.f90 | 346 ++++++++++++++++++++++++++++++++++++++ test/unit/test_gfnff.f90 | 142 +++++++++++++++- 19 files changed, 1051 insertions(+), 207 deletions(-) diff --git a/.param_gfnff.xtb b/.param_gfnff.xtb index 307cad2f3..bcec59f74 100644 --- a/.param_gfnff.xtb +++ b/.param_gfnff.xtb @@ -84,3 +84,20 @@ 84 1.236314 0.076708 0.019269 2.067331 0.192274 0.557520 0.322666 0.415042 0.307481 0.144762 0.767979 85 1.310129 0.000273 0.074803 2.228923 0.127706 0.563373 0.333641 1.259822 0.316447 0.231884 0.536799 86 1.157380 -0.068929 0.016657 1.874218 0.086756 0.484713 0.434163 0.400000 0.400000 0.400000 0.500000 +87 0.789115 0.485375 -0.024950 0.977079 0.150000 0.574626 0.302575 0.032198 0.000000 0.090741 0.000000 +88 0.798704 0.416264 -0.033006 0.770260 0.150000 0.560506 0.163290 0.036663 0.000000 0.076783 0.000000 +89 1.053384 -0.009993 0.065498 0.717658 0.370682 0.682723 0.187645 0.281449 0.078710 0.310896 0.122336 +90 1.056040 -0.010075 0.058290 0.732120 0.368511 0.684824 0.190821 0.280526 0.079266 0.309131 0.131176 +91 1.058772 -0.010138 0.052304 0.743405 0.366339 0.686925 0.193998 0.279603 0.079822 0.307367 0.140015 +92 1.061580 -0.010179 0.047540 0.751515 0.364168 0.689026 0.197174 0.278680 0.080379 0.305602 0.148855 +93 1.064463 -0.010201 0.043999 0.756447 0.361996 0.691127 0.200351 0.277757 0.080935 0.303838 0.157695 +94 1.067422 -0.010201 0.041680 0.758203 0.359825 0.693228 0.203527 0.276834 0.081491 0.302073 0.166534 +95 1.070456 -0.010182 0.040584 0.756783 0.357654 0.695329 0.206703 0.275911 0.082047 0.300309 0.175374 +96 1.073566 -0.010141 0.040710 0.752186 0.355482 0.697430 0.209880 0.274988 0.082603 0.298544 0.184214 +97 1.076751 -0.010080 0.042058 0.744413 0.353311 0.699531 0.213056 0.274065 0.083159 0.296779 0.193053 +98 1.080012 -0.009999 0.044628 0.733463 0.351139 0.701631 0.216233 0.273142 0.083716 0.295015 0.201893 +99 1.083349 -0.009897 0.048421 0.719337 0.348968 0.703732 0.219409 0.272219 0.084272 0.293250 0.210733 +100 1.086761 -0.009775 0.053436 0.702034 0.346797 0.705833 0.222585 0.271296 0.084828 0.291486 0.219572 +101 1.090249 -0.009632 0.059674 0.681554 0.344625 0.707934 0.225762 0.270373 0.085384 0.289721 0.228412 +102 1.093812 -0.009468 0.067134 0.657898 0.342454 0.710035 0.228938 0.269450 0.085940 0.287957 0.237252 +103 1.097451 -0.009284 0.075816 0.631066 0.340282 0.712136 0.232115 0.268528 0.086496 0.286192 0.246091 diff --git a/src/approxrab.f90 b/src/approxrab.f90 index 8ed1969b1..924d200c5 100644 --- a/src/approxrab.f90 +++ b/src/approxrab.f90 @@ -22,71 +22,86 @@ module xtb_approxrab public :: pbc_approx_rab, approx_rab, approx_bonds ! parameter blocks - real(wp),private, dimension(86) :: cnfak - real(wp),private, dimension(86) :: r0 - real(wp),private, dimension(86) :: en + real(wp),private, dimension(103) :: cnfak + real(wp),private, dimension(103) :: r0 + real(wp),private, dimension(103) :: en real(wp),private, dimension(4,2):: p ! START PARAMETER-------------------------------------------------- - data en /& - 2.30085633_wp, 2.78445145_wp, 1.52956084_wp, 1.51714704_wp, 2.20568300_wp,& - 2.49640820_wp, 2.81007174_wp, 4.51078438_wp, 4.67476223_wp, 3.29383610_wp,& - 2.84505365_wp, 2.20047950_wp, 2.31739628_wp, 2.03636974_wp, 1.97558064_wp,& - 2.13446570_wp, 2.91638164_wp, 1.54098156_wp, 2.91656301_wp, 2.26312147_wp,& - 2.25621439_wp, 1.32628677_wp, 2.27050569_wp, 1.86790977_wp, 2.44759456_wp,& - 2.49480042_wp, 2.91545568_wp, 3.25897750_wp, 2.68723778_wp, 1.86132251_wp,& - 2.01200832_wp, 1.97030722_wp, 1.95495427_wp, 2.68920990_wp, 2.84503857_wp,& - 2.61591858_wp, 2.64188286_wp, 2.28442252_wp, 1.33011187_wp, 1.19809388_wp,& - 1.89181390_wp, 2.40186898_wp, 1.89282464_wp, 3.09963488_wp, 2.50677823_wp,& - 2.61196704_wp, 2.09943450_wp, 2.66930105_wp, 1.78349472_wp, 2.09634533_wp,& - 2.00028974_wp, 1.99869908_wp, 2.59072029_wp, 2.54497829_wp, 2.52387890_wp,& - 2.30204667_wp, 1.60119300_wp, 2.00000000_wp, 2.00000000_wp, 2.00000000_wp,& - 2.00000000_wp, 2.00000000_wp, 2.00000000_wp, 2.00000000_wp, 2.00000000_wp,& - 2.00000000_wp, 2.00000000_wp, 2.00000000_wp, 2.00000000_wp, 2.00000000_wp,& - 2.00000000_wp, 2.30089349_wp, 1.75039077_wp, 1.51785130_wp, 2.62972945_wp,& - 2.75372921_wp, 2.62540906_wp, 2.55860939_wp, 3.32492356_wp, 2.65140898_wp,& - 1.52014458_wp, 2.54984804_wp, 1.72021963_wp, 2.69303422_wp, 1.81031095_wp,& - 2.34224386_wp& - / - data r0 /& - 0.55682207_wp, 0.80966997_wp, 2.49092101_wp, 1.91705642_wp, 1.35974851_wp,& - 0.98310699_wp, 0.98423007_wp, 0.76716063_wp, 1.06139799_wp, 1.17736822_wp,& - 2.85570926_wp, 2.56149012_wp, 2.31673425_wp, 2.03181740_wp, 1.82568535_wp,& - 1.73685958_wp, 1.97498207_wp, 2.00136196_wp, 3.58772537_wp, 2.68096221_wp,& - 2.23355957_wp, 2.33135502_wp, 2.15870365_wp, 2.10522128_wp, 2.16376162_wp,& - 2.10804037_wp, 1.96460045_wp, 2.00476257_wp, 2.22628712_wp, 2.43846700_wp,& - 2.39408483_wp, 2.24245792_wp, 2.05751204_wp, 2.15427677_wp, 2.27191920_wp,& - 2.19722638_wp, 3.80910350_wp, 3.26020971_wp, 2.99716916_wp, 2.71707818_wp,& - 2.34950167_wp, 2.11644818_wp, 2.47180659_wp, 2.32198800_wp, 2.32809515_wp,& - 2.15244869_wp, 2.55958313_wp, 2.59141300_wp, 2.62030465_wp, 2.39935278_wp,& - 2.56912355_wp, 2.54374096_wp, 2.56914830_wp, 2.53680807_wp, 4.24537037_wp,& - 3.66542289_wp, 3.19903011_wp, 2.80000000_wp, 2.80000000_wp, 2.80000000_wp,& - 2.80000000_wp, 2.80000000_wp, 2.80000000_wp, 2.80000000_wp, 2.80000000_wp,& - 2.80000000_wp, 2.80000000_wp, 2.80000000_wp, 2.80000000_wp, 2.80000000_wp,& - 2.80000000_wp, 2.34880037_wp, 2.37597108_wp, 2.49067697_wp, 2.14100577_wp,& - 2.33473532_wp, 2.19498900_wp, 2.12678348_wp, 2.34895048_wp, 2.33422774_wp,& - 2.86560827_wp, 2.62488837_wp, 2.88376127_wp, 2.75174124_wp, 2.83054552_wp,& - 2.63264944_wp& - / - data cnfak /& - 0.17957827_wp, 0.25584045_wp,-0.02485871_wp, 0.00374217_wp, 0.05646607_wp,& - 0.10514203_wp, 0.09753494_wp, 0.30470380_wp, 0.23261783_wp, 0.36752208_wp,& - 0.00131819_wp,-0.00368122_wp,-0.01364510_wp, 0.04265789_wp, 0.07583916_wp,& - 0.08973207_wp,-0.00589677_wp, 0.13689929_wp,-0.01861307_wp, 0.11061699_wp,& - 0.10201137_wp, 0.05426229_wp, 0.06014681_wp, 0.05667719_wp, 0.02992924_wp,& - 0.03764312_wp, 0.06140790_wp, 0.08563465_wp, 0.03707679_wp, 0.03053526_wp,& - -0.00843454_wp, 0.01887497_wp, 0.06876354_wp, 0.01370795_wp,-0.01129196_wp,& - 0.07226529_wp, 0.01005367_wp, 0.01541506_wp, 0.05301365_wp, 0.07066571_wp,& - 0.07637611_wp, 0.07873977_wp, 0.02997732_wp, 0.04745400_wp, 0.04582912_wp,& - 0.10557321_wp, 0.02167468_wp, 0.05463616_wp, 0.05370913_wp, 0.05985441_wp,& - 0.02793994_wp, 0.02922983_wp, 0.02220438_wp, 0.03340460_wp,-0.04110969_wp,& - -0.01987240_wp, 0.07260201_wp, 0.07700000_wp, 0.07700000_wp, 0.07700000_wp,& - 0.07700000_wp, 0.07700000_wp, 0.07700000_wp, 0.07700000_wp, 0.07700000_wp,& - 0.07700000_wp, 0.07700000_wp, 0.07700000_wp, 0.07700000_wp, 0.07700000_wp,& - 0.07700000_wp, 0.08379100_wp, 0.07314553_wp, 0.05318438_wp, 0.06799334_wp,& - 0.04671159_wp, 0.06758819_wp, 0.09488437_wp, 0.07556405_wp, 0.13384502_wp,& - 0.03203572_wp, 0.04235009_wp, 0.03153769_wp,-0.00152488_wp, 0.02714675_wp,& - 0.04800662_wp& - / + data en /& + 2.30085633, 2.78445145, 1.52956084, 1.51714704, 2.20568300,& + 2.49640820, 2.81007174, 4.51078438, 4.67476223, 3.29383610,& + 2.84505365, 2.20047950, 2.31739628, 2.03636974, 1.97558064,& + 2.13446570, 2.91638164, 1.54098156, 2.91656301, 2.26312147,& + 2.25621439, 1.32628677, 2.27050569, 1.86790977, 2.44759456,& + 2.49480042, 2.91545568, 3.25897750, 2.68723778, 1.86132251,& + 2.01200832, 1.97030722, 1.95495427, 2.68920990, 2.84503857,& + 2.61591858, 2.64188286, 2.28442252, 1.33011187, 1.19809388,& + 1.89181390, 2.40186898, 1.89282464, 3.09963488, 2.50677823,& + 2.61196704, 2.09943450, 2.66930105, 1.78349472, 2.09634533,& + 2.00028974, 1.99869908, 2.59072029, 2.54497829, 2.52387890,& + 2.30204667, 1.60119300, 2.00000000, 2.00000000, 2.00000000,& ! 60 + 2.00000000, 2.00000000, 2.00000000, 2.00000000, 2.00000000,& + 2.00000000, 2.00000000, 2.00000000, 2.00000000, 2.00000000,& ! 70 + 2.00000000, 2.30089349, 1.75039077, 1.51785130, 2.62972945,& + 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,& + 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 + 2.00000000& + / + data r0 /& + 0.55682207, 0.80966997, 2.49092101, 1.91705642, 1.35974851,& + 0.98310699, 0.98423007, 0.76716063, 1.06139799, 1.17736822,& ! 10 + 2.85570926, 2.56149012, 2.31673425, 2.03181740, 1.82568535,& + 1.73685958, 1.97498207, 2.00136196, 3.58772537, 2.68096221,& ! 20 + 2.23355957, 2.33135502, 2.15870365, 2.10522128, 2.16376162,& + 2.10804037, 1.96460045, 2.00476257, 2.22628712, 2.43846700,& ! 30 + 2.39408483, 2.24245792, 2.05751204, 2.15427677, 2.27191920,& + 2.19722638, 3.80910350, 3.26020971, 2.99716916, 2.71707818,& ! 40 + 2.34950167, 2.11644818, 2.47180659, 2.32198800, 2.32809515,& + 2.15244869, 2.55958313, 2.59141300, 2.62030465, 2.39935278,& ! 50 + 2.56912355, 2.54374096, 2.56914830, 2.53680807, 4.24537037,& + 3.66542289, 3.22480000, 3.21280000, 3.10550000, 3.10200000,& ! 60 + 3.10840000, 3.14030000, 3.06390000, 3.10730000, 3.10000000,& + 3.11910000, 3.10760000, 3.13740000, 3.09740000, 2.92860000,& + 3.05880000, 2.34880037, 2.37597108, 2.49067697, 2.14100577,& + 2.33473532, 2.19498900, 2.12678348, 2.34895048, 2.33422774,& + 2.86560827, 2.62488837, 2.88376127, 2.75174124, 2.83054552,& + 2.63264944,& + 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 + 4.20000000& + / + data cnfak /& + 0.17957827, 0.25584045,-0.02485871, 0.00374217, 0.05646607,& + 0.10514203, 0.09753494, 0.30470380, 0.23261783, 0.36752208,& + 0.00131819,-0.00368122,-0.01364510, 0.04265789, 0.07583916,& + 0.08973207,-0.00589677, 0.13689929,-0.01861307, 0.11061699,& + 0.10201137, 0.05426229, 0.06014681, 0.05667719, 0.02992924,& + 0.03764312, 0.06140790, 0.08563465, 0.03707679, 0.03053526,& + -0.00843454, 0.01887497, 0.06876354, 0.01370795,-0.01129196,& + 0.07226529, 0.01005367, 0.01541506, 0.05301365, 0.07066571,& + 0.07637611, 0.07873977, 0.02997732, 0.04745400, 0.04582912,& + 0.10557321, 0.02167468, 0.05463616, 0.05370913, 0.05985441,& ! 60 + 0.02793994, 0.02922983, 0.02220438, 0.03340460,-0.04110969,& + -0.01987240, 0.07260201, 0.07700000, 0.07700000, 0.07700000,& ! 70 + 0.07700000, 0.07700000, 0.07700000, 0.07700000, 0.07700000,& + 0.07700000, 0.07700000, 0.07700000, 0.07700000, 0.07700000,& ! 80 + 0.07700000, 0.08379100, 0.07314553, 0.05318438, 0.06799334,& + 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,& + 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 + 0.07700000& + / ! END PARAMETER----------------------------------------------------- ! global EN polynomial parameter (NOTE: x 10^3) 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 a797a223b..2335b00d6 100644 --- a/src/gfnff/gfnff_ini.f90 +++ b/src/gfnff/gfnff_ini.f90 @@ -55,7 +55,7 @@ subroutine gfnff_ini(env,pr,makeneighbor,mol,gen,param,topo,neigh,efield,accurac integer ineig,jneig,nrot,bbtyp,ringtyp,nn1,nn2,hybi,hybj,pis,ka,nh,jdum,hcalc,nc integer ringsi,ringsj,ringsk,ringl,npi,nelpi,picount,npiall,maxtors,rings4,nheav integer nm,maxhb,ki,n13,current,ncarbo,mtyp1,mtyp2 - integer ind3(3),sr(20),cr(10,20),niel(86) + integer ind3(3),sr(20),cr(10,20),niel(103) integer qloop_count,nf,nsi,nmet,nhi,nhj,ifrag integer hbA,hbH,Bat,atB,Aat,Hat integer AHB_nr @@ -113,7 +113,7 @@ integer function itabrow6(i) character(len=255) atmp integer :: ich, err real(wp) :: dispthr, cnthr, repthr, hbthr1, hbthr2 - logical :: exitRun, nb_call + logical :: exitRun, nb_call, adjLnAn call gfnff_thresholds(accuracy, dispthr, cnthr, repthr, hbthr1, hbthr2) @@ -167,7 +167,7 @@ integer function itabrow6(i) enddo write(env%unit,'(10x,"Pauling EN used:")') - do i=1,86 + do i=1,103 if(niel(i).gt.0) write(env%unit,'(10x,"Z :",i2," EN :",f6.2)') i,param%en(i) enddo @@ -262,6 +262,14 @@ integer function itabrow6(i) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! neighbor list, hyb and ring info !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! check if Ln or An in coordinates + adjLnAn=.false. + do i=1, mol%n + if ((mol%at(i).ge.57.and.mol%at(i).le.71).or.(mol%at(i).ge.89.and.mol%at(i).le.103)) then + adjLnAn=.true. + exit + endif + enddo topo%qa = 0 qloop_count = 0 @@ -269,7 +277,7 @@ integer function itabrow6(i) !111 continue ! do the loop only if factor is significant - do while (qloop_count.lt.2.and.gen%rqshrink.gt.1.d-3) + do while ((qloop_count.lt.2.and.gen%rqshrink.gt.1.d-3).or.adjLnAn) write(env%unit,'(10x,"----------------------------------------")') write(env%unit,'(10x,"generating topology and atomic info file ...")') @@ -277,6 +285,14 @@ integer function itabrow6(i) & gen%rthr,gen%rthr2,gen%linthr,mchar,topo%hyb,itag,param,topo,mol,neigh,nb_call) nb_call = .true. + ! special treatment for hydrogen bound to Ln or An + if (adjLnAn.and.allocated(topo%qa).and.qloop_count.ne.0) then + call adjust_NB_LnH_AnH(param, mol, topo, neigh) + adjLnAn=.false. + endif + !> gfnff_topo adjustments + call gfnff_topo_changes(env, neigh) + do i=1,mol%n imetal(i)=param%metal(mol%at(i)) ! Sn,Pb,Bi, with small CN are better described as non-metals @@ -644,7 +660,7 @@ integer function itabrow6(i) enddo endif qloop_count=qloop_count+1 - if(qloop_count.lt.2.and.gen%rqshrink.gt.1.d-3) then ! do the loop only if factor is significant + if((qloop_count.lt.2.and.gen%rqshrink.gt.1.d-3).or.adjLnAn) then ! do the loop only if factor is significant deallocate(btyp,pibo,pimvec,neigh%blist) ! goto 111 endif @@ -2175,7 +2191,6 @@ integer function itabrow6(i) call specialTorsList(nn, mol, topo, neigh, topo%sTorsl) endif - if(pr)then write(env%unit,*) write(env%unit,*) 'GFN-FF setup done.' @@ -2184,6 +2199,80 @@ integer function itabrow6(i) contains +! remove hydrogen bonds from topology if bound to Ln or An and topo%qa(H)>-0.0281 +subroutine adjust_NB_LnH_AnH(param, mol, topo, neigh) + type(TGFFData), intent(in) :: param + type(TMolecule), intent(in) :: mol ! # molecule type + type(TGFFTopology), intent(in) :: topo + 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 + ! only apply changes for Ln and An + if ((mol%at(i).ge.57.and.mol%at(i).le.71).or.(mol%at(i).ge.89.and.mol%at(i).le.103)) then + ! loop over all cells + do iTr=1, neigh%numctr + ! loop over all neighbors + nnb = neigh%nb(neigh%numnb, i, iTr) + idx=0 + do count_idx=1, nnb + idx = idx + 1 + inb=neigh%nb(idx,i,iTr) ! atom index of neighbor of i + if (inb.eq.0) cycle + ! 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.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 + nb_tmp(idx) = 0 + ! reset neigh%nb + neigh%nb(1:neigh%numnb-2,i,iTr)=0 + ! update neigh%nb with copied nb list + neigh%nb(neigh%numnb,i,iTr) = nb_tmp(neigh%numnb) ! number neighbors + l=0 + do k=1, neigh%numnb-2 + if (nb_tmp(k).ne.0) then + l = l + 1 + neigh%nb(l,i,iTr) = nb_tmp(k) + end if + enddo + ! setup copy of neighbor list for this H + do iTrH=1, neigh%numctr + do k=1, neigh%nb(neigh%numnb,inb,iTrH) + if (neigh%nb(k,inb,iTrH).eq.i) then + nb_tmp = neigh%nb(:,inb,iTrH) + nb_tmp(neigh%numnb) = neigh%nb(neigh%numnb,inb,iTrH) - 1 + nb_tmp(k) = 0 + ! reset neigh%nb + neigh%nb(1:neigh%numnb-2,inb,iTrH)=0 + ! update neigh%nb with copied nb list + neigh%nb(neigh%numnb,inb,iTrH) = nb_tmp(neigh%numnb) ! number neighbors + l=0 + do m=1, neigh%numnb-2 + if (nb_tmp(m).ne.0) then + l = l + 1 + neigh%nb(l,inb,iTrH) = nb_tmp(m) + end if + enddo + endif + enddo + enddo + ! since the current element was removed from neigh%nb, idx should be reduced by one + idx = idx - 1 + end if ! if topo%qa < -0.0282 + end if + enddo + enddo + end if + enddo + +end subroutine adjust_NB_LnH_AnH + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! special treatment for rotation around carbon triple bonds !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @@ -2298,14 +2387,16 @@ pure elemental function zeta(at,q) real(wp),intent(in) :: q real(wp) :: zeta,qmod - real(wp),parameter :: zeff(86) = (/ & + real(wp),parameter :: zeff(103) = (/ & & 1, 2, & ! H-He & 3, 4, 5, 6, 7, 8, 9,10, & ! Li-Ne & 11,12, 13,14,15,16,17,18, & ! Na-Ar & 19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36, & ! K-Kr & 9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26, & ! Rb-Xe & 9,10,11,30,31,32,33,34,35,36,37,38,39,40,41,42,43, & ! Cs-Lu - & 12,13,14,15,16,17,18,19,20,21,22,23,24,25,26/) ! Hf-Rn + & 12,13,14,15,16,17,18,19,20,21,22,23,24,25,26, & ! Hf-Rn + & 9,10,11,30,31,32,33,34,35,36,37,38,39,40,41,42,43 & ! Fr-Lr + &/) !! Semiempirical Evaluation of the GlobalHardness of the Atoms of 103 !! Elements of the Periodic Table Using the Most Probable Radii as !! their Size Descriptors DULAL C. GHOSH, NAZMUL ISLAM 2009 in @@ -2314,7 +2405,7 @@ pure elemental function zeta(at,q) !! values in the paper multiplied by two because !! (ii:ii)=(IP-EA)=d^2 E/dN^2 but the hardness !! definition they use is 1/2d^2 E/dN^2 (in Eh) - real(wp),parameter :: c(1:86) = (/ & + real(wp),parameter :: c(1:103) = (/ & &0.47259288_wp,0.92203391_wp,0.17452888_wp,0.25700733_wp,0.33949086_wp,0.42195412_wp, & ! H-C &0.50438193_wp,0.58691863_wp,0.66931351_wp,0.75191607_wp,0.17964105_wp,0.22157276_wp, & ! N-Mg &0.26348578_wp,0.30539645_wp,0.34734014_wp,0.38924725_wp,0.43115670_wp,0.47308269_wp, & ! Al-Ar @@ -2329,7 +2420,11 @@ pure elemental function zeta(at,q) &0.25932503_wp,0.27676094_wp,0.29418231_wp,0.31159587_wp,0.32902274_wp,0.34592298_wp, & ! Ho-Hf &0.36388048_wp,0.38130586_wp,0.39877476_wp,0.41614298_wp,0.43364510_wp,0.45104014_wp, & ! Ta-Pt &0.46848986_wp,0.48584550_wp,0.12526730_wp,0.14268677_wp,0.16011615_wp,0.17755889_wp, & ! Au-Po - &0.19497557_wp,0.21240778_wp/) + &0.19497557_wp,0.21240778_wp,& ! At, Rn + &0.07263133_wp,0.09421788_wp,0.09920108_wp,0.10418429_wp,0.14235212_wp,0.16393866_wp, & ! Fr-U + &0.18675998_wp,0.22370039_wp,0.25113742_wp,0.25026279_wp,0.28843797_wp,0.31002451_wp, & ! Np-Cf + &0.33159636_wp,0.35316820_wp,0.36822807_wp,0.39634864_wp,0.40135389_wp & ! Es-Lr + &/) intrinsic :: exp @@ -2345,4 +2440,35 @@ end function zeta end subroutine gfnff_ini +subroutine gfnff_topo_changes(env, neigh) + use xtb_type_environment, only : TEnvironment + use xtb_gfnff_neighbor, only : TNeigh + use xtb_setparam + type(TEnvironment), intent(inout) :: env + type(TNeigh), intent(inout) :: neigh ! main type for introducing PBC + character(len=*), parameter :: source = 'gfnff_ini' + integer :: int_tmp(40) + integer :: i,j,idx,iTr,d1,d2,numnb,nnb + + ! 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.-1) then + d2=size(set%ffnb, dim=2) + do i=1, d2 + 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 + neigh%nb(neigh%numnb,idx,1) = set%ffnb(42,i) + end do + write(env%unit,*) '' + write(env%unit,*) 'The neighborlist has been adjusted according to the input file.' + write(env%unit,*) '' + end if + + +end subroutine gfnff_topo_changes + + end module xtb_gfnff_ini diff --git a/src/gfnff/gfnff_ini2.f90 b/src/gfnff/gfnff_ini2.f90 index e0a70da3e..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) @@ -58,9 +59,9 @@ subroutine gfnff_neigh(env,makeneighbor,natoms,at,xyz,rab,fq,f_in,f2_in,lintr, & real*8 ,allocatable :: cn(:),rtmp(:) integer iat,i,j,k,ni,ii,jj,kk,ll,lin,ati,nb20i,nbdiff,hc_crit,nbmdiff,nnf,nni,nh,nm integer ai,aj,nn,im,ncm,l,no, iTr, iTr2, numnbf, numnbm, numnb, idx, idxdum, idxdum2, numctr - real*8 r,pi,a1,f,f1,phi,f2,rco,fat(86) + real*8 r,pi,a1,f,f1,phi,f2,rco,fat(103) data pi/3.1415926535897932384626433832795029d0/ - data fat / 86 * 1.0d0 / + data fat / 103 * 1.0d0 / ! special hacks fat( 1)=1.02 @@ -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/gfnff_param.f90 b/src/gfnff/gfnff_param.f90 index 3bec740df..d0098523a 100644 --- a/src/gfnff/gfnff_param.f90 +++ b/src/gfnff/gfnff_param.f90 @@ -62,7 +62,7 @@ module xtb_gfnff_param logical :: gff_print = .true. !shows timing of energy + gradient logical :: make_chrg = .true. !generates new eeq chareges based on current geometry - real(wp), parameter :: chi_angewChem2020(86) = [& + real(wp), parameter :: chi_angewChem2020(103) = [& & 1.227054_wp, 1.451412_wp, 0.813363_wp, 1.062841_wp, 1.186499_wp, & & 1.311555_wp, 1.528485_wp, 1.691201_wp, 1.456784_wp, 1.231037_wp, & & 0.772989_wp, 1.199092_wp, 1.221576_wp, 1.245964_wp, 1.248942_wp, & @@ -74,15 +74,19 @@ module xtb_gfnff_param & 0.967863_wp, 1.002901_wp, 1.037940_wp, 1.072978_wp, 1.108017_wp, & & 1.143055_wp, 1.178094_wp, 1.213132_wp, 1.205076_wp, 1.075529_wp, & & 1.206919_wp, 1.303658_wp, 1.332656_wp, 1.179317_wp, 0.789115_wp, & - & 0.798704_wp, 1.127797_wp, 1.127863_wp, 1.127928_wp, 1.127994_wp, & - & 1.128059_wp, 1.128125_wp, 1.128190_wp, 1.128256_wp, 1.128322_wp, & - & 1.128387_wp, 1.128453_wp, 1.128518_wp, 1.128584_wp, 1.128649_wp, & - & 1.128715_wp, 1.128780_wp, 1.129764_wp, 1.130747_wp, 1.131731_wp, & - & 1.132714_wp, 1.133698_wp, 1.134681_wp, 1.135665_wp, 1.136648_wp, & - & 1.061832_wp, 1.053084_wp, 1.207830_wp, 1.236314_wp, 1.310129_wp, & - & 1.157380_wp] - - real(wp), parameter :: gam_angewChem2020(86) = [& + & 0.798704_wp, 0.993208_wp, 0.907847_wp, 0.836114_wp, 0.778008_wp, & ! 60 + & 0.733529_wp, 0.702678_wp, 0.685455_wp, 0.681858_wp, 0.691889_wp, & + & 0.715548_wp, 0.752834_wp, 0.803747_wp, 0.868288_wp, 0.946457_wp, & ! 70 + & 1.038252_wp, 1.128780_wp, 1.129764_wp, 1.130747_wp, 1.131731_wp, & + & 1.132714_wp, 1.133698_wp, 1.134681_wp, 1.135665_wp, 1.136648_wp, & ! 80 + & 1.061832_wp, 1.053084_wp, 1.207830_wp, 1.236314_wp, 1.310129_wp, & + & 1.157380_wp, 0.789115_wp, 0.798704_wp, 1.053384_wp, 1.056040_wp, & ! 90 + & 1.058772_wp, 1.061580_wp, 1.064463_wp, 1.067422_wp, 1.070456_wp, & + & 1.073566_wp, 1.076751_wp, 1.080012_wp, 1.083349_wp, 1.086761_wp, & ! 100 + & 1.090249_wp, 1.093812_wp, 1.097451_wp & ! 103 Lr + &] + + real(wp), parameter :: gam_angewChem2020(103) = [& &-0.448428_wp, 0.131022_wp, 0.571431_wp, 0.334622_wp,-0.089208_wp, & &-0.025895_wp,-0.027280_wp,-0.031236_wp,-0.159892_wp, 0.074198_wp, & & 0.316829_wp, 0.326072_wp, 0.069748_wp,-0.120184_wp,-0.193159_wp, & @@ -94,15 +98,19 @@ module xtb_gfnff_param &-0.010869_wp,-0.000951_wp, 0.008967_wp, 0.018884_wp, 0.028802_wp, & & 0.038720_wp, 0.048638_wp, 0.058556_wp, 0.036488_wp, 0.077711_wp, & & 0.077025_wp, 0.004547_wp, 0.039909_wp, 0.082630_wp, 0.485375_wp, & - & 0.416264_wp,-0.011212_wp,-0.011046_wp,-0.010879_wp,-0.010713_wp, & - &-0.010546_wp,-0.010380_wp,-0.010214_wp,-0.010047_wp,-0.009881_wp, & - &-0.009714_wp,-0.009548_wp,-0.009382_wp,-0.009215_wp,-0.009049_wp, & - &-0.008883_wp,-0.008716_wp,-0.006220_wp,-0.003724_wp,-0.001228_wp, & - & 0.001267_wp, 0.003763_wp, 0.006259_wp, 0.008755_wp, 0.011251_wp, & + & 0.416264_wp,-0.007277_wp,-0.007862_wp,-0.008388_wp,-0.008855_wp, & ! 60 + &-0.009263_wp,-0.009611_wp,-0.009901_wp,-0.010132_wp,-0.010304_wp, & + &-0.010417_wp,-0.010471_wp,-0.010465_wp,-0.010401_wp,-0.010278_wp, & ! 70 + &-0.010096_wp,-0.008716_wp,-0.006220_wp,-0.003724_wp,-0.001228_wp, & + & 0.001267_wp, 0.003763_wp, 0.006259_wp, 0.008755_wp, 0.011251_wp, & ! 80 & 0.020477_wp,-0.056566_wp, 0.051943_wp, 0.076708_wp, 0.000273_wp, & - &-0.068929_wp] + &-0.068929_wp, 0.485375_wp, 0.416264_wp,-0.009993_wp,-0.010075_wp, & ! 90 + &-0.010138_wp,-0.010179_wp,-0.010201_wp,-0.010201_wp,-0.010182_wp, & + &-0.010141_wp,-0.010080_wp,-0.009999_wp,-0.009897_wp,-0.009775_wp, & ! 100 + &-0.009632_wp,-0.009468_wp,-0.009284_wp & ! 103 + &] - real(wp), parameter :: cnf_angewChem2020(86) = [& + real(wp), parameter :: cnf_angewChem2020(103) = [& & 0.008904_wp, 0.004641_wp, 0.048324_wp, 0.080316_wp,-0.051990_wp, & & 0.031779_wp, 0.132184_wp, 0.157353_wp, 0.064120_wp, 0.036540_wp, & &-0.000627_wp, 0.005412_wp, 0.018809_wp, 0.016329_wp, 0.012149_wp, & @@ -114,15 +122,19 @@ module xtb_gfnff_param & 0.004992_wp, 0.009186_wp, 0.013379_wp, 0.017573_wp, 0.021766_wp, & & 0.025960_wp, 0.030153_wp, 0.034347_wp,-0.000052_wp,-0.039776_wp, & & 0.006661_wp, 0.050424_wp, 0.068985_wp, 0.023470_wp,-0.024950_wp, & - &-0.033006_wp, 0.058973_wp, 0.058595_wp, 0.058217_wp, 0.057838_wp, & - & 0.057460_wp, 0.057082_wp, 0.056704_wp, 0.056326_wp, 0.055948_wp, & - & 0.055569_wp, 0.055191_wp, 0.054813_wp, 0.054435_wp, 0.054057_wp, & - & 0.053679_wp, 0.053300_wp, 0.047628_wp, 0.041955_wp, 0.036282_wp, & - & 0.030610_wp, 0.024937_wp, 0.019264_wp, 0.013592_wp, 0.007919_wp, & + &-0.033006_wp, 0.062573_wp, 0.057117_wp, 0.052085_wp, 0.047480_wp, & ! 60 + & 0.043299_wp, 0.039544_wp, 0.036215_wp, 0.033311_wp, 0.030833_wp, & + & 0.028780_wp, 0.027153_wp, 0.025951_wp, 0.025174_wp, 0.024823_wp, & ! 70 + & 0.024898_wp, 0.053300_wp, 0.047628_wp, 0.041955_wp, 0.036282_wp, & + & 0.030610_wp, 0.024937_wp, 0.019264_wp, 0.013592_wp, 0.007919_wp, & ! 80 & 0.006383_wp,-0.089155_wp,-0.001293_wp, 0.019269_wp, 0.074803_wp, & - & 0.016657_wp] + & 0.016657_wp,-0.024950_wp,-0.033006_wp, 0.065498_wp, 0.058290_wp, & ! 90 + & 0.052304_wp, 0.047540_wp, 0.043999_wp, 0.041680_wp, 0.040584_wp, & + & 0.040710_wp, 0.042058_wp, 0.044628_wp, 0.048421_wp, 0.053436_wp, & ! 100 + & 0.059674_wp, 0.067134_wp, 0.075816_wp & ! 103 + &] - real(wp), parameter :: alp_angewChem2020(86) = [& + real(wp), parameter :: alp_angewChem2020(103) = [& & 0.585069_wp, 0.432382_wp, 0.628636_wp, 0.743646_wp, 1.167323_wp, & & 0.903430_wp, 1.278388_wp, 0.905347_wp, 1.067014_wp, 2.941513_wp, & & 0.687680_wp, 0.792170_wp, 1.337040_wp, 1.251409_wp, 1.068295_wp, & @@ -134,35 +146,42 @@ module xtb_gfnff_param & 0.684938_wp, 0.701032_wp, 0.717125_wp, 0.733218_wp, 0.749311_wp, & & 0.765405_wp, 0.781498_wp, 0.797591_wp, 1.296844_wp, 1.534068_wp, & & 1.727781_wp, 1.926871_wp, 2.175548_wp, 2.177702_wp, 0.977079_wp, & - & 0.770260_wp, 0.757372_wp, 0.757352_wp, 0.757332_wp, 0.757313_wp, & - & 0.757293_wp, 0.757273_wp, 0.757253_wp, 0.757233_wp, 0.757213_wp, & - & 0.757194_wp, 0.757174_wp, 0.757154_wp, 0.757134_wp, 0.757114_wp, & - & 0.757095_wp, 0.757075_wp, 0.756778_wp, 0.756480_wp, 0.756183_wp, & - & 0.755886_wp, 0.755589_wp, 0.755291_wp, 0.754994_wp, 0.754697_wp, & + & 0.770260_wp, 0.695335_wp, 0.653051_wp, 0.617284_wp, 0.588034_wp, & ! 60 + & 0.565301_wp, 0.549084_wp, 0.539383_wp, 0.536199_wp, 0.539532_wp, & + & 0.549382_wp, 0.565747_wp, 0.588630_wp, 0.618029_wp, 0.653945_wp, & ! 70 + & 0.696377_wp, 0.757075_wp, 0.756778_wp, 0.756480_wp, 0.756183_wp, & + & 0.755886_wp, 0.755589_wp, 0.755291_wp, 0.754994_wp, 0.754697_wp, & ! 80 & 0.868029_wp, 1.684375_wp, 2.001040_wp, 2.067331_wp, 2.228923_wp, & - & 1.874218_wp] + & 1.874218_wp, 0.977079_wp, 0.770260_wp, 0.717658_wp, 0.732120_wp, & ! 90 + & 0.743405_wp, 0.751515_wp, 0.756447_wp, 0.758203_wp, 0.756783_wp, & + & 0.752186_wp, 0.744413_wp, 0.733463_wp, 0.719337_wp, 0.702034_wp, & ! 100 + & 0.681554_wp, 0.657898_wp, 0.631066_wp & ! 103 + &] - real(wp), parameter :: bond_angewChem2020(86) = [& + real(wp), parameter :: bond_angewChem2020(103) = [& & 0.417997_wp, 0.258490_wp, 0.113608_wp, 0.195935_wp, 0.231217_wp, & - & 0.385248_wp, 0.379257_wp, 0.339249_wp, 0.330706_wp, 0.120319_wp, & + & 0.385248_wp, 0.379257_wp, 0.339249_wp, 0.330706_wp, 0.120319_wp, & ! 10 & 0.127255_wp, 0.173647_wp, 0.183796_wp, 0.273055_wp, 0.249044_wp, & - & 0.290653_wp, 0.218744_wp, 0.034706_wp, 0.136353_wp, 0.192467_wp, & + & 0.290653_wp, 0.218744_wp, 0.034706_wp, 0.136353_wp, 0.192467_wp, & ! 20 & 0.335860_wp, 0.314452_wp, 0.293044_wp, 0.271636_wp, 0.250228_wp, & - & 0.228819_wp, 0.207411_wp, 0.186003_wp, 0.164595_wp, 0.143187_wp, & + & 0.228819_wp, 0.207411_wp, 0.186003_wp, 0.164595_wp, 0.143187_wp, & ! 30 & 0.212434_wp, 0.210451_wp, 0.219870_wp, 0.224618_wp, 0.272206_wp, & - & 0.147864_wp, 0.150000_wp, 0.150000_wp, 0.329501_wp, 0.309632_wp, & + & 0.147864_wp, 0.150000_wp, 0.150000_wp, 0.329501_wp, 0.309632_wp, & ! 40 & 0.289763_wp, 0.269894_wp, 0.250025_wp, 0.230155_wp, 0.210286_wp, & - & 0.190417_wp, 0.170548_wp, 0.150679_wp, 0.192977_wp, 0.173411_wp, & + & 0.190417_wp, 0.170548_wp, 0.150679_wp, 0.192977_wp, 0.173411_wp, & ! 50 & 0.186907_wp, 0.192891_wp, 0.223202_wp, 0.172577_wp, 0.150000_wp, & - & 0.150000_wp, 0.370682_wp, 0.368511_wp, 0.366339_wp, 0.364168_wp, & + & 0.150000_wp, 0.370682_wp, 0.368511_wp, 0.366339_wp, 0.364168_wp, & ! 60 & 0.361996_wp, 0.359825_wp, 0.357654_wp, 0.355482_wp, 0.353311_wp, & - & 0.351139_wp, 0.348968_wp, 0.346797_wp, 0.344625_wp, 0.342454_wp, & + & 0.351139_wp, 0.348968_wp, 0.346797_wp, 0.344625_wp, 0.342454_wp, & ! 70 & 0.340282_wp, 0.338111_wp, 0.305540_wp, 0.272969_wp, 0.240398_wp, & - & 0.207828_wp, 0.175257_wp, 0.142686_wp, 0.110115_wp, 0.077544_wp, & + & 0.207828_wp, 0.175257_wp, 0.142686_wp, 0.110115_wp, 0.077544_wp, & ! 80 & 0.108597_wp, 0.148422_wp, 0.183731_wp, 0.192274_wp, 0.127706_wp, & - & 0.086756_wp] + & 0.086756_wp, 0.150000_wp, 0.150000_wp, 0.370682_wp, 0.368511_wp, & ! 90 + & 0.366339_wp, 0.364168_wp, 0.361996_wp, 0.359825_wp, 0.357654_wp, & + & 0.355482_wp, 0.353311_wp, 0.351139_wp, 0.348968_wp, 0.346797_wp, & ! 100 + & 0.344625_wp, 0.342454_wp, 0.340282_wp] - real(wp), parameter :: repa_angewChem2020(86) = [& + real(wp), parameter :: repa_angewChem2020(103) = [& & 2.639785_wp, 3.575012_wp, 0.732142_wp, 1.159621_wp, 1.561585_wp, & & 1.762895_wp, 2.173015_wp, 2.262269_wp, 2.511112_wp, 3.577220_wp, & & 0.338845_wp, 0.693023_wp, 0.678792_wp, 0.804784_wp, 1.012178_wp, & @@ -174,15 +193,19 @@ module xtb_gfnff_param & 0.768030_wp, 0.840055_wp, 0.912081_wp, 0.984106_wp, 1.056131_wp, & & 1.128156_wp, 1.200181_wp, 1.272206_wp, 0.478807_wp, 0.479759_wp, & & 0.579840_wp, 0.595241_wp, 0.644458_wp, 0.655289_wp, 0.574626_wp, & - & 0.560506_wp, 0.682723_wp, 0.684824_wp, 0.686925_wp, 0.689026_wp, & + & 0.560506_wp, 0.682723_wp, 0.684824_wp, 0.686925_wp, 0.689026_wp, & ! 60 & 0.691127_wp, 0.693228_wp, 0.695329_wp, 0.697430_wp, 0.699531_wp, & - & 0.701631_wp, 0.703732_wp, 0.705833_wp, 0.707934_wp, 0.710035_wp, & + & 0.701631_wp, 0.703732_wp, 0.705833_wp, 0.707934_wp, 0.710035_wp, & ! 70 & 0.712136_wp, 0.714237_wp, 0.745751_wp, 0.777265_wp, 0.808779_wp, & - & 0.840294_wp, 0.871808_wp, 0.903322_wp, 0.934836_wp, 0.966350_wp, & + & 0.840294_wp, 0.871808_wp, 0.903322_wp, 0.934836_wp, 0.966350_wp, & ! 80 & 0.467729_wp, 0.486102_wp, 0.559176_wp, 0.557520_wp, 0.563373_wp, & - & 0.484713_wp] + & 0.484713_wp, 0.574626_wp, 0.560506_wp, 0.682723_wp, 0.684824_wp, & ! 90 + & 0.686925_wp, 0.689026_wp, 0.691127_wp, 0.693228_wp, 0.695329_wp, & + & 0.697430_wp, 0.699531_wp, 0.701631_wp, 0.703732_wp, 0.705833_wp, & + & 0.707934_wp, 0.710035_wp, 0.712136_wp & + &] - real(wp), parameter :: repan_angewChem2020(86) = [& + real(wp), parameter :: repan_angewChem2020(103) = [& & 1.071395_wp, 1.072699_wp, 1.416847_wp, 1.156187_wp, 0.682382_wp, & & 0.556380_wp, 0.746785_wp, 0.847242_wp, 0.997252_wp, 0.873051_wp, & & 0.322503_wp, 0.415554_wp, 0.423946_wp, 0.415776_wp, 0.486773_wp, & @@ -200,9 +223,13 @@ module xtb_gfnff_param & 0.232115_wp, 0.235291_wp, 0.282937_wp, 0.330583_wp, 0.378229_wp, & & 0.425876_wp, 0.473522_wp, 0.521168_wp, 0.568814_wp, 0.616460_wp, & & 0.242521_wp, 0.293680_wp, 0.320931_wp, 0.322666_wp, 0.333641_wp, & - & 0.434163_wp] + & 0.434163_wp, 0.302575_wp, 0.163290_wp, 0.187645_wp, 0.190821_wp, & + & 0.193998_wp, 0.197174_wp, 0.200351_wp, 0.203527_wp, 0.206703_wp, & + & 0.209880_wp, 0.213056_wp, 0.216233_wp, 0.219409_wp, 0.222585_wp, & + & 0.225762_wp, 0.228938_wp, 0.232115_wp & + &] - real(wp), parameter :: angl_angewChem2020(86) = [& + real(wp), parameter :: angl_angewChem2020(103) = [& & 1.661808_wp, 0.300000_wp, 0.018158_wp, 0.029224_wp, 0.572683_wp, & & 0.771055_wp, 1.053577_wp, 2.159889_wp, 1.525582_wp, 0.400000_wp, & & 0.041070_wp, 0.028889_wp, 0.086910_wp, 0.494456_wp, 0.409204_wp, & @@ -220,9 +247,13 @@ module xtb_gfnff_param & 0.268528_wp, 0.267605_wp, 0.253760_wp, 0.239916_wp, 0.226071_wp, & & 0.212227_wp, 0.198382_wp, 0.184538_wp, 0.170693_wp, 0.156849_wp, & & 0.104547_wp, 0.313474_wp, 0.220185_wp, 0.415042_wp, 1.259822_wp, & - & 0.400000_wp] + & 0.400000_wp, 0.032198_wp, 0.036663_wp, 0.281449_wp, 0.280526_wp, & + & 0.279603_wp, 0.278680_wp, 0.277757_wp, 0.276834_wp, 0.275911_wp, & + & 0.274988_wp, 0.274065_wp, 0.273142_wp, 0.272219_wp, 0.271296_wp, & + & 0.270373_wp, 0.269450_wp, 0.268528_wp & + &] - real(wp), parameter :: angl2_angewChem2020(86) = [& + real(wp), parameter :: angl2_angewChem2020(103) = [& & 0.624197_wp, 0.600000_wp, 0.050000_wp, 0.101579_wp, 0.180347_wp, & & 0.755851_wp, 0.761551_wp, 0.813653_wp, 0.791274_wp, 0.400000_wp, & & 0.000000_wp, 0.022706_wp, 0.100000_wp, 0.338514_wp, 0.453023_wp, & @@ -240,9 +271,13 @@ module xtb_gfnff_param & 0.086496_wp, 0.087053_wp, 0.095395_wp, 0.103738_wp, 0.112081_wp, & & 0.120423_wp, 0.128766_wp, 0.137109_wp, 0.145451_wp, 0.153794_wp, & & 0.323570_wp, 0.233450_wp, 0.268137_wp, 0.307481_wp, 0.316447_wp, & - & 0.400000_wp] + & 0.400000_wp, 0.000000_wp, 0.000000_wp, 0.078710_wp, 0.079266_wp, & + & 0.079822_wp, 0.080379_wp, 0.080935_wp, 0.081491_wp, 0.082047_wp, & + & 0.082603_wp, 0.083159_wp, 0.083716_wp, 0.084272_wp, 0.084828_wp, & + & 0.085384_wp, 0.085940_wp, 0.086496_wp & + &] - real(wp), parameter :: tors_angewChem2020(86) = [& + real(wp), parameter :: tors_angewChem2020(103) = [& & 0.100000_wp, 0.100000_wp, 0.100000_wp, 0.000000_wp, 0.121170_wp, & & 0.260028_wp, 0.222546_wp, 0.250620_wp, 0.256328_wp, 0.400000_wp, & & 0.115000_wp, 0.000000_wp, 0.103731_wp, 0.069103_wp, 0.104280_wp, & @@ -256,13 +291,17 @@ module xtb_gfnff_param & 0.280596_wp, 0.229057_wp, 0.300000_wp, 0.423199_wp, 0.090741_wp, & & 0.076783_wp, 0.310896_wp, 0.309131_wp, 0.307367_wp, 0.305602_wp, & & 0.303838_wp, 0.302073_wp, 0.300309_wp, 0.298544_wp, 0.296779_wp, & - & 0.295015_wp, 0.293250_wp, 0.291486_wp, 0.289721_wp, 0.287957_wp, & + & 0.295015_wp, 0.293250_wp, 0.291486_wp, 0.289721_wp, 0.287957_wp, & ! 70 & 0.286192_wp, 0.284427_wp, 0.257959_wp, 0.231490_wp, 0.205022_wp, & - & 0.178553_wp, 0.152085_wp, 0.125616_wp, 0.099147_wp, 0.072679_wp, & + & 0.178553_wp, 0.152085_wp, 0.125616_wp, 0.099147_wp, 0.072679_wp, & ! 80 & 0.203077_wp, 0.169346_wp, 0.090568_wp, 0.144762_wp, 0.231884_wp, & - & 0.400000_wp] + & 0.400000_wp, 0.090741_wp, 0.076783_wp, 0.310896_wp, 0.309131_wp, & ! 90 + & 0.307367_wp, 0.305602_wp, 0.303838_wp, 0.302073_wp, 0.300309_wp, & + & 0.298544_wp, 0.296779_wp, 0.295015_wp, 0.293250_wp, 0.291486_wp, & ! 100 + & 0.289721_wp, 0.287957_wp, 0.286192_wp & + &] - real(wp), parameter :: tors2_angewChem2020(86) = [& + real(wp), parameter :: tors2_angewChem2020(103) = [& & 1.618678_wp, 1.000000_wp, 0.064677_wp, 0.000000_wp, 0.965814_wp, & & 1.324709_wp, 1.079334_wp, 1.478599_wp, 0.304844_wp, 0.500000_wp, & & 0.029210_wp, 0.000000_wp, 0.417423_wp, 0.334275_wp, 0.817008_wp, & @@ -280,7 +319,11 @@ module xtb_gfnff_param & 0.246091_wp, 0.254931_wp, 0.387526_wp, 0.520121_wp, 0.652716_wp, & & 0.785311_wp, 0.917906_wp, 1.050500_wp, 1.183095_wp, 1.315690_wp, & & 0.219729_wp, 0.344830_wp, 0.331862_wp, 0.767979_wp, 0.536799_wp, & - & 0.500000_wp] + & 0.500000_wp, 0.000000_wp, 0.000000_wp, 0.122336_wp, 0.131176_wp, & + & 0.140015_wp, 0.148855_wp, 0.157695_wp, 0.166534_wp, 0.175374_wp, & + & 0.184214_wp, 0.193053_wp, 0.201893_wp, 0.210733_wp, 0.219572_wp, & + & 0.228412_wp, 0.237252_wp, 0.246091_wp & + &] !---------------------------------------------------------------------------------------- @@ -289,7 +332,8 @@ module xtb_gfnff_param !------------------------------------------------------------------------ !Pauling EN - real(wp), parameter :: en(1:86) = [& + ! Period 7: Pauling EN from https://en.wikipedia.org/wiki/Template:Periodic_table_(electronegativity_by_Pauling_scale) + real(wp), parameter :: en(1:103) = [& & 2.200,3.000,0.980,1.570,2.040,2.550,3.040,3.440,3.980 & & ,4.500,0.930,1.310,1.610,1.900,2.190,2.580,3.160,3.500 & & ,0.820,1.000,1.360,1.540,1.630,1.660,1.550,1.830,1.880 & @@ -298,14 +342,18 @@ module xtb_gfnff_param & ,2.200,1.930,1.690,1.780,1.960,2.050,2.100,2.660,2.600 & &,0.79,0.89,1.10,1.12,1.13,1.14,1.15,1.17,1.18,1.20,1.21,1.22 & &,1.23,1.24,1.25,1.26,1.27,1.3,1.5,1.7,1.9,2.1,2.2,2.2,2.2 & ! value of W-Au modified - &,2.00,1.62,2.33,2.02,2.0,2.2,2.2] + &,2.00,1.62,2.33,2.02,2.0,2.2,2.2 & + &,0.79,0.9,1.1,1.3,1.5,1.38,1.36,1.28,1.13,1.28,1.3,1.3,1.3,1.3,1.3,1.3,1.3 & !Fr-Lr + &] - ! COVALENT RADII, used only in neighbor list determination + ! COVALENT RADII, used only in EHB, XB: abhgfnff_eg1,abhgfnff_eg2new, + ! abhgfnff_eg3, and rbxgfnff_eg ! based on "Atomic Radii of the Elements," M. Mantina, R. Valero, C. J. Cramer, and D. G. Truhlar, ! in CRC Handbook of Chemistry and Physics, 91st Edition (2010-2011), ! edited by W. M. Haynes (CRC Press, Boca Raton, FL, 2010), pages 9-49-9-50; ! corrected Nov. 17, 2010 for the 92nd edition. - real(wp), parameter :: rad(1:86) = [& + ! Period 7: Same book, 95th edition + real(wp), parameter :: rad(1:103) = [& &0.32D0,0.37D0,1.30D0,0.99D0,0.84D0,0.75D0,0.71D0,0.64D0,0.60D0,& &0.62D0,1.60D0,1.40D0,1.24D0,1.14D0,1.09D0,1.04D0,1.00D0,1.01D0,& &2.00D0,1.74D0,1.59D0,1.48D0,1.44D0,1.30D0,1.29D0,1.24D0,1.18D0,& @@ -315,41 +363,52 @@ module xtb_gfnff_param &2.38D0,2.06D0,1.94D0,1.84D0,1.90D0,1.88D0,1.86D0,1.85D0,1.83D0,& &1.82D0,1.81D0,1.80D0,1.79D0,1.77D0,1.77D0,1.78D0,1.74D0,1.64D0,& &1.58D0,1.50D0,1.41D0,1.36D0,1.32D0,1.30D0,1.30D0,1.32D0,1.44D0,& - &1.45D0,1.50D0,1.42D0,1.48D0,1.46D0] + &1.45D0,1.50D0,1.42D0,1.48D0,1.46D0, & + &2.42D0,2.11D0,2.01D0,1.90D0,1.84D0,1.83D0,1.80D0,1.80D0,1.73D0,& !Fr-Am + &1.68D0,1.68D0,1.68D0,1.65D0,1.67D0,1.73D0,1.76D0,1.61D0 & !Cm-Lr + &] - integer, parameter :: metal(86) = (/ & + integer, parameter :: metal(103) = (/ & &0, 0,&!He &1,1, 0, 0, 0, 0, 0, 0,&!Ne &1,1, 1, 0, 0, 0, 0, 0,&!Ar &1,1,2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 0, 0, 0, 0, 0,&!Kr &1,2,2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 0, 0, 0, 0,&!Xe &1,2,2, 2,2,2,2,2,2,2,2,2,2,2,2,2,2, & - & 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 0, 0/)!Rn + & 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 0, 0,&!Rn + &1,2,2, 2,2,2,2,2,2,2,2,2,2,2,2,2,2 & !Fr-Lr + &/) ! At is NOT a metal, Po is borderline but slightly better as metal - integer, private, parameter :: group(86) = (/ & + integer, private, parameter :: group(103) = (/ & &1, 8,&!He &1,2, 3, 4, 5, 6, 7, 8,&!Ne &1,2, 3, 4, 5, 6, 7, 8,&!Ar &1,2,-3, -4,-5,-6,-7,-8,-9,-10,-11,-12,3, 4, 5, 6, 7, 8,&!Kr &1,2,-3, -4,-5,-6,-7,-8,-9,-10,-11,-12,3, 4, 5, 6, 7, 8,&!Xe &1,2,-3, -3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3,-3, & - & -4,-5,-6,-7,-8,-9,-10,-11,-12,3, 4, 5, 6, 7, 8/)!Rn - integer, parameter :: normcn(86) = (/ & ! only for non metals well defined + & -4,-5,-6,-7,-8,-9,-10,-11,-12,3, 4, 5, 6, 7, 8,&!Rn + &1,2,-13, -13,-13,-13,-13,-13,-13,-13,-13,-13,-13,-13,-13,-13,-13 & !Fr-Lr + &/) + integer, parameter :: normcn(103) = (/ & ! only for non metals well defined &1, 0,&!He &4,4, 4, 4, 4, 2, 1, 0,&!Ne &4,4, 4, 4, 4, 2, 1, 0,&!Ar &4,4,4, 4, 6, 6, 6, 6, 6, 6, 4, 4, 4, 4, 4, 4, 1, 0,&!Kr &4,4,4, 4, 6, 6, 6, 6, 6, 6, 4, 4, 4, 4, 4, 4, 1, 0,&!Xe &4,4,4, 4,4,4,4,4,4,4,4,4,4,4,4,4,4, & - & 4, 6, 6, 6, 6, 6, 6, 6, 4, 4, 4, 4, 4, 1, 0/) !Rn - real(wp), parameter :: repz(86) = (/ & + & 4, 6, 6, 6, 6, 6, 6, 6, 4, 4, 4, 4, 4, 1, 0,&!Rn + &4,4,4, 4,4,4,4,4,4,4,4,4,4,4,4,4,4 & !Fr-Lr + &/) + real(wp), parameter :: repz(103) = (/ & &1., 2.,&!He &1.,2., 3.,4.,5.,6.,7.,8.,&!Ne &1.,2., 3.,4.,5.,6.,7.,8.,&!Ar &1.,2.,3., 4.,5.,6.,7.,8.,9.,10.,11.,12.,3.,4.,5.,6.,7.,8.,&!Kr &1.,2.,3., 4.,5.,6.,7.,8.,9.,10.,11.,12.,3.,4.,5.,6.,7.,8.,&!Xe &1.,2.,3., 3.,3.,3.,3.,3.,3.,3.,3.,3.,3.,3.,3.,3.,3., & - & 4.,5.,6.,7.,8.,9.,10.,11.,12.,3.,4.,5.,6.,7.,8./) !Rn + & 4.,5.,6.,7.,8.,9.,10.,11.,12.,3.,4.,5.,6.,7.,8.,& !Rn + &1.,2.,3., 3.,3.,3.,3.,3.,3.,3.,3.,3.,3.,3.,3.,3.,3. & !Fr-Lr + &/) contains @@ -431,7 +490,7 @@ subroutine gfnff_set_param(n, gen, param) ! 3B bond prefactors and D3 stuff k=0 - do i=1,86 + do i=1,103 dum=dble(i) param%zb3atm(i)=-dum*gen%batmscal**(1.d0/3.d0) ! inlcude pre-factor do j=1,i @@ -474,11 +533,11 @@ subroutine gfnff_load_param(version, param, exist) exist = .false. - call init(param, 86) + call init(param, 103) param%en(:) = en param%rad(:) = rad - param%rcov(:) = covalentRadD3(1:86) + param%rcov(:) = covalentRadD3(1:103) param%metal(:) = metal param%group(:) = group param%normcn(:) = normcn @@ -522,17 +581,17 @@ subroutine gfnff_read_param(iunit, param) real(wp) :: xx(20) character(len=256) :: atmp - call init(param, 86) + call init(param, 103) param%en(:) = en param%rad(:) = rad - param%rcov(:) = covalentRadD3(1:86) + param%rcov(:) = covalentRadD3(1:103) param%metal(:) = metal param%group(:) = group param%normcn(:) = normcn param%repz(:) = repz - do i=1,86 + do i=1,103 read(iunit,'(a)')atmp call readl(atmp,xx,nn) param%chi(i)=xx(2) @@ -577,7 +636,7 @@ subroutine gfnff_param_alloc(topo, neigh, n) if (.not.allocated(topo%xbatABl)) allocate( topo%xbatABl(5,topo%natxbAB), source = 0 ) if (.not.allocated(neigh%blist)) allocate( neigh%blist(3,neigh%nbond), source = 0 ) - if (.not.allocated(topo%nr_hb)) allocate( topo%nr_hb(topo%nbond_blist), source = 0 ) + if (.not.allocated(neigh%nr_hb)) allocate( neigh%nr_hb(neigh%nbond), source = 0 ) if (.not.allocated(topo%bond_hb_AH)) allocate( topo%bond_hb_AH(4,topo%bond_hb_nr), source = 0 ) if (.not.allocated(topo%bond_hb_B)) allocate( topo%bond_hb_B(2,topo%b_max,topo%bond_hb_nr), source = 0 ) if (.not.allocated(topo%bond_hb_Bn)) allocate( topo%bond_hb_Bn(topo%bond_hb_nr), source = 0 ) diff --git a/src/gfnff/gfnff_rab.f b/src/gfnff/gfnff_rab.f index 79e676367..4e5b46366 100644 --- a/src/gfnff/gfnff_rab.f +++ b/src/gfnff/gfnff_rab.f @@ -36,9 +36,10 @@ integer function itabrow6(i) end function end interface - real*8 cnfak(86),r0(86),en(86) + real*8 cnfak(103),r0(103),en(103) real*8 ra,rb,k1,k2,den,ff real(wp) :: p(6,2) + real*8 scaleF(103,103) c fitted on PBExa-3c geom set by SG, 9/2018 c H B C N O F SI P S CL GE AS SE BR SN SB TE I together (and for glob par) @@ -62,13 +63,18 @@ integer function itabrow6(i) . 1.89181390, 2.40186898, 1.89282464, 3.09963488, 2.50677823, . 2.61196704, 2.09943450, 2.66930105, 1.78349472, 2.09634533, . 2.00028974, 1.99869908, 2.59072029, 2.54497829, 2.52387890, - . 2.30204667, 1.60119300, 2.00000000, 2.00000000, 2.00000000, - . 2.00000000, 2.00000000, 2.00000000, 2.00000000, 2.00000000, + . 2.30204667, 1.60119300, 2.00000000, 2.00000000, 2.00000000, ! 60 . 2.00000000, 2.00000000, 2.00000000, 2.00000000, 2.00000000, + . 2.00000000, 2.00000000, 2.00000000, 2.00000000, 2.00000000, ! 70 . 2.00000000, 2.30089349, 1.75039077, 1.51785130, 2.62972945, - . 2.75372921, 2.62540906, 2.55860939, 3.32492356, 2.65140898, + . 2.75372921, 2.62540906, 2.55860939, 3.32492356, 2.65140898, ! 80 . 1.52014458, 2.54984804, 1.72021963, 2.69303422, 1.81031095, - . 2.34224386 + . 2.34224386, + . 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 + . 2.00000000 ./ data r0 / . 0.55682207, 0.80966997, 2.49092101, 1.91705642, 1.35974851, @@ -82,13 +88,18 @@ integer function itabrow6(i) . 2.34950167, 2.11644818, 2.47180659, 2.32198800, 2.32809515, . 2.15244869, 2.55958313, 2.59141300, 2.62030465, 2.39935278, ! 50 . 2.56912355, 2.54374096, 2.56914830, 2.53680807, 4.24537037, - . 3.66542289, 3.19903011, 2.80000000, 2.80000000, 2.80000000, ! 60 - . 2.80000000, 2.80000000, 2.80000000, 2.80000000, 2.80000000, - . 2.80000000, 2.80000000, 2.80000000, 2.80000000, 2.80000000, - . 2.80000000, 2.34880037, 2.37597108, 2.49067697, 2.14100577, - . 2.33473532, 2.19498900, 2.12678348, 2.34895048, 2.33422774, + . 3.66542289, 3.39000000, 3.22000000, 3.19000000, 3.17000000, ! 60 + . 3.21000000, 3.05000000, 3.08000000, 3.08000000, 3.06000000, + . 2.92000000, 2.99000000, 3.00000000, 2.91000000, 2.92000000, ! 70 + . 2.76000000, 2.34880037, 2.37597108, 2.49067697, 2.14100577, + . 2.33473532, 2.19498900, 2.12678348, 2.34895048, 2.33422774, ! 80 . 2.86560827, 2.62488837, 2.88376127, 2.75174124, 2.83054552, - . 2.63264944 + . 2.63264944, + . 4.24537037, + . 3.66542289,3.440000000,3.110000000,2.750000000,2.690000000, + . 2.97000000,2.860000000,2.860000000,3.010000000,3.180000000, + . 3.12000000,2.620000000,2.670000000,2.670000000,2.750000000, + . 2.62000000 ./ data cnfak / . 0.17957827, 0.25584045,-0.02485871, 0.00374217, 0.05646607, @@ -100,17 +111,44 @@ integer function itabrow6(i) .-0.00843454, 0.01887497, 0.06876354, 0.01370795,-0.01129196, . 0.07226529, 0.01005367, 0.01541506, 0.05301365, 0.07066571, . 0.07637611, 0.07873977, 0.02997732, 0.04745400, 0.04582912, - . 0.10557321, 0.02167468, 0.05463616, 0.05370913, 0.05985441, + . 0.10557321, 0.02167468, 0.05463616, 0.05370913, 0.05985441, ! 50 . 0.02793994, 0.02922983, 0.02220438, 0.03340460,-0.04110969, - .-0.01987240, 0.07260201, 0.07700000, 0.07700000, 0.07700000, - . 0.07700000, 0.07700000, 0.07700000, 0.07700000, 0.07700000, + .-0.01987240, 0.07260201, 0.07700000, 0.07700000, 0.07700000, ! 60 . 0.07700000, 0.07700000, 0.07700000, 0.07700000, 0.07700000, + . 0.07700000, 0.07700000, 0.07700000, 0.07700000, 0.07700000, ! 70 . 0.07700000, 0.08379100, 0.07314553, 0.05318438, 0.06799334, - . 0.04671159, 0.06758819, 0.09488437, 0.07556405, 0.13384502, + . 0.04671159, 0.06758819, 0.09488437, 0.07556405, 0.13384502, ! 80 . 0.03203572, 0.04235009, 0.03153769,-0.00152488, 0.02714675, - . 0.04800662 + . 0.04800662, ! 86 + .-0.04110969, ! 87 + .-0.01987240, 0.07260201, 0.07700000, 0.07700000, 0.07700000, ! 92 + . 0.07700000, 0.07700000, 0.07700000, 0.07700000, 0.07700000, + . 0.07700000, 0.07700000, 0.07700000, 0.07700000, 0.07700000, ! 102 + . 0.07700000 ! 103 ./ + ! setup factor for Fluorine-Ln or F-Ac bonds + scaleF = 1.0 + scaleF(57:71,9) =0.93533568 + scaleF(89:103,9)=0.93533568 + scaleF(9,57:71) =0.93533568 + scaleF(9,89:103)=0.93533568 + ! for Chlorine + scaleF(57:71,17) =1.0190114 + scaleF(89:103,17)=1.0190114 + scaleF(17,57:71) =1.0190114 + scaleF(17,89:103)=1.0190114 + ! for Bromine + scaleF(57:71,35) =1.0425532 + scaleF(89:103,35)=1.0425532 + scaleF(35,57:71) =1.0425532 + scaleF(35,89:103)=1.0425532 + ! for Iodine + scaleF(57:71,53) =1.0551948 + scaleF(89:103,53)=1.0551948 + scaleF(53,57:71) =1.0551948 + scaleF(53,89:103)=1.0551948 + c global EN polynominal parameters x 10^3 p(1,1)= 29.84522887 p(2,1)= -1.70549806 @@ -140,11 +178,11 @@ integer function itabrow6(i) k1=0.005d0*(p(ir,1)+p(jr,1)) k2=0.005d0*(p(ir,2)+p(jr,2)) ff=1.0d0-k1*den-k2*den**2 - rab(k) =(ra+rb+rab(k))*ff - rabdcn(1,k)=cnfak(ati)*ff - rabdcn(2,k)=cnfak(atj)*ff + rab(k) =(ra+rb+rab(k))*ff*scaleF(ati,atj) + rabdcn(1,k)=cnfak(ati)*ff*scaleF(ati,atj) + rabdcn(2,k)=cnfak(atj)*ff*scaleF(ati,atj) do m=1,n - grab(1:3,m,k)=ff*(cnfak(ati)*dcn(1:3,m,i) + grab(1:3,m,k)=scaleF(ati,atj)*ff*(cnfak(ati)*dcn(1:3,m,i) . +cnfak(atj)*dcn(1:3,m,j)) enddo c-------- @@ -188,8 +226,9 @@ subroutine gfnffrab(n,at,cn,rab) integer m,i,j,k,ii,jj,ati,atj,ir,jr INTEGER iTabRow6,lin - real*8 cnfak(86),r0(86),en(86) + real*8 cnfak(103),r0(103),en(103) real*8 ra,rb,k1,k2,den,ff,p(6,2) + real*8 scaleF(103,103) data en / . 2.30085633, 2.78445145, 1.52956084, 1.51714704, 2.20568300, . 2.49640820, 2.81007174, 4.51078438, 4.67476223, 3.29383610, @@ -208,27 +247,37 @@ subroutine gfnffrab(n,at,cn,rab) . 2.00000000, 2.30089349, 1.75039077, 1.51785130, 2.62972945, . 2.75372921, 2.62540906, 2.55860939, 3.32492356, 2.65140898, . 1.52014458, 2.54984804, 1.72021963, 2.69303422, 1.81031095, - . 2.34224386 + . 2.34224386, + . 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 + . 2.00000000 ./ data r0 / . 0.55682207, 0.80966997, 2.49092101, 1.91705642, 1.35974851, - . 0.98310699, 0.98423007, 0.76716063, 1.06139799, 1.17736822, + . 0.98310699, 0.98423007, 0.76716063, 1.06139799, 1.17736822, ! 10 . 2.85570926, 2.56149012, 2.31673425, 2.03181740, 1.82568535, - . 1.73685958, 1.97498207, 2.00136196, 3.58772537, 2.68096221, + . 1.73685958, 1.97498207, 2.00136196, 3.58772537, 2.68096221, ! 20 . 2.23355957, 2.33135502, 2.15870365, 2.10522128, 2.16376162, - . 2.10804037, 1.96460045, 2.00476257, 2.22628712, 2.43846700, + . 2.10804037, 1.96460045, 2.00476257, 2.22628712, 2.43846700, ! 30 . 2.39408483, 2.24245792, 2.05751204, 2.15427677, 2.27191920, - . 2.19722638, 3.80910350, 3.26020971, 2.99716916, 2.71707818, + . 2.19722638, 3.80910350, 3.26020971, 2.99716916, 2.71707818, ! 40 . 2.34950167, 2.11644818, 2.47180659, 2.32198800, 2.32809515, - . 2.15244869, 2.55958313, 2.59141300, 2.62030465, 2.39935278, + . 2.15244869, 2.55958313, 2.59141300, 2.62030465, 2.39935278, ! 50 . 2.56912355, 2.54374096, 2.56914830, 2.53680807, 4.24537037, - . 3.66542289, 3.19903011, 2.80000000, 2.80000000, 2.80000000, - . 2.80000000, 2.80000000, 2.80000000, 2.80000000, 2.80000000, - . 2.80000000, 2.80000000, 2.80000000, 2.80000000, 2.80000000, - . 2.80000000, 2.34880037, 2.37597108, 2.49067697, 2.14100577, - . 2.33473532, 2.19498900, 2.12678348, 2.34895048, 2.33422774, + . 3.66542289, 3.39000000, 3.22000000, 3.19000000, 3.17000000, ! 60 + . 3.21000000, 3.05000000, 3.08000000, 3.08000000, 3.06000000, + . 2.92000000, 2.99000000, 3.00000000, 2.91000000, 2.92000000, ! 70 + . 2.76000000, 2.34880037, 2.37597108, 2.49067697, 2.14100577, + . 2.33473532, 2.19498900, 2.12678348, 2.34895048, 2.33422774, ! 80 . 2.86560827, 2.62488837, 2.88376127, 2.75174124, 2.83054552, - . 2.63264944 + . 2.63264944, + . 4.24537037, + . 3.66542289,3.440000000,3.110000000,2.750000000,2.690000000, + . 2.97000000,2.860000000,2.860000000,3.010000000,3.180000000, + . 3.12000000,2.620000000,2.670000000,2.670000000,2.750000000, + . 2.62000000 ./ data cnfak / . 0.17957827, 0.25584045,-0.02485871, 0.00374217, 0.05646607, @@ -248,9 +297,36 @@ subroutine gfnffrab(n,at,cn,rab) . 0.07700000, 0.08379100, 0.07314553, 0.05318438, 0.06799334, . 0.04671159, 0.06758819, 0.09488437, 0.07556405, 0.13384502, . 0.03203572, 0.04235009, 0.03153769,-0.00152488, 0.02714675, - . 0.04800662 + . 0.04800662, + .-0.04110969, ! 87 + .-0.01987240, 0.07260201, 0.07700000, 0.07700000, 0.07700000, ! 92 + . 0.07700000, 0.07700000, 0.07700000, 0.07700000, 0.07700000, + . 0.07700000, 0.07700000, 0.07700000, 0.07700000, 0.07700000, ! 102 + . 0.07700000 ! 103 ./ + ! setup factor for Fluorine-Ln or F-Ac bonds + scaleF = 1.0 + scaleF(57:71,9) =0.93533568 + scaleF(89:103,9)=0.93533568 + scaleF(9,57:71) =0.93533568 + scaleF(9,89:103)=0.93533568 + ! for Chlorine + scaleF(57:71,17) =1.0190114 + scaleF(89:103,17)=1.0190114 + scaleF(17,57:71) =1.0190114 + scaleF(17,89:103)=1.0190114 + ! for Bromine + scaleF(57:71,35) =1.0425532 + scaleF(89:103,35)=1.0425532 + scaleF(35,57:71) =1.0425532 + scaleF(35,89:103)=1.0425532 + ! for Iodine + scaleF(57:71,53) =1.0551948 + scaleF(89:103,53)=1.0551948 + scaleF(53,57:71) =1.0551948 + scaleF(53,89:103)=1.0551948 + p(1,1)= 29.84522887 p(2,1)= -1.70549806 p(3,1)= 6.54013762 @@ -277,7 +353,7 @@ subroutine gfnffrab(n,at,cn,rab) k1=0.005d0*(p(ir,1)+p(jr,1)) k2=0.005d0*(p(ir,2)+p(jr,2)) ff=1.0d0-k1*den-k2*den**2 - rab(k)=(ra+rb)*ff + rab(k)=(ra+rb)*ff*scaleF(ati,atj) enddo enddo 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/model_hessian.f90 b/src/model_hessian.f90 index 59f4855ea..d20d5872e 100644 --- a/src/model_hessian.f90 +++ b/src/model_hessian.f90 @@ -2686,7 +2686,7 @@ subroutine gff_ddvopt(Cart,nAtoms,Hess,at,s6,param,topo,neigh) Data rkr,rkf,rkt/0.4500D+00,0.3000D+00,0.75000 / ! adjusted to account for redundant angles in old version !cc VDWx-Parameters (Grimme) used for vdw-correction of model hessian - real*8 c6(54),c66 + real*8 c6(55),c66 ! NOT USED ANYMORE (BJ instead) ! data vander & ! H, He @@ -2736,6 +2736,7 @@ subroutine gff_ddvopt(Cart,nAtoms,Hess,at,s6,param,topo,neigh) if(at(i).gt.54) iANr(i)=at(i) - 18 if(at(i).gt.72) iANr(i)=at(i) - 32 if(at(i).gt.56.and.at(i).lt.72) iANr(i)=39 ! set LNs to Y + if(at(i).gt.86) iANr(i)=55 ! use same c6 for all atom types > 86 enddo profile=.false. @@ -2750,11 +2751,13 @@ subroutine gff_ddvopt(Cart,nAtoms,Hess,at,s6,param,topo,neigh) c6(32)=17.10; c6(33)=16.37; c6(34)=12.64; c6(35)=12.47; c6(36)=12.01 c6(37:48)=24.67; c6(49)=37.32; c6(50)=38.71; c6(51)=38.44 c6(52)=31.74; c6(53)=31.50; c6(54)=29.99 + c6(55)=40.0 !hjw threshold reduced rZero=1.0d-10 n3=3*nAtoms -! Hess = 0.0d0 + Hess = 0.0d0 + ! ! Hessian for tension diff --git a/src/optimizer.f90 b/src/optimizer.f90 index e6fdc9a4d..4d5779276 100644 --- a/src/optimizer.f90 +++ b/src/optimizer.f90 @@ -437,6 +437,7 @@ subroutine ancopt(env,ilog,mol,chk,calc, & call env%error("Could not read hessian from file.", source) return endif + ! do not reset the hessian maxmicro = maxopt ex = .true. diff --git a/src/pqn.f90 b/src/pqn.f90 index cb45ec5f6..bb821fb93 100644 --- a/src/pqn.f90 +++ b/src/pqn.f90 @@ -56,5 +56,7 @@ pure elemental integer function ncore(at) ncore=68 elseif(at.le.86)then ncore=78 + elseif(at.le.103)then + ncore=86 endif end function ncore diff --git a/src/set_module.f90 b/src/set_module.f90 index 61e2a65f3..c81703d56 100644 --- a/src/set_module.f90 +++ b/src/set_module.f90 @@ -781,6 +781,7 @@ subroutine rdcontrol(fname,env,copy_file) case('cube' ); call rdblock(env,set_cube, line,id,copy,err,ncount) case('write' ); call rdblock(env,set_write, line,id,copy,err,ncount) case('gfn' ); call rdblock(env,set_gfn, line,id,copy,err,ncount) + case('ffnb' ); call rdblock(env,set_ffnb, line,id,copy,err,ncount) case('scc' ); call rdblock(env,set_scc, line,id,copy,err,ncount) case('oniom' ); call rdblock(env,set_oniom, line,id,copy,err,ncount) case('opt' ); call rdblock(env,set_opt, line,id,copy,err,ncount) @@ -1499,6 +1500,66 @@ subroutine set_gfn(env,key,val) end select end subroutine set_gfn +subroutine set_ffnb(env,key,val) + implicit none + character(len=*), parameter :: source = 'set_ffnb' + type(TEnvironment), intent(inout) :: env + character(len=*),intent(in) :: key + 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.-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 + read(val(1:l), *) set%ffnb(1,j) + i_start=i+1 + exit + endif + enddo + endif + ! read the neighbors into ffnb now + if (val(i:i).eq.",") then + k = k + 1 + i_end=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 + 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 + !> set ONIOM functionality subroutine set_oniom(env,key,val) diff --git a/src/setparam.f90 b/src/setparam.f90 index 8c1a1ec92..13047c724 100644 --- a/src/setparam.f90 +++ b/src/setparam.f90 @@ -516,6 +516,10 @@ module xtb_setparam !! ------------------------------------------------------------------------ !> PTB settings type(TPTBSetup) :: ptbsetup + !> 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) = -1 end type TSet type(TSet) :: set diff --git a/src/type/anc.f90 b/src/type/anc.f90 index 836d70d80..eba5171c1 100644 --- a/src/type/anc.f90 +++ b/src/type/anc.f90 @@ -501,4 +501,4 @@ subroutine get_normal(self,g_cartesian, g_normal) end subroutine get_normal -end module xtb_type_anc \ No newline at end of file +end module xtb_type_anc diff --git a/src/type/molecule.f90 b/src/type/molecule.f90 index 806fcd2c1..715bf0dc1 100644 --- a/src/type/molecule.f90 +++ b/src/type/molecule.f90 @@ -590,6 +590,8 @@ pure elemental integer function ncore(at) ncore=68 elseif(at.le.86)then ncore=78 + elseif(at.le.103)then + ncore=86 endif end function ncore end subroutine mol_set_nuclear_charge 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_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