From b511a75575633a205776a073c544a51d674a1ad4 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Fri, 5 Nov 2010 15:57:43 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@5217 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/create_atoms.cpp | 6 +++--- src/min_fire.cpp | 30 ++++++++++++++++++++-------- src/min_quickmin.cpp | 29 +++++++++++++++++++-------- src/pair_lj_gromacs.cpp | 27 +++++++++++++------------ src/pair_lj_gromacs_coul_gromacs.cpp | 25 ++++++++++++----------- 5 files changed, 73 insertions(+), 44 deletions(-) diff --git a/src/create_atoms.cpp b/src/create_atoms.cpp index 565f50ded2..3d7ac95a73 100644 --- a/src/create_atoms.cpp +++ b/src/create_atoms.cpp @@ -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; diff --git a/src/min_fire.cpp b/src/min_fire.cpp index 67e1dd188b..3bd3efee1f 100644 --- a/src/min_fire.cpp +++ b/src/min_fire.cpp @@ -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; diff --git a/src/min_quickmin.cpp b/src/min_quickmin.cpp index 060dabe5eb..6a7112843d 100644 --- a/src/min_quickmin.cpp +++ b/src/min_quickmin.cpp @@ -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; diff --git a/src/pair_lj_gromacs.cpp b/src/pair_lj_gromacs.cpp index e7abce6af0..6fb2786b08 100644 --- a/src/pair_lj_gromacs.cpp +++ b/src/pair_lj_gromacs.cpp @@ -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; } diff --git a/src/pair_lj_gromacs_coul_gromacs.cpp b/src/pair_lj_gromacs_coul_gromacs.cpp index c7b51d45aa..f5e49a5167 100644 --- a/src/pair_lj_gromacs_coul_gromacs.cpp +++ b/src/pair_lj_gromacs_coul_gromacs.cpp @@ -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;