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

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@9410 f3b2605a-c512-4ea7-a41b-209d697bcdaa

Browse Source
This commit is contained in:
sjplimp
2013-02-08 17:37:46 +00:00
parent d41265166e
commit 3adb315c9e
14 changed files with 37 additions and 15 deletions
Download Patch File Download Diff File Expand all files Collapse all files

5
src/MOLECULE/bond_fene.cpp
View File

@ -231,7 +231,8 @@ void BondFENE::read_restart(FILE *fp)
/* ---------------------------------------------------------------------- */
double BondFENE::single(int type, double rsq, int i, int j)
double BondFENE::single(int type, double rsq, int i, int j,
double &fforce)
{
double r0sq = r0[type] * r0[type];
double rlogarg = 1.0 - rsq/r0sq;
@ -250,11 +251,13 @@ double BondFENE::single(int type, double rsq, int i, int j)
}
double eng = -0.5 * k[type]*r0sq*log(rlogarg);
fforce = -k[type]/rlogarg;
if (rsq < TWO_1_3*sigma[type]*sigma[type]) {
double sr2,sr6;
sr2 = sigma[type]*sigma[type]/rsq;
sr6 = sr2*sr2*sr2;
eng += 4.0*epsilon[type]*sr6*(sr6-1.0) + epsilon[type];
fforce += 48.0*epsilon[type]*sr6*(sr6-0.5)/rsq;
}
return eng;

2
src/MOLECULE/bond_fene.h
View File

