diff --git a/doc/src/compute_mliap.rst b/doc/src/compute_mliap.rst index 5d96de636a..296d3b009a 100644 --- a/doc/src/compute_mliap.rst +++ b/doc/src/compute_mliap.rst @@ -17,10 +17,8 @@ Syntax .. parsed-literal:: - *model* values = style Nelems Nparams + *model* values = style style = *linear* or *quadratic* - Nelems = number of elements - Nparams = number of parameters per element *descriptor* values = style filename style = *sna* filename = name of file containing descriptor definitions @@ -31,7 +29,7 @@ Examples .. code-block:: LAMMPS - compute mliap model linear 2 31 descriptor sna Ta06A.mliap.descriptor + compute mliap model linear descriptor sna Ta06A.mliap.descriptor Description """"""""""" @@ -58,14 +56,8 @@ The compute *mliap* command must be followed by two keywords *model* and *descriptor* in either order. The *model* keyword is followed by a model style, currently limited to -either *linear* or *quadratic*. In both cases, -this is followed by two arguments. *nelems* is the number of elements. -It must be equal to the number of LAMMPS atom types. *nparams* -is the number of parameters per element for this model i.e. -the number of parameter gradients for each element. Note these definitions -are identical to those of *nelems* and *nparams* in the -:doc:`pair_style mliap ` model file. - +either *linear* or *quadratic*. + The *descriptor* keyword is followed by a descriptor style, and additional arguments. Currently the only descriptor style is *sna*, indicating the bispectrum component descriptors used by the Spectral Neighbor Analysis Potential (SNAP) potentials of diff --git a/examples/mliap/in.mliap.quadratic.compute b/examples/mliap/in.mliap.quadratic.compute index 6db1f839c5..929dbf3824 100644 --- a/examples/mliap/in.mliap.quadratic.compute +++ b/examples/mliap/in.mliap.quadratic.compute @@ -78,7 +78,7 @@ thermo 100 # run compute mliap with gradgradflag = 1 -compute snap all mliap descriptor sna compute.mliap.descriptor model quadratic 2 21 gradgradflag 1 +compute snap all mliap descriptor sna compute.mliap.descriptor model quadratic gradgradflag 1 fix snap all ave/time 1 1 1 c_snap[*] file compute.quadratic.gg1.dat mode vector thermo_style custom & @@ -93,7 +93,7 @@ unfix snap # run compute mliap with gradgradflag = 0 -compute snap all mliap descriptor sna compute.mliap.descriptor model quadratic 2 21 gradgradflag 0 +compute snap all mliap descriptor sna compute.mliap.descriptor model quadratic gradgradflag 0 fix snap all ave/time 1 1 1 c_snap[*] file compute.quadratic.gg0.dat mode vector thermo_style custom & diff --git a/examples/mliap/in.mliap.snap.compute b/examples/mliap/in.mliap.snap.compute index 1fb00d90cf..4cfccedbdf 100644 --- a/examples/mliap/in.mliap.snap.compute +++ b/examples/mliap/in.mliap.snap.compute @@ -79,7 +79,7 @@ thermo 100 # run compute mliap with gradgradflag = 1 -compute snap all mliap descriptor sna compute.mliap.descriptor model linear 2 6 gradgradflag 1 +compute snap all mliap descriptor sna compute.mliap.descriptor model linear gradgradflag 1 fix snap all ave/time 1 1 1 c_snap[*] file compute.snap.gg1.dat mode vector thermo_style custom & @@ -94,7 +94,7 @@ unfix snap # run compute mliap with gradgradflag = 0 -compute snap all mliap descriptor sna compute.mliap.descriptor model linear 2 6 gradgradflag 0 +compute snap all mliap descriptor sna compute.mliap.descriptor model linear gradgradflag 0 fix snap all ave/time 1 1 1 c_snap[*] file compute.snap.gg0.dat mode vector thermo_style custom & diff --git a/examples/mliap/log.03Jul20.mliap.quadratic.compute.g++.1 b/examples/mliap/log.03Jul20.mliap.quadratic.compute.g++.1 index dd83662b9e..f4dcb93ff4 100644 --- a/examples/mliap/log.03Jul20.mliap.quadratic.compute.g++.1 +++ b/examples/mliap/log.03Jul20.mliap.quadratic.compute.g++.1 @@ -110,7 +110,7 @@ thermo 100 # run compute mliap with gradgradflag = 1 -compute snap all mliap descriptor sna compute.mliap.descriptor model quadratic 2 21 gradgradflag 1 +compute snap all mliap descriptor sna compute.mliap.descriptor model quadratic gradgradflag 1 SNAP keyword rcutfac 1.0 SNAP keyword twojmax 2 SNAP keyword nelems 2 @@ -160,7 +160,7 @@ Neighbor list info ... pair build: full/bin/atomonly stencil: full/bin/3d bin: standard -Per MPI rank memory allocation (min/avg/max) = 22.45 | 22.45 | 22.45 Mbytes +Per MPI rank memory allocation (min/avg/max) = 22.47 | 22.47 | 22.47 Mbytes PotEng Pxy c_bsum2[20] c_vbsum[220] v_db_2_100 c_snap[1][43] c_snap[13][43] c_snap[1][42] c_snap[12][42] c_snap[6][42] 322.86952 1505558.1 4.2492771e+08 -4952473 28484035 322.86952 1505558.1 4.2492771e+08 -4952473 28484035 Loop time of 1e-06 on 1 procs for 0 steps with 2 atoms @@ -196,7 +196,7 @@ unfix snap # run compute mliap with gradgradflag = 0 -compute snap all mliap descriptor sna compute.mliap.descriptor model quadratic 2 21 gradgradflag 0 +compute snap all mliap descriptor sna compute.mliap.descriptor model quadratic gradgradflag 0 SNAP keyword rcutfac 1.0 SNAP keyword twojmax 2 SNAP keyword nelems 2 @@ -247,12 +247,12 @@ Neighbor list info ... pair build: full/bin/atomonly stencil: full/bin/3d bin: standard -Per MPI rank memory allocation (min/avg/max) = 70.2 | 70.2 | 70.2 Mbytes +Per MPI rank memory allocation (min/avg/max) = 70.24 | 70.24 | 70.24 Mbytes PotEng Pxy c_bsum2[20] c_vbsum[220] v_db_2_100 c_snap[1][43] c_snap[13][43] c_snap[1][42] c_snap[12][42] c_snap[6][42] 322.86952 1505558.1 4.2492771e+08 -4952473 28484035 322.86952 1505558.1 4.2492771e+08 -4952473 28484035 -Loop time of 1e-06 on 1 procs for 0 steps with 2 atoms +Loop time of 2e-06 on 1 procs for 0 steps with 2 atoms -200.0% CPU use with 1 MPI tasks x no OpenMP threads +50.0% CPU use with 1 MPI tasks x no OpenMP threads MPI task timing breakdown: Section | min time | avg time | max time |%varavg| %total @@ -262,7 +262,7 @@ 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 | | 1e-06 | | |100.00 +Other | | 2e-06 | | |100.00 Nlocal: 2.0 ave 2.0 max 2.0 min Histogram: 1 0 0 0 0 0 0 0 0 0 diff --git a/examples/mliap/log.03Jul20.mliap.quadratic.compute.g++.4 b/examples/mliap/log.03Jul20.mliap.quadratic.compute.g++.4 index 5c5787559f..4900b8133d 100644 --- a/examples/mliap/log.03Jul20.mliap.quadratic.compute.g++.4 +++ b/examples/mliap/log.03Jul20.mliap.quadratic.compute.g++.4 @@ -110,7 +110,7 @@ thermo 100 # run compute mliap with gradgradflag = 1 -compute snap all mliap descriptor sna compute.mliap.descriptor model quadratic 2 21 gradgradflag 1 +compute snap all mliap descriptor sna compute.mliap.descriptor model quadratic gradgradflag 1 SNAP keyword rcutfac 1.0 SNAP keyword twojmax 2 SNAP keyword nelems 2 @@ -161,12 +161,12 @@ Neighbor list info ... stencil: full/bin/3d bin: standard WARNING: Proc sub-domain size < neighbor skin, could lead to lost atoms (../domain.cpp:964) -Per MPI rank memory allocation (min/avg/max) = 22.18 | 22.18 | 22.18 Mbytes +Per MPI rank memory allocation (min/avg/max) = 6.429 | 6.43 | 6.432 Mbytes PotEng Pxy c_bsum2[20] c_vbsum[220] v_db_2_100 c_snap[1][43] c_snap[13][43] c_snap[1][42] c_snap[12][42] c_snap[6][42] 322.86952 1505558.1 4.2492771e+08 -4952473 28484035 322.86952 1505558.1 4.2492771e+08 -4952473 28484035 -Loop time of 3.5e-06 on 4 procs for 0 steps with 2 atoms +Loop time of 4.25e-06 on 4 procs for 0 steps with 2 atoms -92.9% CPU use with 4 MPI tasks x no OpenMP threads +82.4% CPU use with 4 MPI tasks x no OpenMP threads MPI task timing breakdown: Section | min time | avg time | max time |%varavg| %total @@ -176,7 +176,7 @@ 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 | | 3.5e-06 | | |100.00 +Other | | 4.25e-06 | | |100.00 Nlocal: 0.5 ave 1.0 max 0.0 min Histogram: 2 0 0 0 0 0 0 0 0 2 @@ -197,7 +197,7 @@ unfix snap # run compute mliap with gradgradflag = 0 -compute snap all mliap descriptor sna compute.mliap.descriptor model quadratic 2 21 gradgradflag 0 +compute snap all mliap descriptor sna compute.mliap.descriptor model quadratic gradgradflag 0 SNAP keyword rcutfac 1.0 SNAP keyword twojmax 2 SNAP keyword nelems 2 @@ -249,12 +249,12 @@ Neighbor list info ... stencil: full/bin/3d bin: standard WARNING: Proc sub-domain size < neighbor skin, could lead to lost atoms (../domain.cpp:964) -Per MPI rank memory allocation (min/avg/max) = 69.68 | 69.81 | 69.93 Mbytes +Per MPI rank memory allocation (min/avg/max) = 53.93 | 54.06 | 54.18 Mbytes PotEng Pxy c_bsum2[20] c_vbsum[220] v_db_2_100 c_snap[1][43] c_snap[13][43] c_snap[1][42] c_snap[12][42] c_snap[6][42] 322.86952 1505558.1 4.2492771e+08 -4952473 28484035 322.86952 1505558.1 4.2492771e+08 -4952473 28484035 -Loop time of 2.75e-06 on 4 procs for 0 steps with 2 atoms +Loop time of 4e-06 on 4 procs for 0 steps with 2 atoms -100.0% CPU use with 4 MPI tasks x no OpenMP threads +106.2% CPU use with 4 MPI tasks x no OpenMP threads MPI task timing breakdown: Section | min time | avg time | max time |%varavg| %total @@ -264,7 +264,7 @@ 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.75e-06 | | |100.00 +Other | | 4e-06 | | |100.00 Nlocal: 0.5 ave 1.0 max 0.0 min Histogram: 2 0 0 0 0 0 0 0 0 2 diff --git a/examples/mliap/log.03Jul20.mliap.snap.compute.g++.1 b/examples/mliap/log.03Jul20.mliap.snap.compute.g++.1 index d314c9c5b5..222479baf3 100644 --- a/examples/mliap/log.03Jul20.mliap.snap.compute.g++.1 +++ b/examples/mliap/log.03Jul20.mliap.snap.compute.g++.1 @@ -111,7 +111,7 @@ thermo 100 # run compute mliap with gradgradflag = 1 -compute snap all mliap descriptor sna compute.mliap.descriptor model linear 2 6 gradgradflag 1 +compute snap all mliap descriptor sna compute.mliap.descriptor model linear gradgradflag 1 SNAP keyword rcutfac 1.0 SNAP keyword twojmax 2 SNAP keyword nelems 2 @@ -161,12 +161,12 @@ Neighbor list info ... pair build: full/bin/atomonly stencil: full/bin/3d bin: standard -Per MPI rank memory allocation (min/avg/max) = 10.73 | 10.73 | 10.73 Mbytes +Per MPI rank memory allocation (min/avg/max) = 10.75 | 10.75 | 10.75 Mbytes PotEng Pxy c_bsum2[1] c_vbsum[55] v_db_2_25 c_snap[1][13] c_snap[13][13] c_snap[1][8] c_snap[12][12] c_snap[6][12] 322.86952 1505558.1 364182.88 -240.25066 1381.7961 322.86952 1505558.1 364182.88 -240.25066 1381.7961 -Loop time of 2e-06 on 1 procs for 0 steps with 2 atoms +Loop time of 0 on 1 procs for 0 steps with 2 atoms -150.0% CPU use with 1 MPI tasks x no OpenMP threads +0.0% CPU use with 1 MPI tasks x no OpenMP threads MPI task timing breakdown: Section | min time | avg time | max time |%varavg| %total @@ -176,7 +176,7 @@ 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 | | 2e-06 | | |100.00 +Other | | 0 | | | 0.00 Nlocal: 2.0 ave 2.0 max 2.0 min Histogram: 1 0 0 0 0 0 0 0 0 0 @@ -197,7 +197,7 @@ unfix snap # run compute mliap with gradgradflag = 0 -compute snap all mliap descriptor sna compute.mliap.descriptor model linear 2 6 gradgradflag 0 +compute snap all mliap descriptor sna compute.mliap.descriptor model linear gradgradflag 0 SNAP keyword rcutfac 1.0 SNAP keyword twojmax 2 SNAP keyword nelems 2 @@ -248,7 +248,7 @@ Neighbor list info ... pair build: full/bin/atomonly stencil: full/bin/3d bin: standard -Per MPI rank memory allocation (min/avg/max) = 22.85 | 22.85 | 22.85 Mbytes +Per MPI rank memory allocation (min/avg/max) = 22.89 | 22.89 | 22.89 Mbytes PotEng Pxy c_bsum2[1] c_vbsum[55] v_db_2_25 c_snap[1][13] c_snap[13][13] c_snap[1][8] c_snap[12][12] c_snap[6][12] 322.86952 1505558.1 364182.88 -240.25066 1381.7961 322.86952 1505558.1 364182.88 -240.25066 1381.7961 Loop time of 1e-06 on 1 procs for 0 steps with 2 atoms diff --git a/examples/mliap/log.03Jul20.mliap.snap.compute.g++.4 b/examples/mliap/log.03Jul20.mliap.snap.compute.g++.4 index 06673e719f..9f4aa7897b 100644 --- a/examples/mliap/log.03Jul20.mliap.snap.compute.g++.4 +++ b/examples/mliap/log.03Jul20.mliap.snap.compute.g++.4 @@ -111,7 +111,7 @@ thermo 100 # run compute mliap with gradgradflag = 1 -compute snap all mliap descriptor sna compute.mliap.descriptor model linear 2 6 gradgradflag 1 +compute snap all mliap descriptor sna compute.mliap.descriptor model linear gradgradflag 1 SNAP keyword rcutfac 1.0 SNAP keyword twojmax 2 SNAP keyword nelems 2 @@ -162,12 +162,12 @@ Neighbor list info ... stencil: full/bin/3d bin: standard WARNING: Proc sub-domain size < neighbor skin, could lead to lost atoms (../domain.cpp:964) -Per MPI rank memory allocation (min/avg/max) = 10.69 | 10.69 | 10.7 Mbytes +Per MPI rank memory allocation (min/avg/max) = 10.7 | 10.7 | 10.71 Mbytes PotEng Pxy c_bsum2[1] c_vbsum[55] v_db_2_25 c_snap[1][13] c_snap[13][13] c_snap[1][8] c_snap[12][12] c_snap[6][12] 322.86952 1505558.1 364182.88 -240.25066 1381.7961 322.86952 1505558.1 364182.88 -240.25066 1381.7961 -Loop time of 5.5e-06 on 4 procs for 0 steps with 2 atoms +Loop time of 3.75e-06 on 4 procs for 0 steps with 2 atoms -104.5% CPU use with 4 MPI tasks x no OpenMP threads +86.7% CPU use with 4 MPI tasks x no OpenMP threads MPI task timing breakdown: Section | min time | avg time | max time |%varavg| %total @@ -177,7 +177,7 @@ 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 | | 5.5e-06 | | |100.00 +Other | | 3.75e-06 | | |100.00 Nlocal: 0.5 ave 1.0 max 0.0 min Histogram: 2 0 0 0 0 0 0 0 0 2 @@ -198,7 +198,7 @@ unfix snap # run compute mliap with gradgradflag = 0 -compute snap all mliap descriptor sna compute.mliap.descriptor model linear 2 6 gradgradflag 0 +compute snap all mliap descriptor sna compute.mliap.descriptor model linear gradgradflag 0 SNAP keyword rcutfac 1.0 SNAP keyword twojmax 2 SNAP keyword nelems 2 @@ -250,12 +250,12 @@ Neighbor list info ... stencil: full/bin/3d bin: standard WARNING: Proc sub-domain size < neighbor skin, could lead to lost atoms (../domain.cpp:964) -Per MPI rank memory allocation (min/avg/max) = 22.57 | 22.69 | 22.82 Mbytes +Per MPI rank memory allocation (min/avg/max) = 22.57 | 22.71 | 22.85 Mbytes PotEng Pxy c_bsum2[1] c_vbsum[55] v_db_2_25 c_snap[1][13] c_snap[13][13] c_snap[1][8] c_snap[12][12] c_snap[6][12] 322.86952 1505558.1 364182.88 -240.25066 1381.7961 322.86952 1505558.1 364182.88 -240.25066 1381.7961 -Loop time of 2.5e-06 on 4 procs for 0 steps with 2 atoms +Loop time of 2e-06 on 4 procs for 0 steps with 2 atoms -120.0% CPU use with 4 MPI tasks x no OpenMP threads +100.0% CPU use with 4 MPI tasks x no OpenMP threads MPI task timing breakdown: Section | min time | avg time | max time |%varavg| %total @@ -265,7 +265,7 @@ 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.5e-06 | | |100.00 +Other | | 2e-06 | | |100.00 Nlocal: 0.5 ave 1.0 max 0.0 min Histogram: 2 0 0 0 0 0 0 0 0 2 diff --git a/examples/mliap/log.21Jun20.mliap.snap.WBe.PRB201.g++.1 b/examples/mliap/log.21Jun20.mliap.snap.WBe.PRB2019.g++.1 similarity index 100% rename from examples/mliap/log.21Jun20.mliap.snap.WBe.PRB201.g++.1 rename to examples/mliap/log.21Jun20.mliap.snap.WBe.PRB2019.g++.1 diff --git a/examples/mliap/log.21Jun20.mliap.snap.WBe.PRB201.g++.4 b/examples/mliap/log.21Jun20.mliap.snap.WBe.PRB2019.g++.4 similarity index 100% rename from examples/mliap/log.21Jun20.mliap.snap.WBe.PRB201.g++.4 rename to examples/mliap/log.21Jun20.mliap.snap.WBe.PRB2019.g++.4 diff --git a/src/MLIAP/compute_mliap.cpp b/src/MLIAP/compute_mliap.cpp index 47d270d071..ca3f473303 100644 --- a/src/MLIAP/compute_mliap.cpp +++ b/src/MLIAP/compute_mliap.cpp @@ -13,7 +13,7 @@ #include #include -#include "mliap.h" +#include "mliap_data.h" #include "mliap_model_linear.h" #include "mliap_model_quadratic.h" #include "mliap_descriptor_snap.h" @@ -61,17 +61,11 @@ ComputeMLIAP::ComputeMLIAP(LAMMPS *lmp, int narg, char **arg) : if (strcmp(arg[iarg],"model") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal compute mliap command"); if (strcmp(arg[iarg+1],"linear") == 0) { - if (iarg+4 > narg) error->all(FLERR,"Illegal compute mliap command"); - int ntmp1 = atoi(arg[iarg+2]); - int ntmp2 = atoi(arg[iarg+3]); - model = new MLIAPModelLinear(lmp,ntmp1,ntmp2); - iarg += 4; + model = new MLIAPModelLinear(lmp); + iarg += 2; } else if (strcmp(arg[iarg+1],"quadratic") == 0) { - if (iarg+4 > narg) error->all(FLERR,"Illegal compute mliap command"); - int ntmp1 = atoi(arg[iarg+2]); - int ntmp2 = atoi(arg[iarg+3]); - model = new MLIAPModelQuadratic(lmp,ntmp1,ntmp2); - iarg += 4; + model = new MLIAPModelQuadratic(lmp); + iarg += 2; } else error->all(FLERR,"Illegal compute mliap command"); modelflag = 1; } else if (strcmp(arg[iarg],"descriptor") == 0) { @@ -95,9 +89,10 @@ ComputeMLIAP::ComputeMLIAP(LAMMPS *lmp, int narg, char **arg) : if (modelflag == 0 || descriptorflag == 0) error->all(FLERR,"Illegal compute_style command"); - ndescriptors = descriptor->ndescriptors; - nparams = model->nparams; - nelements = model->nelements; + // need to tell model how many descriptors + // so it can figure out how many parameters + + model->set_ndescriptors(descriptor->ndescriptors); // create a minimal map, placeholder for more general map @@ -106,14 +101,11 @@ ComputeMLIAP::ComputeMLIAP(LAMMPS *lmp, int narg, char **arg) : for (int i = 1; i <= atom->ntypes; i++) map[i] = i-1; - mliap = new MLIAP(lmp, ndescriptors, nparams, nelements, gradgradflag, map, model, descriptor); + data = new MLIAPData(lmp, gradgradflag, map, model, descriptor); - size_array_rows = mliap->size_array_rows; - size_array_cols = mliap->size_array_cols; + size_array_rows = data->size_array_rows; + size_array_cols = data->size_array_cols; lastcol = size_array_cols-1; - nmax = 0; - natomgamma_max = 0; - printf("Made it to here\n"); } /* ---------------------------------------------------------------------- */ @@ -125,10 +117,9 @@ ComputeMLIAP::~ComputeMLIAP() memory->destroy(mliaparray); memory->destroy(mliaparrayall); - memory->destroy(mliap->gradforce); memory->destroy(map); - delete mliap; + delete data; delete model; delete descriptor; } @@ -162,23 +153,19 @@ void ComputeMLIAP::init() model->init(); descriptor->init(); - mliap->init(); + data->init(); // consistency checks - if (descriptor->ndescriptors != model->ndescriptors) - error->all(FLERR,"Incompatible model and descriptor definitions"); - if (descriptor->nelements != model->nelements) - error->all(FLERR,"Incompatible model and descriptor definitions"); - if (nelements != atom->ntypes) + if (data->nelements != atom->ntypes) error->all(FLERR,"nelements must equal ntypes"); // allocate memory for global array memory->create(mliaparray,size_array_rows,size_array_cols, - "mliap:mliaparray"); + "compute_mliap:mliaparray"); memory->create(mliaparrayall,size_array_rows,size_array_cols, - "mliap:mliaparrayall"); + "compute_mliap:mliaparrayall"); array = mliaparrayall; // find compute for reference energy @@ -223,86 +210,51 @@ void ComputeMLIAP::compute_array() for (int jcol = 0; jcol < size_array_cols; jcol++) mliaparray[irow][jcol] = 0.0; - // grow nmax gradforce array if necessary - - if (atom->nmax > nmax) { - memory->destroy(mliap->gradforce); - nmax = atom->nmax; - memory->create(mliap->gradforce,nmax,mliap->size_gradforce, - "mliap:gradforce"); - } - - // clear gradforce array - - for (int i = 0; i < ntotal; i++) - for (int j = 0; j < mliap->size_gradforce; j++) { - mliap->gradforce[i][j] = 0.0; - } - // invoke full neighbor list (will copy or build if necessary) neighbor->build_one(list); - mliap->generate_neigharrays(list); + data->generate_neighdata(list); // compute descriptors - descriptor->compute_descriptors(mliap->natomdesc, mliap->iatommliap, mliap->ielemmliap, mliap->numneighmliap, - mliap->jatommliap, mliap->jelemmliap, mliap->descriptors); + descriptor->compute_descriptors(data); - if (gradgradflag) { - - // grow gamma arrays if necessary - - const int natomgamma = list->inum; - if (natomgamma_max < natomgamma) { - memory->grow(mliap->gamma_row_index,natomgamma,mliap->gamma_nnz,"ComputeMLIAP:gamma_row_index"); - memory->grow(mliap->gamma_col_index,natomgamma,mliap->gamma_nnz,"ComputeMLIAP:gamma_col_index"); - memory->grow(mliap->gamma,natomgamma,mliap->gamma_nnz,"ComputeMLIAP:gamma"); - natomgamma_max = natomgamma; - } + if (gradgradflag == 1) { // calculate double gradient w.r.t. parameters and descriptors - model->param_gradient(mliap->natomdesc, mliap->iatommliap, mliap->ielemmliap, mliap->descriptors, mliap->gamma_row_index, - mliap->gamma_col_index, mliap->gamma, mliap->egradient); + model->compute_gradgrads(data); - // calculate descriptor gradient contributions to parameter gradients + // calculate gradients of forces w.r.t. parameters - descriptor->compute_gradients(mliap->natomdesc, mliap->iatommliap, mliap->ielemmliap, mliap->numneighmliap, - mliap->jatommliap, mliap->jelemmliap, - mliap->gamma_nnz, mliap->gamma_row_index, - mliap->gamma_col_index, mliap->gamma, mliap->gradforce, - mliap->yoffset, mliap->zoffset); + descriptor->compute_force_gradients(data); - } else { + } else if (gradgradflag == 0) { // calculate descriptor gradients - descriptor->compute_descriptor_gradients(mliap->natomdesc, mliap->iatommliap, mliap->ielemmliap, mliap->numneighmliap, - mliap->jatommliap, mliap->jelemmliap, mliap->graddesc); + descriptor->compute_descriptor_gradients(data); - // calculate force gradients w.r.t. parameters + // calculate gradients of forces w.r.t. parameters - model->compute_force_gradients(mliap->descriptors, mliap->natomdesc, mliap->iatommliap, mliap->ielemmliap, - mliap->numneighmliap, mliap->jatommliap, mliap->jelemmliap, mliap->graddesc, - mliap->yoffset, mliap->zoffset, mliap->gradforce, mliap->egradient); - - } + model->compute_force_gradients(data); + + } else error->all(FLERR,"Invalid value for gradgradflag"); // accumulate descriptor gradient contributions to global array - for (int ielem = 0; ielem < nelements; ielem++) { - const int elemoffset = nparams*ielem; - for (int jparam = 0; jparam < nparams; jparam++) { + for (int ielem = 0; ielem < data->nelements; ielem++) { + const int elemoffset = data->nparams*ielem; + for (int jparam = 0; jparam < data->nparams; jparam++) { int irow = 1; for (int i = 0; i < ntotal; i++) { - double *gradforcei = mliap->gradforce[i]+elemoffset; + double *gradforcei = data->gradforce[i]+elemoffset; int iglobal = atom->tag[i]; int irow = 3*(iglobal-1)+1; mliaparray[irow][jparam+elemoffset] += gradforcei[jparam]; - mliaparray[irow+1][jparam+elemoffset] += gradforcei[jparam+mliap->yoffset]; - mliaparray[irow+2][jparam+elemoffset] += gradforcei[jparam+mliap->zoffset]; + mliaparray[irow+1][jparam+elemoffset] += gradforcei[jparam+data->yoffset]; + mliaparray[irow+2][jparam+elemoffset] += gradforcei[jparam+data->zoffset]; } } } @@ -323,10 +275,10 @@ void ComputeMLIAP::compute_array() // copy energy gradient contributions to global array - for (int ielem = 0; ielem < nelements; ielem++) { - const int elemoffset = nparams*ielem; - for (int jparam = 0; jparam < nparams; jparam++) - mliaparray[0][jparam+elemoffset] = mliap->egradient[jparam+elemoffset]; + for (int ielem = 0; ielem < data->nelements; ielem++) { + const int elemoffset = data->nparams*ielem; + for (int jparam = 0; jparam < data->nparams; jparam++) + mliaparray[0][jparam+elemoffset] = data->egradient[jparam+elemoffset]; } // sum up over all processes @@ -343,7 +295,7 @@ void ComputeMLIAP::compute_array() // switch to Voigt notation c_virial->compute_vector(); - irow += 3*mliap->natoms; + irow += 3*data->natoms_array; mliaparrayall[irow++][lastcol] = c_virial->vector[0]; mliaparrayall[irow++][lastcol] = c_virial->vector[1]; mliaparrayall[irow++][lastcol] = c_virial->vector[2]; @@ -361,20 +313,20 @@ void ComputeMLIAP::compute_array() void ComputeMLIAP::dbdotr_compute() { double **x = atom->x; - int irow0 = 1+mliap->ndims_force*mliap->natoms; + int irow0 = 1+data->ndims_force*data->natoms_array; // sum over bispectrum contributions to forces // on all particles including ghosts int nall = atom->nlocal + atom->nghost; for (int i = 0; i < nall; i++) - for (int ielem = 0; ielem < nelements; ielem++) { - const int elemoffset = nparams*ielem; - double *gradforcei = mliap->gradforce[i]+elemoffset; - for (int jparam = 0; jparam < nparams; jparam++) { + for (int ielem = 0; ielem < data->nelements; ielem++) { + const int elemoffset = data->nparams*ielem; + double *gradforcei = data->gradforce[i]+elemoffset; + for (int jparam = 0; jparam < data->nparams; jparam++) { double dbdx = gradforcei[jparam]; - double dbdy = gradforcei[jparam+mliap->yoffset]; - double dbdz = gradforcei[jparam+mliap->zoffset]; + double dbdy = gradforcei[jparam+data->yoffset]; + double dbdz = gradforcei[jparam+data->zoffset]; int irow = irow0; mliaparray[irow++][jparam+elemoffset] += dbdx*x[i][0]; mliaparray[irow++][jparam+elemoffset] += dbdy*x[i][1]; @@ -394,12 +346,15 @@ double ComputeMLIAP::memory_usage() { double bytes = size_array_rows*size_array_cols * - sizeof(double); // mliaparray + sizeof(double); // mliaparray bytes += size_array_rows*size_array_cols * - sizeof(double); // mliaparrayall - bytes += nmax*mliap->size_gradforce * sizeof(double); // gradforce + sizeof(double); // mliaparrayall int n = atom->ntypes+1; - bytes += n*sizeof(int); // map + bytes += n*sizeof(int); // map + + bytes += descriptor->memory_usage(); // Descriptor object + bytes += model->memory_usage(); // Model object + bytes += data->memory_usage(); // Data object return bytes; } diff --git a/src/MLIAP/compute_mliap.h b/src/MLIAP/compute_mliap.h index 9c3f7d8054..f1ee91d572 100644 --- a/src/MLIAP/compute_mliap.h +++ b/src/MLIAP/compute_mliap.h @@ -39,17 +39,14 @@ class ComputeMLIAP : public Compute { private: double **mliaparray, **mliaparrayall; class NeighList *list; - double **gradforce; int *map; // map types to [0,nelements) int ndescriptors; // number of descriptors int nparams; // number of model parameters per element int nelements; int gradgradflag; // 1 for graddesc, 0 for gamma - int nmax; - int natomgamma_max; // allocated size of atom neighbor arrays class MLIAPModel *model; class MLIAPDescriptor *descriptor; - class MLIAP *mliap; + class MLIAPData *data; Compute *c_pe; Compute *c_virial; diff --git a/src/MLIAP/mliap_data.cpp b/src/MLIAP/mliap_data.cpp new file mode 100644 index 0000000000..5c927ca453 --- /dev/null +++ b/src/MLIAP/mliap_data.cpp @@ -0,0 +1,297 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, 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 +#include +#include "mliap_data.h" +#include "mliap_model_linear.h" +#include "mliap_model_quadratic.h" +#include "mliap_descriptor_snap.h" +#include "compute_mliap.h" +#include "atom.h" +#include "update.h" +#include "modify.h" +#include "neighbor.h" +#include "neigh_list.h" +#include "neigh_request.h" +#include "force.h" +#include "pair.h" +#include "comm.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; + +MLIAPData::MLIAPData(LAMMPS *lmp, + int gradgradflag_in, int *map_in, class MLIAPModel* model_in, + class MLIAPDescriptor* descriptor_in, class PairMLIAP* pairmliap_in) : + Pointers(lmp), + list(NULL), + gradforce(NULL), + betas(NULL), descriptors(NULL), gamma_row_index(NULL), gamma_col_index(NULL), + gamma(NULL), egradient(NULL), model(NULL), descriptor(NULL), + iatoms(NULL), ielems(NULL), numneighs(NULL), + jatoms(NULL), jelems(NULL), rij(NULL), graddesc(NULL) +{ + gradgradflag = gradgradflag_in; + map = map_in; + model = model_in; + descriptor = descriptor_in; + pairmliap = pairmliap_in; + + ndescriptors = descriptor->ndescriptors; + nelements = descriptor->nelements; + nparams = model->get_nparams(); + + gamma_nnz = model->get_gamma_nnz(this); + ndims_force = 3; + ndims_virial = 6; + yoffset = nparams*nelements; + zoffset = 2*yoffset; + natoms_array = atom->natoms; + size_array_rows = 1+ndims_force*natoms_array+ndims_virial; + size_array_cols = nparams*nelements+1; + size_gradforce = ndims_force*nparams*nelements; + + natoms_max = 0; + natomneigh_max = 0; + nneigh_max = 0; + nmax = 0; + natomgamma_max = 0; + +} + +/* ---------------------------------------------------------------------- */ + +MLIAPData::~MLIAPData() +{ + memory->destroy(betas); + memory->destroy(descriptors); + memory->destroy(gamma_row_index); + memory->destroy(gamma_col_index); + memory->destroy(gamma); + memory->destroy(egradient); + memory->destroy(gradforce); + + memory->destroy(iatoms); + memory->destroy(ielems); + memory->destroy(numneighs); + memory->destroy(jatoms); + memory->destroy(jelems); + memory->destroy(rij); + memory->destroy(graddesc); +} + +/* ---------------------------------------------------------------------- */ + +void MLIAPData::init() +{ + memory->create(egradient,nelements*nparams,"MLIAPData:egradient"); +} + +/* ---------------------------------------------------------------------- + generate neighbor arrays +------------------------------------------------------------------------- */ + +void MLIAPData::generate_neighdata(NeighList* list_in, int eflag_in, int vflag_in) +{ + list = list_in; + double **x = atom->x; + int *type = atom->type; + + int *numneigh = list->numneigh; + int **firstneigh = list->firstneigh; + + // grow nmax gradforce array if necessary + + if (atom->nmax > nmax) { + nmax = atom->nmax; + memory->grow(gradforce,nmax,size_gradforce, + "MLIAPData:gradforce"); + } + + // clear gradforce array + + int ntotal = atom->nlocal + atom->nghost; + for (int i = 0; i < ntotal; i++) + for (int j = 0; j < size_gradforce; j++) { + gradforce[i][j] = 0.0; + } + + // grow gamma arrays if necessary + + if (gradgradflag == 1) { + const int natomgamma = list->inum; + if (natomgamma_max < natomgamma) { + memory->grow(gamma_row_index,natomgamma,gamma_nnz,"MLIAPData:gamma_row_index"); + memory->grow(gamma_col_index,natomgamma,gamma_nnz,"MLIAPData:gamma_col_index"); + memory->grow(gamma,natomgamma,gamma_nnz,"MLIAPData:gamma"); + natomgamma_max = natomgamma; + } + } + + // grow arrays if necessary + + natoms = list->inum; + if (natoms_max < natoms) { + memory->grow(betas,natoms,ndescriptors,"MLIAPData:betas"); + memory->grow(descriptors,natoms,ndescriptors,"MLIAPData:descriptors"); + natoms_max = natoms; + } + + grow_neigharrays(); + + int ij = 0; + for (int ii = 0; ii < list->inum; ii++) { + const int i = list->ilist[ii]; + + const double xtmp = x[i][0]; + const double ytmp = x[i][1]; + const double ztmp = x[i][2]; + const int itype = type[i]; + const int ielem = map[itype]; + + int *jlist = firstneigh[i]; + const int jnum = numneigh[i]; + + int ninside = 0; + for (int jj = 0; jj < jnum; jj++) { + int j = jlist[jj]; + j &= NEIGHMASK; + const double delx = x[j][0] - xtmp; + const double dely = x[j][1] - ytmp; + const double delz = x[j][2] - ztmp; + const double rsq = delx*delx + dely*dely + delz*delz; + int jtype = type[j]; + const int jelem = map[jtype]; + + if (rsq < descriptor->cutsq[ielem][jelem]) { + jatoms[ij] = j; + jelems[ij] = jelem; + rij[ij][0] = delx; + rij[ij][1] = dely; + rij[ij][2] = delz; + ij++; + ninside++; + } + } + iatoms[ii] = i; + ielems[ii] = ielem; + numneighs[ii] = ninside; + } + + eflag = eflag_in; + vflag = vflag_in; +} + +/* ---------------------------------------------------------------------- + grow neighbor arrays to handle all neighbors +------------------------------------------------------------------------- */ + +void MLIAPData::grow_neigharrays() +{ + + // grow neighbor atom arrays if necessary + + const int natomneigh = list->inum; + if (natomneigh_max < natomneigh) { + memory->grow(iatoms,natomneigh,"MLIAPData:iatoms"); + memory->grow(ielems,natomneigh,"MLIAPData:ielems"); + memory->grow(numneighs,natomneigh,"MLIAPData:numneighs"); + natomneigh_max = natomneigh; + } + + // grow neighbor arrays if necessary + + int *numneigh = list->numneigh; + int **firstneigh = list->firstneigh; + + int iilast = list->inum-1; + int ilast = list->ilist[iilast]; + int upperbound = firstneigh[ilast] - firstneigh[0] + numneigh[ilast]; + if (nneigh_max >= upperbound) return; + + double **x = atom->x; + int *type = atom->type; + + int nneigh = 0; + for (int ii = 0; ii < list->inum; ii++) { + const int i = list->ilist[ii]; + + const double xtmp = x[i][0]; + const double ytmp = x[i][1]; + const double ztmp = x[i][2]; + const int itype = type[i]; + const int ielem = map[itype]; + + int *jlist = firstneigh[i]; + const int jnum = numneigh[i]; + + int ninside = 0; + for (int jj = 0; jj < jnum; jj++) { + int j = jlist[jj]; + j &= NEIGHMASK; + const double delx = x[j][0] - xtmp; + const double dely = x[j][1] - ytmp; + const double delz = x[j][2] - ztmp; + const double rsq = delx*delx + dely*dely + delz*delz; + int jtype = type[j]; + const int jelem = map[jtype]; + if (rsq < descriptor->cutsq[ielem][jelem]) ninside++; + } + nneigh += ninside; + } + + if (nneigh_max < nneigh) { + memory->grow(jatoms,nneigh,"MLIAPData:jatoms"); + memory->grow(jelems,nneigh,"MLIAPData:jelems"); + memory->grow(rij,nneigh,3,"MLIAPData:rij"); + if (gradgradflag == 0) + memory->grow(graddesc,nneigh,ndescriptors,3,"MLIAPData:graddesc"); + nneigh_max = nneigh; + } +} + +double MLIAPData::memory_usage() +{ + double bytes = 0.0; + + bytes += nelements*nparams*sizeof(double); // egradient + bytes += nmax*size_gradforce*sizeof(double); // gradforce + + if (gradgradflag == 1) { + bytes += natomgamma_max* + gamma_nnz*sizeof(int); //gamma_row_index + bytes += natomgamma_max* + gamma_nnz*sizeof(int); // gamma_col_index + bytes += natomgamma_max* + gamma_nnz*sizeof(double); // gamma + } + + bytes += natoms*ndescriptors*sizeof(int); // betas + bytes += natoms*ndescriptors*sizeof(int); // descriptors + + bytes += natomneigh_max*sizeof(int); // iatoms + bytes += natomneigh_max*sizeof(int); // ielems + bytes += natomneigh_max*sizeof(int); // numneighs + + bytes += nneigh_max*sizeof(int); // jatoms + bytes += nneigh_max*sizeof(int); // jelems + bytes += nneigh_max*3*sizeof(double); // rij" + + if (gradgradflag == 0) + bytes += nneigh_max*ndescriptors*3*sizeof(double);// graddesc + + return bytes; +} + diff --git a/src/MLIAP/mliap_data.h b/src/MLIAP/mliap_data.h new file mode 100644 index 0000000000..31aa52842c --- /dev/null +++ b/src/MLIAP/mliap_data.h @@ -0,0 +1,83 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, 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_MLIAPDATA_H +#define LMP_MLIAPDATA_H + +#include "pointers.h" + +namespace LAMMPS_NS { + +class MLIAPData : protected Pointers { + + public: + MLIAPData(class LAMMPS*, int, int*, class MLIAPModel*, class MLIAPDescriptor*, class PairMLIAP* = NULL); + ~MLIAPData(); + + void init(); + void generate_neighdata(class NeighList *, int=0, int=0); + void grow_neigharrays(); + double memory_usage(); + + int size_array_rows, size_array_cols; + int natoms_array, size_gradforce; + int yoffset, zoffset; + int ndims_force, ndims_virial; + double **gradforce; + double** betas; // betas for all atoms in list + double** descriptors; // descriptors for all atoms in list + int ndescriptors; // number of descriptors + int nparams; // number of model parameters per element + int nelements; // number of elements + + // data structures for grad-grad list (gamma) + + int natomgamma_max; // allocated size of gamma + int gamma_nnz; // number of non-zero entries in gamma + double** gamma; // gamma element + int** gamma_row_index; // row (parameter) index + int** gamma_col_index; // column (descriptor) index + double* egradient; // energy gradient w.r.t. parameters + + // data structures for mliap neighbor list + // only neighbors strictly inside descriptor cutoff + + int natoms; // current number of atoms + int natoms_max; // allocated size of descriptor array + int natomneigh_max; // allocated size of atom neighbor arrays + int *numneighs; // neighbors count for each atom + int *iatoms; // index of each atom + int *ielems; // element of each atom + int nneigh_max; // number of ij neighbors allocated + int *jatoms; // index of each neighbor + int *jelems; // element of each neighbor + double **rij; // distance vector of each neighbor + double ***graddesc; // descriptor gradient w.r.t. each neighbor + int eflag; // indicates if energy is needed + int vflag; // indicates if virial is needed + class PairMLIAP *pairmliap; // access to pair tally functions + + private: + class MLIAPModel* model; + class MLIAPDescriptor* descriptor; + + int nmax; + class NeighList *list; // LAMMPS neighbor list + int *map; // map LAMMPS types to [0,nelements) + int gradgradflag; // 1 for graddesc, 0 for gamma + +}; + +} + +#endif diff --git a/src/MLIAP/mliap_descriptor.h b/src/MLIAP/mliap_descriptor.h index 7e400c544f..03c8652b57 100644 --- a/src/MLIAP/mliap_descriptor.h +++ b/src/MLIAP/mliap_descriptor.h @@ -22,13 +22,10 @@ class MLIAPDescriptor : protected Pointers { public: MLIAPDescriptor(LAMMPS*); ~MLIAPDescriptor(); - virtual void compute_descriptors(int, int*, int*, int*, int*, int*, double**)=0; - virtual void compute_forces(int, int*, int*, int*, - int*, int*, double**, - class PairMLIAP*, int)=0; - virtual void compute_gradients(int, int*, int*, int*, int*, int*, int, int**, int**, double**, - double**, int, int)=0; - virtual void compute_descriptor_gradients(int, int*, int*, int*, int*, int*, double***)=0; + virtual void compute_descriptors(class MLIAPData*)=0; + virtual void compute_forces(class MLIAPData*)=0; + virtual void compute_force_gradients(class MLIAPData*)=0; + virtual void compute_descriptor_gradients(class MLIAPData*)=0; virtual void init()=0; virtual double memory_usage()=0; diff --git a/src/MLIAP/mliap_descriptor_snap.cpp b/src/MLIAP/mliap_descriptor_snap.cpp index 5d57f4516d..f4b7f4fa11 100644 --- a/src/MLIAP/mliap_descriptor_snap.cpp +++ b/src/MLIAP/mliap_descriptor_snap.cpp @@ -13,6 +13,7 @@ #include "mliap_descriptor_snap.h" #include "pair_mliap.h" +#include "mliap_data.h" #include #include #include @@ -73,39 +74,30 @@ MLIAPDescriptorSNAP::~MLIAPDescriptorSNAP() compute descriptors for each atom ---------------------------------------------------------------------- */ -void MLIAPDescriptorSNAP::compute_descriptors(int numlistdesc, - int *iatomdesc, int *ielemdesc, int *numneighdesc, - int *jatomdesc, int *jelemdesc, double **descriptors) +void MLIAPDescriptorSNAP::compute_descriptors(class MLIAPData* data) { double **x = atom->x; int ij = 0; - for (int ii = 0; ii < numlistdesc; ii++) { - const int i = iatomdesc[ii]; - const int ielem = ielemdesc[ii]; - - const double xtmp = x[i][0]; - const double ytmp = x[i][1]; - const double ztmp = x[i][2]; + for (int ii = 0; ii < data->natoms; ii++) { + const int i = data->iatoms[ii]; + const int ielem = data->ielems[ii]; // insure rij, inside, wj, and rcutij are of size jnum - const int jnum = numneighdesc[ii]; + const int jnum = data->numneighs[ii]; snaptr->grow_rij(jnum); int ninside = 0; for (int jj = 0; jj < jnum; jj++) { - const int j = jatomdesc[ij]; - const int jelem = jelemdesc[ij]; + const int j = data->jatoms[ij]; + const int jelem = data->jelems[ij]; + const double *delr = data->rij[ij]; - const double delx = x[j][0] - xtmp; - const double dely = x[j][1] - ytmp; - const double delz = x[j][2] - ztmp; - - snaptr->rij[ninside][0] = delx; - snaptr->rij[ninside][1] = dely; - snaptr->rij[ninside][2] = delz; + snaptr->rij[ninside][0] = delr[0]; + snaptr->rij[ninside][1] = delr[1]; + snaptr->rij[ninside][2] = delr[2]; snaptr->inside[ninside] = j; snaptr->wj[ninside] = wjelem[jelem]; snaptr->rcutij[ninside] = sqrt(cutsq[ielem][jelem]); @@ -125,8 +117,8 @@ void MLIAPDescriptorSNAP::compute_descriptors(int numlistdesc, else snaptr->compute_bi(0); - for (int icoeff = 0; icoeff < ndescriptors; icoeff++) - descriptors[ii][icoeff] = snaptr->blist[icoeff]; + for (int icoeff = 0; icoeff < data->ndescriptors; icoeff++) + data->descriptors[ii][icoeff] = snaptr->blist[icoeff]; } } @@ -135,11 +127,7 @@ void MLIAPDescriptorSNAP::compute_descriptors(int numlistdesc, compute forces for each atom ---------------------------------------------------------------------- */ -void MLIAPDescriptorSNAP::compute_forces(int natomdesc, - int *iatomdesc, int *ielemdesc, int *numneighdesc, - int *jatomdesc, int *jelemdesc, - double **beta, - PairMLIAP *pairmliap, int vflag) +void MLIAPDescriptorSNAP::compute_forces(class MLIAPData* data) { double fij[3]; @@ -147,31 +135,24 @@ void MLIAPDescriptorSNAP::compute_forces(int natomdesc, double **f = atom->f; int ij = 0; - for (int ii = 0; ii < natomdesc; ii++) { - const int i = iatomdesc[ii]; - const int ielem = ielemdesc[ii]; - - const double xtmp = x[i][0]; - const double ytmp = x[i][1]; - const double ztmp = x[i][2]; + for (int ii = 0; ii < data->natoms; ii++) { + const int i = data->iatoms[ii]; + const int ielem = data->ielems[ii]; // insure rij, inside, wj, and rcutij are of size jnum - const int jnum = numneighdesc[ii]; + const int jnum = data->numneighs[ii]; snaptr->grow_rij(jnum); int ninside = 0; for (int jj = 0; jj < jnum; jj++) { - const int j = jatomdesc[ij]; - const int jelem = jelemdesc[ij]; + const int j = data->jatoms[ij]; + const int jelem = data->jelems[ij]; + const double *delr = data->rij[ij]; - const double delx = x[j][0] - xtmp; - const double dely = x[j][1] - ytmp; - const double delz = x[j][2] - ztmp; - - snaptr->rij[ninside][0] = delx; - snaptr->rij[ninside][1] = dely; - snaptr->rij[ninside][2] = delz; + snaptr->rij[ninside][0] = delr[0]; + snaptr->rij[ninside][1] = delr[1]; + snaptr->rij[ninside][2] = delr[2]; snaptr->inside[ninside] = j; snaptr->wj[ninside] = wjelem[jelem]; snaptr->rcutij[ninside] = sqrt(cutsq[ielem][jelem]); @@ -191,7 +172,7 @@ void MLIAPDescriptorSNAP::compute_forces(int natomdesc, // compute Fij = dEi/dRj = -dEi/dRi // add to Fi, subtract from Fj - snaptr->compute_yi(beta[ii]); + snaptr->compute_yi(data->betas[ii]); for (int jj = 0; jj < ninside; jj++) { int j = snaptr->inside[jj]; @@ -214,52 +195,41 @@ void MLIAPDescriptorSNAP::compute_forces(int natomdesc, // add in global and per-atom virial contributions // this is optional and has no effect on force calculation - if (vflag) - pairmliap->v_tally(i,j,fij,snaptr->rij[jj]); + if (data->vflag) + data->pairmliap->v_tally(i,j,fij,snaptr->rij[jj]); } } } /* ---------------------------------------------------------------------- - compute force gradient for each atom + calculate gradients of forces w.r.t. parameters ---------------------------------------------------------------------- */ -void MLIAPDescriptorSNAP::compute_gradients(int natomdesc, - int *iatomdesc, int *ielemdesc, int *numneighdesc, - int *jatomdesc, int *jelemdesc, - int gamma_nnz, int **gamma_row_index, - int **gamma_col_index, double **gamma, double **gradforce, - int yoffset, int zoffset) +void MLIAPDescriptorSNAP::compute_force_gradients(class MLIAPData* data) { double **x = atom->x; int ij = 0; - for (int ii = 0; ii < natomdesc; ii++) { - const int i = iatomdesc[ii]; - const int ielem = ielemdesc[ii]; - - const double xtmp = x[i][0]; - const double ytmp = x[i][1]; - const double ztmp = x[i][2]; + for (int ii = 0; ii < data->natoms; ii++) { + const int i = data->iatoms[ii]; + const int ielem = data->ielems[ii]; // insure rij, inside, wj, and rcutij are of size jnum - const int jnum = numneighdesc[ii]; + const int jnum = data->numneighs[ii]; snaptr->grow_rij(jnum); int ninside = 0; for (int jj = 0; jj < jnum; jj++) { - const int j = jatomdesc[ij]; - const int jelem = jelemdesc[ij]; + const int j = data->jatoms[ij]; + const int jelem = data->jelems[ij]; - const double delx = x[j][0] - xtmp; - const double dely = x[j][1] - ytmp; - const double delz = x[j][2] - ztmp; + const double *delr = data->rij[ij]; - snaptr->rij[ninside][0] = delx; - snaptr->rij[ninside][1] = dely; - snaptr->rij[ninside][2] = delz; + snaptr->rij[ninside][0] = delr[0]; + snaptr->rij[ninside][1] = delr[1]; + snaptr->rij[ninside][2] = delr[2]; snaptr->inside[ninside] = j; snaptr->wj[ninside] = wjelem[jelem]; snaptr->rcutij[ninside] = sqrt(cutsq[ielem][jelem]); @@ -293,15 +263,15 @@ void MLIAPDescriptorSNAP::compute_gradients(int natomdesc, // Accumulate gamma_lk*dB_k/dRi, -gamma_lk**dB_k/dRj - for (int inz = 0; inz < gamma_nnz; inz++) { - const int l = gamma_row_index[ii][inz]; - const int k = gamma_col_index[ii][inz]; - gradforce[i][l] += gamma[ii][inz]*snaptr->dblist[k][0]; - gradforce[i][l+yoffset] += gamma[ii][inz]*snaptr->dblist[k][1]; - gradforce[i][l+zoffset] += gamma[ii][inz]*snaptr->dblist[k][2]; - gradforce[j][l] -= gamma[ii][inz]*snaptr->dblist[k][0]; - gradforce[j][l+yoffset] -= gamma[ii][inz]*snaptr->dblist[k][1]; - gradforce[j][l+zoffset] -= gamma[ii][inz]*snaptr->dblist[k][2]; + for (int inz = 0; inz < data->gamma_nnz; inz++) { + const int l = data->gamma_row_index[ii][inz]; + const int k = data->gamma_col_index[ii][inz]; + data->gradforce[i][l] += data->gamma[ii][inz]*snaptr->dblist[k][0]; + data->gradforce[i][l+data->yoffset] += data->gamma[ii][inz]*snaptr->dblist[k][1]; + data->gradforce[i][l+data->zoffset] += data->gamma[ii][inz]*snaptr->dblist[k][2]; + data->gradforce[j][l] -= data->gamma[ii][inz]*snaptr->dblist[k][0]; + data->gradforce[j][l+data->yoffset] -= data->gamma[ii][inz]*snaptr->dblist[k][1]; + data->gradforce[j][l+data->zoffset] -= data->gamma[ii][inz]*snaptr->dblist[k][2]; } } @@ -313,41 +283,30 @@ void MLIAPDescriptorSNAP::compute_gradients(int natomdesc, compute descriptor gradients for each neighbor atom ---------------------------------------------------------------------- */ -void MLIAPDescriptorSNAP::compute_descriptor_gradients(int numlistdesc, - int *iatomdesc, int *ielemdesc, - int *numneighdesc, - int *jatomdesc, int *jelemdesc, - double ***graddesc) +void MLIAPDescriptorSNAP::compute_descriptor_gradients(class MLIAPData* data) { double **x = atom->x; int ij = 0; - for (int ii = 0; ii < numlistdesc; ii++) { - const int i = iatomdesc[ii]; - const int ielem = ielemdesc[ii]; - - const double xtmp = x[i][0]; - const double ytmp = x[i][1]; - const double ztmp = x[i][2]; + for (int ii = 0; ii < data->natoms; ii++) { + const int i = data->iatoms[ii]; + const int ielem = data->ielems[ii]; // insure rij, inside, wj, and rcutij are of size jnum - const int jnum = numneighdesc[ii]; + const int jnum = data->numneighs[ii]; snaptr->grow_rij(jnum); int ij0 = ij; int ninside = 0; for (int jj = 0; jj < jnum; jj++) { - const int j = jatomdesc[ij]; - const int jelem = jelemdesc[ij]; + const int j = data->jatoms[ij]; + const int jelem = data->jelems[ij]; + const double *delr = data->rij[ij]; - const double delx = x[j][0] - xtmp; - const double dely = x[j][1] - ytmp; - const double delz = x[j][2] - ztmp; - - snaptr->rij[ninside][0] = delx; - snaptr->rij[ninside][1] = dely; - snaptr->rij[ninside][2] = delz; + snaptr->rij[ninside][0] = delr[0]; + snaptr->rij[ninside][1] = delr[1]; + snaptr->rij[ninside][2] = delr[2]; snaptr->inside[ninside] = j; snaptr->wj[ninside] = wjelem[jelem]; snaptr->rcutij[ninside] = sqrt(cutsq[ielem][jelem]); @@ -382,10 +341,10 @@ void MLIAPDescriptorSNAP::compute_descriptor_gradients(int numlistdesc, // Accumulate dB_k^i/dRi, dB_k^i/dRj - for (int k = 0; k < ndescriptors; k++) { - graddesc[ij][k][0] = snaptr->dblist[k][0]; - graddesc[ij][k][1] = snaptr->dblist[k][1]; - graddesc[ij][k][2] = snaptr->dblist[k][2]; + for (int k = 0; k < data->ndescriptors; k++) { + data->graddesc[ij][k][0] = snaptr->dblist[k][0]; + data->graddesc[ij][k][1] = snaptr->dblist[k][1]; + data->graddesc[ij][k][2] = snaptr->dblist[k][2]; } ij++; } @@ -568,7 +527,10 @@ double MLIAPDescriptorSNAP::memory_usage() { double bytes = 0; - bytes += snaptr->memory_usage(); // SNA object + bytes += nelements*sizeof(double); // radelem + bytes += nelements*sizeof(double); // welem + bytes += nelements*nelements*sizeof(int); // cutsq + bytes += snaptr->memory_usage(); // SNA object return bytes; } diff --git a/src/MLIAP/mliap_descriptor_snap.h b/src/MLIAP/mliap_descriptor_snap.h index 2f2c6c3791..40630e5131 100644 --- a/src/MLIAP/mliap_descriptor_snap.h +++ b/src/MLIAP/mliap_descriptor_snap.h @@ -22,18 +22,15 @@ class MLIAPDescriptorSNAP : public MLIAPDescriptor { public: MLIAPDescriptorSNAP(LAMMPS*, char*); ~MLIAPDescriptorSNAP(); - virtual void compute_descriptors(int, int*, int*, int*, int*, int*, double**); - virtual void compute_forces(int, int*, int*, int*, - int*, int*, double**, - class PairMLIAP*, int); - virtual void compute_gradients(int, int*, int*, int*, int*, int*, int, int**, int**, double**, - double**, int, int); - virtual void compute_descriptor_gradients(int, int*, int*, int*, int*, int*, double***); + virtual void compute_descriptors(class MLIAPData*); + virtual void compute_forces(class MLIAPData*); + virtual void compute_force_gradients(class MLIAPData*); + virtual void compute_descriptor_gradients(class MLIAPData*); virtual void init(); virtual double memory_usage(); - double rcutfac; // declared public to workaround gcc 4.9 - // compiler bug, manifest in KOKKOS package + double rcutfac; + protected: class SNA* snaptr; void read_paramfile(char *); diff --git a/src/MLIAP/mliap_model.cpp b/src/MLIAP/mliap_model.cpp index 1f88465304..bc760b7a68 100644 --- a/src/MLIAP/mliap_model.cpp +++ b/src/MLIAP/mliap_model.cpp @@ -33,19 +33,13 @@ using namespace LAMMPS_NS; MLIAPModel::MLIAPModel(LAMMPS* lmp, char* coefffilename) : Pointers(lmp) { - nelements = 0; - coeffelem = NULL; - read_coeffs(coefffilename); - nonlinearflag = 0; -} - -/* ---------------------------------------------------------------------- */ - -MLIAPModel::MLIAPModel(LAMMPS* lmp, int nelements_in, int nparams_in) : Pointers(lmp) -{ - nelements = nelements_in; - nparams = nparams_in; coeffelem = NULL; + if (coefffilename) read_coeffs(coefffilename); + else { + nparams = 0; + nelements = 0; + ndescriptors = 0; + } nonlinearflag = 0; } @@ -64,6 +58,24 @@ void MLIAPModel::init() { } +/* ---------------------------------------------------------------------- + set number of elements + ---------------------------------------------------------------------- */ + +void MLIAPModel::set_nelements(int nelements_in) +{ + nelements = nelements_in; +} + +/* ---------------------------------------------------------------------- + set number of descriptors + ---------------------------------------------------------------------- */ + +void MLIAPModel::set_ndescriptors(int ndescriptors_in) +{ + ndescriptors = ndescriptors_in; +} + /* ---------------------------------------------------------------------- */ void MLIAPModel::read_coeffs(char *coefffilename) diff --git a/src/MLIAP/mliap_model.h b/src/MLIAP/mliap_model.h index 98223c52d2..b2f7312897 100644 --- a/src/MLIAP/mliap_model.h +++ b/src/MLIAP/mliap_model.h @@ -21,13 +21,14 @@ namespace LAMMPS_NS { class MLIAPModel : protected Pointers { public: MLIAPModel(LAMMPS*, char*); - MLIAPModel(LAMMPS*, int, int); ~MLIAPModel(); - virtual void gradient(int, int*, int*, double**, double**, class PairMLIAP*, int)=0; - virtual void param_gradient(int, int*, int*, double**, int**, int**, double**, double*)=0; - virtual int get_gamma_nnz()=0; - virtual void compute_force_gradients(double**, int, int*, int*, int*, int*, - int*, double***, int, int, double**, double*)=0; + void set_ndescriptors(int); + void set_nelements(int); + virtual int get_nparams()=0; + virtual int get_gamma_nnz(class MLIAPData*)=0; + virtual void compute_gradients(class MLIAPData*)=0; + virtual void compute_gradgrads(class MLIAPData*)=0; + virtual void compute_force_gradients(class MLIAPData*)=0; virtual void init(); virtual double memory_usage(); int nelements; // # of unique elements diff --git a/src/MLIAP/mliap_model_linear.cpp b/src/MLIAP/mliap_model_linear.cpp index 55e7e339a3..f5c9c2f504 100644 --- a/src/MLIAP/mliap_model_linear.cpp +++ b/src/MLIAP/mliap_model_linear.cpp @@ -13,6 +13,8 @@ #include "mliap_model_linear.h" #include "pair_mliap.h" +#include "mliap_data.h" +#include "error.h" using namespace LAMMPS_NS; @@ -21,41 +23,46 @@ using namespace LAMMPS_NS; MLIAPModelLinear::MLIAPModelLinear(LAMMPS* lmp, char* coefffilename) : MLIAPModel(lmp, coefffilename) { - ndescriptors = nparams - 1; -} - -/* ---------------------------------------------------------------------- */ - -MLIAPModelLinear::MLIAPModelLinear(LAMMPS* lmp, int nelements_in, int nparams_in) : - MLIAPModel(lmp, nelements_in, nparams_in) -{ - ndescriptors = nparams - 1; + if (nparams > 0) ndescriptors = nparams - 1; } /* ---------------------------------------------------------------------- */ MLIAPModelLinear::~MLIAPModelLinear(){} +/* ---------------------------------------------------------------------- + get number of parameters + ---------------------------------------------------------------------- */ + +int MLIAPModelLinear::get_nparams() +{ + if (nparams == 0) { + if (ndescriptors == 0) error->all(FLERR,"ndescriptors not defined"); + else nparams = ndescriptors + 1; + } + + return nparams; +} + /* ---------------------------------------------------------------------- Calculate model gradients w.r.t descriptors for each atom beta_i = dE(B_i)/dB_i ---------------------------------------------------------------------- */ -void MLIAPModelLinear::gradient(int natomdesc, int *iatomdesc, int *ielemdesc, - double **descriptors, double **beta, PairMLIAP* pairmliap, int eflag) -{ - for (int ii = 0; ii < natomdesc; ii++) { - const int i = iatomdesc[ii]; - const int ielem = ielemdesc[ii]; +void MLIAPModelLinear::compute_gradients(MLIAPData* data) +{ + for (int ii = 0; ii < data->natoms; ii++) { + const int i = data->iatoms[ii]; + const int ielem = data->ielems[ii]; double* coeffi = coeffelem[ielem]; - for (int icoeff = 0; icoeff < ndescriptors; icoeff++) - beta[ii][icoeff] = coeffi[icoeff+1]; + for (int icoeff = 0; icoeff < data->ndescriptors; icoeff++) + data->betas[ii][icoeff] = coeffi[icoeff+1]; // add in contributions to global and per-atom energy // this is optional and has no effect on force calculation - if (eflag) { + if (data->eflag) { // energy of atom I @@ -64,10 +71,10 @@ void MLIAPModelLinear::gradient(int natomdesc, int *iatomdesc, int *ielemdesc, // E_i = beta.B_i - for (int icoeff = 0; icoeff < ndescriptors; icoeff++) - etmp += coeffi[icoeff+1]*descriptors[ii][icoeff]; + for (int icoeff = 0; icoeff < data->ndescriptors; icoeff++) + etmp += coeffi[icoeff+1]*data->descriptors[ii][icoeff]; - pairmliap->e_tally(i,etmp); + data->pairmliap->e_tally(i,etmp); } } } @@ -87,79 +94,67 @@ void MLIAPModelLinear::gradient(int natomdesc, int *iatomdesc, int *ielemdesc, egradient is derivative of energy w.r.t. parameters ---------------------------------------------------------------------- */ -void MLIAPModelLinear::param_gradient(int natommliap, int *iatommliap, int *ielemmliap, - double **descriptors, - int **gamma_row_index, int **gamma_col_index, - double **gamma, double *egradient) +void MLIAPModelLinear::compute_gradgrads(class MLIAPData* data) { // zero out energy gradients - for (int l = 0; l < nelements*nparams; l++) - egradient[l] = 0.0; + for (int l = 0; l < data->nelements*data->nparams; l++) + data->egradient[l] = 0.0; - for (int ii = 0; ii < natommliap; ii++) { - const int i = iatommliap[ii]; - const int ielem = ielemmliap[ii]; - const int elemoffset = nparams*ielem; + for (int ii = 0; ii < data->natoms; ii++) { + const int i = data->iatoms[ii]; + const int ielem = data->ielems[ii]; + const int elemoffset = data->nparams*ielem; int l = elemoffset+1; - for (int icoeff = 0; icoeff < ndescriptors; icoeff++) { - gamma[ii][icoeff] = 1.0; - gamma_row_index[ii][icoeff] = l++; - gamma_col_index[ii][icoeff] = icoeff; + for (int icoeff = 0; icoeff < data->ndescriptors; icoeff++) { + data->gamma[ii][icoeff] = 1.0; + data->gamma_row_index[ii][icoeff] = l++; + data->gamma_col_index[ii][icoeff] = icoeff; } // gradient of energy of atom I w.r.t. parameters l = elemoffset; - egradient[l++] += 1.0; - for (int icoeff = 0; icoeff < ndescriptors; icoeff++) - egradient[l++] += descriptors[ii][icoeff]; + data->egradient[l++] += 1.0; + for (int icoeff = 0; icoeff < data->ndescriptors; icoeff++) + data->egradient[l++] += data->descriptors[ii][icoeff]; } } /* ---------------------------------------------------------------------- - count the number of non-zero entries in gamma matrix + calculate gradients of forces w.r.t. parameters + egradient is derivative of energy w.r.t. parameters ---------------------------------------------------------------------- */ -int MLIAPModelLinear::get_gamma_nnz() +void MLIAPModelLinear::compute_force_gradients(class MLIAPData* data) { - int inz = ndescriptors; - return inz; -} - -void MLIAPModelLinear::compute_force_gradients(double **descriptors, int numlistdesc, - int *iatomdesc, int *ielemdesc, int *numneighdesc, - int *jatomdesc, int *jelemdesc, double ***graddesc, - int yoffset, int zoffset, double **gradforce, - double *egradient) { - // zero out energy gradients - for (int l = 0; l < nelements*nparams; l++) - egradient[l] = 0.0; + for (int l = 0; l < data->nelements*data->nparams; l++) + data->egradient[l] = 0.0; int ij = 0; - for (int ii = 0; ii < numlistdesc; ii++) { - const int i = iatomdesc[ii]; - const int ielem = ielemdesc[ii]; - const int elemoffset = nparams*ielem; + for (int ii = 0; ii < data->natoms; ii++) { + const int i = data->iatoms[ii]; + const int ielem = data->ielems[ii]; + const int elemoffset = data->nparams*ielem; - for (int jj = 0; jj < numneighdesc[ii]; jj++) { - const int j = jatomdesc[ij]; - const int jelem = ielemdesc[ij]; + for (int jj = 0; jj < data->numneighs[ii]; jj++) { + const int j = data->jatoms[ij]; + const int jelem = data->ielems[ij]; int l = elemoffset+1; - for (int icoeff = 0; icoeff < ndescriptors; icoeff++) { - gradforce[i][l] += graddesc[ij][icoeff][0]; - gradforce[i][l+yoffset] += graddesc[ij][icoeff][1]; - gradforce[i][l+zoffset] += graddesc[ij][icoeff][2]; - gradforce[j][l] -= graddesc[ij][icoeff][0]; - gradforce[j][l+yoffset] -= graddesc[ij][icoeff][1]; - gradforce[j][l+zoffset] -= graddesc[ij][icoeff][2]; + for (int icoeff = 0; icoeff < data->ndescriptors; icoeff++) { + data->gradforce[i][l] += data->graddesc[ij][icoeff][0]; + data->gradforce[i][l+data->yoffset] += data->graddesc[ij][icoeff][1]; + data->gradforce[i][l+data->zoffset] += data->graddesc[ij][icoeff][2]; + data->gradforce[j][l] -= data->graddesc[ij][icoeff][0]; + data->gradforce[j][l+data->yoffset] -= data->graddesc[ij][icoeff][1]; + data->gradforce[j][l+data->zoffset] -= data->graddesc[ij][icoeff][2]; l++; } ij++; @@ -168,10 +163,21 @@ void MLIAPModelLinear::compute_force_gradients(double **descriptors, int numlist // gradient of energy of atom I w.r.t. parameters int l = elemoffset; - egradient[l++] += 1.0; - for (int icoeff = 0; icoeff < ndescriptors; icoeff++) - egradient[l++] += descriptors[ii][icoeff]; + data->egradient[l++] += 1.0; + for (int icoeff = 0; icoeff < data->ndescriptors; icoeff++) + data->egradient[l++] += data->descriptors[ii][icoeff]; } } + +/* ---------------------------------------------------------------------- + count the number of non-zero entries in gamma matrix + ---------------------------------------------------------------------- */ + +int MLIAPModelLinear::get_gamma_nnz(class MLIAPData* data) +{ + int inz = data->ndescriptors; + return inz; +} + diff --git a/src/MLIAP/mliap_model_linear.h b/src/MLIAP/mliap_model_linear.h index 11a9e041ff..67ba0985a2 100644 --- a/src/MLIAP/mliap_model_linear.h +++ b/src/MLIAP/mliap_model_linear.h @@ -20,14 +20,13 @@ namespace LAMMPS_NS { class MLIAPModelLinear : public MLIAPModel { public: - MLIAPModelLinear(LAMMPS*, char*); - MLIAPModelLinear(LAMMPS*, int, int); + MLIAPModelLinear(LAMMPS*, char* = NULL); ~MLIAPModelLinear(); - virtual void gradient(int, int*, int*, double**, double**, class PairMLIAP*, int); - virtual void param_gradient(int, int*, int*, double**, int**, int**, double**, double*); - virtual int get_gamma_nnz(); - virtual void compute_force_gradients(double**, int, int*, int*, int*, int*, - int*, double***, int, int, double**, double*); + virtual int get_nparams(); + virtual int get_gamma_nnz(class MLIAPData*); + virtual void compute_gradients(class MLIAPData*); + virtual void compute_gradgrads(class MLIAPData*); + virtual void compute_force_gradients(class MLIAPData*); protected: }; diff --git a/src/MLIAP/mliap_model_quadratic.cpp b/src/MLIAP/mliap_model_quadratic.cpp index f642b07e22..3086ad088f 100644 --- a/src/MLIAP/mliap_model_quadratic.cpp +++ b/src/MLIAP/mliap_model_quadratic.cpp @@ -13,6 +13,8 @@ #include "mliap_model_quadratic.h" #include "pair_mliap.h" +#include "mliap_data.h" +#include "error.h" #include using namespace LAMMPS_NS; @@ -22,48 +24,51 @@ using namespace LAMMPS_NS; MLIAPModelQuadratic::MLIAPModelQuadratic(LAMMPS* lmp, char* coefffilename) : MLIAPModel(lmp, coefffilename) { + if (nparams > 0) ndescriptors = sqrt(2*nparams)-1; nonlinearflag = 1; - ndescriptors = sqrt(2*nparams)-1; -} - -/* ---------------------------------------------------------------------- */ - -MLIAPModelQuadratic::MLIAPModelQuadratic(LAMMPS* lmp, int nelements_in, int nparams_in) : - MLIAPModel(lmp, nelements_in, nparams_in) -{ - nonlinearflag = 1; - ndescriptors = sqrt(2*nparams)-1; } /* ---------------------------------------------------------------------- */ MLIAPModelQuadratic::~MLIAPModelQuadratic(){} +/* ---------------------------------------------------------------------- + get number of parameters + ---------------------------------------------------------------------- */ + +int MLIAPModelQuadratic::get_nparams() +{ + if (nparams == 0) { + if (ndescriptors == 0) error->all(FLERR,"ndescriptors not defined"); + else nparams = ndescriptors + 1 + (ndescriptors*(ndescriptors+1))/2; + } + + return nparams; +} + /* ---------------------------------------------------------------------- Calculate model gradients w.r.t descriptors for each atom dE(B_i)/dB_i ---------------------------------------------------------------------- */ -void MLIAPModelQuadratic::gradient(int natomdesc, int *iatomdesc, int *ielemdesc, - double **descriptors, double **beta, - PairMLIAP* pairmliap, int eflag) +void MLIAPModelQuadratic::compute_gradients(MLIAPData* data) { - for (int ii = 0; ii < natomdesc; ii++) { - const int i = iatomdesc[ii]; - const int ielem = ielemdesc[ii]; + for (int ii = 0; ii < data->natoms; ii++) { + const int i = data->iatoms[ii]; + const int ielem = data->ielems[ii]; double* coeffi = coeffelem[ielem]; - for (int icoeff = 0; icoeff < ndescriptors; icoeff++) - beta[ii][icoeff] = coeffi[icoeff+1]; + for (int icoeff = 0; icoeff < data->ndescriptors; icoeff++) + data->betas[ii][icoeff] = coeffi[icoeff+1]; int k = ndescriptors+1; - for (int icoeff = 0; icoeff < ndescriptors; icoeff++) { - double bveci = descriptors[ii][icoeff]; - beta[ii][icoeff] += coeffi[k]*bveci; + for (int icoeff = 0; icoeff < data->ndescriptors; icoeff++) { + double bveci = data->descriptors[ii][icoeff]; + data->betas[ii][icoeff] += coeffi[k]*bveci; k++; - for (int jcoeff = icoeff+1; jcoeff < ndescriptors; jcoeff++) { - double bvecj = descriptors[ii][jcoeff]; - beta[ii][icoeff] += coeffi[k]*bvecj; - beta[ii][jcoeff] += coeffi[k]*bveci; + for (int jcoeff = icoeff+1; jcoeff < data->ndescriptors; jcoeff++) { + double bvecj = data->descriptors[ii][jcoeff]; + data->betas[ii][icoeff] += coeffi[k]*bvecj; + data->betas[ii][jcoeff] += coeffi[k]*bveci; k++; } } @@ -71,7 +76,7 @@ void MLIAPModelQuadratic::gradient(int natomdesc, int *iatomdesc, int *ielemdesc // add in contributions to global and per-atom energy // this is optional and has no effect on force calculation - if (eflag) { + if (data->eflag) { // energy of atom I @@ -80,21 +85,21 @@ void MLIAPModelQuadratic::gradient(int natomdesc, int *iatomdesc, int *ielemdesc // E_i = beta.B_i + 0.5*B_i^t.alpha.B_i - for (int icoeff = 0; icoeff < ndescriptors; icoeff++) - etmp += coeffi[icoeff+1]*descriptors[ii][icoeff]; + for (int icoeff = 0; icoeff < data->ndescriptors; icoeff++) + etmp += coeffi[icoeff+1]*data->descriptors[ii][icoeff]; // quadratic contributions int k = ndescriptors+1; - for (int icoeff = 0; icoeff < ndescriptors; icoeff++) { - double bveci = descriptors[ii][icoeff]; + for (int icoeff = 0; icoeff < data->ndescriptors; icoeff++) { + double bveci = data->descriptors[ii][icoeff]; etmp += 0.5*coeffi[k++]*bveci*bveci; - for (int jcoeff = icoeff+1; jcoeff < ndescriptors; jcoeff++) { - double bvecj = descriptors[ii][jcoeff]; + for (int jcoeff = icoeff+1; jcoeff < data->ndescriptors; jcoeff++) { + double bvecj = data->descriptors[ii][jcoeff]; etmp += coeffi[k++]*bveci*bvecj; } } - pairmliap->e_tally(i,etmp); + data->pairmliap->e_tally(i,etmp); } } } @@ -114,48 +119,45 @@ void MLIAPModelQuadratic::gradient(int natomdesc, int *iatomdesc, int *ielemdesc egradient is derivative of energy w.r.t. parameters ---------------------------------------------------------------------- */ -void MLIAPModelQuadratic::param_gradient(int natommliap, int *iatommliap, int *ielemmliap, - double **descriptors, - int **gamma_row_index, int **gamma_col_index, - double **gamma, double *egradient) +void MLIAPModelQuadratic::compute_gradgrads(class MLIAPData* data) { // zero out energy gradients - for (int l = 0; l < nelements*nparams; l++) - egradient[l] = 0.0; + for (int l = 0; l < data->nelements*data->nparams; l++) + data->egradient[l] = 0.0; - for (int ii = 0; ii < natommliap; ii++) { - const int i = iatommliap[ii]; - const int ielem = ielemmliap[ii]; - const int elemoffset = nparams*ielem; + for (int ii = 0; ii < data->natoms; ii++) { + const int i = data->iatoms[ii]; + const int ielem = data->ielems[ii]; + const int elemoffset = data->nparams*ielem; // linear contributions int l = elemoffset+1; - for (int icoeff = 0; icoeff < ndescriptors; icoeff++) { - gamma[ii][icoeff] = 1.0; - gamma_row_index[ii][icoeff] = l++; - gamma_col_index[ii][icoeff] = icoeff; + for (int icoeff = 0; icoeff < data->ndescriptors; icoeff++) { + data->gamma[ii][icoeff] = 1.0; + data->gamma_row_index[ii][icoeff] = l++; + data->gamma_col_index[ii][icoeff] = icoeff; } // quadratic contributions - int inz = ndescriptors; - for (int icoeff = 0; icoeff < ndescriptors; icoeff++) { - double bveci = descriptors[ii][icoeff]; - gamma[ii][inz] = bveci; - gamma_row_index[ii][inz] = l++; - gamma_col_index[ii][inz] = icoeff; + int inz = data->ndescriptors; + for (int icoeff = 0; icoeff < data->ndescriptors; icoeff++) { + double bveci = data->descriptors[ii][icoeff]; + data->gamma[ii][inz] = bveci; + data->gamma_row_index[ii][inz] = l++; + data->gamma_col_index[ii][inz] = icoeff; inz++; - for (int jcoeff = icoeff+1; jcoeff < ndescriptors; jcoeff++) { - double bvecj = descriptors[ii][jcoeff]; - gamma[ii][inz] = bvecj; // derivative w.r.t. B[icoeff] - gamma_row_index[ii][inz] = l; - gamma_col_index[ii][inz] = icoeff; + for (int jcoeff = icoeff+1; jcoeff < data->ndescriptors; jcoeff++) { + double bvecj = data->descriptors[ii][jcoeff]; + data->gamma[ii][inz] = bvecj; // derivative w.r.t. B[icoeff] + data->gamma_row_index[ii][inz] = l; + data->gamma_col_index[ii][inz] = icoeff; inz++; - gamma[ii][inz] = bveci; // derivative w.r.t. B[jcoeff] - gamma_row_index[ii][inz] = l; - gamma_col_index[ii][inz] = jcoeff; + data->gamma[ii][inz] = bveci; // derivative w.r.t. B[jcoeff] + data->gamma_row_index[ii][inz] = l; + data->gamma_col_index[ii][inz] = jcoeff; inz++; l++; } @@ -164,18 +166,18 @@ void MLIAPModelQuadratic::param_gradient(int natommliap, int *iatommliap, int *i // gradient of energy of atom I w.r.t. parameters l = elemoffset; - egradient[l++] += 1.0; - for (int icoeff = 0; icoeff < ndescriptors; icoeff++) - egradient[l++] += descriptors[ii][icoeff]; + data->egradient[l++] += 1.0; + for (int icoeff = 0; icoeff < data->ndescriptors; icoeff++) + data->egradient[l++] += data->descriptors[ii][icoeff]; // quadratic contributions - for (int icoeff = 0; icoeff < ndescriptors; icoeff++) { - double bveci = descriptors[ii][icoeff]; - egradient[l++] += 0.5*bveci*bveci; + for (int icoeff = 0; icoeff < data->ndescriptors; icoeff++) { + double bveci = data->descriptors[ii][icoeff]; + data->egradient[l++] += 0.5*bveci*bveci; for (int jcoeff = icoeff+1; jcoeff < ndescriptors; jcoeff++) { - double bvecj = descriptors[ii][jcoeff]; - egradient[l++] += bveci*bvecj; + double bvecj = data->descriptors[ii][jcoeff]; + data->egradient[l++] += bveci*bvecj; } } } @@ -186,12 +188,12 @@ void MLIAPModelQuadratic::param_gradient(int natommliap, int *iatommliap, int *i count the number of non-zero entries in gamma matrix ---------------------------------------------------------------------- */ -int MLIAPModelQuadratic::get_gamma_nnz() +int MLIAPModelQuadratic::get_gamma_nnz(class MLIAPData* data) { int inz = ndescriptors; - for (int icoeff = 0; icoeff < ndescriptors; icoeff++) { + for (int icoeff = 0; icoeff < data->ndescriptors; icoeff++) { inz++; - for (int jcoeff = icoeff+1; jcoeff < ndescriptors; jcoeff++) { + for (int jcoeff = icoeff+1; jcoeff < data->ndescriptors; jcoeff++) { inz++; inz++; } @@ -200,58 +202,53 @@ int MLIAPModelQuadratic::get_gamma_nnz() return inz; } -void MLIAPModelQuadratic::compute_force_gradients(double **descriptors, int numlistdesc, - int *iatomdesc, int *ielemdesc, - int *numneighdesc, int *jatomdesc, - int *jelemdesc, double ***graddesc, - int yoffset, int zoffset, double **gradforce, - double *egradient) { +void MLIAPModelQuadratic::compute_force_gradients(class MLIAPData* data) { // zero out energy gradients - for (int l = 0; l < nelements*nparams; l++) - egradient[l] = 0.0; + for (int l = 0; l < data->nelements*data->nparams; l++) + data->egradient[l] = 0.0; int ij = 0; - for (int ii = 0; ii < numlistdesc; ii++) { - const int i = iatomdesc[ii]; - const int ielem = ielemdesc[ii]; - const int elemoffset = nparams*ielem; + for (int ii = 0; ii < data->natoms; ii++) { + const int i = data->iatoms[ii]; + const int ielem = data->ielems[ii]; + const int elemoffset = data->nparams*ielem; - for (int jj = 0; jj < numneighdesc[ii]; jj++) { - const int j = jatomdesc[ij]; - const int jelem = ielemdesc[ij]; + for (int jj = 0; jj < data->numneighs[ii]; jj++) { + const int j = data->jatoms[ij]; + const int jelem = data->ielems[ij]; const int ij0 = ij; int l = elemoffset+1; - for (int icoeff = 0; icoeff < ndescriptors; icoeff++) { - gradforce[i][l] += graddesc[ij][icoeff][0]; - gradforce[i][l+yoffset] += graddesc[ij][icoeff][1]; - gradforce[i][l+zoffset] += graddesc[ij][icoeff][2]; - gradforce[j][l] -= graddesc[ij][icoeff][0]; - gradforce[j][l+yoffset] -= graddesc[ij][icoeff][1]; - gradforce[j][l+zoffset] -= graddesc[ij][icoeff][2]; + for (int icoeff = 0; icoeff < data->ndescriptors; icoeff++) { + data->gradforce[i][l] += data->graddesc[ij][icoeff][0]; + data->gradforce[i][l+data->yoffset] += data->graddesc[ij][icoeff][1]; + data->gradforce[i][l+data->zoffset] += data->graddesc[ij][icoeff][2]; + data->gradforce[j][l] -= data->graddesc[ij][icoeff][0]; + data->gradforce[j][l+data->yoffset] -= data->graddesc[ij][icoeff][1]; + data->gradforce[j][l+data->zoffset] -= data->graddesc[ij][icoeff][2]; l++; } // quadratic contributions - for (int icoeff = 0; icoeff < ndescriptors; icoeff++) { - double bveci = descriptors[ii][icoeff]; - gradforce[i][l] += graddesc[ij][icoeff][0]*bveci; - gradforce[i][l+yoffset] += graddesc[ij][icoeff][1]*bveci; - gradforce[i][l+zoffset] += graddesc[ij][icoeff][2]*bveci; - gradforce[j][l] -= graddesc[ij][icoeff][0]*bveci; - gradforce[j][l+yoffset] -= graddesc[ij][icoeff][1]*bveci; - gradforce[j][l+zoffset] -= graddesc[ij][icoeff][2]*bveci; + for (int icoeff = 0; icoeff < data->ndescriptors; icoeff++) { + double bveci = data->descriptors[ii][icoeff]; + data->gradforce[i][l] += data->graddesc[ij][icoeff][0]*bveci; + data->gradforce[i][l+data->yoffset] += data->graddesc[ij][icoeff][1]*bveci; + data->gradforce[i][l+data->zoffset] += data->graddesc[ij][icoeff][2]*bveci; + data->gradforce[j][l] -= data->graddesc[ij][icoeff][0]*bveci; + data->gradforce[j][l+data->yoffset] -= data->graddesc[ij][icoeff][1]*bveci; + data->gradforce[j][l+data->zoffset] -= data->graddesc[ij][icoeff][2]*bveci; l++; - for (int jcoeff = icoeff+1; jcoeff < ndescriptors; jcoeff++) { - double bvecj = descriptors[ii][jcoeff]; - gradforce[i][l] += graddesc[ij][icoeff][0]*bvecj + graddesc[ij][jcoeff][0]*bveci; - gradforce[i][l+yoffset] += graddesc[ij][icoeff][1]*bvecj + graddesc[ij][jcoeff][1]*bveci; - gradforce[i][l+zoffset] += graddesc[ij][icoeff][2]*bvecj + graddesc[ij][jcoeff][2]*bveci; - gradforce[j][l] -= graddesc[ij][icoeff][0]*bvecj + graddesc[ij][jcoeff][0]*bveci; - gradforce[j][l+yoffset] -= graddesc[ij][icoeff][1]*bvecj + graddesc[ij][jcoeff][1]*bveci; - gradforce[j][l+zoffset] -= graddesc[ij][icoeff][2]*bvecj + graddesc[ij][jcoeff][2]*bveci; + for (int jcoeff = icoeff+1; jcoeff < data->ndescriptors; jcoeff++) { + double bvecj = data->descriptors[ii][jcoeff]; + data->gradforce[i][l] += data->graddesc[ij][icoeff][0]*bvecj + data->graddesc[ij][jcoeff][0]*bveci; + data->gradforce[i][l+data->yoffset] += data->graddesc[ij][icoeff][1]*bvecj + data->graddesc[ij][jcoeff][1]*bveci; + data->gradforce[i][l+data->zoffset] += data->graddesc[ij][icoeff][2]*bvecj + data->graddesc[ij][jcoeff][2]*bveci; + data->gradforce[j][l] -= data->graddesc[ij][icoeff][0]*bvecj + data->graddesc[ij][jcoeff][0]*bveci; + data->gradforce[j][l+data->yoffset] -= data->graddesc[ij][icoeff][1]*bvecj + data->graddesc[ij][jcoeff][1]*bveci; + data->gradforce[j][l+data->zoffset] -= data->graddesc[ij][icoeff][2]*bvecj + data->graddesc[ij][jcoeff][2]*bveci; l++; } } @@ -261,18 +258,18 @@ void MLIAPModelQuadratic::compute_force_gradients(double **descriptors, int numl // gradient of energy of atom I w.r.t. parameters int l = elemoffset; - egradient[l++] += 1.0; - for (int icoeff = 0; icoeff < ndescriptors; icoeff++) - egradient[l++] += descriptors[ii][icoeff]; + data->egradient[l++] += 1.0; + for (int icoeff = 0; icoeff < data->ndescriptors; icoeff++) + data->egradient[l++] += data->descriptors[ii][icoeff]; // quadratic contributions - for (int icoeff = 0; icoeff < ndescriptors; icoeff++) { - double bveci = descriptors[ii][icoeff]; - egradient[l++] += 0.5*bveci*bveci; - for (int jcoeff = icoeff+1; jcoeff < ndescriptors; jcoeff++) { - double bvecj = descriptors[ii][jcoeff]; - egradient[l++] += bveci*bvecj; + for (int icoeff = 0; icoeff < data->ndescriptors; icoeff++) { + double bveci = data->descriptors[ii][icoeff]; + data->egradient[l++] += 0.5*bveci*bveci; + for (int jcoeff = icoeff+1; jcoeff < data->ndescriptors; jcoeff++) { + double bvecj = data->descriptors[ii][jcoeff]; + data->egradient[l++] += bveci*bvecj; } } diff --git a/src/MLIAP/mliap_model_quadratic.h b/src/MLIAP/mliap_model_quadratic.h index 5024b2e426..188d78661f 100644 --- a/src/MLIAP/mliap_model_quadratic.h +++ b/src/MLIAP/mliap_model_quadratic.h @@ -20,14 +20,13 @@ namespace LAMMPS_NS { class MLIAPModelQuadratic : public MLIAPModel { public: - MLIAPModelQuadratic(LAMMPS*, char*); - MLIAPModelQuadratic(LAMMPS*, int, int); + MLIAPModelQuadratic(LAMMPS*, char* = NULL); ~MLIAPModelQuadratic(); - virtual void gradient(int, int*, int*, double**, double**, class PairMLIAP*, int); - virtual void param_gradient(int, int*, int*, double**, int**, int**, double**, double*); - virtual int get_gamma_nnz(); - virtual void compute_force_gradients(double**, int, int*, int*, int*, int*, - int*, double***, int, int, double**, double*); + virtual int get_nparams(); + virtual int get_gamma_nnz(class MLIAPData*); + virtual void compute_gradients(class MLIAPData*); + virtual void compute_gradgrads(class MLIAPData*); + virtual void compute_force_gradients(class MLIAPData*); protected: }; diff --git a/src/MLIAP/pair_mliap.cpp b/src/MLIAP/pair_mliap.cpp index 04719c749b..a7598a3a03 100644 --- a/src/MLIAP/pair_mliap.cpp +++ b/src/MLIAP/pair_mliap.cpp @@ -15,7 +15,7 @@ #include #include #include -#include "mliap.h" +#include "mliap_data.h" #include "mliap_model_linear.h" #include "mliap_model_quadratic.h" #include "mliap_descriptor_snap.h" @@ -49,7 +49,7 @@ PairMLIAP::~PairMLIAP() delete model; delete descriptor; - delete mliap; + delete data; if (allocated) { memory->destroy(setflag); @@ -67,22 +67,20 @@ void PairMLIAP::compute(int eflag, int vflag) { ev_init(eflag,vflag); - mliap->generate_neigharrays(list); + data->generate_neighdata(list, eflag, vflag); // compute descriptors, if needed if (model->nonlinearflag || eflag) - descriptor->compute_descriptors(mliap->natomdesc, mliap->iatommliap, mliap->ielemmliap, mliap->numneighmliap, - mliap->jatommliap, mliap->jelemmliap, mliap->descriptors); + descriptor->compute_descriptors(data); // compute E_i and beta_i = dE_i/dB_i for all i in list - model->gradient(mliap->natomdesc, mliap->iatommliap, mliap->ielemmliap, mliap->descriptors, mliap->beta, this, eflag); + model->compute_gradients(data); // calculate force contributions beta_i*dB_i/dR_j - descriptor->compute_forces(mliap->natomdesc, mliap->iatommliap, mliap->ielemmliap, mliap->numneighmliap, - mliap->jatommliap, mliap->jelemmliap, mliap->beta, this, vflag); + descriptor->compute_forces(data); // calculate stress @@ -149,10 +147,6 @@ void PairMLIAP::settings(int narg, char ** arg) if (modelflag == 0 || descriptorflag == 0) error->all(FLERR,"Illegal pair_style command"); - ndescriptors = descriptor->ndescriptors; - nparams = model->nparams; - nelements = model->nelements; - } /* ---------------------------------------------------------------------- @@ -209,20 +203,20 @@ void PairMLIAP::coeff(int narg, char **arg) if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); + // set up model, descriptor, and mliap data structures + model->init(); descriptor->init(); + int gradgradflag = -1; + data = new MLIAPData(lmp, gradgradflag, map, model, descriptor, this); + data->init(); // consistency checks - if (ndescriptors != model->ndescriptors) + if (data->ndescriptors != model->ndescriptors) error->all(FLERR,"Incompatible model and descriptor definitions"); - if (descriptor->nelements != model->nelements) + if (data->nelements != model->nelements) error->all(FLERR,"Incompatible model and descriptor definitions"); - - int gradgradflag = 0; - mliap = new MLIAP(lmp, ndescriptors, nparams, nelements, gradgradflag, map, model, descriptor); - mliap->init(); - } /* ---------------------------------------------------------------------- @@ -315,12 +309,12 @@ double PairMLIAP::memory_usage() double bytes = Pair::memory_usage(); int n = atom->ntypes+1; - bytes += n*n*sizeof(int); // setflag - bytes += mliap->natomdesc_max*ndescriptors*sizeof(double); // descriptors - bytes += mliap->natomdesc_max*ndescriptors*sizeof(double); // beta - + bytes += n*n*sizeof(int); // setflag + bytes += n*n*sizeof(int); // cutsq + bytes += n*sizeof(int); // map bytes += descriptor->memory_usage(); // Descriptor object bytes += model->memory_usage(); // Model object + bytes += data->memory_usage(); // Data object return bytes; } diff --git a/src/MLIAP/pair_mliap.h b/src/MLIAP/pair_mliap.h index e1aed97026..61acbbb278 100644 --- a/src/MLIAP/pair_mliap.h +++ b/src/MLIAP/pair_mliap.h @@ -41,13 +41,9 @@ public: protected: virtual void allocate(); - int ndescriptors; // number of descriptors - int nparams; // number of model parameters per element - int nelements; - class MLIAPModel* model; class MLIAPDescriptor* descriptor; - class MLIAP *mliap; + class MLIAPData *data; }; }