diff --git a/src/USER-DPD/atom_vec_dpd.cpp b/src/USER-DPD/atom_vec_dpd.cpp index 0b1866e65c..92028b60e6 100644 --- a/src/USER-DPD/atom_vec_dpd.cpp +++ b/src/USER-DPD/atom_vec_dpd.cpp @@ -35,9 +35,9 @@ AtomVecDPD::AtomVecDPD(LAMMPS *lmp) : AtomVec(lmp) mass_type = 1; comm_x_only = comm_f_only = 0; // we communicate not only x forward but also dpdTheta - size_forward = 6; // 3 + dpdTheta + uCond + uMech + size_forward = 7; // 3 + dpdTheta + uCond + uMech + uChem size_reverse = 3; // 3 - size_border = 9; // 6 + dpdTheta + uCond + uMech + size_border = 12; // 6 + dpdTheta + uCond + uMech + uChem + uCG + uCGnew size_velocity = 3; size_data_atom = 6; // we read id + type + dpdTheta + x + y + z size_data_vel = 4; @@ -73,8 +73,12 @@ void AtomVecDPD::grow(int n) dpdTheta = memory->grow(atom->dpdTheta, nmax, "atom:dpdTheta"); uCond = memory->grow(atom->uCond,nmax,"atom:uCond"); uMech = memory->grow(atom->uMech,nmax,"atom:uMech"); + uChem = memory->grow(atom->uChem,nmax,"atom:uChem"); + uCG = memory->grow(atom->uCG,nmax,"atom:uCG"); + uCGnew = memory->grow(atom->uCGnew,nmax,"atom:uCGnew"); duCond = memory->grow(atom->duCond,nmax,"atom:duCond"); duMech = memory->grow(atom->duMech,nmax,"atom:duMech"); + duChem = memory->grow(atom->duChem,nmax,"atom:duChem"); if (atom->nextra_grow) for (int iextra = 0; iextra < atom->nextra_grow; iextra++) @@ -94,8 +98,12 @@ void AtomVecDPD::grow_reset() dpdTheta = atom->dpdTheta; uCond = atom->uCond; uMech = atom->uMech; + uChem = atom->uChem; + uCG = atom->uCG; + uCGnew = atom->uCGnew; duCond = atom->duCond; duMech = atom->duMech; + duChem = atom->duChem; } /* ---------------------------------------------------------------------- @@ -117,6 +125,9 @@ void AtomVecDPD::copy(int i, int j, int delflag) dpdTheta[j] = dpdTheta[i]; uCond[j] = uCond[i]; uMech[j] = uMech[i]; + uChem[j] = uChem[i]; + uCG[j] = uCG[i]; + uCGnew[j] = uCGnew[i]; if (atom->nextra_grow) for (int iextra = 0; iextra < atom->nextra_grow; iextra++) @@ -141,6 +152,7 @@ int AtomVecDPD::pack_comm(int n, int *list, double *buf, buf[m++] = dpdTheta[j]; buf[m++] = uCond[j]; buf[m++] = uMech[j]; + buf[m++] = uChem[j]; } } else { if (domain->triclinic == 0) { @@ -160,6 +172,7 @@ int AtomVecDPD::pack_comm(int n, int *list, double *buf, buf[m++] = dpdTheta[j]; buf[m++] = uCond[j]; buf[m++] = uMech[j]; + buf[m++] = uChem[j]; } } return m; @@ -186,6 +199,7 @@ int AtomVecDPD::pack_comm_vel(int n, int *list, double *buf, buf[m++] = dpdTheta[j]; buf[m++] = uCond[j]; buf[m++] = uMech[j]; + buf[m++] = uChem[j]; } } else { if (domain->triclinic == 0) { @@ -209,6 +223,7 @@ int AtomVecDPD::pack_comm_vel(int n, int *list, double *buf, buf[m++] = dpdTheta[j]; buf[m++] = uCond[j]; buf[m++] = uMech[j]; + buf[m++] = uChem[j]; } } else { dvx = pbc[0]*h_rate[0] + pbc[5]*h_rate[5] + pbc[4]*h_rate[4]; @@ -231,6 +246,7 @@ int AtomVecDPD::pack_comm_vel(int n, int *list, double *buf, buf[m++] = dpdTheta[j]; buf[m++] = uCond[j]; buf[m++] = uMech[j]; + buf[m++] = uChem[j]; } } } @@ -252,6 +268,7 @@ void AtomVecDPD::unpack_comm(int n, int first, double *buf) dpdTheta[i] = buf[m++]; uCond[i] = buf[m++]; uMech[i] = buf[m++]; + uChem[i] = buf[m++]; } } @@ -273,6 +290,7 @@ void AtomVecDPD::unpack_comm_vel(int n, int first, double *buf) dpdTheta[i] = buf[m++]; uCond[i] = buf[m++]; uMech[i] = buf[m++]; + uChem[i] = buf[m++]; } } @@ -328,6 +346,9 @@ int AtomVecDPD::pack_border(int n, int *list, double *buf, buf[m++] = dpdTheta[j]; buf[m++] = uCond[j]; buf[m++] = uMech[j]; + buf[m++] = uChem[j]; + buf[m++] = uCG[j]; + buf[m++] = uCGnew[j]; } } else { if (domain->triclinic == 0) { @@ -350,6 +371,9 @@ int AtomVecDPD::pack_border(int n, int *list, double *buf, buf[m++] = dpdTheta[j]; buf[m++] = uCond[j]; buf[m++] = uMech[j]; + buf[m++] = uChem[j]; + buf[m++] = uCG[j]; + buf[m++] = uCGnew[j]; } } @@ -384,6 +408,9 @@ int AtomVecDPD::pack_border_vel(int n, int *list, double *buf, buf[m++] = dpdTheta[j]; buf[m++] = uCond[j]; buf[m++] = uMech[j]; + buf[m++] = uChem[j]; + buf[m++] = uCG[j]; + buf[m++] = uCGnew[j]; } } else { if (domain->triclinic == 0) { @@ -410,6 +437,9 @@ int AtomVecDPD::pack_border_vel(int n, int *list, double *buf, buf[m++] = dpdTheta[j]; buf[m++] = uCond[j]; buf[m++] = uMech[j]; + buf[m++] = uChem[j]; + buf[m++] = uCG[j]; + buf[m++] = uCGnew[j]; } } else { dvx = pbc[0]*h_rate[0] + pbc[5]*h_rate[5] + pbc[4]*h_rate[4]; @@ -435,6 +465,9 @@ int AtomVecDPD::pack_border_vel(int n, int *list, double *buf, buf[m++] = dpdTheta[j]; buf[m++] = uCond[j]; buf[m++] = uMech[j]; + buf[m++] = uChem[j]; + buf[m++] = uCG[j]; + buf[m++] = uCGnew[j]; } } } @@ -458,6 +491,9 @@ int AtomVecDPD::pack_comm_hybrid(int n, int *list, double *buf) buf[m++] = dpdTheta[j]; buf[m++] = uCond[j]; buf[m++] = uMech[j]; + buf[m++] = uChem[j]; + buf[m++] = uCG[j]; + buf[m++] = uCGnew[j]; } return m; } @@ -474,6 +510,9 @@ int AtomVecDPD::pack_border_hybrid(int n, int *list, double *buf) buf[m++] = dpdTheta[j]; buf[m++] = uCond[j]; buf[m++] = uMech[j]; + buf[m++] = uChem[j]; + buf[m++] = uCG[j]; + buf[m++] = uCGnew[j]; } return m; } @@ -497,6 +536,9 @@ void AtomVecDPD::unpack_border(int n, int first, double *buf) dpdTheta[i] = buf[m++]; uCond[i] = buf[m++]; uMech[i] = buf[m++]; + uChem[i] = buf[m++]; + uCG[i] = buf[m++]; + uCGnew[i] = buf[m++]; } if (atom->nextra_border) @@ -527,6 +569,9 @@ void AtomVecDPD::unpack_border_vel(int n, int first, double *buf) dpdTheta[i] = buf[m++]; uCond[i] = buf[m++]; uMech[i] = buf[m++]; + uChem[i] = buf[m++]; + uCG[i] = buf[m++]; + uCGnew[i] = buf[m++]; } if (atom->nextra_border) @@ -547,6 +592,7 @@ int AtomVecDPD::unpack_comm_hybrid(int n, int first, double *buf) dpdTheta[i] = buf[m++]; uCond[i] = buf[m++]; uMech[i] = buf[m++]; + uChem[i] = buf[m++]; } return m; } @@ -563,6 +609,9 @@ int AtomVecDPD::unpack_border_hybrid(int n, int first, double *buf) dpdTheta[i] = buf[m++]; uCond[i] = buf[m++]; uMech[i] = buf[m++]; + uChem[i] = buf[m++]; + uCG[i] = buf[m++]; + uCGnew[i] = buf[m++]; } return m; } @@ -588,6 +637,9 @@ int AtomVecDPD::pack_exchange(int i, double *buf) buf[m++] = dpdTheta[i]; buf[m++] = uCond[i]; buf[m++] = uMech[i]; + buf[m++] = uChem[i]; + buf[m++] = uCG[i]; + buf[m++] = uCGnew[i]; if (atom->nextra_grow) for (int iextra = 0; iextra < atom->nextra_grow; iextra++) @@ -618,6 +670,9 @@ int AtomVecDPD::unpack_exchange(double *buf) dpdTheta[nlocal] = buf[m++]; uCond[nlocal] = buf[m++]; uMech[nlocal] = buf[m++]; + uChem[nlocal] = buf[m++]; + uCG[nlocal] = buf[m++]; + uCGnew[nlocal] = buf[m++]; if (atom->nextra_grow) for (int iextra = 0; iextra < atom->nextra_grow; iextra++) @@ -638,7 +693,7 @@ int AtomVecDPD::size_restart() int i; int nlocal = atom->nlocal; - int n = 14 * nlocal; // 11 + dpdTheta + uCond + uMech + int n = 15 * nlocal; // 11 + dpdTheta + uCond + uMech + uChem if (atom->nextra_restart) for (int iextra = 0; iextra < atom->nextra_restart; iextra++) @@ -670,6 +725,7 @@ int AtomVecDPD::pack_restart(int i, double *buf) buf[m++] = dpdTheta[i]; buf[m++] = uCond[i]; buf[m++] = uMech[i]; + buf[m++] = uChem[i]; if (atom->nextra_restart) for (int iextra = 0; iextra < atom->nextra_restart; iextra++) @@ -706,6 +762,7 @@ int AtomVecDPD::unpack_restart(double *buf) dpdTheta[nlocal] = buf[m++]; uCond[nlocal] = buf[m++]; uMech[nlocal] = buf[m++]; + uChem[nlocal] = buf[m++]; double **extra = atom->extra; if (atom->nextra_store) { @@ -742,8 +799,12 @@ void AtomVecDPD::create_atom(int itype, double *coord) dpdTheta[nlocal] = 0.0; uCond[nlocal] = 0.0; uMech[nlocal] = 0.0; + uChem[nlocal] = 0.0; + uCG[nlocal] = 0.0; + uCGnew[nlocal] = 0.0; duCond[nlocal] = 0.0; duMech[nlocal] = 0.0; + duChem[nlocal] = 0.0; atom->nlocal++; } @@ -781,6 +842,9 @@ void AtomVecDPD::data_atom(double *coord, tagint imagetmp, char **values) rho[nlocal] = 0.0; uCond[nlocal] = 0.0; uMech[nlocal] = 0.0; + uChem[nlocal] = 0.0; + uCG[nlocal] = 0.0; + uCGnew[nlocal] = 0.0; atom->nlocal++; } @@ -870,8 +934,12 @@ bigint AtomVecDPD::memory_usage() if (atom->memcheck("dpdTheta")) bytes += memory->usage(dpdTheta,nmax); if (atom->memcheck("uCond")) bytes += memory->usage(uCond,nmax); if (atom->memcheck("uMech")) bytes += memory->usage(uMech,nmax); + if (atom->memcheck("uChem")) bytes += memory->usage(uChem,nmax); + if (atom->memcheck("uCG")) bytes += memory->usage(uCG,nmax); + if (atom->memcheck("uCGnew")) bytes += memory->usage(uCGnew,nmax); if (atom->memcheck("duCond")) bytes += memory->usage(duCond,nmax); if (atom->memcheck("duMech")) bytes += memory->usage(duMech,nmax); + if (atom->memcheck("duChem")) bytes += memory->usage(duChem,nmax); return bytes; } diff --git a/src/USER-DPD/atom_vec_dpd.h b/src/USER-DPD/atom_vec_dpd.h index 753705e763..71ef48d053 100644 --- a/src/USER-DPD/atom_vec_dpd.h +++ b/src/USER-DPD/atom_vec_dpd.h @@ -58,8 +58,8 @@ class AtomVecDPD : public AtomVec { void write_data(FILE *, int, double **); int write_data_hybrid(FILE *, double *); bigint memory_usage(); - double *uCond,*uMech,*dpdTheta,*rho; - double *duCond,*duMech; + double *uCond,*uMech,*uChem,*uCG,*uCGnew,*rho,*dpdTheta; + double *duCond,*duMech,*duChem; protected: tagint *tag; diff --git a/src/USER-DPD/compute_dpd.cpp b/src/USER-DPD/compute_dpd.cpp index 87eab0baac..59716362c7 100644 --- a/src/USER-DPD/compute_dpd.cpp +++ b/src/USER-DPD/compute_dpd.cpp @@ -59,6 +59,7 @@ void ComputeDpd::compute_vector() double *uCond = atom->uCond; double *uMech = atom->uMech; + double *uChem = atom->uChem; double *dpdTheta = atom->dpdTheta; int nlocal = atom->nlocal; int *mask = atom->mask; @@ -72,7 +73,7 @@ void ComputeDpd::compute_vector() if (mask[i] & groupbit){ dpdU[0] += uCond[i]; dpdU[1] += uMech[i]; - dpdU[2] += uCond[i] + uMech[i]; + dpdU[2] += uChem[i]; dpdU[3] += 1.0 / dpdTheta[i]; dpdU[4]++; } diff --git a/src/USER-DPD/compute_dpd_atom.cpp b/src/USER-DPD/compute_dpd_atom.cpp index 8d944a292f..10f5d8203f 100644 --- a/src/USER-DPD/compute_dpd_atom.cpp +++ b/src/USER-DPD/compute_dpd_atom.cpp @@ -40,7 +40,7 @@ ComputeDpdAtom::ComputeDpdAtom(LAMMPS *lmp, int narg, char **arg) : if (narg != 3) error->all(FLERR,"Illegal compute dpd/atom command"); peratom_flag = 1; - size_peratom_cols = 3; + size_peratom_cols = 4; nmax = 0; dpdAtom = NULL; @@ -77,21 +77,23 @@ void ComputeDpdAtom::compute_peratom() double *uCond = atom->uCond; double *uMech = atom->uMech; + double *uChem = atom->uChem; double *dpdTheta = atom->dpdTheta; + int nlocal = atom->nlocal; int *mask = atom->mask; - if (atom->nmax > nmax) { + if (nlocal > nmax) { memory->destroy(dpdAtom); nmax = atom->nmax; memory->create(dpdAtom,nmax,size_peratom_cols,"dpd/atom:dpdAtom"); array_atom = dpdAtom; } - const int nlocal = atom->nlocal; for (int i = 0; i < nlocal; i++){ if (mask[i] & groupbit){ dpdAtom[i][0] = uCond[i]; dpdAtom[i][1] = uMech[i]; - dpdAtom[i][2] = dpdTheta[i]; + dpdAtom[i][2] = uChem[i]; + dpdAtom[i][3] = dpdTheta[i]; } } }