diff --git a/doc/src/Commands_compute.rst b/doc/src/Commands_compute.rst index 61c5e83eda..bb0c3e9c2a 100644 --- a/doc/src/Commands_compute.rst +++ b/doc/src/Commands_compute.rst @@ -138,6 +138,8 @@ KOKKOS, o = OPENMP, t = OPT. * :doc:`smd/vol ` * :doc:`snap ` * :doc:`sna/atom ` + * :doc:`sna/grid ` + * :doc:`sna/grid/local ` * :doc:`snad/atom ` * :doc:`snav/atom ` * :doc:`sph/e/atom ` diff --git a/doc/src/Errors_warnings.rst b/doc/src/Errors_warnings.rst index 2d588a1b77..ab06ac523c 100644 --- a/doc/src/Errors_warnings.rst +++ b/doc/src/Errors_warnings.rst @@ -470,6 +470,12 @@ This will most likely cause errors in kinetic fluctuations. *More than one compute sna/atom* Self-explanatory. +*More than one compute sna/grid* + Self-explanatory. + +*More than one compute sna/grid/local* + Self-explanatory. + *More than one compute snad/atom* Self-explanatory. diff --git a/doc/src/Packages_details.rst b/doc/src/Packages_details.rst index b5d704190e..65bab0e5a2 100644 --- a/doc/src/Packages_details.rst +++ b/doc/src/Packages_details.rst @@ -1801,6 +1801,8 @@ computes which analyze attributes of the potential. * src/ML-SNAP: filenames -> commands * :doc:`pair_style snap ` * :doc:`compute sna/atom ` +* :doc:`compute sna/grid ` +* :doc:`compute sna/grid/local ` * :doc:`compute snad/atom ` * :doc:`compute snav/atom ` * examples/snap diff --git a/doc/src/compute.rst b/doc/src/compute.rst index 6fdedbbb95..508c440e78 100644 --- a/doc/src/compute.rst +++ b/doc/src/compute.rst @@ -284,6 +284,8 @@ The individual style names on the :doc:`Commands compute ` pag * :doc:`smd/vol ` - per-particle volumes and their sum in Smooth Mach Dynamics * :doc:`snap ` - gradients of SNAP energy and forces w.r.t. linear coefficients and related quantities for fitting SNAP potentials * :doc:`sna/atom ` - bispectrum components for each atom +* :doc:`sna/grid ` - global array of bispectrum components on a regular grid +* :doc:`sna/grid/local ` - local array of bispectrum components on a regular grid * :doc:`snad/atom ` - derivative of bispectrum components for each atom * :doc:`snav/atom ` - virial contribution from bispectrum components for each atom * :doc:`sph/e/atom ` - per-atom internal energy of Smooth-Particle Hydrodynamics atoms diff --git a/doc/src/compute_sna_atom.rst b/doc/src/compute_sna_atom.rst index 54a6df02a2..deb717d84d 100644 --- a/doc/src/compute_sna_atom.rst +++ b/doc/src/compute_sna_atom.rst @@ -2,6 +2,8 @@ .. index:: compute snad/atom .. index:: compute snav/atom .. index:: compute snap +.. index:: compute sna/grid +.. index:: compute sna/grid/local compute sna/atom command ======================== @@ -15,6 +17,12 @@ compute snav/atom command compute snap command ==================== +compute sna/grid command +======================== + +compute sna/grid/local command +============================== + Syntax """""" @@ -24,6 +32,9 @@ Syntax compute ID group-ID snad/atom rcutfac rfac0 twojmax R_1 R_2 ... w_1 w_2 ... keyword values ... compute ID group-ID snav/atom rcutfac rfac0 twojmax R_1 R_2 ... w_1 w_2 ... keyword values ... compute ID group-ID snap rcutfac rfac0 twojmax R_1 R_2 ... w_1 w_2 ... keyword values ... + compute ID group-ID snap rcutfac rfac0 twojmax R_1 R_2 ... w_1 w_2 ... keyword values ... + compute ID group-ID sna/grid nx ny nz rcutfac rfac0 twojmax R_1 R_2 ... w_1 w_2 ... keyword values ... + compute ID group-ID sna/grid/local nx ny nz rcutfac rfac0 twojmax R_1 R_2 ... w_1 w_2 ... keyword values ... * ID, group-ID are documented in :doc:`compute ` command * sna/atom = style name of this compute command @@ -32,6 +43,7 @@ Syntax * twojmax = band limit for bispectrum components (non-negative integer) * R_1, R_2,... = list of cutoff radii, one for each type (distance units) * w_1, w_2,... = list of neighbor weights, one for each type +* nx, ny, nz = number of grid points in x, y, and z directions (positive integer) * zero or more keyword/value pairs may be appended * keyword = *rmin0* or *switchflag* or *bzeroflag* or *quadraticflag* or *chem* or *bnormflag* or *wselfallflag* or *bikflag* or *switchinnerflag* or *sinner* or *dinner* @@ -78,6 +90,7 @@ Examples compute snap all snap 1.4 0.95 6 2.0 1.0 compute snap all snap 1.0 0.99363 6 3.81 3.83 1.0 0.93 chem 2 0 1 compute snap all snap 1.0 0.99363 6 3.81 3.83 1.0 0.93 switchinnerflag 1 sinner 1.35 1.6 dinner 0.25 0.3 + compute bgrid all sna/grid/local 200 200 200 1.4 0.95 6 2.0 1.0 Description """"""""""" @@ -212,6 +225,46 @@ command: See section below on output for a detailed explanation of the data layout in the global array. +The compute *sna/grid* and *sna/grid/local* commands calculate +bispectrum components for a regular grid of points. +These are calculated from the local density of nearby atoms *i'* +around each grid point, as if there was a central atom *i* +at the grid point. This is useful for characterizing fine-scale +structure in a configuration of atoms, and it is used +in the `MALA package `_ +to build machine-learning surrogates for finite-temperature Kohn-Sham +density functional theory (:ref:`Ellis et al. `) +Neighbor atoms not in the group do not contribute to the +bispectrum components of the grid points. The distance cutoff :math:`R_{ii'}` +assumes that *i* has the same type as the neighbor atom *i'*. + +Compute *sna/grid* calculates a global array containing bispectrum +components for a regular grid of points. +The grid is aligned with the current box dimensions, with the +first point at the box origin, and forming a regular 3d array with +*nx*, *ny*, and *nz* points in the x, y, and z directions. For triclinic +boxes, the array is congruent with the periodic lattice vectors +a, b, and c. The array contains one row for each of the +:math:`nx \times ny \times nz` grid points, looping over the index for *ix* fastest, +then *iy*, and *iz* slowest. Each row of the array contains the *x*, *y*, +and *z* coordinates of the grid point, followed by the bispectrum +components. See section below on output for a detailed explanation of the data +layout in the global array. + +Compute *sna/grid/local* calculates bispectrum components of a regular +grid of points similarly to compute *sna/grid* described above. +However, because the array is local, it contains only rows for grid points +that are local to the processor sub-domain. The global grid +of :math:`nx \times ny \times nz` points is still laid out in space the same as for *sna/grid*, +but grid points are strictly partitioned, so that every grid point appears in +one and only one local array. The array contains one row for each of the +local grid points, looping over the global index *ix* fastest, +then *iy*, and *iz* slowest. Each row of the array contains +the global indexes *ix*, *iy*, and *iz* first, followed by the *x*, *y*, +and *z* coordinates of the grid point, followed by the bispectrum +components. See section below on output for a detailed explanation of the data +layout in the global array. + The value of all bispectrum components will be zero for atoms not in the group. Neighbor atoms not in the group do not contribute to the bispectrum of atoms in the group. @@ -414,6 +467,21 @@ number of columns in the global array generated by *snap* are 31, and 931, respectively, while the number of rows is 1+3\*\ *N*\ +6, where *N* is the total number of atoms. +Compute *sna/grid* evaluates a global array. +The array contains one row for each of the +:math:`nx \times ny \times nz` grid points, looping over the index for *ix* fastest, +then *iy*, and *iz* slowest. Each row of the array contains the *x*, *y*, +and *z* coordinates of the grid point, followed by the bispectrum +components. + +Compute *sna/grid/local* evaluates a local array. +The array contains one row for each of the +local grid points, looping over the global index *ix* fastest, +then *iy*, and *iz* slowest. Each row of the array contains +the global indexes *ix*, *iy*, and *iz* first, followed by the *x*, *y*, +and *z* coordinates of the grid point, followed by the bispectrum +components. + If the *quadratic* keyword value is set to 1, then additional columns are generated, corresponding to the products of all distinct pairs of bispectrum components. If the number of bispectrum components is *K*, @@ -464,8 +532,7 @@ The optional keyword defaults are *rmin0* = 0, .. _Thompson20141: -**(Thompson)** Thompson, Swiler, Trott, Foiles, Tucker, under review, preprint -available at `arXiv:1409.3880 `_ +**(Thompson)** Thompson, Swiler, Trott, Foiles, Tucker, J Comp Phys, 285, 316, (2015). .. _Bartok20101: @@ -486,4 +553,8 @@ of Angular Momentum, World Scientific, Singapore (1987). .. _Cusentino2020: -**(Cusentino)** Cusentino, Wood, and Thompson, J Phys Chem A, xxx, xxxxx, (2020) +**(Cusentino)** Cusentino, Wood, Thompson, J Phys Chem A, 124, 5456, (2020) + +.. _Ellis2021: + +**(Ellis)** Ellis, Fiedler, Popoola, Modine, Stephens, Thompson, Cangi, Rajamanickam, Phys Rev B, 104, 035120, (2021) diff --git a/doc/utils/sphinx-config/false_positives.txt b/doc/utils/sphinx-config/false_positives.txt index cce6e92d05..abc6b47ccd 100644 --- a/doc/utils/sphinx-config/false_positives.txt +++ b/doc/utils/sphinx-config/false_positives.txt @@ -373,6 +373,7 @@ Caltech Camilloni Camiloni Campana +Cangi Cao Capolungo Caro @@ -2681,6 +2682,7 @@ polyhedra Polym polymorphism popen +Popoola Popov popstore Poresag @@ -2824,6 +2826,7 @@ radians Rafferty rahman Rahman +Rajamanickam Ralf Raman ramped diff --git a/examples/snap/in.grid.snap b/examples/snap/in.grid.snap new file mode 100644 index 0000000000..08c95a004f --- /dev/null +++ b/examples/snap/in.grid.snap @@ -0,0 +1,94 @@ +# Demonstrate calculation of SNAP bispectrum descriptors on a grid + +# CORRECTNESS: The two atom positions coincide with two of +# the gridpoints, so c_b[2][1-5] should match c_mygrid[8][4-8]. +# The same is true for compute grid/local c_mygridlocal[8][4-11]. +# Local arrays can not be access directly in the script, +# but they are printed out to file dump.blocal + +variable nrep index 1 +variable a index 3.316 +variable ngrid index 2 + +units metal +atom_modify map hash + +# generate the box and atom positions using a BCC lattice + +variable nx equal ${nrep} +variable ny equal ${nrep} +variable nz equal ${nrep} + +boundary p p p + +lattice custom $a & + a1 1 0 0 & + a2 0 1 0 & + a3 0 0 1 & + basis 0 0 0 & + basis 0.5 0.5 0.5 & + +region box block 0 ${nx} 0 ${ny} 0 ${nz} +create_box 1 box +create_atoms 1 box + +mass 1 180.88 + +# define atom compute and grid compute + +group snapgroup type 1 +variable twojmax equal 2 +variable rcutfac equal 4.67637 +variable rfac0 equal 0.99363 +variable rmin0 equal 0 +variable wj equal 1 +variable radelem equal 0.5 +variable bzero equal 0 +variable quadratic equal 0 +variable switch equal 1 + +variable snap_options string & + "${rcutfac} ${rfac0} ${twojmax} ${radelem} & + ${wj} rmin0 ${rmin0} quadraticflag ${quadratic} & + bzeroflag ${bzero} switchflag ${switch}" + +# build zero potential to satisfy compute sna/atom + +pair_style zero ${rcutfac} +pair_coeff * * + +# define atom and grid computes + +compute b all sna/atom ${snap_options} +compute mygrid all sna/grid grid ${ngrid} ${ngrid} ${ngrid} & + ${snap_options} +compute mygridlocal all sna/grid/local grid ${ngrid} ${ngrid} ${ngrid} & + ${snap_options} + +# define output + +variable B5atom equal c_b[2][5] +variable B5grid equal c_mygrid[8][8] + +variable rmse_global equal "sqrt( & + (c_mygrid[8][1] - x[2])^2 + & + (c_mygrid[8][2] - y[2])^2 + & + (c_mygrid[8][3] - z[2])^2 + & + (c_mygrid[8][4] - c_b[2][1])^2 + & + (c_mygrid[8][5] - c_b[2][2])^2 + & + (c_mygrid[8][6] - c_b[2][3])^2 + & + (c_mygrid[8][7] - c_b[2][4])^2 + & + (c_mygrid[8][8] - c_b[2][5])^2 & + )" + +thermo_style custom step v_B5atom v_B5grid v_rmse_global + +# this is the only way to view the local grid + +dump 1 all local 1000 dump.blocal c_mygridlocal[*] +dump 2 all custom 1000 dump.batom id x y z c_b[*] + +# run + +run 0 + diff --git a/examples/snap/in.grid.tri b/examples/snap/in.grid.tri new file mode 100644 index 0000000000..5283957eb8 --- /dev/null +++ b/examples/snap/in.grid.tri @@ -0,0 +1,114 @@ +# Demonstrate calculation of SNAP bispectrum +# descriptors on a grid for triclinic cell + +# This triclinic cell has 6 times the volume of the single +# unit cell used by in.grid +# and contains 12 atoms. It is a 3x2x1 supercell +# with each unit cell containing 2 atoms and the +# reduced lattice vectors are [1 0 0], [1 1 0], and [1 1 1]. +# The grid is listed in x-fastest order + +# CORRECTNESS: The atom positions coincide with certain +# gridpoints, so c_b[1][1-5] should match c_mygrid[1][4-8] +# and c_b[7][1-5] should match c_mygrid[13][4-8]. +# Local arrays can not be access directly in the script, +# but they are printed out to file dump.blocal.tri + +# Initialize simulation + +variable nrep index 1 +variable a index 3.316 +variable ngrid index 2 + +variable nrepx equal 3*${nrep} +variable nrepy equal 2*${nrep} +variable nrepz equal 1*${nrep} + +variable ngridx equal 3*${ngrid} +variable ngridy equal 2*${ngrid} +variable ngridz equal 1*${ngrid} + +units metal +atom_modify map hash sort 0 0 + +# generate the box and atom positions using a BCC lattice + +variable nx equal ${nrepx} +variable ny equal ${nrepy} +variable nz equal ${nrepz} + +boundary p p p + +lattice custom $a & + a1 1 0 0 & + a2 1 1 0 & + a3 1 1 1 & + basis 0 0 0 & + basis 0.0 0.0 0.5 & + spacing 1 1 1 + +box tilt large +region box prism 0 ${nx} 0 ${ny} 0 ${nz} ${ny} ${nz} ${nz} +create_box 1 box +create_atoms 1 box + +mass 1 180.88 + +# define atom compute and grid compute + +group snapgroup type 1 +variable twojmax equal 2 +variable rcutfac equal 4.67637 +variable rfac0 equal 0.99363 +variable rmin0 equal 0 +variable wj equal 1 +variable radelem equal 0.5 +variable bzero equal 0 +variable quadratic equal 0 +variable switch equal 1 + +variable snap_options string & + "${rcutfac} ${rfac0} ${twojmax} ${radelem} & + ${wj} rmin0 ${rmin0} quadraticflag ${quadratic} & + bzeroflag ${bzero} switchflag ${switch}" + +# build zero potential to satisfy compute sna/atom + +pair_style zero ${rcutfac} +pair_coeff * * + +# define atom and grid computes + +compute b all sna/atom ${snap_options} +compute mygrid all sna/grid grid ${ngridx} ${ngridy} ${ngridz} & + ${snap_options} +compute mygridlocal all sna/grid/local grid ${ngridx} ${ngridy} ${ngridz} & + ${snap_options} + +# define output + +variable B5atom equal c_b[7][5] +variable B5grid equal c_mygrid[13][8] + +# do not compare x,y,z because assignment of ids +# to atoms is not unnique for different processor grids + +variable rmse_global equal "sqrt( & + (c_mygrid[13][4] - c_b[7][1])^2 + & + (c_mygrid[13][5] - c_b[7][2])^2 + & + (c_mygrid[13][6] - c_b[7][3])^2 + & + (c_mygrid[13][7] - c_b[7][4])^2 + & + (c_mygrid[13][8] - c_b[7][5])^2 & + )" + +thermo_style custom step v_B5atom v_B5grid v_rmse_global + +# this is the only way to view the local grid + +dump 1 all local 1000 dump.blocal.tri c_mygridlocal[*] +dump 2 all custom 1000 dump.batom.tri id x y z c_b[*] + +# run + +run 0 + diff --git a/examples/snap/log.15Jun22.grid.snap.g++.1 b/examples/snap/log.15Jun22.grid.snap.g++.1 new file mode 100644 index 0000000000..ec2026b16e --- /dev/null +++ b/examples/snap/log.15Jun22.grid.snap.g++.1 @@ -0,0 +1,159 @@ +LAMMPS (2 Jun 2022) + using 1 OpenMP thread(s) per MPI task +# Demonstrate calculation of SNAP bispectrum descriptors on a grid + +# CORRECTNESS: The two atom positions coincide with two of +# the gridpoints, so c_b[2][1-5] should match c_mygrid[8][4-8]. +# The same is true for compute grid/local c_mygridlocal[8][4-11]. +# Local arrays can not be access directly in the script, +# but they are printed out to file dump.blocal + +variable nrep index 1 +variable a index 3.316 +variable ngrid index 2 + +units metal +atom_modify map hash + +# generate the box and atom positions using a BCC lattice + +variable nx equal ${nrep} +variable nx equal 1 +variable ny equal ${nrep} +variable ny equal 1 +variable nz equal ${nrep} +variable nz equal 1 + +boundary p p p + +lattice custom $a a1 1 0 0 a2 0 1 0 a3 0 0 1 basis 0 0 0 basis 0.5 0.5 0.5 +lattice custom 3.316 a1 1 0 0 a2 0 1 0 a3 0 0 1 basis 0 0 0 basis 0.5 0.5 0.5 +Lattice spacing in x,y,z = 3.316 3.316 3.316 +region box block 0 ${nx} 0 ${ny} 0 ${nz} +region box block 0 1 0 ${ny} 0 ${nz} +region box block 0 1 0 1 0 ${nz} +region box block 0 1 0 1 0 1 +create_box 1 box +Created orthogonal box = (0 0 0) to (3.316 3.316 3.316) + 1 by 1 by 1 MPI processor grid +create_atoms 1 box +Created 2 atoms + using lattice units in orthogonal box = (0 0 0) to (3.316 3.316 3.316) + create_atoms CPU = 0.000 seconds + +mass 1 180.88 + +# define atom compute and grid compute + +group snapgroup type 1 +2 atoms in group snapgroup +variable twojmax equal 2 +variable rcutfac equal 4.67637 +variable rfac0 equal 0.99363 +variable rmin0 equal 0 +variable wj equal 1 +variable radelem equal 0.5 +variable bzero equal 0 +variable quadratic equal 0 +variable switch equal 1 + +variable snap_options string "${rcutfac} ${rfac0} ${twojmax} ${radelem} ${wj} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}" +4.67637 ${rfac0} ${twojmax} ${radelem} ${wj} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch} +4.67637 0.99363 ${twojmax} ${radelem} ${wj} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch} +4.67637 0.99363 2 ${radelem} ${wj} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch} +4.67637 0.99363 2 0.5 ${wj} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch} +4.67637 0.99363 2 0.5 1 rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch} +4.67637 0.99363 2 0.5 1 rmin0 0 quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch} +4.67637 0.99363 2 0.5 1 rmin0 0 quadraticflag 0 bzeroflag ${bzero} switchflag ${switch} +4.67637 0.99363 2 0.5 1 rmin0 0 quadraticflag 0 bzeroflag 0 switchflag ${switch} +4.67637 0.99363 2 0.5 1 rmin0 0 quadraticflag 0 bzeroflag 0 switchflag 1 + +# build zero potential to satisfy compute sna/atom + +pair_style zero ${rcutfac} +pair_style zero 4.67637 +pair_coeff * * + +# define atom and grid computes + +compute b all sna/atom ${snap_options} +compute b all sna/atom 4.67637 0.99363 2 0.5 1 rmin0 0 quadraticflag 0 bzeroflag 0 switchflag 1 +compute mygrid all sna/grid grid ${ngrid} ${ngrid} ${ngrid} ${snap_options} +compute mygrid all sna/grid grid 2 ${ngrid} ${ngrid} ${snap_options} +compute mygrid all sna/grid grid 2 2 ${ngrid} ${snap_options} +compute mygrid all sna/grid grid 2 2 2 ${snap_options} +compute mygrid all sna/grid grid 2 2 2 4.67637 0.99363 2 0.5 1 rmin0 0 quadraticflag 0 bzeroflag 0 switchflag 1 +compute mygridlocal all sna/grid/local grid ${ngrid} ${ngrid} ${ngrid} ${snap_options} +compute mygridlocal all sna/grid/local grid 2 ${ngrid} ${ngrid} ${snap_options} +compute mygridlocal all sna/grid/local grid 2 2 ${ngrid} ${snap_options} +compute mygridlocal all sna/grid/local grid 2 2 2 ${snap_options} +compute mygridlocal all sna/grid/local grid 2 2 2 4.67637 0.99363 2 0.5 1 rmin0 0 quadraticflag 0 bzeroflag 0 switchflag 1 + +# define output + +variable B5atom equal c_b[2][5] +variable B5grid equal c_mygrid[8][8] + +variable rmse_global equal "sqrt( (c_mygrid[8][1] - x[2])^2 + (c_mygrid[8][2] - y[2])^2 + (c_mygrid[8][3] - z[2])^2 + (c_mygrid[8][4] - c_b[2][1])^2 + (c_mygrid[8][5] - c_b[2][2])^2 + (c_mygrid[8][6] - c_b[2][3])^2 + (c_mygrid[8][7] - c_b[2][4])^2 + (c_mygrid[8][8] - c_b[2][5])^2 )" + +thermo_style custom step v_B5atom v_B5grid v_rmse_global + +# this is the only way to view the local grid + +dump 1 all local 1000 dump.blocal c_mygridlocal[*] +dump 2 all custom 1000 dump.batom id x y z c_b[*] + +# run + +run 0 +WARNING: No fixes with time integration, atoms won't move (src/verlet.cpp:60) +Generated 0 of 0 mixed pair_coeff terms from geometric mixing rule +Neighbor list info ... + update every 1 steps, delay 10 steps, check yes + max neighbors/atom: 2000, page size: 100000 + master list distance cutoff = 6.67637 + ghost atom cutoff = 6.67637 + binsize = 3.338185, bins = 1 1 1 + 2 neighbor lists, perpetual/occasional/extra = 1 1 0 + (1) pair zero, perpetual + attributes: half, newton on + pair build: half/bin/atomonly/newton + stencil: half/bin/3d + bin: standard + (2) compute sna/atom, occasional + attributes: full, newton on + pair build: full/bin/atomonly + stencil: full/bin/3d + bin: standard +Per MPI rank memory allocation (min/avg/max) = 7.127 | 7.127 | 7.127 Mbytes + Step v_B5atom v_B5grid v_rmse_global + 0 1.0427295 1.0427295 9.1551336e-16 +Loop time of 1.43e-06 on 1 procs for 0 steps with 2 atoms + +139.9% CPU use with 1 MPI tasks x 1 OpenMP threads + +MPI task timing breakdown: +Section | min time | avg time | max time |%varavg| %total +--------------------------------------------------------------- +Pair | 0 | 0 | 0 | 0.0 | 0.00 +Neigh | 0 | 0 | 0 | 0.0 | 0.00 +Comm | 0 | 0 | 0 | 0.0 | 0.00 +Output | 0 | 0 | 0 | 0.0 | 0.00 +Modify | 0 | 0 | 0 | 0.0 | 0.00 +Other | | 1.43e-06 | | |100.00 + +Nlocal: 2 ave 2 max 2 min +Histogram: 1 0 0 0 0 0 0 0 0 0 +Nghost: 339 ave 339 max 339 min +Histogram: 1 0 0 0 0 0 0 0 0 0 +Neighs: 64 ave 64 max 64 min +Histogram: 1 0 0 0 0 0 0 0 0 0 +FullNghs: 128 ave 128 max 128 min +Histogram: 1 0 0 0 0 0 0 0 0 0 + +Total # of neighbors = 128 +Ave neighs/atom = 64 +Neighbor list builds = 0 +Dangerous builds = 0 + +Total wall time: 0:00:00 diff --git a/examples/snap/log.15Jun22.grid.snap.g++.4 b/examples/snap/log.15Jun22.grid.snap.g++.4 new file mode 100644 index 0000000000..5be17ada7d --- /dev/null +++ b/examples/snap/log.15Jun22.grid.snap.g++.4 @@ -0,0 +1,160 @@ +LAMMPS (2 Jun 2022) + using 1 OpenMP thread(s) per MPI task +# Demonstrate calculation of SNAP bispectrum descriptors on a grid + +# CORRECTNESS: The two atom positions coincide with two of +# the gridpoints, so c_b[2][1-5] should match c_mygrid[8][4-8]. +# The same is true for compute grid/local c_mygridlocal[8][4-11]. +# Local arrays can not be access directly in the script, +# but they are printed out to file dump.blocal + +variable nrep index 1 +variable a index 3.316 +variable ngrid index 2 + +units metal +atom_modify map hash + +# generate the box and atom positions using a BCC lattice + +variable nx equal ${nrep} +variable nx equal 1 +variable ny equal ${nrep} +variable ny equal 1 +variable nz equal ${nrep} +variable nz equal 1 + +boundary p p p + +lattice custom $a a1 1 0 0 a2 0 1 0 a3 0 0 1 basis 0 0 0 basis 0.5 0.5 0.5 +lattice custom 3.316 a1 1 0 0 a2 0 1 0 a3 0 0 1 basis 0 0 0 basis 0.5 0.5 0.5 +Lattice spacing in x,y,z = 3.316 3.316 3.316 +region box block 0 ${nx} 0 ${ny} 0 ${nz} +region box block 0 1 0 ${ny} 0 ${nz} +region box block 0 1 0 1 0 ${nz} +region box block 0 1 0 1 0 1 +create_box 1 box +Created orthogonal box = (0 0 0) to (3.316 3.316 3.316) + 1 by 2 by 2 MPI processor grid +create_atoms 1 box +Created 2 atoms + using lattice units in orthogonal box = (0 0 0) to (3.316 3.316 3.316) + create_atoms CPU = 0.001 seconds + +mass 1 180.88 + +# define atom compute and grid compute + +group snapgroup type 1 +2 atoms in group snapgroup +variable twojmax equal 2 +variable rcutfac equal 4.67637 +variable rfac0 equal 0.99363 +variable rmin0 equal 0 +variable wj equal 1 +variable radelem equal 0.5 +variable bzero equal 0 +variable quadratic equal 0 +variable switch equal 1 + +variable snap_options string "${rcutfac} ${rfac0} ${twojmax} ${radelem} ${wj} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}" +4.67637 ${rfac0} ${twojmax} ${radelem} ${wj} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch} +4.67637 0.99363 ${twojmax} ${radelem} ${wj} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch} +4.67637 0.99363 2 ${radelem} ${wj} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch} +4.67637 0.99363 2 0.5 ${wj} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch} +4.67637 0.99363 2 0.5 1 rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch} +4.67637 0.99363 2 0.5 1 rmin0 0 quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch} +4.67637 0.99363 2 0.5 1 rmin0 0 quadraticflag 0 bzeroflag ${bzero} switchflag ${switch} +4.67637 0.99363 2 0.5 1 rmin0 0 quadraticflag 0 bzeroflag 0 switchflag ${switch} +4.67637 0.99363 2 0.5 1 rmin0 0 quadraticflag 0 bzeroflag 0 switchflag 1 + +# build zero potential to satisfy compute sna/atom + +pair_style zero ${rcutfac} +pair_style zero 4.67637 +pair_coeff * * + +# define atom and grid computes + +compute b all sna/atom ${snap_options} +compute b all sna/atom 4.67637 0.99363 2 0.5 1 rmin0 0 quadraticflag 0 bzeroflag 0 switchflag 1 +compute mygrid all sna/grid grid ${ngrid} ${ngrid} ${ngrid} ${snap_options} +compute mygrid all sna/grid grid 2 ${ngrid} ${ngrid} ${snap_options} +compute mygrid all sna/grid grid 2 2 ${ngrid} ${snap_options} +compute mygrid all sna/grid grid 2 2 2 ${snap_options} +compute mygrid all sna/grid grid 2 2 2 4.67637 0.99363 2 0.5 1 rmin0 0 quadraticflag 0 bzeroflag 0 switchflag 1 +compute mygridlocal all sna/grid/local grid ${ngrid} ${ngrid} ${ngrid} ${snap_options} +compute mygridlocal all sna/grid/local grid 2 ${ngrid} ${ngrid} ${snap_options} +compute mygridlocal all sna/grid/local grid 2 2 ${ngrid} ${snap_options} +compute mygridlocal all sna/grid/local grid 2 2 2 ${snap_options} +compute mygridlocal all sna/grid/local grid 2 2 2 4.67637 0.99363 2 0.5 1 rmin0 0 quadraticflag 0 bzeroflag 0 switchflag 1 + +# define output + +variable B5atom equal c_b[2][5] +variable B5grid equal c_mygrid[8][8] + +variable rmse_global equal "sqrt( (c_mygrid[8][1] - x[2])^2 + (c_mygrid[8][2] - y[2])^2 + (c_mygrid[8][3] - z[2])^2 + (c_mygrid[8][4] - c_b[2][1])^2 + (c_mygrid[8][5] - c_b[2][2])^2 + (c_mygrid[8][6] - c_b[2][3])^2 + (c_mygrid[8][7] - c_b[2][4])^2 + (c_mygrid[8][8] - c_b[2][5])^2 )" + +thermo_style custom step v_B5atom v_B5grid v_rmse_global + +# this is the only way to view the local grid + +dump 1 all local 1000 dump.blocal c_mygridlocal[*] +dump 2 all custom 1000 dump.batom id x y z c_b[*] + +# run + +run 0 +WARNING: No fixes with time integration, atoms won't move (src/verlet.cpp:60) +Generated 0 of 0 mixed pair_coeff terms from geometric mixing rule +Neighbor list info ... + update every 1 steps, delay 10 steps, check yes + max neighbors/atom: 2000, page size: 100000 + master list distance cutoff = 6.67637 + ghost atom cutoff = 6.67637 + binsize = 3.338185, bins = 1 1 1 + 2 neighbor lists, perpetual/occasional/extra = 1 1 0 + (1) pair zero, perpetual + attributes: half, newton on + pair build: half/bin/atomonly/newton + stencil: half/bin/3d + bin: standard + (2) compute sna/atom, occasional + attributes: full, newton on + pair build: full/bin/atomonly + stencil: full/bin/3d + bin: standard +WARNING: Proc sub-domain size < neighbor skin, could lead to lost atoms (src/domain.cpp:970) +Per MPI rank memory allocation (min/avg/max) = 6.123 | 6.631 | 7.139 Mbytes + Step v_B5atom v_B5grid v_rmse_global + 0 1.0427295 1.0427295 1.6316879e-15 +Loop time of 2.57125e-06 on 4 procs for 0 steps with 2 atoms + +107.0% CPU use with 4 MPI tasks x 1 OpenMP threads + +MPI task timing breakdown: +Section | min time | avg time | max time |%varavg| %total +--------------------------------------------------------------- +Pair | 0 | 0 | 0 | 0.0 | 0.00 +Neigh | 0 | 0 | 0 | 0.0 | 0.00 +Comm | 0 | 0 | 0 | 0.0 | 0.00 +Output | 0 | 0 | 0 | 0.0 | 0.00 +Modify | 0 | 0 | 0 | 0.0 | 0.00 +Other | | 2.571e-06 | | |100.00 + +Nlocal: 0.5 ave 1 max 0 min +Histogram: 2 0 0 0 0 0 0 0 0 2 +Nghost: 274.5 ave 275 max 274 min +Histogram: 2 0 0 0 0 0 0 0 0 2 +Neighs: 16 ave 40 max 0 min +Histogram: 2 0 0 0 0 0 1 0 0 1 +FullNghs: 32 ave 64 max 0 min +Histogram: 2 0 0 0 0 0 0 0 0 2 + +Total # of neighbors = 128 +Ave neighs/atom = 64 +Neighbor list builds = 0 +Dangerous builds = 0 + +Total wall time: 0:00:00 diff --git a/examples/snap/log.15Jun22.grid.tri.g++.1 b/examples/snap/log.15Jun22.grid.tri.g++.1 new file mode 100644 index 0000000000..e26315235b --- /dev/null +++ b/examples/snap/log.15Jun22.grid.tri.g++.1 @@ -0,0 +1,192 @@ +LAMMPS (2 Jun 2022) + using 1 OpenMP thread(s) per MPI task +# Demonstrate calculation of SNAP bispectrum +# descriptors on a grid for triclinic cell + +# This triclinic cell has 6 times the volume of the single +# unit cell used by in.grid +# and contains 12 atoms. It is a 3x2x1 supercell +# with each unit cell containing 2 atoms and the +# reduced lattice vectors are [1 0 0], [1 1 0], and [1 1 1]. +# The grid is listed in x-fastest order + +# CORRECTNESS: The atom positions coincide with certain +# gridpoints, so c_b[1][1-5] should match c_mygrid[1][4-8] +# and c_b[7][1-5] should match c_mygrid[13][4-8]. +# Local arrays can not be access directly in the script, +# but they are printed out to file dump.blocal.tri + +# Initialize simulation + +variable nrep index 1 +variable a index 3.316 +variable ngrid index 2 + +variable nrepx equal 3*${nrep} +variable nrepx equal 3*1 +variable nrepy equal 2*${nrep} +variable nrepy equal 2*1 +variable nrepz equal 1*${nrep} +variable nrepz equal 1*1 + +variable ngridx equal 3*${ngrid} +variable ngridx equal 3*2 +variable ngridy equal 2*${ngrid} +variable ngridy equal 2*2 +variable ngridz equal 1*${ngrid} +variable ngridz equal 1*2 + +units metal +atom_modify map hash sort 0 0 + +# generate the box and atom positions using a BCC lattice + +variable nx equal ${nrepx} +variable nx equal 3 +variable ny equal ${nrepy} +variable ny equal 2 +variable nz equal ${nrepz} +variable nz equal 1 + +boundary p p p + +lattice custom $a a1 1 0 0 a2 1 1 0 a3 1 1 1 basis 0 0 0 basis 0.0 0.0 0.5 spacing 1 1 1 +lattice custom 3.316 a1 1 0 0 a2 1 1 0 a3 1 1 1 basis 0 0 0 basis 0.0 0.0 0.5 spacing 1 1 1 +Lattice spacing in x,y,z = 3.316 3.316 3.316 + +box tilt large +region box prism 0 ${nx} 0 ${ny} 0 ${nz} ${ny} ${nz} ${nz} +region box prism 0 3 0 ${ny} 0 ${nz} ${ny} ${nz} ${nz} +region box prism 0 3 0 2 0 ${nz} ${ny} ${nz} ${nz} +region box prism 0 3 0 2 0 1 ${ny} ${nz} ${nz} +region box prism 0 3 0 2 0 1 2 ${nz} ${nz} +region box prism 0 3 0 2 0 1 2 1 ${nz} +region box prism 0 3 0 2 0 1 2 1 1 +create_box 1 box +Created triclinic box = (0 0 0) to (9.948 6.632 3.316) with tilt (6.632 3.316 3.316) +WARNING: Triclinic box skew is large (src/domain.cpp:224) + 1 by 1 by 1 MPI processor grid +create_atoms 1 box +Created 12 atoms + using lattice units in triclinic box = (0 0 0) to (9.948 6.632 3.316) with tilt (6.632 3.316 3.316) + create_atoms CPU = 0.000 seconds + +mass 1 180.88 + +# define atom compute and grid compute + +group snapgroup type 1 +12 atoms in group snapgroup +variable twojmax equal 2 +variable rcutfac equal 4.67637 +variable rfac0 equal 0.99363 +variable rmin0 equal 0 +variable wj equal 1 +variable radelem equal 0.5 +variable bzero equal 0 +variable quadratic equal 0 +variable switch equal 1 + +variable snap_options string "${rcutfac} ${rfac0} ${twojmax} ${radelem} ${wj} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}" +4.67637 ${rfac0} ${twojmax} ${radelem} ${wj} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch} +4.67637 0.99363 ${twojmax} ${radelem} ${wj} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch} +4.67637 0.99363 2 ${radelem} ${wj} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch} +4.67637 0.99363 2 0.5 ${wj} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch} +4.67637 0.99363 2 0.5 1 rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch} +4.67637 0.99363 2 0.5 1 rmin0 0 quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch} +4.67637 0.99363 2 0.5 1 rmin0 0 quadraticflag 0 bzeroflag ${bzero} switchflag ${switch} +4.67637 0.99363 2 0.5 1 rmin0 0 quadraticflag 0 bzeroflag 0 switchflag ${switch} +4.67637 0.99363 2 0.5 1 rmin0 0 quadraticflag 0 bzeroflag 0 switchflag 1 + +# build zero potential to satisfy compute sna/atom + +pair_style zero ${rcutfac} +pair_style zero 4.67637 +pair_coeff * * + +# define atom and grid computes + +compute b all sna/atom ${snap_options} +compute b all sna/atom 4.67637 0.99363 2 0.5 1 rmin0 0 quadraticflag 0 bzeroflag 0 switchflag 1 +compute mygrid all sna/grid grid ${ngridx} ${ngridy} ${ngridz} ${snap_options} +compute mygrid all sna/grid grid 6 ${ngridy} ${ngridz} ${snap_options} +compute mygrid all sna/grid grid 6 4 ${ngridz} ${snap_options} +compute mygrid all sna/grid grid 6 4 2 ${snap_options} +compute mygrid all sna/grid grid 6 4 2 4.67637 0.99363 2 0.5 1 rmin0 0 quadraticflag 0 bzeroflag 0 switchflag 1 +compute mygridlocal all sna/grid/local grid ${ngridx} ${ngridy} ${ngridz} ${snap_options} +compute mygridlocal all sna/grid/local grid 6 ${ngridy} ${ngridz} ${snap_options} +compute mygridlocal all sna/grid/local grid 6 4 ${ngridz} ${snap_options} +compute mygridlocal all sna/grid/local grid 6 4 2 ${snap_options} +compute mygridlocal all sna/grid/local grid 6 4 2 4.67637 0.99363 2 0.5 1 rmin0 0 quadraticflag 0 bzeroflag 0 switchflag 1 + +# define output + +variable B5atom equal c_b[7][5] +variable B5grid equal c_mygrid[13][8] + +# do not compare x,y,z because assignment of ids +# to atoms is not unnique for different processor grids + +variable rmse_global equal "sqrt( (c_mygrid[13][4] - c_b[7][1])^2 + (c_mygrid[13][5] - c_b[7][2])^2 + (c_mygrid[13][6] - c_b[7][3])^2 + (c_mygrid[13][7] - c_b[7][4])^2 + (c_mygrid[13][8] - c_b[7][5])^2 )" + +thermo_style custom step v_B5atom v_B5grid v_rmse_global + +# this is the only way to view the local grid + +dump 1 all local 1000 dump.blocal.tri c_mygridlocal[*] +dump 2 all custom 1000 dump.batom.tri id x y z c_b[*] + +# run + +run 0 +WARNING: No fixes with time integration, atoms won't move (src/verlet.cpp:60) +Generated 0 of 0 mixed pair_coeff terms from geometric mixing rule +Neighbor list info ... + update every 1 steps, delay 10 steps, check yes + max neighbors/atom: 2000, page size: 100000 + master list distance cutoff = 6.67637 + ghost atom cutoff = 6.67637 + binsize = 3.338185, bins = 6 3 1 + 2 neighbor lists, perpetual/occasional/extra = 1 1 0 + (1) pair zero, perpetual + attributes: half, newton on + pair build: half/bin/newton/tri + stencil: half/bin/3d/tri + bin: standard + (2) compute sna/atom, occasional + attributes: full, newton on + pair build: full/bin/atomonly + stencil: full/bin/3d + bin: standard +Per MPI rank memory allocation (min/avg/max) = 7.183 | 7.183 | 7.183 Mbytes + Step v_B5atom v_B5grid v_rmse_global + 0 1.0427295 1.0427295 7.2262471e-14 +Loop time of 1.414e-06 on 1 procs for 0 steps with 12 atoms + +70.7% CPU use with 1 MPI tasks x 1 OpenMP threads + +MPI task timing breakdown: +Section | min time | avg time | max time |%varavg| %total +--------------------------------------------------------------- +Pair | 0 | 0 | 0 | 0.0 | 0.00 +Neigh | 0 | 0 | 0 | 0.0 | 0.00 +Comm | 0 | 0 | 0 | 0.0 | 0.00 +Output | 0 | 0 | 0 | 0.0 | 0.00 +Modify | 0 | 0 | 0 | 0.0 | 0.00 +Other | | 1.414e-06 | | |100.00 + +Nlocal: 12 ave 12 max 12 min +Histogram: 1 0 0 0 0 0 0 0 0 0 +Nghost: 604 ave 604 max 604 min +Histogram: 1 0 0 0 0 0 0 0 0 0 +Neighs: 384 ave 384 max 384 min +Histogram: 1 0 0 0 0 0 0 0 0 0 +FullNghs: 768 ave 768 max 768 min +Histogram: 1 0 0 0 0 0 0 0 0 0 + +Total # of neighbors = 768 +Ave neighs/atom = 64 +Neighbor list builds = 0 +Dangerous builds = 0 + +Total wall time: 0:00:00 diff --git a/examples/snap/log.15Jun22.grid.tri.g++.4 b/examples/snap/log.15Jun22.grid.tri.g++.4 new file mode 100644 index 0000000000..cee3ce7f12 --- /dev/null +++ b/examples/snap/log.15Jun22.grid.tri.g++.4 @@ -0,0 +1,192 @@ +LAMMPS (2 Jun 2022) + using 1 OpenMP thread(s) per MPI task +# Demonstrate calculation of SNAP bispectrum +# descriptors on a grid for triclinic cell + +# This triclinic cell has 6 times the volume of the single +# unit cell used by in.grid +# and contains 12 atoms. It is a 3x2x1 supercell +# with each unit cell containing 2 atoms and the +# reduced lattice vectors are [1 0 0], [1 1 0], and [1 1 1]. +# The grid is listed in x-fastest order + +# CORRECTNESS: The atom positions coincide with certain +# gridpoints, so c_b[1][1-5] should match c_mygrid[1][4-8] +# and c_b[7][1-5] should match c_mygrid[13][4-8]. +# Local arrays can not be access directly in the script, +# but they are printed out to file dump.blocal.tri + +# Initialize simulation + +variable nrep index 1 +variable a index 3.316 +variable ngrid index 2 + +variable nrepx equal 3*${nrep} +variable nrepx equal 3*1 +variable nrepy equal 2*${nrep} +variable nrepy equal 2*1 +variable nrepz equal 1*${nrep} +variable nrepz equal 1*1 + +variable ngridx equal 3*${ngrid} +variable ngridx equal 3*2 +variable ngridy equal 2*${ngrid} +variable ngridy equal 2*2 +variable ngridz equal 1*${ngrid} +variable ngridz equal 1*2 + +units metal +atom_modify map hash sort 0 0 + +# generate the box and atom positions using a BCC lattice + +variable nx equal ${nrepx} +variable nx equal 3 +variable ny equal ${nrepy} +variable ny equal 2 +variable nz equal ${nrepz} +variable nz equal 1 + +boundary p p p + +lattice custom $a a1 1 0 0 a2 1 1 0 a3 1 1 1 basis 0 0 0 basis 0.0 0.0 0.5 spacing 1 1 1 +lattice custom 3.316 a1 1 0 0 a2 1 1 0 a3 1 1 1 basis 0 0 0 basis 0.0 0.0 0.5 spacing 1 1 1 +Lattice spacing in x,y,z = 3.316 3.316 3.316 + +box tilt large +region box prism 0 ${nx} 0 ${ny} 0 ${nz} ${ny} ${nz} ${nz} +region box prism 0 3 0 ${ny} 0 ${nz} ${ny} ${nz} ${nz} +region box prism 0 3 0 2 0 ${nz} ${ny} ${nz} ${nz} +region box prism 0 3 0 2 0 1 ${ny} ${nz} ${nz} +region box prism 0 3 0 2 0 1 2 ${nz} ${nz} +region box prism 0 3 0 2 0 1 2 1 ${nz} +region box prism 0 3 0 2 0 1 2 1 1 +create_box 1 box +Created triclinic box = (0 0 0) to (9.948 6.632 3.316) with tilt (6.632 3.316 3.316) +WARNING: Triclinic box skew is large (src/domain.cpp:224) + 2 by 2 by 1 MPI processor grid +create_atoms 1 box +Created 12 atoms + using lattice units in triclinic box = (0 0 0) to (9.948 6.632 3.316) with tilt (6.632 3.316 3.316) + create_atoms CPU = 0.000 seconds + +mass 1 180.88 + +# define atom compute and grid compute + +group snapgroup type 1 +12 atoms in group snapgroup +variable twojmax equal 2 +variable rcutfac equal 4.67637 +variable rfac0 equal 0.99363 +variable rmin0 equal 0 +variable wj equal 1 +variable radelem equal 0.5 +variable bzero equal 0 +variable quadratic equal 0 +variable switch equal 1 + +variable snap_options string "${rcutfac} ${rfac0} ${twojmax} ${radelem} ${wj} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}" +4.67637 ${rfac0} ${twojmax} ${radelem} ${wj} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch} +4.67637 0.99363 ${twojmax} ${radelem} ${wj} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch} +4.67637 0.99363 2 ${radelem} ${wj} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch} +4.67637 0.99363 2 0.5 ${wj} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch} +4.67637 0.99363 2 0.5 1 rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch} +4.67637 0.99363 2 0.5 1 rmin0 0 quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch} +4.67637 0.99363 2 0.5 1 rmin0 0 quadraticflag 0 bzeroflag ${bzero} switchflag ${switch} +4.67637 0.99363 2 0.5 1 rmin0 0 quadraticflag 0 bzeroflag 0 switchflag ${switch} +4.67637 0.99363 2 0.5 1 rmin0 0 quadraticflag 0 bzeroflag 0 switchflag 1 + +# build zero potential to satisfy compute sna/atom + +pair_style zero ${rcutfac} +pair_style zero 4.67637 +pair_coeff * * + +# define atom and grid computes + +compute b all sna/atom ${snap_options} +compute b all sna/atom 4.67637 0.99363 2 0.5 1 rmin0 0 quadraticflag 0 bzeroflag 0 switchflag 1 +compute mygrid all sna/grid grid ${ngridx} ${ngridy} ${ngridz} ${snap_options} +compute mygrid all sna/grid grid 6 ${ngridy} ${ngridz} ${snap_options} +compute mygrid all sna/grid grid 6 4 ${ngridz} ${snap_options} +compute mygrid all sna/grid grid 6 4 2 ${snap_options} +compute mygrid all sna/grid grid 6 4 2 4.67637 0.99363 2 0.5 1 rmin0 0 quadraticflag 0 bzeroflag 0 switchflag 1 +compute mygridlocal all sna/grid/local grid ${ngridx} ${ngridy} ${ngridz} ${snap_options} +compute mygridlocal all sna/grid/local grid 6 ${ngridy} ${ngridz} ${snap_options} +compute mygridlocal all sna/grid/local grid 6 4 ${ngridz} ${snap_options} +compute mygridlocal all sna/grid/local grid 6 4 2 ${snap_options} +compute mygridlocal all sna/grid/local grid 6 4 2 4.67637 0.99363 2 0.5 1 rmin0 0 quadraticflag 0 bzeroflag 0 switchflag 1 + +# define output + +variable B5atom equal c_b[7][5] +variable B5grid equal c_mygrid[13][8] + +# do not compare x,y,z because assignment of ids +# to atoms is not unnique for different processor grids + +variable rmse_global equal "sqrt( (c_mygrid[13][4] - c_b[7][1])^2 + (c_mygrid[13][5] - c_b[7][2])^2 + (c_mygrid[13][6] - c_b[7][3])^2 + (c_mygrid[13][7] - c_b[7][4])^2 + (c_mygrid[13][8] - c_b[7][5])^2 )" + +thermo_style custom step v_B5atom v_B5grid v_rmse_global + +# this is the only way to view the local grid + +dump 1 all local 1000 dump.blocal.tri c_mygridlocal[*] +dump 2 all custom 1000 dump.batom.tri id x y z c_b[*] + +# run + +run 0 +WARNING: No fixes with time integration, atoms won't move (src/verlet.cpp:60) +Generated 0 of 0 mixed pair_coeff terms from geometric mixing rule +Neighbor list info ... + update every 1 steps, delay 10 steps, check yes + max neighbors/atom: 2000, page size: 100000 + master list distance cutoff = 6.67637 + ghost atom cutoff = 6.67637 + binsize = 3.338185, bins = 6 3 1 + 2 neighbor lists, perpetual/occasional/extra = 1 1 0 + (1) pair zero, perpetual + attributes: half, newton on + pair build: half/bin/newton/tri + stencil: half/bin/3d/tri + bin: standard + (2) compute sna/atom, occasional + attributes: full, newton on + pair build: full/bin/atomonly + stencil: full/bin/3d + bin: standard +Per MPI rank memory allocation (min/avg/max) = 7.15 | 7.15 | 7.15 Mbytes + Step v_B5atom v_B5grid v_rmse_global + 0 1.0427295 1.0427295 1.9367585e-14 +Loop time of 2.65825e-06 on 4 procs for 0 steps with 12 atoms + +84.6% CPU use with 4 MPI tasks x 1 OpenMP threads + +MPI task timing breakdown: +Section | min time | avg time | max time |%varavg| %total +--------------------------------------------------------------- +Pair | 0 | 0 | 0 | 0.0 | 0.00 +Neigh | 0 | 0 | 0 | 0.0 | 0.00 +Comm | 0 | 0 | 0 | 0.0 | 0.00 +Output | 0 | 0 | 0 | 0.0 | 0.00 +Modify | 0 | 0 | 0 | 0.0 | 0.00 +Other | | 2.658e-06 | | |100.00 + +Nlocal: 3 ave 4 max 2 min +Histogram: 2 0 0 0 0 0 0 0 0 2 +Nghost: 459 ave 460 max 458 min +Histogram: 2 0 0 0 0 0 0 0 0 2 +Neighs: 96 ave 128 max 64 min +Histogram: 2 0 0 0 0 0 0 0 0 2 +FullNghs: 192 ave 256 max 128 min +Histogram: 2 0 0 0 0 0 0 0 0 2 + +Total # of neighbors = 768 +Ave neighs/atom = 64 +Neighbor list builds = 0 +Dangerous builds = 0 + +Total wall time: 0:00:00 diff --git a/src/.gitignore b/src/.gitignore index 220cd621db..4c1cf49b3b 100644 --- a/src/.gitignore +++ b/src/.gitignore @@ -173,12 +173,20 @@ /pair_tdpd.cpp /pair_tdpd.h +/compute_grid.cpp +/compute_grid.h +/compute_grid_local.cpp +/compute_grid_local.h /compute_sna_atom.cpp /compute_sna_atom.h /compute_snad_atom.cpp /compute_snad_atom.h /compute_snav_atom.cpp /compute_snav_atom.h +/compute_sna_grid.cpp +/compute_sna_grid.h +/compute_sna_grid_local.cpp +/compute_sna_grid_local.h /compute_snap.cpp /compute_snap.h /openmp_snap.h diff --git a/src/ML-SNAP/compute_grid.cpp b/src/ML-SNAP/compute_grid.cpp new file mode 100644 index 0000000000..7a5ffe1e24 --- /dev/null +++ b/src/ML-SNAP/compute_grid.cpp @@ -0,0 +1,242 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/ Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#include "compute_grid.h" + +#include "atom.h" +#include "comm.h" +#include "domain.h" +#include "error.h" +#include "force.h" +#include "memory.h" +#include "modify.h" +#include "update.h" + +#include + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +ComputeGrid::ComputeGrid(LAMMPS *lmp, int narg, char **arg) : + Compute(lmp, narg, arg), grid(nullptr), gridall(nullptr), gridlocal(nullptr) +{ + if (narg < 6) error->all(FLERR, "Illegal compute grid command"); + + array_flag = 1; + size_array_cols = 0; + size_array_rows = 0; + extarray = 0; + + int iarg0 = 3; + int iarg = iarg0; + if (strcmp(arg[iarg], "grid") == 0) { + if (iarg + 4 > narg) error->all(FLERR, "Illegal compute grid command"); + nx = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); + ny = utils::inumeric(FLERR, arg[iarg + 2], false, lmp); + nz = utils::inumeric(FLERR, arg[iarg + 3], false, lmp); + if (nx <= 0 || ny <= 0 || nz <= 0) error->all(FLERR, "All grid dimensions must be positive"); + iarg += 4; + } else + error->all(FLERR, "Illegal compute grid command"); + + nargbase = iarg - iarg0; + + size_array_rows = nx * ny * nz; + size_array_cols_base = 3; + gridlocal_allocated = 0; +} + +/* ---------------------------------------------------------------------- */ + +ComputeGrid::~ComputeGrid() +{ + deallocate(); +} + +/* ---------------------------------------------------------------------- */ + +void ComputeGrid::setup() +{ + deallocate(); + set_grid_global(); + set_grid_local(); + allocate(); +} + +/* ---------------------------------------------------------------------- + convert global array index to box coords +------------------------------------------------------------------------- */ + +void ComputeGrid::grid2x(int igrid, double *x) +{ + int iz = igrid / (nx * ny); + igrid -= iz * (nx * ny); + int iy = igrid / nx; + igrid -= iy * nx; + int ix = igrid; + + x[0] = ix * delx; + x[1] = iy * dely; + x[2] = iz * delz; + + if (triclinic) domain->lamda2x(x, x); +} + +/* ---------------------------------------------------------------------- + copy coords to global array +------------------------------------------------------------------------- */ + +void ComputeGrid::assign_coords_all() +{ + double x[3]; + for (int igrid = 0; igrid < size_array_rows; igrid++) { + grid2x(igrid, x); + gridall[igrid][0] = x[0]; + gridall[igrid][1] = x[1]; + gridall[igrid][2] = x[2]; + } +} + +/* ---------------------------------------------------------------------- + create arrays +------------------------------------------------------------------------- */ + +void ComputeGrid::allocate() +{ + // allocate arrays + + memory->create(grid, size_array_rows, size_array_cols, "grid:grid"); + memory->create(gridall, size_array_rows, size_array_cols, "grid:gridall"); + if (nxlo <= nxhi && nylo <= nyhi && nzlo <= nzhi) { + gridlocal_allocated = 1; + memory->create4d_offset(gridlocal, size_array_cols, nzlo, nzhi, nylo, nyhi, nxlo, nxhi, + "grid:gridlocal"); + } + array = gridall; +} + +/* ---------------------------------------------------------------------- + free arrays +------------------------------------------------------------------------- */ + +void ComputeGrid::deallocate() +{ + memory->destroy(grid); + memory->destroy(gridall); + if (gridlocal_allocated) { + gridlocal_allocated = 0; + memory->destroy4d_offset(gridlocal, nzlo, nylo, nxlo); + } + array = nullptr; +} + +/* ---------------------------------------------------------------------- + set global grid +------------------------------------------------------------------------- */ + +void ComputeGrid::set_grid_global() +{ + // calculate grid layout + + triclinic = domain->triclinic; + + if (triclinic == 0) { + prd = domain->prd; + boxlo = domain->boxlo; + sublo = domain->sublo; + subhi = domain->subhi; + } else { + prd = domain->prd_lamda; + boxlo = domain->boxlo_lamda; + sublo = domain->sublo_lamda; + subhi = domain->subhi_lamda; + } + + double xprd = prd[0]; + double yprd = prd[1]; + double zprd = prd[2]; + + delxinv = nx / xprd; + delyinv = ny / yprd; + delzinv = nz / zprd; + + delx = 1.0 / delxinv; + dely = 1.0 / delyinv; + delz = 1.0 / delzinv; +} + +/* ---------------------------------------------------------------------- + set local subset of grid that I own + n xyz lo/hi = 3d brick that I own (inclusive) +------------------------------------------------------------------------- */ + +void ComputeGrid::set_grid_local() +{ + // nx,ny,nz = extent of global grid + // indices into the global grid range from 0 to N-1 in each dim + // if grid point is inside my sub-domain I own it, + // this includes sub-domain lo boundary but excludes hi boundary + // ixyz lo/hi = inclusive lo/hi bounds of global grid sub-brick I own + // if proc owns no grid cells in a dim, then ilo > ihi + // if 2 procs share a boundary a grid point is exactly on, + // the 2 equality if tests insure a consistent decision + // as to which proc owns it + + double xfraclo, xfrachi, yfraclo, yfrachi, zfraclo, zfrachi; + + if (comm->layout != Comm::LAYOUT_TILED) { + xfraclo = comm->xsplit[comm->myloc[0]]; + xfrachi = comm->xsplit[comm->myloc[0] + 1]; + yfraclo = comm->ysplit[comm->myloc[1]]; + yfrachi = comm->ysplit[comm->myloc[1] + 1]; + zfraclo = comm->zsplit[comm->myloc[2]]; + zfrachi = comm->zsplit[comm->myloc[2] + 1]; + } else { + xfraclo = comm->mysplit[0][0]; + xfrachi = comm->mysplit[0][1]; + yfraclo = comm->mysplit[1][0]; + yfrachi = comm->mysplit[1][1]; + zfraclo = comm->mysplit[2][0]; + zfrachi = comm->mysplit[2][1]; + } + + nxlo = static_cast(xfraclo * nx); + if (1.0 * nxlo != xfraclo * nx) nxlo++; + nxhi = static_cast(xfrachi * nx); + if (1.0 * nxhi == xfrachi * nx) nxhi--; + + nylo = static_cast(yfraclo * ny); + if (1.0 * nylo != yfraclo * ny) nylo++; + nyhi = static_cast(yfrachi * ny); + if (1.0 * nyhi == yfrachi * ny) nyhi--; + + nzlo = static_cast(zfraclo * nz); + if (1.0 * nzlo != zfraclo * nz) nzlo++; + nzhi = static_cast(zfrachi * nz); + if (1.0 * nzhi == zfrachi * nz) nzhi--; + + ngridlocal = (nxhi - nxlo + 1) * (nyhi - nylo + 1) * (nzhi - nzlo + 1); +} + +/* ---------------------------------------------------------------------- + memory usage of local data +------------------------------------------------------------------------- */ + +double ComputeGrid::memory_usage() +{ + double nbytes = size_array_rows * size_array_cols * sizeof(double); // grid + nbytes += size_array_rows * size_array_cols * sizeof(double); // gridall + nbytes += size_array_cols * ngridlocal * sizeof(double); // gridlocal + return nbytes; +} diff --git a/src/ML-SNAP/compute_grid.h b/src/ML-SNAP/compute_grid.h new file mode 100644 index 0000000000..84b3bc4e98 --- /dev/null +++ b/src/ML-SNAP/compute_grid.h @@ -0,0 +1,58 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/ Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifndef LMP_COMPUTE_GRID_H +#define LMP_COMPUTE_GRID_H + +#include "compute.h" + +namespace LAMMPS_NS { + +class ComputeGrid : public Compute { + public: + ComputeGrid(class LAMMPS *, int, char **); + ~ComputeGrid() override; + void setup() override; + void compute_array() override = 0; + + double memory_usage() override; + + protected: + int nx, ny, nz; // global grid dimensions + int nxlo, nxhi, nylo, nyhi, nzlo, nzhi; // local grid bounds, inclusive + int ngridlocal; // number of local grid points + int nvalues; // number of values per grid point + double **grid; // global grid + double **gridall; // global grid summed over procs + double ****gridlocal; // local grid + int triclinic; // triclinic flag + double *boxlo, *prd; // box info (units real/ortho or reduced/tri) + double *sublo, *subhi; // subdomain info (units real/ortho or reduced/tri) + double delxinv, delyinv, delzinv; // inverse grid spacing + double delx, dely, delz; // grid spacing + int nargbase; // number of base class args + double cutmax; // largest cutoff distance + int size_array_cols_base; // number of columns used for coords, etc. + int gridlocal_allocated; // shows if gridlocal allocated + + void allocate(); // create arrays + void deallocate(); // free arrays + void grid2x(int, double *); // convert grid point to coord + void assign_coords_all(); // assign coords for global grid + void set_grid_global(); // set global grid + void set_grid_local(); // set bounds for local grid +}; + +} // namespace LAMMPS_NS + +#endif diff --git a/src/ML-SNAP/compute_grid_local.cpp b/src/ML-SNAP/compute_grid_local.cpp new file mode 100644 index 0000000000..48714962ac --- /dev/null +++ b/src/ML-SNAP/compute_grid_local.cpp @@ -0,0 +1,270 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/ Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#include "compute_grid_local.h" + +#include "atom.h" +#include "comm.h" +#include "domain.h" +#include "error.h" +#include "force.h" +#include "memory.h" +#include "modify.h" +#include "update.h" + +#include + +// For the subdomain test below; grid-points and subdomain boundaries +// sometimes differ by minimal amounts (in the order of 2e-17). +static constexpr double EPSILON = 1.0e-10; + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +ComputeGridLocal::ComputeGridLocal(LAMMPS *lmp, int narg, char **arg) : + Compute(lmp, narg, arg), alocal(nullptr) +{ + if (narg < 6) error->all(FLERR, "Illegal compute grid/local command"); + + local_flag = 1; + size_local_cols = 0; + size_local_rows = 0; + extarray = 0; + + int iarg0 = 3; + int iarg = iarg0; + if (strcmp(arg[iarg], "grid") == 0) { + if (iarg + 4 > narg) error->all(FLERR, "Illegal compute grid/local command"); + nx = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); + ny = utils::inumeric(FLERR, arg[iarg + 2], false, lmp); + nz = utils::inumeric(FLERR, arg[iarg + 3], false, lmp); + if (nx <= 0 || ny <= 0 || nz <= 0) + error->all(FLERR, "All grid/local dimensions must be positive"); + iarg += 4; + } else + error->all(FLERR, "Illegal compute grid/local command"); + + nargbase = iarg - iarg0; + + size_local_cols_base = 6; + gridlocal_allocated = 0; +} + +/* ---------------------------------------------------------------------- */ + +ComputeGridLocal::~ComputeGridLocal() +{ + deallocate(); +} + +/* ---------------------------------------------------------------------- */ + +void ComputeGridLocal::setup() +{ + deallocate(); + set_grid_global(); + set_grid_local(); + allocate(); + assign_coords(); +} + +/* ---------------------------------------------------------------------- + convert global array indexes to box coords +------------------------------------------------------------------------- */ + +void ComputeGridLocal::grid2x(int ix, int iy, int iz, double *x) +{ + x[0] = ix * delx; + x[1] = iy * dely; + x[2] = iz * delz; + + if (triclinic) domain->lamda2x(x, x); +} + +/* ---------------------------------------------------------------------- + convert global array indexes to lamda coords; for orthorombic + cells defaults to grid2x. +------------------------------------------------------------------------- */ + +void ComputeGridLocal::grid2lamda(int ix, int iy, int iz, double *x) +{ + x[0] = ix * delx; + x[1] = iy * dely; + x[2] = iz * delz; +} + +/* ---------------------------------------------------------------------- + create arrays +------------------------------------------------------------------------- */ + +void ComputeGridLocal::allocate() +{ + if (nxlo <= nxhi && nylo <= nyhi && nzlo <= nzhi) { + gridlocal_allocated = 1; + memory->create(alocal, size_local_rows, size_local_cols, "compute/grid/local:alocal"); + array_local = alocal; + } +} + +/* ---------------------------------------------------------------------- + free arrays +------------------------------------------------------------------------- */ + +void ComputeGridLocal::deallocate() +{ + if (gridlocal_allocated) { + gridlocal_allocated = 0; + memory->destroy(alocal); + } + array_local = nullptr; +} + +/* ---------------------------------------------------------------------- + set global grid +------------------------------------------------------------------------- */ + +void ComputeGridLocal::set_grid_global() +{ + // calculate grid layout + + triclinic = domain->triclinic; + + if (triclinic == 0) { + prd = domain->prd; + boxlo = domain->boxlo; + sublo = domain->sublo; + subhi = domain->subhi; + } else { + prd = domain->prd_lamda; + boxlo = domain->boxlo_lamda; + sublo = domain->sublo_lamda; + subhi = domain->subhi_lamda; + } + + double xprd = prd[0]; + double yprd = prd[1]; + double zprd = prd[2]; + + delxinv = nx / xprd; + delyinv = ny / yprd; + delzinv = nz / zprd; + + delx = 1.0 / delxinv; + dely = 1.0 / delyinv; + delz = 1.0 / delzinv; +} + +/* ---------------------------------------------------------------------- + set local subset of grid that I own + n xyz lo/hi = 3d brick that I own (inclusive) +------------------------------------------------------------------------- */ + +void ComputeGridLocal::set_grid_local() +{ + // nx,ny,nz = extent of global grid + // indices into the global grid range from 0 to N-1 in each dim + // if grid point is inside my sub-domain I own it, + // this includes sub-domain lo boundary but excludes hi boundary + // ixyz lo/hi = inclusive lo/hi bounds of global grid sub-brick I own + // if proc owns no grid cells in a dim, then ilo > ihi + // if 2 procs share a boundary a grid point is exactly on, + // the 2 equality if tests insure a consistent decision + // as to which proc owns it + + double xfraclo, xfrachi, yfraclo, yfrachi, zfraclo, zfrachi; + + if (comm->layout != Comm::LAYOUT_TILED) { + xfraclo = comm->xsplit[comm->myloc[0]]; + xfrachi = comm->xsplit[comm->myloc[0] + 1]; + yfraclo = comm->ysplit[comm->myloc[1]]; + yfrachi = comm->ysplit[comm->myloc[1] + 1]; + zfraclo = comm->zsplit[comm->myloc[2]]; + zfrachi = comm->zsplit[comm->myloc[2] + 1]; + } else { + xfraclo = comm->mysplit[0][0]; + xfrachi = comm->mysplit[0][1]; + yfraclo = comm->mysplit[1][0]; + yfrachi = comm->mysplit[1][1]; + zfraclo = comm->mysplit[2][0]; + zfrachi = comm->mysplit[2][1]; + } + + nxlo = static_cast(xfraclo * nx); + if (1.0 * nxlo != xfraclo * nx) nxlo++; + nxhi = static_cast(xfrachi * nx); + if (1.0 * nxhi == xfrachi * nx) nxhi--; + + nylo = static_cast(yfraclo * ny); + if (1.0 * nylo != yfraclo * ny) nylo++; + nyhi = static_cast(yfrachi * ny); + if (1.0 * nyhi == yfrachi * ny) nyhi--; + + nzlo = static_cast(zfraclo * nz); + if (1.0 * nzlo != zfraclo * nz) nzlo++; + nzhi = static_cast(zfrachi * nz); + if (1.0 * nzhi == zfrachi * nz) nzhi--; + + size_local_rows = (nxhi - nxlo + 1) * (nyhi - nylo + 1) * (nzhi - nzlo + 1); +} + +/* ---------------------------------------------------------------------- + copy coords to local array +------------------------------------------------------------------------- */ + +void ComputeGridLocal::assign_coords() +{ + int igrid = 0; + for (int iz = nzlo; iz <= nzhi; iz++) + for (int iy = nylo; iy <= nyhi; iy++) + for (int ix = nxlo; ix <= nxhi; ix++) { + alocal[igrid][0] = ix; + alocal[igrid][1] = iy; + alocal[igrid][2] = iz; + double xgrid[3]; + + // for triclinic: create gridpoint in lamda coordinates and transform after check. + // for orthorombic: create gridpoint in box coordinates. + + if (triclinic) + grid2lamda(ix, iy, iz, xgrid); + else + grid2x(ix, iy, iz, xgrid); + + // ensure gridpoint is not strictly outside subdomain + + if ((sublo[0] - xgrid[0]) > EPSILON || (xgrid[0] - subhi[0]) > EPSILON || + (sublo[1] - xgrid[1]) > EPSILON || (xgrid[1] - subhi[1]) > EPSILON || + (sublo[2] - xgrid[2]) > EPSILON || (xgrid[2] - subhi[2]) > EPSILON) + error->one(FLERR, "Invalid gridpoint position in compute grid/local"); + + // convert lamda to x, y, z, after sudomain check + + if (triclinic) domain->lamda2x(xgrid, xgrid); + + alocal[igrid][3] = xgrid[0]; + alocal[igrid][4] = xgrid[1]; + alocal[igrid][5] = xgrid[2]; + igrid++; + } +} + +/* ---------------------------------------------------------------------- + memory usage of local data +------------------------------------------------------------------------- */ + +double ComputeGridLocal::memory_usage() +{ + int nbytes = size_local_rows * size_local_cols * sizeof(double); // gridlocal + return nbytes; +} diff --git a/src/ML-SNAP/compute_grid_local.h b/src/ML-SNAP/compute_grid_local.h new file mode 100644 index 0000000000..5bc6f2082a --- /dev/null +++ b/src/ML-SNAP/compute_grid_local.h @@ -0,0 +1,56 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/ Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifndef LMP_COMPUTE_GRID_LOCAL_H +#define LMP_COMPUTE_GRID_LOCAL_H + +#include "compute.h" + +namespace LAMMPS_NS { + +class ComputeGridLocal : public Compute { + public: + ComputeGridLocal(class LAMMPS *, int, char **); + ~ComputeGridLocal() override; + void setup() override; + void compute_local() override = 0; + + double memory_usage() override; + + protected: + int nx, ny, nz; // global grid dimensions + int nxlo, nxhi, nylo, nyhi, nzlo, nzhi; // local grid bounds, inclusive + int nvalues; // number of values per grid point + double **alocal; // pointer to Compute::array_local + int triclinic; // triclinic flag + double *boxlo, *prd; // box info (units real/ortho or reduced/tri) + double *sublo, *subhi; // subdomain info (units real/ortho or reduced/tri) + double delxinv, delyinv, delzinv; // inverse grid spacing + double delx, dely, delz; // grid spacing + int nargbase; // number of base class args + double cutmax; // largest cutoff distance + int size_local_cols_base; // number of columns used for coords, etc. + int gridlocal_allocated; // shows if gridlocal allocated + + void allocate(); // create arrays + void deallocate(); // free arrays + void grid2x(int, int, int, double *); // convert global indices to coordinates + void grid2lamda(int, int, int, double *); // convert global indices to lamda coordinates + void set_grid_global(); // set global grid + void set_grid_local(); // set bounds for local grid + void assign_coords(); // assign coords for grid +}; + +} // namespace LAMMPS_NS + +#endif diff --git a/src/ML-SNAP/compute_sna_atom.cpp b/src/ML-SNAP/compute_sna_atom.cpp index 5f24ebeaa6..1d2190d50d 100644 --- a/src/ML-SNAP/compute_sna_atom.cpp +++ b/src/ML-SNAP/compute_sna_atom.cpp @@ -35,20 +35,21 @@ ComputeSNAAtom::ComputeSNAAtom(LAMMPS *lmp, int narg, char **arg) : radelem(nullptr), wjelem(nullptr), sinnerelem(nullptr), dinnerelem(nullptr) { - double rmin0, rfac0; + // begin code common to all SNAP computes + + double rfac0, rmin0; int twojmax, switchflag, bzeroflag, bnormflag, wselfallflag; int ntypes = atom->ntypes; - int nargmin = 6+2*ntypes; + int nargmin = 6 + 2 * ntypes; - if (narg < nargmin) error->all(FLERR,"Illegal compute sna/atom command"); + if (narg < nargmin) error->all(FLERR, "Illegal compute {} command", style); // default values rmin0 = 0.0; switchflag = 1; bzeroflag = 1; - bnormflag = 0; quadraticflag = 0; chemflag = 0; bnormflag = 0; @@ -56,32 +57,34 @@ ComputeSNAAtom::ComputeSNAAtom(LAMMPS *lmp, int narg, char **arg) : switchinnerflag = 0; nelements = 1; - // offset by 1 to match up with types + // process required arguments - memory->create(radelem,ntypes+1,"sna/atom:radelem"); - memory->create(wjelem,ntypes+1,"sna/atom:wjelem"); + memory->create(radelem, ntypes + 1, "sna/atom:radelem"); // offset by 1 to match up with types + memory->create(wjelem, ntypes + 1, "sna/atom:wjelem"); - rcutfac = atof(arg[3]); - rfac0 = atof(arg[4]); - twojmax = atoi(arg[5]); + rcutfac = utils::numeric(FLERR, arg[3], false, lmp); + rfac0 = utils::numeric(FLERR, arg[4], false, lmp); + twojmax = utils::inumeric(FLERR, arg[5], false, lmp); for (int i = 0; i < ntypes; i++) - radelem[i+1] = atof(arg[6+i]); + radelem[i + 1] = + utils::numeric(FLERR, arg[6 + i], false, lmp); for (int i = 0; i < ntypes; i++) - wjelem[i+1] = atof(arg[6+ntypes+i]); + wjelem[i + 1] = + utils::numeric(FLERR, arg[6 + ntypes + i], false, lmp); // construct cutsq double cut; cutmax = 0.0; - memory->create(cutsq,ntypes+1,ntypes+1,"sna/atom:cutsq"); + memory->create(cutsq, ntypes + 1, ntypes + 1, "sna/atom:cutsq"); for (int i = 1; i <= ntypes; i++) { - cut = 2.0*radelem[i]*rcutfac; + cut = 2.0 * radelem[i] * rcutfac; if (cut > cutmax) cutmax = cut; - cutsq[i][i] = cut*cut; - for (int j = i+1; j <= ntypes; j++) { - cut = (radelem[i]+radelem[j])*rcutfac; - cutsq[i][j] = cutsq[j][i] = cut*cut; + cutsq[i][i] = cut * cut; + for (int j = i + 1; j <= ntypes; j++) { + cut = (radelem[i] + radelem[j]) * rcutfac; + cutsq[i][j] = cutsq[j][i] = cut * cut; } } @@ -95,89 +98,87 @@ ComputeSNAAtom::ComputeSNAAtom(LAMMPS *lmp, int narg, char **arg) : int iarg = nargmin; while (iarg < narg) { - if (strcmp(arg[iarg],"rmin0") == 0) { - if (iarg+2 > narg) - error->all(FLERR,"Illegal compute sna/atom command"); - rmin0 = atof(arg[iarg+1]); + if (strcmp(arg[iarg], "rmin0") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + rmin0 = utils::numeric(FLERR, arg[iarg + 1], false, lmp); iarg += 2; - } else if (strcmp(arg[iarg],"switchflag") == 0) { - if (iarg+2 > narg) - error->all(FLERR,"Illegal compute sna/atom command"); - switchflag = atoi(arg[iarg+1]); + } else if (strcmp(arg[iarg], "switchflag") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + switchflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); iarg += 2; - } else if (strcmp(arg[iarg],"bzeroflag") == 0) { - if (iarg+2 > narg) - error->all(FLERR,"Illegal compute sna/atom command"); - bzeroflag = atoi(arg[iarg+1]); + } else if (strcmp(arg[iarg], "bzeroflag") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + bzeroflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); iarg += 2; - } else if (strcmp(arg[iarg],"quadraticflag") == 0) { - if (iarg+2 > narg) - error->all(FLERR,"Illegal compute sna/atom command"); - quadraticflag = atoi(arg[iarg+1]); + } else if (strcmp(arg[iarg], "quadraticflag") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + quadraticflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); iarg += 2; - } else if (strcmp(arg[iarg],"chem") == 0) { - if (iarg+2+ntypes > narg) - error->all(FLERR,"Illegal compute sna/atom command"); + } else if (strcmp(arg[iarg], "chem") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); chemflag = 1; - memory->create(map,ntypes+1,"compute_sna_atom:map"); - nelements = utils::inumeric(FLERR,arg[iarg+1],false,lmp); + memory->create(map, ntypes + 1, "compute_sna_grid:map"); + nelements = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); for (int i = 0; i < ntypes; i++) { - int jelem = utils::inumeric(FLERR,arg[iarg+2+i],false,lmp); - if (jelem < 0 || jelem >= nelements) - error->all(FLERR,"Illegal compute sna/atom command"); - map[i+1] = jelem; + int jelem = utils::inumeric(FLERR, arg[iarg + 2 + i], false, lmp); + if (jelem < 0 || jelem >= nelements) error->all(FLERR, "Illegal compute {} command", style); + map[i + 1] = jelem; } - iarg += 2+ntypes; - } else if (strcmp(arg[iarg],"bnormflag") == 0) { - if (iarg+2 > narg) - error->all(FLERR,"Illegal compute sna/atom command"); - bnormflag = atoi(arg[iarg+1]); + iarg += 2 + ntypes; + } else if (strcmp(arg[iarg], "bnormflag") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + bnormflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); iarg += 2; - } else if (strcmp(arg[iarg],"wselfallflag") == 0) { - if (iarg+2 > narg) - error->all(FLERR,"Illegal compute sna/atom command"); - wselfallflag = atoi(arg[iarg+1]); + } else if (strcmp(arg[iarg], "wselfallflag") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + wselfallflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); iarg += 2; - } else if (strcmp(arg[iarg],"switchinnerflag") == 0) { - if (iarg+2 > narg) - error->all(FLERR,"Illegal compute sna/atom command"); - switchinnerflag = atoi(arg[iarg+1]); + } else if (strcmp(arg[iarg], "switchinnerflag") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + switchinnerflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); iarg += 2; - } else if (strcmp(arg[iarg],"sinner") == 0) { + } else if (strcmp(arg[iarg], "sinner") == 0) { iarg++; - if (iarg+ntypes > narg) - error->all(FLERR,"Illegal compute sna/atom command"); - memory->create(sinnerelem,ntypes+1,"sna/atom:sinnerelem"); + if (iarg + ntypes > narg) error->all(FLERR, "Illegal compute {} command", style); + memory->create(sinnerelem, ntypes + 1, "snap:sinnerelem"); for (int i = 0; i < ntypes; i++) - sinnerelem[i+1] = utils::numeric(FLERR,arg[iarg+i],false,lmp); + sinnerelem[i + 1] = utils::numeric(FLERR, arg[iarg + i], false, lmp); sinnerflag = 1; iarg += ntypes; - } else if (strcmp(arg[iarg],"dinner") == 0) { + } else if (strcmp(arg[iarg], "dinner") == 0) { iarg++; - if (iarg+ntypes > narg) - error->all(FLERR,"Illegal compute sna/atom command"); - memory->create(dinnerelem,ntypes+1,"sna/atom:dinnerelem"); + if (iarg + ntypes > narg) error->all(FLERR, "Illegal compute {} command", style); + memory->create(dinnerelem, ntypes + 1, "snap:dinnerelem"); for (int i = 0; i < ntypes; i++) - dinnerelem[i+1] = utils::numeric(FLERR,arg[iarg+i],false,lmp); + dinnerelem[i + 1] = utils::numeric(FLERR, arg[iarg + i], false, lmp); dinnerflag = 1; iarg += ntypes; - } else error->all(FLERR,"Illegal compute sna/atom command"); + } else + error->all(FLERR, "Illegal compute {} command", style); } if (switchinnerflag && !(sinnerflag && dinnerflag)) - error->all(FLERR,"Illegal compute sna/atom command: switchinnerflag = 1, missing sinner/dinner keyword"); + error->all( + FLERR, + "Illegal compute {} command: switchinnerflag = 1, missing sinner/dinner keyword", + style); if (!switchinnerflag && (sinnerflag || dinnerflag)) - error->all(FLERR,"Illegal compute sna/atom command: switchinnerflag = 0, unexpected sinner/dinner keyword"); + error->all( + FLERR, + "Illegal compute {} command: switchinnerflag = 0, unexpected sinner/dinner keyword", + style); - snaptr = new SNA(lmp, rfac0, twojmax, - rmin0, switchflag, bzeroflag, - chemflag, bnormflag, wselfallflag, - nelements, switchinnerflag); + snaptr = new SNA(lmp, rfac0, twojmax, rmin0, switchflag, bzeroflag, chemflag, bnormflag, + wselfallflag, nelements, switchinnerflag); ncoeff = snaptr->ncoeff; - size_peratom_cols = ncoeff; - if (quadraticflag) size_peratom_cols += (ncoeff*(ncoeff+1))/2; + nvalues = ncoeff; + if (quadraticflag) nvalues += (ncoeff * (ncoeff + 1)) / 2; + + // end code common to all SNAP computes + + size_peratom_cols = nvalues; peratom_flag = 1; nmax = 0; diff --git a/src/ML-SNAP/compute_sna_atom.h b/src/ML-SNAP/compute_sna_atom.h index 7d1ebfa2f8..059062670f 100644 --- a/src/ML-SNAP/compute_sna_atom.h +++ b/src/ML-SNAP/compute_sna_atom.h @@ -50,6 +50,7 @@ class ComputeSNAAtom : public Compute { class SNA *snaptr; double cutmax; int quadraticflag; + int nvalues; }; } // namespace LAMMPS_NS diff --git a/src/ML-SNAP/compute_sna_grid.cpp b/src/ML-SNAP/compute_sna_grid.cpp new file mode 100644 index 0000000000..497390f99a --- /dev/null +++ b/src/ML-SNAP/compute_sna_grid.cpp @@ -0,0 +1,320 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/ Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#include "compute_sna_grid.h" + +#include "atom.h" +#include "comm.h" +#include "error.h" +#include "force.h" +#include "memory.h" +#include "modify.h" +#include "sna.h" +#include "update.h" + +#include +#include + +using namespace LAMMPS_NS; + +ComputeSNAGrid::ComputeSNAGrid(LAMMPS *lmp, int narg, char **arg) : + ComputeGrid(lmp, narg, arg), cutsq(nullptr), radelem(nullptr), wjelem(nullptr) +{ + // skip over arguments used by base class + // so that argument positions are identical to + // regular per-atom compute + + arg += nargbase; + narg -= nargbase; + + // begin code common to all SNAP computes + + double rfac0, rmin0; + int twojmax, switchflag, bzeroflag, bnormflag, wselfallflag; + + int ntypes = atom->ntypes; + int nargmin = 6 + 2 * ntypes; + + if (narg < nargmin) error->all(FLERR, "Illegal compute {} command", style); + + // default values + + rmin0 = 0.0; + switchflag = 1; + bzeroflag = 1; + quadraticflag = 0; + chemflag = 0; + bnormflag = 0; + wselfallflag = 0; + switchinnerflag = 0; + nelements = 1; + + // process required arguments + + memory->create(radelem, ntypes + 1, "sna/atom:radelem"); // offset by 1 to match up with types + memory->create(wjelem, ntypes + 1, "sna/atom:wjelem"); + + rcutfac = utils::numeric(FLERR, arg[3], false, lmp); + rfac0 = utils::numeric(FLERR, arg[4], false, lmp); + twojmax = utils::inumeric(FLERR, arg[5], false, lmp); + + for (int i = 0; i < ntypes; i++) radelem[i + 1] = utils::numeric(FLERR, arg[6 + i], false, lmp); + for (int i = 0; i < ntypes; i++) + wjelem[i + 1] = utils::numeric(FLERR, arg[6 + ntypes + i], false, lmp); + + // construct cutsq + + double cut; + cutmax = 0.0; + memory->create(cutsq, ntypes + 1, ntypes + 1, "sna/atom:cutsq"); + for (int i = 1; i <= ntypes; i++) { + cut = 2.0 * radelem[i] * rcutfac; + if (cut > cutmax) cutmax = cut; + cutsq[i][i] = cut * cut; + for (int j = i + 1; j <= ntypes; j++) { + cut = (radelem[i] + radelem[j]) * rcutfac; + cutsq[i][j] = cutsq[j][i] = cut * cut; + } + } + + // set local input checks + + int sinnerflag = 0; + int dinnerflag = 0; + + // process optional args + + int iarg = nargmin; + + while (iarg < narg) { + if (strcmp(arg[iarg], "rmin0") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + rmin0 = utils::numeric(FLERR, arg[iarg + 1], false, lmp); + iarg += 2; + } else if (strcmp(arg[iarg], "switchflag") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + switchflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); + iarg += 2; + } else if (strcmp(arg[iarg], "bzeroflag") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + bzeroflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); + iarg += 2; + } else if (strcmp(arg[iarg], "quadraticflag") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + quadraticflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); + iarg += 2; + } else if (strcmp(arg[iarg], "chem") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + chemflag = 1; + memory->create(map, ntypes + 1, "compute_sna_grid:map"); + nelements = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); + for (int i = 0; i < ntypes; i++) { + int jelem = utils::inumeric(FLERR, arg[iarg + 2 + i], false, lmp); + if (jelem < 0 || jelem >= nelements) error->all(FLERR, "Illegal compute {} command", style); + map[i + 1] = jelem; + } + iarg += 2 + ntypes; + } else if (strcmp(arg[iarg], "bnormflag") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + bnormflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); + iarg += 2; + } else if (strcmp(arg[iarg], "wselfallflag") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + wselfallflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); + iarg += 2; + } else if (strcmp(arg[iarg], "switchinnerflag") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + switchinnerflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); + iarg += 2; + } else if (strcmp(arg[iarg], "sinner") == 0) { + iarg++; + if (iarg + ntypes > narg) error->all(FLERR, "Illegal compute {} command", style); + memory->create(sinnerelem, ntypes + 1, "snap:sinnerelem"); + for (int i = 0; i < ntypes; i++) + sinnerelem[i + 1] = utils::numeric(FLERR, arg[iarg + i], false, lmp); + sinnerflag = 1; + iarg += ntypes; + } else if (strcmp(arg[iarg], "dinner") == 0) { + iarg++; + if (iarg + ntypes > narg) error->all(FLERR, "Illegal compute {} command", style); + memory->create(dinnerelem, ntypes + 1, "snap:dinnerelem"); + for (int i = 0; i < ntypes; i++) + dinnerelem[i + 1] = utils::numeric(FLERR, arg[iarg + i], false, lmp); + dinnerflag = 1; + iarg += ntypes; + } else + error->all(FLERR, "Illegal compute {} command", style); + } + + if (switchinnerflag && !(sinnerflag && dinnerflag)) + error->all(FLERR, + "Illegal compute {} command: switchinnerflag = 1, missing sinner/dinner keyword", + style); + + if (!switchinnerflag && (sinnerflag || dinnerflag)) + error->all(FLERR, + "Illegal compute {} command: switchinnerflag = 0, unexpected sinner/dinner keyword", + style); + + snaptr = new SNA(lmp, rfac0, twojmax, rmin0, switchflag, bzeroflag, chemflag, bnormflag, + wselfallflag, nelements, switchinnerflag); + + ncoeff = snaptr->ncoeff; + nvalues = ncoeff; + if (quadraticflag) nvalues += (ncoeff * (ncoeff + 1)) / 2; + + // end code common to all SNAP computes + + size_array_cols = size_array_cols_base + nvalues; + array_flag = 1; +} + +/* ---------------------------------------------------------------------- */ + +ComputeSNAGrid::~ComputeSNAGrid() +{ + memory->destroy(radelem); + memory->destroy(wjelem); + memory->destroy(cutsq); + delete snaptr; + + if (chemflag) memory->destroy(map); +} + +/* ---------------------------------------------------------------------- */ + +void ComputeSNAGrid::init() +{ + if ((modify->get_compute_by_style("^sna/grid$").size() > 1) && (comm->me == 0)) + error->warning(FLERR, "More than one instance of compute sna/grid"); + snaptr->init(); +} + +/* ---------------------------------------------------------------------- */ + +void ComputeSNAGrid::compute_array() +{ + invoked_array = update->ntimestep; + + // compute sna for each gridpoint + + double **const x = atom->x; + const int *const mask = atom->mask; + int *const type = atom->type; + const int ntotal = atom->nlocal + atom->nghost; + + // insure rij, inside, and typej are of size jnum + + snaptr->grow_rij(ntotal); + + for (int iz = nzlo; iz <= nzhi; iz++) + for (int iy = nylo; iy <= nyhi; iy++) + for (int ix = nxlo; ix <= nxhi; ix++) { + double xgrid[3]; + const int igrid = iz * (nx * ny) + iy * nx + ix; + grid2x(igrid, xgrid); + const double xtmp = xgrid[0]; + const double ytmp = xgrid[1]; + const double ztmp = xgrid[2]; + + // currently, all grid points are type 1 + // not clear what a better choice would be + + const int itype = 1; + int ielem = 0; + if (chemflag) ielem = map[itype]; + + // rij[][3] = displacements between atom I and those neighbors + // inside = indices of neighbors of I within cutoff + // typej = types of neighbors of I within cutoff + + int ninside = 0; + for (int j = 0; j < ntotal; j++) { + + // check that j is in compute group + + if (!(mask[j] & groupbit)) continue; + + const double delx = xtmp - x[j][0]; + const double dely = ytmp - x[j][1]; + const double delz = ztmp - x[j][2]; + const double rsq = delx * delx + dely * dely + delz * delz; + int jtype = type[j]; + int jelem = 0; + if (chemflag) jelem = map[jtype]; + + if (rsq < cutsq[jtype][jtype] && rsq > 1e-20) { + snaptr->rij[ninside][0] = delx; + snaptr->rij[ninside][1] = dely; + snaptr->rij[ninside][2] = delz; + snaptr->inside[ninside] = j; + snaptr->wj[ninside] = wjelem[jtype]; + snaptr->rcutij[ninside] = 2.0 * radelem[jtype] * rcutfac; + if (switchinnerflag) { + snaptr->sinnerij[ninside] = sinnerelem[jelem]; + snaptr->dinnerij[ninside] = dinnerelem[jelem]; + } + if (chemflag) snaptr->element[ninside] = jelem; + ninside++; + } + } + + snaptr->compute_ui(ninside, ielem); + snaptr->compute_zi(); + snaptr->compute_bi(ielem); + + // linear contributions + + for (int icoeff = 0; icoeff < ncoeff; icoeff++) + gridlocal[size_array_cols_base + icoeff][iz][iy][ix] = snaptr->blist[icoeff]; + + // quadratic contributions + + if (quadraticflag) { + int ncount = ncoeff; + for (int icoeff = 0; icoeff < ncoeff; icoeff++) { + double bveci = snaptr->blist[icoeff]; + gridlocal[size_array_cols_base + ncount++][iz][iy][ix] = 0.5 * bveci * bveci; + for (int jcoeff = icoeff + 1; jcoeff < ncoeff; jcoeff++) + gridlocal[size_array_cols_base + ncount++][iz][iy][ix] = + bveci * snaptr->blist[jcoeff]; + } + } + } + + memset(&grid[0][0], 0, size_array_rows * size_array_cols * sizeof(double)); + + for (int iz = nzlo; iz <= nzhi; iz++) + for (int iy = nylo; iy <= nyhi; iy++) + for (int ix = nxlo; ix <= nxhi; ix++) { + const int igrid = iz * (nx * ny) + iy * nx + ix; + for (int j = 0; j < nvalues; j++) + grid[igrid][size_array_cols_base + j] = gridlocal[size_array_cols_base + j][iz][iy][ix]; + } + MPI_Allreduce(&grid[0][0], &gridall[0][0], size_array_rows * size_array_cols, MPI_DOUBLE, MPI_SUM, + world); + assign_coords_all(); +} + +/* ---------------------------------------------------------------------- + memory usage +------------------------------------------------------------------------- */ + +double ComputeSNAGrid::memory_usage() +{ + double nbytes = snaptr->memory_usage(); // SNA object + int n = atom->ntypes + 1; + nbytes += (double) n * sizeof(int); // map + + return nbytes; +} diff --git a/src/ML-SNAP/compute_sna_grid.h b/src/ML-SNAP/compute_sna_grid.h new file mode 100644 index 0000000000..839003dc2e --- /dev/null +++ b/src/ML-SNAP/compute_sna_grid.h @@ -0,0 +1,54 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/ Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef COMPUTE_CLASS +// clang-format off +ComputeStyle(sna/grid,ComputeSNAGrid); +// clang-format on +#else + +#ifndef LMP_COMPUTE_SNA_GRID_H +#define LMP_COMPUTE_SNA_GRID_H + +#include "compute_grid.h" + +namespace LAMMPS_NS { + +class ComputeSNAGrid : public ComputeGrid { + public: + ComputeSNAGrid(class LAMMPS *, int, char **); + ~ComputeSNAGrid() override; + void init() override; + void compute_array() override; + double memory_usage() override; + + private: + int ncoeff; + double **cutsq; + double rcutfac; + double *radelem; + double *wjelem; + int *map; // map types to [0,nelements) + int nelements, chemflag; + int switchinnerflag; + double *sinnerelem; + double *dinnerelem; + class SNA *snaptr; + double cutmax; + int quadraticflag; +}; + +} // namespace LAMMPS_NS + +#endif +#endif diff --git a/src/ML-SNAP/compute_sna_grid_local.cpp b/src/ML-SNAP/compute_sna_grid_local.cpp new file mode 100644 index 0000000000..80a1baddab --- /dev/null +++ b/src/ML-SNAP/compute_sna_grid_local.cpp @@ -0,0 +1,306 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/ Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#include "compute_sna_grid_local.h" + +#include "atom.h" +#include "comm.h" +#include "error.h" +#include "force.h" +#include "memory.h" +#include "modify.h" +#include "sna.h" +#include "update.h" + +#include +#include + +using namespace LAMMPS_NS; + +ComputeSNAGridLocal::ComputeSNAGridLocal(LAMMPS *lmp, int narg, char **arg) : + ComputeGridLocal(lmp, narg, arg), cutsq(nullptr), radelem(nullptr), wjelem(nullptr) +{ + // skip over arguments used by base class + // so that argument positions are identical to + // regular per-atom compute + + arg += nargbase; + narg -= nargbase; + + // begin code common to all SNAP computes + + double rfac0, rmin0; + int twojmax, switchflag, bzeroflag, bnormflag, wselfallflag; + + int ntypes = atom->ntypes; + int nargmin = 6 + 2 * ntypes; + + if (narg < nargmin) error->all(FLERR, "Illegal compute {} command", style); + + // default values + + rmin0 = 0.0; + switchflag = 1; + bzeroflag = 1; + quadraticflag = 0; + chemflag = 0; + bnormflag = 0; + wselfallflag = 0; + switchinnerflag = 0; + nelements = 1; + + // process required arguments + + memory->create(radelem, ntypes + 1, "sna/atom:radelem"); // offset by 1 to match up with types + memory->create(wjelem, ntypes + 1, "sna/atom:wjelem"); + + rcutfac = utils::numeric(FLERR, arg[3], false, lmp); + rfac0 = utils::numeric(FLERR, arg[4], false, lmp); + twojmax = utils::inumeric(FLERR, arg[5], false, lmp); + + for (int i = 0; i < ntypes; i++) radelem[i + 1] = utils::numeric(FLERR, arg[6 + i], false, lmp); + for (int i = 0; i < ntypes; i++) + wjelem[i + 1] = utils::numeric(FLERR, arg[6 + ntypes + i], false, lmp); + + // construct cutsq + + double cut; + cutmax = 0.0; + memory->create(cutsq, ntypes + 1, ntypes + 1, "sna/atom:cutsq"); + for (int i = 1; i <= ntypes; i++) { + cut = 2.0 * radelem[i] * rcutfac; + if (cut > cutmax) cutmax = cut; + cutsq[i][i] = cut * cut; + for (int j = i + 1; j <= ntypes; j++) { + cut = (radelem[i] + radelem[j]) * rcutfac; + cutsq[i][j] = cutsq[j][i] = cut * cut; + } + } + + // set local input checks + + int sinnerflag = 0; + int dinnerflag = 0; + + // process optional args + + int iarg = nargmin; + + while (iarg < narg) { + if (strcmp(arg[iarg], "rmin0") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + rmin0 = utils::numeric(FLERR, arg[iarg + 1], false, lmp); + iarg += 2; + } else if (strcmp(arg[iarg], "switchflag") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + switchflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); + iarg += 2; + } else if (strcmp(arg[iarg], "bzeroflag") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + bzeroflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); + iarg += 2; + } else if (strcmp(arg[iarg], "quadraticflag") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + quadraticflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); + iarg += 2; + } else if (strcmp(arg[iarg], "chem") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + chemflag = 1; + memory->create(map, ntypes + 1, "compute_sna_grid:map"); + nelements = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); + for (int i = 0; i < ntypes; i++) { + int jelem = utils::inumeric(FLERR, arg[iarg + 2 + i], false, lmp); + if (jelem < 0 || jelem >= nelements) error->all(FLERR, "Illegal compute {} command", style); + map[i + 1] = jelem; + } + iarg += 2 + ntypes; + } else if (strcmp(arg[iarg], "bnormflag") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + bnormflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); + iarg += 2; + } else if (strcmp(arg[iarg], "wselfallflag") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + wselfallflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); + iarg += 2; + } else if (strcmp(arg[iarg], "switchinnerflag") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + switchinnerflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); + iarg += 2; + } else if (strcmp(arg[iarg], "sinner") == 0) { + iarg++; + if (iarg + ntypes > narg) error->all(FLERR, "Illegal compute {} command", style); + memory->create(sinnerelem, ntypes + 1, "snap:sinnerelem"); + for (int i = 0; i < ntypes; i++) + sinnerelem[i + 1] = utils::numeric(FLERR, arg[iarg + i], false, lmp); + sinnerflag = 1; + iarg += ntypes; + } else if (strcmp(arg[iarg], "dinner") == 0) { + iarg++; + if (iarg + ntypes > narg) error->all(FLERR, "Illegal compute {} command", style); + memory->create(dinnerelem, ntypes + 1, "snap:dinnerelem"); + for (int i = 0; i < ntypes; i++) + dinnerelem[i + 1] = utils::numeric(FLERR, arg[iarg + i], false, lmp); + dinnerflag = 1; + iarg += ntypes; + } else + error->all(FLERR, "Illegal compute {} command", style); + } + + if (switchinnerflag && !(sinnerflag && dinnerflag)) + error->all(FLERR, + "Illegal compute {} command: switchinnerflag = 1, missing sinner/dinner keyword", + style); + + if (!switchinnerflag && (sinnerflag || dinnerflag)) + error->all(FLERR, + "Illegal compute {} command: switchinnerflag = 0, unexpected sinner/dinner keyword", + style); + + snaptr = new SNA(lmp, rfac0, twojmax, rmin0, switchflag, bzeroflag, chemflag, bnormflag, + wselfallflag, nelements, switchinnerflag); + + ncoeff = snaptr->ncoeff; + nvalues = ncoeff; + if (quadraticflag) nvalues += (ncoeff * (ncoeff + 1)) / 2; + + // end code common to all SNAP computes + + size_local_cols = size_local_cols_base + nvalues; +} + +/* ---------------------------------------------------------------------- */ + +ComputeSNAGridLocal::~ComputeSNAGridLocal() +{ + memory->destroy(radelem); + memory->destroy(wjelem); + memory->destroy(cutsq); + delete snaptr; + + if (chemflag) memory->destroy(map); +} + +/* ---------------------------------------------------------------------- */ + +void ComputeSNAGridLocal::init() +{ + if ((modify->get_compute_by_style("^sna/grid/local$").size() > 1) && (comm->me == 0)) + error->warning(FLERR, "More than one instance of compute sna/grid/local"); + snaptr->init(); +} + +/* ---------------------------------------------------------------------- */ + +void ComputeSNAGridLocal::compute_local() +{ + invoked_array = update->ntimestep; + + // compute sna for each gridpoint + + double **const x = atom->x; + const int *const mask = atom->mask; + int *const type = atom->type; + const int ntotal = atom->nlocal + atom->nghost; + + // insure rij, inside, and typej are of size jnum + + snaptr->grow_rij(ntotal); + + int igrid = 0; + for (int iz = nzlo; iz <= nzhi; iz++) + for (int iy = nylo; iy <= nyhi; iy++) + for (int ix = nxlo; ix <= nxhi; ix++) { + double xgrid[3]; + grid2x(ix, iy, iz, xgrid); + const double xtmp = xgrid[0]; + const double ytmp = xgrid[1]; + const double ztmp = xgrid[2]; + + // currently, all grid points are type 1 + // not clear what a better choice would be + + const int itype = 1; + int ielem = 0; + if (chemflag) ielem = map[itype]; + + // rij[][3] = displacements between atom I and those neighbors + // inside = indices of neighbors of I within cutoff + // typej = types of neighbors of I within cutoff + + int ninside = 0; + for (int j = 0; j < ntotal; j++) { + + // check that j is in compute group + + if (!(mask[j] & groupbit)) continue; + + const double delx = xtmp - x[j][0]; + const double dely = ytmp - x[j][1]; + const double delz = ztmp - x[j][2]; + const double rsq = delx * delx + dely * dely + delz * delz; + int jtype = type[j]; + int jelem = 0; + if (chemflag) jelem = map[jtype]; + if (rsq < cutsq[jtype][jtype] && rsq > 1e-20) { + snaptr->rij[ninside][0] = delx; + snaptr->rij[ninside][1] = dely; + snaptr->rij[ninside][2] = delz; + snaptr->inside[ninside] = j; + snaptr->wj[ninside] = wjelem[jtype]; + snaptr->rcutij[ninside] = 2.0 * radelem[jtype] * rcutfac; + if (switchinnerflag) { + snaptr->sinnerij[ninside] = sinnerelem[jelem]; + snaptr->dinnerij[ninside] = dinnerelem[jelem]; + } + if (chemflag) + snaptr->element[ninside] = jelem; // element index for multi-element snap + ninside++; + } + } + + snaptr->compute_ui(ninside, ielem); + snaptr->compute_zi(); + snaptr->compute_bi(ielem); + + // linear contributions + + for (int icoeff = 0; icoeff < ncoeff; icoeff++) + alocal[igrid][size_local_cols_base + icoeff] = snaptr->blist[icoeff]; + + // quadratic contributions + + if (quadraticflag) { + int ncount = ncoeff; + for (int icoeff = 0; icoeff < ncoeff; icoeff++) { + double bveci = snaptr->blist[icoeff]; + alocal[igrid][size_local_cols_base + ncount++] = 0.5 * bveci * bveci; + for (int jcoeff = icoeff + 1; jcoeff < ncoeff; jcoeff++) + alocal[igrid][size_local_cols_base + ncount++] = bveci * snaptr->blist[jcoeff]; + } + } + igrid++; + } +} + +/* ---------------------------------------------------------------------- + memory usage +------------------------------------------------------------------------- */ + +double ComputeSNAGridLocal::memory_usage() +{ + double nbytes = snaptr->memory_usage(); // SNA object + int n = atom->ntypes + 1; + nbytes += (double) n * sizeof(int); // map + + return nbytes; +} diff --git a/src/ML-SNAP/compute_sna_grid_local.h b/src/ML-SNAP/compute_sna_grid_local.h new file mode 100644 index 0000000000..60f93b92b0 --- /dev/null +++ b/src/ML-SNAP/compute_sna_grid_local.h @@ -0,0 +1,54 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/ Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef COMPUTE_CLASS +// clang-format off +ComputeStyle(sna/grid/local,ComputeSNAGridLocal); +// clang-format on +#else + +#ifndef LMP_COMPUTE_SNA_GRID_LOCAL_H +#define LMP_COMPUTE_SNA_GRID_LOCAL_H + +#include "compute_grid_local.h" + +namespace LAMMPS_NS { + +class ComputeSNAGridLocal : public ComputeGridLocal { + public: + ComputeSNAGridLocal(class LAMMPS *, int, char **); + ~ComputeSNAGridLocal() override; + void init() override; + void compute_local() override; + double memory_usage() override; + + private: + int ncoeff; + double **cutsq; + double rcutfac; + double *radelem; + double *wjelem; + int *map; // map types to [0,nelements) + int nelements, chemflag; + int switchinnerflag; + double *sinnerelem; + double *dinnerelem; + class SNA *snaptr; + double cutmax; + int quadraticflag; +}; + +} // namespace LAMMPS_NS + +#endif +#endif diff --git a/src/ML-SNAP/compute_snad_atom.cpp b/src/ML-SNAP/compute_snad_atom.cpp index 838d8a85c9..12839629a9 100644 --- a/src/ML-SNAP/compute_snad_atom.cpp +++ b/src/ML-SNAP/compute_snad_atom.cpp @@ -34,20 +34,22 @@ ComputeSNADAtom::ComputeSNADAtom(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg), cutsq(nullptr), list(nullptr), snad(nullptr), radelem(nullptr), wjelem(nullptr), sinnerelem(nullptr), dinnerelem(nullptr) { + + // begin code common to all SNAP computes + double rfac0, rmin0; int twojmax, switchflag, bzeroflag, bnormflag, wselfallflag; int ntypes = atom->ntypes; - int nargmin = 6+2*ntypes; + int nargmin = 6 + 2 * ntypes; - if (narg < nargmin) error->all(FLERR,"Illegal compute snad/atom command"); + if (narg < nargmin) error->all(FLERR, "Illegal compute {} command", style); // default values rmin0 = 0.0; switchflag = 1; bzeroflag = 1; - bnormflag = 0; quadraticflag = 0; chemflag = 0; bnormflag = 0; @@ -57,28 +59,32 @@ ComputeSNADAtom::ComputeSNADAtom(LAMMPS *lmp, int narg, char **arg) : // process required arguments - memory->create(radelem,ntypes+1,"snad/atom:radelem"); // offset by 1 to match up with types - memory->create(wjelem,ntypes+1,"snad/atom:wjelem"); - rcutfac = atof(arg[3]); - rfac0 = atof(arg[4]); - twojmax = atoi(arg[5]); + memory->create(radelem, ntypes + 1, "sna/atom:radelem"); // offset by 1 to match up with types + memory->create(wjelem, ntypes + 1, "sna/atom:wjelem"); + + rcutfac = utils::numeric(FLERR, arg[3], false, lmp); + rfac0 = utils::numeric(FLERR, arg[4], false, lmp); + twojmax = utils::inumeric(FLERR, arg[5], false, lmp); + for (int i = 0; i < ntypes; i++) - radelem[i+1] = atof(arg[6+i]); + radelem[i + 1] = + utils::numeric(FLERR, arg[6 + i], false, lmp); for (int i = 0; i < ntypes; i++) - wjelem[i+1] = atof(arg[6+ntypes+i]); + wjelem[i + 1] = + utils::numeric(FLERR, arg[6 + ntypes + i], false, lmp); // construct cutsq double cut; cutmax = 0.0; - memory->create(cutsq,ntypes+1,ntypes+1,"snad/atom:cutsq"); + memory->create(cutsq, ntypes + 1, ntypes + 1, "sna/atom:cutsq"); for (int i = 1; i <= ntypes; i++) { - cut = 2.0*radelem[i]*rcutfac; + cut = 2.0 * radelem[i] * rcutfac; if (cut > cutmax) cutmax = cut; - cutsq[i][i] = cut*cut; - for (int j = i+1; j <= ntypes; j++) { - cut = (radelem[i]+radelem[j])*rcutfac; - cutsq[i][j] = cutsq[j][i] = cut*cut; + cutsq[i][i] = cut * cut; + for (int j = i + 1; j <= ntypes; j++) { + cut = (radelem[i] + radelem[j]) * rcutfac; + cutsq[i][j] = cutsq[j][i] = cut * cut; } } @@ -92,93 +98,89 @@ ComputeSNADAtom::ComputeSNADAtom(LAMMPS *lmp, int narg, char **arg) : int iarg = nargmin; while (iarg < narg) { - if (strcmp(arg[iarg],"rmin0") == 0) { - if (iarg+2 > narg) - error->all(FLERR,"Illegal compute snad/atom command"); - rmin0 = atof(arg[iarg+1]); + if (strcmp(arg[iarg], "rmin0") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + rmin0 = utils::numeric(FLERR, arg[iarg + 1], false, lmp); iarg += 2; - } else if (strcmp(arg[iarg],"bzeroflag") == 0) { - if (iarg+2 > narg) - error->all(FLERR,"Illegal compute snad/atom command"); - bzeroflag = atoi(arg[iarg+1]); + } else if (strcmp(arg[iarg], "switchflag") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + switchflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); iarg += 2; - } else if (strcmp(arg[iarg],"switchflag") == 0) { - if (iarg+2 > narg) - error->all(FLERR,"Illegal compute snad/atom command"); - switchflag = atoi(arg[iarg+1]); + } else if (strcmp(arg[iarg], "bzeroflag") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + bzeroflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); iarg += 2; - } else if (strcmp(arg[iarg],"quadraticflag") == 0) { - if (iarg+2 > narg) - error->all(FLERR,"Illegal compute snad/atom command"); - quadraticflag = atoi(arg[iarg+1]); + } else if (strcmp(arg[iarg], "quadraticflag") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + quadraticflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); iarg += 2; - } else if (strcmp(arg[iarg],"chem") == 0) { - if (iarg+2+ntypes > narg) - error->all(FLERR,"Illegal compute snad/atom command"); + } else if (strcmp(arg[iarg], "chem") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); chemflag = 1; - memory->create(map,ntypes+1,"compute_snad_atom:map"); - nelements = utils::inumeric(FLERR,arg[iarg+1],false,lmp); + memory->create(map, ntypes + 1, "compute_sna_grid:map"); + nelements = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); for (int i = 0; i < ntypes; i++) { - int jelem = utils::inumeric(FLERR,arg[iarg+2+i],false,lmp); - if (jelem < 0 || jelem >= nelements) - error->all(FLERR,"Illegal compute snad/atom command"); - map[i+1] = jelem; + int jelem = utils::inumeric(FLERR, arg[iarg + 2 + i], false, lmp); + if (jelem < 0 || jelem >= nelements) error->all(FLERR, "Illegal compute {} command", style); + map[i + 1] = jelem; } - iarg += 2+ntypes; - } else if (strcmp(arg[iarg],"bnormflag") == 0) { - if (iarg+2 > narg) - error->all(FLERR,"Illegal compute snad/atom command"); - bnormflag = atoi(arg[iarg+1]); + iarg += 2 + ntypes; + } else if (strcmp(arg[iarg], "bnormflag") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + bnormflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); iarg += 2; - } else if (strcmp(arg[iarg],"wselfallflag") == 0) { - if (iarg+2 > narg) - error->all(FLERR,"Illegal compute snad/atom command"); - wselfallflag = atoi(arg[iarg+1]); + } else if (strcmp(arg[iarg], "wselfallflag") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + wselfallflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); iarg += 2; - } else if (strcmp(arg[iarg],"switchinnerflag") == 0) { - if (iarg+2 > narg) - error->all(FLERR,"Illegal compute snad/atom command"); - switchinnerflag = atoi(arg[iarg+1]); + } else if (strcmp(arg[iarg], "switchinnerflag") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + switchinnerflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); iarg += 2; - } else if (strcmp(arg[iarg],"sinner") == 0) { + } else if (strcmp(arg[iarg], "sinner") == 0) { iarg++; - if (iarg+ntypes > narg) - error->all(FLERR,"Illegal compute snad/atom command"); - memory->create(sinnerelem,ntypes+1,"snad/atom:sinnerelem"); + if (iarg + ntypes > narg) error->all(FLERR, "Illegal compute {} command", style); + memory->create(sinnerelem, ntypes + 1, "snap:sinnerelem"); for (int i = 0; i < ntypes; i++) - sinnerelem[i+1] = utils::numeric(FLERR,arg[iarg+i],false,lmp); + sinnerelem[i + 1] = utils::numeric(FLERR, arg[iarg + i], false, lmp); sinnerflag = 1; iarg += ntypes; - } else if (strcmp(arg[iarg],"dinner") == 0) { + } else if (strcmp(arg[iarg], "dinner") == 0) { iarg++; - if (iarg+ntypes > narg) - error->all(FLERR,"Illegal compute snad/atom command"); - memory->create(dinnerelem,ntypes+1,"snad/atom:dinnerelem"); + if (iarg + ntypes > narg) error->all(FLERR, "Illegal compute {} command", style); + memory->create(dinnerelem, ntypes + 1, "snap:dinnerelem"); for (int i = 0; i < ntypes; i++) - dinnerelem[i+1] = utils::numeric(FLERR,arg[iarg+i],false,lmp); + dinnerelem[i + 1] = utils::numeric(FLERR, arg[iarg + i], false, lmp); dinnerflag = 1; iarg += ntypes; - } else error->all(FLERR,"Illegal compute snad/atom command"); + } else + error->all(FLERR, "Illegal compute {} command", style); } if (switchinnerflag && !(sinnerflag && dinnerflag)) - error->all(FLERR,"Illegal compute snad/atom command: switchinnerflag = 1, missing sinner/dinner keyword"); + error->all( + FLERR, + "Illegal compute {} command: switchinnerflag = 1, missing sinner/dinner keyword", + style); if (!switchinnerflag && (sinnerflag || dinnerflag)) - error->all(FLERR,"Illegal compute snad/atom command: switchinnerflag = 0, unexpected sinner/dinner keyword"); + error->all( + FLERR, + "Illegal compute {} command: switchinnerflag = 0, unexpected sinner/dinner keyword", + style); - - snaptr = new SNA(lmp, rfac0, twojmax, - rmin0, switchflag, bzeroflag, - chemflag, bnormflag, wselfallflag, - nelements, switchinnerflag); + snaptr = new SNA(lmp, rfac0, twojmax, rmin0, switchflag, bzeroflag, chemflag, bnormflag, + wselfallflag, nelements, switchinnerflag); ncoeff = snaptr->ncoeff; - nperdim = ncoeff; - if (quadraticflag) nperdim += (ncoeff*(ncoeff+1))/2; - yoffset = nperdim; - zoffset = 2*nperdim; - size_peratom_cols = 3*nperdim*atom->ntypes; + nvalues = ncoeff; + if (quadraticflag) nvalues += (ncoeff * (ncoeff + 1)) / 2; + + // end code common to all SNAP computes + + yoffset = nvalues; + zoffset = 2*nvalues; + size_peratom_cols = 3*nvalues*atom->ntypes; comm_reverse = size_peratom_cols; peratom_flag = 1; @@ -289,7 +291,7 @@ void ComputeSNADAtom::compute_peratom() // const int typeoffset = threencoeff*(atom->type[i]-1); // const int quadraticoffset = threencoeff*atom->ntypes + // threencoeffq*(atom->type[i]-1); - const int typeoffset = 3*nperdim*(atom->type[i]-1); + const int typeoffset = 3*nvalues*(atom->type[i]-1); // insure rij, inside, and typej are of size jnum diff --git a/src/ML-SNAP/compute_snad_atom.h b/src/ML-SNAP/compute_snad_atom.h index 0cdb075daf..ac27e10769 100644 --- a/src/ML-SNAP/compute_snad_atom.h +++ b/src/ML-SNAP/compute_snad_atom.h @@ -37,7 +37,7 @@ class ComputeSNADAtom : public Compute { private: int nmax; - int ncoeff, nperdim, yoffset, zoffset; + int ncoeff, nvalues, yoffset, zoffset; double **cutsq; class NeighList *list; double **snad; diff --git a/src/ML-SNAP/compute_snap.cpp b/src/ML-SNAP/compute_snap.cpp index 83f56e338c..d1f75f8382 100644 --- a/src/ML-SNAP/compute_snap.cpp +++ b/src/ML-SNAP/compute_snap.cpp @@ -41,13 +41,15 @@ ComputeSnap::ComputeSnap(LAMMPS *lmp, int narg, char **arg) : array_flag = 1; extarray = 0; + // begin code common to all SNAP computes + double rfac0, rmin0; int twojmax, switchflag, bzeroflag, bnormflag, wselfallflag; int ntypes = atom->ntypes; - int nargmin = 6+2*ntypes; + int nargmin = 6 + 2 * ntypes; - if (narg < nargmin) error->all(FLERR,"Illegal compute snap command"); + if (narg < nargmin) error->all(FLERR, "Illegal compute {} command", style); // default values @@ -55,7 +57,6 @@ ComputeSnap::ComputeSnap(LAMMPS *lmp, int narg, char **arg) : switchflag = 1; bzeroflag = 1; quadraticflag = 0; - bikflag = 0; chemflag = 0; bnormflag = 0; wselfallflag = 0; @@ -64,28 +65,32 @@ ComputeSnap::ComputeSnap(LAMMPS *lmp, int narg, char **arg) : // process required arguments - memory->create(radelem,ntypes+1,"snap:radelem"); // offset by 1 to match up with types - memory->create(wjelem,ntypes+1,"snap:wjelem"); - rcutfac = atof(arg[3]); - rfac0 = atof(arg[4]); - twojmax = atoi(arg[5]); + memory->create(radelem, ntypes + 1, "sna/atom:radelem"); // offset by 1 to match up with types + memory->create(wjelem, ntypes + 1, "sna/atom:wjelem"); + + rcutfac = utils::numeric(FLERR, arg[3], false, lmp); + rfac0 = utils::numeric(FLERR, arg[4], false, lmp); + twojmax = utils::inumeric(FLERR, arg[5], false, lmp); + for (int i = 0; i < ntypes; i++) - radelem[i+1] = atof(arg[6+i]); + radelem[i + 1] = + utils::numeric(FLERR, arg[6 + i], false, lmp); for (int i = 0; i < ntypes; i++) - wjelem[i+1] = atof(arg[6+ntypes+i]); + wjelem[i + 1] = + utils::numeric(FLERR, arg[6 + ntypes + i], false, lmp); // construct cutsq double cut; cutmax = 0.0; - memory->create(cutsq,ntypes+1,ntypes+1,"snap:cutsq"); + memory->create(cutsq, ntypes + 1, ntypes + 1, "sna/atom:cutsq"); for (int i = 1; i <= ntypes; i++) { - cut = 2.0*radelem[i]*rcutfac; + cut = 2.0 * radelem[i] * rcutfac; if (cut > cutmax) cutmax = cut; - cutsq[i][i] = cut*cut; - for (int j = i+1; j <= ntypes; j++) { - cut = (radelem[i]+radelem[j])*rcutfac; - cutsq[i][j] = cutsq[j][i] = cut*cut; + cutsq[i][i] = cut * cut; + for (int j = i + 1; j <= ntypes; j++) { + cut = (radelem[i] + radelem[j]) * rcutfac; + cutsq[i][j] = cutsq[j][i] = cut * cut; } } @@ -99,107 +104,99 @@ ComputeSnap::ComputeSnap(LAMMPS *lmp, int narg, char **arg) : int iarg = nargmin; while (iarg < narg) { - if (strcmp(arg[iarg],"rmin0") == 0) { - if (iarg+2 > narg) - error->all(FLERR,"Illegal compute snap command"); - rmin0 = atof(arg[iarg+1]); + if (strcmp(arg[iarg], "rmin0") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + rmin0 = utils::numeric(FLERR, arg[iarg + 1], false, lmp); iarg += 2; - } else if (strcmp(arg[iarg],"bzeroflag") == 0) { - if (iarg+2 > narg) - error->all(FLERR,"Illegal compute snap command"); - bzeroflag = atoi(arg[iarg+1]); + } else if (strcmp(arg[iarg], "switchflag") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + switchflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); iarg += 2; - } else if (strcmp(arg[iarg],"switchflag") == 0) { - if (iarg+2 > narg) - error->all(FLERR,"Illegal compute snap command"); - switchflag = atoi(arg[iarg+1]); + } else if (strcmp(arg[iarg], "bzeroflag") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + bzeroflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); iarg += 2; - } else if (strcmp(arg[iarg],"quadraticflag") == 0) { - if (iarg+2 > narg) - error->all(FLERR,"Illegal compute snap command"); - quadraticflag = atoi(arg[iarg+1]); + } else if (strcmp(arg[iarg], "quadraticflag") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + quadraticflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); iarg += 2; - } else if (strcmp(arg[iarg],"chem") == 0) { - if (iarg+2+ntypes > narg) - error->all(FLERR,"Illegal compute snap command"); + } else if (strcmp(arg[iarg], "chem") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); chemflag = 1; - memory->create(map,ntypes+1,"compute_snap:map"); - nelements = utils::inumeric(FLERR,arg[iarg+1],false,lmp); + memory->create(map, ntypes + 1, "compute_sna_grid:map"); + nelements = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); for (int i = 0; i < ntypes; i++) { - int jelem = utils::inumeric(FLERR,arg[iarg+2+i],false,lmp); - if (jelem < 0 || jelem >= nelements) - error->all(FLERR,"Illegal compute snap command"); - map[i+1] = jelem; + int jelem = utils::inumeric(FLERR, arg[iarg + 2 + i], false, lmp); + if (jelem < 0 || jelem >= nelements) error->all(FLERR, "Illegal compute {} command", style); + map[i + 1] = jelem; } - iarg += 2+ntypes; - } else if (strcmp(arg[iarg],"bnormflag") == 0) { - if (iarg+2 > narg) - error->all(FLERR,"Illegal compute snap command"); - bnormflag = atoi(arg[iarg+1]); + iarg += 2 + ntypes; + } else if (strcmp(arg[iarg], "bnormflag") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + bnormflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); iarg += 2; - } else if (strcmp(arg[iarg],"wselfallflag") == 0) { - if (iarg+2 > narg) - error->all(FLERR,"Illegal compute snap command"); - wselfallflag = atoi(arg[iarg+1]); + } else if (strcmp(arg[iarg], "wselfallflag") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + wselfallflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); iarg += 2; - } else if (strcmp(arg[iarg],"bikflag") == 0) { - if (iarg+2 > narg) - error->all(FLERR,"Illegal compute snap command"); - bikflag = atoi(arg[iarg+1]); + } else if (strcmp(arg[iarg], "switchinnerflag") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + switchinnerflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); iarg += 2; - } else if (strcmp(arg[iarg],"switchinnerflag") == 0) { - if (iarg+2 > narg) - error->all(FLERR,"Illegal compute snap command"); - switchinnerflag = atoi(arg[iarg+1]); - iarg += 2; - } else if (strcmp(arg[iarg],"sinner") == 0) { + } else if (strcmp(arg[iarg], "sinner") == 0) { iarg++; - if (iarg+ntypes > narg) - error->all(FLERR,"Illegal compute snap command"); - memory->create(sinnerelem,ntypes+1,"snap:sinnerelem"); + if (iarg + ntypes > narg) error->all(FLERR, "Illegal compute {} command", style); + memory->create(sinnerelem, ntypes + 1, "snap:sinnerelem"); for (int i = 0; i < ntypes; i++) - sinnerelem[i+1] = utils::numeric(FLERR,arg[iarg+i],false,lmp); + sinnerelem[i + 1] = utils::numeric(FLERR, arg[iarg + i], false, lmp); sinnerflag = 1; iarg += ntypes; - } else if (strcmp(arg[iarg],"dinner") == 0) { + } else if (strcmp(arg[iarg], "dinner") == 0) { iarg++; - if (iarg+ntypes > narg) - error->all(FLERR,"Illegal compute snap command"); - memory->create(dinnerelem,ntypes+1,"snap:dinnerelem"); + if (iarg + ntypes > narg) error->all(FLERR, "Illegal compute {} command", style); + memory->create(dinnerelem, ntypes + 1, "snap:dinnerelem"); for (int i = 0; i < ntypes; i++) - dinnerelem[i+1] = utils::numeric(FLERR,arg[iarg+i],false,lmp); + dinnerelem[i + 1] = utils::numeric(FLERR, arg[iarg + i], false, lmp); dinnerflag = 1; iarg += ntypes; - } else error->all(FLERR,"Illegal compute snap command"); + } else + error->all(FLERR, "Illegal compute {} command", style); } if (switchinnerflag && !(sinnerflag && dinnerflag)) - error->all(FLERR,"Illegal compute snap command: switchinnerflag = 1, missing sinner/dinner keyword"); + error->all( + FLERR, + "Illegal compute {} command: switchinnerflag = 1, missing sinner/dinner keyword", + style); if (!switchinnerflag && (sinnerflag || dinnerflag)) - error->all(FLERR,"Illegal compute snap command: switchinnerflag = 0, unexpected sinner/dinner keyword"); + error->all( + FLERR, + "Illegal compute {} command: switchinnerflag = 0, unexpected sinner/dinner keyword", + style); - snaptr = new SNA(lmp, rfac0, twojmax, - rmin0, switchflag, bzeroflag, - chemflag, bnormflag, wselfallflag, - nelements, switchinnerflag); + snaptr = new SNA(lmp, rfac0, twojmax, rmin0, switchflag, bzeroflag, chemflag, bnormflag, + wselfallflag, nelements, switchinnerflag); ncoeff = snaptr->ncoeff; - nperdim = ncoeff; - if (quadraticflag) nperdim += (ncoeff*(ncoeff+1))/2; + nvalues = ncoeff; + if (quadraticflag) nvalues += (ncoeff * (ncoeff + 1)) / 2; + + // end code common to all SNAP computes + ndims_force = 3; ndims_virial = 6; - yoffset = nperdim; - zoffset = 2*nperdim; + yoffset = nvalues; + zoffset = 2*nvalues; natoms = atom->natoms; bik_rows = 1; if (bikflag) bik_rows = natoms; size_array_rows = bik_rows+ndims_force*natoms+ndims_virial; - size_array_cols = nperdim*atom->ntypes+1; + size_array_cols = nvalues*atom->ntypes+1; lastcol = size_array_cols-1; ndims_peratom = ndims_force; - size_peratom = ndims_peratom*nperdim*atom->ntypes; + size_peratom = ndims_peratom*nvalues*atom->ntypes; nmax = 0; } @@ -341,8 +338,8 @@ void ComputeSnap::compute_array() const double radi = radelem[itype]; const int* const jlist = firstneigh[i]; const int jnum = numneigh[i]; - const int typeoffset_local = ndims_peratom*nperdim*(itype-1); - const int typeoffset_global = nperdim*(itype-1); + const int typeoffset_local = ndims_peratom*nvalues*(itype-1); + const int typeoffset_global = nvalues*(itype-1); // insure rij, inside, and typej are of size jnum @@ -481,9 +478,9 @@ void ComputeSnap::compute_array() // accumulate bispectrum force contributions to global array for (int itype = 0; itype < atom->ntypes; itype++) { - const int typeoffset_local = ndims_peratom*nperdim*itype; - const int typeoffset_global = nperdim*itype; - for (int icoeff = 0; icoeff < nperdim; icoeff++) { + const int typeoffset_local = ndims_peratom*nvalues*itype; + const int typeoffset_global = nvalues*itype; + for (int icoeff = 0; icoeff < nvalues; icoeff++) { for (int i = 0; i < ntotal; i++) { double *snadi = snap_peratom[i]+typeoffset_local; int iglobal = atom->tag[i]; @@ -549,10 +546,10 @@ void ComputeSnap::dbdotr_compute() int nall = atom->nlocal + atom->nghost; for (int i = 0; i < nall; i++) for (int itype = 0; itype < atom->ntypes; itype++) { - const int typeoffset_local = ndims_peratom*nperdim*itype; - const int typeoffset_global = nperdim*itype; + const int typeoffset_local = ndims_peratom*nvalues*itype; + const int typeoffset_global = nvalues*itype; double *snadi = snap_peratom[i]+typeoffset_local; - for (int icoeff = 0; icoeff < nperdim; icoeff++) { + for (int icoeff = 0; icoeff < nvalues; icoeff++) { double dbdx = snadi[icoeff]; double dbdy = snadi[icoeff+yoffset]; double dbdz = snadi[icoeff+zoffset]; diff --git a/src/ML-SNAP/compute_snap.h b/src/ML-SNAP/compute_snap.h index bc0670e2c7..9f4162d6c4 100644 --- a/src/ML-SNAP/compute_snap.h +++ b/src/ML-SNAP/compute_snap.h @@ -35,7 +35,7 @@ class ComputeSnap : public Compute { private: int natoms, nmax, size_peratom, lastcol; - int ncoeff, nperdim, yoffset, zoffset; + int ncoeff, nvalues, yoffset, zoffset; int ndims_peratom, ndims_force, ndims_virial; double **cutsq; class NeighList *list; diff --git a/src/ML-SNAP/compute_snav_atom.cpp b/src/ML-SNAP/compute_snav_atom.cpp index 0baa4bc110..949dda8e5c 100644 --- a/src/ML-SNAP/compute_snav_atom.cpp +++ b/src/ML-SNAP/compute_snav_atom.cpp @@ -21,6 +21,7 @@ #include "neighbor.h" #include "neigh_list.h" #include "force.h" +#include "pair.h" #include "comm.h" #include "memory.h" #include "error.h" @@ -33,20 +34,22 @@ ComputeSNAVAtom::ComputeSNAVAtom(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg), cutsq(nullptr), list(nullptr), snav(nullptr), radelem(nullptr), wjelem(nullptr), sinnerelem(nullptr), dinnerelem(nullptr) { + + // begin code common to all SNAP computes + double rfac0, rmin0; int twojmax, switchflag, bzeroflag, bnormflag, wselfallflag; int ntypes = atom->ntypes; - int nargmin = 6+2*ntypes; + int nargmin = 6 + 2 * ntypes; - if (narg < nargmin) error->all(FLERR,"Illegal compute snav/atom command"); + if (narg < nargmin) error->all(FLERR, "Illegal compute {} command", style); // default values rmin0 = 0.0; switchflag = 1; bzeroflag = 1; - bnormflag = 0; quadraticflag = 0; chemflag = 0; bnormflag = 0; @@ -56,24 +59,32 @@ ComputeSNAVAtom::ComputeSNAVAtom(LAMMPS *lmp, int narg, char **arg) : // process required arguments - memory->create(radelem,ntypes+1,"snav/atom:radelem"); // offset by 1 to match up with types - memory->create(wjelem,ntypes+1,"snav/atom:wjelem"); - rcutfac = atof(arg[3]); - rfac0 = atof(arg[4]); - twojmax = atoi(arg[5]); + memory->create(radelem, ntypes + 1, "sna/atom:radelem"); // offset by 1 to match up with types + memory->create(wjelem, ntypes + 1, "sna/atom:wjelem"); + + rcutfac = utils::numeric(FLERR, arg[3], false, lmp); + rfac0 = utils::numeric(FLERR, arg[4], false, lmp); + twojmax = utils::inumeric(FLERR, arg[5], false, lmp); + for (int i = 0; i < ntypes; i++) - radelem[i+1] = atof(arg[6+i]); + radelem[i + 1] = + utils::numeric(FLERR, arg[6 + i], false, lmp); for (int i = 0; i < ntypes; i++) - wjelem[i+1] = atof(arg[6+ntypes+i]); + wjelem[i + 1] = + utils::numeric(FLERR, arg[6 + ntypes + i], false, lmp); + // construct cutsq + double cut; - memory->create(cutsq,ntypes+1,ntypes+1,"snav/atom:cutsq"); + cutmax = 0.0; + memory->create(cutsq, ntypes + 1, ntypes + 1, "sna/atom:cutsq"); for (int i = 1; i <= ntypes; i++) { - cut = 2.0*radelem[i]*rcutfac; - cutsq[i][i] = cut*cut; - for (int j = i+1; j <= ntypes; j++) { - cut = (radelem[i]+radelem[j])*rcutfac; - cutsq[i][j] = cutsq[j][i] = cut*cut; + cut = 2.0 * radelem[i] * rcutfac; + if (cut > cutmax) cutmax = cut; + cutsq[i][i] = cut * cut; + for (int j = i + 1; j <= ntypes; j++) { + cut = (radelem[i] + radelem[j]) * rcutfac; + cutsq[i][j] = cutsq[j][i] = cut * cut; } } @@ -87,90 +98,87 @@ ComputeSNAVAtom::ComputeSNAVAtom(LAMMPS *lmp, int narg, char **arg) : int iarg = nargmin; while (iarg < narg) { - if (strcmp(arg[iarg],"rmin0") == 0) { - if (iarg+2 > narg) - error->all(FLERR,"Illegal compute snav/atom command"); - rmin0 = atof(arg[iarg+1]); + if (strcmp(arg[iarg], "rmin0") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + rmin0 = utils::numeric(FLERR, arg[iarg + 1], false, lmp); iarg += 2; - } else if (strcmp(arg[iarg],"switchflag") == 0) { - if (iarg+2 > narg) - error->all(FLERR,"Illegal compute snav/atom command"); - switchflag = atoi(arg[iarg+1]); + } else if (strcmp(arg[iarg], "switchflag") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + switchflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); iarg += 2; - } else if (strcmp(arg[iarg],"bzeroflag") == 0) { - if (iarg+2 > narg) - error->all(FLERR,"Illegal compute snav/atom command"); - bzeroflag = atoi(arg[iarg+1]); + } else if (strcmp(arg[iarg], "bzeroflag") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + bzeroflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); iarg += 2; - } else if (strcmp(arg[iarg],"quadraticflag") == 0) { - if (iarg+2 > narg) - error->all(FLERR,"Illegal compute snav/atom command"); - quadraticflag = atoi(arg[iarg+1]); + } else if (strcmp(arg[iarg], "quadraticflag") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + quadraticflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); iarg += 2; - } else if (strcmp(arg[iarg],"chem") == 0) { - if (iarg+2+ntypes > narg) - error->all(FLERR,"Illegal compute snav/atom command"); + } else if (strcmp(arg[iarg], "chem") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); chemflag = 1; - memory->create(map,ntypes+1,"compute_sna_atom:map"); - nelements = utils::inumeric(FLERR,arg[iarg+1],false,lmp); + memory->create(map, ntypes + 1, "compute_sna_grid:map"); + nelements = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); for (int i = 0; i < ntypes; i++) { - int jelem = utils::inumeric(FLERR,arg[iarg+2+i],false,lmp); - if (jelem < 0 || jelem >= nelements) - error->all(FLERR,"Illegal compute snav/atom command"); - map[i+1] = jelem; + int jelem = utils::inumeric(FLERR, arg[iarg + 2 + i], false, lmp); + if (jelem < 0 || jelem >= nelements) error->all(FLERR, "Illegal compute {} command", style); + map[i + 1] = jelem; } - iarg += 2+ntypes; - } else if (strcmp(arg[iarg],"bnormflag") == 0) { - if (iarg+2 > narg) - error->all(FLERR,"Illegal compute snav/atom command"); - bnormflag = atoi(arg[iarg+1]); + iarg += 2 + ntypes; + } else if (strcmp(arg[iarg], "bnormflag") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + bnormflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); iarg += 2; - } else if (strcmp(arg[iarg],"wselfallflag") == 0) { - if (iarg+2 > narg) - error->all(FLERR,"Illegal compute snav/atom command"); - wselfallflag = atoi(arg[iarg+1]); + } else if (strcmp(arg[iarg], "wselfallflag") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + wselfallflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); iarg += 2; - } else if (strcmp(arg[iarg],"switchinnerflag") == 0) { - if (iarg+2 > narg) - error->all(FLERR,"Illegal compute snav/atom command"); - switchinnerflag = atoi(arg[iarg+1]); + } else if (strcmp(arg[iarg], "switchinnerflag") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style); + switchinnerflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); iarg += 2; - } else if (strcmp(arg[iarg],"sinner") == 0) { + } else if (strcmp(arg[iarg], "sinner") == 0) { iarg++; - if (iarg+ntypes > narg) - error->all(FLERR,"Illegal compute snav/atom command"); - memory->create(sinnerelem,ntypes+1,"snav/atom:sinnerelem"); + if (iarg + ntypes > narg) error->all(FLERR, "Illegal compute {} command", style); + memory->create(sinnerelem, ntypes + 1, "snap:sinnerelem"); for (int i = 0; i < ntypes; i++) - sinnerelem[i+1] = utils::numeric(FLERR,arg[iarg+i],false,lmp); + sinnerelem[i + 1] = utils::numeric(FLERR, arg[iarg + i], false, lmp); sinnerflag = 1; iarg += ntypes; - } else if (strcmp(arg[iarg],"dinner") == 0) { + } else if (strcmp(arg[iarg], "dinner") == 0) { iarg++; - if (iarg+ntypes > narg) - error->all(FLERR,"Illegal compute snav/atom command"); - memory->create(dinnerelem,ntypes+1,"snav/atom:dinnerelem"); + if (iarg + ntypes > narg) error->all(FLERR, "Illegal compute {} command", style); + memory->create(dinnerelem, ntypes + 1, "snap:dinnerelem"); for (int i = 0; i < ntypes; i++) - dinnerelem[i+1] = utils::numeric(FLERR,arg[iarg+i],false,lmp); + dinnerelem[i + 1] = utils::numeric(FLERR, arg[iarg + i], false, lmp); dinnerflag = 1; iarg += ntypes; - } else error->all(FLERR,"Illegal compute snav/atom command"); + } else + error->all(FLERR, "Illegal compute {} command", style); } if (switchinnerflag && !(sinnerflag && dinnerflag)) - error->all(FLERR,"Illegal compute snav/atom command: switchinnerflag = 1, missing sinner/dinner keyword"); + error->all( + FLERR, + "Illegal compute {} command: switchinnerflag = 1, missing sinner/dinner keyword", + style); if (!switchinnerflag && (sinnerflag || dinnerflag)) - error->all(FLERR,"Illegal compute snav/atom command: switchinnerflag = 0, unexpected sinner/dinner keyword"); + error->all( + FLERR, + "Illegal compute {} command: switchinnerflag = 0, unexpected sinner/dinner keyword", + style); - snaptr = new SNA(lmp, rfac0, twojmax, - rmin0, switchflag, bzeroflag, - chemflag, bnormflag, wselfallflag, - nelements, switchinnerflag); + snaptr = new SNA(lmp, rfac0, twojmax, rmin0, switchflag, bzeroflag, chemflag, bnormflag, + wselfallflag, nelements, switchinnerflag); ncoeff = snaptr->ncoeff; - nperdim = ncoeff; - if (quadraticflag) nperdim += (ncoeff*(ncoeff+1))/2; - size_peratom_cols = 6*nperdim*atom->ntypes; + nvalues = ncoeff; + if (quadraticflag) nvalues += (ncoeff * (ncoeff + 1)) / 2; + + // end code common to all SNAP computes + + size_peratom_cols = 6*nvalues*atom->ntypes; comm_reverse = size_peratom_cols; peratom_flag = 1; @@ -203,10 +211,9 @@ void ComputeSNAVAtom::init() { if (force->pair == nullptr) error->all(FLERR,"Compute snav/atom requires a pair style be defined"); - // TODO: Not sure what to do with this error check since cutoff radius is not - // a single number - //if (sqrt(cutsq) > force->pair->cutforce) - // error->all(FLERR,"Compute snav/atom cutoff is longer than pairwise cutoff"); + + if (cutmax > force->pair->cutforce) + error->all(FLERR,"Compute snav/atom cutoff is longer than pairwise cutoff"); // need an occasional full neighbor list @@ -280,7 +287,7 @@ void ComputeSNAVAtom::compute_peratom() const int* const jlist = firstneigh[i]; const int jnum = numneigh[i]; - const int typeoffset = 6*nperdim*(atom->type[i]-1); + const int typeoffset = 6*nvalues*(atom->type[i]-1); // insure rij, inside, and typej are of size jnum @@ -339,17 +346,17 @@ void ComputeSNAVAtom::compute_peratom() for (int icoeff = 0; icoeff < ncoeff; icoeff++) { snavi[icoeff] += snaptr->dblist[icoeff][0]*xtmp; - snavi[icoeff+nperdim] += snaptr->dblist[icoeff][1]*ytmp; - snavi[icoeff+2*nperdim] += snaptr->dblist[icoeff][2]*ztmp; - snavi[icoeff+3*nperdim] += snaptr->dblist[icoeff][1]*ztmp; - snavi[icoeff+4*nperdim] += snaptr->dblist[icoeff][0]*ztmp; - snavi[icoeff+5*nperdim] += snaptr->dblist[icoeff][0]*ytmp; + snavi[icoeff+nvalues] += snaptr->dblist[icoeff][1]*ytmp; + snavi[icoeff+2*nvalues] += snaptr->dblist[icoeff][2]*ztmp; + snavi[icoeff+3*nvalues] += snaptr->dblist[icoeff][1]*ztmp; + snavi[icoeff+4*nvalues] += snaptr->dblist[icoeff][0]*ztmp; + snavi[icoeff+5*nvalues] += snaptr->dblist[icoeff][0]*ytmp; snavj[icoeff] -= snaptr->dblist[icoeff][0]*x[j][0]; - snavj[icoeff+nperdim] -= snaptr->dblist[icoeff][1]*x[j][1]; - snavj[icoeff+2*nperdim] -= snaptr->dblist[icoeff][2]*x[j][2]; - snavj[icoeff+3*nperdim] -= snaptr->dblist[icoeff][1]*x[j][2]; - snavj[icoeff+4*nperdim] -= snaptr->dblist[icoeff][0]*x[j][2]; - snavj[icoeff+5*nperdim] -= snaptr->dblist[icoeff][0]*x[j][1]; + snavj[icoeff+nvalues] -= snaptr->dblist[icoeff][1]*x[j][1]; + snavj[icoeff+2*nvalues] -= snaptr->dblist[icoeff][2]*x[j][2]; + snavj[icoeff+3*nvalues] -= snaptr->dblist[icoeff][1]*x[j][2]; + snavj[icoeff+4*nvalues] -= snaptr->dblist[icoeff][0]*x[j][2]; + snavj[icoeff+5*nvalues] -= snaptr->dblist[icoeff][0]*x[j][1]; } if (quadraticflag) { @@ -369,17 +376,17 @@ void ComputeSNAVAtom::compute_peratom() double dbytmp = bi*biy; double dbztmp = bi*biz; snavi[ncount] += dbxtmp*xtmp; - snavi[ncount+nperdim] += dbytmp*ytmp; - snavi[ncount+2*nperdim] += dbztmp*ztmp; - snavi[ncount+3*nperdim] += dbytmp*ztmp; - snavi[ncount+4*nperdim] += dbxtmp*ztmp; - snavi[ncount+5*nperdim] += dbxtmp*ytmp; + snavi[ncount+nvalues] += dbytmp*ytmp; + snavi[ncount+2*nvalues] += dbztmp*ztmp; + snavi[ncount+3*nvalues] += dbytmp*ztmp; + snavi[ncount+4*nvalues] += dbxtmp*ztmp; + snavi[ncount+5*nvalues] += dbxtmp*ytmp; snavj[ncount] -= dbxtmp*x[j][0]; - snavj[ncount+nperdim] -= dbytmp*x[j][1]; - snavj[ncount+2*nperdim] -= dbztmp*x[j][2]; - snavj[ncount+3*nperdim] -= dbytmp*x[j][2]; - snavj[ncount+4*nperdim] -= dbxtmp*x[j][2]; - snavj[ncount+5*nperdim] -= dbxtmp*x[j][1]; + snavj[ncount+nvalues] -= dbytmp*x[j][1]; + snavj[ncount+2*nvalues] -= dbztmp*x[j][2]; + snavj[ncount+3*nvalues] -= dbytmp*x[j][2]; + snavj[ncount+4*nvalues] -= dbxtmp*x[j][2]; + snavj[ncount+5*nvalues] -= dbxtmp*x[j][1]; ncount++; // upper-triangular elements of quadratic matrix @@ -392,17 +399,17 @@ void ComputeSNAVAtom::compute_peratom() double dbztmp = bi*snaptr->dblist[jcoeff][2] + biz*snaptr->blist[jcoeff]; snavi[ncount] += dbxtmp*xtmp; - snavi[ncount+nperdim] += dbytmp*ytmp; - snavi[ncount+2*nperdim] += dbztmp*ztmp; - snavi[ncount+3*nperdim] += dbytmp*ztmp; - snavi[ncount+4*nperdim] += dbxtmp*ztmp; - snavi[ncount+5*nperdim] += dbxtmp*ytmp; + snavi[ncount+nvalues] += dbytmp*ytmp; + snavi[ncount+2*nvalues] += dbztmp*ztmp; + snavi[ncount+3*nvalues] += dbytmp*ztmp; + snavi[ncount+4*nvalues] += dbxtmp*ztmp; + snavi[ncount+5*nvalues] += dbxtmp*ytmp; snavj[ncount] -= dbxtmp*x[j][0]; - snavj[ncount+nperdim] -= dbytmp*x[j][1]; - snavj[ncount+2*nperdim] -= dbztmp*x[j][2]; - snavj[ncount+3*nperdim] -= dbytmp*x[j][2]; - snavj[ncount+4*nperdim] -= dbxtmp*x[j][2]; - snavj[ncount+5*nperdim] -= dbxtmp*x[j][1]; + snavj[ncount+nvalues] -= dbytmp*x[j][1]; + snavj[ncount+2*nvalues] -= dbztmp*x[j][2]; + snavj[ncount+3*nvalues] -= dbytmp*x[j][2]; + snavj[ncount+4*nvalues] -= dbxtmp*x[j][2]; + snavj[ncount+5*nvalues] -= dbxtmp*x[j][1]; ncount++; } } diff --git a/src/ML-SNAP/compute_snav_atom.h b/src/ML-SNAP/compute_snav_atom.h index 1eb84d8df7..2eac0b7b28 100644 --- a/src/ML-SNAP/compute_snav_atom.h +++ b/src/ML-SNAP/compute_snav_atom.h @@ -37,7 +37,7 @@ class ComputeSNAVAtom : public Compute { private: int nmax; - int ncoeff, nperdim; + int ncoeff, nvalues; double **cutsq; class NeighList *list; double **snav; @@ -50,6 +50,7 @@ class ComputeSNAVAtom : public Compute { double *sinnerelem; double *dinnerelem; class SNA *snaptr; + double cutmax; int quadraticflag; };