diff --git a/examples/snap/compute.snap.dat b/examples/snap/compute.snap.dat new file mode 100644 index 0000000000..60aa40fdd6 --- /dev/null +++ b/examples/snap/compute.snap.dat @@ -0,0 +1,59 @@ +# Time-averaged data for fix snap +# TimeStep Number-of-rows +# Row c_snap[1] c_snap[2] c_snap[3] c_snap[4] c_snap[5] c_snap[6] c_snap[7] c_snap[8] c_snap[9] c_snap[10] c_snap[11] +0 55 +1 0 0 0 0 0 3.12659e+06 1.91282e+06 1.01756e+06 1.18149e+06 419003 2775.75 +2 0 0 0 0 0 0 -2617.97 -11804.8 -32003.5 -14156.5 -126.705 +3 0 0 0 0 0 0 -2414.16 -4239.67 -6275.15 -3852.23 -118.927 +4 0 0 0 0 0 0 2529.98 3883.7 6245.75 2522.89 103.66 +5 0 0 0 0 0 0 411.847 604.579 57.0959 1095.67 -188.806 +6 0 0 0 0 0 0 1541.86 4697.43 11841.7 5519.43 275.079 +7 0 0 0 0 0 0 -2870.68 -1447.5 4412.24 1032.92 -63.9586 +8 0 0 0 0 0 0 1193.62 7012.92 20475.9 9007.1 230.377 +9 0 0 0 0 0 0 4848.36 11241.9 22593.7 11630.3 42.8991 +10 0 0 0 0 0 0 -1770.07 -2679.25 -3788.5 -2555.62 -135.264 +11 0 0 0 0 0 0 -4969.62 -8016.32 -11201.8 -7220.33 -85.5022 +12 0 0 0 0 0 0 1641.76 3596.16 7806.47 3219.57 40.8509 +13 0 0 0 0 0 0 325.571 4349.75 13049 5826.43 27.2534 +14 0 0 0 0 0 0 5920.17 5611.27 846.546 2245.23 83.7477 +15 0 0 0 0 0 0 -888.529 -848.965 -1874.49 -290.268 -68.0047 +16 0 0 0 0 0 0 -1916.74 67.9945 4784.3 2143.56 -39.6058 +17 0 0 0 0 0 0 -4098.57 -10375.2 -22007.6 -10355 -200.101 +18 0 0 0 0 0 0 -2284.58 -6551.33 -15184.8 -7117.19 -67.4731 +19 0 0 0 0 0 0 -2737.86 -632.669 6669.64 2094.01 52.5289 +20 0 0 0 0 0 0 -2329.4 -41.9068 7566.17 1913.97 100.188 +21 0 0 0 0 0 0 -444.112 -2754.7 -8428.65 -3849.65 -122.932 +22 0 0 0 0 0 0 -70.5051 111.212 854.264 255.733 65.2259 +23 0 0 0 0 0 0 3554.61 12874.2 31397 14566.8 47.5973 +24 0 0 0 0 0 0 1865.24 2108.07 1180.27 1465.26 91.3443 +25 0 0 0 0 0 0 -889.973 2561.32 11256.4 4537.35 77.4022 +26 0 0 0 0 0 0 3550.36 106.913 -9710.14 -2944.98 144.241 +27 0 0 0 0 0 0 -4712.47 -8838.63 -14464.9 -8091.56 -224.069 +28 0 0 0 0 0 0 -2024.94 -4432.38 -9505.05 -4018.8 -207.602 +29 0 0 0 0 0 0 2379.69 4724.47 7670.76 5006.86 -23.6309 +30 0 0 0 0 0 0 376.992 1771.26 5976.85 2024.35 134.961 +31 0 0 0 0 0 0 1237.27 -1519.65 -9085.33 -3530.88 -43.4288 +32 0 0 0 0 0 0 583.161 6064.47 18404.5 7643.32 243.05 +33 0 0 0 0 0 0 -2538.86 -2021.15 691.987 -389.262 -141.239 +34 0 0 0 0 0 0 2885.38 5612.51 9715.93 5772.93 193.908 +35 0 0 0 0 0 0 -6048.23 -11209.3 -18774.1 -10567.4 -252.412 +36 0 0 0 0 0 0 -1418.32 -3619.88 -5764.64 -4231.84 203.031 +37 0 0 0 0 0 0 3007.44 1474.23 -3713.21 -994.284 140.462 +38 0 0 0 0 0 0 4888.42 4654.63 805.35 2190.37 43.3575 +39 0 0 0 0 0 0 969.58 3277.56 6218.65 3924.82 -58.9942 +40 0 0 0 0 0 0 2987.73 4234.51 5529.54 3085.54 43.2781 +41 0 0 0 0 0 0 810.067 -1872.94 -8730.18 -3125.43 -210.33 +42 0 0 0 0 0 0 2844.79 2986.48 1115.95 1588.01 123.161 +43 0 0 0 0 0 0 134.538 -4097.82 -14380.1 -6204.27 -19.7911 +44 0 0 0 0 0 0 -2999.2 -2447.09 1548.16 -1098.43 162.086 +45 0 0 0 0 0 0 -2288.5 -5930.54 -12773.2 -6503.71 -200.232 +46 0 0 0 0 0 0 -2625.62 -6290.98 -12970.9 -6562.73 -182.126 +47 0 0 0 0 0 0 -228.949 4114.07 13655.9 5798.77 32.8425 +48 0 0 0 0 0 0 2900.97 5126.05 7340.27 4953.94 90.5452 +49 0 0 0 0 0 0 1798.49 -1194.98 -9074.02 -3404.76 -11.9431 +50 0 0 0 0 0 0 -3.09692e+06 -3.518e+06 -4.33318e+06 -2.30338e+06 1.32116e+08 +51 0 0 0 0 0 0 -3.10721e+06 -3.53165e+06 -4.34977e+06 -2.31581e+06 1.28785e+08 +52 0 0 0 0 0 0 -3.10871e+06 -3.53788e+06 -4.36295e+06 -2.32103e+06 1.4248e+08 +53 0 0 0 0 0 0 3585.35 6805.98 11450.9 6458.62 914589 +54 0 0 0 0 0 0 -6674.27 -11551.6 -17884.1 -10474.7 -2.08251e+06 +55 0 0 0 0 0 0 -11913.9 -22733.1 -38858.2 -21261 -7.73337e+06 diff --git a/examples/snap/in.snap.compute b/examples/snap/in.snap.compute index b0c7314882..005ba8b1a3 100644 --- a/examples/snap/in.snap.compute +++ b/examples/snap/in.snap.compute @@ -3,7 +3,7 @@ # initialize simulation variable nsteps index 0 -variable nrep equal 1 +variable nrep equal 2 variable a equal 2.0 units metal @@ -81,7 +81,7 @@ thermo 100 # test output: 1: total potential energy # 2: xy component of stress tensor -# 3: Sum(B_{000}^i, all i of type 2) +# 3: Sum(B_{000}^i, all i of type 2) # 4: xz component of Sum(Sum(r_j*dB_{222}^i/dR[j]), all i of type 2), all j) # 5: y component of -Sum(d(B_{222}^i)/dR[2]), all i of type 2) # @@ -89,7 +89,7 @@ thermo 100 thermo_style custom & pe pxy c_bsum2[1] c_vbsum[55] v_db_2_25 & - c_snap[1][11] c_snap[13][11] c_snap[1][6] c_snap[12][10] c_snap[6][10] + c_snap[1][11] c_snap[13][11] c_snap[1][6] c_snap[12][10] c_snap[6][10] thermo_modify norm no # dump mydump_db all custom 1000 dump_db id c_db[*] diff --git a/python/lammps/mliap/__init__.py b/python/lammps/mliap/__init__.py index 57fe97d803..bbb5bb51a5 100644 --- a/python/lammps/mliap/__init__.py +++ b/python/lammps/mliap/__init__.py @@ -4,7 +4,8 @@ # try to improperly start up a new interpreter. import sysconfig import ctypes -library = sysconfig.get_config_vars('INSTSONAME')[0] +#library = sysconfig.get_config_vars('INSTSONAME')[0] +library="/usr/local/Cellar/python@3.10/3.10.2/Frameworks/Python.framework/Versions/3.10/Python" try: pylib = ctypes.CDLL(library) except OSError as e: diff --git a/src/ML-SNAP/compute_snap.cpp b/src/ML-SNAP/compute_snap.cpp index 83f56e338c..70464420c2 100644 --- a/src/ML-SNAP/compute_snap.cpp +++ b/src/ML-SNAP/compute_snap.cpp @@ -56,6 +56,7 @@ ComputeSnap::ComputeSnap(LAMMPS *lmp, int narg, char **arg) : bzeroflag = 1; quadraticflag = 0; bikflag = 0; + dbirjflag = 0; chemflag = 0; bnormflag = 0; wselfallflag = 0; @@ -81,6 +82,7 @@ ComputeSnap::ComputeSnap(LAMMPS *lmp, int narg, char **arg) : memory->create(cutsq,ntypes+1,ntypes+1,"snap:cutsq"); for (int i = 1; i <= ntypes; i++) { cut = 2.0*radelem[i]*rcutfac; + printf("cut: %f\n", cut); if (cut > cutmax) cutmax = cut; cutsq[i][i] = cut*cut; for (int j = i+1; j <= ntypes; j++) { @@ -147,6 +149,11 @@ ComputeSnap::ComputeSnap(LAMMPS *lmp, int narg, char **arg) : error->all(FLERR,"Illegal compute snap command"); bikflag = atoi(arg[iarg+1]); iarg += 2; + } else if (strcmp(arg[iarg],"dbirjflag") == 0) { + if (iarg+2 > narg) + error->all(FLERR,"Illegal compute snap command"); + dbirjflag = atoi(arg[iarg+1]); + iarg += 2; } else if (strcmp(arg[iarg],"switchinnerflag") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal compute snap command"); @@ -194,7 +201,9 @@ ComputeSnap::ComputeSnap(LAMMPS *lmp, int narg, char **arg) : natoms = atom->natoms; bik_rows = 1; if (bikflag) bik_rows = natoms; - size_array_rows = bik_rows+ndims_force*natoms+ndims_virial; + //size_array_rows = bik_rows+ndims_force*natoms+ndims_virial; + dbirj_rows = ndims_force*natoms; + size_array_rows = bik_rows+dbirj_rows+ndims_virial; size_array_cols = nperdim*atom->ntypes+1; lastcol = size_array_cols-1; @@ -222,6 +231,14 @@ ComputeSnap::~ComputeSnap() memory->destroy(sinnerelem); memory->destroy(dinnerelem); } + + if (dbirjflag){ + printf("dbirj_rows: %d\n", dbirj_rows); + memory->destroy(dbirj); + memory->destroy(nneighs); + memory->destroy(neighsum); + memory->destroy(icounter); + } } /* ---------------------------------------------------------------------- */ @@ -231,8 +248,10 @@ void ComputeSnap::init() if (force->pair == nullptr) error->all(FLERR,"Compute snap requires a pair style be defined"); - if (cutmax > force->pair->cutforce) + if (cutmax > force->pair->cutforce){ + //printf("----- cutmax cutforce: %f %f\n", cutmax, force->pair->cutforce); error->all(FLERR,"Compute snap cutoff is longer than pairwise cutoff"); + } // need an occasional full neighbor list @@ -243,7 +262,7 @@ void ComputeSnap::init() snaptr->init(); // allocate memory for global array - + printf("----- dbirjflag: %d\n", dbirjflag); memory->create(snap,size_array_rows,size_array_cols, "snap:snap"); memory->create(snapall,size_array_rows,size_array_cols, @@ -283,6 +302,13 @@ void ComputeSnap::init_list(int /*id*/, NeighList *ptr) void ComputeSnap::compute_array() { + if (dbirjflag){ + printf("----- dbirjflag true.\n"); + get_dbirj_length(); + printf("----- got dbirj_length\n"); + } + printf("----- cutmax cutforce: %f %f\n", cutmax, force->pair->cutforce); + //else{ int ntotal = atom->nlocal + atom->nghost; invoked_array = update->ntimestep; @@ -295,22 +321,22 @@ void ComputeSnap::compute_array() memory->create(snap_peratom,nmax,size_peratom, "snap:snap_peratom"); } - // clear global array - - for (int irow = 0; irow < size_array_rows; irow++) - for (int icoeff = 0; icoeff < size_array_cols; icoeff++) + printf("size_array_rows: %d\n", size_array_rows); + for (int irow = 0; irow < size_array_rows; irow++){ + for (int icoeff = 0; icoeff < size_array_cols; icoeff++){ + //printf("%d %d\n", irow, icoeff); snap[irow][icoeff] = 0.0; + } + } // clear local peratom array - for (int i = 0; i < ntotal; i++) for (int icoeff = 0; icoeff < size_peratom; icoeff++) { snap_peratom[i][icoeff] = 0.0; } // invoke full neighbor list (will copy or build if necessary) - neighbor->build_one(list); const int inum = list->inum; @@ -324,11 +350,17 @@ void ComputeSnap::compute_array() double** const x = atom->x; const int* const mask = atom->mask; - + //printf("----- inum: %d\n", inum); + //printf("----- NEIGHMASK: %d\n", NEIGHMASK); + int ninside; + int numneigh_sum = 0; for (int ii = 0; ii < inum; ii++) { int irow = 0; if (bikflag) irow = atom->tag[ilist[ii] & NEIGHMASK]-1; + printf("----- ii, ilist, itag, irow: %d %d %d %d\n", ii, ilist[ii] & NEIGHMASK, atom->tag[ilist[ii]], irow); const int i = ilist[ii]; + //printf("----- ii, i: %d %d\n", ii, i); + //printf("----- mask[i] groupbit: %d %d\n", mask[i], groupbit); if (mask[i] & groupbit) { const double xtmp = x[i][0]; @@ -353,7 +385,13 @@ void ComputeSnap::compute_array() // typej = types of neighbors of I within cutoff // note Rij sign convention => dU/dRij = dU/dRj = -dU/dRi - int ninside = 0; + /* + This loop assigns quantities in snaptr. + snaptr is a SNA class instance, see sna.h + + */ + //int ninside = 0; + ninside=0; for (int jj = 0; jj < jnum; jj++) { int j = jlist[jj]; j &= NEIGHMASK; @@ -367,6 +405,7 @@ void ComputeSnap::compute_array() if (chemflag) jelem = map[jtype]; if (rsq < cutsq[itype][jtype]&&rsq>1e-20) { + //printf("cutsq: %f\n", cutsq[itype][jtype]); snaptr->rij[ninside][0] = delx; snaptr->rij[ninside][1] = dely; snaptr->rij[ninside][2] = delz; @@ -382,27 +421,94 @@ void ComputeSnap::compute_array() } } + /* + Now that we have assigned neighbor quantities with previous loop, we are ready to compute things. + Here we compute the wigner functions (U), Z is some other quantity, and bi is bispectrum. + */ snaptr->compute_ui(ninside, ielem); snaptr->compute_zi(); snaptr->compute_bi(ielem); + /* + Looks like this loop computes derivatives. + How does snaptr know what atom I we're dealing with? + I think it only needs neighbor info, and then it goes from there. + */ + //printf("----- Derivative loop - looping over neighbors j.\n"); + printf("----- ninside: %d\n", ninside); // numneighs of I within cutoff for (int jj = 0; jj < ninside; jj++) { + //printf("----- jj: %d\n", jj); const int j = snaptr->inside[jj]; + //printf("----- jj, j, jtag: %d %d %d\n", jj, j, atom->tag[j]); + //int dbirj_row_indx = 3*neighsum[i] + 3*jj ; // THIS IS WRONG, SEE NEXT LINE. + //int dbirj_row_indx = 3*neighsum[j] + 3*i_indx; // need to get i_indx. + // How to get i_indx? + /* + i_indx must start at zero and end at (nneighs[j]-1). + We can start a counter for each atom j. + Maybe this icounter can serve as an index for i as a neighbor of j. + icounter starts at zero and ends at (nneighs[j]-1). + */ + //icounter[atom->tag[j]-1] += 1; + int dbirj_row_indx = 3*neighsum[atom->tag[j]-1] + 3*icounter[atom->tag[j]-1] ; // THIS IS WRONG, SEE NEXT VAR. + printf("jtag, icounter, dbirj_row_indx: %d, %d, %d %d %d\n", atom->tag[j], icounter[atom->tag[j]-1], dbirj_row_indx+0, dbirj_row_indx+1, dbirj_row_indx+2); + icounter[atom->tag[j]-1] += 1; + /* + j is an atom index starting from 0. + Use atom->tag[j] to get the atom in the box (index starts at 1). + Need to make sure that the order of these ij pairs is the same when multiplying by dE/dD later. + */ + //printf("----- jj, j, jtag: %d %d %d\n", jj, j, atom->tag[j]); snaptr->compute_duidrj(jj); snaptr->compute_dbidrj(); // Accumulate dBi/dRi, -dBi/dRj + /* + snap_peratom[i] has type double * because each atom index has indices for descriptors. + */ double *snadi = snap_peratom[i]+typeoffset_local; double *snadj = snap_peratom[j]+typeoffset_local; - + //printf("----- ncoeff: %d\n", ncoeff); + //printf("snadi: %f %f %f %f %f\n", snadi[0], snadi[1], snadi[2], snadi[3], snadi[4]); + //printf("----- typeoffset_local: %d\n", typeoffset_local); + //printf("snadi: "); + //for (int s=0; s<(ncoeff*3); s++){ + // printf("%f ", snadi[s]); + //} + /* + printf("snadj: "); + for (int s=0; s<(ncoeff*3); s++){ + printf("%f ", snadj[s]); + } + */ + //printf("\n"); for (int icoeff = 0; icoeff < ncoeff; icoeff++) { + //printf("----- dblist[icoeff]: %f %f %f\n", snaptr->dblist[icoeff][0], snaptr->dblist[icoeff][1], snaptr->dblist[icoeff][2]); + /* + I think these are the descriptor derivatives. + Desriptor derivatives wrt atom i. + What exactly is being summed here? + This is a loop over descriptors or coeff k. + + */ snadi[icoeff] += snaptr->dblist[icoeff][0]; snadi[icoeff+yoffset] += snaptr->dblist[icoeff][1]; snadi[icoeff+zoffset] += snaptr->dblist[icoeff][2]; + /* + Descriptor derivatives wrt atom j + */ snadj[icoeff] -= snaptr->dblist[icoeff][0]; snadj[icoeff+yoffset] -= snaptr->dblist[icoeff][1]; snadj[icoeff+zoffset] -= snaptr->dblist[icoeff][2]; + + + if (dbirjflag){ + dbirj[dbirj_row_indx+0][icoeff] = snaptr->dblist[icoeff][0]; + dbirj[dbirj_row_indx+1][icoeff] = snaptr->dblist[icoeff][1]; + dbirj[dbirj_row_indx+2][icoeff] = snaptr->dblist[icoeff][2]; + } + } if (quadraticflag) { @@ -453,9 +559,11 @@ void ComputeSnap::compute_array() } } - } + } // for (int jj = 0; jj < ninside; jj++) + //printf("---- irow after jj loop: %d\n", irow); // Accumulate Bi + //printf("----- Accumulate Bi.\n"); // linear contributions @@ -476,18 +584,38 @@ void ComputeSnap::compute_array() } } } + + numneigh_sum += ninside; + } // for (int ii = 0; ii < inum; ii++) + + // Check icounter. + for (int ii = 0; ii < inum; ii++) { + const int i = ilist[ii]; + if (mask[i] & groupbit) { + printf("icounter[i]: %d\n", icounter[i]); + } } + //printf("----- Accumulate bispecturm force contributions to global array.\n"); // accumulate bispectrum force contributions to global array - + //printf("----- ntotal, nmax, natoms: %d %d %d\n", ntotal, nmax, atom->natoms); for (int itype = 0; itype < atom->ntypes; itype++) { const int typeoffset_local = ndims_peratom*nperdim*itype; const int typeoffset_global = nperdim*itype; + //printf("----- nperdim: %d\n", nperdim); + /*nperdim = ncoeff set previsouly*/ for (int icoeff = 0; icoeff < nperdim; icoeff++) { + //printf("----- icoeff: %d\n", icoeff); for (int i = 0; i < ntotal; i++) { double *snadi = snap_peratom[i]+typeoffset_local; int iglobal = atom->tag[i]; + if (icoeff==4){ + if ( (snadi[icoeff] != 0.0) || (snadi[icoeff+yoffset] != 0.0) || (snadi[icoeff+zoffset] != 0.0) ){ + //printf("%d %d %f %f %f\n", i, iglobal, snadi[icoeff], snadi[icoeff+yoffset], snadi[icoeff+zoffset]); + } + } int irow = 3*(iglobal-1)+bik_rows; + //printf("----- snadi[icoeff]: %f\n", snadi[icoeff]); snap[irow++][icoeff+typeoffset_global] += snadi[icoeff]; snap[irow++][icoeff+typeoffset_global] += snadi[icoeff+yoffset]; snap[irow][icoeff+typeoffset_global] += snadi[icoeff+zoffset]; @@ -495,13 +623,20 @@ void ComputeSnap::compute_array() } } + //printf("----- Accumulate forces to global array.\n"); + /* + These are the last columns. + */ // accumulate forces to global array for (int i = 0; i < atom->nlocal; i++) { int iglobal = atom->tag[i]; int irow = 3*(iglobal-1)+bik_rows; + //printf("---- irow: %d\n", irow); snap[irow++][lastcol] = atom->f[i][0]; + //printf("---- irow: %d\n", irow); snap[irow++][lastcol] = atom->f[i][1]; + //printf("---- irow: %d\n", irow); snap[irow][lastcol] = atom->f[i][2]; } @@ -531,6 +666,9 @@ void ComputeSnap::compute_array() snapall[irow++][lastcol] = c_virial->vector[4]; snapall[irow][lastcol] = c_virial->vector[3]; + //}// else + + printf("----- End of compute_array.\n"); } /* ---------------------------------------------------------------------- @@ -567,6 +705,111 @@ void ComputeSnap::dbdotr_compute() } } +/* ---------------------------------------------------------------------- + compute dbirj length +------------------------------------------------------------------------- */ + +void ComputeSnap::get_dbirj_length() +{ + // invoke full neighbor list (will copy or build if necessary) + neighbor->build_one(list); + //memory->destroy(snap); + //memory->destroy(snapall); + dbirj_rows = 0; + const int inum = list->inum; + const int* const ilist = list->ilist; + const int* const numneigh = list->numneigh; + int** const firstneigh = list->firstneigh; + int * const type = atom->type; + const int* const mask = atom->mask; + double** const x = atom->x; + //printf("----- inum: %d\n", inum); + memory->create(neighsum, inum, "snap:neighsum"); + memory->create(nneighs, inum, "snap:nneighs"); + memory->create(icounter, inum, "snap:icounter"); + for (int ii = 0; ii < inum; ii++) { + const int i = ilist[ii]; + if (mask[i] & groupbit) { + icounter[i]=0; + neighsum[i] = 0; + nneighs[i] = 0; + const double xtmp = x[i][0]; + const double ytmp = x[i][1]; + const double ztmp = x[i][2]; + const int itype = type[i]; + const int* const jlist = firstneigh[i]; + const int jnum = numneigh[i]; + //printf("----- jnum: %d\n", jnum); + int jnum_cutoff = 0; + for (int jj = 0; jj < jnum; jj++) { + int j = jlist[jj]; + j &= NEIGHMASK; + + const double delx = x[j][0] - xtmp; + const double dely = x[j][1] - ytmp; + const double delz = x[j][2] - ztmp; + const double rsq = delx*delx + dely*dely + delz*delz; + int jtype = type[j]; + + if (rsq < cutsq[itype][jtype]&&rsq>1e-20) { + dbirj_rows += 1; //jnum + 1; + jnum_cutoff += 1; + nneighs[i]+=1; + } + } + //printf("----- jnum_cutoff: %d\n", jnum_cutoff); + } + } + + dbirj_rows *= ndims_force; + + // Loop over all atoms again to calculate neighsum. + for (int ii = 0; ii < inum; ii++) { + const int i = ilist[ii]; + if (mask[i] & groupbit) { + //printf("nneighs[i]: %d\n", nneighs[i]); + //neighsum[i] = 0; + //printf("i nneighs[i]: %d %d\n", i, nneighs[i]); + if (i==0){ + neighsum[i]=0; + } + else{ + for (int jj=0; jj < ii; jj++){ + const int j = ilist[jj]; + if (mask[j] & groupbit) { + //printf(" j nneighs[j-1]: %d %d\n", j, nneighs[j]); + neighsum[i] += nneighs[j]; + } + } + } + } + //printf("%d\n", neighsum[i]); + } + + memory->create(dbirj, dbirj_rows, ncoeff, "snap:dbirj"); + // Set size array rows which now depends on dbirj_rows. + //size_array_rows = bik_rows+dbirj_rows+ndims_virial; + //printf("----- dbirj_rows: %d\n", dbirj_rows); + //printf("----- end of dbirj length.\n"); + /* + memory->create(snap,size_array_rows,size_array_cols, + "snap:snap"); + memory->create(snapall,size_array_rows,size_array_cols, + "snap:snapall"); + array = snapall; + */ +} + +/* ---------------------------------------------------------------------- + compute array length +------------------------------------------------------------------------- */ + +double ComputeSnap::compute_scalar() +{ + if (dbirjflag) get_dbirj_length(); + return size_array_rows; +} + /* ---------------------------------------------------------------------- memory usage ------------------------------------------------------------------------- */ diff --git a/src/ML-SNAP/compute_snap.h b/src/ML-SNAP/compute_snap.h index bc0670e2c7..d077b8f463 100644 --- a/src/ML-SNAP/compute_snap.h +++ b/src/ML-SNAP/compute_snap.h @@ -31,6 +31,7 @@ class ComputeSnap : public Compute { void init() override; void init_list(int, class NeighList *) override; void compute_array() override; + double compute_scalar() override; double memory_usage() override; private: @@ -52,13 +53,19 @@ class ComputeSnap : public Compute { class SNA *snaptr; double cutmax; int quadraticflag; - int bikflag; - int bik_rows; + //int bikflag; + //int bik_rows; + int bikflag, bik_rows, dbirjflag, dbirj_rows; + double **dbirj; + int *nneighs; // number of neighs inside the snap cutoff. + int *neighsum; + int *icounter; // counting atoms i for each j. Compute *c_pe; Compute *c_virial; void dbdotr_compute(); + void get_dbirj_length(); }; } // namespace LAMMPS_NS diff --git a/src/ML-SNAP/sna.cpp b/src/ML-SNAP/sna.cpp index 33937b9c45..d71eeaa7cb 100644 --- a/src/ML-SNAP/sna.cpp +++ b/src/ML-SNAP/sna.cpp @@ -776,6 +776,7 @@ void SNA::compute_dbidrj() int elem3 = elem_duarray; + //printf("----- idxb_max: %d\n", idxb_max); for (int jjb = 0; jjb < idxb_max; jjb++) { const int j1 = idxb[jjb].j1; const int j2 = idxb[jjb].j2; @@ -1334,6 +1335,9 @@ double SNA::memory_usage() void SNA::create_twojmax_arrays() { + + //printf("----- idxb_max: %d\n", idxb_max); + //printf("----- ntriples: %d\n", ntriples); int jdimpq = twojmax + 2; memory->create(rootpqarray, jdimpq, jdimpq, "sna:rootpqarray"); @@ -1595,4 +1599,3 @@ double SNA::compute_dsfac(double r, double rcut, double sinner, double dinner) return dsfac; } - diff --git a/src/fix_move.cpp b/src/fix_move.cpp index 37e8647671..7cf5567028 100644 --- a/src/fix_move.cpp +++ b/src/fix_move.cpp @@ -35,11 +35,16 @@ #include #include +#include +#include +#include +#include + using namespace LAMMPS_NS; using namespace FixConst; using namespace MathConst; -enum { LINEAR, WIGGLE, ROTATE, VARIABLE, TRANSROT }; +enum { LINEAR, WIGGLE, ROTATE, VARIABLE, TRANSROT, WIGGLE_EIGEN }; enum { EQUAL, ATOM }; #define INERTIA 0.2 // moment of inertia prefactor for ellipsoid @@ -187,7 +192,80 @@ FixMove::FixMove(LAMMPS *lmp, int narg, char **arg) : } else error->all(FLERR, "Illegal fix move command"); - } else + } + // Fix move wiggle_eigen + else if (strcmp(arg[3], "wiggle_eigen") == 0) { + int natoms = atom->natoms; + int nlocal = atom->nlocal; + printf("----- fix move wiggle_eigen initialize -----\n"); + printf("natoms: %d\n", natoms); + memory->create(emat,natoms*3,natoms*3,"move:emat"); + memory->create(freq,natoms*3,"move:freq"); + + // Read EMAT + std::ifstream readfile2; + for (int i = 0; i < natoms*3; i++) { + for (int j = 0; j < natoms*3; j++) { + emat[i][j] = 0.0; + } + } + readfile2.open("../EMAT"); + //readfile2.open("ev_real.txt"); + if (!readfile2.is_open()) { + //printf("ASDFASDF"); + printf("Unable to open ../EMAT\n"); + exit(1); + } + //printf("natoms: %d\n", natoms); + for (int i=0;i<3*natoms;i++){ + for (int j=0;j<3*natoms;j++){ + readfile2>>emat[i][j]; + } + } + + // Read FREQUENCIES + std::ifstream readfile3; + for (int i = 0; i < natoms*3; i++) { + freq[i]=0.0; + } + readfile3.open("FREQUENCIES"); + //readfile2.open("ev_real.txt"); + + if (!readfile3.is_open()) { + printf("Unable to open FREQUENCIES.\n"); + exit(1); + } + //printf("natoms: %d\n", natoms); + for (int i=0;i<3*natoms;i++){ + readfile3>>freq[i]; + } + + + if (narg < 8) error->all(FLERR, "Illegal fix move command"); + iarg = 8; + mstyle = WIGGLE_EIGEN; + if (strcmp(arg[4], "NULL") == 0) + axflag = 0; + else { + axflag = 1; + ax = utils::numeric(FLERR, arg[4], false, lmp); + } + if (strcmp(arg[5], "NULL") == 0) + ayflag = 0; + else { + ayflag = 1; + ay = utils::numeric(FLERR, arg[5], false, lmp); + } + if (strcmp(arg[6], "NULL") == 0) + azflag = 0; + else { + azflag = 1; + az = utils::numeric(FLERR, arg[6], false, lmp); + } + period = utils::numeric(FLERR, arg[7], false, lmp); + if (period <= 0.0) error->all(FLERR, "Illegal fix move command"); + } + else error->all(FLERR, "Illegal fix move command"); // optional args @@ -236,6 +314,10 @@ FixMove::FixMove(LAMMPS *lmp, int narg, char **arg) : if (axflag) ax *= xscale; if (ayflag) ay *= yscale; if (azflag) az *= zscale; + } else if (mstyle == WIGGLE_EIGEN) { + if (axflag) ax *= xscale; + if (ayflag) ay *= yscale; + if (azflag) az *= zscale; } else if (mstyle == ROTATE) { point[0] *= xscale; point[1] *= yscale; @@ -252,7 +334,7 @@ FixMove::FixMove(LAMMPS *lmp, int narg, char **arg) : // set omega_rotate from period - if ((mstyle == WIGGLE) || (mstyle == ROTATE) || (mstyle == TRANSROT)) + if ((mstyle == WIGGLE) || (mstyle == ROTATE) || (mstyle == TRANSROT) || (mstyle == WIGGLE_EIGEN)) omega_rotate = MY_2PI / period; // runit = unit vector along rotation axis @@ -295,10 +377,10 @@ FixMove::FixMove(LAMMPS *lmp, int narg, char **arg) : // AtomVec pointers to retrieve per-atom storage of extra quantities - avec_ellipsoid = dynamic_cast(atom->style_match("ellipsoid")); - avec_line = dynamic_cast(atom->style_match("line")); - avec_tri = dynamic_cast(atom->style_match("tri")); - avec_body = dynamic_cast(atom->style_match("body")); + avec_ellipsoid = dynamic_cast( atom->style_match("ellipsoid")); + avec_line = dynamic_cast( atom->style_match("line")); + avec_tri = dynamic_cast( atom->style_match("tri")); + avec_body = dynamic_cast( atom->style_match("body")); // xoriginal = initial unwrapped positions of atoms // toriginal = initial theta of lines @@ -380,6 +462,11 @@ FixMove::~FixMove() memory->destroy(displace); memory->destroy(velocity); + if (mstyle==WIGGLE_EIGEN){ + memory->destroy(emat); + memory->destroy(freq); + } + delete[] xvarstr; delete[] yvarstr; delete[] zvarstr; @@ -495,7 +582,7 @@ void FixMove::init() velocity = nullptr; if (utils::strmatch(update->integrate_style, "^respa")) - nlevels_respa = (dynamic_cast(update->integrate))->nlevels; + nlevels_respa = (dynamic_cast( update->integrate))->nlevels; } /* ---------------------------------------------------------------------- @@ -527,6 +614,7 @@ void FixMove::initial_integrate(int /*vflag*/) int *tri = atom->tri; int *body = atom->body; int *mask = atom->mask; + int *tag = atom->tag; int nlocal = atom->nlocal; @@ -652,6 +740,142 @@ void FixMove::initial_integrate(int /*vflag*/) // X = P + C + A cos(w*dt) + B sin(w*dt) // V = w R0 cross (A cos(w*dt) + B sin(w*dt)) +} else if (mstyle == WIGGLE_EIGEN) { + //printf("----- Wiggling!-----\n"); + int modeindx = 10; + //printf("%f %f\n", omega_rotate, MY_2PI*freq[modeindx]); + omega_rotate = MY_2PI*freq[modeindx]; + double arg = omega_rotate * delta; + double sine = sin(arg); + double cosine = cos(arg); + + //printf("arg: %f\n", arg); + for (int alpha=0; alpha<3; alpha++){ + for (int i=0; iremap_near(x[i], xold); + } + + } + + + } + + /* + for (int i=0; iremap_near(x[i], xold); + } + + } + */ + + /* + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + xold[0] = x[i][0]; + xold[1] = x[i][1]; + xold[2] = x[i][2]; + + if (axflag) { + v[i][0] = ax * omega_rotate * cosine; + x[i][0] = xoriginal[i][0] + ax * sine; + } else if (rmass) { + dtfm = dtf / rmass[i]; + v[i][0] += dtfm * f[i][0]; + x[i][0] += dtv * v[i][0]; + } else { + dtfm = dtf / mass[type[i]]; + v[i][0] += dtfm * f[i][0]; + x[i][0] += dtv * v[i][0]; + } + + if (ayflag) { + v[i][1] = ay * omega_rotate * cosine; + x[i][1] = xoriginal[i][1] + ay * sine; + } else if (rmass) { + dtfm = dtf / rmass[i]; + v[i][1] += dtfm * f[i][1]; + x[i][1] += dtv * v[i][1]; + } else { + dtfm = dtf / mass[type[i]]; + v[i][1] += dtfm * f[i][1]; + x[i][1] += dtv * v[i][1]; + } + + if (azflag) { + v[i][2] = az * omega_rotate * cosine; + x[i][2] = xoriginal[i][2] + az * sine; + } else if (rmass) { + dtfm = dtf / rmass[i]; + v[i][2] += dtfm * f[i][2]; + x[i][2] += dtv * v[i][2]; + } else { + dtfm = dtf / mass[type[i]]; + v[i][2] += dtfm * f[i][2]; + x[i][2] += dtv * v[i][2]; + } + + domain->remap_near(x[i], xold); + } + } + */ + // for rotate by right-hand rule around omega: + // P = point = vector = point of rotation + // R = vector = axis of rotation + // w = omega of rotation (from period) + // X0 = xoriginal = initial coord of atom + // R0 = runit = unit vector for R + // D = X0 - P = vector from P to X0 + // C = (D dot R0) R0 = projection of atom coord onto R line + // A = D - C = vector from R line to X0 + // B = R0 cross A = vector perp to A in plane of rotation + // A,B define plane of circular rotation around R line + // X = P + C + A cos(w*dt) + B sin(w*dt) + // V = w R0 cross (A cos(w*dt) + B sin(w*dt)) + } else if (mstyle == ROTATE) { double arg = omega_rotate * delta; double cosine = cos(arg); @@ -1303,6 +1527,13 @@ void FixMove::set_arrays(int i) if (ayflag) xoriginal[i][1] -= ay * sine; if (azflag) xoriginal[i][2] -= az * sine; +} else if (mstyle == WIGGLE_EIGEN) { + double arg = omega_rotate * delta; + double sine = sin(arg); + if (axflag) xoriginal[i][0] -= ax * sine; + if (ayflag) xoriginal[i][1] -= ay * sine; + if (azflag) xoriginal[i][2] -= az * sine; + } else if (mstyle == ROTATE) { double a[3], b[3], c[3], d[3], disp[3], ddotr; double arg = -omega_rotate * delta; diff --git a/src/fix_move.h b/src/fix_move.h index e6b253a2a2..8023a66c1b 100644 --- a/src/fix_move.h +++ b/src/fix_move.h @@ -73,6 +73,11 @@ class FixMove : public Fix { int maxatom; double **displace, **velocity; + // fix move wiggle_eigen variables + + double **emat; + double *freq; + class AtomVecEllipsoid *avec_ellipsoid; class AtomVecLine *avec_line; class AtomVecTri *avec_tri; @@ -83,3 +88,53 @@ class FixMove : public Fix { #endif #endif + +/* ERROR/WARNING messages: + +E: Illegal ... command + +Self-explanatory. Check the input script syntax and compare to the +documentation for the command. You can use -echo screen as a +command-line option when running LAMMPS to see the offending line. + +E: Fix move cannot set linear z motion for 2d problem + +Self-explanatory. + +E: Fix move cannot set wiggle z motion for 2d problem + +Self-explanatory. + +E: Fix move cannot rotate around non z-axis for 2d problem + +UNDOCUMENTED + +E: Fix move cannot define z or vz variable for 2d problem + +Self-explanatory. + +E: Zero length rotation vector with fix move + +Self-explanatory. + +E: Variable name for fix move does not exist + +Self-explanatory. + +E: Variable for fix move is invalid style + +Only equal-style variables can be used. + +E: Cannot add atoms to fix move variable + +Atoms can not be added afterwards to this fix option. + +E: Resetting timestep size is not allowed with fix move + +This is because fix move is moving atoms based on elapsed time. + +U: Fix move cannot rotate aroung non z-axis for 2d problem + +Self-explanatory. + +*/