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

enable and apply clang-format

Browse Source
This commit is contained in:
Axel Kohlmeyer
2024-08-06 15:34:16 -04:00
parent ce3b3b2950
commit 2344b7bcc3
1 changed files with 283 additions and 271 deletions
Download Patch File Download Diff File Expand all files Collapse all files

554
src/EXTRA-COMPUTE/compute_basal_atom.cpp
View File

@ -1,4 +1,3 @@
// clang-format off
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories https://www.lammps.org/, Sandia National Laboratories
@ -37,10 +36,9 @@ using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
ComputeBasalAtom::ComputeBasalAtom(LAMMPS *lmp, int narg, char **arg) : ComputeBasalAtom::ComputeBasalAtom(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg)
Compute(lmp, narg, arg)
{ {
if (narg != 3) error->all(FLERR,"Illegal compute basal/atom command"); if (narg != 3) error->all(FLERR, "Illegal compute basal/atom command");
peratom_flag = 1; peratom_flag = 1;
size_peratom_cols = 3; size_peratom_cols = 3;
@ -74,7 +72,7 @@ void ComputeBasalAtom::init()
neighbor->add_request(this, NeighConst::REQ_FULL | NeighConst::REQ_OCCASIONAL); neighbor->add_request(this, NeighConst::REQ_FULL | NeighConst::REQ_OCCASIONAL);
if ((modify->get_compute_by_style("basal/atom").size() > 1) && (comm->me == 0)) if ((modify->get_compute_by_style("basal/atom").size() > 1) && (comm->me == 0))
error->warning(FLERR,"More than one compute basal/atom"); error->warning(FLERR, "More than one compute basal/atom");
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
@ -88,16 +86,16 @@ void ComputeBasalAtom::init_list(int /*id*/, NeighList *ptr)
void ComputeBasalAtom::compute_peratom() void ComputeBasalAtom::compute_peratom()
{ {
int i,j,ii,jj,k,n,inum,jnum; int i, j, ii, jj, k, n, inum, jnum;
double xtmp,ytmp,ztmp,delx,dely,delz,rsq,var5,var6,var7; double xtmp, ytmp, ztmp, delx, dely, delz, rsq, var5, var6, var7;
int *ilist,*jlist,*numneigh,**firstneigh; int *ilist, *jlist, *numneigh, **firstneigh;
int chi[8]; int chi[8];
int value; int value;
int count; int count;
int k2[3]; int k2[3];
int j1[3]; int j1[3];
double x4[3],y4[3],z4[3],x5[3],y5[3],z5[3],x6[3],y6[3],z6[3]; double x4[3], y4[3], z4[3], x5[3], y5[3], z5[3], x6[3], y6[3], z6[3];
double x7[3],y7[3],z7[3]; double x7[3], y7[3], z7[3];
invoked_peratom = update->ntimestep; invoked_peratom = update->ntimestep;
@ -106,7 +104,7 @@ void ComputeBasalAtom::compute_peratom()
if (atom->nmax > nmax) { if (atom->nmax > nmax) {
memory->destroy(BPV); memory->destroy(BPV);
nmax = atom->nmax; nmax = atom->nmax;
memory->create(BPV,nmax,3,"basal/atom:basal"); memory->create(BPV, nmax, 3, "basal/atom:basal");
array_atom = BPV; array_atom = BPV;
} }
@ -143,10 +141,10 @@ void ComputeBasalAtom::compute_peratom()
memory->destroy(nearest_n0); memory->destroy(nearest_n0);
memory->destroy(nearest_n1); memory->destroy(nearest_n1);
maxneigh = jnum; maxneigh = jnum;
memory->create(distsq,maxneigh,"compute/basal/atom:distsq"); memory->create(distsq, maxneigh, "compute/basal/atom:distsq");
memory->create(nearest,maxneigh,"compute/basal/atom:nearest"); memory->create(nearest, maxneigh, "compute/basal/atom:nearest");
memory->create(nearest_n0,maxneigh,"compute/basal/atom:nearest_n0"); memory->create(nearest_n0, maxneigh, "compute/basal/atom:nearest_n0");
memory->create(nearest_n1,maxneigh,"compute/basal/atom:nearest_n1"); memory->create(nearest_n1, maxneigh, "compute/basal/atom:nearest_n1");
} }
// neighbor selection is identical to ackland/atom algorithm // neighbor selection is identical to ackland/atom algorithm
@ -162,7 +160,7 @@ void ComputeBasalAtom::compute_peratom()
delx = xtmp - x[j][0]; delx = xtmp - x[j][0];
dely = ytmp - x[j][1]; dely = ytmp - x[j][1];
delz = ztmp - x[j][2]; delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz; rsq = delx * delx + dely * dely + delz * delz;
if (rsq < cutsq) { if (rsq < cutsq) {
distsq[n] = rsq; distsq[n] = rsq;
nearest[n++] = j; nearest[n++] = j;
@ -171,7 +169,7 @@ void ComputeBasalAtom::compute_peratom()
// Select 6 nearest neighbors // Select 6 nearest neighbors
select2(6,n,distsq,nearest); select2(6, n, distsq, nearest);
// Mean squared separation // Mean squared separation
@ -182,248 +180,258 @@ void ComputeBasalAtom::compute_peratom()
// n0 near neighbors with: distsq<1.45*r0_sq // n0 near neighbors with: distsq<1.45*r0_sq
// n1 near neighbors with: distsq<1.55*r0_sq // n1 near neighbors with: distsq<1.55*r0_sq
double n0_dist_sq = 1.45*r0_sq, double n0_dist_sq = 1.45 * r0_sq, n1_dist_sq = 1.55 * r0_sq;
n1_dist_sq = 1.55*r0_sq;
int n0 = 0, n1 = 0; int n0 = 0, n1 = 0;
for (j = 0; j < n; j++) { for (j = 0; j < n; j++) {
if (distsq[j] < n1_dist_sq) { if (distsq[j] < n1_dist_sq) {
nearest_n1[n1++] = nearest[j]; nearest_n1[n1++] = nearest[j];
if (distsq[j] < n0_dist_sq) { if (distsq[j] < n0_dist_sq) nearest_n0[n0++] = nearest[j];
nearest_n0[n0++] = nearest[j]; }
}
}
} }
// Evaluate all angles <(r_ij,rik) forall n0 particles with: distsq<1.45*r0_sq // Evaluate all angles <(r_ij,rik) forall n0 particles with: distsq<1.45*r0_sq
double bond_angle; double bond_angle;
double norm_j, norm_k; double norm_j, norm_k;
chi[0] = chi[1] = chi[2] = chi[3] = chi[4] = chi[5] = chi[6] = chi[7] = 0; chi[0] = chi[1] = chi[2] = chi[3] = chi[4] = chi[5] = chi[6] = chi[7] = 0;
double x_ij, y_ij, z_ij, x_ik, y_ik, z_ik, xmean5, ymean5, zmean5, double x_ij, y_ij, z_ij, x_ik, y_ik, z_ik, xmean5, ymean5, zmean5, xmean6, ymean6, zmean6,
xmean6, ymean6, zmean6, xmean7, ymean7, zmean7; xmean7, ymean7, zmean7;
auto x3 = new double[n0]; auto x3 = new double[n0];
auto y3 = new double[n0]; auto y3 = new double[n0];
auto z3 = new double[n0]; auto z3 = new double[n0];
for (j = 0; j < n0; j++) { for (j = 0; j < n0; j++) {
x_ij = x[i][0]-x[nearest_n0[j]][0]; x_ij = x[i][0] - x[nearest_n0[j]][0];
y_ij = x[i][1]-x[nearest_n0[j]][1]; y_ij = x[i][1] - x[nearest_n0[j]][1];
z_ij = x[i][2]-x[nearest_n0[j]][2]; z_ij = x[i][2] - x[nearest_n0[j]][2];
norm_j = sqrt (x_ij*x_ij + y_ij*y_ij + z_ij*z_ij); norm_j = sqrt(x_ij * x_ij + y_ij * y_ij + z_ij * z_ij);
if (norm_j <= 0.) {continue;} if (norm_j <= 0.0) continue;
for (k = j+1; k < n0; k++) { for (k = j + 1; k < n0; k++) {
x_ik = x[i][0]-x[nearest_n0[k]][0]; x_ik = x[i][0] - x[nearest_n0[k]][0];
y_ik = x[i][1]-x[nearest_n0[k]][1]; y_ik = x[i][1] - x[nearest_n0[k]][1];
z_ik = x[i][2]-x[nearest_n0[k]][2]; z_ik = x[i][2] - x[nearest_n0[k]][2];
norm_k = sqrt (x_ik*x_ik + y_ik*y_ik + z_ik*z_ik); norm_k = sqrt(x_ik * x_ik + y_ik * y_ik + z_ik * z_ik);
if (norm_k <= 0.) {continue;} if (norm_k <= 0.) { continue; }
bond_angle = (x_ij*x_ik + y_ij*y_ik + z_ij*z_ik) / (norm_j*norm_k); bond_angle = (x_ij * x_ik + y_ij * y_ik + z_ij * z_ik) / (norm_j * norm_k);
//find all bond angles that are about 180 degrees //find all bond angles that are about 180 degrees
if (-1. <= bond_angle && bond_angle < -0.945) { if (-1. <= bond_angle && bond_angle < -0.945) {
x3[chi[0]] = x_ik - x_ij; x3[chi[0]] = x_ik - x_ij;
y3[chi[0]] = y_ik - y_ij; y3[chi[0]] = y_ik - y_ij;
z3[chi[0]] = z_ik - z_ij; z3[chi[0]] = z_ik - z_ij;
chi[0]++; chi[0]++;
} }
} }
} }
// for atoms that have 2 or 3 ~180 bond angles: // for atoms that have 2 or 3 ~180 bond angles:
if (2 == chi[0] || 3 == chi[0]) { if (2 == chi[0] || 3 == chi[0]) {
count = 0; count = 0;
if (chi[0] == 2) { if (chi[0] == 2) {
k2[0] = 0; k2[0] = 0;
j1[0] = 1; j1[0] = 1;
} } else {
else {
k2[0] = 0;
k2[1] = 0;
k2[2] = 1;
j1[0]=1;
j1[1]=2;
j1[2]=2;
}
xmean5 = ymean5 = zmean5 = xmean6 = ymean6 = zmean6 = xmean7 = ymean7 = zmean7 = 0.0;
for (j = 0; j < chi[0]; j++) {
for (k = j+1; k < chi[0]; k++) {
//get cross products
x4[count] = y3[j1[count]]*z3[k2[count]]-y3[k2[count]]*z3[j1[count]];
y4[count] = z3[j1[count]]*x3[k2[count]]-z3[k2[count]]*x3[j1[count]];
z4[count] = x3[j1[count]]*y3[k2[count]]-x3[k2[count]]*y3[j1[count]];
//get all sign combinations of cross products
x5[count] = x4[count]*copysign(1.0,x4[count]);
y5[count] = y4[count]*copysign(1.0,x4[count]);
z5[count] = z4[count]*copysign(1.0,x4[count]);
x6[count] = x4[count]*copysign(1.0,y4[count]);
y6[count] = y4[count]*copysign(1.0,y4[count]);
z6[count] = z4[count]*copysign(1.0,y4[count]);
x7[count] = x4[count]*copysign(1.0,z4[count]);
y7[count] = y4[count]*copysign(1.0,z4[count]);
z7[count] = z4[count]*copysign(1.0,z4[count]);
//get average cross products
xmean5 += x5[count];
ymean5 += y5[count];
zmean5 += z5[count];
xmean6 += x6[count];
ymean6 += y6[count];
zmean6 += z6[count];
xmean7 += x7[count];
ymean7 += y7[count];
zmean6 += z7[count];
count++;
}
}
if (count > 0) {
xmean5 /= count;
xmean6 /= count;
xmean7 /= count;
ymean5 /= count;
ymean6 /= count;
ymean7 /= count;
zmean5 /= count;
zmean6 /= count;
zmean7 /= count;
}
var5 = var6 = var7 = 0.0;
//find standard deviations
for (j=0;j<count;j++) {
var5 = var5 + x5[j]*x5[j]-2*x5[j]*xmean5+xmean5*xmean5+y5[j]*y5[j]-2*y5[j]*ymean5+ymean5*ymean5+z5[j]*z5[j]-2*z5[j]*zmean5+zmean5*zmean5;
var6 = var6 + x6[j]*x6[j]-2*x6[j]*xmean6+xmean6*xmean6+y6[j]*y6[j]-2*y6[j]*ymean6+ymean6*ymean6+z6[j]*z6[j]-2*z6[j]*zmean6+zmean6*zmean6;
var7 = var7 + x7[j]*x7[j]-2*x7[j]*xmean7+xmean7*xmean7+y7[j]*y7[j]-2*y7[j]*ymean7+ymean7*ymean7+z7[j]*z7[j]-2*z7[j]*zmean7+zmean7*zmean7;
}
//select sign combination with minimum standard deviation
if (var5 < var6) {
if (var5 < var7) { value = 0;}
else {value = 2;}
}
else if (var6 < var7) {value = 1;}
else {value = 2;}
//BPV is average of cross products of all neighbor vectors which are part of 180 degree angles
BPV[i][0] = 0;
BPV[i][1] = 0;
BPV[i][2] = 0;
for (k=0;k<count;k++) {
if (value == 0) {
BPV[i][0] = BPV[i][0]+x5[k];
BPV[i][1] = BPV[i][1]+y5[k];
BPV[i][2] = BPV[i][2]+z5[k];
}
else if (value == 1) {
BPV[i][0] = BPV[i][0]+x6[k];
BPV[i][1] = BPV[i][1]+y6[k];
BPV[i][2] = BPV[i][2]+z6[k];
}
else {
BPV[i][0] = BPV[i][0]+x7[k];
BPV[i][1] = BPV[i][1]+y7[k];
BPV[i][2] = BPV[i][2]+z7[k];
}
}
}
//for atoms with more than three 180 degree bond angles:
else if (chi[0] > 3) {
double x44[3], y44[3], z44[3], S0;
int l, m;
count = 0;
S0 = 100000;
k2[0] = 0; k2[0] = 0;
k2[1] = 0; k2[1] = 0;
k2[2] = 1; k2[2] = 1;
j1[0]=1; j1[0] = 1;
j1[1]=2; j1[1] = 2;
j1[2]=2; j1[2] = 2;
//algorithm is as above, but now all combinations of three 180 degree angles are compared, and the combination with minimum standard deviation is chosen }
for (j=0; j<chi[0]; j++) { xmean5 = ymean5 = zmean5 = xmean6 = ymean6 = zmean6 = xmean7 = ymean7 = zmean7 = 0.0;
for (k=j+1; k<chi[0]; k++) { for (j = 0; j < chi[0]; j++) {
for (l=k+1; l<chi[0]; l++) { for (k = j + 1; k < chi[0]; k++) {
if (k >= chi[0] || l >= chi[0]) continue; //get cross products
//get unique combination of three neighbor vectors x4[count] = y3[j1[count]] * z3[k2[count]] - y3[k2[count]] * z3[j1[count]];
x4[0] = x3[j]; y4[count] = z3[j1[count]] * x3[k2[count]] - z3[k2[count]] * x3[j1[count]];
x4[1] = x3[k]; z4[count] = x3[j1[count]] * y3[k2[count]] - x3[k2[count]] * y3[j1[count]];
x4[2] = x3[l]; //get all sign combinations of cross products
y4[0] = y3[j]; x5[count] = x4[count] * copysign(1.0, x4[count]);
y4[1] = y3[k]; y5[count] = y4[count] * copysign(1.0, x4[count]);
y4[2] = y3[l]; z5[count] = z4[count] * copysign(1.0, x4[count]);
z4[0] = z3[j]; x6[count] = x4[count] * copysign(1.0, y4[count]);
z4[1] = z3[k]; y6[count] = y4[count] * copysign(1.0, y4[count]);
z4[2] = z3[l]; z6[count] = z4[count] * copysign(1.0, y4[count]);
xmean5 = ymean5 = zmean5 = xmean6 = ymean6 = zmean6 = xmean7 = ymean7 = zmean7 = 0; x7[count] = x4[count] * copysign(1.0, z4[count]);
for (m=0;m<3;m++) { y7[count] = y4[count] * copysign(1.0, z4[count]);
//get cross products z7[count] = z4[count] * copysign(1.0, z4[count]);
x44[m] = y4[j1[m]]*z4[k2[m]]-y4[k2[m]]*z4[j1[m]]; //get average cross products
y44[m] = z4[j1[m]]*x4[k2[m]]-z4[k2[m]]*x4[j1[m]]; xmean5 += x5[count];
z44[m] = x4[j1[m]]*y4[k2[m]]-x4[k2[m]]*y4[j1[m]]; ymean5 += y5[count];
x5[m] = x44[m]*copysign(1.0,x44[m]); zmean5 += z5[count];
y5[m] = y44[m]*copysign(1.0,x44[m]); xmean6 += x6[count];
z5[m] = z44[m]*copysign(1.0,x44[m]); ymean6 += y6[count];
x6[m] = x44[m]*copysign(1.0,y44[m]); zmean6 += z6[count];
y6[m] = y44[m]*copysign(1.0,y44[m]); xmean7 += x7[count];
z6[m] = z44[m]*copysign(1.0,y44[m]); ymean7 += y7[count];
x7[m] = x44[m]*copysign(1.0,z44[m]); zmean6 += z7[count];
y7[m] = y44[m]*copysign(1.0,z44[m]); count++;
z7[m] = z44[m]*copysign(1.0,z44[m]);
//get average cross products
xmean5 = xmean5 + x5[m];
ymean5 = ymean5 + y5[m];
zmean5 = zmean5 + z5[m];
xmean6 = xmean6 + x6[m];
ymean6 = ymean6 + y6[m];
zmean6 = zmean6 + z6[m];
xmean7 = xmean7 + x7[m];
ymean7 = ymean7 + y7[m];
zmean6 = zmean6 + z7[m];
}
xmean5 = xmean5/3;
xmean6 = xmean6/3;
xmean7 = xmean7/3;
ymean5 = ymean5/3;
ymean6 = ymean6/3;
ymean7 = ymean7/3;
zmean5 = zmean5/3;
zmean6 = zmean6/3;
zmean7 = zmean7/3;
var5 = var6 = var7 = 0;
//get standard deviations
for (m=0;m<3;m++) {
var5 = var5 + x5[m]*x5[m]-2*x5[m]*xmean5+xmean5*xmean5+y5[m]*y5[m]-2*y5[m]*ymean5+ymean5*ymean5+z5[m]*z5[m]-2*z5[m]*zmean5+zmean5*zmean5;
var6 = var6 + x6[m]*x6[m]-2*x6[m]*xmean6+xmean6*xmean6+y6[m]*y6[m]-2*y6[m]*ymean6+ymean6*ymean6+z6[m]*z6[m]-2*z6[m]*zmean6+zmean6*zmean6;
var7 = var7 + x7[m]*x7[m]-2*x7[m]*xmean7+xmean7*xmean7+y7[m]*y7[m]-2*y7[m]*ymean7+ymean7*ymean7+z7[m]*z7[m]-2*z7[m]*zmean7+zmean7*zmean7;
}
//choose minimum standard deviation
if (var5 < S0) {
S0 = var5;
BPV[i][0] = (x5[0]+x5[1]+x5[2])/3;
BPV[i][1] = (y5[0]+y5[1]+x5[2])/3;
BPV[i][2] = (z5[0]+z5[1]+z5[2])/3;
}
if (var6 < S0) {
S0 = var6;
BPV[i][0] = (x6[0]+x6[1]+x6[2])/3;
BPV[i][1] = (y6[0]+y6[1]+x6[2])/3;
BPV[i][2] = (z6[0]+z6[1]+z6[2])/3;
}
if (var7 < S0) {
S0 = var7;
BPV[i][0] = (x7[0]+x7[1]+x7[2])/3;
BPV[i][1] = (y7[0]+y7[1]+x7[2])/3;
BPV[i][2] = (z7[0]+z7[1]+z7[2])/3;
}
}
}
} }
//if there are less than two ~180 degree bond angles, the algorithm returns null }
} else BPV[i][0] = BPV[i][1] = BPV[i][2] = 0.0; if (count > 0) {
xmean5 /= count;
xmean6 /= count;
xmean7 /= count;
ymean5 /= count;
ymean6 /= count;
ymean7 /= count;
zmean5 /= count;
zmean6 /= count;
zmean7 /= count;
}
var5 = var6 = var7 = 0.0;
//find standard deviations
for (j = 0; j < count; j++) {
var5 = var5 + x5[j] * x5[j] - 2 * x5[j] * xmean5 + xmean5 * xmean5 + y5[j] * y5[j] -
2 * y5[j] * ymean5 + ymean5 * ymean5 + z5[j] * z5[j] - 2 * z5[j] * zmean5 +
zmean5 * zmean5;
var6 = var6 + x6[j] * x6[j] - 2 * x6[j] * xmean6 + xmean6 * xmean6 + y6[j] * y6[j] -
2 * y6[j] * ymean6 + ymean6 * ymean6 + z6[j] * z6[j] - 2 * z6[j] * zmean6 +
zmean6 * zmean6;
var7 = var7 + x7[j] * x7[j] - 2 * x7[j] * xmean7 + xmean7 * xmean7 + y7[j] * y7[j] -
2 * y7[j] * ymean7 + ymean7 * ymean7 + z7[j] * z7[j] - 2 * z7[j] * zmean7 +
zmean7 * zmean7;
}
//select sign combination with minimum standard deviation
if (var5 < var6) {
if (var5 < var7)
value = 0;
else
value = 2;
} else if (var6 < var7) {
value = 1;
} else
value = 2;
//BPV is average of cross products of all neighbor vectors which are part of 180 degree angles
BPV[i][0] = 0;
BPV[i][1] = 0;
BPV[i][2] = 0;
for (k = 0; k < count; k++) {
if (value == 0) {
BPV[i][0] = BPV[i][0] + x5[k];
BPV[i][1] = BPV[i][1] + y5[k];
BPV[i][2] = BPV[i][2] + z5[k];
} else if (value == 1) {
BPV[i][0] = BPV[i][0] + x6[k];
BPV[i][1] = BPV[i][1] + y6[k];
BPV[i][2] = BPV[i][2] + z6[k];
} else {
BPV[i][0] = BPV[i][0] + x7[k];
BPV[i][1] = BPV[i][1] + y7[k];
BPV[i][2] = BPV[i][2] + z7[k];
}
}
}
//for atoms with more than three 180 degree bond angles:
else if (chi[0] > 3) {
double x44[3], y44[3], z44[3], S0;
int l, m;
count = 0;
S0 = 100000;
k2[0] = 0;
k2[1] = 0;
k2[2] = 1;
j1[0] = 1;
j1[1] = 2;
j1[2] = 2;
//algorithm is as above, but now all combinations of three 180 degree angles are compared, and the combination with minimum standard deviation is chosen
for (j = 0; j < chi[0]; j++) {
for (k = j + 1; k < chi[0]; k++) {
for (l = k + 1; l < chi[0]; l++) {
if (k >= chi[0] || l >= chi[0]) continue;
//get unique combination of three neighbor vectors
x4[0] = x3[j];
x4[1] = x3[k];
x4[2] = x3[l];
y4[0] = y3[j];
y4[1] = y3[k];
y4[2] = y3[l];
z4[0] = z3[j];
z4[1] = z3[k];
z4[2] = z3[l];
xmean5 = ymean5 = zmean5 = xmean6 = ymean6 = zmean6 = xmean7 = ymean7 = zmean7 = 0;
for (m = 0; m < 3; m++) {
//get cross products
x44[m] = y4[j1[m]] * z4[k2[m]] - y4[k2[m]] * z4[j1[m]];
y44[m] = z4[j1[m]] * x4[k2[m]] - z4[k2[m]] * x4[j1[m]];
z44[m] = x4[j1[m]] * y4[k2[m]] - x4[k2[m]] * y4[j1[m]];
x5[m] = x44[m] * copysign(1.0, x44[m]);
y5[m] = y44[m] * copysign(1.0, x44[m]);
z5[m] = z44[m] * copysign(1.0, x44[m]);
x6[m] = x44[m] * copysign(1.0, y44[m]);
y6[m] = y44[m] * copysign(1.0, y44[m]);
z6[m] = z44[m] * copysign(1.0, y44[m]);
x7[m] = x44[m] * copysign(1.0, z44[m]);
y7[m] = y44[m] * copysign(1.0, z44[m]);
z7[m] = z44[m] * copysign(1.0, z44[m]);
//get average cross products
xmean5 = xmean5 + x5[m];
ymean5 = ymean5 + y5[m];
zmean5 = zmean5 + z5[m];
xmean6 = xmean6 + x6[m];
ymean6 = ymean6 + y6[m];
zmean6 = zmean6 + z6[m];
xmean7 = xmean7 + x7[m];
ymean7 = ymean7 + y7[m];
zmean6 = zmean6 + z7[m];
}
xmean5 = xmean5 / 3;
xmean6 = xmean6 / 3;
xmean7 = xmean7 / 3;
ymean5 = ymean5 / 3;
ymean6 = ymean6 / 3;
ymean7 = ymean7 / 3;
zmean5 = zmean5 / 3;
zmean6 = zmean6 / 3;
zmean7 = zmean7 / 3;
var5 = var6 = var7 = 0;
//get standard deviations
for (m = 0; m < 3; m++) {
var5 = var5 + x5[m] * x5[m] - 2 * x5[m] * xmean5 + xmean5 * xmean5 + y5[m] * y5[m] -
2 * y5[m] * ymean5 + ymean5 * ymean5 + z5[m] * z5[m] - 2 * z5[m] * zmean5 +
zmean5 * zmean5;
var6 = var6 + x6[m] * x6[m] - 2 * x6[m] * xmean6 + xmean6 * xmean6 + y6[m] * y6[m] -
2 * y6[m] * ymean6 + ymean6 * ymean6 + z6[m] * z6[m] - 2 * z6[m] * zmean6 +
zmean6 * zmean6;
var7 = var7 + x7[m] * x7[m] - 2 * x7[m] * xmean7 + xmean7 * xmean7 + y7[m] * y7[m] -
2 * y7[m] * ymean7 + ymean7 * ymean7 + z7[m] * z7[m] - 2 * z7[m] * zmean7 +
zmean7 * zmean7;
}
//choose minimum standard deviation
if (var5 < S0) {
S0 = var5;
BPV[i][0] = (x5[0] + x5[1] + x5[2]) / 3;
BPV[i][1] = (y5[0] + y5[1] + x5[2]) / 3;
BPV[i][2] = (z5[0] + z5[1] + z5[2]) / 3;
}
if (var6 < S0) {
S0 = var6;
BPV[i][0] = (x6[0] + x6[1] + x6[2]) / 3;
BPV[i][1] = (y6[0] + y6[1] + x6[2]) / 3;
BPV[i][2] = (z6[0] + z6[1] + z6[2]) / 3;
}
if (var7 < S0) {
S0 = var7;
BPV[i][0] = (x7[0] + x7[1] + x7[2]) / 3;
BPV[i][1] = (y7[0] + y7[1] + x7[2]) / 3;
BPV[i][2] = (z7[0] + z7[1] + z7[2]) / 3;
}
}
}
}
//if there are less than two ~180 degree bond angles, the algorithm returns null
} else
BPV[i][0] = BPV[i][1] = BPV[i][2] = 0.0;
delete[] x3; delete[] x3;
delete[] y3; delete[] y3;
delete[] z3; delete[] z3;
//normalize BPV: //normalize BPV:
double Mag = sqrt(BPV[i][0]*BPV[i][0] + double Mag = sqrt(BPV[i][0] * BPV[i][0] + BPV[i][1] * BPV[i][1] + BPV[i][2] * BPV[i][2]);
BPV[i][1]*BPV[i][1] + BPV[i][2]*BPV[i][2]);
if (Mag > 0) { if (Mag > 0) {
BPV[i][0] = BPV[i][0]/Mag; BPV[i][0] = BPV[i][0] / Mag;
BPV[i][1] = BPV[i][1]/Mag; BPV[i][1] = BPV[i][1] / Mag;
BPV[i][2] = BPV[i][2]/Mag; BPV[i][2] = BPV[i][2] / Mag;
} }
} else BPV[i][0] = BPV[i][1] = BPV[i][2] = 0.0; } else
BPV[i][0] = BPV[i][1] = BPV[i][2] = 0.0;
} }
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
@ -433,35 +441,37 @@ void ComputeBasalAtom::compute_peratom()
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void ComputeBasalAtom::select(int k, int n, double *arr) void ComputeBasalAtom::select(int k, int n, double *arr)
{ {
int i,ir,j,l,mid; int i, ir, j, l, mid;
double a; double a;
arr--; arr--;
l = 1; l = 1;
ir = n; ir = n;
while (true) { while (true) {
if (ir <= l+1) { if (ir <= l + 1) {
if (ir == l+1 && arr[ir] < arr[l]) std::swap(arr[l],arr[ir]); if (ir == l + 1 && arr[ir] < arr[l]) std::swap(arr[l], arr[ir]);
return; return;
} else { } else {
mid=(l+ir) >> 1; mid = (l + ir) >> 1;
std::swap(arr[mid],arr[l+1]); std::swap(arr[mid], arr[l + 1]);
if (arr[l] > arr[ir]) std::swap(arr[l],arr[ir]); if (arr[l] > arr[ir]) std::swap(arr[l], arr[ir]);
if (arr[l+1] > arr[ir]) std::swap(arr[l+1],arr[ir]); if (arr[l + 1] > arr[ir]) std::swap(arr[l + 1], arr[ir]);
if (arr[l] > arr[l+1]) std::swap(arr[l],arr[l+1]); if (arr[l] > arr[l + 1]) std::swap(arr[l], arr[l + 1]);
i = l+1; i = l + 1;
j = ir; j = ir;
a = arr[l+1]; a = arr[l + 1];
while (true) { while (true) {
do i++; while (arr[i] < a); do i++;
do j--; while (arr[j] > a); while (arr[i] < a);
do j--;
while (arr[j] > a);
if (j < i) break; if (j < i) break;
std::swap(arr[i],arr[j]); std::swap(arr[i], arr[j]);
} }
arr[l+1] = arr[j]; arr[l + 1] = arr[j];
arr[j] = a; arr[j] = a;
if (j >= k) ir = j-1; if (j >= k) ir = j - 1;
if (j <= k) l = i; if (j <= k) l = i;
} }
} }
@ -471,7 +481,7 @@ void ComputeBasalAtom::select(int k, int n, double *arr)
void ComputeBasalAtom::select2(int k, int n, double *arr, int *iarr) void ComputeBasalAtom::select2(int k, int n, double *arr, int *iarr)
{ {
int i,ir,j,l,mid,ia; int i, ir, j, l, mid, ia;
double a; double a;
arr--; arr--;
@ -479,44 +489,46 @@ void ComputeBasalAtom::select2(int k, int n, double *arr, int *iarr)
l = 1; l = 1;
ir = n; ir = n;
while (true) { while (true) {
if (ir <= l+1) { if (ir <= l + 1) {
if (ir == l+1 && arr[ir] < arr[l]) { if (ir == l + 1 && arr[ir] < arr[l]) {
std::swap(arr[l],arr[ir]); std::swap(arr[l], arr[ir]);
std::swap(iarr[l],iarr[ir]); std::swap(iarr[l], iarr[ir]);
} }
return; return;
} else { } else {
mid=(l+ir) >> 1; mid = (l + ir) >> 1;
std::swap(arr[mid],arr[l+1]); std::swap(arr[mid], arr[l + 1]);
std::swap(iarr[mid],iarr[l+1]); std::swap(iarr[mid], iarr[l + 1]);
if (arr[l] > arr[ir]) { if (arr[l] > arr[ir]) {
std::swap(arr[l],arr[ir]); std::swap(arr[l], arr[ir]);
std::swap(iarr[l],iarr[ir]); std::swap(iarr[l], iarr[ir]);
} }
if (arr[l+1] > arr[ir]) { if (arr[l + 1] > arr[ir]) {
std::swap(arr[l+1],arr[ir]); std::swap(arr[l + 1], arr[ir]);
std::swap(iarr[l+1],iarr[ir]); std::swap(iarr[l + 1], iarr[ir]);
} }
if (arr[l] > arr[l+1]) { if (arr[l] > arr[l + 1]) {
std::swap(arr[l],arr[l+1]); std::swap(arr[l], arr[l + 1]);
std::swap(iarr[l],iarr[l+1]); std::swap(iarr[l], iarr[l + 1]);
} }
i = l+1; i = l + 1;
j = ir; j = ir;
a = arr[l+1]; a = arr[l + 1];
ia = iarr[l+1]; ia = iarr[l + 1];
while (true) { while (true) {
do i++; while (arr[i] < a); do i++;
do j--; while (arr[j] > a); while (arr[i] < a);
do j--;
while (arr[j] > a);
if (j < i) break; if (j < i) break;
std::swap(arr[i],arr[j]); std::swap(arr[i], arr[j]);
std::swap(iarr[i],iarr[j]); std::swap(iarr[i], iarr[j]);
} }
arr[l+1] = arr[j]; arr[l + 1] = arr[j];
arr[j] = a; arr[j] = a;
iarr[l+1] = iarr[j]; iarr[l + 1] = iarr[j];
iarr[j] = ia; iarr[j] = ia;
if (j >= k) ir = j-1; if (j >= k) ir = j - 1;
if (j <= k) l = i; if (j <= k) l = i;
} }
} }
@ -528,6 +540,6 @@ void ComputeBasalAtom::select2(int k, int n, double *arr, int *iarr)
double ComputeBasalAtom::memory_usage() double ComputeBasalAtom::memory_usage()
{ {
double bytes = 3*nmax * sizeof(double); double bytes = 3 * nmax * sizeof(double);
return bytes; return bytes;
} }