diff --git a/tools/msi2lmp/src/Forcefield.h b/tools/msi2lmp/src/Forcefield.h index d6ce024f8e..085e6c2c6b 100644 --- a/tools/msi2lmp/src/Forcefield.h +++ b/tools/msi2lmp/src/Forcefield.h @@ -6,7 +6,7 @@ * It defines the data structures used to store the force field in memory */ -#define MAX_NO_MEMS 5 +#define MAX_NO_MEMS 6 #define MAX_NO_PARAMS 8 struct FrcFieldData { diff --git a/tools/msi2lmp/src/GetParameters.c b/tools/msi2lmp/src/GetParameters.c index a0729ea22c..22e3524507 100644 --- a/tools/msi2lmp/src/GetParameters.c +++ b/tools/msi2lmp/src/GetParameters.c @@ -73,7 +73,7 @@ void GetParameters(int Forcefield) if (k < 0) { get_equivs(1,potential_types,equiv_types); - if (pflag > 2) printf("Using equivalences for VDW %s -> %s\n", + if (pflag > 2) printf(" Using equivalences for VDW %s -> %s\n", potential_types[0],equiv_types[0]); k = find_match(1,equiv_types,ff_vdw,&backwards); @@ -124,7 +124,7 @@ void GetParameters(int Forcefield) get_equivs(2,potential_types,equiv_types); if (pflag > 2) { - printf("Using equivalences for bond %s %s -> %s %s\n", + printf(" Using equivalences for bond %s %s -> %s %s\n", potential_types[0],potential_types[1], equiv_types[0],equiv_types[1]); } @@ -149,7 +149,7 @@ void GetParameters(int Forcefield) printf("\n Bond Types and Parameters\n"); for (i=0; i < no_bond_types; i++) { for (j=0; j < 2; j++) - printf("%-3s",atomtypes[bondtypes[i].types[j]].potential); + printf(" %-3s",atomtypes[bondtypes[i].types[j]].potential); for (j=0; j < 4; j++) printf(" %8.4f",bondtypes[i].params[j]); printf("\n"); @@ -177,7 +177,7 @@ void GetParameters(int Forcefield) if (k < 0) { get_equivs(3,potential_types,equiv_types); if (pflag > 2) { - printf("Using equivalences for angle %s %s %s -> %s %s %s\n", + printf(" Using equivalences for angle %s %s %s -> %s %s %s\n", potential_types[0],potential_types[1], potential_types[2], equiv_types[0],equiv_types[1], @@ -201,7 +201,7 @@ void GetParameters(int Forcefield) if (Forcefield > 1) { get_equivs(3,potential_types,equiv_types); if (pflag > 2) { - printf("Using equivalences for 3 body cross terms %s %s %s -> %s %s %s\n", + printf(" Using equivalences for 3 body cross terms %s %s %s -> %s %s %s\n", potential_types[0],potential_types[1],potential_types[2], equiv_types[0],equiv_types[1],equiv_types[2]); } @@ -260,7 +260,7 @@ void GetParameters(int Forcefield) printf("\n BondBond Types and Parameters\n"); for (i=0; i < no_angle_types; i++) { for (j=0; j < 3; j++) - printf("%-3s",atomtypes[angletypes[i].types[j]].potential); + printf(" %-3s",atomtypes[angletypes[i].types[j]].potential); for (j=0; j < 3; j++) printf(" %8.4f",angletypes[i].bondbond_cross_term[j]); printf("\n"); @@ -268,7 +268,7 @@ void GetParameters(int Forcefield) printf("\n BondAngle Types and Parameters\n"); for (i=0; i < no_angle_types; i++) { for (j=0; j < 3; j++) - printf("%-3s",atomtypes[angletypes[i].types[j]].potential); + printf(" %-3s",atomtypes[angletypes[i].types[j]].potential); for (j=0; j < 4; j++) printf(" %8.4f",angletypes[i].bondangle_cross_term[j]); printf("\n"); @@ -302,7 +302,7 @@ void GetParameters(int Forcefield) get_equivs(4,potential_types,equiv_types); if (pflag > 2) { - printf("Using equivalences for dihedral %s %s %s %s -> %s %s %s %s\n", + printf(" Using equivalences for dihedral %s %s %s %s -> %s %s %s %s\n", potential_types[0],potential_types[1], potential_types[2],potential_types[3], equiv_types[0],equiv_types[1], @@ -333,7 +333,7 @@ void GetParameters(int Forcefield) else if (ff_tor.data[k].ff_param[2] == 180.0) dihedraltypes[i].params[1] = -1.0; else { - printf("Non planar phi0 for %s %s %s %s\n", + printf(" Non planar phi0 for %s %s %s %s\n", potential_types[0],potential_types[1], potential_types[2],potential_types[3]); dihedraltypes[i].params[1] = 0.0; @@ -349,7 +349,7 @@ void GetParameters(int Forcefield) if (Forcefield > 1) { get_equivs(4,potential_types,equiv_types); if (pflag > 2) { - printf("Using equivalences for linear 4 body cross terms %s %s %s %s -> %s %s %s %s\n", + printf(" Using equivalences for linear 4 body cross terms %s %s %s %s -> %s %s %s %s\n", potential_types[0],potential_types[1], potential_types[2],potential_types[3], equiv_types[0],equiv_types[1], @@ -529,7 +529,7 @@ void GetParameters(int Forcefield) printf("\n Dihedral Types and Parameters\n"); for (i=0; i < no_dihedral_types; i++) { for (j=0; j < 4; j++) - printf("%-3s",atomtypes[dihedraltypes[i].types[j]].potential); + printf(" %-3s",atomtypes[dihedraltypes[i].types[j]].potential); for (j=0; j < 6; j++) printf(" %8.4f",dihedraltypes[i].params[j]); printf("\n"); @@ -540,7 +540,7 @@ void GetParameters(int Forcefield) printf("\n EndBondDihedral Types and Parameters\n"); for (i=0; i < no_dihedral_types; i++) { for (j=0; j < 4; j++) - printf("%-3s",atomtypes[dihedraltypes[i].types[j]].potential); + printf(" %-3s",atomtypes[dihedraltypes[i].types[j]].potential); for (j=0; j < 8; j++) printf(" %8.4f",dihedraltypes[i].endbonddihedral_cross_term[j]); printf("\n"); @@ -557,7 +557,7 @@ void GetParameters(int Forcefield) printf("\n AngleDihedral Types and Parameters\n"); for (i=0; i < no_dihedral_types; i++) { for (j=0; j < 4; j++) - printf("%-3s",atomtypes[dihedraltypes[i].types[j]].potential); + printf(" %-3s",atomtypes[dihedraltypes[i].types[j]].potential); for (j=0; j < 8; j++) printf(" %8.4f",dihedraltypes[i].angledihedral_cross_term[j]); printf("\n"); @@ -568,7 +568,7 @@ void GetParameters(int Forcefield) for (j=0; j < 4; j++) printf(" %-3s",atomtypes[dihedraltypes[i].types[j]].potential); for (j=0; j < 3; j++) - printf("%8.4f",dihedraltypes[i].angleangledihedral_cross_term[j]); + printf(" %8.4f",dihedraltypes[i].angleangledihedral_cross_term[j]); printf("\n"); } @@ -618,7 +618,7 @@ void GetParameters(int Forcefield) get_equivs(5,potential_types,equiv_types); if (pflag > 2) { - printf("Using equivalences for oop %s %s %s %s -> %s %s %s %s\n", + printf(" Using equivalences for oop %s %s %s %s -> %s %s %s %s\n", potential_types[0],potential_types[1], potential_types[2],potential_types[3], equiv_types[0],equiv_types[1], @@ -638,7 +638,7 @@ void GetParameters(int Forcefield) else if (ff_oop.data[k].ff_param[2] == 180.0) ooptypes[i].params[1] = -1.0; else { - printf("Non planar phi0 for %s %s %s %s\n", + printf(" Non planar phi0 for %s %s %s %s\n", potential_types[0],potential_types[1], potential_types[2],potential_types[3]); ooptypes[i].params[1] = 0.0; @@ -658,7 +658,7 @@ void GetParameters(int Forcefield) if (k < 0) { get_equivs(5,potential_types,equiv_types); if (pflag > 2) { - printf("Using equivalences for oop %s %s %s %s -> %s %s %s %s\n", + printf(" Using equivalences for oop %s %s %s %s -> %s %s %s %s\n", potential_types[0],potential_types[1], potential_types[2],potential_types[3], equiv_types[0],equiv_types[1], @@ -682,7 +682,7 @@ void GetParameters(int Forcefield) printf("\n OOP Types and Parameters\n"); for (i=0; i < no_oop_types; i++) { for (j=0; j < 4; j++) - printf("%-3s",atomtypes[ooptypes[i].types[j]].potential); + printf(" %-3s",atomtypes[ooptypes[i].types[j]].potential); for (j=0; j < 3; j++) printf(" %8.4f",ooptypes[i].params[j]); printf("\n"); @@ -733,7 +733,7 @@ void GetParameters(int Forcefield) if (k < 0) { get_equivs(5,potential_types,equiv_types); if (pflag > 2) { - printf("Using equivalences for angleangle %s %s %s %s -> %s %s %s %s\n", + printf(" Using equivalences for angleangle %s %s %s %s -> %s %s %s %s\n", potential_types[0],potential_types[1], potential_types[2],potential_types[3], equiv_types[0],equiv_types[1], @@ -803,7 +803,7 @@ void GetParameters(int Forcefield) printf("\n AngleAngle Types and Parameters\n"); for (i=0; i < no_oop_types; i++) { for (j=0; j < 4; j++) - printf("%-3s",atomtypes[ooptypes[i].types[j]].potential); + printf(" %-3s",atomtypes[ooptypes[i].types[j]].potential); for (j=0; j < 6; j++) printf(" %8.4f",ooptypes[i].angleangle_params[j]); printf("\n"); @@ -1193,7 +1193,7 @@ double get_r0(int typei,int typej) } if (match == 0) - printf("Unable to find r0 for types %d %d\n",typei,typej); + printf(" Unable to find r0 for types %d %d\n",typei,typej); return r; } @@ -1279,13 +1279,12 @@ void get_equivs(int ic,char potential_types[][5],char equiv_types[][5]) for (i=0; i < 4; i++) { k = find_equiv_type(potential_types[i]); if (k > -1) - /* XXX: this leads to an out-of-bounds access. - the change below does not seem to make a difference. - strncpy(equiv_types[i],equivalence.data[k].ff_types[5],5); */ - strncpy(equiv_types[i],equivalence.data[k].ff_types[4],5); + strncpy(equiv_types[i],equivalence.data[k].ff_types[5],5); } break; default: + printf(" Requesting equivalences of unsupported type: %d\n",ic); + condexit(26); break; } return; diff --git a/tools/msi2lmp/src/Makefile b/tools/msi2lmp/src/Makefile index bd12849f35..26ced8c31a 100644 --- a/tools/msi2lmp/src/Makefile +++ b/tools/msi2lmp/src/Makefile @@ -10,8 +10,7 @@ SRCS = msi2lmp.c \ SearchAndFill.c \ GetParameters.c \ CheckLists.c \ - WriteDataFile05.c \ - WriteDataFile01.c + WriteDataFile.c OBJS = $(SRCS:.c=.o) @@ -42,5 +41,4 @@ ReadCarFile.o: ReadCarFile.c msi2lmp.h ReadFrcFile.o: ReadFrcFile.c msi2lmp.h Forcefield.h ReadMdfFile.o: ReadMdfFile.c msi2lmp.h SearchAndFill.o: SearchAndFill.c msi2lmp.h Forcefield.h -WriteDataFile01.o: WriteDataFile01.c msi2lmp.h Forcefield.h -WriteDataFile05.o: WriteDataFile05.c msi2lmp.h Forcefield.h +WriteDataFile.o: WriteDataFile.c msi2lmp.h Forcefield.h diff --git a/tools/msi2lmp/src/ReadCarFile.c b/tools/msi2lmp/src/ReadCarFile.c index 5edb407d47..44bdabda69 100644 --- a/tools/msi2lmp/src/ReadCarFile.c +++ b/tools/msi2lmp/src/ReadCarFile.c @@ -12,17 +12,18 @@ void ReadCarFile(void) { char line[MAX_LINE_LENGTH]; /* Stores lines as they are read in */ - int k,m,n; /* counters */ - int skip; /* lines to skip at beginning of file */ - double lowest, highest; /* temp coordinate finding variables */ + int k,m,n; /* counters */ + int skip; /* lines to skip at beginning of file */ + double lowest, highest; /* temp coordinate finding variables */ double total_q; double sq_c; - double cos_alpha; // Added by SLTM Sept 13, 2010 + double cos_alpha; /* Added by SLTM Sept 13, 2010 */ double cos_gamma; double sin_gamma; double cos_beta; double sin_beta; double A, B, C; + double center[3]; /* Open .car file for reading */ @@ -51,7 +52,7 @@ void ReadCarFile(void) fscanf(CarF,"%*s %lf %lf %lf %lf %lf %lf %*s", &pbc[0],&pbc[1],&pbc[2],&pbc[3],&pbc[4],&pbc[5]); - // Added triclinic flag for non-orthogonal boxes Oct 5, 2010 SLTM + /* Added triclinic flag for non-orthogonal boxes Oct 5, 2010 SLTM */ if(pbc[3] != 90.0 || pbc[4] != 90.0 || pbc[5] != 90.0) { TriclinicFlag = 1; } else TriclinicFlag = 0; @@ -115,7 +116,7 @@ void ReadCarFile(void) } /* Third pass through file -- Read+Parse Car File */ - + center[0] = center[1] = center[2] = 0.0; rewind(CarF); for(n=0; n < skip; n++) fgets(line,MAX_LINE_LENGTH,CarF); @@ -136,12 +137,21 @@ void ReadCarFile(void) atoms[k].potential, atoms[k].element, &(atoms[k].q)); + if (centerflag) { + center[0] += atoms[k].x[0]; + center[1] += atoms[k].x[1]; + center[2] += atoms[k].x[2]; + } } fgets(line,MAX_LINE_LENGTH,CarF); fgets(line,MAX_LINE_LENGTH,CarF); } /* End m (molecule) loop */ + center[0] /= (double) total_no_atoms; + center[1] /= (double) total_no_atoms; + center[2] /= (double) total_no_atoms; + for (total_q=0.0,k=0; k < total_no_atoms; k++) total_q += atoms[k].q; @@ -154,8 +164,10 @@ void ReadCarFile(void) /* Search coordinates to find lowest and highest for x, y, and z */ if (periodic == 0) { - // Added if/else statment STLM Oct 5 2010 + /* Added if/else statment STLM Oct 5 2010 */ if (TriclinicFlag == 0) { + /* no need to re-center the box, if we use min/max values */ + center[0] = center[1] = center[2] = 0.0; for ( k = 0; k < 3; k++) { lowest = atoms[0].x[k]; highest = atoms[0].x[k]; @@ -164,8 +176,9 @@ void ReadCarFile(void) if (atoms[m].x[k] < lowest) lowest = atoms[m].x[k]; if (atoms[m].x[k] > highest) highest = atoms[m].x[k]; } - pbc[k] = lowest; - pbc[k+3] = highest; + box[0][k] = lowest - 0.5; + box[1][k] = highest + 0.5; + box[2][k] = 0.0; } } else { printf("This tool only works for periodic systems with triclinic boxes"); @@ -173,11 +186,12 @@ void ReadCarFile(void) } } else { - // Modified lines 176 - 201 Oct 5th 2010 + /* Modified lines 176 - 201 Oct 5th 2010 */ if (TriclinicFlag == 0) { for (k=0; k < 3; k++) { - pbc[k+3] = pbc[k]; - pbc[k] = 0.0; + box[0][k] = -0.5*pbc[k] + center[k]; + box[1][k] = 0.5*pbc[k] + center[k]; + box[2][k] = 0.0; } } else { sq_c = pbc[2]*pbc[2]; @@ -186,21 +200,24 @@ void ReadCarFile(void) sin_gamma = sin(pbc[5]*PI_180); cos_beta = cos(pbc[4]*PI_180); sin_beta = sin(pbc[4]*PI_180); - if (pflag > 1) { - printf("pbc[3] %f pbc[4] %f pbc[5] %f\n", pbc[3] ,pbc[4] ,pbc[5]); - printf("cos_alpha %f cos_beta %f cos_gamma %f\n", cos_alpha ,cos_beta ,cos_gamma); + if (pflag > 2) { + printf(" pbc[3] %f pbc[4] %f pbc[5] %f\n", pbc[3] ,pbc[4] ,pbc[5]); + printf(" cos_alpha %f cos_beta %f cos_gamma %f\n", cos_alpha ,cos_beta ,cos_gamma); } A = pbc[0]; B = pbc[1]; C = pbc[2]; - pbc[0] = A; - pbc[1] = B*sin_gamma; - pbc[2] = sqrt(sq_c * sin_beta*sin_beta - C*(cos_alpha-cos_gamma*cos_beta)/sin_gamma); - pbc[3] = B * cos_gamma; // This is xy SLTM - pbc[4] = C * cos_beta; // This is xz SLTM - pbc[5] = C*(cos_alpha-cos_gamma*cos_beta)/sin_gamma; // This is yz SLTM + box[0][0] = -0.5*A + center[0]; + box[1][0] = 0.5*A + center[0]; + box[0][1] = -0.5*B*sin_gamma + center[1]; + box[1][1] = 0.5*B*sin_gamma + center[1]; + box[0][2] = -0.5*sqrt(sq_c * sin_beta*sin_beta - C*(cos_alpha-cos_gamma*cos_beta)/sin_gamma) + center[2]; + box[1][2] = 0.5*sqrt(sq_c * sin_beta*sin_beta - C*(cos_alpha-cos_gamma*cos_beta)/sin_gamma) + center[2]; + box[2][0] = B * cos_gamma; /* This is xy SLTM */ + box[2][1] = C * cos_beta; /* This is xz SLTM */ + box[2][2] = C*(cos_alpha-cos_gamma*cos_beta)/sin_gamma; /* This is yz SLTM */ } } diff --git a/tools/msi2lmp/src/msi2lmp.c b/tools/msi2lmp/src/msi2lmp.c index 6e26ad8791..b07d883894 100644 --- a/tools/msi2lmp/src/msi2lmp.c +++ b/tools/msi2lmp/src/msi2lmp.c @@ -1,6 +1,6 @@ /* * -* msi2lmp.exe V3.6 +* msi2lmp.exe V3.8 * * v3.6 KLA - Changes to output to either lammps 2001 (F90 version) or to * lammps 2005 (C++ version) @@ -33,10 +33,11 @@ * The program is started by supplying information at the command prompt * according to the usage described below. * -* USAGE: msi2lmp3 ROOTNAME {-print #} {-class #} {-frc FRC_FILE} -2001 +* USAGE: msi2lmp3 ROOTNAME {-print #} {-class #} {-frc FRC_FILE} {-ignore} {-nocenter} * * -- msi2lmp3 is the name of the executable * -- ROOTNAME is the base name of the .car and .mdf files +* -- all opther flags are optional and can be abbreviated (e.g. -p instead of -print) * * -- -print * # is the print level: 0 - silent except for errors @@ -48,7 +49,12 @@ * (II = Class II e.g., CFFx ) * default is -class I * -* -- -frc - specifies name of the forcefield file (e.g., cff91) +* -- -ignore - tells msi2lmp to ignore warnings and errors and keep going +* +* -- -nocenter - tells msi2lmp to not center the box around the (geometrical) +* center of the atoms, but around the origin +* +* -- -frc - specifies name of the forcefield file (e.g., cff91) * * If the name includes a hard wired directory (i.e., if the name * starts with . or /), then the name is used alone. Otherwise, @@ -70,10 +76,7 @@ * * By default, the program uses $BIOSYM_LIBRARY/cvff.frc * -* -- -2001 will output a data file for the FORTRAN 90 version of LAMMPS (2001) -* By default, the program will output for the C++ version of LAMMPS. -* -* -- output is written to a file called ROOTNAME.lammps{01/05} +* -- output is written to a file called ROOTNAME.data * * **************************************************************** @@ -116,10 +119,12 @@ /* global variables */ char *rootname; -double pbc[9]; +double pbc[6]; +double box[3][3]; int periodic = 1; int TriclinicFlag = 0; int forcefield = 0; +int centerflag = 1; int pflag; int iflag; @@ -175,12 +180,10 @@ static int check_arg(char **arg, const char *flag, int num, int argc) int main (int argc, char *argv[]) { - int n,i,found_sep; /* Counter */ - int outv; + int n,i,found_sep; const char *frc_dir_name = NULL; const char *frc_file_name = NULL; - outv = 2005; pflag = 1; iflag = 0; forcefield = 1; @@ -188,7 +191,7 @@ int main (int argc, char *argv[]) frc_dir_name = getenv("BIOSYM_LIBRARY"); if (argc < 2) { - printf("usage: %s [-class ] [-frc ] [-p #] [-i]\n",argv[0]); + printf("usage: %s [-class ] [-frc ] [-print #] [-ignore] [-nocenter]\n",argv[0]); return 1; } else { /* rootname was supplied as first argument, copy to rootname */ int len = strlen(argv[1]) + 1; @@ -210,10 +213,6 @@ int main (int argc, char *argv[]) printf("Unrecognized Forcefield class: %s\n",argv[n]); return 3; } - } else if (strcmp(argv[n],"-2001") == 0) { - outv = 2001; - } else if (strcmp(argv[n],"-2005") == 0) { - outv = 2005; } else if (strncmp(argv[n],"-f",2) == 0) { n++; if (check_arg(argv,"-frc",n,argc)) @@ -221,6 +220,8 @@ int main (int argc, char *argv[]) frc_file_name = argv[n]; } else if (strncmp(argv[n],"-i",2) == 0 ) { iflag = 1; + } else if (strncmp(argv[n],"-n",2) == 0 ) { + centerflag = 0; } else if (strncmp(argv[n],"-p",2) == 0) { n++; if (check_arg(argv,"-print",n,argc)) @@ -333,12 +334,8 @@ int main (int argc, char *argv[]) printf("\n Check parameters for internal consistency\n"); CheckLists(); - if (outv == 2001) { - WriteDataFile01(rootname,forcefield); - - } else if (outv == 2005) { - WriteDataFile05(rootname,forcefield); - } + /* Write out the final data */ + WriteDataFile(rootname,forcefield); free(rootname); if (pflag > 0) diff --git a/tools/msi2lmp/src/msi2lmp.h b/tools/msi2lmp/src/msi2lmp.h index f8c83202fe..6339432bc2 100644 --- a/tools/msi2lmp/src/msi2lmp.h +++ b/tools/msi2lmp/src/msi2lmp.h @@ -24,12 +24,14 @@ * and to make the program fully dynamic. The second version used * fixed dimension arrays for the internal coordinates. * -* John Carpenter can be contacted by sending email to -* jec374@earthlink.net -* -* November 2000 -* +* The thrid version was revised in Fall 2011 by +* Stephanie Teich-McGoldrick to add support non-orthogonal cells. * +* The next revision was done in Summer 2013 by +* Axel Kohlmeyer to improve portability to Windows compilers, +* clean up command line parsing and improve compatibility with +* the then current LAMMPS versions. This revision removes +* compatibility with the obsolete LAMMPS version written in Fortran 90. */ # include @@ -153,10 +155,12 @@ struct Atom { extern char *rootname; extern char *FrcFileName; -extern double pbc[9]; +extern double pbc[6]; /* A, B, C, alpha, beta, gamma */ +extern double box[3][3]; /* hi/lo for x/y/z and xy, xz, yz for triclinic */ extern int periodic; /* 0= nonperiodic 1= 3-D periodic */ extern int TriclinicFlag; /* 0= Orthogonal 1= Triclinic */ extern int forcefield; /* 1= ClassI 2= ClassII */ +extern int centerflag; /* 1= center box 0= keep box */ extern int pflag; /* print level: 0, 1, 2, 3 */ extern int iflag; /* 0 stop at errors 1 = ignore errors */ extern int *no_atoms; @@ -200,5 +204,4 @@ extern void ReadFrcFile(); extern void MakeLists(); extern void GetParameters(int); extern void CheckLists(); -extern void WriteDataFile01(char *,int); -extern void WriteDataFile05(char *,int); +extern void WriteDataFile(char *,int);