jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

New ILP paramters

This commit is contained in:
oywg11 2020-01-16 16:46:04 +02:00
parent f1c79fb914
commit b2b28015c4
8 changed files with 83 additions and 81 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

17
doc/txt/pair_ilp_graphene_hbn.txt
View File

@ -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).

11
doc/txt/pair_kolmogorov_crespi_full.txt
View File

@ -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).

6
potentials/BNCH-old.ILP
View File

@ -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

4
potentials/BNCH.ILP
View File

@ -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 ++++++++++++++++************

16
potentials/BNC_MBD_bulk.ILP Normal file
View File

@ -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

16
potentials/BNC_TS_bulk.ILP Normal file
View File

@ -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

40
src/USER-MISC/pair_ilp_graphene_hbn.cpp
View File

@ -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];

54
src/USER-MISC/pair_kolmogorov_crespi_full.cpp
View File

@ -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];