diff --git a/doc/src/pair_ilp_tmd.rst b/doc/src/pair_ilp_tmd.rst index 3ac22e5140..b2fd405a39 100644 --- a/doc/src/pair_ilp_tmd.rst +++ b/doc/src/pair_ilp_tmd.rst @@ -30,7 +30,7 @@ Description """"""""""" The *ilp/tmd* style computes the registry-dependent interlayer -potential (ILP) potential for itransition metal dichalcogenide (TMD) +potential (ILP) potential for itransition metal dichalcogenide (TMD) as described in :ref:`(Ouyang3) `. .. math:: @@ -57,11 +57,11 @@ calculating the normals. .. note:: - Since each MX2 (M = Mo, W and X = S, Se Te) layer contains two + Since each MX2 (M = Mo, W and X = S, Se Te) layer contains two sublayers of X atoms and one sublayer of M atoms, the definition of the normal vectors used for graphene and h-BN is no longer valid for TMDs. - In :ref:`(Ouyang3) `, a new definition is proposed, where for - each atom `i`, its six nearest neighboring atoms belonging to the same + In :ref:`(Ouyang3) `, a new definition is proposed, where for + each atom `i`, its six nearest neighboring atoms belonging to the same sublayer are chosen to define the normal vector `{\bf n}_i`. The parameter file (e.g. TMD.ILP), is intended for use with *metal* @@ -75,7 +75,7 @@ list for calculating the normals for each atom pair. The parameters presented in the parameter file (e.g. BNCH.ILP), are fitted with taper function by setting the cutoff equal to 16.0 Angstrom. Using different cutoff or taper function should be careful. - These parameters provide a good description in both short- and long-range + These parameters provide a good description in both short- and long-range interaction regimes. 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 diff --git a/src/INTERLAYER/pair_ilp_graphene_hbn.cpp b/src/INTERLAYER/pair_ilp_graphene_hbn.cpp index 8464c44ddf..a13a3464ca 100644 --- a/src/INTERLAYER/pair_ilp_graphene_hbn.cpp +++ b/src/INTERLAYER/pair_ilp_graphene_hbn.cpp @@ -58,8 +58,7 @@ static const char cite_ilp[] = /* ---------------------------------------------------------------------- */ -PairILPGrapheneHBN::PairILPGrapheneHBN(LAMMPS *lmp) - : Pair(lmp), variant(ILP_GrhBN) +PairILPGrapheneHBN::PairILPGrapheneHBN(LAMMPS *lmp) : Pair(lmp), variant(ILP_GrhBN) { restartinfo = 0; one_coeff = 1; @@ -212,20 +211,20 @@ void PairILPGrapheneHBN::read_file(char *filename) if (comm->me == 0) { std::string potential_name; switch (variant) { - case ILP_GrhBN: - potential_name = "ilp/graphene/hbn"; - break; + case ILP_GrhBN: + potential_name = "ilp/graphene/hbn"; + break; - case ILP_TMD: - potential_name = "ilp/tmd"; - break; + case ILP_TMD: + potential_name = "ilp/tmd"; + break; - case SAIP_METAL: - potential_name = "saip/metal"; - break; + case SAIP_METAL: + potential_name = "saip/metal"; + break; - default: - error->one(FLERR,"Unknown ILP style variant {}",variant); + default: + error->one(FLERR, "Unknown ILP style variant {}", variant); } PotentialFileReader reader(lmp, filename, potential_name, unit_convert_flag); char *line; @@ -342,10 +341,8 @@ void PairILPGrapheneHBN::read_file(char *filename) void PairILPGrapheneHBN::init_style() { - if (force->newton_pair == 0) - error->all(FLERR, "Pair style ilp/* requires newton pair on"); - if (!atom->molecule_flag) - error->all(FLERR, "Pair style ilp/* requires atom attribute molecule"); + if (force->newton_pair == 0) error->all(FLERR, "Pair style ilp/* requires newton pair on"); + if (!atom->molecule_flag) error->all(FLERR, "Pair style ilp/* requires atom attribute molecule"); // need a full neighbor list, including neighbors of ghosts diff --git a/src/INTERLAYER/pair_ilp_graphene_hbn.h b/src/INTERLAYER/pair_ilp_graphene_hbn.h index deae47fbee..171c30f1bf 100644 --- a/src/INTERLAYER/pair_ilp_graphene_hbn.h +++ b/src/INTERLAYER/pair_ilp_graphene_hbn.h @@ -68,8 +68,8 @@ class PairILPGrapheneHBN : public Pair { double ***dnormdri; double ****dnormal; -// adds for ilp/tmd - int Nnei; // max # of nearest neighbors for one atom + // adds for ilp/tmd + int Nnei; // max # of nearest neighbors for one atom double **dnn; double **vect; double **pvet; diff --git a/src/INTERLAYER/pair_ilp_tmd.cpp b/src/INTERLAYER/pair_ilp_tmd.cpp index df78c8a31e..943d90e6a5 100644 --- a/src/INTERLAYER/pair_ilp_tmd.cpp +++ b/src/INTERLAYER/pair_ilp_tmd.cpp @@ -42,16 +42,17 @@ using namespace InterLayer; #define DELTA 4 #define PGDELTA 1 -static const char cite_ilp_tmd[] = - "ilp/tmd potential doi/10.1021/acs.jctc.1c00782\n" - "@Article{Ouyang2021\n" - " author = {W. Ouyang, R. Sofer, X. Gao, J. Hermann, A. Tkatchenko, L. Kronik, M. Urbakh, and O. Hod},\n" - " title = {Anisotropic Interlayer Force Field for Transition Metal Dichalcogenides: The Case of Molybdenum Disulfide},\n" - " journal = {J. Chem. Theory Comput.},\n" - " volume = 17,\n" - " pages = {7237–7245}\n" - " year = 2021,\n" - "}\n\n"; +static const char cite_ilp_tmd[] = "ilp/tmd potential doi/10.1021/acs.jctc.1c00782\n" + "@Article{Ouyang2021\n" + " author = {W. Ouyang, R. Sofer, X. Gao, J. Hermann, A. " + "Tkatchenko, L. Kronik, M. Urbakh, and O. Hod},\n" + " title = {Anisotropic Interlayer Force Field for Transition " + "Metal Dichalcogenides: The Case of Molybdenum Disulfide},\n" + " journal = {J. Chem. Theory Comput.},\n" + " volume = 17,\n" + " pages = {7237–7245}\n" + " year = 2021,\n" + "}\n\n"; /* ---------------------------------------------------------------------- */ @@ -79,21 +80,20 @@ void PairILPTMD::settings(int narg, char **arg) if (narg == 2) tap_flag = utils::numeric(FLERR, arg[1], false, lmp); } -/* ---------------------------------------------------------------------- +/* ---------------------------------------------------------------------- Repulsive forces and energy ------------------------------------------------------------------------- */ void PairILPTMD::calc_FRep(int eflag, int /* vflag */) { - int i,j,ii,jj,inum,jnum,itype,jtype,k,kk; - double prodnorm1,fkcx,fkcy,fkcz; - double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair,fpair1; - double rsq,r,Rcut,rhosq1,exp0,exp1,r2inv,r6inv,r8inv,Tap,dTap,Vilp; - double frho1,TSvdw,TSvdw2inv,Erep,fsum,rdsq1; - int *ilist,*jlist,*numneigh,**firstneigh; + int i, j, ii, jj, inum, jnum, itype, jtype, k, kk; + double prodnorm1, fkcx, fkcy, fkcz; + double xtmp, ytmp, ztmp, delx, dely, delz, evdwl, fpair, fpair1; + double rsq, r, Rcut, rhosq1, exp0, exp1, r2inv, r6inv, r8inv, Tap, dTap, Vilp; + double frho1, TSvdw, TSvdw2inv, Erep, fsum, rdsq1; + int *ilist, *jlist, *numneigh, **firstneigh; int *ILP_neighs_i; - evdwl = 0.0; double **x = atom->x; @@ -131,58 +131,64 @@ void PairILPTMD::calc_FRep(int eflag, int /* vflag */) delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; delz = ztmp - x[j][2]; - rsq = delx*delx + dely*dely + delz*delz; + rsq = delx * delx + dely * dely + delz * delz; // only include the interation between different layers if (rsq < cutsq[itype][jtype] && atom->molecule[i] != atom->molecule[j]) { int iparam_ij = elem2param[map[itype]][map[jtype]]; - Param& p = params[iparam_ij]; + Param &p = params[iparam_ij]; r = sqrt(rsq); // turn on/off taper function if (tap_flag) { Rcut = sqrt(cutsq[itype][jtype]); - Tap = calc_Tap(r,Rcut); - dTap = calc_dTap(r,Rcut); - } else {Tap = 1.0; dTap = 0.0;} + Tap = calc_Tap(r, Rcut); + dTap = calc_dTap(r, Rcut); + } else { + Tap = 1.0; + dTap = 0.0; + } // Calculate the transverse distance // note that rho_ij does not equal to rho_ji except when normals are all along z - prodnorm1 = normal[i][0]*delx + normal[i][1]*dely + normal[i][2]*delz; - rhosq1 = rsq - prodnorm1*prodnorm1; // rho_ij - rdsq1 = rhosq1*p.delta2inv; // (rho_ij/delta)^2 + prodnorm1 = normal[i][0] * delx + normal[i][1] * dely + normal[i][2] * delz; + rhosq1 = rsq - prodnorm1 * prodnorm1; // rho_ij + rdsq1 = rhosq1 * p.delta2inv; // (rho_ij/delta)^2 // store exponents - exp0 = exp(-p.lambda*(r-p.z0)); + exp0 = exp(-p.lambda * (r - p.z0)); exp1 = exp(-rdsq1); - frho1 = exp1*p.C; - Erep = 0.5*p.epsilon + frho1; - Vilp = exp0*Erep; + frho1 = exp1 * p.C; + Erep = 0.5 * p.epsilon + frho1; + Vilp = exp0 * Erep; // derivatives - fpair = p.lambda*exp0/r*Erep; - fpair1 = 2.0*exp0*frho1*p.delta2inv; + fpair = p.lambda * exp0 / r * Erep; + fpair1 = 2.0 * exp0 * frho1 * p.delta2inv; fsum = fpair + fpair1; // derivatives of the product of rij and ni, the result is a vector - dprodnorm1[0] = dnormdri[i][0][0]*delx + dnormdri[i][1][0]*dely + dnormdri[i][2][0]*delz; - dprodnorm1[1] = dnormdri[i][0][1]*delx + dnormdri[i][1][1]*dely + dnormdri[i][2][1]*delz; - dprodnorm1[2] = dnormdri[i][0][2]*delx + dnormdri[i][1][2]*dely + dnormdri[i][2][2]*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; + dprodnorm1[0] = + dnormdri[i][0][0] * delx + dnormdri[i][1][0] * dely + dnormdri[i][2][0] * delz; + dprodnorm1[1] = + dnormdri[i][0][1] * delx + dnormdri[i][1][1] * dely + dnormdri[i][2][1] * delz; + dprodnorm1[2] = + dnormdri[i][0][2] * delx + dnormdri[i][1][2] * dely + dnormdri[i][2][2] * 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; - fkcx = (delx*fsum - fp1[0])*Tap - Vilp*dTap*delx/r; - fkcy = (dely*fsum - fp1[1])*Tap - Vilp*dTap*dely/r; - fkcz = (delz*fsum - fp1[2])*Tap - Vilp*dTap*delz/r; + fkcx = (delx * fsum - fp1[0]) * Tap - Vilp * dTap * delx / r; + fkcy = (dely * fsum - fp1[1]) * Tap - Vilp * dTap * dely / r; + fkcz = (delz * fsum - fp1[2]) * Tap - Vilp * dTap * delz / r; - f[i][0] += fkcx - fprod1[0]*Tap; - f[i][1] += fkcy - fprod1[1]*Tap; - f[i][2] += fkcz - fprod1[2]*Tap; + f[i][0] += fkcx - fprod1[0] * Tap; + f[i][1] += fkcy - fprod1[1] * Tap; + f[i][2] += fkcz - fprod1[2] * Tap; f[j][0] -= fkcx; f[j][1] -= fkcy; f[j][2] -= fkcz; @@ -193,23 +199,29 @@ void PairILPTMD::calc_FRep(int eflag, int /* vflag */) k = ILP_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[i][0][kk][0]*delx + dnormal[i][1][kk][0]*dely + dnormal[i][2][kk][0]*delz; - dprodnorm1[1] = dnormal[i][0][kk][1]*delx + dnormal[i][1][kk][1]*dely + dnormal[i][2][kk][1]*delz; - dprodnorm1[2] = dnormal[i][0][kk][2]*delx + dnormal[i][1][kk][2]*dely + dnormal[i][2][kk][2]*delz; - fk[0] = (-prodnorm1*dprodnorm1[0]*fpair1)*Tap; - fk[1] = (-prodnorm1*dprodnorm1[1]*fpair1)*Tap; - fk[2] = (-prodnorm1*dprodnorm1[2]*fpair1)*Tap; + dprodnorm1[0] = dnormal[i][0][kk][0] * delx + dnormal[i][1][kk][0] * dely + + dnormal[i][2][kk][0] * delz; + dprodnorm1[1] = dnormal[i][0][kk][1] * delx + dnormal[i][1][kk][1] * dely + + dnormal[i][2][kk][1] * delz; + dprodnorm1[2] = dnormal[i][0][kk][2] * delx + dnormal[i][1][kk][2] * dely + + dnormal[i][2][kk][2] * delz; + fk[0] = (-prodnorm1 * dprodnorm1[0] * fpair1) * Tap; + fk[1] = (-prodnorm1 * dprodnorm1[1] * fpair1) * Tap; + fk[2] = (-prodnorm1 * dprodnorm1[2] * fpair1) * Tap; f[k][0] += fk[0]; f[k][1] += fk[1]; f[k][2] += fk[2]; delki[0] = x[k][0] - x[i][0]; delki[1] = x[k][1] - x[i][1]; delki[2] = x[k][2] - x[i][2]; - if (evflag) ev_tally_xyz(k,j,nlocal,newton_pair,0.0,0.0,fk[0],fk[1],fk[2],delki[0],delki[1],delki[2]); + if (evflag) + ev_tally_xyz(k, j, nlocal, newton_pair, 0.0, 0.0, fk[0], fk[1], fk[2], delki[0], + delki[1], delki[2]); } - if (eflag) pvector[1] += evdwl = Tap*Vilp; - if (evflag) ev_tally_xyz(i,j,nlocal,newton_pair,evdwl,0.0,fkcx,fkcy,fkcz,delx,dely,delz); + if (eflag) pvector[1] += evdwl = Tap * Vilp; + if (evflag) + ev_tally_xyz(i, j, nlocal, newton_pair, evdwl, 0.0, fkcx, fkcy, fkcz, delx, dely, delz); } } } @@ -221,11 +233,11 @@ void PairILPTMD::calc_FRep(int eflag, int /* vflag */) void PairILPTMD::ILP_neigh() { - int i,j,l,ii,jj,ll,n,allnum,jnum,itype,jtype,ltype,imol,jmol,count; - double xtmp,ytmp,ztmp,delx,dely,delz,deljx,deljy,deljz,rsq,rsqlj; - int *ilist,*jlist,*numneigh,**firstneigh; + int i, j, l, ii, jj, ll, n, allnum, jnum, itype, jtype, ltype, imol, jmol, count; + double xtmp, ytmp, ztmp, delx, dely, delz, deljx, deljy, deljz, rsq, rsqlj; + int *ilist, *jlist, *numneigh, **firstneigh; int *neighsort; - int neighptr[10],check[10]; + int neighptr[10], check[10]; double **x = atom->x; int *type = atom->type; @@ -234,8 +246,8 @@ void PairILPTMD::ILP_neigh() maxlocal = atom->nmax; memory->destroy(ILP_numneigh); memory->sfree(ILP_firstneigh); - memory->create(ILP_numneigh,maxlocal,"ILPTMD:numneigh"); - ILP_firstneigh = (int **) memory->smalloc(maxlocal*sizeof(int *),"ILPTMD:firstneigh"); + memory->create(ILP_numneigh, maxlocal, "ILPTMD:numneigh"); + ILP_firstneigh = (int **) memory->smalloc(maxlocal * sizeof(int *), "ILPTMD:firstneigh"); } allnum = list->inum + list->gnum; @@ -254,7 +266,7 @@ void PairILPTMD::ILP_neigh() //initialize varibles n = 0; neighsort = ipage->vget(); - for (ll = 0; ll < 10; ll++){ + for (ll = 0; ll < 10; ll++) { neighptr[ll] = -1; check[ll] = -1; } @@ -275,67 +287,61 @@ void PairILPTMD::ILP_neigh() delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; delz = ztmp - x[j][2]; - rsq = delx*delx + dely*dely + delz*delz; - + rsq = delx * delx + dely * dely + delz * delz; + // check if the atom i is TMD, i.e., Mo/S/W/Se - if(strcmp(elements[itype],"Mo") == 0 || strcmp(elements[itype],"W") == 0 - || strcmp(elements[itype],"S") == 0 || strcmp(elements[itype],"Se") == 0 - || strcmp(elements[itype],"Te") == 0) { + if (strcmp(elements[itype], "Mo") == 0 || strcmp(elements[itype], "W") == 0 || + strcmp(elements[itype], "S") == 0 || strcmp(elements[itype], "Se") == 0 || + strcmp(elements[itype], "Te") == 0) { if (rsq != 0 && rsq < cutILPsq[itype][jtype] && imol == jmol && type[i] == type[j]) { neighptr[n++] = j; } + } else { // atom i is C, B, N or H. + if (rsq != 0 && rsq < cutILPsq[itype][jtype] && imol == jmol) { neighptr[n++] = j; } } - else {// atom i is C, B, N or H. - if (rsq != 0 && rsq < cutILPsq[itype][jtype] && imol == jmol) { - neighptr[n++] = j; - } - } - } // loop over jj + } // loop over jj // if atom i is Mo/W/S/Se/Te, then sorting the orders of neighbors - if(strcmp(elements[itype],"Mo") == 0 || strcmp(elements[itype],"W") == 0 - || strcmp(elements[itype],"S") == 0 || strcmp(elements[itype],"Se") == 0 - || strcmp(elements[itype],"Te") == 0) { + if (strcmp(elements[itype], "Mo") == 0 || strcmp(elements[itype], "W") == 0 || + strcmp(elements[itype], "S") == 0 || strcmp(elements[itype], "Se") == 0 || + strcmp(elements[itype], "Te") == 0) { // initialize neighsort - for (ll = 0; ll < n; ll++){ + for (ll = 0; ll < n; ll++) { neighsort[ll] = neighptr[ll]; check[ll] = neighptr[ll]; } // select the first neighbor of atomi - if(n == Nnei) { + if (n == Nnei) { neighsort[0] = neighptr[0]; check[0] = -1; - } - else if (n < Nnei && n > 0) { - for (jj = 0; jj < n; jj++) { //identify the first neighbor + } else if (n < Nnei && n > 0) { + for (jj = 0; jj < n; jj++) { //identify the first neighbor j = neighptr[jj]; jtype = map[type[j]]; count = 0; for (ll = 0; ll < n; ll++) { l = neighptr[ll]; ltype = map[type[l]]; - if(l == j) continue; + if (l == j) continue; deljx = x[l][0] - x[j][0]; deljy = x[l][1] - x[j][1]; deljz = x[l][2] - x[j][2]; - rsqlj = deljx*deljx + deljy*deljy + deljz*deljz; - if(rsqlj != 0 && rsqlj < cutILPsq[ltype][jtype]) { - count++; - } + rsqlj = deljx * deljx + deljy * deljy + deljz * deljz; + if (rsqlj != 0 && rsqlj < cutILPsq[ltype][jtype]) { count++; } } - if(count == 1) { + if (count == 1) { neighsort[0] = neighptr[jj]; check[jj] = -1; break; } - } // end of idenfying the first neighbor - } - else if (n > Nnei) { - fprintf (screen,"Molecule ID = %d\n",imol); - fprintf (screen,"Atom Type = %d\n",type[i]); - fprintf (screen,"Neinum = %d\n",n); - error->one(FLERR,"There are too many neighbors for TMD atoms, please check your configuration"); + } // end of idenfying the first neighbor + } else if (n > Nnei) { + fprintf(screen, "Molecule ID = %d\n", imol); + fprintf(screen, "Atom Type = %d\n", type[i]); + fprintf(screen, "Neinum = %d\n", n); + error->one(FLERR, + "There are too many neighbors for TMD atoms, please check your configuration"); } // sort the order of neighbors of atomi @@ -343,36 +349,38 @@ void PairILPTMD::ILP_neigh() j = neighsort[jj]; jtype = map[type[j]]; ll = 0; - while(ll < n) { + while (ll < n) { l = neighptr[ll]; - if(check[ll] == -1) {ll++; continue;} + if (check[ll] == -1) { + ll++; + continue; + } ltype = map[type[l]]; deljx = x[l][0] - x[j][0]; deljy = x[l][1] - x[j][1]; deljz = x[l][2] - x[j][2]; - rsqlj = deljx*deljx + deljy*deljy + deljz*deljz; - if(rsqlj != 0 && rsqlj < cutILPsq[ltype][jtype]) { - neighsort[jj+1] = l; + rsqlj = deljx * deljx + deljy * deljy + deljz * deljz; + if (rsqlj != 0 && rsqlj < cutILPsq[ltype][jtype]) { + neighsort[jj + 1] = l; check[ll] = -1; break; } ll++; } - } // end of sorting the order of neighbors - } - else {// for B/N/C/H atoms - if (n > 3) error->one(FLERR,"There are too many neighbors for B/N/C/H atoms, please check your configuration"); - for (ll = 0; ll < n; ll++) { - neighsort[ll] = neighptr[ll]; - } + } // end of sorting the order of neighbors + } else { // for B/N/C/H atoms + if (n > 3) + error->one( + FLERR, + "There are too many neighbors for B/N/C/H atoms, please check your configuration"); + for (ll = 0; ll < n; ll++) { neighsort[ll] = neighptr[ll]; } } ILP_firstneigh[i] = neighsort; ILP_numneigh[i] = n; ipage->vgot(n); - if (ipage->status()) - error->one(FLERR,"Neighbor list overflow, boost neigh_modify one"); + if (ipage->status()) error->one(FLERR, "Neighbor list overflow, boost neigh_modify one"); } } @@ -381,11 +389,11 @@ void PairILPTMD::ILP_neigh() ------------------------------------------------------------------------- */ void PairILPTMD::calc_normal() { - int i,j,ii,jj,inum,jnum; - int cont,id,ip,m,k,itype; - double nn,xtp,ytp,ztp,delx,dely,delz,nn2; - int *ilist,*jlist; - double Nave[3],dni[3],dpvdri[3][3]; + int i, j, ii, jj, inum, jnum; + int cont, id, ip, m, k, itype; + double nn, xtp, ytp, ztp, delx, dely, delz, nn2; + int *ilist, *jlist; + double Nave[3], dni[3], dpvdri[3][3]; double **x = atom->x; int *type = atom->type; @@ -396,12 +404,12 @@ void PairILPTMD::calc_normal() memory->destroy(dpvet1); memory->destroy(dpvet2); memory->destroy(dNave); - memory->create(dnn,Nnei,3,"ILPTMD:dnn"); - memory->create(vect,Nnei,3,"ILPTMD:vect"); - memory->create(pvet,Nnei,3,"ILPTMD:pvet"); - memory->create(dpvet1,Nnei,3,3,"ILPTMD:dpvet1"); - memory->create(dpvet2,Nnei,3,3,"ILPTMD:dpvet2"); - memory->create(dNave,3,Nnei,3,"ILPTMD:dNave"); + memory->create(dnn, Nnei, 3, "ILPTMD:dnn"); + memory->create(vect, Nnei, 3, "ILPTMD:vect"); + memory->create(pvet, Nnei, 3, "ILPTMD:pvet"); + memory->create(dpvet1, Nnei, 3, 3, "ILPTMD:dpvet1"); + memory->create(dpvet2, Nnei, 3, 3, "ILPTMD:dpvet2"); + memory->create(dNave, 3, Nnei, 3, "ILPTMD:dNave"); // grow normal array if necessary @@ -410,9 +418,9 @@ void PairILPTMD::calc_normal() memory->destroy(dnormal); memory->destroy(dnormdri); nmax = atom->nmax; - memory->create(normal,nmax,3,"ILPTMD:normal"); - memory->create(dnormdri,nmax,3,3,"ILPTMD:dnormdri"); - memory->create(dnormal,nmax,3,Nnei,3,"ILPTMD:dnormal"); + memory->create(normal, nmax, 3, "ILPTMD:normal"); + memory->create(dnormdri, nmax, 3, 3, "ILPTMD:dnormdri"); + memory->create(dnormal, nmax, 3, Nnei, 3, "ILPTMD:dnormal"); } inum = list->inum; @@ -426,14 +434,14 @@ void PairILPTMD::calc_normal() ztp = x[i][2]; // Initialize the arrays - for (id = 0; id < 3; id++){ + for (id = 0; id < 3; id++) { Nave[id] = 0.0; dni[id] = 0.0; normal[i][id] = 0.0; - for (ip = 0; ip < 3; ip++){ + for (ip = 0; ip < 3; ip++) { dpvdri[ip][id] = 0.0; dnormdri[i][ip][id] = 0.0; - for (m = 0; m < Nnei; m++){ + for (m = 0; m < Nnei; m++) { dnn[m][id] = 0.0; vect[m][id] = 0.0; pvet[m][id] = 0.0; @@ -461,89 +469,85 @@ void PairILPTMD::calc_normal() cont++; } -//############################ For the dangling atoms ############################ + //############################ For the dangling atoms ############################ if (cont <= 1) { normal[i][0] = 0.0; normal[i][1] = 0.0; normal[i][2] = 1.0; - for (id = 0; id < 3; id++){ - for (ip = 0; ip < 3; ip++){ + for (id = 0; id < 3; id++) { + for (ip = 0; ip < 3; ip++) { dnormdri[i][id][ip] = 0.0; - for (m = 0; m < Nnei; m++){ - dnormal[i][id][m][ip] = 0.0; - } + for (m = 0; m < Nnei; m++) { dnormal[i][id][m][ip] = 0.0; } } } } -//############################ For the edge atoms of TMD ################################ + //############################ For the edge atoms of TMD ################################ else if (cont > 1 && cont < Nnei) { - if(strcmp(elements[itype],"Mo") == 0|| strcmp(elements[itype],"W") == 0 || strcmp(elements[itype],"S") == 0|| strcmp(elements[itype],"Se") == 0) { + if (strcmp(elements[itype], "Mo") == 0 || strcmp(elements[itype], "W") == 0 || + strcmp(elements[itype], "S") == 0 || strcmp(elements[itype], "Se") == 0) { // derivatives of Ni[l] respect to the cont neighbors - for (k = 0; k < cont-1; k++){ - for (ip = 0; ip < 3; ip++){ - pvet[k][ip] = vect[k][modulo(ip+1,3)]*vect[k+1][modulo(ip+2,3)] - vect[k][modulo(ip+2,3)]*vect[k+1][modulo(ip+1,3)]; + for (k = 0; k < cont - 1; k++) { + for (ip = 0; ip < 3; ip++) { + pvet[k][ip] = vect[k][modulo(ip + 1, 3)] * vect[k + 1][modulo(ip + 2, 3)] - + vect[k][modulo(ip + 2, 3)] * vect[k + 1][modulo(ip + 1, 3)]; } - // dpvet1[k][l][ip]: the derivatve of the k (=0,...cont-1)th Nik respect to the ip component of atom l - // derivatives respect to atom l - // dNik,x/drl - dpvet1[k][0][0] = 0.0; - dpvet1[k][0][1] = vect[modulo(k+1,Nnei)][2]; - dpvet1[k][0][2] = -vect[modulo(k+1,Nnei)][1]; - // dNik,y/drl - dpvet1[k][1][0] = -vect[modulo(k+1,Nnei)][2]; - dpvet1[k][1][1] = 0.0; - dpvet1[k][1][2] = vect[modulo(k+1,Nnei)][0]; - // dNik,z/drl - dpvet1[k][2][0] = vect[modulo(k+1,Nnei)][1]; - dpvet1[k][2][1] = -vect[modulo(k+1,Nnei)][0]; - dpvet1[k][2][2] = 0.0; + // dpvet1[k][l][ip]: the derivatve of the k (=0,...cont-1)th Nik respect to the ip component of atom l + // derivatives respect to atom l + // dNik,x/drl + dpvet1[k][0][0] = 0.0; + dpvet1[k][0][1] = vect[modulo(k + 1, Nnei)][2]; + dpvet1[k][0][2] = -vect[modulo(k + 1, Nnei)][1]; + // dNik,y/drl + dpvet1[k][1][0] = -vect[modulo(k + 1, Nnei)][2]; + dpvet1[k][1][1] = 0.0; + dpvet1[k][1][2] = vect[modulo(k + 1, Nnei)][0]; + // dNik,z/drl + dpvet1[k][2][0] = vect[modulo(k + 1, Nnei)][1]; + dpvet1[k][2][1] = -vect[modulo(k + 1, Nnei)][0]; + dpvet1[k][2][2] = 0.0; - // dpvet2[k][l][ip]: the derivatve of the k (=0,...cont-1)th Nik respect to the ip component of atom l+1 - // derivatives respect to atom l+1 - // dNik,x/drl+1 - dpvet2[k][0][0] = 0.0; - dpvet2[k][0][1] = -vect[modulo(k,Nnei)][2]; - dpvet2[k][0][2] = vect[modulo(k,Nnei)][1]; - // dNik,y/drl+1 - dpvet2[k][1][0] = vect[modulo(k,Nnei)][2]; - dpvet2[k][1][1] = 0.0; - dpvet2[k][1][2] = -vect[modulo(k,Nnei)][0]; - // dNik,z/drl+1 - dpvet2[k][2][0] = -vect[modulo(k,Nnei)][1]; - dpvet2[k][2][1] = vect[modulo(k,Nnei)][0]; - dpvet2[k][2][2] = 0.0; + // dpvet2[k][l][ip]: the derivatve of the k (=0,...cont-1)th Nik respect to the ip component of atom l+1 + // derivatives respect to atom l+1 + // dNik,x/drl+1 + dpvet2[k][0][0] = 0.0; + dpvet2[k][0][1] = -vect[modulo(k, Nnei)][2]; + dpvet2[k][0][2] = vect[modulo(k, Nnei)][1]; + // dNik,y/drl+1 + dpvet2[k][1][0] = vect[modulo(k, Nnei)][2]; + dpvet2[k][1][1] = 0.0; + dpvet2[k][1][2] = -vect[modulo(k, Nnei)][0]; + // dNik,z/drl+1 + dpvet2[k][2][0] = -vect[modulo(k, Nnei)][1]; + dpvet2[k][2][1] = vect[modulo(k, Nnei)][0]; + dpvet2[k][2][2] = 0.0; } // average the normal vectors by using the Nnei neighboring planes - for (ip = 0; ip < 3; ip++){ + for (ip = 0; ip < 3; ip++) { Nave[ip] = 0.0; - for (k = 0; k < cont-1; k++){ - Nave[ip] += pvet[k][ip]; - } - Nave[ip] /= (cont - 1); + for (k = 0; k < cont - 1; k++) { Nave[ip] += pvet[k][ip]; } + Nave[ip] /= (cont - 1); } // the magnitude of the normal vector - nn2 = Nave[0]*Nave[0] + Nave[1]*Nave[1] + Nave[2]*Nave[2]; + nn2 = Nave[0] * Nave[0] + Nave[1] * Nave[1] + Nave[2] * Nave[2]; nn = sqrt(nn2); - if (nn == 0) error->one(FLERR,"The magnitude of the normal vector is zero"); + if (nn == 0) error->one(FLERR, "The magnitude of the normal vector is zero"); // the unit normal vector - normal[i][0] = Nave[0]/nn; - normal[i][1] = Nave[1]/nn; - normal[i][2] = Nave[2]/nn; + normal[i][0] = Nave[0] / nn; + normal[i][1] = Nave[1] / nn; + normal[i][2] = Nave[2] / nn; // derivatives of non-normalized normal vector, dNave:3xcontx3 array // dNave[id][m][ip]: the derivatve of the id component of Nave respect to the ip component of atom m - for (id = 0; id < 3; id++){ - for (ip = 0; ip < 3; ip++){ - for (m = 0; m < cont; m++){ - if(m == 0) { - dNave[id][m][ip] = dpvet1[m][id][ip]/(cont-1); - } - else if(m == cont - 1) { - dNave[id][m][ip] = dpvet2[m-1][id][ip]/(cont-1); - } - else { // sum of the derivatives of the mth and (m-1)th normal vector respect to the atom m - dNave[id][m][ip] = (dpvet1[m][id][ip] + dpvet2[m-1][id][ip])/(cont-1); + for (id = 0; id < 3; id++) { + for (ip = 0; ip < 3; ip++) { + for (m = 0; m < cont; m++) { + if (m == 0) { + dNave[id][m][ip] = dpvet1[m][id][ip] / (cont - 1); + } else if (m == cont - 1) { + dNave[id][m][ip] = dpvet2[m - 1][id][ip] / (cont - 1); + } else { // sum of the derivatives of the mth and (m-1)th normal vector respect to the atom m + dNave[id][m][ip] = (dpvet1[m][id][ip] + dpvet2[m - 1][id][ip]) / (cont - 1); } } } @@ -551,102 +555,99 @@ void PairILPTMD::calc_normal() // derivatives of nn, dnn:contx3 vector // dnn[m][id]: the derivative of nn respect to r[m][id], m=0,...Nnei-1; id=0,1,2 // r[m][id]: the id's component of atom m - for (m = 0; m < cont; m++){ - for (id = 0; id < 3; id++){ - dnn[m][id] = (Nave[0]*dNave[0][m][id] + Nave[1]*dNave[1][m][id] + Nave[2]*dNave[2][m][id])/nn; + for (m = 0; m < cont; m++) { + for (id = 0; id < 3; id++) { + dnn[m][id] = (Nave[0] * dNave[0][m][id] + Nave[1] * dNave[1][m][id] + + Nave[2] * dNave[2][m][id]) / + nn; } } // dnormal[i][id][m][ip]: the derivative of normal[i][id] respect to r[m][ip], id,ip=0,1,2. // for atom m, which is a neighbor atom of atom i, m = 0,...,Nnei-1 - for (m = 0; m < cont; m++){ - for (id = 0; id < 3; id++){ - for (ip = 0; ip < 3; ip++){ - dnormal[i][id][m][ip] = dNave[id][m][ip]/nn - Nave[id]*dnn[m][ip]/nn2; + for (m = 0; m < cont; m++) { + for (id = 0; id < 3; id++) { + for (ip = 0; ip < 3; ip++) { + dnormal[i][id][m][ip] = dNave[id][m][ip] / nn - Nave[id] * dnn[m][ip] / nn2; } } } // Calculte dNave/dri, defined as dpvdri - for (id = 0; id < 3; id++){ - for (ip = 0; ip < 3; ip++){ + for (id = 0; id < 3; id++) { + for (ip = 0; ip < 3; ip++) { dpvdri[id][ip] = 0.0; - for (k = 0; k < cont; k++){ - dpvdri[id][ip] -= dNave[id][k][ip]; - } + for (k = 0; k < cont; k++) { dpvdri[id][ip] -= dNave[id][k][ip]; } } } // derivatives of nn, dnn:3x1 vector - dni[0] = (Nave[0]*dpvdri[0][0] + Nave[1]*dpvdri[1][0] + Nave[2]*dpvdri[2][0])/nn; - dni[1] = (Nave[0]*dpvdri[0][1] + Nave[1]*dpvdri[1][1] + Nave[2]*dpvdri[2][1])/nn; - dni[2] = (Nave[0]*dpvdri[0][2] + Nave[1]*dpvdri[1][2] + Nave[2]*dpvdri[2][2])/nn; + dni[0] = (Nave[0] * dpvdri[0][0] + Nave[1] * dpvdri[1][0] + Nave[2] * dpvdri[2][0]) / nn; + dni[1] = (Nave[0] * dpvdri[0][1] + Nave[1] * dpvdri[1][1] + Nave[2] * dpvdri[2][1]) / nn; + dni[2] = (Nave[0] * dpvdri[0][2] + Nave[1] * dpvdri[1][2] + Nave[2] * dpvdri[2][2]) / nn; // derivatives of unit vector ni respect to ri, the result is 3x3 matrix - for (id = 0; id < 3; id++){ - for (ip = 0; ip < 3; ip++){ - dnormdri[i][id][ip] = dpvdri[id][ip]/nn - Nave[id]*dni[ip]/nn2; + for (id = 0; id < 3; id++) { + for (ip = 0; ip < 3; ip++) { + dnormdri[i][id][ip] = dpvdri[id][ip] / nn - Nave[id] * dni[ip] / nn2; } } - } // for TMD -//############################ For the edge & bulk atoms of GrhBN ################################ + } // for TMD + //############################ For the edge & bulk atoms of GrhBN ################################ else { - if(cont == 2) { - for (ip = 0; ip < 3; ip++){ - pvet[0][ip] = vect[0][modulo(ip+1,3)]*vect[1][modulo(ip+2,3)] - vect[0][modulo(ip+2,3)]*vect[1][modulo(ip+1,3)]; + if (cont == 2) { + for (ip = 0; ip < 3; ip++) { + pvet[0][ip] = vect[0][modulo(ip + 1, 3)] * vect[1][modulo(ip + 2, 3)] - + vect[0][modulo(ip + 2, 3)] * vect[1][modulo(ip + 1, 3)]; } // dpvet1[k][l][ip]: the derivatve of the k (=0,...cont-1)th Nik respect to the ip component of atom l // derivatives respect to atom l // dNik,x/drl - dpvet1[0][0][0] = 0.0; - dpvet1[0][0][1] = vect[1][2]; + dpvet1[0][0][0] = 0.0; + dpvet1[0][0][1] = vect[1][2]; dpvet1[0][0][2] = -vect[1][1]; // dNi0,y/drl dpvet1[0][1][0] = -vect[1][2]; - dpvet1[0][1][1] = 0.0; - dpvet1[0][1][2] = vect[1][0]; + dpvet1[0][1][1] = 0.0; + dpvet1[0][1][2] = vect[1][0]; // dNi0,z/drl - dpvet1[0][2][0] = vect[1][1]; + dpvet1[0][2][0] = vect[1][1]; dpvet1[0][2][1] = -vect[1][0]; - dpvet1[0][2][2] = 0.0; + dpvet1[0][2][2] = 0.0; // dpvet2[0][l][ip]: the derivatve of the 0 (=0,...cont-1)th Ni0 respect to the ip component of atom l+1 // derivatives respect to atom l+1 // dNi0,x/drl+1 - dpvet2[0][0][0] = 0.0; + dpvet2[0][0][0] = 0.0; dpvet2[0][0][1] = -vect[0][2]; - dpvet2[0][0][2] = vect[0][1]; + dpvet2[0][0][2] = vect[0][1]; // dNi0,y/drl+1 - dpvet2[0][1][0] = vect[0][2]; - dpvet2[0][1][1] = 0.0; + dpvet2[0][1][0] = vect[0][2]; + dpvet2[0][1][1] = 0.0; dpvet2[0][1][2] = -vect[0][0]; // dNi0,z/drl+1 dpvet2[0][2][0] = -vect[0][1]; - dpvet2[0][2][1] = vect[0][0]; - dpvet2[0][2][2] = 0.0; + dpvet2[0][2][1] = vect[0][0]; + dpvet2[0][2][2] = 0.0; - for (ip = 0; ip < 3; ip++){ - Nave[ip] += pvet[0][ip]; - } + for (ip = 0; ip < 3; ip++) { Nave[ip] += pvet[0][ip]; } // the magnitude of the normal vector - nn2 = Nave[0]*Nave[0] + Nave[1]*Nave[1] + Nave[2]*Nave[2]; + nn2 = Nave[0] * Nave[0] + Nave[1] * Nave[1] + Nave[2] * Nave[2]; nn = sqrt(nn2); - if (nn == 0) error->one(FLERR,"The magnitude of the normal vector is zero"); + if (nn == 0) error->one(FLERR, "The magnitude of the normal vector is zero"); // the unit normal vector - normal[i][0] = Nave[0]/nn; - normal[i][1] = Nave[1]/nn; - normal[i][2] = Nave[2]/nn; + normal[i][0] = Nave[0] / nn; + normal[i][1] = Nave[1] / nn; + normal[i][2] = Nave[2] / nn; // derivatives of non-normalized normal vector, dNave:3xcontx3 array // dNave[id][m][ip]: the derivatve of the id component of Nave respect to the ip component of atom m - for (id = 0; id < 3; id++){ - for (ip = 0; ip < 3; ip++){ - for (m = 0; m < cont; m++){ - if(m == 0) { - dNave[id][m][ip] = dpvet1[m][id][ip]/(cont-1); - } - else if(m == cont - 1) { - dNave[id][m][ip] = dpvet2[m-1][id][ip]/(cont-1); - } - else { // sum of the derivatives of the mth and (m-1)th normal vector respect to the atom m - dNave[id][m][ip] = (dpvet1[m][id][ip] + dpvet2[m-1][id][ip])/(cont-1); + for (id = 0; id < 3; id++) { + for (ip = 0; ip < 3; ip++) { + for (m = 0; m < cont; m++) { + if (m == 0) { + dNave[id][m][ip] = dpvet1[m][id][ip] / (cont - 1); + } else if (m == cont - 1) { + dNave[id][m][ip] = dpvet2[m - 1][id][ip] / (cont - 1); + } else { // sum of the derivatives of the mth and (m-1)th normal vector respect to the atom m + dNave[id][m][ip] = (dpvet1[m][id][ip] + dpvet2[m - 1][id][ip]) / (cont - 1); } } } @@ -654,224 +655,230 @@ void PairILPTMD::calc_normal() // derivatives of nn, dnn:contx3 vector // dnn[m][id]: the derivative of nn respect to r[m][id], m=0,...Nnei-1; id=0,1,2 // r[m][id]: the id's component of atom m - for (m = 0; m < cont; m++){ - for (id = 0; id < 3; id++){ - dnn[m][id] = (Nave[0]*dNave[0][m][id] + Nave[1]*dNave[1][m][id] + Nave[2]*dNave[2][m][id])/nn; + for (m = 0; m < cont; m++) { + for (id = 0; id < 3; id++) { + dnn[m][id] = (Nave[0] * dNave[0][m][id] + Nave[1] * dNave[1][m][id] + + Nave[2] * dNave[2][m][id]) / + nn; } } // dnormal[i][id][m][ip]: the derivative of normal[i][id] respect to r[m][ip], id,ip=0,1,2. // for atom m, which is a neighbor atom of atom i, m = 0,...,Nnei-1 - for (m = 0; m < cont; m++){ - for (id = 0; id < 3; id++){ - for (ip = 0; ip < 3; ip++){ - dnormal[i][id][m][ip] = dNave[id][m][ip]/nn - Nave[id]*dnn[m][ip]/nn2; + for (m = 0; m < cont; m++) { + for (id = 0; id < 3; id++) { + for (ip = 0; ip < 3; ip++) { + dnormal[i][id][m][ip] = dNave[id][m][ip] / nn - Nave[id] * dnn[m][ip] / nn2; } } } // Calculte dNave/dri, defined as dpvdri - for (id = 0; id < 3; id++){ - for (ip = 0; ip < 3; ip++){ + for (id = 0; id < 3; id++) { + for (ip = 0; ip < 3; ip++) { dpvdri[id][ip] = 0.0; - for (k = 0; k < cont; k++){ - dpvdri[id][ip] -= dNave[id][k][ip]; - } + for (k = 0; k < cont; k++) { dpvdri[id][ip] -= dNave[id][k][ip]; } } } // derivatives of nn, dnn:3x1 vector - dni[0] = (Nave[0]*dpvdri[0][0] + Nave[1]*dpvdri[1][0] + Nave[2]*dpvdri[2][0])/nn; - dni[1] = (Nave[0]*dpvdri[0][1] + Nave[1]*dpvdri[1][1] + Nave[2]*dpvdri[2][1])/nn; - dni[2] = (Nave[0]*dpvdri[0][2] + Nave[1]*dpvdri[1][2] + Nave[2]*dpvdri[2][2])/nn; + dni[0] = (Nave[0] * dpvdri[0][0] + Nave[1] * dpvdri[1][0] + Nave[2] * dpvdri[2][0]) / nn; + dni[1] = (Nave[0] * dpvdri[0][1] + Nave[1] * dpvdri[1][1] + Nave[2] * dpvdri[2][1]) / nn; + dni[2] = (Nave[0] * dpvdri[0][2] + Nave[1] * dpvdri[1][2] + Nave[2] * dpvdri[2][2]) / nn; // derivatives of unit vector ni respect to ri, the result is 3x3 matrix - for (id = 0; id < 3; id++){ - for (ip = 0; ip < 3; ip++){ - dnormdri[i][id][ip] = dpvdri[id][ip]/nn - Nave[id]*dni[ip]/nn2; + for (id = 0; id < 3; id++) { + for (ip = 0; ip < 3; ip++) { + dnormdri[i][id][ip] = dpvdri[id][ip] / nn - Nave[id] * dni[ip] / nn2; } } - } // end of cont == 2 - else if(cont == 3) { + } // end of cont == 2 + else if (cont == 3) { // derivatives of Ni[l] respect to the 3 neighbors - for (k = 0; k < 3; k++){ - for (ip = 0; ip < 3; ip++){ - pvet[k][ip] = vect[modulo(k,3)][modulo(ip+1,3)]*vect[modulo(k+1,3)][modulo(ip+2,3)] \ - -vect[modulo(k,3)][modulo(ip+2,3)]*vect[modulo(k+1,3)][modulo(ip+1,3)]; + for (k = 0; k < 3; k++) { + for (ip = 0; ip < 3; ip++) { + pvet[k][ip] = vect[modulo(k, 3)][modulo(ip + 1, 3)] * + vect[modulo(k + 1, 3)][modulo(ip + 2, 3)] - + vect[modulo(k, 3)][modulo(ip + 2, 3)] * vect[modulo(k + 1, 3)][modulo(ip + 1, 3)]; } - // dpvet1[k][l][ip]: the derivatve of the k (=0,...cont-1)th Nik respect to the ip component of atom l - // derivatives respect to atom l - // dNik,x/drl - dpvet1[k][0][0] = 0.0; - dpvet1[k][0][1] = vect[modulo(k+1,3)][2]; - dpvet1[k][0][2] = -vect[modulo(k+1,3)][1]; - // dNik,y/drl - dpvet1[k][1][0] = -vect[modulo(k+1,3)][2]; - dpvet1[k][1][1] = 0.0; - dpvet1[k][1][2] = vect[modulo(k+1,3)][0]; - // dNik,z/drl - dpvet1[k][2][0] = vect[modulo(k+1,3)][1]; - dpvet1[k][2][1] = -vect[modulo(k+1,3)][0]; - dpvet1[k][2][2] = 0.0; + // dpvet1[k][l][ip]: the derivatve of the k (=0,...cont-1)th Nik respect to the ip component of atom l + // derivatives respect to atom l + // dNik,x/drl + dpvet1[k][0][0] = 0.0; + dpvet1[k][0][1] = vect[modulo(k + 1, 3)][2]; + dpvet1[k][0][2] = -vect[modulo(k + 1, 3)][1]; + // dNik,y/drl + dpvet1[k][1][0] = -vect[modulo(k + 1, 3)][2]; + dpvet1[k][1][1] = 0.0; + dpvet1[k][1][2] = vect[modulo(k + 1, 3)][0]; + // dNik,z/drl + dpvet1[k][2][0] = vect[modulo(k + 1, 3)][1]; + dpvet1[k][2][1] = -vect[modulo(k + 1, 3)][0]; + dpvet1[k][2][2] = 0.0; - // dpvet2[k][l][ip]: the derivatve of the k (=0,...cont-1)th Nik respect to the ip component of atom l+1 - // derivatives respect to atom l+1 - // dNik,x/drl+1 - dpvet2[k][0][0] = 0.0; - dpvet2[k][0][1] = -vect[modulo(k,3)][2]; - dpvet2[k][0][2] = vect[modulo(k,3)][1]; - // dNik,y/drl+1 - dpvet2[k][1][0] = vect[modulo(k,3)][2]; - dpvet2[k][1][1] = 0.0; - dpvet2[k][1][2] = -vect[modulo(k,3)][0]; - // dNik,z/drl+1 - dpvet2[k][2][0] = -vect[modulo(k,3)][1]; - dpvet2[k][2][1] = vect[modulo(k,3)][0]; - dpvet2[k][2][2] = 0.0; + // dpvet2[k][l][ip]: the derivatve of the k (=0,...cont-1)th Nik respect to the ip component of atom l+1 + // derivatives respect to atom l+1 + // dNik,x/drl+1 + dpvet2[k][0][0] = 0.0; + dpvet2[k][0][1] = -vect[modulo(k, 3)][2]; + dpvet2[k][0][2] = vect[modulo(k, 3)][1]; + // dNik,y/drl+1 + dpvet2[k][1][0] = vect[modulo(k, 3)][2]; + dpvet2[k][1][1] = 0.0; + dpvet2[k][1][2] = -vect[modulo(k, 3)][0]; + // dNik,z/drl+1 + dpvet2[k][2][0] = -vect[modulo(k, 3)][1]; + dpvet2[k][2][1] = vect[modulo(k, 3)][0]; + dpvet2[k][2][2] = 0.0; } // average the normal vectors by using the 3 neighboring planes - for (ip = 0; ip < 3; ip++){ + for (ip = 0; ip < 3; ip++) { Nave[ip] = 0.0; - for (k = 0; k < 3; k++){ - Nave[ip] += pvet[k][ip]; - } - Nave[ip] /= 3; + for (k = 0; k < 3; k++) { Nave[ip] += pvet[k][ip]; } + Nave[ip] /= 3; } // the magnitude of the normal vector - nn2 = Nave[0]*Nave[0] + Nave[1]*Nave[1] + Nave[2]*Nave[2]; + nn2 = Nave[0] * Nave[0] + Nave[1] * Nave[1] + Nave[2] * Nave[2]; nn = sqrt(nn2); - if (nn == 0) error->one(FLERR,"The magnitude of the normal vector is zero"); + if (nn == 0) error->one(FLERR, "The magnitude of the normal vector is zero"); // the unit normal vector - normal[i][0] = Nave[0]/nn; - normal[i][1] = Nave[1]/nn; - normal[i][2] = Nave[2]/nn; + normal[i][0] = Nave[0] / nn; + normal[i][1] = Nave[1] / nn; + normal[i][2] = Nave[2] / nn; // for the central atoms, dnormdri is always zero - for (id = 0; id < 3; id++){ - for (ip = 0; ip < 3; ip++){ - dnormdri[i][id][ip] = 0.0; - } + for (id = 0; id < 3; id++) { + for (ip = 0; ip < 3; ip++) { dnormdri[i][id][ip] = 0.0; } } // derivatives of non-normalized normal vector, dNave:3x3x3 array // dNave[id][m][ip]: the derivatve of the id component of Nave respect to the ip component of atom m - for (id = 0; id < 3; id++){ - for (ip = 0; ip < 3; ip++){ - for (m = 0; m < 3; m++){ // sum of the derivatives of the mth and (m-1)th normal vector respect to the atom m - dNave[id][m][ip] = (dpvet1[modulo(m,3)][id][ip] + dpvet2[modulo(m-1,3)][id][ip])/3; + for (id = 0; id < 3; id++) { + for (ip = 0; ip < 3; ip++) { + for ( + m = 0; m < 3; + m++) { // sum of the derivatives of the mth and (m-1)th normal vector respect to the atom m + dNave[id][m][ip] = + (dpvet1[modulo(m, 3)][id][ip] + dpvet2[modulo(m - 1, 3)][id][ip]) / 3; } } } // derivatives of nn, dnn:3x3 vector // dnn[m][id]: the derivative of nn respect to r[m][id], m=0,...3-1; id=0,1,2 // r[m][id]: the id's component of atom m - for (m = 0; m < 3; m++){ - for (id = 0; id < 3; id++){ - dnn[m][id] = (Nave[0]*dNave[0][m][id] + Nave[1]*dNave[1][m][id] + Nave[2]*dNave[2][m][id])/nn; + for (m = 0; m < 3; m++) { + for (id = 0; id < 3; id++) { + dnn[m][id] = (Nave[0] * dNave[0][m][id] + Nave[1] * dNave[1][m][id] + + Nave[2] * dNave[2][m][id]) / + nn; } } // dnormal[i][id][m][ip]: the derivative of normal[i][id] respect to r[m][ip], id,ip=0,1,2. // for atom m, which is a neighbor atom of atom i, m = 0,...,3-1 - for (m = 0; m < 3; m++){ - for (id = 0; id < 3; id++){ - for (ip = 0; ip < 3; ip++){ - dnormal[i][id][m][ip] = dNave[id][m][ip]/nn - Nave[id]*dnn[m][ip]/nn2; + for (m = 0; m < 3; m++) { + for (id = 0; id < 3; id++) { + for (ip = 0; ip < 3; ip++) { + dnormal[i][id][m][ip] = dNave[id][m][ip] / nn - Nave[id] * dnn[m][ip] / nn2; } } } - } // end of cont == 3 - else error->one(FLERR,"There are too many neighbors for calculating normals of B/N/C/H atoms"); - } // for B/N/C/H - } // end of if(contone(FLERR, + "There are too many neighbors for calculating normals of B/N/C/H atoms"); + } // for B/N/C/H + } // end of if(contone(FLERR,"The magnitude of the normal vector is zero"); + if (nn == 0.0) error->one(FLERR, "The magnitude of the normal vector is zero"); // the unit normal vector - normal[i][0] = Nave[0]/nn; - normal[i][1] = Nave[1]/nn; - normal[i][2] = Nave[2]/nn; + normal[i][0] = Nave[0] / nn; + normal[i][1] = Nave[1] / nn; + normal[i][2] = Nave[2] / nn; // for the central atoms, dnormdri is always zero - for (id = 0; id < 3; id++){ - for (ip = 0; ip < 3; ip++){ - dnormdri[i][id][ip] = 0.0; - } + for (id = 0; id < 3; id++) { + for (ip = 0; ip < 3; ip++) { dnormdri[i][id][ip] = 0.0; } } // derivatives of non-normalized normal vector, dNave:3xNneix3 array // dNave[id][m][ip]: the derivatve of the id component of Nave respect to the ip component of atom m - for (id = 0; id < 3; id++){ - for (ip = 0; ip < 3; ip++){ - for (m = 0; m < Nnei; m++){ // sum of the derivatives of the mth and (m-1)th normal vector respect to the atom m - dNave[id][m][ip] = (dpvet1[modulo(m,Nnei)][id][ip] + dpvet2[modulo(m-1,Nnei)][id][ip])/Nnei; + for (id = 0; id < 3; id++) { + for (ip = 0; ip < 3; ip++) { + for ( + m = 0; m < Nnei; + m++) { // sum of the derivatives of the mth and (m-1)th normal vector respect to the atom m + dNave[id][m][ip] = + (dpvet1[modulo(m, Nnei)][id][ip] + dpvet2[modulo(m - 1, Nnei)][id][ip]) / Nnei; } } } // derivatives of nn, dnn:Nneix3 vector // dnn[m][id]: the derivative of nn respect to r[m][id], m=0,...Nnei-1; id=0,1,2 // r[m][id]: the id's component of atom m - for (m = 0; m < Nnei; m++){ - for (id = 0; id < 3; id++){ - dnn[m][id] = (Nave[0]*dNave[0][m][id] + Nave[1]*dNave[1][m][id] + Nave[2]*dNave[2][m][id])/nn; + for (m = 0; m < Nnei; m++) { + for (id = 0; id < 3; id++) { + dnn[m][id] = + (Nave[0] * dNave[0][m][id] + Nave[1] * dNave[1][m][id] + Nave[2] * dNave[2][m][id]) / + nn; } } // dnormal[i][id][m][ip]: the derivative of normal[i][id] respect to r[m][ip], id,ip=0,1,2. // for atom m, which is a neighbor atom of atom i, m = 0,...,Nnei-1 - for (m = 0; m < Nnei; m++){ - for (id = 0; id < 3; id++){ - for (ip = 0; ip < 3; ip++){ - dnormal[i][id][m][ip] = dNave[id][m][ip]/nn - Nave[id]*dnn[m][ip]/nn2; + for (m = 0; m < Nnei; m++) { + for (id = 0; id < 3; id++) { + for (ip = 0; ip < 3; ip++) { + dnormal[i][id][m][ip] = dNave[id][m][ip] / nn - Nave[id] * dnn[m][ip] / nn2; } } } - } - else { - error->one(FLERR,"There are too many neighbors for calculating normals of TMD atoms"); - } // end of four cases of cont - } // end of i loop + } else { + error->one(FLERR, "There are too many neighbors for calculating normals of TMD atoms"); + } // end of four cases of cont + } // end of i loop } diff --git a/src/INTERLAYER/pair_ilp_tmd.h b/src/INTERLAYER/pair_ilp_tmd.h index ac8d86a253..7ada9d59fd 100644 --- a/src/INTERLAYER/pair_ilp_tmd.h +++ b/src/INTERLAYER/pair_ilp_tmd.h @@ -43,7 +43,6 @@ class PairILPTMD : public PairILPGrapheneHBN { return k % range; } /**************************************************************/ - }; } // namespace LAMMPS_NS diff --git a/src/INTERLAYER/pair_kolmogorov_crespi_z.cpp b/src/INTERLAYER/pair_kolmogorov_crespi_z.cpp index 96b9c5b2a8..bd2ce680d5 100644 --- a/src/INTERLAYER/pair_kolmogorov_crespi_z.cpp +++ b/src/INTERLAYER/pair_kolmogorov_crespi_z.cpp @@ -28,8 +28,8 @@ #include "error.h" #include "force.h" #include "memory.h" -#include "neighbor.h" #include "neigh_list.h" +#include "neighbor.h" #include "potential_file_reader.h" #include @@ -160,15 +160,15 @@ void PairKolmogorovCrespiZ::compute(int eflag, int vflag) if (evflag) { ev_tally(i, j, nlocal, newton_pair, evdwl, 0.0, fpair, delx, dely, delz); if (vflag_either) { - double fi[3],fj[3]; + double fi[3], fj[3]; fi[0] = delx * fpair1; fi[1] = dely * fpair1; fi[2] = 0; fj[0] = -delx * fpair1; fj[1] = -dely * fpair1; fj[2] = 0; - v_tally2_newton(i,fi,x[i]); - v_tally2_newton(j,fj,x[j]); + v_tally2_newton(i, fi, x[i]); + v_tally2_newton(j, fj, x[j]); } } } @@ -246,9 +246,9 @@ void PairKolmogorovCrespiZ::coeff(int narg, char **arg) void PairKolmogorovCrespiZ::init_style() { if (force->newton_pair == 0) - error->all(FLERR,"Pair style kolmogorov/crespi/z requires newton pair on"); + error->all(FLERR, "Pair style kolmogorov/crespi/z requires newton pair on"); - neighbor->request(this,instance_me); + neighbor->request(this, instance_me); } /* ---------------------------------------------------------------------- diff --git a/src/INTERLAYER/pair_saip_metal.cpp b/src/INTERLAYER/pair_saip_metal.cpp index 02a5c9f323..f18b2b8ee4 100644 --- a/src/INTERLAYER/pair_saip_metal.cpp +++ b/src/INTERLAYER/pair_saip_metal.cpp @@ -52,7 +52,8 @@ static const char cite_saip[] = /* ---------------------------------------------------------------------- */ -PairSAIPMETAL::PairSAIPMETAL(LAMMPS *lmp) : PairILPGrapheneHBN(lmp) { +PairSAIPMETAL::PairSAIPMETAL(LAMMPS *lmp) : PairILPGrapheneHBN(lmp) +{ variant = SAIP_METAL; if (lmp->citeme) lmp->citeme->add(cite_saip); } @@ -71,19 +72,18 @@ void PairSAIPMETAL::settings(int narg, char **arg) if (narg == 2) tap_flag = utils::numeric(FLERR, arg[1], false, lmp); } - /* ---------------------------------------------------------------------- Repulsive forces and energy ------------------------------------------------------------------------- */ void PairSAIPMETAL::calc_FRep(int eflag, int /* vflag */) { - int i,j,ii,jj,inum,jnum,itype,jtype,k,kk; - double prodnorm1,fkcx,fkcy,fkcz,filp; - double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair,fpair1; - double rsq,r,Rcut,rhosq1,exp0,exp1,Tap,dTap,Vilp; - double frho1,Erep,fsum,rdsq1; - int *ilist,*jlist,*numneigh,**firstneigh; + int i, j, ii, jj, inum, jnum, itype, jtype, k, kk; + double prodnorm1, fkcx, fkcy, fkcz, filp; + double xtmp, ytmp, ztmp, delx, dely, delz, evdwl, fpair, fpair1; + double rsq, r, Rcut, rhosq1, exp0, exp1, Tap, dTap, Vilp; + double frho1, Erep, fsum, rdsq1; + int *ilist, *jlist, *numneigh, **firstneigh; int *ILP_neighs_i; evdwl = 0.0; @@ -123,35 +123,38 @@ void PairSAIPMETAL::calc_FRep(int eflag, int /* vflag */) delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; delz = ztmp - x[j][2]; - rsq = delx*delx + dely*dely + delz*delz; + rsq = delx * delx + dely * dely + delz * delz; // only include the interaction between different layers if (rsq < cutsq[itype][jtype] && atom->molecule[i] != atom->molecule[j]) { int iparam_ij = elem2param[map[itype]][map[jtype]]; - Param& p = params[iparam_ij]; + Param &p = params[iparam_ij]; r = sqrt(rsq); // turn on/off taper function if (tap_flag) { Rcut = sqrt(cutsq[itype][jtype]); - Tap = calc_Tap(r,Rcut); - dTap = calc_dTap(r,Rcut); - } else {Tap = 1.0; dTap = 0.0;} + Tap = calc_Tap(r, Rcut); + dTap = calc_dTap(r, Rcut); + } else { + Tap = 1.0; + dTap = 0.0; + } // for atoms in bulk materials - if (strcmp(elements[map[itype]],"C") != 0 && strcmp(elements[map[itype]],"H") != 0 && - strcmp(elements[map[itype]],"B") != 0 && strcmp(elements[map[itype]],"N") != 0) { + if (strcmp(elements[map[itype]], "C") != 0 && strcmp(elements[map[itype]], "H") != 0 && + strcmp(elements[map[itype]], "B") != 0 && strcmp(elements[map[itype]], "N") != 0) { // Set ni along rij - exp0 = exp(-p.lambda*(r-p.z0)); - frho1 = 1.0*p.C; - Erep = 0.5*p.epsilon + frho1; - Vilp = exp0*Erep; + exp0 = exp(-p.lambda * (r - p.z0)); + frho1 = 1.0 * p.C; + Erep = 0.5 * p.epsilon + frho1; + Vilp = exp0 * Erep; // derivatives - fpair = p.lambda*exp0/r*Erep; - filp = fpair*Tap - Vilp*dTap/r; - fkcx = delx*filp; - fkcy = dely*filp; - fkcz = delz*filp; + fpair = p.lambda * exp0 / r * Erep; + filp = fpair * Tap - Vilp * dTap / r; + fkcx = delx * filp; + fkcy = dely * filp; + fkcz = delz * filp; f[i][0] += fkcx; f[i][1] += fkcy; @@ -159,46 +162,48 @@ void PairSAIPMETAL::calc_FRep(int eflag, int /* vflag */) f[j][0] -= fkcx; f[j][1] -= fkcy; f[j][2] -= fkcz; - - if (eflag) pvector[1] += evdwl = Tap*Vilp; - if (evflag) ev_tally(i,j,nlocal,newton_pair,evdwl,0.0,filp,delx,dely,delz); - } - else { // for atoms in 2D materials + + if (eflag) pvector[1] += evdwl = Tap * Vilp; + if (evflag) ev_tally(i, j, nlocal, newton_pair, evdwl, 0.0, filp, delx, dely, delz); + } else { // for atoms in 2D materials // Calculate the transverse distance - prodnorm1 = normal[i][0]*delx + normal[i][1]*dely + normal[i][2]*delz; - rhosq1 = rsq - prodnorm1*prodnorm1; // rho_ij - rdsq1 = rhosq1*p.delta2inv; // (rho_ij/delta)^2 + prodnorm1 = normal[i][0] * delx + normal[i][1] * dely + normal[i][2] * delz; + rhosq1 = rsq - prodnorm1 * prodnorm1; // rho_ij + rdsq1 = rhosq1 * p.delta2inv; // (rho_ij/delta)^2 // store exponents - exp0 = exp(-p.lambda*(r-p.z0)); + exp0 = exp(-p.lambda * (r - p.z0)); exp1 = exp(-rdsq1); - frho1 = exp1*p.C; - Erep = 0.5*p.epsilon + frho1; - Vilp = exp0*Erep; + frho1 = exp1 * p.C; + Erep = 0.5 * p.epsilon + frho1; + Vilp = exp0 * Erep; // derivatives - fpair = p.lambda*exp0/r*Erep; - fpair1 = 2.0*exp0*frho1*p.delta2inv; + fpair = p.lambda * exp0 / r * Erep; + fpair1 = 2.0 * exp0 * frho1 * p.delta2inv; fsum = fpair + fpair1; // derivatives of the product of rij and ni, the result is a vector - 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; + 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; - fkcx = (delx*fsum - fp1[0])*Tap - Vilp*dTap*delx/r; - fkcy = (dely*fsum - fp1[1])*Tap - Vilp*dTap*dely/r; - fkcz = (delz*fsum - fp1[2])*Tap - Vilp*dTap*delz/r; + fkcx = (delx * fsum - fp1[0]) * Tap - Vilp * dTap * delx / r; + fkcy = (dely * fsum - fp1[1]) * Tap - Vilp * dTap * dely / r; + fkcz = (delz * fsum - fp1[2]) * Tap - Vilp * dTap * delz / r; - f[i][0] += fkcx - fprod1[0]*Tap; - f[i][1] += fkcy - fprod1[1]*Tap; - f[i][2] += fkcz - fprod1[2]*Tap; + f[i][0] += fkcx - fprod1[0] * Tap; + f[i][1] += fkcy - fprod1[1] * Tap; + f[i][2] += fkcz - fprod1[2] * Tap; f[j][0] -= fkcx; f[j][1] -= fkcy; f[j][2] -= fkcz; @@ -209,25 +214,31 @@ void PairSAIPMETAL::calc_FRep(int eflag, int /* vflag */) k = ILP_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; - dprodnorm1[1] = dnormal[0][1][kk][i]*delx + dnormal[1][1][kk][i]*dely + dnormal[2][1][kk][i]*delz; - dprodnorm1[2] = dnormal[0][2][kk][i]*delx + dnormal[1][2][kk][i]*dely + dnormal[2][2][kk][i]*delz; - fk[0] = (-prodnorm1*dprodnorm1[0]*fpair1)*Tap; - fk[1] = (-prodnorm1*dprodnorm1[1]*fpair1)*Tap; - fk[2] = (-prodnorm1*dprodnorm1[2]*fpair1)*Tap; + dprodnorm1[0] = dnormal[0][0][kk][i] * delx + dnormal[1][0][kk][i] * dely + + dnormal[2][0][kk][i] * delz; + dprodnorm1[1] = dnormal[0][1][kk][i] * delx + dnormal[1][1][kk][i] * dely + + dnormal[2][1][kk][i] * delz; + dprodnorm1[2] = dnormal[0][2][kk][i] * delx + dnormal[1][2][kk][i] * dely + + dnormal[2][2][kk][i] * delz; + fk[0] = (-prodnorm1 * dprodnorm1[0] * fpair1) * Tap; + fk[1] = (-prodnorm1 * dprodnorm1[1] * fpair1) * Tap; + fk[2] = (-prodnorm1 * dprodnorm1[2] * fpair1) * Tap; f[k][0] += fk[0]; f[k][1] += fk[1]; f[k][2] += fk[2]; delki[0] = x[k][0] - x[i][0]; delki[1] = x[k][1] - x[i][1]; delki[2] = x[k][2] - x[i][2]; - if (evflag) ev_tally_xyz(k,j,nlocal,newton_pair,0.0,0.0,fk[0],fk[1],fk[2],delki[0],delki[1],delki[2]); + if (evflag) + ev_tally_xyz(k, j, nlocal, newton_pair, 0.0, 0.0, fk[0], fk[1], fk[2], delki[0], + delki[1], delki[2]); } - - if (eflag) pvector[1] += evdwl = Tap*Vilp; - if (evflag) ev_tally_xyz(i,j,nlocal,newton_pair,evdwl,0.0,fkcx,fkcy,fkcz,delx,dely,delz); - } // end of speration atoms for bulk materials - } // end of rsq < cutoff - } // loop over jj - } // loop over ii + + if (eflag) pvector[1] += evdwl = Tap * Vilp; + if (evflag) + ev_tally_xyz(i, j, nlocal, newton_pair, evdwl, 0.0, fkcx, fkcy, fkcz, delx, dely, delz); + } // end of speration atoms for bulk materials + } // end of rsq < cutoff + } // loop over jj + } // loop over ii }