diff --git a/examples/snap/README.md b/examples/snap/README.md new file mode 100644 index 0000000000..7c34fe7021 --- /dev/null +++ b/examples/snap/README.md @@ -0,0 +1,9 @@ +See `compute_snap_dgrad.py` for a test that compares the dBi/dRj from compute snap (`dbirjflag=1`) to the sum of dBi/dRj from usual compute snap (`dbirjflag=0`). + +The format of the global array from `dbirjflag=1` is as follows. + +The first N rows belong to bispectrum components for each atom, if we use `bikflag=1`. The first `K` columns correspond to each bispectrum coefficient. The final 3 columns contain reference force components for each atom. + +The rows after the first N rows contain dBi/dRj values for all pairs. These values are arranged in row chunks for each atom `j`, where all the rows in a chunk are associated with the neighbors `i` of `j`, as well as the self-terms where `i=j`. So for atom `j`, the number of rows is equal to the number of atoms within the SNAP cutoff, plus 1 for the `i=j` terms, times 3 for each Cartesian component. The total number of dBi/dRj rows is therefore equal to `N*(Nneigh+1)*3`, and `Nneigh` may be different for each atom. To facilitate with determining which row belong to which atom pair `ij`, the last 3 columns contain indices; the 3rd to last column contains global indices of atoms `i` (the neighbors), the 2nd to last column contains global indices of atoms `j`, and the last column contains an index 0,1,2 for the Cartesian component. Like the `bik` rows, the first `K` columns correspond to each bispectrum coefficient. + +Finally, the first column of the last row contains the reference energy. diff --git a/examples/snap/compute_snap_dgrad.py b/examples/snap/compute_snap_dgrad.py new file mode 100644 index 0000000000..dc06ae1576 --- /dev/null +++ b/examples/snap/compute_snap_dgrad.py @@ -0,0 +1,130 @@ +""" +compute_snap_dgrad.py +Purpose: Demonstrate extraction of descriptor gradient (dB/dR) array from compute snap. + Show that dBi/dRj components summed over neighbors i yields same output as regular compute snap with dbirjflag=0. + This shows that the dBi/dRj components extracted with dbirjflag=1 are correct. +Serial syntax: + python compute_snap_dgrad.py +Parallel syntax: + mpirun -np 2 python compute_snap_dgrad.py +""" + +from __future__ import print_function +import sys +import ctypes +import numpy as np + +# uncomment this if running in parallel via mpi4py +#me = 0 +#from mpi4py import MPI +#me = MPI.COMM_WORLD.Get_rank() +#nprocs = MPI.COMM_WORLD.Get_size() + +from lammps import lammps, LMP_TYPE_ARRAY, LMP_STYLE_GLOBAL +cmds = ["-screen", "none", "-log", "none"] +lmp = lammps(cmdargs=cmds) + +def run_lammps(dbirjflag): + lmp.command("clear") + lmp.command("units metal") + lmp.command("boundary p p p") + lmp.command("atom_modify map hash") + lmp.command(f"lattice bcc {latparam}") + lmp.command(f"region box block 0 {nx} 0 {ny} 0 {nz}") + lmp.command(f"create_box {ntypes} box") + lmp.command(f"create_atoms {ntypes} box") + lmp.command("mass * 180.88") + lmp.command("displace_atoms all random 0.01 0.01 0.01 123456") + # Pair style + snap_options=f'{rcutfac} {rfac0} {twojmax} {radelem1} {radelem2} {wj1} {wj2} rmin0 {rmin0} quadraticflag {quadratic} bzeroflag {bzero} switchflag {switch} bikflag {bikflag} dbirjflag {dbirjflag}' + lmp.command(f"pair_style zero {rcutfac}") + lmp.command(f"pair_coeff * *") + lmp.command(f"pair_style zbl {zblcutinner} {zblcutouter}") + lmp.command(f"pair_coeff * * {zblz} {zblz}") + # set up compute snap generating global array + lmp.command(f"compute snap all snap {snap_options}") + # Run + lmp.command(f"thermo 100") + lmp.command(f"run {nsteps}") + +# Declare simulation/structure variables +nsteps=0 +nrep=2 +latparam=2.0 +ntypes=2 +nx=nrep +ny=nrep +nz=nrep + +# Declare compute snap variables +twojmax=8 +m = (twojmax/2)+1 +rcutfac=1.0 +rfac0=0.99363 +rmin0=0 +radelem1=2.3 +radelem2=2.0 +wj1=1.0 +wj2=0.96 +quadratic=0 +bzero=0 +switch=0 +bikflag=1 +dbirjflag=1 + +# Declare reference potential variables +zblcutinner=4.0 +zblcutouter=4.8 +zblz=73 + +# Number of descriptors +if (twojmax % 2 == 0): + nd = int(m*(m+1)*(2*m+1)/6) +else: + nd = int(m*(m+1)*(m+2)/3) +print(f"Number of descriptors based on twojmax : {nd}") + +# Run lammps with dbirjflag on +run_lammps(1) + +# Get global snap array +lmp_snap = lmp.numpy.extract_compute("snap",0, 2) + +# Extract dBj/dRi (includes dBi/dRi) +natoms = lmp.get_natoms() +fref1 = lmp_snap[0:natoms,-3:].flatten() +eref1 = lmp_snap[-6,0] +dbdr_length = np.shape(lmp_snap)[0]-(natoms)-6 # Length of neighborlist pruned dbdr array +dBdR = lmp_snap[natoms:(natoms+dbdr_length),:] +force_indices = lmp_snap[natoms:(natoms+dbdr_length),-3:].astype(np.int32) + +# Sum over neighbors j for each atom i, like dbirjflag=0 does. +array1 = np.zeros((3*natoms,nd)) +a = 0 +for k in range(0,nd): + for l in range(0,dbdr_length): + j = force_indices[l,0] + i = force_indices[l,1] + #array1[3*(i-1)+a,k] += dBdR[l,k] + array1[3*(i)+a,k] += dBdR[l,k] + a = a+1 + if (a>2): + a=0 + +# Run lammps with dbirjflag off +run_lammps(0) + +# Get global snap array +lmp_snap = lmp.numpy.extract_compute("snap",0, 2) +natoms = lmp.get_natoms() +fref2 = lmp_snap[natoms:(natoms+3*natoms),-1] +eref2 = lmp_snap[0,-1] +array2 = lmp_snap[natoms:natoms+(3*natoms), nd:-1] + +# Sum the arrays obtained from dbirjflag on and off. +summ = array1 + array2 +#np.savetxt("sum.dat", summ) + +print(f"Maximum difference in descriptor sums: {np.max(summ)}") +print(f"Maximum difference in reference forces: {np.max(fref1-fref2)}") +print(f"Difference in reference energy: {np.max(eref1-eref2)}") diff --git a/python/examples/in.lj b/python/examples/in.lj new file mode 100644 index 0000000000..3dc80bdfea --- /dev/null +++ b/python/examples/in.lj @@ -0,0 +1,25 @@ +# 3d Lennard-Jones melt + +units lj +atom_style atomic +atom_modify map array + +lattice fcc 0.8442 +region box block 0 4 0 4 0 4 +create_box 1 box +create_atoms 1 box +mass 1 1.0 + +velocity all create 1.44 87287 loop geom + +pair_style lj/cut 2.5 +pair_coeff 1 1 1.0 1.0 2.5 + +neighbor 0.3 bin +neigh_modify delay 0 every 20 check no + +fix 1 all nve + +variable fx atom fx + +run 10 diff --git a/python/examples/simple.py b/python/examples/simple.py index 8925ce48c0..2f9c191165 100755 --- a/python/examples/simple.py +++ b/python/examples/simple.py @@ -5,10 +5,10 @@ # Purpose: mimic operation of examples/COUPLE/simple/simple.cpp via Python # Serial syntax: simple.py in.lammps -# in.lammps = LAMMPS input script +# in.simple = LAMMPS input script -# Parallel syntax: mpirun -np 4 simple.py in.lammps -# in.lammps = LAMMPS input script +# Parallel syntax: mpirun -np 4 python simple.py in.simple +# in.simple = LAMMPS input script # also need to uncomment mpi4py sections below from __future__ import print_function @@ -27,9 +27,9 @@ infile = sys.argv[1] me = 0 # uncomment this if running in parallel via mpi4py -#from mpi4py import MPI -#me = MPI.COMM_WORLD.Get_rank() -#nprocs = MPI.COMM_WORLD.Get_size() +from mpi4py import MPI +me = MPI.COMM_WORLD.Get_rank() +nprocs = MPI.COMM_WORLD.Get_size() from lammps import lammps lmp = lammps() @@ -122,10 +122,10 @@ if me == 0: print("Gather post scatter subset:", boxlo,boxhi,xy,yz,xz,periodicity,box_change = lmp.extract_box() if me == 0: print("Box info",boxlo,boxhi,xy,yz,xz,periodicity,box_change) -lmp.reset_box([0,0,0],[10,10,8],0,0,0) +#lmp.reset_box([0,0,0],[10,10,8],0,0,0) -boxlo,boxhi,xy,yz,xz,periodicity,box_change = lmp.extract_box() -if me == 0: print("Box info",boxlo,boxhi,xy,yz,xz,periodicity,box_change) +#boxlo,boxhi,xy,yz,xz,periodicity,box_change = lmp.extract_box() +#if me == 0: print("Box info",boxlo,boxhi,xy,yz,xz,periodicity,box_change) # uncomment if running in parallel via mpi4py -#print("Proc %d out of %d procs has" % (me,nprocs), lmp) +print("Proc %d out of %d procs has" % (me,nprocs), lmp) diff --git a/src/ML-SNAP/compute_snap.cpp b/src/ML-SNAP/compute_snap.cpp index 3927703f01..d1b0980dfb 100644 --- a/src/ML-SNAP/compute_snap.cpp +++ b/src/ML-SNAP/compute_snap.cpp @@ -13,7 +13,6 @@ ------------------------------------------------------------------------- */ #include "compute_snap.h" - #include "sna.h" #include "atom.h" #include "update.h" @@ -25,6 +24,7 @@ #include "comm.h" #include "memory.h" #include "error.h" +#include "universe.h" // For MPI #include @@ -82,7 +82,6 @@ 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++) { @@ -201,10 +200,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; dbirj_rows = ndims_force*natoms; size_array_rows = bik_rows+dbirj_rows+ndims_virial; - if (dbirjflag) size_array_cols = nperdim+3; // plus 3 for tag[i], tag[2], and cartesian index + if (dbirjflag) size_array_cols = nperdim+3; // plus 3 for tag[i], tag[j], and cartesian index else size_array_cols = nperdim*atom->ntypes+1; lastcol = size_array_cols-1; @@ -220,41 +218,25 @@ ComputeSnap::~ComputeSnap() { memory->destroy(snap); - //printf("---------------- 1\n"); memory->destroy(snapall); - //printf("---------------- 2\n"); memory->destroy(snap_peratom); - //printf("---------------- 3\n"); memory->destroy(radelem); - //printf("---------------- 4\n"); memory->destroy(wjelem); - //printf("---------------- 5\n"); memory->destroy(cutsq); - //printf("---------------- 6\n"); delete snaptr; - //printf("---------------- 7\n"); if (chemflag) memory->destroy(map); if (switchinnerflag) { memory->destroy(sinnerelem); memory->destroy(dinnerelem); } - //printf("---------------- 8\n"); if (dbirjflag){ - //printf("dbirj_rows: %d\n", dbirj_rows); - //printf("----- destroy dbirj\n"); memory->destroy(dbirj); - //printf("----- 1-1-1-1-1-1\n"); memory->destroy(nneighs); - //printf("----- 2-1-1-1-1-1\n"); memory->destroy(neighsum); - //printf("----- 3-1-1-1-1-1\n"); memory->destroy(icounter); - //printf("----- 4-1-1-1-1-1\n"); memory->destroy(dbiri); - //printf("----- 5-1-1-1-1-1\n"); } - //printf("---------------- 9\n"); } /* ---------------------------------------------------------------------- */ @@ -265,7 +247,6 @@ void ComputeSnap::init() error->all(FLERR,"Compute snap requires a pair style be defined"); 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"); } @@ -279,7 +260,6 @@ void ComputeSnap::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, @@ -319,18 +299,10 @@ void ComputeSnap::init_list(int /*id*/, NeighList *ptr) void ComputeSnap::compute_array() { - //printf(" -----2 dbirjflag: %d\n", dbirjflag); if (dbirjflag){ - //printf("----- dbirjflag true.\n"); get_dbirj_length(); - //printf("----- got dbirj_length\n"); } - //else{ - // printf("----- dbirjflag false.\n"); - //} - //printf("----- cutmax cutforce: %f %f\n", cutmax, force->pair->cutforce); - //else{ int ntotal = atom->nlocal + atom->nghost; invoked_array = update->ntimestep; @@ -343,11 +315,10 @@ void ComputeSnap::compute_array() memory->create(snap_peratom,nmax,size_peratom, "snap:snap_peratom"); } + // clear global array - //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; } } @@ -372,18 +343,13 @@ 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; int dbirj_row_indx; for (int ii = 0; ii < inum; ii++) { int irow = 0; if (bikflag) irow = atom->tag[ilist[ii] & NEIGHMASK]-1; - //printf("----- i, itag: %d %d\n", ilist[ii] & NEIGHMASK, atom->tag[ilist[ii]]); 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]; @@ -408,12 +374,8 @@ void ComputeSnap::compute_array() // typej = types of neighbors of I within cutoff // note Rij sign convention => dU/dRij = dU/dRj = -dU/dRi - /* - This loop assigns quantities in snaptr. - snaptr is a SNA class instance, see sna.h + // assign quantities in snaptr - */ - //int ninside = 0; ninside=0; for (int jj = 0; jj < jnum; jj++) { int j = jlist[jj]; @@ -428,7 +390,6 @@ 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; @@ -444,88 +405,36 @@ 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. - */ + // compute bispectrum for atom i + 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 + // loop over neighbors for descriptors derivatives + 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; + if (dbirjflag){ - 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); + dbirj_row_indx = 3*neighsum[atom->tag[j]-1] + 3*icounter[atom->tag[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. - */ + for (int icoeff = 0; icoeff < ncoeff; icoeff++) { + 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]; @@ -546,11 +455,11 @@ void ComputeSnap::compute_array() dbirj[dbirj_row_indx+2][ncoeff+1] = atom->tag[j]-1; dbirj[dbirj_row_indx+2][ncoeff+2] = 2; } - // Accumulate dBi/dRi = sum (-dBi/dRj) for neighbors j of if i. + // Accumulate dBi/dRi = sum (-dBi/dRj) for neighbors j of if i dbiri[3*(atom->tag[i]-1)+0][icoeff] -= snaptr->dblist[icoeff][0]; dbiri[3*(atom->tag[i]-1)+1][icoeff] -= snaptr->dblist[icoeff][1]; dbiri[3*(atom->tag[i]-1)+2][icoeff] -= snaptr->dblist[icoeff][2]; - // Get last columns + // Get last columns which are i, j, and Cartesian index if (icoeff==(ncoeff-1)){ dbiri[3*(atom->tag[i]-1)+0][ncoeff] = atom->tag[i]-1; dbiri[3*(atom->tag[i]-1)+0][ncoeff+1] = atom->tag[i]-1; @@ -617,18 +526,13 @@ 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 int k; if (dbirjflag) k = 0; else k = typeoffset_global; for (int icoeff = 0; icoeff < ncoeff; icoeff++){ - //printf("----- %d %f\n", icoeff, snaptr->blist[icoeff]); snap[irow][k++] += snaptr->blist[icoeff]; } @@ -649,85 +553,25 @@ void ComputeSnap::compute_array() numneigh_sum += ninside; } // for (int ii = 0; ii < inum; ii++) - //printf("----- bik_rows: %d\n", bik_rows); - //printf("----- bikflag: %d\n", bikflag); - - // 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]); - } - } - */ - - // Sum all the derivatives we calculated to check usual compute snap output. - /* - if (dbirjflag){ - fh_d = fopen("DEBUG", "w"); - int row_indx=0; - double sum; - for (int ii=0; iintypes); - //for (int itype = 0; itype < atom->ntypes; itype++) { int dbiri_indx; int irow; for (int itype=0; itype<1; itype++){ const int typeoffset_local = ndims_peratom*nperdim*itype; const int typeoffset_global = nperdim*itype; - //printf("----- nperdim: %d\n", nperdim); for (int icoeff = 0; icoeff < nperdim; icoeff++) { - //printf("----- icoeff: %d\n", icoeff); dbiri_indx=0; for (int i = 0; i < atom->nlocal; i++) { - //printf("i: %d\n", i); - //int dbiri_indx = 0; - //int irow; for (int jj=0; jjtag[i]-1] + 3*jj; int snap_row_indx = 3*neighsum[atom->tag[i]-1] + 3*(atom->tag[i]-1) + 3*jj; - //printf("snap_row_indx: %d\n", snap_row_indx); - //int irow = dbirj_row_indx+bik_rows; irow = snap_row_indx + bik_rows; - //printf(" row_indx, irow: %d %d\n", dbirj_row_indx, irow); + // x-coordinate snap[irow][icoeff+typeoffset_global] += dbirj[dbirj_row_indx+0][icoeff]; if (icoeff==(ncoeff-1)){ @@ -736,7 +580,7 @@ void ComputeSnap::compute_array() snap[irow][ncoeff+2] += dbirj[dbirj_row_indx+0][ncoeff+2]; } irow++; - //printf(" irow: %d\n", irow); + // y-coordinate snap[irow][icoeff+typeoffset_global] += dbirj[dbirj_row_indx+1][icoeff]; if (icoeff==(ncoeff-1)){ @@ -745,7 +589,7 @@ void ComputeSnap::compute_array() snap[irow][ncoeff+2] += dbirj[dbirj_row_indx+1][ncoeff+2]; } irow++; - //printf(" irow: %d\n", irow); + // z-coordinate snap[irow][icoeff+typeoffset_global] += dbirj[dbirj_row_indx+2][icoeff]; if (icoeff==(ncoeff-1)){ @@ -755,11 +599,11 @@ void ComputeSnap::compute_array() } dbiri_indx = dbiri_indx+3; } + // Put dBi/dRi at end of each dBj/dRi chunk. - //int dbiri_row_indx; + irow = dbiri_indx + bik_rows; - //printf("dbiri_indx: %d\n", dbiri_indx); - //printf("dbiri: %f %f %f\n", dbiri[3*i+0][icoeff], dbiri[3*i+1][icoeff], dbiri[3*i+2][icoeff]); + // x-coordinate snap[irow][icoeff+typeoffset_global] += dbiri[3*i+0][icoeff]; if (icoeff==(ncoeff-1)){ @@ -768,7 +612,7 @@ void ComputeSnap::compute_array() snap[irow][ncoeff+2] += dbiri[3*i+0][ncoeff+2]; } irow++; - //printf(" irow: %d\n", irow); + // y-coordinate snap[irow][icoeff+typeoffset_global] += dbiri[3*i+1][icoeff]; if (icoeff==(ncoeff-1)){ @@ -777,7 +621,8 @@ void ComputeSnap::compute_array() snap[irow][ncoeff+2] += dbiri[3*i+1][ncoeff+2]; } irow++; - //printf(" irow: %d\n", irow); + + // z-coordinate snap[irow][icoeff+typeoffset_global] += dbiri[3*i+2][icoeff]; if (icoeff==(ncoeff-1)){ snap[irow][ncoeff] += dbiri[3*i+2][ncoeff]; @@ -789,31 +634,21 @@ void ComputeSnap::compute_array() } } } - //printf(" Accumulated to global array.\n"); } else{ - //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]; @@ -822,15 +657,11 @@ void ComputeSnap::compute_array() } } - //printf("----- Accumulate forces to global array.\n"); - /* - These are the last columns. - */ // accumulate forces to global array + if (dbirjflag){ + // for dbirjflag=1, put forces at last 3 columns of bik rows for (int i=0; inlocal; i++){ - - int iglobal = atom->tag[i]; snap[iglobal-1][ncoeff+0] = atom->f[i][0]; snap[iglobal-1][ncoeff+1] = atom->f[i][1]; @@ -842,11 +673,8 @@ void ComputeSnap::compute_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]; } } @@ -854,7 +682,7 @@ void ComputeSnap::compute_array() // accumulate bispectrum virial contributions to global array if (dbirjflag){ - + // no virial terms for dbirj yet } else{ dbdotr_compute(); @@ -865,8 +693,8 @@ void ComputeSnap::compute_array() // assign energy to last column if (dbirjflag){ // Assign reference energy right after the dbirj rows, first column. - int irow = bik_rows + dbirj_rows + 3*natoms; // add 3N for the dBi/dRi rows - //printf("irow bik_rows dbirj_rows: %d %d %d\n", irow, bik_rows, dbirj_rows); + // Add 3N for the dBi/dRi rows. + int irow = bik_rows + dbirj_rows + 3*natoms; double reference_energy = c_pe->compute_scalar(); snapall[irow][0] = reference_energy; } @@ -882,7 +710,7 @@ void ComputeSnap::compute_array() c_virial->compute_vector(); if (dbirjflag){ - + // no virial terms for dbirj yet } else{ int irow = 3*natoms+bik_rows; @@ -894,7 +722,6 @@ void ComputeSnap::compute_array() snapall[irow][lastcol] = c_virial->vector[3]; } - //}// else } /* ---------------------------------------------------------------------- @@ -906,14 +733,7 @@ void ComputeSnap::dbdotr_compute() { double **x = atom->x; - int irow0; - if (dbirjflag){ - irow0 = bik_rows+dbirj_rows+(3*natoms); - } - else{ - irow0 = bik_rows+ndims_force*natoms; - } - //int irow0 = bik_rows+ndims_force*natoms; + int irow0 = bik_rows+ndims_force*natoms; // sum over bispectrum contributions to forces // on all particles including ghosts @@ -946,9 +766,11 @@ void ComputeSnap::dbdotr_compute() void ComputeSnap::get_dbirj_length() { + int rank = universe->me; // for MPI debugging + memory->destroy(snap); memory->destroy(snapall); - // invoke full neighbor list (will copy or build if necessary) + // invoke full neighbor list neighbor->build_one(list); dbirj_rows = 0; const int inum = list->inum; @@ -958,12 +780,20 @@ void ComputeSnap::get_dbirj_length() 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"); memory->create(dbiri, 3*atom->nlocal,ncoeff+3, "snap:dbiri"); - for (int ii=0; ii<3*atom->nlocal; ii++){ + */ + memory->create(neighsum, natoms, "snap:neighsum"); + memory->create(nneighs, natoms, "snap:nneighs"); + memory->create(icounter, natoms, "snap:icounter"); + memory->create(dbiri, 3*natoms,ncoeff+3, "snap:dbiri"); + if (atom->nlocal != natoms){ + error->all(FLERR,"Compute snap dbirjflag=1 does not support parallelism."); + } + for (int ii=0; ii<3*natoms; ii++){ for (int icoeff=0; icoeffcreate(dbirj, dbirj_rows, ncoeff+3, "snap:dbirj"); for (int i=0; inlocal; // Add 3*N for dBi/dRi - //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");