Homec4science
Diffusion lammps (lammps-icms) Browse

lammps/tools/replicate.c
60c7f745b9b7lammps-icms

Last Change
Raw File

replicate.c
View Options

/* replicate tool
read LAMMPS data file as input
create new LAMMPS data file as output
new file can be translated, dilated, mapped to periodic box
and/or replicated to create a larger system
create new molecule tags as needed
Syntax: replicate [options] < infile > outfile
Options:
-style atom_style
use the LAMMPS atom style that corresponds to this file
e.g. atomic or bond or angle or full or granular
will be used for reading and writing files
this option must be specified
for atom_style = hybrid, append subsequent styles
e.g. -style hybrid bond charge
-repeat nx ny nz
use the input data as a unit cell and replicate the unit
cell by (nx,ny,nz) in the +x, +y, +z directions
this creates new atoms, bonds, angles, etc
if not specified, nx = ny = nz = 1
-inbox xlo xhi ylo yhi zlo zhi
use these values for the bounding box of the input system
and for the unit cell used for replicating purposes
if not specified use the box bounds read in from infile
-outbox xlo xhi ylo yhi zlo zhi
use these values for the bounding box of the output system
if not specified the output bounding box is computed from the input
bounding box by holding the "lo" values constant and creating
new "hi" values by multiplying by the replicating factors
-unmap
unmap all input atom positions using input image flags
and inbox in a periodic sense
this enables molecular bonds to be replicated
-map
map all final atom positions into the output bounding box
in a periodic sense
works no matter how far away from the box the atoms are
this operation is independent of unmap option
-imageout
image flags will be written to outfile, else they won't be
if unmapping/mapping is done, image flags are modified accordingly
Notes:
works with "Velocities" entry by giving same velocity to
all replicated images of an atom, which is probably not ideal
can perform system dilation/contraction by using "-repeat 1 1 1"
and specifying a slightly larger/smaller output box than input
box.
can add/delete image flags by using -imageout
*/
#include <string.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
struct box {
double xlo,xhi,xsize;
double ylo,yhi,ysize;
double zlo,zhi,zsize;
};
#define MAXLINE 1000
main(int argc, char *argv[])
{
char line[MAXLINE]; /* strings for reading/parsing input file */
char linetoken[MAXLINE];
char *token,*ptr;
/* numbers of various quantities as read in */
int natoms,nbonds,nangles,ndihedrals,nimpropers;
int ntypes,nbondtypes,nangletypes,ndihedtypes,nimprotypes;
int nmolecules; /* total # of mols specified by mol tags */
/* vectors of various quantities as read in */
int *tag,*molecule,*type;
double *q,*radius,*density,*x,*y,*z,*mux,*muy,*muz;
int *imagex,*imagey,*imagez;
int *vtag;
double *vx,*vy,*vz,*vphix,*vphiy,*vphiz;
int *btype,*atype,*dtype,*itype;
int *bond1,*bond2,*angle1,*angle2,*angle3;
int *dihed1,*dihed2,*dihed3,*dihed4,*impro1,*impro2,*impro3,*impro4;
int atomstyle; /* atom styles and settings */
int style_angle,style_atomic,style_bond,style_charge,style_dipole;
int style_dpd,style_full,style_granular,style_molecular;
int charge,molecular,dipole;
int nx,ny,nz; /* replication factors */
int inbox,outbox; /* flags for whether in/out boxes were explicitly
specified */
int unmap; /* 0 = no unmapping, 1 = unmap from input box */
int map; /* 0 = no mapping, 1 = map to output box */
int imageout; /* 0/1 = no/yes output of image flags */
int imageflag; /* 0/1 = no/yes if image flags exist in input */
struct box in,out; /* input and output bounding boxes */
int ntotal; /* total replication factor = nx*ny*nz */
double newx,newy,newz; /* replicated atom position */
int size_atom_valid; /* values that should be on atom line */
int size_atom_actual; /* values that are on atom line */
int i,j,k,m,n,tmp,del,tag_offset,mol_offset; /* temp variables */
char **values;
char *gettoken(char *, int); /* function defs */
int count_words(char *);
void copyline(int);
void scale(struct box, int, int, int,
struct box, double *, double *, double *);
void pbc(struct box, double *, double *, double *, int *, int *, int *);
/* default input values */
atomstyle = 0;
style_angle = style_atomic = style_bond = style_charge = style_dipole =
style_dpd = style_full = style_granular = style_molecular = 0;
nx = ny = nz = 1;
inbox = outbox = 0;
unmap = 0;
map = 0;
imageout = 0;
/* read input options from command line
should do more error checking for missing args */
i = 1;
while (i < argc) {
if (!strcmp(argv[i],"-style")) {
atomstyle = 1;
if (strcmp(argv[i+1],"angle") == 0) style_angle = 1;
else if (strcmp(argv[i+1],"atomic") == 0) style_atomic = 1;
else if (strcmp(argv[i+1],"bond") == 0) style_bond = 1;
else if (strcmp(argv[i+1],"charge") == 0) style_charge = 1;
else if (strcmp(argv[i+1],"dipole") == 0) style_dipole = 1;
else if (strcmp(argv[i+1],"dpd") == 0) style_dpd = 1;
else if (strcmp(argv[i+1],"full") == 0) style_full = 1;
else if (strcmp(argv[i+1],"granular") == 0) style_granular = 1;
else if (strcmp(argv[i+1],"molecular") == 0) style_molecular = 1;
else if (strcmp(argv[i+1],"hybrid") == 0) {
n = i + 2;
while (n < argc && argv[n][0] != '-') {
if (strcmp(argv[n],"angle") == 0) style_angle = 1;
else if (strcmp(argv[n],"atomic") == 0) style_atomic = 1;
else if (strcmp(argv[n],"bond") == 0) style_bond = 1;
else if (strcmp(argv[n],"charge") == 0) style_charge = 1;
else if (strcmp(argv[n],"dipole") == 0) style_dipole = 1;
else if (strcmp(argv[n],"dpd") == 0) style_dpd = 1;
else if (strcmp(argv[n],"full") == 0) style_full = 1;
else if (strcmp(argv[n],"granular") == 0) style_granular = 1;
else if (strcmp(argv[n],"molecular") == 0) style_molecular = 1;
else {
fprintf(stderr,"Error with command-line arg style\n");
exit(1);
}
n++;
}
if (i == n - 2) {
fprintf(stderr,"Error with command-line arg style\n");
exit(1);
}
i = n - 2;
} else {
fprintf(stderr,"Error with command-line arg style\n");
exit(1);
}
i += 2;
} else if (!strcmp(argv[i],"-repeat")) {
sscanf(argv[i+1],"%d",&nx);
sscanf(argv[i+2],"%d",&ny);
sscanf(argv[i+3],"%d",&nz);
i += 4;
} else if (!strcmp(argv[i],"-inbox")) {
inbox = 1;
sscanf(argv[i+1],"%lg",&in.xlo);
sscanf(argv[i+2],"%lg",&in.xhi);
sscanf(argv[i+3],"%lg",&in.ylo);
sscanf(argv[i+4],"%lg",&in.yhi);
sscanf(argv[i+5],"%lg",&in.zlo);
sscanf(argv[i+6],"%lg",&in.zhi);
i += 7;
} else if (!strcmp(argv[i],"-outbox")) {
outbox = 1;
sscanf(argv[i+1],"%lg",&out.xlo);
sscanf(argv[i+2],"%lg",&out.xhi);
sscanf(argv[i+3],"%lg",&out.ylo);
sscanf(argv[i+4],"%lg",&out.yhi);
sscanf(argv[i+5],"%lg",&out.zlo);
sscanf(argv[i+6],"%lg",&out.zhi);
i += 7;
} else if (!strcmp(argv[i],"-unmap")) {
unmap = 1;
i += 1;
} else if (!strcmp(argv[i],"-map")) {
map = 1;
i += 1;
} else if (!strcmp(argv[i],"-imageout")) {
imageout = 1;
i += 1;
} else {
fprintf(stderr,"Syntax error: replicate [options] < infile > outfile\n");
exit(1);
}
}
/* set low-level info from style flags */
charge = 0;
if (style_charge || style_full || style_dipole) charge = 1;
molecular = 0;
if (style_angle || style_bond || style_full || style_molecular)
molecular = 1;
dipole = 0;
if (style_dipole) dipole = 1;
size_atom_valid = 5;
if (charge) size_atom_valid += 1;
if (style_dipole) size_atom_valid += 3;
if (style_granular) size_atom_valid += 2;
if (molecular) size_atom_valid += 1;
/* error checks */
if (atomstyle == 0) {
fprintf(stderr,"ERROR: must use -style to set atom style\n");
exit(1);
}
/* ntotal = total replication factor */
ntotal = nx*ny*nz;
/* defaults */
natoms = nbonds = nangles = ndihedrals = nimpropers = 0;
ntypes = nbondtypes = nangletypes = ndihedtypes = nimprotypes = 0;
/* read/write header */
copyline(1); // echo 1st line
while (1) {
fgets(line,MAXLINE,stdin);
// trim anything from '#' onward
// if line is blank, continue
if (ptr = strchr(line,'#')) *ptr = '\0';
if (strspn(line," \t\n\r") == strlen(line)) {
printf("\n");
continue;
}
if (strstr(line,"atoms")) {
sscanf(line,"%d",&natoms);
printf("%d atoms\n",natoms*ntotal);
} else if (strstr(line,"bonds")) {
sscanf(line,"%d",&nbonds);
printf("%d bonds\n",nbonds*ntotal);
} else if (strstr(line,"angles")) {
sscanf(line,"%d",&nangles);
printf("%d angles\n",nangles*ntotal);
} else if (strstr(line,"dihedrals")) {
sscanf(line,"%d",&ndihedrals);
printf("%d dihedrals\n",ndihedrals*ntotal);
} else if (strstr(line,"impropers")) {
sscanf(line,"%d",&nimpropers);
printf("%d impropers\n",nimpropers*ntotal);
} else if (strstr(line,"atom types")) {
sscanf(line,"%d",&ntypes);
printf("%d atom types\n",ntypes);
} else if (strstr(line,"bond types")) {
sscanf(line,"%d",&nbondtypes);
printf("%d bond types\n",nbondtypes);
} else if (strstr(line,"angle types")) {
sscanf(line,"%d",&nangletypes);
printf("%d angle types\n",nangletypes);
} else if (strstr(line,"dihedral types")) {
sscanf(line,"%d",&ndihedtypes);
printf("%d dihedral types\n",ndihedtypes);
} else if (strstr(line,"improper types")) {
sscanf(line,"%d",&nimprotypes);
printf("%d improper types\n",nimprotypes);
} else if (strstr(line,"xlo xhi")) {
if (!inbox) sscanf(line,"%lg %lg",&in.xlo,&in.xhi);
in.xsize = in.xhi - in.xlo;
if (!outbox) {
out.xlo = in.xlo;
if (nx > 1) out.xhi = out.xlo + (in.xhi-in.xlo)*nx;
else out.xhi = in.xhi;
}
out.xsize = out.xhi - out.xlo;
printf("%g %g xlo xhi\n",out.xlo,out.xhi);
} else if (strstr(line,"ylo yhi")) {
if (!inbox) sscanf(line,"%lg %lg",&in.ylo,&in.yhi);
in.ysize = in.yhi - in.ylo;
if (!outbox) {
out.ylo = in.ylo;
if (ny > 1) out.yhi = out.ylo + (in.yhi-in.ylo)*ny;
else out.yhi = in.yhi;
}
out.ysize = out.yhi - out.ylo;
printf("%g %g ylo yhi\n",out.ylo,out.yhi);
} else if (strstr(line,"zlo zhi")) {
if (!inbox) sscanf(line,"%lg %lg",&in.zlo,&in.zhi);
in.zsize = in.zhi - in.zlo;
if (!outbox) {
out.zlo = in.zlo;
if (nz > 1) out.zhi = out.zlo + (in.zhi-in.zlo)*nz;
else out.zhi = in.zhi;
}
out.zsize = out.zhi - out.zlo;
printf("%g %g zlo zhi\n",out.zlo,out.zhi);
} else break;
}
/* malloc space for input atoms and topology */
tag = (int *) malloc(natoms*sizeof(int));
type = (int *) malloc(natoms*sizeof(int));
x = (double *) malloc(natoms*sizeof(double));
y = (double *) malloc(natoms*sizeof(double));
z = (double *) malloc(natoms*sizeof(double));
imagex = (int *) malloc(natoms*sizeof(int));
imagey = (int *) malloc(natoms*sizeof(int));
imagez = (int *) malloc(natoms*sizeof(int));
vtag = (int *) malloc(natoms*sizeof(int));
vx = (double *) malloc(natoms*sizeof(double));
vy = (double *) malloc(natoms*sizeof(double));
vz = (double *) malloc(natoms*sizeof(double));
if (tag == NULL || type == NULL ||
x == NULL || y == NULL || z == NULL ||
imagex == NULL || imagey == NULL || imagez == NULL ||
vtag == NULL || vx == NULL || vy == NULL || vz == NULL) {
fprintf(stderr,"Error in atom malloc - no space\n");
exit(1);
}
if (charge) {
q = (double *) malloc(natoms*sizeof(double));
if (q == NULL) {
fprintf(stderr,"Error in atom malloc - no space\n");
exit(1);
}
}
if (molecular) {
molecule = (int *) malloc(natoms*sizeof(int));
if (molecule == NULL) {
fprintf(stderr,"Error in atom malloc - no space\n");
exit(1);
}
}
if (dipole) {
mux = (double *) malloc(natoms*sizeof(double));
muy = (double *) malloc(natoms*sizeof(double));
muz = (double *) malloc(natoms*sizeof(double));
if (mux == NULL || muy == NULL || muz == NULL) {
fprintf(stderr,"Error in atom malloc - no space\n");
exit(1);
}
}
if (style_granular) {
radius = (double *) malloc(natoms*sizeof(double));
density = (double *) malloc(natoms*sizeof(double));
vphix = (double *) malloc(natoms*sizeof(double));
vphiy = (double *) malloc(natoms*sizeof(double));
vphiz = (double *) malloc(natoms*sizeof(double));
if (radius == NULL || density == NULL ||
vphix == NULL || vphiy == NULL || vphiz == NULL) {
fprintf(stderr,"Error in atom malloc - no space\n");
exit(1);
}
}
if (nbonds) {
btype = (int *) malloc(nbonds*sizeof(int));
bond1 = (int *) malloc(nbonds*sizeof(int));
bond2 = (int *) malloc(nbonds*sizeof(int));
if (btype == NULL || bond1 == NULL || bond2 == NULL) {
fprintf(stderr,"Error in bond malloc - no space\n");
exit(1);
}
}
if (nangles) {
atype = (int *) malloc(nangles*sizeof(int));
angle1 = (int *) malloc(nangles*sizeof(int));
angle2 = (int *) malloc(nangles*sizeof(int));
angle3 = (int *) malloc(nangles*sizeof(int));
if (atype == NULL || angle1 == NULL || angle2 == NULL || angle3 == NULL) {
fprintf(stderr,"Error in angle malloc - no space\n");
exit(1);
}
}
if (ndihedrals) {
dtype = (int *) malloc(ndihedrals*sizeof(int));
dihed1 = (int *) malloc(ndihedrals*sizeof(int));
dihed2 = (int *) malloc(ndihedrals*sizeof(int));
dihed3 = (int *) malloc(ndihedrals*sizeof(int));
dihed4 = (int *) malloc(ndihedrals*sizeof(int));
if (dtype == NULL || dihed1 == NULL || dihed2 == NULL ||
dihed3 == NULL || dihed4 == NULL) {
fprintf(stderr,"Error in dihedral malloc - no space\n");
exit(1);
}
}
if (nimpropers) {
itype = (int *) malloc(nimpropers*sizeof(int));
impro1 = (int *) malloc(nimpropers*sizeof(int));
impro2 = (int *) malloc(nimpropers*sizeof(int));
impro3 = (int *) malloc(nimpropers*sizeof(int));
impro4 = (int *) malloc(nimpropers*sizeof(int));
if (itype == NULL || impro1 == NULL || impro2 == NULL ||
impro3 == NULL || impro4 == NULL) {
fprintf(stderr,"Error in improper malloc - no space\n");
exit(1);
}
}
/* read identifier strings one by one in free-form part of data file */
token = gettoken(line,1);
while (token) {
/* read atoms */
if (!strcmp(token,"Atoms")) {
for (i = 0; i < natoms; i++) {
fgets(line,MAXLINE,stdin);
if (i == 0) {
size_atom_actual = count_words(line);
if (size_atom_actual == size_atom_valid) imageflag = 0;
else if (size_atom_actual == size_atom_valid + 3) imageflag = 1;
else {
fprintf(stderr,"Error in atom lines\n");
exit(1);
}
if (unmap == 1 && imageflag == 0) {
fprintf(stderr,
"If unmap is requested, file must have image flags\n");
exit(1);
}
values = (char **) malloc(size_atom_actual*sizeof(char *));
}
values[0] = strtok(line," \t\n");
for (m = 1; m < size_atom_actual; m++)
values[m] = strtok(NULL," \t\n");
m = 0;
tag[i] = atoi(values[m++]);
if (molecular) molecule[i] = atoi(values[m++]);
type[i] = atoi(values[m++]);
if (charge) q[i] = atof(values[m++]);
if (style_granular) {
radius[i] = 0.5 * atof(values[m++]);
density[i] = atof(values[m++]);
}
x[i] = atof(values[m++]);
y[i] = atof(values[m++]);
z[i] = atof(values[m++]);
if (style_dipole) {
mux[i] = atof(values[m++]);
muy[i] = atof(values[m++]);
muz[i] = atof(values[m++]);
}
if (imageflag) {
imagex[i] = atoi(values[m++]);
imagey[i] = atoi(values[m++]);
imagez[i] = atoi(values[m++]);
} else imagex[i] = imagey[i] = imagez[i] = 0;
}
/* unmap atoms if requested */
if (unmap) {
for (i = 0; i < natoms; i++) {
x[i] = x[i] + imagex[i]*in.xsize;
y[i] = y[i] + imagey[i]*in.ysize;
z[i] = z[i] + imagez[i]*in.zsize;
imagex[i] = imagey[i] = imagez[i] = 0;
}
}
/* find number of molecules */
nmolecules = 0;
if (molecular)
for (i = 0; i < natoms; i++)
if (molecule[i] > nmolecules) nmolecules = molecule[i];
/* replicate set of N atoms as many times as requested
generate new tags, mol number, coords, image flags as needed */
for (k = 0; k < nz; k++) {
for (j = 0; j < ny; j++) {
for (i = 0; i < nx; i++) {
for (m = 0; m < natoms; m++) {
newx = x[m] + i*in.xsize;
newy = y[m] + j*in.ysize;
newz = z[m] + k*in.zsize;
if (outbox) scale(in,nx,ny,nz,out,&newx,&newy,&newz);
if (map) pbc(out,&newx,&newy,&newz,
&imagex[m],&imagey[m],&imagez[m]);
tag_offset = k*ny*nx*natoms + j*nx*natoms + i*natoms;
mol_offset = k*ny*nx*nmolecules + j*nx*nmolecules + i*nmolecules;
line[0] = '\0';
sprintf(&line[strlen(line)],"%d",tag[m]+tag_offset);
if (molecular) sprintf(&line[strlen(line)]," %d",
molecule[m]+mol_offset);
sprintf(&line[strlen(line)]," %d",type[m]);
if (charge) sprintf(&line[strlen(line)]," %g",q[m]);
if (style_granular) sprintf(&line[strlen(line)]," %g %g",
radius[m],density[m]);
sprintf(&line[strlen(line)]," %g %g %g",newx,newy,newz);
if (style_dipole) sprintf(&line[strlen(line)]," %g %g %g",
mux[m],muy[m],muz[m]);
if (imageout) sprintf(&line[strlen(line)]," %d %d %d",
imagex[m],imagey[m],imagez[m]);
printf("%s\n",line);
}
}
}
}
}
/* read velocities and replicate */
else if (!strcmp(token,"Velocities")) {
for (i = 0; i < natoms; i++) {
fgets(line,MAXLINE,stdin);
if (!style_granular)
sscanf(line,"%d %lg %lg %lg",&vtag[i],&vx[i],&vy[i],&vz[i]);
else
sscanf(line,"%d %lg %lg %lg %lg %lg %lg",&vtag[i],
&vx[i],&vy[i],&vz[i],&vphix[i],&vphiy[i],&vphiz[i]);
}
for (k = 0; k < nz; k++) {
for (j = 0; j < ny; j++) {
for (i = 0; i < nx; i++) {
for (m = 0; m < natoms; m++) {
tag_offset = k*ny*nx*natoms + j*nx*natoms + i*natoms;
if (!style_granular)
printf("%d %g %g %g\n",vtag[m]+tag_offset,vx[m],vy[m],vz[m]);
else
printf("%d %g %g %g %g %g %g\n",vtag[m]+tag_offset,
vx[m],vy[m],vz[m],vphix[m],vphiy[m],vphiz[m]);
}
}
}
}
}
/* read bonds and replicate */
else if (!strcmp(token,"Bonds")) {
for (i = 0; i < nbonds; i++) {
fgets(line,MAXLINE,stdin);
sscanf(line,"%d %d %d %d",&n,&btype[i],&bond1[i],&bond2[i]);
}
n = 0;
for (k = 0; k < nz; k++)
for (j = 0; j < ny; j++)
for (i = 0; i < nx; i++)
for (m = 0; m < nbonds; m++) {
n++;
del = k*ny*nx*natoms + j*nx*natoms + i*natoms;
printf("%d %d %d %d\n",n,btype[m],bond1[m]+del,bond2[m]+del);
}
}
/* read angles and replicate */
else if (!strcmp(token,"Angles")) {
for (i = 0; i < nangles; i++) {
fgets(line,MAXLINE,stdin);
sscanf(line,"%d %d %d %d %d",&n,&atype[i],
&angle1[i],&angle2[i],&angle3[i]);
}
n = 0;
for (k = 0; k < nz; k++)
for (j = 0; j < ny; j++)
for (i = 0; i < nx; i++)
for (m = 0; m < nangles; m++) {
n++;
del = k*ny*nx*natoms + j*nx*natoms + i*natoms;
printf("%d %d %d %d %d\n",n,atype[m],
angle1[m]+del,angle2[m]+del,angle3[m]+del);
}
}
/* read dihedrals and replicate */
else if (!strcmp(token,"Dihedrals")) {
for (i = 0; i < ndihedrals; i++) {
fgets(line,MAXLINE,stdin);
sscanf(line,"%d %d %d %d %d %d",&n,&dtype[i],&dihed1[i],&dihed2[i],
&dihed3[i],&dihed4[i]);
}
n = 0;
for (k = 0; k < nz; k++)
for (j = 0; j < ny; j++)
for (i = 0; i < nx; i++)
for (m = 0; m < ndihedrals; m++) {
n++;
del = k*ny*nx*natoms + j*nx*natoms + i*natoms;
printf("%d %d %d %d %d %d\n",n,dtype[m],
dihed1[m]+del,dihed2[m]+del,dihed3[m]+del,dihed4[m]+del);
}
}
/* read impropers and replicate */
else if (!strcmp(token,"Impropers")) {
for (i = 0; i < nimpropers; i++) {
fgets(line,MAXLINE,stdin);
sscanf(line,"%d %d %d %d %d %d",&n,&itype[i],&impro1[i],&impro2[i],
&impro3[i],&impro4[i]);
}
n = 0;
for (k = 0; k < nz; k++)
for (j = 0; j < ny; j++)
for (i = 0; i < nx; i++)
for (m = 0; m < nimpropers; m++) {
n++;
del = k*ny*nx*natoms + j*nx*natoms + i*natoms;
printf("%d %d %d %d %d %d\n",n,itype[m],
impro1[m]+del,impro2[m]+del,impro3[m]+del,impro4[m]+del);
}
}
/* non-replicated sections
just copy proper number of lines from in to out */
else if (!strcmp(token,"Masses"))
copyline(ntypes);
else if (!strcmp(token,"Dipoles"))
copyline(ntypes);
else if (!strcmp(token,"Pair Coeffs"))
copyline(ntypes);
else if (!strcmp(token,"Bond Coeffs"))
copyline(nbondtypes);
else if (!strcmp(token,"Angle Coeffs"))
copyline(nangletypes);
else if (!strcmp(token,"Dihedral Coeffs"))
copyline(ndihedtypes);
else if (!strcmp(token,"Improper Coeffs"))
copyline(nimprotypes);
else if (!strcmp(token,"BondBond Coeffs"))
copyline(nangletypes);
else if (!strcmp(token,"BondAngle Coeffs"))
copyline(nangletypes);
else if (!strcmp(token,"MiddleBondTorsion Coeffs"))
copyline(ndihedtypes);
else if (!strcmp(token,"EndBondTorsion Coeffs"))
copyline(ndihedtypes);
else if (!strcmp(token,"AngleTorsion Coeffs"))
copyline(ndihedtypes);
else if (!strcmp(token,"AngleAngleTorsion Coeffs"))
copyline(ndihedtypes);
else if (!strcmp(token,"BondBond13 Coeffs"))
copyline(ndihedtypes);
else if (!strcmp(token,"AngleAngle Coeffs"))
copyline(nimprotypes);
else {
fprintf(stderr,
"Error in input data file - unknown identifier %s\n",token);
exit(1);
}
token = gettoken(line,0);
}
}
/* ------------------------------------------------------------------- */
/* return a LAMMPS keyword from data file
if first is 1, non-blank line with token is in line
else read until find a non-blank line
keyword is all text on line w/out leading & trailing white space
read one additional line after non-blank line (assumed blank)
return ptr to keyword
return NULL if end of file */
char *gettoken(char *line, int first)
{
char buffer[MAXLINE];
// read upto non-blank line plus 1 following line
// eof is set to 1 if any read hits end-of-file
int eof = 0;
if (!first) if (fgets(line,MAXLINE,stdin) == NULL) eof = 1;
while (eof == 0 && strspn(line," \t\n\r") == strlen(line))
if (fgets(line,MAXLINE,stdin) == NULL) eof = 1;
if (fgets(buffer,MAXLINE,stdin) == NULL) eof = 1;
// if eof, return NULL
if (eof) return NULL;
// bracket non-whitespace portion of line
int start = strspn(line," \t\n\r");
int stop = strlen(line) - 1;
while (line[stop] == ' ' || line[stop] == '\t'
|| line[stop] == '\n' || line[stop] == '\r') stop--;
line[stop+1] = '\0';
// print keyword into output file
if (first) printf("%s\n\n",&line[start]);
else printf("\n%s\n\n",&line[start]);
// return ptr to keyword
return &line[start];
}
/* count and return words in a single line
make copy of line before using strtok so as not to change line
trim anything from '#' onward */
int count_words(char *line)
{
int n = strlen(line) + 1;
char *copy = (char *) malloc(n*sizeof(char));
strcpy(copy,line);
char *ptr;
if (ptr = strchr(copy,'#')) *ptr = '\0';
if (strtok(copy," \t\n") == NULL) return 0;
n = 1;
while (strtok(NULL," \t\n")) n++;
free(copy);
return n;
}
/* copy n lines from stdin to stdout */
void copyline(int n)
{
char line[1000];
while (n) {
fgets(line,MAXLINE,stdin);
printf("%s",line);
n--;
}
}
/* point (x,y,z) is somewhere in (or near) the replicated box "in"
rescale (x,y,z) so it is in the same relative location in box "out"
do this by 1st converting to 0-1 coord system in "in" space */
void scale(struct box in, int nx, int ny, int nz,
struct box out, double *x, double *y, double *z)
{
double rel;
rel = (*x-in.xlo) / (nx*in.xsize);
*x = out.xlo + rel*out.xsize;
rel = (*y-in.ylo) / (ny*in.ysize);
*y = out.ylo + rel*out.ysize;
rel = (*z-in.zlo) / (nz*in.zsize);
*z = out.zlo + rel*out.zsize;
}
/* remap the point (x,y,z) so it is guaranteed to be inside the box
using standard periodic imaging, update image flags */
void pbc(struct box box, double *x, double *y, double *z,
int *imagex, int *imagey, int *imagez)
{
while (*x < box.xlo) {
*x += box.xsize;
*imagex -= 1;
}
while (*x >= box.xhi) {
*x -= box.xsize;
*imagex += 1;
}
while (*y < box.ylo) {
*y += box.ysize;
*imagey -= 1;
}
while (*y >= box.yhi) {
*y -= box.ysize;
*imagey += 1;
}
while (*z < box.zlo) {
*z += box.zsize;
*imagez -= 1;
}
while (*z >= box.zhi) {
*z -= box.zsize;
*imagez += 1;
}
}
c4science · Help