diff --git a/src/atom.cpp b/src/atom.cpp index 0249547253..2c2ebd911f 100644 --- a/src/atom.cpp +++ b/src/atom.cpp @@ -1181,7 +1181,7 @@ void Atom::data_atoms(int n, char *buf, tagint id_offset, tagint mol_offset, // convert atom coords from general triclinic to restricted triclinic - if (triclinic_general) domain->general_to_restricted(xdata); + if (triclinic_general) domain->general_to_restricted_coords(xdata); // apply shift if requested by read_data command diff --git a/src/create_atoms.cpp b/src/create_atoms.cpp index ac29387ff3..f7347f9ad1 100644 --- a/src/create_atoms.cpp +++ b/src/create_atoms.cpp @@ -1236,31 +1236,31 @@ void CreateAtoms::add_lattice() double point[3]; point[0] = bboxlo[0]; point[1] = bboxlo[1]; point[2] = bboxlo[2]; - domain->restricted_to_general(point); + domain->restricted_to_general_coords(point); domain->lattice->bbox(1, point[0], point[1], point[2], xmin, ymin, zmin, xmax, ymax, zmax); point[0] = bboxhi[0]; point[1] = bboxlo[1]; point[2] = bboxlo[2]; - domain->restricted_to_general(point); + domain->restricted_to_general_coords(point); domain->lattice->bbox(1, point[0], point[1], point[2], xmin, ymin, zmin, xmax, ymax, zmax); point[0] = bboxlo[0]; point[1] = bboxhi[1]; point[2] = bboxlo[2]; - domain->restricted_to_general(point); + domain->restricted_to_general_coords(point); domain->lattice->bbox(1, point[0], point[1], point[2], xmin, ymin, zmin, xmax, ymax, zmax); point[0] = bboxhi[0]; point[1] = bboxhi[1]; point[2] = bboxlo[2]; - domain->restricted_to_general(point); + domain->restricted_to_general_coords(point); domain->lattice->bbox(1, point[0], point[1], point[2], xmin, ymin, zmin, xmax, ymax, zmax); point[0] = bboxlo[0]; point[1] = bboxlo[1]; point[2] = bboxhi[2]; - domain->restricted_to_general(point); + domain->restricted_to_general_coords(point); domain->lattice->bbox(1, point[0], point[1], point[2], xmin, ymin, zmin, xmax, ymax, zmax); point[0] = bboxhi[0]; point[1] = bboxlo[1]; point[2] = bboxhi[2]; - domain->restricted_to_general(point); + domain->restricted_to_general_coords(point); domain->lattice->bbox(1, point[0], point[1], point[2], xmin, ymin, zmin, xmax, ymax, zmax); point[0] = bboxlo[0]; point[1] = bboxhi[1]; point[2] = bboxhi[2]; - domain->restricted_to_general(point); + domain->restricted_to_general_coords(point); domain->lattice->bbox(1, point[0], point[1], point[2], xmin, ymin, zmin, xmax, ymax, zmax); point[0] = bboxhi[0]; point[1] = bboxhi[1]; point[2] = bboxhi[2]; - domain->restricted_to_general(point); + domain->restricted_to_general_coords(point); domain->lattice->bbox(1, point[0], point[1], point[2], xmin, ymin, zmin, xmax, ymax, zmax); } @@ -1373,7 +1373,7 @@ void CreateAtoms::loop_lattice(int action) // convert from general to restricted triclinic coords - if (triclinic_general) domain->general_to_restricted(x); + if (triclinic_general) domain->general_to_restricted_coords(x); // if a region was specified, test if atom is in it diff --git a/src/create_box.cpp b/src/create_box.cpp index c534bf63f8..6a04744826 100644 --- a/src/create_box.cpp +++ b/src/create_box.cpp @@ -96,7 +96,7 @@ void CreateBox::command(int narg, char **arg) // setup general triclinic box (with no region) // read next box extent arguments to create ABC edge vectors + origin - // setup_general_triclinic() converts + // define_general_triclinic() converts // ABC edge vectors + origin to restricted triclinic } else if (triclinic_general) { @@ -146,9 +146,9 @@ void CreateBox::command(int narg, char **arg) cvec[1] = py - origin[1]; cvec[2] = pz - origin[2]; - // setup general triclinic box within Domain + // define general triclinic box within Domain - domain->setup_general_triclinic(avec,bvec,cvec,origin); + domain->define_general_triclinic(avec,bvec,cvec,origin); } // if molecular, zero out topology info diff --git a/src/domain.cpp b/src/domain.cpp index 29e0e2cd89..810814a80f 100644 --- a/src/domain.cpp +++ b/src/domain.cpp @@ -290,6 +290,29 @@ void Domain::set_global_box() boxhi_bound[1] = MAX(boxhi[1],boxhi[1]+yz); boxhi_bound[2] = boxhi[2]; } + + // update general triclinic box if defined + // reset ABC edge vectors from restricted triclinic box + // boxlo = lower left corner of general triclinic box + + if (triclinic_general) { + double aprime[3],bprime[3],cprime[3]; + + aprime[0] = boxhi[0] - boxlo[0]; + aprime[1] = aprime[2] = 0.0; + bprime[0] = xy; + bprime[1] = boxhi[1] - boxlo[1]; + bprime[2] = 0.0; + cprime[0] = xz; + cprime[1] = yz; + cprime[2] = boxhi[2] - boxlo[2]; + + // transform restricted A'B'C' to general triclinic A,B,C + + MathExtra::matvec(rotate_r2g,aprime,avec); + MathExtra::matvec(rotate_r2g,bprime,bvec); + MathExtra::matvec(rotate_r2g,cprime,cvec); + } } /* ---------------------------------------------------------------------- @@ -527,8 +550,8 @@ void Domain::reset_box() create rotation matrices for general <--> restricted transformations ------------------------------------------------------------------------- */ -void Domain::setup_general_triclinic(double *avec_caller, double *bvec_caller, - double *cvec_caller, double *origin_caller) +void Domain::define_general_triclinic(double *avec_caller, double *bvec_caller, + double *cvec_caller, double *origin_caller) { if (triclinic || triclinic_general) error->all(FLERR,"General triclinic box edge vectors are already set"); @@ -547,9 +570,9 @@ void Domain::setup_general_triclinic(double *avec_caller, double *bvec_caller, cvec[1] = cvec_caller[1]; cvec[2] = cvec_caller[2]; - gtri_origin[0] = origin_caller[0]; - gtri_origin[1] = origin_caller[1]; - gtri_origin[2] = origin_caller[2]; + boxlo[0] = origin_caller[0]; + boxlo[1] = origin_caller[1]; + boxlo[2] = origin_caller[2]; // error check on cvec for 2d systems @@ -623,9 +646,9 @@ void Domain::setup_general_triclinic(double *avec_caller, double *bvec_caller, } MathExtra::transpose3(rotate_g2r,rotate_r2g); - - // A',B',C' = transformation of A,B,C to restricted triclinic - + + // transform general ABC to restricted triclinic A'B'C' + double aprime[3],bprime[3],cprime[3]; MathExtra::matvec(rotate_g2r,avec,aprime); MathExtra::matvec(rotate_g2r,bvec,bprime); @@ -633,10 +656,6 @@ void Domain::setup_general_triclinic(double *avec_caller, double *bvec_caller, // set restricted triclinic boxlo, boxhi, and tilt factors - boxlo[0] = gtri_origin[0]; - boxlo[1] = gtri_origin[1]; - boxlo[2] = gtri_origin[2]; - boxhi[0] = boxlo[0] + aprime[0]; boxhi[1] = boxlo[1] + bprime[1]; boxhi[2] = boxlo[2] + cprime[2]; @@ -667,28 +686,72 @@ void Domain::setup_general_triclinic(double *avec_caller, double *bvec_caller, transform atom coords from general triclinic to restricted triclinic ------------------------------------------------------------------------- */ -void Domain::general_to_restricted(double *x) +void Domain::general_to_restricted_coords(double *x) { - double xnew[3]; + double xshift[3],xnew[3]; - MathExtra::matvec(rotate_g2r,x,xnew); - x[0] = xnew[0] + gtri_origin[0]; - x[1] = xnew[1] + gtri_origin[1]; - x[2] = xnew[2] + gtri_origin[2]; + xshift[0] = x[0] - boxlo[0]; + xshift[1] = x[1] - boxlo[1]; + xshift[2] = x[2] - boxlo[2]; + MathExtra::matvec(rotate_g2r,xshift,xnew); + x[0] = xnew[0] + boxlo[0]; + x[1] = xnew[1] + boxlo[1]; + x[2] = xnew[2] + boxlo[2]; } /* ---------------------------------------------------------------------- transform atom coords from restricted triclinic to general triclinic ------------------------------------------------------------------------- */ -void Domain::restricted_to_general(double *x) +void Domain::restricted_to_general_coords(double *x) { double xshift[3],xnew[3]; - xshift[0] = x[0] - gtri_origin[0]; - xshift[1] = x[1] - gtri_origin[1]; - xshift[2] = x[2] - gtri_origin[2]; + xshift[0] = x[0] - boxlo[0]; + xshift[1] = x[1] - boxlo[1]; + xshift[2] = x[2] - boxlo[2]; MathExtra::matvec(rotate_r2g,xshift,xnew); + x[0] = xnew[0] + boxlo[0]; + x[1] = xnew[1] + boxlo[1]; + x[2] = xnew[2] + boxlo[2]; +} + +void Domain::restricted_to_general_coords(double *x, double *xnew) +{ + double xshift[3]; + + xshift[0] = x[0] - boxlo[0]; + xshift[1] = x[1] - boxlo[1]; + xshift[2] = x[2] - boxlo[2]; + MathExtra::matvec(rotate_r2g,xshift,xnew); + xnew[0] += boxlo[0]; + xnew[1] += boxlo[1]; + xnew[2] += boxlo[2]; +} + +/* ---------------------------------------------------------------------- + transform atom vector from general triclinic to restricted triclinic +------------------------------------------------------------------------- */ + +void Domain::general_to_restricted_vector(double *x) +{ + double xnew[3]; + + MathExtra::matvec(rotate_g2r,x,xnew); + x[0] = xnew[0]; + x[1] = xnew[1]; + x[2] = xnew[2]; +} + +/* ---------------------------------------------------------------------- + transform atom vector from restricted triclinic to general triclinic +------------------------------------------------------------------------- */ + +void Domain::restricted_to_general_vector(double *x) +{ + double xnew[3]; + + MathExtra::matvec(rotate_r2g,x,xnew); x[0] = xnew[0]; x[1] = xnew[1]; x[2] = xnew[2]; diff --git a/src/domain.h b/src/domain.h index d09f4cf70d..7b306a9426 100644 --- a/src/domain.h +++ b/src/domain.h @@ -89,9 +89,9 @@ class Domain : protected Pointers { double h_rate[6], h_ratelo[3]; // rate of box size/shape change // general triclinic box - + // boxlo = lower left corner + double avec[3], bvec[3], cvec[3]; // ABC edge vectors of general triclinic box - double gtri_origin[3]; // origin of general triclinic box double rotate_g2r[3][3]; // rotation matrix from general --> restricted tri double rotate_r2g[3][3]; // rotation matrix from restricted --> general tri @@ -141,10 +141,13 @@ class Domain : protected Pointers { void image_flip(int, int, int); int ownatom(int, double *, imageint *, int); - void setup_general_triclinic(double *, double *, double *, double *); - void general_to_restricted(double *); - void restricted_to_general(double *); - + void define_general_triclinic(double *, double *, double *, double *); + void general_to_restricted_coords(double *); + void restricted_to_general_coords(double *); + void restricted_to_general_coords(double *, double *); + void general_to_restricted_vector(double *); + void restricted_to_general_vector(double *); + void set_lattice(int, char **); void add_region(int, char **); void delete_region(Region *); diff --git a/src/dump_custom.cpp b/src/dump_custom.cpp index 1e60295bbe..d96800554d 100644 --- a/src/dump_custom.cpp +++ b/src/dump_custom.cpp @@ -34,7 +34,7 @@ using namespace LAMMPS_NS; // customize by adding keyword -// also customize compute_atom_property.cpp +// also customize compute_property_atom.cpp enum{ID,MOL,PROC,PROCP1,TYPE,ELEMENT,MASS, X,Y,Z,XS,YS,ZS,XSTRI,YSTRI,ZSTRI,XU,YU,ZU,XUTRI,YUTRI,ZUTRI, @@ -88,6 +88,7 @@ DumpCustom::DumpCustom(LAMMPS *lmp, int narg, char **arg) : buffer_allow = 1; buffer_flag = 1; + triclinic_general = 0; nthresh = 0; nthreshlast = 0; @@ -294,10 +295,14 @@ void DumpCustom::init_style() if (binary && domain->triclinic == 0) header_choice = &DumpCustom::header_binary; + else if (binary && triclinic_general == 1) + header_choice = &DumpCustom::header_binary_triclinic_general; else if (binary && domain->triclinic == 1) header_choice = &DumpCustom::header_binary_triclinic; else if (!binary && domain->triclinic == 0) header_choice = &DumpCustom::header_item; + else if (!binary && triclinic_general == 1) + header_choice = &DumpCustom::header_item_triclinic_general; else if (!binary && domain->triclinic == 1) header_choice = &DumpCustom::header_item_triclinic; @@ -489,6 +494,31 @@ void DumpCustom::header_binary_triclinic(bigint ndump) /* ---------------------------------------------------------------------- */ +void DumpCustom::header_binary_triclinic_general(bigint ndump) +{ + header_format_binary(); + + fwrite(&update->ntimestep,sizeof(bigint),1,fp); + fwrite(&ndump,sizeof(bigint),1,fp); + fwrite(&domain->triclinic,sizeof(int),1,fp); + fwrite(&domain->triclinic_general,sizeof(int),1,fp); + fwrite(&domain->boundary[0][0],6*sizeof(int),1,fp); + fwrite(domain->avec,3*sizeof(double),1,fp); + fwrite(domain->bvec,3*sizeof(double),1,fp); + fwrite(domain->cvec,3*sizeof(double),1,fp); + fwrite(domain->boxlo,3*sizeof(double),1,fp); + fwrite(&nfield,sizeof(int),1,fp); + + header_unit_style_binary(); + header_time_binary(); + header_columns_binary(); + + if (multiproc) fwrite(&nclusterprocs,sizeof(int),1,fp); + else fwrite(&nprocs,sizeof(int),1,fp); +} + +/* ---------------------------------------------------------------------- */ + void DumpCustom::header_item(bigint ndump) { if (unit_flag && !unit_count) { @@ -535,6 +565,32 @@ void DumpCustom::header_item_triclinic(bigint ndump) /* ---------------------------------------------------------------------- */ +void DumpCustom::header_item_triclinic_general(bigint ndump) +{ + if (unit_flag && !unit_count) { + ++unit_count; + fmt::print(fp,"ITEM: UNITS\n{}\n",update->unit_style); + } + if (time_flag) fmt::print(fp,"ITEM: TIME\n{:.16}\n",compute_time()); + + fmt::print(fp,"ITEM: TIMESTEP\n{}\n" + "ITEM: NUMBER OF ATOMS\n{}\n", + update->ntimestep, ndump); + + fmt::print(fp,"ITEM: BOX BOUNDS abc origin {}\n" + "{:>1.16e} {:>1.16e} {:>1.16e} {:>1.16e}\n" + "{:>1.16e} {:>1.16e} {:>1.16e} {:>1.16e}\n" + "{:>1.16e} {:>1.16e} {:>1.16e} {:>1.16e}\n", + boundstr, + domain->avec[0],domain->avec[1],domain->avec[2],domain->boxlo[0], + domain->bvec[0],domain->bvec[1],domain->bvec[2],domain->boxlo[1], + domain->cvec[0],domain->cvec[1],domain->cvec[2],domain->boxlo[2]); + + fmt::print(fp,"ITEM: ATOMS {}\n",columns); +} + +/* ---------------------------------------------------------------------- */ + int DumpCustom::count() { int i; @@ -1286,13 +1342,16 @@ int DumpCustom::parse_fields(int narg, char **arg) vtype[iarg] = Dump::DOUBLE; } else if (strcmp(arg[iarg],"x") == 0) { - pack_choice[iarg] = &DumpCustom::pack_x; + if (triclinic_general) pack_choice[iarg] = &DumpCustom::pack_x_triclinic_general; + else pack_choice[iarg] = &DumpCustom::pack_x; vtype[iarg] = Dump::DOUBLE; } else if (strcmp(arg[iarg],"y") == 0) { - pack_choice[iarg] = &DumpCustom::pack_y; + if (triclinic_general) pack_choice[iarg] = &DumpCustom::pack_y_triclinic_general; + else pack_choice[iarg] = &DumpCustom::pack_y; vtype[iarg] = Dump::DOUBLE; - } else if (strcmp(arg[iarg],"z") == 0) { - pack_choice[iarg] = &DumpCustom::pack_z; + } else if (strcmp(arg[iarg],"z") == 0) { + if (triclinic_general) pack_choice[iarg] = &DumpCustom::pack_z_triclinic_general; + else pack_choice[iarg] = &DumpCustom::pack_z; vtype[iarg] = Dump::DOUBLE; } else if (strcmp(arg[iarg],"xs") == 0) { if (domain->triclinic) pack_choice[iarg] = &DumpCustom::pack_xs_triclinic; @@ -1307,15 +1366,18 @@ int DumpCustom::parse_fields(int narg, char **arg) else pack_choice[iarg] = &DumpCustom::pack_zs; vtype[iarg] = Dump::DOUBLE; } else if (strcmp(arg[iarg],"xu") == 0) { - if (domain->triclinic) pack_choice[iarg] = &DumpCustom::pack_xu_triclinic; + if (triclinic_general) pack_choice[iarg] = &DumpCustom::pack_xu_triclinic_general; + else if (domain->triclinic) pack_choice[iarg] = &DumpCustom::pack_xu_triclinic; else pack_choice[iarg] = &DumpCustom::pack_xu; vtype[iarg] = Dump::DOUBLE; } else if (strcmp(arg[iarg],"yu") == 0) { - if (domain->triclinic) pack_choice[iarg] = &DumpCustom::pack_yu_triclinic; + if (triclinic_general) pack_choice[iarg] = &DumpCustom::pack_yu_triclinic_general; + else if (domain->triclinic) pack_choice[iarg] = &DumpCustom::pack_yu_triclinic; else pack_choice[iarg] = &DumpCustom::pack_yu; vtype[iarg] = Dump::DOUBLE; } else if (strcmp(arg[iarg],"zu") == 0) { - if (domain->triclinic) pack_choice[iarg] = &DumpCustom::pack_zu_triclinic; + if (triclinic_general) pack_choice[iarg] = &DumpCustom::pack_zu_triclinic_general; + else if (domain->triclinic) pack_choice[iarg] = &DumpCustom::pack_zu_triclinic; else pack_choice[iarg] = &DumpCustom::pack_zu; vtype[iarg] = Dump::DOUBLE; } else if (strcmp(arg[iarg],"xsu") == 0) { @@ -1341,22 +1403,28 @@ int DumpCustom::parse_fields(int narg, char **arg) vtype[iarg] = Dump::INT; } else if (strcmp(arg[iarg],"vx") == 0) { - pack_choice[iarg] = &DumpCustom::pack_vx; + if (triclinic_general) pack_choice[iarg] = &DumpCustom::pack_vx_triclinic_general; + else pack_choice[iarg] = &DumpCustom::pack_vx; vtype[iarg] = Dump::DOUBLE; } else if (strcmp(arg[iarg],"vy") == 0) { - pack_choice[iarg] = &DumpCustom::pack_vy; + if (triclinic_general) pack_choice[iarg] = &DumpCustom::pack_vy_triclinic_general; + else pack_choice[iarg] = &DumpCustom::pack_vy; vtype[iarg] = Dump::DOUBLE; } else if (strcmp(arg[iarg],"vz") == 0) { - pack_choice[iarg] = &DumpCustom::pack_vz; + if (triclinic_general) pack_choice[iarg] = &DumpCustom::pack_vz_triclinic_general; + else pack_choice[iarg] = &DumpCustom::pack_vz; vtype[iarg] = Dump::DOUBLE; } else if (strcmp(arg[iarg],"fx") == 0) { - pack_choice[iarg] = &DumpCustom::pack_fx; + if (triclinic_general) pack_choice[iarg] = &DumpCustom::pack_fx_triclinic_general; + else pack_choice[iarg] = &DumpCustom::pack_fx; vtype[iarg] = Dump::DOUBLE; } else if (strcmp(arg[iarg],"fy") == 0) { - pack_choice[iarg] = &DumpCustom::pack_fy; + if (triclinic_general) pack_choice[iarg] = &DumpCustom::pack_fy_triclinic_general; + else pack_choice[iarg] = &DumpCustom::pack_fy; vtype[iarg] = Dump::DOUBLE; } else if (strcmp(arg[iarg],"fz") == 0) { - pack_choice[iarg] = &DumpCustom::pack_fz; + if (triclinic_general) pack_choice[iarg] = &DumpCustom::pack_fz_triclinic_general; + else pack_choice[iarg] = &DumpCustom::pack_fz; vtype[iarg] = Dump::DOUBLE; } else if (strcmp(arg[iarg],"q") == 0) { @@ -1395,6 +1463,7 @@ int DumpCustom::parse_fields(int narg, char **arg) error->all(FLERR,"Dumping an atom property that isn't allocated"); pack_choice[iarg] = &DumpCustom::pack_diameter; vtype[iarg] = Dump::DOUBLE; + } else if (strcmp(arg[iarg],"heatflow") == 0) { if (!atom->heatflow_flag) error->all(FLERR,"Dumping an atom property that isn't allocated"); @@ -1405,6 +1474,7 @@ int DumpCustom::parse_fields(int narg, char **arg) error->all(FLERR,"Dumping an atom property that isn't allocated"); pack_choice[iarg] = &DumpCustom::pack_temperature; vtype[iarg] = Dump::DOUBLE; + } else if (strcmp(arg[iarg],"omegax") == 0) { if (!atom->omega_flag) error->all(FLERR,"Dumping an atom property that isn't allocated"); @@ -1702,6 +1772,14 @@ int DumpCustom::modify_param(int narg, char **arg) return 2; } + if (strcmp(arg[0],"triclinic/general") == 0) { + if (narg < 2) error->all(FLERR,"Illegal dump_modify command"); + triclinic_general = utils::logical(FLERR,arg[1],false,lmp); + if (triclinic_general && !domain->triclinic_general) + error->all(FLERR,"Dump_modify triclinic/general invalid b/c simulation box is not"); + return 2; + } + if (strcmp(arg[0],"format") == 0) { if (narg < 2) utils::missing_cmd_args(FLERR, "dump_modify format", error); @@ -2283,6 +2361,48 @@ void DumpCustom::pack_z(int n) /* ---------------------------------------------------------------------- */ +void DumpCustom::pack_x_triclinic_general(int n) +{ + double **x = atom->x; + double xtri[3]; + + for (int i = 0; i < nchoose; i++) { + domain->restricted_to_general_coords(x[clist[i]],xtri); + buf[n] = xtri[0]; + n += size_one; + } +} + +/* ---------------------------------------------------------------------- */ + +void DumpCustom::pack_y_triclinic_general(int n) +{ + double **x = atom->x; + double xtri[3]; + + for (int i = 0; i < nchoose; i++) { + domain->restricted_to_general_coords(x[clist[i]],xtri); + buf[n] = xtri[1]; + n += size_one; + } +} + +/* ---------------------------------------------------------------------- */ + +void DumpCustom::pack_z_triclinic_general(int n) +{ + double **x = atom->x; + double xtri[3]; + + for (int i = 0; i < nchoose; i++) { + domain->restricted_to_general_coords(x[clist[i]],xtri); + buf[n] = xtri[2]; + n += size_one; + } +} + +/* ---------------------------------------------------------------------- */ + void DumpCustom::pack_xs(int n) { double **x = atom->x; @@ -2489,6 +2609,84 @@ void DumpCustom::pack_zu_triclinic(int n) /* ---------------------------------------------------------------------- */ +void DumpCustom::pack_xu_triclinic_general(int n) +{ + int j; + double **x = atom->x; + imageint *image = atom->image; + + double *h = domain->h; + double xu[3]; + int xbox,ybox,zbox; + + for (int i = 0; i < nchoose; i++) { + j = clist[i]; + xbox = (image[j] & IMGMASK) - IMGMAX; + ybox = (image[j] >> IMGBITS & IMGMASK) - IMGMAX; + zbox = (image[j] >> IMG2BITS) - IMGMAX; + xu[0] = x[j][0] + h[0]*xbox + h[5]*ybox + h[4]*zbox; + xu[1] = x[j][1] + h[1]*ybox + h[3]*zbox; + xu[2] = x[j][2] + h[2]*zbox; + domain->restricted_to_general_coords(xu); + buf[n] = xu[0]; + n += size_one; + } +} + +/* ---------------------------------------------------------------------- */ + +void DumpCustom::pack_yu_triclinic_general(int n) +{ + int j; + double **x = atom->x; + imageint *image = atom->image; + + double *h = domain->h; + double xu[3]; + int xbox,ybox,zbox; + + for (int i = 0; i < nchoose; i++) { + j = clist[i]; + xbox = (image[j] & IMGMASK) - IMGMAX; + ybox = (image[j] >> IMGBITS & IMGMASK) - IMGMAX; + zbox = (image[j] >> IMG2BITS) - IMGMAX; + xu[0] = x[j][0] + h[0]*xbox + h[5]*ybox + h[4]*zbox; + xu[1] = x[j][1] + h[1]*ybox + h[3]*zbox; + xu[2] = x[j][2] + h[2]*zbox; + domain->restricted_to_general_coords(xu); + buf[n] = xu[1]; + n += size_one; + } +} + +/* ---------------------------------------------------------------------- */ + +void DumpCustom::pack_zu_triclinic_general(int n) +{ + int j; + double **x = atom->x; + imageint *image = atom->image; + + double *h = domain->h; + double xu[3]; + int xbox,ybox,zbox; + + for (int i = 0; i < nchoose; i++) { + j = clist[i]; + xbox = (image[j] & IMGMASK) - IMGMAX; + ybox = (image[j] >> IMGBITS & IMGMASK) - IMGMAX; + zbox = (image[j] >> IMG2BITS) - IMGMAX; + xu[0] = x[j][0] + h[0]*xbox + h[5]*ybox + h[4]*zbox; + xu[1] = x[j][1] + h[1]*ybox + h[3]*zbox; + xu[2] = x[j][2] + h[2]*zbox; + domain->restricted_to_general_coords(xu); + buf[n] = xu[2]; + n += size_one; + } +} + +/* ---------------------------------------------------------------------- */ + void DumpCustom::pack_xsu(int n) { int j; @@ -2671,6 +2869,48 @@ void DumpCustom::pack_vz(int n) /* ---------------------------------------------------------------------- */ +void DumpCustom::pack_vx_triclinic_general(int n) +{ + double **v = atom->v; + double vtri[3]; + + for (int i = 0; i < nchoose; i++) { + domain->restricted_to_general_coords(v[clist[i]],vtri); + buf[n] = vtri[0]; + n += size_one; + } +} + +/* ---------------------------------------------------------------------- */ + +void DumpCustom::pack_vy_triclinic_general(int n) +{ + double **v = atom->v; + double vtri[3]; + + for (int i = 0; i < nchoose; i++) { + domain->restricted_to_general_coords(v[clist[i]],vtri); + buf[n] = vtri[1]; + n += size_one; + } +} + +/* ---------------------------------------------------------------------- */ + +void DumpCustom::pack_vz_triclinic_general(int n) +{ + double **v = atom->v; + double vtri[3]; + + for (int i = 0; i < nchoose; i++) { + domain->restricted_to_general_coords(v[clist[i]],vtri); + buf[n] = vtri[2]; + n += size_one; + } +} + +/* ---------------------------------------------------------------------- */ + void DumpCustom::pack_fx(int n) { double **f = atom->f; @@ -2707,6 +2947,48 @@ void DumpCustom::pack_fz(int n) /* ---------------------------------------------------------------------- */ +void DumpCustom::pack_fx_triclinic_general(int n) +{ + double **f = atom->f; + double ftri[3]; + + for (int i = 0; i < nchoose; i++) { + domain->restricted_to_general_coords(f[clist[i]],ftri); + buf[n] = ftri[0]; + n += size_one; + } +} + +/* ---------------------------------------------------------------------- */ + +void DumpCustom::pack_fy_triclinic_general(int n) +{ + double **f = atom->f; + double ftri[3]; + + for (int i = 0; i < nchoose; i++) { + domain->restricted_to_general_coords(f[clist[i]],ftri); + buf[n] = ftri[1]; + n += size_one; + } +} + +/* ---------------------------------------------------------------------- */ + +void DumpCustom::pack_fz_triclinic_general(int n) +{ + double **f = atom->f; + double ftri[3]; + + for (int i = 0; i < nchoose; i++) { + domain->restricted_to_general_coords(f[clist[i]],ftri); + buf[n] = ftri[2]; + n += size_one; + } +} + +/* ---------------------------------------------------------------------- */ + void DumpCustom::pack_q(int n) { double *q = atom->q; diff --git a/src/dump_custom.h b/src/dump_custom.h index 2b04944ec3..83adbe0d7f 100644 --- a/src/dump_custom.h +++ b/src/dump_custom.h @@ -37,6 +37,8 @@ class DumpCustom : public Dump { int nevery; // dump frequency for output char *idregion; // region ID, nullptr if no region + int triclinic_general; // set by dump_modify + int nthresh; // # of defined thresholds int nthreshlast; // # of defined thresholds with value = LAST // @@ -124,8 +126,10 @@ class DumpCustom : public Dump { FnPtrHeader header_choice; // ptr to write header functions void header_binary(bigint); void header_binary_triclinic(bigint); + void header_binary_triclinic_general(bigint); void header_item(bigint); void header_item_triclinic(bigint); + void header_item_triclinic_general(bigint); typedef void (DumpCustom::*FnPtrWrite)(int, double *); FnPtrWrite write_choice; // ptr to write data functions @@ -153,34 +157,52 @@ class DumpCustom : public Dump { void pack_x(int); void pack_y(int); void pack_z(int); + void pack_x_triclinic_general(int); + void pack_y_triclinic_general(int); + void pack_z_triclinic_general(int); + void pack_xs(int); void pack_ys(int); void pack_zs(int); void pack_xs_triclinic(int); void pack_ys_triclinic(int); void pack_zs_triclinic(int); + void pack_xu(int); void pack_yu(int); void pack_zu(int); void pack_xu_triclinic(int); void pack_yu_triclinic(int); void pack_zu_triclinic(int); + void pack_xu_triclinic_general(int); + void pack_yu_triclinic_general(int); + void pack_zu_triclinic_general(int); + void pack_xsu(int); void pack_ysu(int); void pack_zsu(int); void pack_xsu_triclinic(int); void pack_ysu_triclinic(int); void pack_zsu_triclinic(int); + void pack_ix(int); void pack_iy(int); void pack_iz(int); void pack_vx(int); void pack_vy(int); - void pack_vz(int); + void pack_vz(int); + void pack_vx_triclinic_general(int); + void pack_vy_triclinic_general(int); + void pack_vz_triclinic_general(int); + void pack_fx(int); void pack_fy(int); void pack_fz(int); + void pack_fx_triclinic_general(int); + void pack_fy_triclinic_general(int); + void pack_fz_triclinic_general(int); + void pack_q(int); void pack_mux(int); void pack_muy(int); @@ -188,6 +210,7 @@ class DumpCustom : public Dump { void pack_mu(int); void pack_radius(int); void pack_diameter(int); + void pack_heatflow(int); void pack_temperature(int); diff --git a/src/lmprestart.h b/src/lmprestart.h index e5e1da7edf..b3982ac8c1 100644 --- a/src/lmprestart.h +++ b/src/lmprestart.h @@ -29,6 +29,7 @@ enum{VERSION,SMALLINT,TAGINT,BIGINT, NDIHEDRALS,NDIHEDRALTYPES,DIHEDRAL_PER_ATOM, NIMPROPERS,NIMPROPERTYPES,IMPROPER_PER_ATOM, TRICLINIC,BOXLO,BOXHI,XY,XZ,YZ, + TRICLINIC_GENERAL,ROTATE_G2R,ROTATE_R2G, SPECIAL_LJ,SPECIAL_COUL, MASS,PAIR,BOND,ANGLE,DIHEDRAL,IMPROPER, MULTIPROC,MPIIO,PROCSPERFILE,PERPROC, diff --git a/src/pair_lj_cut.cpp b/src/pair_lj_cut.cpp index 0e4398f676..fae308c57f 100644 --- a/src/pair_lj_cut.cpp +++ b/src/pair_lj_cut.cpp @@ -111,8 +111,6 @@ void PairLJCut::compute(int eflag, int vflag) jtype = type[j]; if (rsq < cutsq[itype][jtype]) { - printf("AAA tags %d %d ij %d %d nlocal %d dist %g\n", - atom->tag[i],atom->tag[j],i,j,atom->nlocal,sqrt(rsq)); r2inv = 1.0 / rsq; r6inv = r2inv * r2inv * r2inv; forcelj = r6inv * (lj1[itype][jtype] * r6inv - lj2[itype][jtype]); diff --git a/src/read_data.cpp b/src/read_data.cpp index 75d1c23b75..c19250c2aa 100644 --- a/src/read_data.cpp +++ b/src/read_data.cpp @@ -480,7 +480,7 @@ void ReadData::command(int narg, char **arg) ellipsoidflag = lineflag = triflag = bodyflag = 0; xloxhi_flag = yloyhi_flag = zlozhi_flag = tilt_flag = 0; - avec_flag = bvec_flag = cvec_flag = gtri_origin_flag = 0; + avec_flag = bvec_flag = cvec_flag = abc_origin_flag = 0; // values in this data file @@ -495,7 +495,7 @@ void ReadData::command(int narg, char **arg) bvec[0] = bvec[1] = bvec[2] = 0.0; cvec[0] = cvec[1] = cvec[2] = 0.0; avec[0] = bvec[1] = cvec[2] = 1.0; - gtri_origin[0] = gtri_origin[1] = gtri_origin[2] = 0.0; + abc_origin[0] = abc_origin[1] = abc_origin[2] = 0.0; keyword[0] = '\0'; @@ -518,7 +518,7 @@ void ReadData::command(int narg, char **arg) // check if simulation box specified consistently - if (!avec_flag && !bvec_flag && !cvec_flag && !gtri_origin_flag) { + if (!avec_flag && !bvec_flag && !cvec_flag && !abc_origin_flag) { triclinic = triclinic_general = 0; if (tilt_flag) triclinic = 1; } else { @@ -579,11 +579,11 @@ void ReadData::command(int narg, char **arg) } // general triclinic box - // setup_general_triclinic() converts - // ABC edge vectors + gtri_origin to restricted triclinic + // define_general_triclinic() converts + // ABC edge vectors + abc_origin to restricted triclinic } else if (triclinic_general) { - domain->setup_general_triclinic(avec,bvec,cvec,gtri_origin); + domain->define_general_triclinic(avec,bvec,cvec,abc_origin); } } @@ -596,7 +596,7 @@ void ReadData::command(int narg, char **arg) // general triclinic // first data file must also be general triclinic - // avec,bvec,vec must match first data file + // avec,bvec,vec and origin must match first data file // shift not allowed if (triclinic_general) { @@ -609,8 +609,8 @@ void ReadData::command(int narg, char **arg) errflag = 1; if (cvec[0] != domain->cvec[0] || cvec[1] != domain->cvec[1] || cvec[2] != domain->cvec[2]) errflag = 1; - if (gtri_origin[0] != domain->gtri_origin[0] || gtri_origin[1] != domain->gtri_origin[1] || - gtri_origin[2] != domain->gtri_origin[2]) + if (abc_origin[0] != domain->boxlo[0] || abc_origin[1] != domain->boxlo[1] || + abc_origin[2] != domain->boxlo[2]) errflag = 1; if (errflag) error->all(FLERR,"Read_data subsequent file ABC vectors must be same as first file"); @@ -1432,11 +1432,11 @@ void ReadData::header(int firstpass) cvec[1] = utils::numeric(FLERR, words[1], false, lmp); cvec[2] = utils::numeric(FLERR, words[2], false, lmp); - } else if (utils::strmatch(line, "^\\s*\\f+\\s+\\f+\\s+\\f+\\s+\\gtri\\s+origin\\s")) { - gtri_origin_flag = 1; - gtri_origin[0] = utils::numeric(FLERR, words[0], false, lmp); - gtri_origin[1] = utils::numeric(FLERR, words[1], false, lmp); - gtri_origin[2] = utils::numeric(FLERR, words[2], false, lmp); + } else if (utils::strmatch(line, "^\\s*\\f+\\s+\\f+\\s+\\f+\\s+\\abc\\s+origin\\s")) { + abc_origin_flag = 1; + abc_origin[0] = utils::numeric(FLERR, words[0], false, lmp); + abc_origin[1] = utils::numeric(FLERR, words[1], false, lmp); + abc_origin[2] = utils::numeric(FLERR, words[2], false, lmp); } else break; diff --git a/src/read_data.h b/src/read_data.h index e3c269e685..28b277860a 100644 --- a/src/read_data.h +++ b/src/read_data.h @@ -63,10 +63,11 @@ class ReadData : public Command { double boxlo[3], boxhi[3]; double xy, xz, yz; - double avec[3], bvec[3], cvec[3], gtri_origin[3]; + double avec[3], bvec[3], cvec[3]; + double abc_origin[3]; int triclinic, triclinic_general; int xloxhi_flag, yloyhi_flag, zlozhi_flag, tilt_flag; - int avec_flag, bvec_flag, cvec_flag, gtri_origin_flag; + int avec_flag, bvec_flag, cvec_flag, abc_origin_flag; // optional args diff --git a/src/read_restart.cpp b/src/read_restart.cpp index 6925bd6096..372d1bcfe8 100644 --- a/src/read_restart.cpp +++ b/src/read_restart.cpp @@ -783,6 +783,13 @@ void ReadRestart::header() } else if (flag == YZ) { domain->yz = read_double(); + } else if (flag == TRICLINIC_GENERAL) { + domain->triclinic_general = read_int(); + } else if (flag == ROTATE_G2R) { + read_double_vec(9,&domain->rotate_g2r[0][0]); + } else if (flag == ROTATE_R2G) { + read_double_vec(9,&domain->rotate_r2g[0][0]); + } else if (flag == SPECIAL_LJ) { read_int(); read_double_vec(3,&force->special_lj[1]); diff --git a/src/write_data.cpp b/src/write_data.cpp index ad0c508d15..960a4ff1cf 100644 --- a/src/write_data.cpp +++ b/src/write_data.cpp @@ -227,7 +227,7 @@ void WriteData::write(const std::string &file) memory->create(xstore,nlocal,3,"write_data:xstore"); if (nlocal) memcpy(&xstore[0][0],&x[0][0],3*nlocal*sizeof(double)); for (int i = 0; i < nlocal; i++) - domain->restricted_to_general(x[i]); + domain->restricted_to_general_coords(x[i]); } if (natoms) atoms(); @@ -335,8 +335,8 @@ void WriteData::header() domain->avec[0],domain->avec[1],domain->avec[2], domain->bvec[0],domain->bvec[1],domain->bvec[2], domain->cvec[0],domain->cvec[1],domain->cvec[2]); - fmt::print(fp,"{} {} {} gtri origin\n", - domain->gtri_origin[0],domain->gtri_origin[1],domain->gtri_origin[2]); + fmt::print(fp,"{} {} {} abc origin\n", + domain->boxlo[0],domain->boxlo[1],domain->boxhi[2]); } } diff --git a/src/write_restart.cpp b/src/write_restart.cpp index a996532687..fb5b41eb02 100644 --- a/src/write_restart.cpp +++ b/src/write_restart.cpp @@ -448,6 +448,12 @@ void WriteRestart::header() write_double(XZ,domain->xz); write_double(YZ,domain->yz); + write_int(TRICLINIC_GENERAL,domain->triclinic_general); + if (domain->triclinic_general) { + write_double_vec(ROTATE_G2R,9,&domain->rotate_g2r[0][0]); + write_double_vec(ROTATE_R2G,9,&domain->rotate_r2g[0][0]); + } + write_double_vec(SPECIAL_LJ,3,&force->special_lj[1]); write_double_vec(SPECIAL_COUL,3,&force->special_coul[1]);