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

Build dbidrj array.

Browse Source
This commit is contained in:
rohskopf
2022-05-24 15:33:40 -06:00
parent 223aebe3fb
commit 43048811dd
8 changed files with 628 additions and 29 deletions
Download Patch File Download Diff File Expand all files Collapse all files

59
examples/snap/compute.snap.dat Normal file
View File

@ -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

2
examples/snap/in.snap.compute
View File

@ -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

3
python/lammps/mliap/__init__.py
View File

@ -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:

271
src/ML-SNAP/compute_snap.cpp
View File

@ -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
------------------------------------------------------------------------- */

11
src/ML-SNAP/compute_snap.h
View File

@ -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

5
src/ML-SNAP/sna.cpp
View File

@ -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;
}

247
src/fix_move.cpp
View File

@ -35,11 +35,16 @@
#include <cmath>
#include <cstring>
#include <sstream>
#include <vector>
#include <fstream>
#include <string>
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<AtomVecEllipsoid *>(atom->style_match("ellipsoid"));
avec_line = dynamic_cast<AtomVecLine *>(atom->style_match("line"));
avec_tri = dynamic_cast<AtomVecTri *>(atom->style_match("tri"));
avec_body = dynamic_cast<AtomVecBody *>(atom->style_match("body"));
avec_ellipsoid = dynamic_cast<AtomVecEllipsoid *>( atom->style_match("ellipsoid"));
avec_line = dynamic_cast<AtomVecLine *>( atom->style_match("line"));
avec_tri = dynamic_cast<AtomVecTri *>( atom->style_match("tri"));
avec_body = dynamic_cast<AtomVecBody *>( 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<Respa *>(update->integrate))->nlevels;
nlevels_respa = (dynamic_cast<Respa *>( 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; i<nlocal; i++){
if(mask[i] & groupbit){
xold[0] = x[i][0];
xold[1] = x[i][1];
xold[2] = x[i][2];
if (axflag) {
//int modeindx = 3;
//int alpha = 0;
ax = emat[3*(tag[i]-1)+alpha][modeindx];
v[i][alpha] = ax * omega_rotate * cosine;
x[i][alpha] = xoriginal[i][alpha] + ax * sine;
} else if (rmass) {
dtfm = dtf / rmass[i];
v[i][alpha] += dtfm * f[i][alpha];
x[i][alpha] += dtv * v[i][alpha];
} else {
dtfm = dtf / mass[type[i]];
v[i][alpha] += dtfm * f[i][alpha];
x[i][alpha] += dtv * v[i][alpha];
}
domain->remap_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) {
//int modeindx = 3;
int alpha = 0;
ax = emat[3*(tag[i]-1)+alpha][modeindx];
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];
}
domain->remap_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;

55
src/fix_move.h
View File

@ -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.
*/