ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

fix permissions, error docs, URLs, whitespace

Browse Source
This commit is contained in:
Axel Kohlmeyer
2022-10-18 11:53:38 -04:00
parent 72e859acf8
commit 39020379da
6 changed files with 68 additions and 887133 deletions
Download Patch File Download Diff File Expand all files Collapse all files

12
doc/src/pair_ylz.rst Executable file → Normal file
View File

@ -2,7 +2,7 @@
.. index:: pair_style ylz/gpu .. index:: pair_style ylz/gpu
.. index:: pair_style ylz/intel .. index:: pair_style ylz/intel
.. index:: pair_style ylz/omp .. index:: pair_style ylz/omp
pair_style ylc command pair_style ylc command
=========================== ===========================
@ -24,18 +24,18 @@ Examples
.. code-block:: LAMMPS .. code-block:: LAMMPS
pair_style ylz 2.6 pair_style ylz 2.6
pair_coeff * * 1.0 1.0 4 3 0.0 2.6 pair_coeff * * 1.0 1.0 4 3 0.0 2.6
Description Description
""""""""""" """""""""""
The *ylz* (Yuan-Li-Zhang) The *ylz* (Yuan-Li-Zhang)
:ref:`(Yuan) <Yuan>` style computes anisotropic interactions between pairs of particles considering the relative particle orientations via the formulas :ref:`(Yuan) <Yuan>` style computes anisotropic interactions between pairs of particles considering the relative particle orientations via the formulas
.. math:: .. math::
U ( \mathbf{r}_{ij}, \mathbf{n}_i, \mathbf{n}_j ) =\left\{\begin{matrix} \mathbf{u}_R(r)+\left [ 1-\phi (\mathbf{r\hat{}}_{ij}, \mathbf{n}_i, \mathbf{n}_j ) \right ]\epsilon, ~~ r<\mathbf{r}_{min} \\ \mathbf{u}_A(r)\phi (\mathbf{r\hat{}}_{ij}, \mathbf{n}_i, \mathbf{n}_j ),~~ \mathbf{r}_{min}<r<\mathbf{r}_{c} \\ \end{matrix}\right.\\\\ \phi (\mathbf{r\hat{}}_{ij}, \mathbf{n}_i, \mathbf{n}_j )=1+\mu (a(\mathbf{r\hat{}}_{ij}, \mathbf{n}_i, \mathbf{n}_j )-1) \\\\ a(\mathbf{r\hat{}}_{ij}, \mathbf{n}_i, \mathbf{n}_j )=(\mathbf{n}_i\times\mathbf{r\hat{}}_{ij} )\cdot (\mathbf{n}_j\times\mathbf{r\hat{}}_{ij} )+sin\mathbf{\theta}_0(\mathbf{n}_i-\mathbf{n}_j)\cdot \mathbf{r\hat{}}_{ij}\\\\ \mathbf{u}_R(r)=\epsilon \left [ \left ( \frac{{r}_{min}}{r} \right )^{4}-2\left ( \frac{{r}_{min}}{r}\right )^{2} \right ] \\\\ \mathbf{u}_A(r)=-\epsilon\;cos^{2\zeta }\left ( \frac{\pi}{2}\frac{\left ( {r}-{r}_{min} \right )}{\left ( {r}_{c}-{r}_{min} \right )} \right ) \\ U ( \mathbf{r}_{ij}, \mathbf{n}_i, \mathbf{n}_j ) =\left\{\begin{matrix} \mathbf{u}_R(r)+\left [ 1-\phi (\mathbf{r\hat{}}_{ij}, \mathbf{n}_i, \mathbf{n}_j ) \right ]\epsilon, ~~ r<\mathbf{r}_{min} \\ \mathbf{u}_A(r)\phi (\mathbf{r\hat{}}_{ij}, \mathbf{n}_i, \mathbf{n}_j ),~~ \mathbf{r}_{min}<r<\mathbf{r}_{c} \\ \end{matrix}\right.\\\\ \phi (\mathbf{r\hat{}}_{ij}, \mathbf{n}_i, \mathbf{n}_j )=1+\mu (a(\mathbf{r\hat{}}_{ij}, \mathbf{n}_i, \mathbf{n}_j )-1) \\\\ a(\mathbf{r\hat{}}_{ij}, \mathbf{n}_i, \mathbf{n}_j )=(\mathbf{n}_i\times\mathbf{r\hat{}}_{ij} )\cdot (\mathbf{n}_j\times\mathbf{r\hat{}}_{ij} )+sin\mathbf{\theta}_0(\mathbf{n}_i-\mathbf{n}_j)\cdot \mathbf{r\hat{}}_{ij}\\\\ \mathbf{u}_R(r)=\epsilon \left [ \left ( \frac{{r}_{min}}{r} \right )^{4}-2\left ( \frac{{r}_{min}}{r}\right )^{2} \right ] \\\\ \mathbf{u}_A(r)=-\epsilon\;cos^{2\zeta }\left ( \frac{\pi}{2}\frac{\left ( {r}-{r}_{min} \right )}{\left ( {r}_{c}-{r}_{min} \right )} \right ) \\
where :math:`r_{i}` and :math:`r_{j}` are the center position vectors of particles i and j, respectively, :math:`r_{ij}=r_{i}-r_{j}` is the inter-particle distance vector, :math:`r=\left|r_{ij} \right|` and :math:`{r\hat{}}_{ij}=\mathbf{r}_{ij}/r`. The unit vectors :math:`n_{i}` and :math:`n_{j}` represent the axes of symmetry of particles i and j, respectively. :math:`u_R` and :math:`u_A` are the repulsive and attractive potentials. :math:`\phi` is an angular function which depends on the relative orientation between pair particles. :math:`\mu` is the parameter related to bending rigidity, :math:`\theta_{0}` is the parameter related to the spontaneous curvature, and :math:`\epsilon` is the energy unit, respectively. The :math:`\zeta` controls the slope of the attractive branch and :math:`{r}_{c}`is the cutoff radius. :math:`r_{min}` is the distance which minimizes the potential energy :math:`u_{A}(r)`and :math:`r_{min}=2^{1/6}\sigma`, where :math:`\sigma` is the length unit. where :math:`r_{i}` and :math:`r_{j}` are the center position vectors of particles i and j, respectively, :math:`r_{ij}=r_{i}-r_{j}` is the inter-particle distance vector, :math:`r=\left|r_{ij} \right|` and :math:`{r\hat{}}_{ij}=\mathbf{r}_{ij}/r`. The unit vectors :math:`n_{i}` and :math:`n_{j}` represent the axes of symmetry of particles i and j, respectively. :math:`u_R` and :math:`u_A` are the repulsive and attractive potentials. :math:`\phi` is an angular function which depends on the relative orientation between pair particles. :math:`\mu` is the parameter related to bending rigidity, :math:`\theta_{0}` is the parameter related to the spontaneous curvature, and :math:`\epsilon` is the energy unit, respectively. The :math:`\zeta` controls the slope of the attractive branch and :math:`{r}_{c}`is the cutoff radius. :math:`r_{min}` is the distance which minimizes the potential energy :math:`u_{A}(r)`and :math:`r_{min}=2^{1/6}\sigma`, where :math:`\sigma` is the length unit.

0
examples/ASPHERE/flat_membrane/in_flat_membrane Executable file → Normal file
View File

887047
examples/ASPHERE/vesicle/onlymembrane2.lammpstrj
View File

File diff suppressed because it is too large Load Diff

0
examples/ASPHERE/vesicle/read_data.vesicle1026 Executable file → Normal file
View File

116
src/ASPHERE/pair_ylz.cpp Executable file → Normal file
View File

@ -86,7 +86,7 @@ void PairYLZ::compute(int eflag, int vflag)
evdwl = 0.0; evdwl = 0.0;
ev_init(eflag , vflag); ev_init(eflag , vflag);
AtomVecEllipsoid::Bonus *bonus = avec->bonus; AtomVecEllipsoid::Bonus *bonus = avec->bonus;
int *ellipsoid = atom->ellipsoid; int *ellipsoid = atom->ellipsoid;
double **x = atom->x; double **x = atom->x;
@ -111,7 +111,7 @@ void PairYLZ::compute(int eflag, int vflag)
iquat = bonus[ellipsoid[i]].quat; iquat = bonus[ellipsoid[i]].quat;
MathExtra::quat_to_mat_trans(iquat,a1); MathExtra::quat_to_mat_trans(iquat,a1);
jlist = firstneigh[i]; jlist = firstneigh[i];
jnum = numneigh[i]; jnum = numneigh[i];
@ -133,7 +133,7 @@ void PairYLZ::compute(int eflag, int vflag)
if (rsq < cutsq[itype][jtype]) { if (rsq < cutsq[itype][jtype]) {
jquat = bonus[ellipsoid[j]].quat; jquat = bonus[ellipsoid[j]].quat;
MathExtra::quat_to_mat_trans(jquat,a2); MathExtra::quat_to_mat_trans(jquat,a2);
one_eng = ylz_analytic(i,j,a1,a2,r12,rsq,fforce,ttor,rtor); one_eng = ylz_analytic(i,j,a1,a2,r12,rsq,fforce,ttor,rtor);
fforce[0] *= factor_lj; fforce[0] *= factor_lj;
@ -172,8 +172,8 @@ void PairYLZ::compute(int eflag, int vflag)
} }
if (vflag_fdotr) virial_fdotr_compute(); if (vflag_fdotr) virial_fdotr_compute();
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
@ -197,7 +197,7 @@ void PairYLZ::allocate()
memory->create(zeta,n+1,n+1,"pair:zeta"); memory->create(zeta,n+1,n+1,"pair:zeta");
memory->create(mu,n+1,n+1,"pair:mu"); memory->create(mu,n+1,n+1,"pair:mu");
memory->create(beta,n+1,n+1,"pair:beta"); memory->create(beta,n+1,n+1,"pair:beta");
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
@ -207,7 +207,7 @@ void PairYLZ::allocate()
void PairYLZ::settings(int narg, char **arg) void PairYLZ::settings(int narg, char **arg)
{ {
if (narg != 1) error->all(FLERR,"Illegal pair_style command"); if (narg != 1) error->all(FLERR,"Illegal pair_style command");
cut_global = utils::numeric(FLERR,arg[0],false, lmp); cut_global = utils::numeric(FLERR,arg[0],false, lmp);
} }
@ -232,7 +232,7 @@ void PairYLZ::coeff(int narg, char **arg)
double mu_one = utils::numeric(FLERR,arg[5], false, lmp); double mu_one = utils::numeric(FLERR,arg[5], false, lmp);
double beta_one = utils::numeric(FLERR,arg[6], false, lmp); double beta_one = utils::numeric(FLERR,arg[6], false, lmp);
double cut_one = utils::numeric(FLERR,arg[7], false, lmp); double cut_one = utils::numeric(FLERR,arg[7], false, lmp);
int count = 0; int count = 0;
for (int i = ilo; i <= ihi; i++) { for (int i = ilo; i <= ihi; i++) {
for (int j = MAX(jlo,i); j <= jhi; j++) { for (int j = MAX(jlo,i); j <= jhi; j++) {
@ -262,7 +262,7 @@ void PairYLZ::init_style()
if (!avec) error->all(FLERR,"Pair YLZ requires atom style ellipsoid"); if (!avec) error->all(FLERR,"Pair YLZ requires atom style ellipsoid");
neighbor->request(this,instance_me); neighbor->request(this,instance_me);
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
@ -299,7 +299,7 @@ void PairYLZ::write_restart(FILE *fp)
int i,j; int i,j;
for (i = 1; i <= atom->ntypes; i++) { for (i = 1; i <= atom->ntypes; i++) {
for (j = i; j <= atom->ntypes; j++) { for (j = i; j <= atom->ntypes; j++) {
fwrite(&setflag[i][j],sizeof(int),1,fp); fwrite(&setflag[i][j],sizeof(int),1,fp);
if (setflag[i][j]) { if (setflag[i][j]) {
@ -325,10 +325,10 @@ void PairYLZ::read_restart(FILE *fp)
int i,j; int i,j;
int me = comm->me; int me = comm->me;
for (i = 1; i <= atom->ntypes; i++) { for (i = 1; i <= atom->ntypes; i++) {
for (j = i; j <= atom->ntypes; j++) { for (j = i; j <= atom->ntypes; j++) {
if (me == 0) utils::sfread(FLERR, &setflag[i][j],sizeof(int),1,fp, nullptr, error); if (me == 0) utils::sfread(FLERR, &setflag[i][j],sizeof(int),1,fp, nullptr, error);
MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world); MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world);
@ -340,17 +340,17 @@ void PairYLZ::read_restart(FILE *fp)
utils::sfread(FLERR, &zeta[i][j],sizeof(double),1,fp, nullptr, error); utils::sfread(FLERR, &zeta[i][j],sizeof(double),1,fp, nullptr, error);
utils::sfread(FLERR, &mu[i][j],sizeof(double),1,fp, nullptr, error); utils::sfread(FLERR, &mu[i][j],sizeof(double),1,fp, nullptr, error);
utils::sfread(FLERR, &beta[i][j],sizeof(double),1,fp, nullptr, error); utils::sfread(FLERR, &beta[i][j],sizeof(double),1,fp, nullptr, error);
} }
MPI_Bcast(&epsilon[i][j],1,MPI_DOUBLE,0,world); MPI_Bcast(&epsilon[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&sigma[i][j],1,MPI_DOUBLE,0,world); MPI_Bcast(&sigma[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world); MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&zeta[i][j],1,MPI_DOUBLE,0,world); MPI_Bcast(&zeta[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&mu[i][j],1,MPI_DOUBLE,0,world); MPI_Bcast(&mu[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&beta[i][j],1,MPI_DOUBLE,0,world); MPI_Bcast(&beta[i][j],1,MPI_DOUBLE,0,world);
} }
} }
} }
@ -376,13 +376,13 @@ void PairYLZ::read_restart_settings(FILE *fp)
{ {
int me = comm->me; int me = comm->me;
if (me == 0) { if (me == 0) {
// utils::sfread(FLERR, &model_type,sizeof(int),1,fp, nullptr, error); // utils::sfread(FLERR, &model_type,sizeof(int),1,fp, nullptr, error);
utils::sfread(FLERR, &cut_global,sizeof(double),1,fp, nullptr, error); utils::sfread(FLERR, &cut_global,sizeof(double),1,fp, nullptr, error);
utils::sfread(FLERR, &offset_flag,sizeof(int),1,fp, nullptr, error); utils::sfread(FLERR, &offset_flag,sizeof(int),1,fp, nullptr, error);
utils::sfread(FLERR, &mix_flag,sizeof(int),1,fp, nullptr, error); utils::sfread(FLERR, &mix_flag,sizeof(int),1,fp, nullptr, error);
} }
//MPI_Bcast(&model_type,1,MPI_DOUBLE,0,world); //MPI_Bcast(&model_type,1,MPI_DOUBLE,0,world);
MPI_Bcast(&cut_global,1,MPI_DOUBLE,0,world); MPI_Bcast(&cut_global,1,MPI_DOUBLE,0,world);
MPI_Bcast(&offset_flag,1,MPI_INT,0,world); MPI_Bcast(&offset_flag,1,MPI_INT,0,world);
@ -417,7 +417,7 @@ void PairYLZ::write_data_all(FILE *fp)
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
compute analytic energy, force (fforce), and torque (ttor & rtor) compute analytic energy, force (fforce), and torque (ttor & rtor)
based on rotation matrices a based on rotation matrices a
if newton is off, rtor is not calculated for ghost atoms if newton is off, rtor is not calculated for ghost atoms
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
@ -426,7 +426,7 @@ double PairYLZ::ylz_analytic(const int i,const int j,double a1[3][3],
const double rsq, double *fforce, const double rsq, double *fforce,
double *ttor, double *rtor) double *ttor, double *rtor)
{ {
int *type = atom->type; int *type = atom->type;
int newton_pair = force->newton_pair; int newton_pair = force->newton_pair;
int nlocal = atom->nlocal; int nlocal = atom->nlocal;
@ -436,43 +436,43 @@ double PairYLZ::ylz_analytic(const int i,const int j,double a1[3][3],
double r = sqrt(rsq); double r = sqrt(rsq);
double ni1[3],nj1[3],dphi_drhat[3],dUdrhat[3],dUdni1[3],dUdnj1[3],dphi_dni1[3],dphi_dnj1[3] ; double ni1[3],nj1[3],dphi_drhat[3],dUdrhat[3],dUdni1[3],dUdnj1[3],dphi_dni1[3],dphi_dnj1[3] ;
double t,t1,t2,t4,cos_t,U,uR,uA,dUdr,dUdphi; double t,t1,t2,t4,cos_t,U,uR,uA,dUdr,dUdphi;
double pi = 3.141592653589793, pow2_1by6 = 1.122462048309373; double pi = 3.141592653589793, pow2_1by6 = 1.122462048309373;
double energy_well = epsilon[type[i]][type[j]]; double energy_well = epsilon[type[i]][type[j]];
double rmin = pow2_1by6*sigma[type[i]][type[j]]; double rmin = pow2_1by6*sigma[type[i]][type[j]];
double rcut = cut[type[i]][type[j]]; double rcut = cut[type[i]][type[j]];
double zt = zeta[type[i]][type[j]]; double zt = zeta[type[i]][type[j]];
double muu = mu[type[i]][type[j]]; double muu = mu[type[i]][type[j]];
double sint = beta[type[i]][type[j]]; double sint = beta[type[i]][type[j]];
ni1[0]=a1[0][0]; ni1[0]=a1[0][0];
ni1[1]=a1[0][1]; ni1[1]=a1[0][1];
ni1[2]=a1[0][2]; ni1[2]=a1[0][2];
nj1[0]=a2[0][0]; nj1[0]=a2[0][0];
nj1[1]=a2[0][1]; nj1[1]=a2[0][1];
nj1[2]=a2[0][2]; nj1[2]=a2[0][2];
double ninj = MathExtra::dot3(ni1,nj1); double ninj = MathExtra::dot3(ni1,nj1);
double ni1rhat = MathExtra::dot3(ni1,r12hat); double ni1rhat = MathExtra::dot3(ni1,r12hat);
double nj1rhat = MathExtra::dot3(nj1,r12hat); double nj1rhat = MathExtra::dot3(nj1,r12hat);
double a = ninj + (sint-ni1rhat)*(sint+nj1rhat) - 2.0*sint*sint; double a = ninj + (sint-ni1rhat)*(sint+nj1rhat) - 2.0*sint*sint;
double phi = 1.0 + (a-1.0)*muu; double phi = 1.0 + (a-1.0)*muu;
dphi_drhat[0] = muu*( (sint-ni1rhat)*nj1[0]-ni1[0]*(sint+nj1rhat) ); dphi_drhat[0] = muu*( (sint-ni1rhat)*nj1[0]-ni1[0]*(sint+nj1rhat) );
dphi_drhat[1] = muu*( (sint-ni1rhat)*nj1[1]-ni1[1]*(sint+nj1rhat) ); dphi_drhat[1] = muu*( (sint-ni1rhat)*nj1[1]-ni1[1]*(sint+nj1rhat) );
dphi_drhat[2] = muu*( (sint-ni1rhat)*nj1[2]-ni1[2]*(sint+nj1rhat) ); dphi_drhat[2] = muu*( (sint-ni1rhat)*nj1[2]-ni1[2]*(sint+nj1rhat) );
dphi_dni1[0] = muu*(nj1[0]-r12hat[0]*(sint+nj1rhat)); dphi_dni1[0] = muu*(nj1[0]-r12hat[0]*(sint+nj1rhat));
dphi_dni1[1] = muu*(nj1[1]-r12hat[1]*(sint+nj1rhat)); dphi_dni1[1] = muu*(nj1[1]-r12hat[1]*(sint+nj1rhat));
dphi_dni1[2] = muu*(nj1[2]-r12hat[2]*(sint+nj1rhat)); dphi_dni1[2] = muu*(nj1[2]-r12hat[2]*(sint+nj1rhat));
dphi_dnj1[0] = muu*(ni1[0]+r12hat[0]*(sint-ni1rhat)); dphi_dnj1[0] = muu*(ni1[0]+r12hat[0]*(sint-ni1rhat));
dphi_dnj1[1] = muu*(ni1[1]+r12hat[1]*(sint-ni1rhat)); dphi_dnj1[1] = muu*(ni1[1]+r12hat[1]*(sint-ni1rhat));
dphi_dnj1[2] = muu*(ni1[2]+r12hat[2]*(sint-ni1rhat)); dphi_dnj1[2] = muu*(ni1[2]+r12hat[2]*(sint-ni1rhat));
if (r < rmin) if (r < rmin)
{ {
@ -480,71 +480,71 @@ double PairYLZ::ylz_analytic(const int i,const int j,double a1[3][3],
t2=t*t; t2=t*t;
t4=t2*t2; t4=t2*t2;
uR = (t4 - 2.0*t2)*energy_well; uR = (t4 - 2.0*t2)*energy_well;
U = uR +(1.0 - phi)*energy_well; U = uR +(1.0 - phi)*energy_well;
dUdr = 4.0*(t2-t4)/r*energy_well; dUdr = 4.0*(t2-t4)/r*energy_well;
dUdphi = -energy_well; dUdphi = -energy_well;
} }
else else
{ {
t = pi/2.0*(r-rmin)/(rcut-rmin); t = pi/2.0*(r-rmin)/(rcut-rmin);
cos_t = cos(t); cos_t = cos(t);
t1=cos_t; t1=cos_t;
for (int k = 1; k <= 2*zt-2; k++) {t1 *= cos_t;} // get cos()^(2zt-1) for (int k = 1; k <= 2*zt-2; k++) {t1 *= cos_t;} // get cos()^(2zt-1)
uA = -energy_well*t1*cos_t; uA = -energy_well*t1*cos_t;
U = uA*phi; U = uA*phi;
dUdr = 4.0*pi/(rcut-rmin)*(t1)*sin(t)*phi*energy_well; dUdr = 4.0*pi/(rcut-rmin)*(t1)*sin(t)*phi*energy_well;
dUdphi = uA; dUdphi = uA;
} }
dUdrhat[0] = dUdphi*dphi_drhat[0]; dUdrhat[0] = dUdphi*dphi_drhat[0];
dUdrhat[1] = dUdphi*dphi_drhat[1]; dUdrhat[1] = dUdphi*dphi_drhat[1];
dUdrhat[2] = dUdphi*dphi_drhat[2]; dUdrhat[2] = dUdphi*dphi_drhat[2];
double dUdrhatrhat = MathExtra::dot3(dUdrhat,r12hat); double dUdrhatrhat = MathExtra::dot3(dUdrhat,r12hat);
fforce[0] = dUdr*r12hat[0] + (dUdrhat[0]-dUdrhatrhat*r12hat[0])/r; fforce[0] = dUdr*r12hat[0] + (dUdrhat[0]-dUdrhatrhat*r12hat[0])/r;
fforce[1] = dUdr*r12hat[1] + (dUdrhat[1]-dUdrhatrhat*r12hat[1])/r; fforce[1] = dUdr*r12hat[1] + (dUdrhat[1]-dUdrhatrhat*r12hat[1])/r;
fforce[2] = dUdr*r12hat[2] + (dUdrhat[2]-dUdrhatrhat*r12hat[2])/r; fforce[2] = dUdr*r12hat[2] + (dUdrhat[2]-dUdrhatrhat*r12hat[2])/r;
// torque i // torque i
dUdni1[0] = dUdphi*dphi_dni1[0]; dUdni1[0] = dUdphi*dphi_dni1[0];
dUdni1[1] = dUdphi*dphi_dni1[1]; dUdni1[1] = dUdphi*dphi_dni1[1];
dUdni1[2] = dUdphi*dphi_dni1[2]; dUdni1[2] = dUdphi*dphi_dni1[2];
MathExtra::cross3(dUdni1,ni1,ttor); //minus sign is replace by swapping ni1 and dUdni1 MathExtra::cross3(dUdni1,ni1,ttor); //minus sign is replace by swapping ni1 and dUdni1
if (newton_pair || j < nlocal) { if (newton_pair || j < nlocal) {
dUdnj1[0] = dUdphi*dphi_dnj1[0]; dUdnj1[0] = dUdphi*dphi_dnj1[0];
dUdnj1[1] = dUdphi*dphi_dnj1[1]; dUdnj1[1] = dUdphi*dphi_dnj1[1];
dUdnj1[2] = dUdphi*dphi_dnj1[2]; dUdnj1[2] = dUdphi*dphi_dnj1[2];
MathExtra::cross3(dUdnj1,nj1,rtor); //minus sign is replace by swapping ni1 and dUdni1 MathExtra::cross3(dUdnj1,nj1,rtor); //minus sign is replace by swapping ni1 and dUdni1
} }
// output energy, force, torque, only for the in_two_particles input file, for checking the model implementation // output energy, force, torque, only for the in_two_particles input file, for checking the model implementation
/* /*
fprintf(screen,"energy = %f \n",U); fprintf(screen,"energy = %f \n",U);
fprintf(screen,"force_i = %f %f %f \n",fforce[0],fforce[1],fforce[2]); fprintf(screen,"force_i = %f %f %f \n",fforce[0],fforce[1],fforce[2]);
fprintf(screen,"torque_i = %f %f %f \n",ttor[0],ttor[1],ttor[2]); fprintf(screen,"torque_i = %f %f %f \n",ttor[0],ttor[1],ttor[2]);
fprintf(screen,"torque_j = %f %f %f \n",rtor[0],rtor[1],rtor[2]); fprintf(screen,"torque_j = %f %f %f \n",rtor[0],rtor[1],rtor[2]);
*/ */
return U; return U;
} }

26
src/ASPHERE/pair_ylz.h Executable file → Normal file
View File

@ -1,6 +1,6 @@
/* -*- c++ -*- ---------------------------------------------------------- /* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories https://www.lammps.org/ Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract Copyright (2003) Sandia Corporation. Under the terms of Contract
@ -40,12 +40,12 @@ class PairYLZ : public Pair {
void read_restart_settings(FILE *); void read_restart_settings(FILE *);
void write_data(FILE *); void write_data(FILE *);
void write_data_all(FILE *); void write_data_all(FILE *);
protected: protected:
double cut_global; double cut_global;
double **epsilon,**sigma,**cut,**zeta,**mu,**beta; // model parameter values for atom-type pairs double **epsilon,**sigma,**cut,**zeta,**mu,**beta; // model parameter values for atom-type pairs
class AtomVecEllipsoid *avec; class AtomVecEllipsoid *avec;
@ -63,21 +63,3 @@ class PairYLZ : public Pair {
#endif #endif
#endif #endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Incorrect args for pair coefficients
Self-explanatory. Check the input script or data file.
E: Pair ylz requires atom style ellipsoid
Self-explanatory.
*/