add general tri options to read and write data commands and Domain class
This commit is contained in:
Steve Plimpton
27
src/atom.cpp
27
src/atom.cpp
@ -1041,7 +1041,8 @@ void Atom::deallocate_topology()
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void Atom::data_atoms(int n, char *buf, tagint id_offset, tagint mol_offset,
|
||||
int type_offset, int shiftflag, double *shift,
|
||||
int type_offset, int triclinic_general,
|
||||
int shiftflag, double *shift,
|
||||
int labelflag, int *ilabel)
|
||||
{
|
||||
int xptr,iptr;
|
||||
@ -1053,6 +1054,7 @@ void Atom::data_atoms(int n, char *buf, tagint id_offset, tagint mol_offset,
|
||||
auto location = "Atoms section of data file";
|
||||
|
||||
// use the first line to detect and validate the number of words/tokens per line
|
||||
|
||||
next = strchr(buf,'\n');
|
||||
if (!next) error->all(FLERR, "Missing data in {}", location);
|
||||
*next = '\0';
|
||||
@ -1069,6 +1071,7 @@ void Atom::data_atoms(int n, char *buf, tagint id_offset, tagint mol_offset,
|
||||
error->all(FLERR,"Incorrect format in {}: {}", location, utils::trim(buf));
|
||||
|
||||
*next = '\n';
|
||||
|
||||
// set bounds for my proc
|
||||
// if periodic and I am lo/hi proc, adjust bounds by EPSILON
|
||||
// ensures all data atoms will be owned even with round-off
|
||||
@ -1143,11 +1146,19 @@ void Atom::data_atoms(int n, char *buf, tagint id_offset, tagint mol_offset,
|
||||
*next = '\0';
|
||||
auto values = Tokenizer(buf).as_vector();
|
||||
int nvalues = values.size();
|
||||
|
||||
// skip comment lines
|
||||
|
||||
if ((nvalues == 0) || (utils::strmatch(values[0],"^#.*"))) {
|
||||
// skip over empty or comment lines
|
||||
|
||||
// check that line has correct # of words
|
||||
|
||||
} else if ((nvalues < nwords) ||
|
||||
((nvalues > nwords) && (!utils::strmatch(values[nwords],"^#")))) {
|
||||
error->all(FLERR, "Incorrect format in {}: {}", location, utils::trim(buf));
|
||||
|
||||
// extract the atom coords and image flags (if they exist)
|
||||
|
||||
} else {
|
||||
int imx = 0, imy = 0, imz = 0;
|
||||
if (imageflag) {
|
||||
@ -1167,13 +1178,25 @@ void Atom::data_atoms(int n, char *buf, tagint id_offset, tagint mol_offset,
|
||||
xdata[0] = utils::numeric(FLERR,values[xptr],false,lmp);
|
||||
xdata[1] = utils::numeric(FLERR,values[xptr+1],false,lmp);
|
||||
xdata[2] = utils::numeric(FLERR,values[xptr+2],false,lmp);
|
||||
|
||||
// convert atom coords from general triclinic to restricted triclinic
|
||||
|
||||
if (triclinic_general) domain->general_to_restricted(xdata);
|
||||
|
||||
// apply shift if requested by read_data command
|
||||
|
||||
if (shiftflag) {
|
||||
xdata[0] += shift[0];
|
||||
xdata[1] += shift[1];
|
||||
xdata[2] += shift[2];
|
||||
}
|
||||
|
||||
// map atom into simulation box for periodic dimensions
|
||||
|
||||
domain->remap(xdata,imagedata);
|
||||
|
||||
// determine if this proc owns the atom
|
||||
|
||||
if (triclinic) {
|
||||
domain->x2lamda(xdata,lamda);
|
||||
coord = lamda;
|
||||
|
||||
@ -328,7 +328,8 @@ class Atom : protected Pointers {
|
||||
|
||||
void deallocate_topology();
|
||||
|
||||
void data_atoms(int, char *, tagint, tagint, int, int, double *, int, int *);
|
||||
void data_atoms(int, char *, tagint, tagint, int, int,
|
||||
int, double *, int, int *);
|
||||
void data_vels(int, char *, tagint);
|
||||
void data_bonds(int, char *, int *, tagint, int, int, int *);
|
||||
void data_angles(int, char *, int *, tagint, int, int, int *);
|
||||
|
||||
181
src/domain.cpp
181
src/domain.cpp
@ -28,6 +28,7 @@
|
||||
#include "force.h"
|
||||
#include "kspace.h"
|
||||
#include "lattice.h"
|
||||
#include "math_extra.h"
|
||||
#include "memory.h"
|
||||
#include "modify.h"
|
||||
#include "molecule.h"
|
||||
@ -41,6 +42,7 @@
|
||||
#include <cmath>
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
using namespace MathExtra;
|
||||
|
||||
#define BIG 1.0e20
|
||||
#define SMALL 1.0e-4
|
||||
@ -81,7 +83,7 @@ Domain::Domain(LAMMPS *lmp) : Pointers(lmp)
|
||||
minylo = minyhi = 0.0;
|
||||
minzlo = minzhi = 0.0;
|
||||
|
||||
triclinic = 0;
|
||||
triclinic = triclinic_general = 0;
|
||||
|
||||
boxlo[0] = boxlo[1] = boxlo[2] = -0.5;
|
||||
boxhi[0] = boxhi[1] = boxhi[2] = 0.5;
|
||||
@ -516,6 +518,180 @@ void Domain::reset_box()
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
define and store a general triclinic simulation box
|
||||
3 edge vectors of box = avec/bvec/cvec caller
|
||||
origin of edge vectors = origin_caller = lower left corner of box
|
||||
create mapping to restricted triclinic box
|
||||
set boxlo[3], boxhi[3] and 3 tilt factors
|
||||
create rotation matrices for general <--> restricted transformations
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void Domain::set_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");
|
||||
|
||||
triclinic = triclinic_general = 1;
|
||||
|
||||
avec[0] = avec_caller[0];
|
||||
avec[1] = avec_caller[1];
|
||||
avec[2] = avec_caller[2];
|
||||
|
||||
bvec[0] = bvec_caller[0];
|
||||
bvec[1] = bvec_caller[1];
|
||||
bvec[2] = bvec_caller[2];
|
||||
|
||||
cvec[0] = cvec_caller[0];
|
||||
cvec[1] = cvec_caller[1];
|
||||
cvec[2] = cvec_caller[2];
|
||||
|
||||
tri_origin[0] = origin_caller[0];
|
||||
tri_origin[1] = origin_caller[1];
|
||||
tri_origin[2] = origin_caller[2];
|
||||
|
||||
// error check for co-planar A,B,C
|
||||
|
||||
double abcross[3];
|
||||
MathExtra::cross3(avec,bvec,abcross);
|
||||
double dot = MathExtra::dot3(abcross,cvec);
|
||||
if (dot == 0.0)
|
||||
error->all(FLERR,"General triclinic box edge vectors are co-planar");
|
||||
|
||||
// quat1 = convert A into A' along x-axis
|
||||
// rot1 = unit vector to rotate A around
|
||||
// theta1 = angle of rotation calculated from
|
||||
// A dot xunit = Ax = |A| cos(theta1)
|
||||
|
||||
double rot1[3],quat1[4];
|
||||
double xaxis[3] = {1.0, 0.0, 0.0};
|
||||
|
||||
double avec_len = MathExtra::len3(avec);
|
||||
MathExtra::cross3(avec,xaxis,rot1);
|
||||
MathExtra::norm3(rot1);
|
||||
double theta1 = acos(avec[0]/avec_len);
|
||||
MathExtra::axisangle_to_quat(rot1,theta1,quat1);
|
||||
|
||||
// rotmat1 = rotation matrix associated with quat1
|
||||
|
||||
double rotmat1[3][3];
|
||||
MathExtra::quat_to_mat(quat1,rotmat1);
|
||||
|
||||
// B1 = rotation of B by quat1 rotation matrix
|
||||
|
||||
double bvec1[3];
|
||||
|
||||
MathExtra::matvec(rotmat1,bvec,bvec1);
|
||||
|
||||
// quat2 = rotation to convert B1 into B' in xy plane
|
||||
// Byz1 = projection of B1 into yz plane
|
||||
// rot2 = unit vector to rotate B1 around = -x axis
|
||||
// theta2 = angle of rotation calculated from
|
||||
// Byz1 dot yunit = B1y = |Byz1} cos(theta2)
|
||||
|
||||
double byzvec1[3],quat2[4];
|
||||
MathExtra::copy3(bvec1,byzvec1);
|
||||
byzvec1[0] = 0.0;
|
||||
double byzvec1_len = MathExtra::len3(byzvec1);
|
||||
double rot2[3] = {-1.0, 0.0, 0.0};
|
||||
double theta2 = acos(bvec1[1]/byzvec1_len);
|
||||
MathExtra::axisangle_to_quat(rot2,theta2,quat2);
|
||||
|
||||
// quat = product of quat2 * quat1 = transformation via single quat
|
||||
// rotate_g2r = general to restricted transformation matrix
|
||||
// rotate_r2g = restricted to general transformation matrix
|
||||
// if A x B not in direction of C, flip sign of z component of transform
|
||||
// done by flipping sign of 3rd row of rotate_g2r matrix
|
||||
|
||||
double quat[4];
|
||||
MathExtra::quatquat(quat2,quat1,quat);
|
||||
MathExtra::quat_to_mat(quat,rotate_g2r);
|
||||
|
||||
if (dot < 0.0) {
|
||||
rotate_g2r[2][0] = -rotate_g2r[2][0];
|
||||
rotate_g2r[2][1] = -rotate_g2r[2][1];
|
||||
rotate_g2r[2][2] = -rotate_g2r[2][2];
|
||||
}
|
||||
|
||||
MathExtra::transpose3(rotate_g2r,rotate_r2g);
|
||||
|
||||
// A',B',C' = transformation of A,B,C to restricted triclinic
|
||||
|
||||
double aprime[3],bprime[3],cprime[3];
|
||||
MathExtra::matvec(rotate_g2r,avec,aprime);
|
||||
MathExtra::matvec(rotate_g2r,bvec,bprime);
|
||||
MathExtra::matvec(rotate_g2r,cvec,cprime);
|
||||
|
||||
// set restricted triclinic boxlo, boxhi, and tilt factors
|
||||
|
||||
boxlo[0] = tri_origin[0];
|
||||
boxlo[1] = tri_origin[1];
|
||||
boxlo[2] = tri_origin[2];
|
||||
|
||||
boxhi[0] = boxlo[0] + aprime[0];
|
||||
boxhi[1] = boxlo[1] + bprime[1];
|
||||
boxhi[2] = boxlo[2] + cprime[2];
|
||||
|
||||
xy = bprime[1];
|
||||
xz = cprime[0];
|
||||
yz = cprime[1];
|
||||
|
||||
// debug
|
||||
|
||||
/*
|
||||
printf("Quat: %g %g %g %g\n",quat[0],quat[1],quat[2],quat[3]);
|
||||
double angle = 2.0*acos(quat[0]);
|
||||
printf("Theta: %g\n",angle);
|
||||
printf("Rotvec: %g %g %g\n",quat[1]/sin(0.5*angle),quat[2]/sin(0.5*angle),quat[3]/sin(0.5*angle));
|
||||
printf("Aprime: %g %g %g\n",aprime[0],aprime[1],aprime[2]);
|
||||
printf("Bprime: %g %g %g\n",bprime[0],bprime[1],bprime[2]);
|
||||
printf("Cprime: %g %g %g\n",cprime[0],cprime[1],cprime[2]);
|
||||
printf("Length A: %g %g\n",MathExtra::len3(avec),MathExtra::len3(aprime));
|
||||
printf("Length B: %g %g\n",MathExtra::len3(bvec),MathExtra::len3(bprime));
|
||||
printf("Length C: %g %g\n",MathExtra::len3(cvec),MathExtra::len3(cprime));
|
||||
|
||||
double coord1[3] = {0.5,0.0,0.0};
|
||||
double coord2[3] = {0.5,0.0,0.3};
|
||||
double newcoord[3];
|
||||
MathExtra::matvec(rotate_g2r,coord1,newcoord);
|
||||
printf("Atom1: %g %g %g\n",newcoord[0],newcoord[1],newcoord[2]);
|
||||
MathExtra::matvec(rotate_g2r,coord2,newcoord);
|
||||
printf("Atom2: %g %g %g\n",newcoord[0],newcoord[1],newcoord[2]);
|
||||
*/
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
transform one atom's coords from general triclinic to restricted triclinic
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void Domain::general_to_restricted(double *x)
|
||||
{
|
||||
double xnew[3];
|
||||
|
||||
MathExtra::matvec(rotate_g2r,x,xnew);
|
||||
x[0] = xnew[0] + tri_origin[0];
|
||||
x[1] = xnew[1] + tri_origin[1];
|
||||
x[2] = xnew[2] + tri_origin[2];
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
transform one atom's coords from restricted triclinic to general triclinic
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void Domain::restricted_to_general(double *x)
|
||||
{
|
||||
double xshift[3],xnew[3];
|
||||
|
||||
xshift[0] = x[0] - tri_origin[0];
|
||||
xshift[1] = x[1] - tri_origin[1];
|
||||
xshift[2] = x[2] - tri_origin[2];
|
||||
MathExtra::matvec(rotate_r2g,xshift,xnew);
|
||||
x[0] = xnew[0];
|
||||
x[1] = xnew[1];
|
||||
x[2] = xnew[2];
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
enforce PBC and modify box image flags for each atom
|
||||
called every reneighboring and by other commands that change atoms
|
||||
@ -676,9 +852,7 @@ int Domain::inside(double* x)
|
||||
lamda[1] < lo[1] || lamda[1] >= hi[1] ||
|
||||
lamda[2] < lo[2] || lamda[2] >= hi[2]) return 0;
|
||||
else return 1;
|
||||
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
@ -713,7 +887,6 @@ int Domain::inside_nonperiodic(double* x)
|
||||
if (!zperiodic && (lamda[2] < lo[2] || lamda[2] >= hi[2])) return 0;
|
||||
return 1;
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
|
||||
30
src/domain.h
30
src/domain.h
@ -39,7 +39,8 @@ class Domain : protected Pointers {
|
||||
// 2 = shrink-wrap non-periodic
|
||||
// 3 = shrink-wrap non-per w/ min
|
||||
|
||||
int triclinic; // 0 = orthog box, 1 = triclinic
|
||||
int triclinic; // 0 = orthog box, 1 = triclinic (restricted or general)
|
||||
int triclinic_general; // 1 if mapping to/from general triclinic is stored, 0 if not
|
||||
|
||||
// orthogonal box
|
||||
|
||||
@ -48,8 +49,8 @@ class Domain : protected Pointers {
|
||||
double prd[3]; // array form of dimensions
|
||||
double prd_half[3]; // array form of half dimensions
|
||||
|
||||
// triclinic box
|
||||
// xyzprd,xyzprd_half and prd,prd_half = same as if untilted
|
||||
// restricted triclinic box
|
||||
// xyzprd,xyzprd_half and prd,prd_half = same as if not tilted
|
||||
|
||||
double prd_lamda[3]; // lamda box = (1,1,1)
|
||||
double prd_half_lamda[3]; // lamda half box = (0.5,0.5,0.5)
|
||||
@ -58,14 +59,14 @@ class Domain : protected Pointers {
|
||||
|
||||
double boxlo[3], boxhi[3];
|
||||
|
||||
// triclinic box
|
||||
// boxlo/hi = same as if untilted
|
||||
// restricted triclinic box
|
||||
// boxlo/hi = same as if not tilted
|
||||
|
||||
double boxlo_lamda[3], boxhi_lamda[3]; // lamda box = (0,1)
|
||||
double boxlo_bound[3], boxhi_bound[3]; // bounding box of tilted domain
|
||||
double corners[8][3]; // 8 corner points
|
||||
|
||||
// orthogonal box & triclinic box
|
||||
// orthogonal box & restricted triclinic box
|
||||
|
||||
double minxlo, minxhi; // minimum size of global box
|
||||
double minylo, minyhi; // when shrink-wrapping
|
||||
@ -75,18 +76,27 @@ class Domain : protected Pointers {
|
||||
|
||||
double sublo[3], subhi[3]; // sub-box bounds on this proc
|
||||
|
||||
// triclinic box
|
||||
// restricted triclinic box
|
||||
// sublo/hi = undefined
|
||||
|
||||
double sublo_lamda[3], subhi_lamda[3]; // bounds of subbox in lamda
|
||||
|
||||
// triclinic box
|
||||
// restricted triclinic box
|
||||
|
||||
double xy, xz, yz; // 3 tilt factors
|
||||
double h[6], h_inv[6]; // shape matrix in Voigt ordering
|
||||
// Voigt = xx,yy,zz,yz,xz,xy
|
||||
double h_rate[6], h_ratelo[3]; // rate of box size/shape change
|
||||
|
||||
// general triclinic box
|
||||
|
||||
double avec[3], bvec[3], cvec[3]; // ABC edge vectors of general triclinic box
|
||||
double tri_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
|
||||
|
||||
// box flags
|
||||
|
||||
int box_change; // 1 if any of next 3 flags are set, else 0
|
||||
int box_change_size; // 1 if box size changes, 0 if not
|
||||
int box_change_shape; // 1 if box shape changes, 0 if not
|
||||
@ -131,6 +141,10 @@ class Domain : protected Pointers {
|
||||
void image_flip(int, int, int);
|
||||
int ownatom(int, double *, imageint *, int);
|
||||
|
||||
void set_general_triclinic(double *, double *, double *, double *);
|
||||
void general_to_restricted(double *);
|
||||
void restricted_to_general(double *);
|
||||
|
||||
void set_lattice(int, char **);
|
||||
void add_region(int, char **);
|
||||
void delete_region(Region *);
|
||||
|
||||
@ -47,6 +47,7 @@ inline void cross3(const double *v1, const double *v2, double *ans);
|
||||
|
||||
inline void zeromat3(double m[3][3]);
|
||||
inline void zeromat3(double **m);
|
||||
inline void transpose3(const double m[3][3], double ans[3][3]);
|
||||
|
||||
inline void col2mat(const double *ex, const double *ey, const double *ez, double m[3][3]);
|
||||
inline double det3(const double mat[3][3]);
|
||||
@ -763,6 +764,21 @@ inline void MathExtra::zeromat3(double **m)
|
||||
m[2][0] = m[2][1] = m[2][2] = 0.0;
|
||||
}
|
||||
|
||||
// transpose a matrix
|
||||
|
||||
inline void MathExtra::transpose3(const double m[3][3], double ans[3][3])
|
||||
{
|
||||
ans[0][0] = m[0][0];
|
||||
ans[0][1] = m[1][0];
|
||||
ans[0][2] = m[2][0];
|
||||
ans[1][0] = m[0][1];
|
||||
ans[1][1] = m[1][1];
|
||||
ans[1][2] = m[2][1];
|
||||
ans[2][0] = m[0][2];
|
||||
ans[2][1] = m[1][2];
|
||||
ans[2][2] = m[2][2];
|
||||
}
|
||||
|
||||
// add two matrices
|
||||
|
||||
inline void MathExtra::plus3(const double m[3][3], double **m2, double **ans)
|
||||
|
||||
@ -479,6 +479,9 @@ void ReadData::command(int narg, char **arg)
|
||||
bondflag = angleflag = dihedralflag = improperflag = 0;
|
||||
ellipsoidflag = lineflag = triflag = bodyflag = 0;
|
||||
|
||||
xloxhi_flag = yloyhi_flag = zlozhi_flag = tilt_flag = 0;
|
||||
avec_flag = bvec_flag = cvec_flag = 0;
|
||||
|
||||
// values in this data file
|
||||
|
||||
natoms = 0;
|
||||
@ -488,7 +491,12 @@ void ReadData::command(int narg, char **arg)
|
||||
|
||||
boxlo[0] = boxlo[1] = boxlo[2] = -0.5;
|
||||
boxhi[0] = boxhi[1] = boxhi[2] = 0.5;
|
||||
triclinic = 0;
|
||||
avec[0] = avec[1] = avec[2] = 0.0;
|
||||
bvec[0] = bvec[1] = bvec[2] = 0.0;
|
||||
cvec[0] = cvec[1] = cvec[2] = 0.0;
|
||||
avec[0] = bvec[1] = cvec[2] = 1.0;
|
||||
tri_origin[0] = tri_origin[1] = tri_origin[2] = 0.0;
|
||||
|
||||
keyword[0] = '\0';
|
||||
|
||||
nlocal_previous = atom->nlocal;
|
||||
@ -508,6 +516,17 @@ void ReadData::command(int narg, char **arg)
|
||||
|
||||
header(firstpass);
|
||||
|
||||
// check if simulation box specified consistently
|
||||
|
||||
if (!avec_flag && !bvec_flag && !cvec_flag) {
|
||||
triclinic = triclinic_general = 0;
|
||||
if (tilt_flag) triclinic = 1;
|
||||
} else {
|
||||
if (xloxhi_flag || yloyhi_flag || zlozhi_flag || tilt_flag)
|
||||
error->all(FLERR,"Read_data header cannot specify simulation box lo/hi/tilt and ABC vectors");
|
||||
triclinic = triclinic_general = 1;
|
||||
}
|
||||
|
||||
// problem setup using info from header
|
||||
// only done once, if firstpass and first data file
|
||||
// apply extra settings before grow(), even if no topology in file
|
||||
@ -536,6 +555,11 @@ void ReadData::command(int narg, char **arg)
|
||||
n = static_cast<int>(nbig);
|
||||
atom->avec->grow(n);
|
||||
|
||||
// setup simulation box
|
||||
// 3 options: orthogonal, restricted triclinic, general triclinic
|
||||
// for general triclinic: convect general ABC edge vectors to LAMMPS restricted triclinic
|
||||
|
||||
if (!triclinic_general) {
|
||||
domain->boxlo[0] = boxlo[0];
|
||||
domain->boxhi[0] = boxhi[0];
|
||||
domain->boxlo[1] = boxlo[1];
|
||||
@ -550,17 +574,50 @@ void ReadData::command(int narg, char **arg)
|
||||
domain->yz = yz;
|
||||
}
|
||||
|
||||
domain->print_box(" ");
|
||||
domain->set_initial_box();
|
||||
domain->set_global_box();
|
||||
comm->set_proc_grid();
|
||||
domain->set_local_box();
|
||||
} else if (triclinic_general) {
|
||||
domain->set_general_triclinic(avec,bvec,cvec,tri_origin);
|
||||
}
|
||||
}
|
||||
|
||||
// change simulation box to be union of existing box and new box + shift
|
||||
// only done if firstpass and not first data file
|
||||
// shift not allowed for general triclinic
|
||||
|
||||
if (firstpass && addflag != NONE) {
|
||||
|
||||
// general triclinic
|
||||
// first data file must also be general triclinic
|
||||
// avec,bvec,vec must match first data file
|
||||
// shift not allowed
|
||||
|
||||
if (triclinic_general) {
|
||||
if (!domain->triclinic_general)
|
||||
error->all(FLERR,"Read_data subsequent file cannot switch to general triclinic");
|
||||
int errflag = 0;
|
||||
if (avec[0] != domain->avec[0] || avec[1] != domain->avec[1] || avec[2] != domain->avec[2])
|
||||
errflag = 1;
|
||||
if (bvec[0] != domain->bvec[0] || bvec[1] != domain->bvec[1] || bvec[2] != domain->bvec[2])
|
||||
errflag = 1;
|
||||
if (cvec[0] != domain->cvec[0] || cvec[1] != domain->cvec[1] || cvec[2] != domain->cvec[2])
|
||||
errflag = 1;
|
||||
if (tri_origin[0] != domain->tri_origin[0] || tri_origin[1] != domain->tri_origin[1] ||
|
||||
tri_origin[2] != domain->tri_origin[2])
|
||||
errflag = 1;
|
||||
if (errflag)
|
||||
error->all(FLERR,"Read_data subsequent file ABC vectors must be same as first file");
|
||||
if (shift[0] != 0.0 || shift[1] != 0.0 || shift[2] != 0.0)
|
||||
error->all(FLERR,"Read_data subsequent file with ABC vectors cannot define shift");
|
||||
|
||||
// restricted triclinic
|
||||
// tilt factors must match first data file
|
||||
|
||||
} else if (triclinic) {
|
||||
if (!domain->triclinic || domain->triclinic_general)
|
||||
error->all(FLERR,"Read_data subsequent file cannot switch to restricted triclinic");
|
||||
if (xy != domain->xy || xz != domain->xz || yz != domain->yz)
|
||||
error->all(FLERR,"Read_data subsequent file tilt factors must be same as first file");
|
||||
}
|
||||
|
||||
double oldboxlo[3] = {domain->boxlo[0], domain->boxlo[1], domain->boxlo[2]};
|
||||
double oldboxhi[3] = {domain->boxhi[0], domain->boxhi[1], domain->boxhi[2]};
|
||||
domain->boxlo[0] = MIN(domain->boxlo[0], boxlo[0] + shift[0]);
|
||||
@ -570,7 +627,9 @@ void ReadData::command(int narg, char **arg)
|
||||
domain->boxlo[2] = MIN(domain->boxlo[2], boxlo[2] + shift[2]);
|
||||
domain->boxhi[2] = MAX(domain->boxhi[2], boxhi[2] + shift[2]);
|
||||
|
||||
// check of box has changed. If yes, warn about non-zero image flags
|
||||
// check if box has changed
|
||||
// if yes, warn about non-zero image flags
|
||||
|
||||
if ((oldboxlo[0] != domain->boxlo[0]) || (oldboxlo[1] != domain->boxlo[1]) ||
|
||||
(oldboxlo[2] != domain->boxlo[2]) || (oldboxhi[0] != domain->boxhi[0]) ||
|
||||
(oldboxhi[1] != domain->boxhi[1]) || (oldboxhi[2] != domain->boxhi[2])) {
|
||||
@ -588,18 +647,15 @@ void ReadData::command(int narg, char **arg)
|
||||
if ((flag_all > 0) && (comm->me == 0))
|
||||
error->warning(FLERR, "Non-zero image flags with growing box leads to bad coordinates");
|
||||
}
|
||||
}
|
||||
|
||||
// NOTE: not sure what to do about tilt value in subsequent data files
|
||||
//if (triclinic) {
|
||||
// domain->xy = xy; domain->xz = xz; domain->yz = yz;
|
||||
// }
|
||||
// setup simulation box and paritioning in Domain and Comm classes
|
||||
|
||||
domain->print_box(" ");
|
||||
domain->set_initial_box();
|
||||
domain->set_global_box();
|
||||
comm->set_proc_grid();
|
||||
domain->set_local_box();
|
||||
}
|
||||
|
||||
// allocate space for type label map
|
||||
|
||||
@ -608,8 +664,9 @@ void ReadData::command(int narg, char **arg)
|
||||
lmap = new LabelMap(lmp, ntypes, nbondtypes, nangletypes, ndihedraltypes, nimpropertypes);
|
||||
}
|
||||
|
||||
// rest of data file is Sections
|
||||
// read in any order, except where error checks
|
||||
// customize for new sections
|
||||
// read rest of file in free format
|
||||
|
||||
while (strlen(keyword)) {
|
||||
|
||||
@ -1140,7 +1197,8 @@ void ReadData::header(int firstpass)
|
||||
|
||||
// initialize type counts by the "extra" numbers so they get counted
|
||||
// in case the corresponding "types" line is missing and thus the extra
|
||||
// value will not be processed.
|
||||
// value will not be processed
|
||||
|
||||
if (addflag == NONE) {
|
||||
atom->ntypes = extra_atom_types;
|
||||
atom->nbondtypes = extra_bond_types;
|
||||
@ -1156,6 +1214,7 @@ void ReadData::header(int firstpass)
|
||||
if (eof == nullptr) error->one(FLERR, "Unexpected end of data file");
|
||||
|
||||
// check for units keyword in first line and print warning on mismatch
|
||||
|
||||
auto units = Tokenizer(utils::strfind(line, "units = \\w+")).as_vector();
|
||||
if (units.size() > 2) {
|
||||
if (units[2] != update->unit_style)
|
||||
@ -1326,23 +1385,47 @@ void ReadData::header(int firstpass)
|
||||
// so can treat differently for first vs subsequent data files
|
||||
|
||||
} else if (utils::strmatch(line, "^\\s*\\f+\\s+\\f+\\s+xlo\\s+xhi\\s")) {
|
||||
xloxhi_flag = 1;
|
||||
boxlo[0] = utils::numeric(FLERR, words[0], false, lmp);
|
||||
boxhi[0] = utils::numeric(FLERR, words[1], false, lmp);
|
||||
|
||||
} else if (utils::strmatch(line, "^\\s*\\f+\\s+\\f+\\s+ylo\\s+yhi\\s")) {
|
||||
yloyhi_flag = 1;
|
||||
boxlo[1] = utils::numeric(FLERR, words[0], false, lmp);
|
||||
boxhi[1] = utils::numeric(FLERR, words[1], false, lmp);
|
||||
|
||||
} else if (utils::strmatch(line, "^\\s*\\f+\\s+\\f+\\s+zlo\\s+zhi\\s")) {
|
||||
zlozhi_flag = 1;
|
||||
boxlo[2] = utils::numeric(FLERR, words[0], false, lmp);
|
||||
boxhi[2] = utils::numeric(FLERR, words[1], false, lmp);
|
||||
|
||||
} else if (utils::strmatch(line, "^\\s*\\f+\\s+\\f+\\s+\\f+\\s+xy\\s+xz\\s+yz\\s")) {
|
||||
triclinic = 1;
|
||||
tilt_flag = 1;
|
||||
xy = utils::numeric(FLERR, words[0], false, lmp);
|
||||
xz = utils::numeric(FLERR, words[1], false, lmp);
|
||||
yz = utils::numeric(FLERR, words[2], false, lmp);
|
||||
|
||||
} else if (utils::strmatch(line, "^\\s*\\f+\\s+\\f+\\s+\\f+\\s+\\f+\\s+\\avec\\s")) {
|
||||
avec_flag = 1;
|
||||
avec[0] = utils::numeric(FLERR, words[0], false, lmp);
|
||||
avec[1] = utils::numeric(FLERR, words[1], false, lmp);
|
||||
avec[2] = utils::numeric(FLERR, words[2], false, lmp);
|
||||
tri_origin[0] = utils::numeric(FLERR, words[3], false, lmp);
|
||||
|
||||
} else if (utils::strmatch(line, "^\\s*\\f+\\s+\\f+\\s+\\f+\\s+\\f+\\s+\\bvec\\s")) {
|
||||
bvec_flag = 1;
|
||||
bvec[0] = utils::numeric(FLERR, words[0], false, lmp);
|
||||
bvec[1] = utils::numeric(FLERR, words[1], false, lmp);
|
||||
bvec[2] = utils::numeric(FLERR, words[2], false, lmp);
|
||||
tri_origin[1] = utils::numeric(FLERR, words[3], false, lmp);
|
||||
|
||||
} else if (utils::strmatch(line, "^\\s*\\f+\\s+\\f+\\s+\\f+\\s+\\f+\\s+\\cvec\\s")) {
|
||||
cvec_flag = 1;
|
||||
cvec[0] = utils::numeric(FLERR, words[0], false, lmp);
|
||||
cvec[1] = utils::numeric(FLERR, words[1], false, lmp);
|
||||
cvec[2] = utils::numeric(FLERR, words[2], false, lmp);
|
||||
tri_origin[2] = utils::numeric(FLERR, words[3], false, lmp);
|
||||
|
||||
} else
|
||||
break;
|
||||
}
|
||||
@ -1407,8 +1490,8 @@ void ReadData::atoms()
|
||||
if (eof) error->all(FLERR, "Unexpected end of data file");
|
||||
if (tlabelflag && !lmap->is_complete(Atom::ATOM))
|
||||
error->all(FLERR, "Label map is incomplete: all types must be assigned a unique type label");
|
||||
atom->data_atoms(nchunk, buffer, id_offset, mol_offset, toffset, shiftflag, shift, tlabelflag,
|
||||
lmap->lmap2lmap.atom);
|
||||
atom->data_atoms(nchunk, buffer, id_offset, mol_offset, toffset, triclinic_general,
|
||||
shiftflag, shift, tlabelflag, lmap->lmap2lmap.atom);
|
||||
nread += nchunk;
|
||||
}
|
||||
|
||||
|
||||
@ -63,7 +63,10 @@ class ReadData : public Command {
|
||||
|
||||
double boxlo[3], boxhi[3];
|
||||
double xy, xz, yz;
|
||||
int triclinic;
|
||||
double avec[3], bvec[3], cvec[3], tri_origin[3];
|
||||
int triclinic, triclinic_general;
|
||||
int xloxhi_flag, yloyhi_flag, zlozhi_flag, tilt_flag;
|
||||
int avec_flag, bvec_flag, cvec_flag;
|
||||
|
||||
// optional args
|
||||
|
||||
|
||||
@ -72,8 +72,11 @@ void WriteData::command(int narg, char **arg)
|
||||
pairflag = II;
|
||||
coeffflag = 1;
|
||||
fixflag = 1;
|
||||
triclinic_general = 0;
|
||||
lmapflag = 1;
|
||||
// store current (default) setting since we may change it.
|
||||
|
||||
// store current (default) setting since we may change it
|
||||
|
||||
int types_style = atom->types_style;
|
||||
int noinit = 0;
|
||||
|
||||
@ -94,6 +97,9 @@ void WriteData::command(int narg, char **arg)
|
||||
} else if (strcmp(arg[iarg],"nofix") == 0) {
|
||||
fixflag = 0;
|
||||
iarg++;
|
||||
} else if (strcmp(arg[iarg],"triclinic") == 0) {
|
||||
triclinic_general = 1;
|
||||
iarg++;
|
||||
} else if (strcmp(arg[iarg],"nolabelmap") == 0) {
|
||||
lmapflag = 0;
|
||||
iarg++;
|
||||
@ -135,7 +141,9 @@ void WriteData::command(int narg, char **arg)
|
||||
}
|
||||
|
||||
write(file);
|
||||
|
||||
// restore saved setting
|
||||
|
||||
atom->types_style = types_style;
|
||||
}
|
||||
|
||||
@ -206,8 +214,32 @@ void WriteData::write(const std::string &file)
|
||||
}
|
||||
|
||||
// per atom info in Atoms and Velocities sections
|
||||
// if general triclinic:
|
||||
// save restricted triclinic atom coords
|
||||
// transform atom coords from restricted to general
|
||||
// restore save atom coords after output
|
||||
// NOTE: do same for velocities as well ?
|
||||
|
||||
double **xstore = nullptr;
|
||||
|
||||
if (triclinic_general) {
|
||||
double **x = atom->x;
|
||||
int nlocal = atom->nlocal;
|
||||
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]);
|
||||
}
|
||||
|
||||
if (natoms) atoms();
|
||||
|
||||
if (triclinic_general) {
|
||||
double **x = atom->x;
|
||||
int nlocal = atom->nlocal;
|
||||
if (nlocal) memcpy(&x[0][0],&xstore[0][0],3*nlocal*sizeof(double));
|
||||
memory->destroy(xstore);
|
||||
}
|
||||
|
||||
if (natoms) velocities();
|
||||
|
||||
// molecular topology info if defined
|
||||
@ -289,15 +321,22 @@ void WriteData::header()
|
||||
for (int m = 0; m < ifix->wd_header; m++)
|
||||
ifix->write_data_header(fp,m);
|
||||
|
||||
// box info
|
||||
// box info: orthogonal, restricted triclinic, or general triclinic (if requested)
|
||||
|
||||
auto box = fmt::format("\n{} {} xlo xhi\n{} {} ylo yhi\n{} {} zlo zhi\n",
|
||||
if (!triclinic_general) {
|
||||
fmt::print(fp,"\n{} {} xlo xhi\n{} {} ylo yhi\n{} {} zlo zhi\n",
|
||||
domain->boxlo[0],domain->boxhi[0],
|
||||
domain->boxlo[1],domain->boxhi[1],
|
||||
domain->boxlo[2],domain->boxhi[2]);
|
||||
if (domain->triclinic)
|
||||
box += fmt::format("{} {} {} xy xz yz\n",domain->xy,domain->xz,domain->yz);
|
||||
fputs(box.c_str(),fp);
|
||||
fmt::print(fp,"{} {} {} xy xz yz\n",domain->xy,domain->xz,domain->yz);
|
||||
|
||||
} else if (triclinic_general) {
|
||||
fmt::print(fp,"\n{} {} {} {} avec\n{} {} {} {} bvec\n{} {} {} {} cvec\n",
|
||||
domain->avec[0],domain->avec[1],domain->avec[2],domain->tri_origin[0],
|
||||
domain->bvec[0],domain->bvec[1],domain->bvec[2],domain->tri_origin[1],
|
||||
domain->cvec[0],domain->cvec[1],domain->cvec[2],domain->tri_origin[2]);
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
|
||||
@ -35,6 +35,7 @@ class WriteData : public Command {
|
||||
int pairflag;
|
||||
int coeffflag;
|
||||
int fixflag;
|
||||
int triclinic_general;
|
||||
int lmapflag;
|
||||
FILE *fp;
|
||||
bigint nbonds_local, nbonds;
|
||||
|
||||
Reference in New Issue
Block a user