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

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

This commit is contained in:
sjplimp 2010-11-05 15:57:43 +00:00
parent 236cf58328
commit b511a75575
5 changed files with 73 additions and 44 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
src/create_atoms.cpp
View File

@ -309,9 +309,9 @@ void CreateAtoms::add_random()
if (domain->triclinic) {
domain->x2lamda(xone,lamda);
coord = lamda;
if (coord[0] >= boxlo[0] && coord[0] < boxhi[0] &&
coord[1] >= boxlo[1] && coord[1] < boxhi[1] &&
coord[2] >= boxlo[2] && coord[2] < boxhi[2]) valid = 0;
if (coord[0] < boxlo[0] || coord[0] >= boxhi[0] ||
coord[1] < boxlo[1] || coord[1] >= boxhi[1] ||
coord[2] < boxlo[2] || coord[2] >= boxhi[2]) valid = 0;
} else coord = xone;
if (valid) break;

30
src/min_fire.cpp
View File

@ -149,6 +149,7 @@ int MinFire::iterate(int maxiter)
// limit timestep so no particle moves further than dmax
double *rmass = atom->rmass;
double *mass = atom->mass;
int *type = atom->type;
@ -165,14 +166,27 @@ int MinFire::iterate(int maxiter)
double **x = atom->x;
for (int i = 0; i < nlocal; i++) {
dtfm = dtv / mass[type[i]];
x[i][0] += dtv * v[i][0];
x[i][1] += dtv * v[i][1];
x[i][2] += dtv * v[i][2];
v[i][0] += dtfm * f[i][0];
v[i][1] += dtfm * f[i][1];
v[i][2] += dtfm * f[i][2];
if (rmass) {
for (int i = 0; i < nlocal; i++) {
dtfm = dtv / rmass[i];
x[i][0] += dtv * v[i][0];
x[i][1] += dtv * v[i][1];
x[i][2] += dtv * v[i][2];
v[i][0] += dtfm * f[i][0];
v[i][1] += dtfm * f[i][1];
v[i][2] += dtfm * f[i][2];
}
} else {
for (int i = 0; i < nlocal; i++) {
dtfm = dtv / mass[type[i]];
x[i][0] += dtv * v[i][0];
x[i][1] += dtv * v[i][1];
x[i][2] += dtv * v[i][2];
v[i][0] += dtfm * f[i][0];
v[i][1] += dtfm * f[i][1];
v[i][2] += dtfm * f[i][2];
}
}
eprevious = ecurrent;

29
src/min_quickmin.cpp
View File

@ -119,6 +119,7 @@ int MinQuickMin::iterate(int maxiter)
// limit timestep so no particle moves further than dmax
double *rmass = atom->rmass;
double *mass = atom->mass;
int *type = atom->type;
@ -135,14 +136,26 @@ int MinQuickMin::iterate(int maxiter)
double **x = atom->x;
for (int i = 0; i < nlocal; i++) {
dtfm = dtv / mass[type[i]];
x[i][0] += dtv * v[i][0];
x[i][1] += dtv * v[i][1];
x[i][2] += dtv * v[i][2];
v[i][0] += dtfm * f[i][0];
v[i][1] += dtfm * f[i][1];
v[i][2] += dtfm * f[i][2];
if (rmass) {
for (int i = 0; i < nlocal; i++) {
dtfm = dtv / rmass[i];
x[i][0] += dtv * v[i][0];
x[i][1] += dtv * v[i][1];
x[i][2] += dtv * v[i][2];
v[i][0] += dtfm * f[i][0];
v[i][1] += dtfm * f[i][1];
v[i][2] += dtfm * f[i][2];
}
} else {
for (int i = 0; i < nlocal; i++) {
dtfm = dtv / mass[type[i]];
x[i][0] += dtv * v[i][0];
x[i][1] += dtv * v[i][1];
x[i][2] += dtv * v[i][2];
v[i][0] += dtfm * f[i][0];
v[i][1] += dtfm * f[i][1];
v[i][2] += dtfm * f[i][2];
}
}
eprevious = ecurrent;

27
src/pair_lj_gromacs.cpp
View File

@ -138,10 +138,10 @@ void PairLJGromacs::compute(int eflag, int vflag)
}
if (eflag) {
evdwl = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) +
ljsw5[itype][jtype];
evdwl = r6inv * (lj3[itype][jtype]*r6inv - lj4[itype][jtype]);
if (rsq > cut_inner_sq[itype][jtype]) {
eswitch = t*t*t*(ljsw3[itype][jtype] + ljsw4[itype][jtype]*t);
eswitch = t*t*t*(ljsw3[itype][jtype] + ljsw4[itype][jtype]*t) +
ljsw5[itype][jtype];
evdwl += eswitch;
}
evdwl *= factor_lj;
@ -283,17 +283,18 @@ double PairLJGromacs::init_one(int i, int j)
double t2inv = 1.0/(t*t);
double t3inv = t2inv/t;
double t3 = 1.0/t3inv;
double a6 = ( 7.0*cut_inner[i][j] - 10.0*cut[i][j])*r8inv*t2inv;
double b6 = ( 9.0*cut[i][j] - 7.0*cut_inner[i][j])*r8inv*t3inv;
double a6 = (7.0*cut_inner[i][j] - 10.0*cut[i][j])*r8inv*t2inv;
double b6 = (9.0*cut[i][j] - 7.0*cut_inner[i][j])*r8inv*t3inv;
double a12 = (13.0*cut_inner[i][j] - 16.0*cut[i][j])*r6inv*r8inv*t2inv;
double b12 = (15.0*cut[i][j] - 13.0*cut_inner[i][j])*r6inv*r8inv*t3inv;
double c6 = r6inv - t3*(a6/3.0 + b6*t/4.0);
double c12 = r6inv*r6inv - t3*(a12/3.0 + b12*t/4.0);
double c6 = r6inv - t3*(6.0*a6/3.0 + 6.0*b6*t/4.0);
double c12 = r6inv*r6inv - t3*(12.0*a12/3.0 + 12.0*b12*t/4.0);
ljsw1[i][j] = lj1[i][j]*a12 - lj2[i][j]*a6;
ljsw2[i][j] = lj1[i][j]*b12 - lj2[i][j]*b6;
ljsw3[i][j] =-lj3[i][j]*a12/3.0 + lj4[i][j]*a6/3.0;
ljsw4[i][j] =-lj3[i][j]*b12/4.0 + lj4[i][j]*b6/4.0;
ljsw5[i][j] =-lj3[i][j]*c12 + lj4[i][j]*c6;
ljsw3[i][j] = -lj3[i][j]*12.0*a12/3.0 + lj4[i][j]*6.0*a6/3.0;
ljsw4[i][j] = -lj3[i][j]*12.0*b12/4.0 + lj4[i][j]*6.0*b6/4.0;
ljsw5[i][j] = -lj3[i][j]*c12 + lj4[i][j]*c6;
cut_inner[j][i] = cut_inner[i][j];
cut_inner_sq[j][i] = cut_inner_sq[i][j];
@ -413,10 +414,10 @@ double PairLJGromacs::single(int i, int j, int itype, int jtype,
}
fforce = factor_lj*forcelj*r2inv;
philj = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) +
ljsw5[itype][jtype];
philj = r6inv * (lj3[itype][jtype]*r6inv - lj4[itype][jtype]);
if (rsq > cut_inner_sq[itype][jtype]) {
phiswitch = t*t*t*(ljsw3[itype][jtype] + ljsw4[itype][jtype]*t);
phiswitch = t*t*t*(ljsw3[itype][jtype] + ljsw4[itype][jtype]*t) +
ljsw5[itype][jtype];
philj += phiswitch;
}

