diff --git a/doc/txt/pair_ilp_graphene_hbn.txt b/doc/txt/pair_ilp_graphene_hbn.txt index 6d922e20f9..d4c2fb2c60 100644 --- a/doc/txt/pair_ilp_graphene_hbn.txt +++ b/doc/txt/pair_ilp_graphene_hbn.txt @@ -72,15 +72,15 @@ a screened Coulomb interaction "coul/shield"_pair_coul_shield.html. Therefore, to simulated the properties of h-BN correctly, this potential must be used in combination with the pair style "coul/shield"_pair_coul_shield.html. -NOTE: Two new sets of parameters of ILP for two-dimensional hexagonal -Materials are presented in "(Ouyang)"_#Ouyang. These parameters provide -a good description in both short- and long-range interaction regimes. +NOTE: Four new sets of parameters of ILP for 2D layered Materials with bilayer and +bulk configurations are presented in "(Ouyang1)"_#Ouyang1 and "(Ouyang2)"_#Ouyang2, respectively. +These parameters provide a good description in both short- and long-range interaction regimes. While the old ILP parameters published in "(Leven2)"_#Leven2 and "(Maaravi)"_#Maaravi2 are only suitable for long-range interaction regime. This feature is essential for simulations in high pressure regime (i.e., the interlayer distance is smaller than the equilibrium -distance). The benchmark tests and comparison of these parameters can -be found in "(Ouyang)"_#Ouyang. +distance). The benchmark tests and comparison of these parameters can +be found in "(Ouyang1)"_#Ouyang1 and "(Ouyang2)"_#Ouyang2. This potential must be used in combination with hybrid/overlay. Other interactions can be set to zero using pair_style {none}. @@ -155,5 +155,8 @@ units, if your simulation does not use {metal} units. :link(Kolmogorov2) [(Kolmogorov)] A. N. Kolmogorov, V. H. Crespi, Phys. Rev. B 71, 235415 (2005). -:link(Ouyang) -[(Ouyang)] W. Ouyang, D. Mandelli, M. Urbakh and O. Hod, Nano Lett. 18, 6009-6016 (2018). +:link(Ouyang1) +[(Ouyang1)] W. Ouyang, D. Mandelli, M. Urbakh and O. Hod, Nano Lett. 18, 6009-6016 (2018). + +:link(Ouyang2) +[(Ouyang2)] W. Ouyang et al., J. Chem. Theory Comput. 16(1), 666-676 (2020). diff --git a/doc/txt/pair_kolmogorov_crespi_full.txt b/doc/txt/pair_kolmogorov_crespi_full.txt index b42027aada..f33b301d28 100644 --- a/doc/txt/pair_kolmogorov_crespi_full.txt +++ b/doc/txt/pair_kolmogorov_crespi_full.txt @@ -61,7 +61,7 @@ list for calculating the normals for each atom pair. NOTE: Two new sets of parameters of KC potential for hydrocarbons, CH.KC (without the taper function) and CH_taper.KC (with the taper function) -are presented in "(Ouyang)"_#Ouyang1. The energy for the KC potential +are presented in "(Ouyang)"_#Ouyang3. The energy for the KC potential with the taper function goes continuously to zero at the cutoff. The parameters in both CH.KC and CH_taper.KC provide a good description in both short- and long-range interaction regimes. While the original @@ -69,7 +69,7 @@ parameters (CC.KC) published in "(Kolmogorov)"_#Kolmogorov1 are only suitable for long-range interaction regime. This feature is essential for simulations in high pressure regime (i.e., the interlayer distance is smaller than the equilibrium distance). The benchmark tests and -comparison of these parameters can be found in "(Ouyang)"_#Ouyang1. +comparison of these parameters can be found in "(Ouyang1)"_#Ouyang3" and (Ouyang2)"_#Ouyang4. This potential must be used in combination with hybrid/overlay. Other interactions can be set to zero using pair_style {none}. @@ -134,5 +134,8 @@ units. :link(Kolmogorov1) [(Kolmogorov)] A. N. Kolmogorov, V. H. Crespi, Phys. Rev. B 71, 235415 (2005) -:link(Ouyang1) -[(Ouyang)] W. Ouyang, D. Mandelli, M. Urbakh and O. Hod, Nano Lett. 18, 6009-6016 (2018). +:link(Ouyang3) +[(Ouyang1)] W. Ouyang, D. Mandelli, M. Urbakh and O. Hod, Nano Lett. 18, 6009-6016 (2018). + +:link(Ouyang4) +[(Ouyang2)] W. Ouyang et al., J. Chem. Theory Comput. 16(1), 666-676 (2020). diff --git a/potentials/BNCH-old.ILP b/potentials/BNCH-old.ILP index e3e834a400..30219e2a1e 100644 --- a/potentials/BNCH-old.ILP +++ b/potentials/BNCH-old.ILP @@ -1,8 +1,8 @@ -# Interlayer Potential for graphitic and boron nitride systems -# +# Interlayer Potential (ILP) for bilayer graphene/graphene, graphene/hBN and hBN/hBN junctions +# The parameters below are fitted against the HSE + MBD DFT referece data from 3.1 A to 15 A. # Cite J. Chem.Theory Comput. 2016, 12, 2896-905 and J. Phys. Chem. C 2017, 121, 22826-22835. -# beta alpha delta epsilon C d sR reff C6 S rcut +# beta alpha delta epsilon C d sR reff C6 S rcut C C 3.22 9.200 1.20 0.010 0.800 15.0 0.704 3.586 522.915 43.363442016573508 2.0 B B 3.10 8.000 1.60 0.460 0.450 15.0 0.800 3.786 1037.322 43.363442016573508 2.0 N N 3.34 8.000 1.20 0.210 0.680 15.0 0.800 3.365 310.433 43.363442016573508 2.0 diff --git a/potentials/BNCH.ILP b/potentials/BNCH.ILP index 6ca56ca186..ed58adf946 100644 --- a/potentials/BNCH.ILP +++ b/potentials/BNCH.ILP @@ -1,5 +1,5 @@ -# Interlayer Potential (ILP) for graphene/graphene, graphene/hBN and hBN/hBN junctions -# +# Interlayer Potential (ILP) for bilayer graphene/graphene, graphene/hBN and hBN/hBN junctions +# The parameters below are fitted against the HSE + MBD DFT referece data from 2.5 A to 15 A. # Cite as W. Ouyang, D. Mandelli, M. Urbakh and O. Hod, Nano Letters 18, 6009-6016 (2018). # # ----------------- Repulsion Potential ------------------++++++++++++++ Vdw Potential ++++++++++++++++************ diff --git a/potentials/BNC_MBD_bulk.ILP b/potentials/BNC_MBD_bulk.ILP new file mode 100644 index 0000000000..9502e0a89a --- /dev/null +++ b/potentials/BNC_MBD_bulk.ILP @@ -0,0 +1,16 @@ +# Interlayer Potential (ILP) for graphite, bulk-hBN and their heterojunctions +# The parameters below are fitted against the HSE + MBD DFT referece data from 2 A to 10 A. +# Cite as W. Ouyang et al., J. Chem. Theory Comput. 16(1), 666-676 (2020). +# +# ------------------------------ Repulsion Potential --------------------++++++++++++++ Vdw Potential ++++++++++++++++************ +# MBD-HSE beta(A) alpha delta(A) epsilon(meV) C(meV) d sR reff(A) C6(meV*A^6) S rcut +C C 3.1894274136 8.2113165501 1.2600313066 0.0106237125 38.9820878926 10.9736146687 0.7869029010 3.4578620004 25249.6185284695 1.0000 2.0 +B B 3.2146562020 7.1651845022 1.7458546494 11.0735774589 15.4819142891 15.4815063183 0.8550308760 3.4423900567 49498.4383474008 1.0000 2.0 +N N 3.3006476373 6.9225730132 1.4844871606 7.9907993147 46.6114968784 16.9081462104 0.7584806746 3.3265576243 14810.6448568309 1.0000 2.0 +B N 3.1708595336 8.5168240743 2.8657479230 5.4561495348 2.5548134497 13.5321053144 0.8863432069 3.4553049811 24670.8164462408 1.0000 2.0 +C B 3.1007371653 5.1145801996 3.8386588076 18.2345048230 1.1901887968 10.2155326647 0.7686152602 3.5030009241 39262.8518949659 1.0000 2.0 +C N 3.3172548125 10.3496923621 1.3792655319 16.3162761182 19.5690538017 15.7748377566 0.5645056777 3.2659337344 19963.0795570299 1.0000 2.0 +# Symmetric part +B C 3.1007371653 5.1145801996 3.8386588076 18.2345048230 1.1901887968 10.2155326647 0.7686152602 3.5030009241 39262.8518949659 1.0000 2.0 +N C 3.3172548125 10.3496923621 1.3792655319 16.3162761182 19.5690538017 15.7748377566 0.5645056777 3.2659337344 19963.0795570299 1.0000 2.0 +N B 3.1708595336 8.5168240743 2.8657479230 5.4561495348 2.5548134497 13.5321053144 0.8863432069 3.4553049811 24670.8164462408 1.0000 2.0 diff --git a/potentials/BNC_TS_bulk.ILP b/potentials/BNC_TS_bulk.ILP new file mode 100644 index 0000000000..1baca2ace0 --- /dev/null +++ b/potentials/BNC_TS_bulk.ILP @@ -0,0 +1,16 @@ +# Interlayer Potential (ILP) for graphite, bulk-hBN and their heterojunctions +# The parameters below are fitted against the HSE + TS DFT referece data from 2 A to 10 A. +# Cite as W. Ouyang et al., J. Chem. Theory Comput. 16(1), 666-676 (2020). +# +# ------------------------------ Repulsion Potential ------------------++++++++++++++ Vdw Potential ++++++++++++++++************ +# TS-HSE beta(A) alpha delta(A) epsilon(meV) C(meV) d sR reff(A) C6(meV*A^6) S rcut +C C 3.1911991861 8.8422960372 1.1312263335 0.0863236828 33.4354373112 10.0195636456 0.9251350921 3.4842048950 32402.5447674022 1.0000 2.0 +B B 3.5386170858 5.1268088040 2.2006291426 12.8752511690 27.5894275824 13.3599605541 0.8414408912 3.6431051884 99513.2942026427 1.0000 2.0 +N N 3.5915052274 3.2218485106 1.4354352315 6.6765695916 73.1025964821 13.0709665964 0.7465905646 3.3082847788 74823.5568235260 1.0000 2.0 +B N 3.9928670174 7.8553264966 2.5853334572 4.5784769626 2.3283590129 16.2664653527 0.8669270315 3.9824166141 84699.9699356515 1.0000 2.0 +C B 3.0183281153 9.8126192181 3.6974442648 22.1591263112 0.8264956969 11.1782507988 0.9510337752 3.8465295183 40165.3497011995 1.0000 2.0 +C N 3.4895869665 10.1614317973 1.1614801250 4.2614500372 11.1810783861 11.0390855050 0.9257256781 3.2512446017 29066.8955087607 1.0000 2.0 +# Symmetric part +B C 3.0183281153 9.8126192181 3.6974442648 22.1591263112 0.8264956969 11.1782507988 0.9510337752 3.8465295183 40165.3497011995 1.0000 2.0 +N C 3.4895869665 10.1614317973 1.1614801250 4.2614500372 11.1810783861 11.0390855050 0.9257256781 3.2512446017 29066.8955087607 1.0000 2.0 +N B 3.9928670174 7.8553264966 2.5853334572 4.5784769626 2.3283590129 16.2664653527 0.8669270315 3.9824166141 84699.9699356515 1.0000 2.0 diff --git a/src/USER-MISC/pair_ilp_graphene_hbn.cpp b/src/USER-MISC/pair_ilp_graphene_hbn.cpp index e998abf005..bdfbbaa5b7 100644 --- a/src/USER-MISC/pair_ilp_graphene_hbn.cpp +++ b/src/USER-MISC/pair_ilp_graphene_hbn.cpp @@ -442,7 +442,7 @@ void PairILPGrapheneHBN::compute(int eflag, int vflag) if (vflag_fdotr) virial_fdotr_compute(); } -/* ---------------------------------------------------------------------- +/* ---------------------------------------------------------------------- van der Waals forces and energy ------------------------------------------------------------------------- */ @@ -523,7 +523,7 @@ void PairILPGrapheneHBN::calc_FvdW(int eflag, int /* vflag */) // derivatives fpair = -6.0*p.C6*r8inv/TSvdw + p.C6*p.d/p.seff*(TSvdw-1.0)*TSvdw2inv*r8inv*r; - fsum = fpair*Tap - Vilp*dTap/r; + fsum = fpair*Tap - Vilp*dTap/r; f[i][0] += fsum*delx; f[i][1] += fsum*dely; @@ -540,7 +540,7 @@ void PairILPGrapheneHBN::calc_FvdW(int eflag, int /* vflag */) } } -/* ---------------------------------------------------------------------- +/* ---------------------------------------------------------------------- Repulsive forces and energy ------------------------------------------------------------------------- */ @@ -576,9 +576,6 @@ void PairILPGrapheneHBN::calc_FRep(int eflag, int /* vflag */) // loop over neighbors of owned atoms for (ii = 0; ii < inum; ii++) { i = ilist[ii]; - if (ILP_numneigh[i] == -1) { - continue; - } xtmp = x[i][0]; ytmp = x[i][1]; ztmp = x[i][2]; @@ -589,9 +586,6 @@ void PairILPGrapheneHBN::calc_FRep(int eflag, int /* vflag */) for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; j &= NEIGHMASK; - if (ILP_numneigh[j] == -1) { - continue; - } jtype = type[j]; delx = xtmp - x[j][0]; @@ -634,12 +628,12 @@ void PairILPGrapheneHBN::calc_FRep(int eflag, int /* vflag */) dprodnorm1[0] = dnormdri[0][0][i]*delx + dnormdri[1][0][i]*dely + dnormdri[2][0][i]*delz; dprodnorm1[1] = dnormdri[0][1][i]*delx + dnormdri[1][1][i]*dely + dnormdri[2][1][i]*delz; dprodnorm1[2] = dnormdri[0][2][i]*delx + dnormdri[1][2][i]*dely + dnormdri[2][2][i]*delz; - fp1[0] = prodnorm1*normal[i][0]*fpair1; - fp1[1] = prodnorm1*normal[i][1]*fpair1; - fp1[2] = prodnorm1*normal[i][2]*fpair1; - fprod1[0] = prodnorm1*dprodnorm1[0]*fpair1; - fprod1[1] = prodnorm1*dprodnorm1[1]*fpair1; - fprod1[2] = prodnorm1*dprodnorm1[2]*fpair1; + fp1[0] = prodnorm1*normal[i][0]*fpair1; + fp1[1] = prodnorm1*normal[i][1]*fpair1; + fp1[2] = prodnorm1*normal[i][2]*fpair1; + fprod1[0] = prodnorm1*dprodnorm1[0]*fpair1; + fprod1[1] = prodnorm1*dprodnorm1[1]*fpair1; + fprod1[2] = prodnorm1*dprodnorm1[2]*fpair1; fkcx = (delx*fsum - fp1[0])*Tap - Vilp*dTap*delx/r; fkcy = (dely*fsum - fp1[1])*Tap - Vilp*dTap*dely/r; @@ -741,17 +735,8 @@ void PairILPGrapheneHBN::ILP_neigh() } // loop over jj ILP_firstneigh[i] = neighptr; - if (n == 3) { - ILP_numneigh[i] = n; - } - else if (n < 3) { - if (i < inum) { - ILP_numneigh[i] = n; - } else { - ILP_numneigh[i] = -1; - } - } - else if (n > 3) error->one(FLERR,"There are too many neighbors for some atoms, please check your configuration"); + ILP_numneigh[i] = n; + if (n > 3) error->one(FLERR,"There are too many neighbors for some atoms, please check your configuration"); ipage->vgot(n); if (ipage->status()) @@ -814,9 +799,6 @@ void PairILPGrapheneHBN::calc_normal() } } - if (ILP_numneigh[i] == -1) { - continue; - } xtp = x[i][0]; ytp = x[i][1]; ztp = x[i][2]; diff --git a/src/USER-MISC/pair_kolmogorov_crespi_full.cpp b/src/USER-MISC/pair_kolmogorov_crespi_full.cpp index eea0b1261c..da6e892d54 100644 --- a/src/USER-MISC/pair_kolmogorov_crespi_full.cpp +++ b/src/USER-MISC/pair_kolmogorov_crespi_full.cpp @@ -444,7 +444,7 @@ void PairKolmogorovCrespiFull::compute(int eflag, int vflag) if (vflag_fdotr) virial_fdotr_compute(); } -/* ---------------------------------------------------------------------- +/* ---------------------------------------------------------------------- van der Waals forces and energy ------------------------------------------------------------------------- */ @@ -519,11 +519,11 @@ void PairKolmogorovCrespiFull::calc_FvdW(int eflag, int /* vflag */) dTap = calc_dTap(r,Rcut); } else {Tap = 1.0; dTap = 0.0;} - Vkc = -p.A*p.z06*r6inv; + Vkc = -p.A*p.z06*r6inv; // derivatives fpair = -6.0*p.A*p.z06*r8inv; - fsum = fpair*Tap - Vkc*dTap/r; + fsum = fpair*Tap - Vkc*dTap/r; f[i][0] += fsum*delx; f[i][1] += fsum*dely; @@ -540,7 +540,7 @@ void PairKolmogorovCrespiFull::calc_FvdW(int eflag, int /* vflag */) } } -/* ---------------------------------------------------------------------- +/* ---------------------------------------------------------------------- Repulsive forces and energy ------------------------------------------------------------------------- */ @@ -576,9 +576,6 @@ void PairKolmogorovCrespiFull::calc_FRep(int eflag, int /* vflag */) // loop over neighbors of owned atoms for (ii = 0; ii < inum; ii++) { i = ilist[ii]; - if (KC_numneigh[i] == -1) { - continue; - } xtmp = x[i][0]; ytmp = x[i][1]; ztmp = x[i][2]; @@ -589,9 +586,6 @@ void PairKolmogorovCrespiFull::calc_FRep(int eflag, int /* vflag */) for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; j &= NEIGHMASK; - if (KC_numneigh[j] == -1) { - continue; - } jtype = type[j]; delx = xtmp - x[j][0]; @@ -607,11 +601,11 @@ void PairKolmogorovCrespiFull::calc_FRep(int eflag, int /* vflag */) r = sqrt(rsq); - // turn on/off taper function - if (tap_flag) { - Tap = calc_Tap(r,sqrt(cutsq[itype][jtype])); - dTap = calc_dTap(r,sqrt(cutsq[itype][jtype])); - } else {Tap = 1.0; dTap = 0.0;} + // turn on/off taper function + if (tap_flag) { + Tap = calc_Tap(r,sqrt(cutsq[itype][jtype])); + dTap = calc_dTap(r,sqrt(cutsq[itype][jtype])); + } else {Tap = 1.0; dTap = 0.0;} // Calculate the transverse distance prodnorm1 = normal[i][0]*delx + normal[i][1]*dely + normal[i][2]*delz; @@ -626,7 +620,7 @@ void PairKolmogorovCrespiFull::calc_FRep(int eflag, int /* vflag */) sumC11 = (p.C2 + 2.0*p.C4*rho_ij)*p.delta2inv; frho_ij = exp1*sumC1; sumCff = 0.5*p.C + frho_ij; - Vkc = exp0*sumCff; + Vkc = exp0*sumCff; // derivatives fpair = p.lambda*exp0/r*sumCff; @@ -653,10 +647,10 @@ void PairKolmogorovCrespiFull::calc_FRep(int eflag, int /* vflag */) f[j][1] -= fkcy; f[j][2] -= fkcz; - // calculate the forces acted on the neighbors of atom i from atom j - KC_neighs_i = KC_firstneigh[i]; - for (kk = 0; kk < KC_numneigh[i]; kk++) { - k = KC_neighs_i[kk]; + // calculate the forces acted on the neighbors of atom i from atom j + KC_neighs_i = KC_firstneigh[i]; + for (kk = 0; kk < KC_numneigh[i]; kk++) { + k = KC_neighs_i[kk]; if (k == i) continue; // derivatives of the product of rij and ni respect to rk, k=0,1,2, where atom k is the neighbors of atom i dprodnorm1[0] = dnormal[0][0][kk][i]*delx + dnormal[1][0][kk][i]*dely + dnormal[2][0][kk][i]*delz; @@ -672,7 +666,7 @@ void PairKolmogorovCrespiFull::calc_FRep(int eflag, int /* vflag */) delkj[1] = x[k][1] - x[j][1]; delkj[2] = x[k][2] - x[j][2]; if (evflag) ev_tally_xyz(k,j,nlocal,newton_pair,0.0,0.0,fk[0],fk[1],fk[2],delkj[0],delkj[1],delkj[2]); - } + } if (eflag) { if (tap_flag) pvector[1] += evdwl = Tap*Vkc; @@ -746,17 +740,8 @@ void PairKolmogorovCrespiFull::KC_neigh() } KC_firstneigh[i] = neighptr; - if (n == 3) { - KC_numneigh[i] = n; - } - else if (n < 3) { - if (i < inum) { - KC_numneigh[i] = n; - } else { - KC_numneigh[i] = -1; - } - } - else if (n > 3) error->one(FLERR,"There are too many neighbors for some atoms, please check your configuration"); + KC_numneigh[i] = n; + if (n > 3) error->one(FLERR,"There are too many neighbors for some atoms, please check your configuration"); ipage->vgot(n); if (ipage->status()) @@ -790,7 +775,7 @@ void PairKolmogorovCrespiFull::calc_normal() memory->create(dnormal,3,3,3,nmax,"KolmogorovCrespiFull:dnormal"); } - inum = list->inum; + inum = list->inum; ilist = list->ilist; //Calculate normals for (ii = 0; ii < inum; ii++) { @@ -819,9 +804,6 @@ void PairKolmogorovCrespiFull::calc_normal() } } - if (KC_numneigh[i] == -1) { - continue; - } xtp = x[i][0]; ytp = x[i][1]; ztp = x[i][2];