@ -35,7 +35,7 @@ class BondFENE : public Bond {
double equilibrium_distance(int);
void write_restart(FILE *);
void read_restart(FILE *);
double single(int, double, int, int);
double single(int, double, int, int, double &);
protected:
double TWO_1_3;

5
src/MOLECULE/bond_fene_expand.cpp
View File

@ -242,7 +242,8 @@ void BondFENEExpand::read_restart(FILE *fp)
/* ---------------------------------------------------------------------- */
double BondFENEExpand::single(int type, double rsq, int i, int j)
double BondFENEExpand::single(int type, double rsq, int i, int j,
double &fforce)
{
double r = sqrt(rsq);
double rshift = r - shift[type];
@ -264,11 +265,13 @@ double BondFENEExpand::single(int type, double rsq, int i, int j)
}
double eng = -0.5 * k[type]*r0sq*log(rlogarg);
fforce = -k[type]*rshift/rlogarg/r;
if (rshiftsq < TWO_1_3*sigma[type]*sigma[type]) {
double sr2,sr6;
sr2 = sigma[type]*sigma[type]/rshiftsq;
sr6 = sr2*sr2*sr2;
eng += 4.0*epsilon[type]*sr6*(sr6-1.0) + epsilon[type];
fforce += 48.0*epsilon[type]*sr6*(sr6-0.5)/rshift/r;
}
return eng;

2
src/MOLECULE/bond_fene_expand.h
View File

@ -35,7 +35,7 @@ class BondFENEExpand : public Bond {
double equilibrium_distance(int);
void write_restart(FILE *);
void read_restart(FILE *);
double single(int, double, int, int);
double single(int, double, int, int, double &);
protected:
double TWO_1_3;

5
src/MOLECULE/bond_harmonic.cpp
View File

@ -176,10 +176,13 @@ void BondHarmonic::read_restart(FILE *fp)
/* ---------------------------------------------------------------------- */
double BondHarmonic::single(int type, double rsq, int i, int j)
double BondHarmonic::single(int type, double rsq, int i, int j,
double &fforce)
{
double r = sqrt(rsq);
double dr = r - r0[type];
double rk = k[type] * dr;
fforce = 0;
if (r > 0.0) fforce = -2.0*rk/r;
return rk*dr;
}

2
src/MOLECULE/bond_harmonic.h
View File

@ -34,7 +34,7 @@ class BondHarmonic : public Bond {
double equilibrium_distance(int);
void write_restart(FILE *);
void read_restart(FILE *);
double single(int, double, int, int);
double single(int, double, int, int, double &);
protected:
double *k,*r0;

5
src/MOLECULE/bond_morse.cpp
View File

@ -186,10 +186,13 @@ void BondMorse::read_restart(FILE *fp)
/* ---------------------------------------------------------------------- */
double BondMorse::single(int type, double rsq, int i, int j)
double BondMorse::single(int type, double rsq, int i, int j,
double &fforce)
{
double r = sqrt(rsq);
double dr = r - r0[type];
double ralpha = exp(-alpha[type]*dr);
fforce = 0;
if (r > 0.0) fforce = -2.0*d0[type]*alpha[type]*(1-ralpha)*ralpha/r;
return d0[type]*(1-ralpha)*(1-ralpha);
}

2
src/MOLECULE/bond_morse.h
View File

@ -34,7 +34,7 @@ class BondMorse : public Bond {
double equilibrium_distance(int);
void write_restart(FILE *);
void read_restart(FILE *);
double single(int, double, int, int);
double single(int, double, int, int, double &);
protected:
double *d0,*alpha,*r0;

5
src/MOLECULE/bond_nonlinear.cpp
View File

@ -181,12 +181,15 @@ void BondNonlinear::read_restart(FILE *fp)
/* ---------------------------------------------------------------------- */
double BondNonlinear::single(int type, double rsq, int i, int j)
double BondNonlinear::single(int type, double rsq, int i, int j,
double &fforce)
{
double r = sqrt(rsq);
double dr = r - r0[type];
double drsq = dr*dr;
double lamdasq = lamda[type]*lamda[type];
double denom = lamdasq - drsq;
double denomsq = denom*denom;
fforce = -epsilon[type]/r * 2.0*dr*lamdasq/denomsq;
return epsilon[type] * drsq / denom;
}

2
src/MOLECULE/bond_nonlinear.h
View File

@ -34,7 +34,7 @@ class BondNonlinear : public Bond {
double equilibrium_distance(int);
void write_restart(FILE *);
void read_restart(FILE *);
double single(int, double, int, int);
double single(int, double, int, int, double &);
protected:
double *epsilon,*r0,*lamda;

6
src/MOLECULE/bond_quartic.cpp
View File

@ -296,7 +296,8 @@ void BondQuartic::read_restart(FILE *fp)
/* ---------------------------------------------------------------------- */
double BondQuartic::single(int type, double rsq, int i, int j)
double BondQuartic::single(int type, double rsq, int i, int j,
double &fforce)
{
double r,dr,r2,ra,rb,sr2,sr6;
@ -325,10 +326,13 @@ double BondQuartic::single(int type, double rsq, int i, int j)
rb = dr - b2[type];
eng += k[type]*r2*ra*rb + u0[type];
fforce = -k[type]/r * (r2*(ra+rb) + 2.0*dr*ra*rb);
if (rsq < TWO_1_3) {
sr2 = 1.0/rsq;
sr6 = sr2*sr2*sr2;
eng += 4.0*sr6*(sr6-1.0) + 1.0;
fforce += 48.0*sr6*(sr6-0.5)/rsq;
}
return eng;

2
src/MOLECULE/bond_quartic.h
View File

@ -35,7 +35,7 @@ class BondQuartic : public Bond {
double equilibrium_distance(int);
void write_restart(FILE *);
void read_restart(FILE *);
double single(int, double, int, int);
double single(int, double, int, int, double &);
protected:
double TWO_1_3;

7
src/MOLECULE/bond_table.cpp
View File

@ -241,11 +241,14 @@ void BondTable::read_restart(FILE *fp)
/* ---------------------------------------------------------------------- */
double BondTable::single(int type, double rsq, int i, int j)
double BondTable::single(int type, double rsq, int i, int j,
double &fforce)
{
double r = sqrt(rsq);
double u;
u_lookup(type,r,u);
double mdu;
uf_lookup(type,r,u,mdu);
fforce = mdu/r;
return u;
}

2
src/MOLECULE/bond_table.h
View File

@ -35,7 +35,7 @@ class BondTable : public Bond {
double equilibrium_distance(int);
void write_restart(FILE *);
void read_restart(FILE *);
double single(int, double, int, int);
double single(int, double, int, int, double &);
protected:
int tabstyle,tablength;