diff --git a/src/MLIAP/compute_mliap.cpp b/src/MLIAP/compute_mliap.cpp index 46c062ecaa..b32078bdc1 100644 --- a/src/MLIAP/compute_mliap.cpp +++ b/src/MLIAP/compute_mliap.cpp @@ -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++) { diff --git a/src/MLIAP/mliap_descriptor_snap.cpp b/src/MLIAP/mliap_descriptor_snap.cpp index ffb4ab65ff..bdbeaff43e 100644 --- a/src/MLIAP/mliap_descriptor_snap.cpp +++ b/src/MLIAP/mliap_descriptor_snap.cpp @@ -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; qidxu_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]); } } } diff --git a/src/thermo.cpp b/src/thermo.cpp index ea476dda04..339b641bb2 100644 --- a/src/thermo.cpp +++ b/src/thermo.cpp @@ -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); }