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

Correctly reproduced examples/in.snap.compute, not yet for quadratic case

Browse Source
This commit is contained in:
Aidan Thompson
2020-07-02 14:32:47 -06:00
parent 2a4e51fa38
commit bc36511767
3 changed files with 21 additions and 27 deletions
Download Patch File Download Diff File Expand all files Collapse all files

12
src/MLIAP/compute_mliap.cpp
View File

@ -93,8 +93,8 @@ ComputeMLIAP::ComputeMLIAP(LAMMPS *lmp, int narg, char **arg) :
nperdim = nparams;
ndims_force = 3;
ndims_virial = 6;
yoffset = nperdim;
zoffset = 2*nperdim;
yoffset = nperdim*atom->ntypes;
zoffset = 2*yoffset;
natoms = atom->natoms;
size_array_rows = 1+ndims_force*natoms+ndims_virial;
size_array_cols = nperdim*atom->ntypes+1;
@ -169,14 +169,12 @@ void ComputeMLIAP::init()
// allocate memory for global array
printf("size_array_rows = %d size_array_cols = %d\n",size_array_rows,size_array_cols);
memory->create(mliap,size_array_rows,size_array_cols,
"mliap:mliap");
memory->create(mliapall,size_array_rows,size_array_cols,
"mliap:mliapall");
array = mliapall;
printf("nelements = %d nparams = %d\n",nelements,nparams);
memory->create(egradient,nelements*nparams,"ComputeMLIAP:egradient");
// find compute for reference energy
@ -227,7 +225,6 @@ void ComputeMLIAP::compute_array()
if (atom->nmax > nmax) {
memory->destroy(mliap_peratom);
nmax = atom->nmax;
printf("nmax = %d size_peratom = %d\n",nmax,size_peratom);
memory->create(mliap_peratom,nmax,size_peratom,
"mliap:mliap_peratom");
}
@ -256,7 +253,6 @@ void ComputeMLIAP::compute_array()
memory->grow(gamma,list->inum,gamma_nnz,"ComputeMLIAP:gamma");
gamma_max = list->inum;
}
printf("gamma_max %d %d %d %d %d %d %p\n",gamma_max, ndescriptors, gamma_nnz, nelements, nparams, list->inum, list);
// compute descriptors, if needed
@ -283,7 +279,7 @@ void ComputeMLIAP::compute_array()
// accumulate descriptor gradient contributions to global array
for (int itype = 0; itype < atom->ntypes; itype++) {
const int typeoffset_local = ndims_peratom*nperdim*itype;
const int typeoffset_local = nperdim*itype;
const int typeoffset_global = nperdim*itype;
for (int icoeff = 0; icoeff < nperdim; icoeff++) {
int irow = 1;
@ -360,7 +356,7 @@ void ComputeMLIAP::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_local = nperdim*itype;
const int typeoffset_global = nperdim*itype;
double *snadi = mliap_peratom[i]+typeoffset_local;
for (int icoeff = 0; icoeff < nperdim; icoeff++) {

35
src/MLIAP/mliap_descriptor_snap.cpp
View File

@ -228,9 +228,6 @@ void MLIAPDescriptorSNAP::backward(PairMLIAP* pairmliap, NeighList* list, double
// add to Fi, subtract from Fj
snaptr->compute_yi(beta[ii]);
//for (int q=0; q<snaptr->idxu_max*2; q++){
// fprintf(screen, "%i %f\n",q, snaptr->ylist_r[q]);
//}
for (int jj = 0; jj < ninside; jj++) {
int j = snaptr->inside[jj];
@ -294,11 +291,11 @@ void MLIAPDescriptorSNAP::param_backward(int *map, NeighList* list,
const double ytmp = x[i][1];
const double ztmp = x[i][2];
const int itype = type[i];
int ielem = 0;
if (chemflag)
int ielem = itype-1;
// element map not yet implemented
// if (chemflag)
// element map not yet implemented
// ielem = map[itype];
const int ielem = itype-1;
jlist = firstneigh[i];
jnum = numneigh[i];
@ -322,14 +319,17 @@ void MLIAPDescriptorSNAP::param_backward(int *map, NeighList* list,
delz = x[j][2] - ztmp;
rsq = delx*delx + dely*dely + delz*delz;
int jtype = type[j];
int jelem = jtype-1;
// element map not yet implemented
// int jelem = map[jtype];
const int jelem = jtype-1;
// if (chemflag)
// jelem = map[jtype];
if (rsq < cutsq[ielem][jelem]) {
snaptr->rij[ninside][0] = delx;
snaptr->rij[ninside][1] = dely;
snaptr->rij[ninside][2] = delz;
snaptr->inside[ninside] = j;
// element map not yet implemented
snaptr->wj[ninside] = wjelem[jelem];
snaptr->rcutij[ninside] = sqrt(cutsq[ielem][jelem]);
snaptr->element[ninside] = jelem; // element index for chem snap
@ -337,14 +337,13 @@ void MLIAPDescriptorSNAP::param_backward(int *map, NeighList* list,
}
}
if(chemflag)
if (chemflag)
snaptr->compute_ui(ninside, ielem);
else
snaptr->compute_ui(ninside, 0);
snaptr->compute_zi();
if(chemflag)
if (chemflag)
snaptr->compute_bi(ielem);
else
snaptr->compute_bi(0);
@ -352,12 +351,12 @@ void MLIAPDescriptorSNAP::param_backward(int *map, NeighList* list,
for (int jj = 0; jj < ninside; jj++) {
const int j = snaptr->inside[jj];
if(chemflag)
snaptr->compute_duidrj(snaptr->rij[jj], snaptr->wj[jj],
snaptr->rcutij[jj], jj, snaptr->element[jj]);
else
snaptr->compute_duidrj(snaptr->rij[jj], snaptr->wj[jj],
snaptr->rcutij[jj], jj, 0);
if(chemflag)
snaptr->compute_duidrj(snaptr->rij[jj], snaptr->wj[jj],
snaptr->rcutij[jj],jj, snaptr->element[jj]);
else
snaptr->compute_duidrj(snaptr->rij[jj], snaptr->wj[jj],
snaptr->rcutij[jj],jj, 0);
snaptr->compute_dbidrj();
@ -541,11 +540,9 @@ void MLIAPDescriptorSNAP::read_paramfile(char *paramfilename)
cut = 2.0*radelem[ielem]*rcutfac;
if (cut > cutmax) cutmax = cut;
cutsq[ielem][ielem] = cut*cut;
printf("ielem = %d ielem = %d cutsq[i][i] = %g\n",ielem, ielem, cutsq[ielem][ielem]);
for(int jelem = ielem+1; jelem < nelements; jelem++) {
cut = (radelem[ielem]+radelem[jelem])*rcutfac;
cutsq[ielem][jelem] = cutsq[jelem][ielem] = cut*cut;
printf("ielem = %d jelem = %d cutsq[i][j] = %g\n",ielem, jelem, cutsq[ielem][jelem]);
}
}
}

1
src/thermo.cpp
View File

@ -315,6 +315,7 @@ void Thermo::header()
std::string hdr;
for (int i = 0; i < nfield; i++) hdr += keyword[i] + std::string(" ");
hdr += "\n";
if (me == 0) utils::logmesg(lmp,hdr);
}