From e13c488f9653222d3c990aa3a0002d714d044357 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Fri, 18 Dec 2009 16:57:20 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@3573 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/compute_angle_local.cpp | 25 +++++++++++++++---------- src/compute_angle_local.h | 2 +- src/compute_bond_local.cpp | 8 ++++++-- src/fix_rigid.cpp | 21 +++++++++------------ src/fix_rigid.h | 2 +- 5 files changed, 32 insertions(+), 26 deletions(-) diff --git a/src/compute_angle_local.cpp b/src/compute_angle_local.cpp index b57cdc049b..ebb710c62b 100644 --- a/src/compute_angle_local.cpp +++ b/src/compute_angle_local.cpp @@ -47,14 +47,14 @@ ComputeAngleLocal::ComputeAngleLocal(LAMMPS *lmp, int narg, char **arg) : which = new int[nvalues]; pack_choice = new FnPtrPack[nvalues]; - dflag = eflag = 0; + tflag = eflag = 0; int i; for (int iarg = 3; iarg < narg; iarg++) { i = iarg-3; if (strcmp(arg[iarg],"theta") == 0) { - dflag = 1; + tflag = 1; which[i] = THETA; pack_choice[i] = &ComputeAngleLocal::pack_theta; } else if (strcmp(arg[iarg],"energy") == 0) { @@ -124,6 +124,7 @@ void ComputeAngleLocal::compute_local() if angle is deleted (type = 0), do not count if angle is turned off (type < 0), still count if flag is set, compute requested info about angle + if angle is turned off (type < 0), energy = 0.0 ------------------------------------------------------------------------- */ int ComputeAngleLocal::compute_angles(int flag) @@ -143,6 +144,7 @@ int ComputeAngleLocal::compute_angles(int flag) int nlocal = atom->nlocal; Angle *angle = force->angle; + double PI = 4.0*atan(1.0); int m = 0; for (atom2 = 0; atom2 < nlocal; atom2++) { @@ -153,10 +155,10 @@ int ComputeAngleLocal::compute_angles(int flag) if (atom1 < 0 || !(mask[atom1] & groupbit)) continue; atom3 = atom->map(angle_atom3[atom2][i]); if (atom3 < 0 || !(mask[atom3] & groupbit)) continue; + if (angle_type[atom2][i] == 0) continue; if (flag) { - - if (dflag) { + if (tflag) { delx1 = x[atom1][0] - x[atom2][0]; dely1 = x[atom1][1] - x[atom2][1]; delz1 = x[atom1][2] - x[atom2][2]; @@ -179,11 +181,14 @@ int ComputeAngleLocal::compute_angles(int flag) c /= r1*r2; if (c > 1.0) c = 1.0; if (c < -1.0) c = -1.0; - theta[m] = acos(c); + theta[m] = 180.0*acos(c)/PI; } - if (eflag) - energy[m] = angle->single(angle_type[atom1][i],atom1,atom2,atom3); + if (eflag) { + if (angle_type[atom2][i] > 0) + energy[m] = angle->single(angle_type[atom2][i],atom1,atom2,atom3); + else energy[m] = 0.0; + } } m++; @@ -221,7 +226,7 @@ void ComputeAngleLocal::reallocate(int n) while (nmax < n) nmax += DELTA; - if (dflag) { + if (tflag) { memory->sfree(theta); theta = (double *) memory->smalloc(nmax*sizeof(double), "bond/local:theta"); @@ -233,7 +238,7 @@ void ComputeAngleLocal::reallocate(int n) } if (nvalues == 1) { - if (dflag) vector_local = theta; + if (tflag) vector_local = theta; if (eflag) vector_local = energy; } else { memory->destroy_2d_double_array(array); @@ -250,7 +255,7 @@ void ComputeAngleLocal::reallocate(int n) double ComputeAngleLocal::memory_usage() { double bytes = 0.0; - if (dflag) bytes += nmax * sizeof(double); + if (tflag) bytes += nmax * sizeof(double); if (eflag) bytes += nmax * sizeof(double); if (nvalues > 1) bytes += nmax*nvalues * sizeof(double); return bytes; diff --git a/src/compute_angle_local.h b/src/compute_angle_local.h index bc1bf849f4..c82a4f5a72 100644 --- a/src/compute_angle_local.h +++ b/src/compute_angle_local.h @@ -27,7 +27,7 @@ class ComputeAngleLocal : public Compute { double memory_usage(); private: - int nvalues,dflag,eflag; + int nvalues,tflag,eflag; int *which; int ncount; diff --git a/src/compute_bond_local.cpp b/src/compute_bond_local.cpp index a0179d29c2..429aa63228 100644 --- a/src/compute_bond_local.cpp +++ b/src/compute_bond_local.cpp @@ -124,6 +124,7 @@ void ComputeBondLocal::compute_local() if bond is deleted (type = 0), do not count if bond is turned off (type < 0), still count if flag is set, compute requested info about bond + if bond is turned off (type < 0), energy = 0.0 ------------------------------------------------------------------------- */ int ComputeBondLocal::compute_bonds(int flag) @@ -158,8 +159,11 @@ int ComputeBondLocal::compute_bonds(int flag) domain->minimum_image(delx,dely,delz); rsq = delx*delx + dely*dely + delz*delz; if (dflag) distance[m] = sqrt(rsq); - if (eflag) - energy[m] = bond->single(bond_type[atom1][i],rsq,atom1,atom2); + if (eflag) { + if (bond_type[atom1][i] > 0) + energy[m] = bond->single(bond_type[atom1][i],rsq,atom1,atom2); + else energy[m] = 0.0; + } } m++; diff --git a/src/fix_rigid.cpp b/src/fix_rigid.cpp index 15d41bba11..7af5e1d268 100644 --- a/src/fix_rigid.cpp +++ b/src/fix_rigid.cpp @@ -203,10 +203,11 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) : // initialize force/torque flags to default = 1.0 - vector_flag = 1; - size_vector = 12*nbody; + array_flag = 1; + size_array_rows = nbody; + size_array_cols = 12; global_freq = 1; - extvector = 0; + extarray = 0; for (i = 0; i < nbody; i++) { fflag[i][0] = fflag[i][1] = fflag[i][2] = 1.0; @@ -2139,14 +2140,10 @@ void FixRigid::reset_dt() xcm = 1,2,3; vcm = 4,5,6; fcm = 7,8,9; torque = 10,11,12 ------------------------------------------------------------------------- */ -double FixRigid::compute_vector(int n) +double FixRigid::compute_array(int i, int j) { - int ibody = n/12; - int iattribute = (n % 12) / 3; - int index = n % 3; - - if (iattribute == 0) return xcm[ibody][index]; - else if (iattribute == 1) return vcm[ibody][index]; - else if (iattribute == 2) return fcm[ibody][index]; - else return torque[ibody][index]; + if (j <= 3) return xcm[i][j]; + if (j <= 6) return vcm[i][j-3]; + if (j <= 9) return fcm[i][j-6]; + return torque[i][j-9]; } diff --git a/src/fix_rigid.h b/src/fix_rigid.h index 5617ce774b..209b5df018 100644 --- a/src/fix_rigid.h +++ b/src/fix_rigid.h @@ -41,7 +41,7 @@ class FixRigid : public Fix { int dof(int); void deform(int); void reset_dt(); - double compute_vector(int); + double compute_array(int, int); private: double dtv,dtf,dtq;