git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@10275 f3b2605a-c512-4ea7-a41b-209d697bcdaa

This commit is contained in:
sjplimp 2013-07-22 22:49:21 +00:00
parent 9034be4e5e
commit 268c782771
6 changed files with 95 additions and 81 deletions

View File

@ -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 {

View File

@ -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;

View File

@ -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

View File

@ -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 */
}
}

View File

@ -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 <rootname> [-class <I|1|II|2>] [-frc <path to frc file>] [-p #] [-i]\n",argv[0]);
printf("usage: %s <rootname> [-class <I|1|II|2>] [-frc <path to frc file>] [-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)

View File

@ -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 <stdio.h>
@ -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);