// PROBLEM APPEARS TO BE THAT THE SCALED COORDS DON'T RUN FROM -.5 to +.5 // Something is slightly wrong with the transformations. #include #include #include #include #include #include #include #include"global.h" #define PI 3.1415926536 using namespace std; void boxptrsetup(int i, int j, int k, int l, int m, int n); int main (int argc, const char *argv[]) { N = atoi(argv[1]); Nframes = atoi(argv[2]); R = atoi(argv[3]); S = atoi(argv[4]); T = atoi(argv[5]); const char *infname = argv[6]; if (R < 3 || S < 3 || T < 3) { cout << "Need to use at least three boxes in each direction for code to work. No speedup unless use at least 4 boxes in some direction." << endl; exit(EXIT_FAILURE); } x.resize(N); y.resize(N); z.resize(N); ix.resize(N); iy.resize(N); iz.resize(N); type.resize(N); // matrices for scaling and unscaling coords -- needed to properly apply triclinic PBCs / min-img convention double h[3][3], hinv[3][3]; double det; // determinant of h double scaledxq, scaledyq, scaledzq, scaledxr, scaledyr, scaledzr; double scaleddelx, scaleddely, scaleddelz; double scaledxlo, scaledylo, scaledzlo, scaledxhi, scaledyhi, scaledzhi; double xlo_bound, xhi_bound, ylo_bound, yhi_bound, zlo_bound, zhi_bound; // various helper variables double delx, dely, delz, delsq, del; int f, i, j, k, l, m, n, o, p, q, r, t, u, v; int index, dum; int tox,toy,toz; // used to place atoms in bins double Rd,Sd,Td; // double versions of R,S,T Rd = R; Sd = S; Td = T; // LAMMPS tilt factors: see https://lammps.sandia.gov/doc/dump.html and https://docs.lammps.org/Howto_triclinic.html double xy, xz, yz; // for measuring density fluctuations vector rho; rho.resize(Nframes); // overall number density vector rhosq; rhosq.resize(Nframes); double cellvol, rholocal, phi; // for voronoi analyses -- use only if voronoi volumes are written to the config.dump file double vorov; int nface; // vector avgvorov; avgvorov.resize(Nframes); // vector avgvorovsq; avgvorovsq.resize(Nframes); // for contact analyses vector Nc; Nc.resize(Nframes); // NEXT UP: add scaled avg. contact distance? vector scavgcontdist; scavgcontdist.resize(Nframes); // input, output files char sponge[150]; ifstream fin; fin.open(infname); if (!fin.is_open()) { cout << "Error: couldn't open configs file" << endl; exit(1); } ofstream fout; // timers time_t tstart, tend, tmid; tstart = time(0); // Here do initial box setup vector boxx, boxy, boxz; boxx.resize(N); boxy.resize(N); boxz.resize(N); // set up bins vector< vector< vector > > supercell(R,vector< vector >(S,vector(T))); vector< vector< vector > > boxpops(R,vector< vector >(S,vector(T))); // link bins; each box has a pointer to each of its neighbors, 3D-pbcs assumed for (i = 0; i < R; i++) for (j = 0; j < S; j++) for (k = 0; k < T; k++) { // each box has a pointer to each of its neighbors in 3D for (l = -1; l < 2; l++) for (m = -1; m < 2; m++) for (n = -1; n < 2; n++) { boxptrsetup(i,j,k,l,m,n); supercell[i][j][k].neighbor[l+1][m+1][n+1] = &supercell[pointvec[0]][pointvec[1]][pointvec[2]]; } } // main loop over frames for (f = 0; f < Nframes; f++) { tmid = time(0); // cout << "Frame = " << f << ", time = " << difftime(tmid, tstart) << endl; // read in frame header for (i = 0; i < 5; i++) fin.getline(sponge,150); // note edits made for sheared-triclinic-cell format // DO THESE NEED TO BE CHANGED TO XLO_BOUND, ETC? YES! See https://docs.lammps.org/Howto_triclinic.html // fin >> xlo >> xhi >> xy >> ylo >> yhi >> xz >> zlo >> zhi >> yz; fin >> xlo_bound >> xhi_bound >> xy >> ylo_bound >> yhi_bound >> xz >> zlo_bound >> zhi_bound >> yz; for (i = 0; i < 2; i++) fin.getline(sponge,150); // Set cell edge lengths, cell volume, avg. number density // First define corrected box bounds; see https://docs.lammps.org/Howto_triclinic.html vector shiftxbounds{0,xy,xz,xy+xz}; sort(shiftxbounds.begin(),shiftxbounds.end()); xlo = xlo_bound - shiftxbounds[0]; xhi = xhi_bound - shiftxbounds[3]; ylo = ylo_bound - min(0.0,yz); yhi = yhi_bound - max(0.0,yz); zlo = zlo_bound; zhi = zhi_bound; DX = xhi - xlo; DY = yhi - ylo; DZ = zhi - zlo; // cout << DX << " " << xhi_bound - xlo_bound << " " << DY << " " << yhi_bound - ylo_bound << " " << DZ << " " << zhi_bound - zlo_bound << endl; cellvol = (DX/R)*(DY/S)*(DZ/T); rho[f] = N/(DX*DY*DZ); // estimated phi assumes 50:50 binary 1:1.4 mixture phi = rho[f]*.5*(PI/6)*(1 + 1.4*1.4*1.4); // Next set cell side vectors -- used for PBCs: see https://docs.lammps.org/Howto_triclinic.html // Are these the correct vectors? // Vertical orientation of the vectors a, b, c, i.e. making them the COLUMN vectors of the shape matrix h[0][0] = DX; h[0][1] = xy; h[0][2] = xz; h[1][0] = 0; h[1][1] = DY; h[1][2] = yz; h[2][0] = 0; h[2][1] = 0; h[2][2] = DZ; cout << "Double-precision h matrix: " << endl; cout.precision(15); cout << DX << " " << xy << " " << xz << endl << 0 << " " << DY << " " << yz << endl << 0 << " " << 0 << " " << DZ << endl; det = h[0][0]*(h[1][1]*h[2][2] - h[2][1]*h[1][2]) - h[0][1]*(h[1][0]*h[2][2] - h[2][0]*h[1][2]) + h[0][2]*(h[1][0]*h[2][1] - h[2][0]*h[1][1]); phi = N*(PI/6)*(1 + 1.4*1.4*1.4)/det; // cout << "Frame = " << f << " phi = " << phi << " time = " << difftime(tmid, tstart) << endl; hinv[0][0] = h[1][1]*h[2][2] - h[2][1]*h[1][2]; hinv[0][1] = h[0][2]*h[2][1] - h[0][1]*h[2][2]; hinv[0][2] = h[0][1]*h[1][2] - h[0][2]*h[1][1]; hinv[1][0] = h[1][2]*h[2][0] - h[1][0]*h[2][2]; hinv[1][1] = h[0][0]*h[2][2] - h[0][2]*h[2][0]; hinv[1][2] = h[1][0]*h[0][2] - h[0][0]*h[1][2]; hinv[2][0] = h[1][0]*h[2][1] - h[2][0]*h[1][1]; hinv[2][1] = h[2][0]*h[0][1] - h[0][0]*h[2][1]; hinv[2][2] = h[0][0]*h[1][1] - h[1][0]*h[0][1]; for (i = 0; i < 3; i++) for (j = 0; j < 3; j++) hinv[i][j] /= det; // reads in positions for this frame -- tailored to match Kai's sheared-cohesive-grain dump files for (i = 0; i < N; i++) { fin >> index; // >> dum; index -= 1; fin >> type[index] >> x[index] >> y[index] >> z[index] >> ix[index] >> iy[index] >> iz[index] >> vorov >> nface; // avgvorov[f] += vorov/N; // avgvorovsq[f] += vorov*vorov/N; } fin.getline(sponge,150); // first step of putting atoms in boxes: set atoms' boxx,boxy,boxz, box populations // code needed to be modified to handle triclinic boxes, use scaled coordinates // Q: DO THESE NEED TO BE XLO_BOUND/XHI_BOUND, YLO_BOUND/YHI_BOUND, ZLO_BOUND/ZHI_BOUND instead? // ARE THE SCALED COORDINATES ALL IN THE RANGE (-.5, .5)? // NO!!! For example, one config. gives "Scaled box bounds: -0.500497 0.500772 -0.500138 0.500138 -0.5 0.5" // Fixed one bug, now gives "Scaled box bounds: -0.500634 0.500634 -0.500138 0.500138 -0.5 0.5" // Scaled box is symmetric about the origin but slightly too large along the x, y directions scaledxlo = hinv[0][0]*xlo + hinv[0][1]*ylo + hinv[0][2]*zlo; scaledylo = hinv[1][0]*xlo + hinv[1][1]*ylo + hinv[1][2]*zlo; scaledzlo = hinv[2][0]*xlo + hinv[2][1]*ylo + hinv[2][2]*zlo; scaledxhi = hinv[0][0]*xhi + hinv[0][1]*yhi + hinv[0][2]*zhi; scaledyhi = hinv[1][0]*xhi + hinv[1][1]*yhi + hinv[1][2]*zhi; scaledzhi = hinv[2][0]*xhi + hinv[2][1]*yhi + hinv[2][2]*zhi; cout << "Scaled box bounds: xlo " << scaledxlo << " xhi " << scaledxhi << " ylo " << scaledylo << " yhi " << scaledyhi << " zlo " << scaledzlo << " zhi " << scaledzhi << endl; for (i = 0; i < R; i++) for (j = 0; j < S; j++) for (k = 0; k < T; k++) boxpops[i][j][k] = 0; for (q = 0; q < N; q++) { scaledxq = hinv[0][0]*x[q] + hinv[0][1]*y[q] + hinv[0][2]*z[q]; scaledyq = hinv[1][0]*x[q] + hinv[1][1]*y[q] + hinv[1][2]*z[q]; scaledzq = hinv[2][0]*x[q] + hinv[2][1]*y[q] + hinv[2][2]*z[q]; if (scaledxq < -.5 || scaledxq > .5 || scaledyq < -.5 || scaledyq > .5 || scaledzq < -.5 || scaledzq > .5) { // cout << "q = " << q << " scaled coords " << scaledxq << " " << scaledyq << " " << scaledzq << endl; // exit(1); } tox = static_cast ((scaledxq-scaledxlo)/(1.0/Rd)); if (tox < 0) { tox = R-1; x[u] += DX; } else if (tox >= R) { tox = 0; x[u] -= DX; } boxx[q] = tox; toy = static_cast ((scaledyq-scaledylo)/(1.0/Sd)); if (toy < 0) { toy = S-1; y[u] += DY; } else if (toy >= S) { toy = 0; y[u] -= DY; } boxy[q] = toy; toz = static_cast ((scaledzq-scaledzlo)/(1.0/Td)); if (toz < 0) { toz = T-1; z[u] += DZ; } else if (toz >= T) { toz = 0; z[u] -= DZ; } boxz[q] = toz; boxpops[tox][toy][toz] += 1; } // set up atom lists of length boxpops for (i = 0; i < R; i++) for (j = 0; j < S; j++) for (k = 0; k < T; k++) { supercell[i][j][k].atomlist = new int[boxpops[i][j][k]]; supercell[i][j][k].pos = 0; } // second step of putting atoms in boxes for (u = 0; u < N; u++) { supercell[boxx[u]][boxy[u]][boxz[u]].atomlist[supercell[boxx[u]][boxy[u]][boxz[u]].pos] = u; supercell[boxx[u]][boxy[u]][boxz[u]].pos += 1; } // Calculate density fluctuations for (i = 0; i < R; i++) for (j = 0; j < S; j++) for (k = 0; k < T; k++) { rholocal = boxpops[i][j][k]/cellvol; rhosq[f] += rholocal*rholocal/(R*S*T); } // Calculate Z. Assumes type=1 atoms have radius .5 and type 2 have radius .7 // PBC routine follows Tuckerman textbook for (i = 0; i < R; i++) for (j = 0; j < S; j++) for (k = 0; k < T; k++) for (o = 0; o < boxpops[i][j][k]; o++) { q = supercell[i][j][k].atomlist[o]; // index of the atom I'm currently mapping contacts for // calculate scaled coordinates scaledxq = hinv[0][0]*x[q] + hinv[0][1]*y[q] + hinv[0][2]*z[q]; scaledyq = hinv[1][0]*x[q] + hinv[1][1]*y[q] + hinv[1][2]*z[q]; scaledzq = hinv[2][0]*x[q] + hinv[2][1]*y[q] + hinv[2][2]*z[q]; // loop over neighbors for (l = 0; l < 3; l++) for (m = 0; m < 3; m++) for (n = 0; n < 3; n++) for (p = 0; p < supercell[i][j][k].neighbor[l][m][n]->pos; p++) { r = supercell[i][j][k].neighbor[l][m][n]->atomlist[p]; // index of the atom testing whether in range of atom q if (r > q) { // calculate scaled coordinates scaledxr = hinv[0][0]*x[r] + hinv[0][1]*y[r] + hinv[0][2]*z[r]; scaledyr = hinv[1][0]*x[r] + hinv[1][1]*y[r] + hinv[1][2]*z[r]; scaledzr = hinv[2][0]*x[r] + hinv[2][1]*y[r] + hinv[2][2]*z[r]; scaleddelx = scaledxr - scaledxq; scaleddely = scaledyr - scaledyq; scaleddelz = scaledzr - scaledzq; // min. image. conv in scaled coords. if (scaleddelx > 1.0/2) scaleddelx -= 1; else if (scaleddelx < -1.0/2) scaleddelx += 1; if (scaleddely > 1.0/2) scaleddely -= 1; else if (scaleddely < -1.0/2) scaleddely += 1; if (scaleddelz > 1.0/2) scaleddelz -= 1; else if (scaleddelz < -1.0/2) scaleddelz += 1; // now convert back to unscaled coords to obtain the correct "del" delx = h[0][0]*scaleddelx + h[0][1]*scaleddely + h[0][2]*scaleddelz; dely = h[1][0]*scaleddelx + h[1][1]*scaleddely + h[1][2]*scaleddelz; delz = h[2][0]*scaleddelx + h[2][1]*scaleddely + h[2][2]*scaleddelz; delsq = delx*delx + dely*dely + delz*delz; del = sqrt(delsq); // contact-counting if (del < 1.4) { if (type[q] == 2 && type[r] == 2) Nc[f] += 2; else if ((type[q] == 1 && type[r] == 2) || (type[q] == 2 && type[r] == 1)) if (del < 1.2) Nc[f] += 2; else if (type[q] == 1 && type[r] == 1) if (del < 1) Nc[f] += 2; } } // end of if (r > q) conditional } // end of the l-m-n-p loop } // end of the i-j-k-o loop // delete atoms in linked cells for (i = 0; i < R; i++) for (j = 0; j < S; j++) for (k = 0; k < T; k++) delete [] supercell[i][j][k].atomlist; } // end of f-loop fin.close(); // output Z char ofname[50]; sprintf(ofname,"Zvsstrain"); fout.open(ofname); for (f = 0; f < Nframes; f++) fout << f << " " << (1.0*Nc[f])/N << endl; fout.close(); // timer tend = time(0); cout << "Execution took " << difftime(tend, tstart) << " seconds." << endl; return 0; } // assignment of pointers to neighboring boxes (respecting PBCs) void boxptrsetup(int i, int j, int k, int l, int m, int n) { int o, p, q; // PBCs expressed here: (l, m, n) is the vector from the pointer if ((i + l) == -1) o = R-1; else if ((i + l) == R) o = 0; // handles x rollover else o = i + l; pointvec[0] = o; if ((j + m) == -1) p = S-1; else if ((j + m) == S) p = 0; // handles y rollover else p = j + m; pointvec[1] = p; if ((k + n) == -1) q = T-1; // handles z rollover else if ((k + n) == T) q = 0; else q = k + n; pointvec[2] = q; return; }