diff --git a/examples/PACKAGES/pace/README.md b/examples/PACKAGES/pace/README.md new file mode 100644 index 0000000000..ab2d3a54f5 --- /dev/null +++ b/examples/PACKAGES/pace/README.md @@ -0,0 +1,9 @@ +# This folder contains examples for pace in LAMMPS + + +## Compute pace usage +compute/latte_cell_0.data # lammps data file with C-H-O structure +compute/latte_cell_0.xyz # xyz file with C-H-O structure +compute/coupling_coefficients.yace # .yace file containing coupling coefficients (or ACE potential parameters) +compute/in.lammps # input file for calling `compute pace` + diff --git a/examples/PACKAGES/pace/compute/README.md b/examples/PACKAGES/pace/compute/README.md new file mode 100644 index 0000000000..e69de29bb2 diff --git a/examples/pace/multi-elem/coupling_coefficients.yace b/examples/PACKAGES/pace/compute/coupling_coefficients.yace similarity index 100% rename from examples/pace/multi-elem/coupling_coefficients.yace rename to examples/PACKAGES/pace/compute/coupling_coefficients.yace diff --git a/examples/PACKAGES/pace/compute/in.lammps b/examples/PACKAGES/pace/compute/in.lammps new file mode 100644 index 0000000000..a0ef25f606 --- /dev/null +++ b/examples/PACKAGES/pace/compute/in.lammps @@ -0,0 +1,22 @@ +#info all out log +units metal +atom_style atomic +boundary p p p +atom_modify map hash +boundary p p p +read_data latte_cell_0.data +mass 1 1.00 +mass 2 14.00 +mass 3 15.999 + + # potential settings + +pair_style zero 5.7 +pair_coeff * * + +compute pace all pace coupling_coefficients.yace 1 0 + +thermo 1 +thermo_style custom step temp c_pace[1][183] + +run 0 diff --git a/examples/pace/multi-elem/latte_cell_0.data b/examples/PACKAGES/pace/compute/latte_cell_0.data similarity index 100% rename from examples/pace/multi-elem/latte_cell_0.data rename to examples/PACKAGES/pace/compute/latte_cell_0.data diff --git a/examples/pace/multi-elem/latte_cell_0.xyz b/examples/PACKAGES/pace/compute/latte_cell_0.xyz similarity index 100% rename from examples/pace/multi-elem/latte_cell_0.xyz rename to examples/PACKAGES/pace/compute/latte_cell_0.xyz diff --git a/examples/pace/multi-elem/README.md b/examples/pace/multi-elem/README.md deleted file mode 100644 index e2ab583d10..0000000000 --- a/examples/pace/multi-elem/README.md +++ /dev/null @@ -1,5 +0,0 @@ -Check that seg fault doesn't occur by doing: - - ./loop.sh - -which will loop through many runs of `python test_en.py`. diff --git a/examples/pace/multi-elem/latte_cell_0_chi_i.npy b/examples/pace/multi-elem/latte_cell_0_chi_i.npy deleted file mode 100644 index 27f0878145..0000000000 Binary files a/examples/pace/multi-elem/latte_cell_0_chi_i.npy and /dev/null differ diff --git a/examples/pace/multi-elem/loop.sh b/examples/pace/multi-elem/loop.sh deleted file mode 100755 index 613636a30e..0000000000 --- a/examples/pace/multi-elem/loop.sh +++ /dev/null @@ -1,8 +0,0 @@ -for (( ; ; )) -do - python test_en.py - # terminate loop if seg fault - if [[ $? -eq 139 ]]; then - break - fi -done diff --git a/examples/pace/multi-elem/test_en.py b/examples/pace/multi-elem/test_en.py deleted file mode 100644 index 17c1924080..0000000000 --- a/examples/pace/multi-elem/test_en.py +++ /dev/null @@ -1,86 +0,0 @@ -from __future__ import print_function -import sys, os -import ctypes -import numpy as np -from ase.io import read,write -from lammps import lammps, LMP_TYPE_ARRAY, LMP_STYLE_GLOBAL - -# get MPI settings from LAMMPS -def run_struct(f): - file_prefix = f.split('.')[0] - atoms = read(f) - lmp = lammps() - - me = lmp.extract_setting("world_rank") - nprocs = lmp.extract_setting("world_size") - - write('%s.data' % file_prefix,atoms,format='lammps-data') - - cmds = ["-screen", "none", "-log", "none"] - lmp = lammps(cmdargs = cmds) - - print("Made LAMMPS instance") - - def run_lammps(dgradflag): - - # simulation settings - fname = file_prefix - lmp.command("clear") - lmp.command("info all out log") - lmp.command('units metal') - lmp.command('atom_style atomic') - lmp.command("boundary p p p") - lmp.command("atom_modify map hash") - lmp.command('neighbor 2.3 bin') - # boundary - lmp.command('boundary p p p') - # read atoms - lmp.command('read_data %s.data' % fname ) - lmp.command('mass 1 1.00') - lmp.command('mass 2 14.00') - lmp.command('mass 3 15.999') - - # potential settings - - lmp.command(f"pair_style zero 5.7") - lmp.command(f"pair_coeff * *") - - - if dgradflag: - lmp.command(f"compute pace all pace coupling_coefficients.yace 1 1") - else: - lmp.command(f"compute pace all pace coupling_coefficients.yace 1 0 ") - - # run - - lmp.command(f"thermo 100") - lmp.command(f"run {nsteps}") - - # declare simulation/structure variables - - nsteps = 0 - ntypes = 3 - - # declare compute pace variables - - - bikflag = 1 - - # NUMBER of descriptors - nd = 91 - - dgradflag = 0 - - run_lammps(dgradflag) - - lmp_pace = lmp.numpy.extract_compute("pace", LMP_STYLE_GLOBAL, LMP_TYPE_ARRAY) - print ('global shape',np.shape(lmp_pace)) - np.save('%s_chi_i.npy' % file_prefix,lmp_pace) - lmp.close() - del lmp - return None - -import glob -for f in sorted(glob.glob('*.xyz')): - print ('running %s' % f) - run_struct(f) diff --git a/src/ML-PACE/compute_pace.cpp b/src/ML-PACE/compute_pace.cpp index f9e5ff8a30..4d9e761d26 100644 --- a/src/ML-PACE/compute_pace.cpp +++ b/src/ML-PACE/compute_pace.cpp @@ -77,6 +77,7 @@ ComputePACE::ComputePACE(LAMMPS *lmp, int narg, char **arg) : //read in file with CG coefficients or c_tilde coefficients auto potential_file_name = utils::get_potential_file_path(arg[3]); + delete acecimpl -> basis_set; acecimpl -> basis_set = new ACECTildeBasisSet(potential_file_name); double cut = acecimpl -> basis_set->cutoffmax; cutmax = acecimpl -> basis_set->cutoffmax; @@ -138,6 +139,7 @@ ComputePACE::ComputePACE(LAMMPS *lmp, int narg, char **arg) : ComputePACE::~ComputePACE() { + delete acecimpl; memory->destroy(pace); memory->destroy(paceall); memory->destroy(cutsq); @@ -244,10 +246,6 @@ void ComputePACE::compute_array() NS_TYPE *ns; LS_TYPE *ls; - int n_r1, n_rp = 0; - n_r1 = acecimpl -> basis_set->total_basis_size_rank1[0]; - n_rp = acecimpl -> basis_set->total_basis_size[0]; - const int inum = list->inum; const int* const ilist = list->ilist; const int* const numneigh = list->numneigh; @@ -285,15 +283,18 @@ void ComputePACE::compute_array() const int typeoffset_local = ndims_peratom*nvalues*(itype-1); const int typeoffset_global = nvalues*(itype-1); + delete acecimpl -> ace; acecimpl -> ace = new ACECTildeEvaluator(*acecimpl -> basis_set); + acecimpl -> ace->compute_projections = 1; + acecimpl -> ace->compute_b_grad = 1; int n_r1, n_rp = 0; - n_r1 = basis_set->total_basis_size_rank1[0]; - n_rp = basis_set->total_basis_size[0]; + n_r1 = acecimpl -> basis_set->total_basis_size_rank1[0]; + n_rp = acecimpl -> basis_set->total_basis_size[0]; int ncoeff = n_r1 + n_rp; acecimpl -> ace->element_type_mapping.init(ntypes+1); for (int ik = 1; ik <= ntypes; ik++) { - for(int mu = 0; mu < basis_set->nelements; mu++){ + for(int mu = 0; mu < acecimpl -> basis_set ->nelements; mu++){ if (mu != -1) { if (mu == ik - 1) { map[ik] = mu; @@ -336,7 +337,7 @@ void ComputePACE::compute_array() // resize the neighbor cache after setting the basis acecimpl -> ace->resize_neighbours_cache(max_jnum); acecimpl -> ace->compute_atom(i, atom->x, atom->type, list->numneigh[i], list->firstneigh[i]); - Array1D Bs =ace->B_all; + Array1D Bs = acecimpl -> ace -> projections; for (int jj = 0; jj < jnum; jj++) { const int j = jlist[jj]; @@ -345,14 +346,13 @@ void ComputePACE::compute_array() double *pacedi = pace_peratom[i]+typeoffset_local; double *pacedj = pace_peratom[j]+typeoffset_local; - Array3D fs = acecimpl -> ace->neighbours_dB; //force array in (func_ind,neighbour_ind,xyz_ind) format // dimension: (n_descriptors,max_jnum,3) //example to access entries for neighbour jj after running compute_atom for atom i: for (int func_ind =0; func_ind < n_r1 + n_rp; func_ind++){ - DOUBLE_TYPE fx_dB = fs(func_ind,jj,0); - DOUBLE_TYPE fy_dB = fs(func_ind,jj,1); - DOUBLE_TYPE fz_dB = fs(func_ind,jj,2); + DOUBLE_TYPE fx_dB = acecimpl -> ace -> neighbours_dB(func_ind,jj,0); + DOUBLE_TYPE fy_dB = acecimpl -> ace -> neighbours_dB(func_ind,jj,1); + DOUBLE_TYPE fz_dB = acecimpl -> ace -> neighbours_dB(func_ind,jj,2); pacedi[func_ind] += fx_dB; pacedi[func_ind+yoffset] += fy_dB; pacedi[func_ind+zoffset] += fz_dB; @@ -362,15 +362,14 @@ void ComputePACE::compute_array() } } else { //printf("inside dBi/dRj logical : ncoeff = %d \n", ncoeff); - Array3D fs = acecimpl -> ace->neighbours_dB; for (int iicoeff = 0; iicoeff < ncoeff; iicoeff++) { // add to pace array for this proc //printf("inside dBi/dRj loop\n"); // dBi/dRj - DOUBLE_TYPE fx_dB = fs(iicoeff,jj,0); - DOUBLE_TYPE fy_dB = fs(iicoeff,jj,1); - DOUBLE_TYPE fz_dB = fs(iicoeff,jj,2); + DOUBLE_TYPE fx_dB = acecimpl -> ace -> neighbours_dB(iicoeff,jj,0); + DOUBLE_TYPE fy_dB = acecimpl -> ace -> neighbours_dB(iicoeff,jj,1); + DOUBLE_TYPE fz_dB = acecimpl -> ace -> neighbours_dB(iicoeff,jj,2); pace[bik_rows + ((atom->tag[j]-1)*3*natoms) + 3*(atom->tag[i]-1) + 0][iicoeff+3] -= fx_dB; pace[bik_rows + ((atom->tag[j]-1)*3*natoms) + 3*(atom->tag[i]-1) + 1][iicoeff+3] -= fy_dB; pace[bik_rows + ((atom->tag[j]-1)*3*natoms) + 3*(atom->tag[i]-1) + 2][iicoeff+3] -= fz_dB; @@ -395,10 +394,8 @@ void ComputePACE::compute_array() pace[irow][k++] += Bs(icoeff); } } - delete acecimpl -> ace; } //group bit } // for ii loop - delete acecimpl -> basis_set; // accumulate force contributions to global array if (!dgradflag){ for (int itype = 0; itype < atom->ntypes; itype++) {