25
src/pair_lj_gromacs_coul_gromacs.cpp
View File

@ -166,11 +166,11 @@ void PairLJGromacsCoulGromacs::compute(int eflag, int vflag)
ecoul *= factor_coul;
} else ecoul = 0.0;
if (rsq < cut_ljsq) {
evdwl = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) +
ljsw5[itype][jtype];
evdwl = r6inv * (lj3[itype][jtype]*r6inv - lj4[itype][jtype]);
if (rsq > cut_lj_innersq) {
eswitch = tlj*tlj*tlj * (ljsw3[itype][jtype] +
ljsw4[itype][jtype]*tlj);
eswitch = tlj*tlj*tlj *
(ljsw3[itype][jtype] + ljsw4[itype][jtype]*tlj) +
ljsw5[itype][jtype];
evdwl += eswitch;
}
evdwl *= factor_lj;
@ -316,12 +316,13 @@ double PairLJGromacsCoulGromacs::init_one(int i, int j)
double b6 = (9.0*cut_lj - 7.0*cut_lj_inner)*r8inv*t3inv;
double a12 = (13.0*cut_lj_inner - 16.0*cut_lj)*r6inv*r8inv*t2inv;
double b12 = (15.0*cut_lj - 13.0*cut_lj_inner)*r6inv*r8inv*t3inv;
double c6 = r6inv - t3*(a6/3.0 + b6*t/4.0);
double c12 = r6inv*r6inv - t3*(a12/3.0 + b12*t/4.0);
double c6 = r6inv - t3*(6.0*a6/3.0 + 6.0*b6*t/4.0);
double c12 = r6inv*r6inv - t3*(12.0*a12/3.0 + 12.0*b12*t/4.0);
ljsw1[i][j] = lj1[i][j]*a12 - lj2[i][j]*a6;
ljsw2[i][j] = lj1[i][j]*b12 - lj2[i][j]*b6;
ljsw3[i][j] = -lj3[i][j]*a12/3.0 + lj4[i][j]*a6/3.0;
ljsw4[i][j] = -lj3[i][j]*b12/4.0 + lj4[i][j]*b6/4.0;
ljsw3[i][j] = -lj3[i][j]*12.0*a12/3.0 + lj4[i][j]*6.0*a6/3.0;
ljsw4[i][j] = -lj3[i][j]*12.0*b12/4.0 + lj4[i][j]*6.0*b6/4.0;
ljsw5[i][j] = -lj3[i][j]*c12 + lj4[i][j]*c6;
double r3inv = 1.0/pow(cut_coul,3.0);
@ -477,11 +478,11 @@ double PairLJGromacsCoulGromacs::single(int i, int j, int itype, int jtype,
}
if (rsq < cut_ljsq) {
philj = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) +
ljsw5[itype][jtype];
philj = r6inv * (lj3[itype][jtype]*r6inv - lj4[itype][jtype]);
if (rsq > cut_lj_innersq) {
phiswitch = tlj*tlj*tlj * (ljsw3[itype][jtype] +
ljsw4[itype][jtype]*tlj);
phiswitch = tlj*tlj*tlj *
(ljsw3[itype][jtype] + ljsw4[itype][jtype]*tlj) +
ljsw5[itype][jtype];
philj += phiswitch;
}
eng += factor_lj*philj;