Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F74204422
ReadCarFile.c
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Fri, Jul 26, 11:25
Size
9 KB
Mime Type
text/x-c
Expires
Sun, Jul 28, 11:25 (1 d, 22 h)
Engine
blob
Format
Raw Data
Handle
19354406
Attached To
rLAMMPS lammps
ReadCarFile.c
View Options
/*
* This function opens the .car file and extracts coordinate information
* into the atoms Atom structure
*/
#include "msi2lmp.h"
#include <stdlib.h>
#include <string.h>
#include <math.h>
/* ----------------------------------------------------------------------
set global box params
assumes boxlo/hi and triclinic tilts are already set
------------------------------------------------------------------------- */
void set_box(double box[3][3], double *h, double *h_inv)
{
h[0] = box[1][0] - box[0][0];
h[1] = box[1][1] - box[0][1];
h[2] = box[1][2] - box[0][2];
h_inv[0] = 1.0/h[0];
h_inv[1] = 1.0/h[1];
h_inv[2] = 1.0/h[2];
h[3] = box[2][2];
h[4] = box[2][1];
h[5] = box[2][0];
h_inv[3] = -h[3] / (h[1]*h[2]);
h_inv[4] = (h[3]*h[5] - h[1]*h[4]) / (h[0]*h[1]*h[2]);
h_inv[5] = -h[5] / (h[0]*h[1]);
}
/* ----------------------------------------------------------------------
convert triclinic 0-1 lamda coords to box coords for one atom
x = H lamda + x0;
lamda and x can point to same 3-vector
------------------------------------------------------------------------- */
void lamda2x(double *lamda, double *x, double *h, double *boxlo)
{
x[0] = h[0]*lamda[0] + h[5]*lamda[1] + h[4]*lamda[2] + boxlo[0];
x[1] = h[1]*lamda[1] + h[3]*lamda[2] + boxlo[1];
x[2] = h[2]*lamda[2] + boxlo[2];
}
/* ----------------------------------------------------------------------
convert box coords to triclinic 0-1 lamda coords for one atom
lamda = H^-1 (x - x0)
x and lamda can point to same 3-vector
------------------------------------------------------------------------- */
void x2lamda(double *x, double *lamda, double *h_inv, double *boxlo)
{
double delta[3];
delta[0] = x[0] - boxlo[0];
delta[1] = x[1] - boxlo[1];
delta[2] = x[2] - boxlo[2];
lamda[0] = h_inv[0]*delta[0] + h_inv[5]*delta[1] + h_inv[4]*delta[2];
lamda[1] = h_inv[1]*delta[1] + h_inv[3]*delta[2];
lamda[2] = h_inv[2]*delta[2];
}
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 */
double total_q;
double sq_c;
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];
double hmat[6];
double hinv[6];
double lamda[3];
/* Open .car file for reading */
sprintf(line,"%s.car",rootname);
if (pflag > 0) printf(" Reading car file: %s\n",line);
if( (CarF = fopen(line,"r")) == NULL ) {
printf("Cannot open %s\n",line);
exit(33);
}
/* Determine Number of molecules & atoms */
rewind(CarF);
no_molecules = -1; /* Set to -1 because counter will be incremented an
extra time at the end of the file */
fgets(line,MAX_LINE_LENGTH,CarF); /* Read header line */
/* Check for periodicity, if present, read cell constants */
if( strncmp(fgets(line,MAX_LINE_LENGTH,CarF),"PBC=ON",6) == 0) {
periodic = 1;
skip = 5; /* Data starts 5 lines from beginning of file */
fgets(line,MAX_LINE_LENGTH,CarF); /* Comment line */
fgets(line,MAX_LINE_LENGTH,CarF); /* Date stamp */
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 */
if(pbc[3] != 90.0 || pbc[4] != 90.0 || pbc[5] != 90.0) {
TriclinicFlag = 1;
} else TriclinicFlag = 0;
} else {
periodic = 0;
skip = 4;
if (pflag > 1) {
printf(" %s is not a periodic system\n", rootname);
printf(" Assigning cell parameters based on coordinates\n");
}
fgets(line,MAX_LINE_LENGTH, CarF); /* Comment line */
fgets(line,MAX_LINE_LENGTH, CarF); /* Date Stamp */
}
/* First pass through file -- Count molecules */
while(fgets(line,MAX_LINE_LENGTH,CarF) != NULL )
if( strncmp(line,"end",3) == 0 )
no_molecules++;
/* Allocate space to keep track of the number of atoms within a molecule */
no_atoms = (int *) calloc(no_molecules,sizeof(int));
if ( no_atoms == NULL ) {
printf("Could not allocate memory for no_atoms\n");
exit(32);
}
/* Second pass through file -- Count atoms */
rewind(CarF);
for(n=0; n < skip; n++) /* Skip beginning lines */
fgets(line,MAX_LINE_LENGTH,CarF);
for(n=0; n < no_molecules; n++)
while( strncmp(fgets(line,MAX_LINE_LENGTH,CarF),"end",3) )
no_atoms[n]++;
for( total_no_atoms=0, n=0; n < no_molecules; n++ )
total_no_atoms += no_atoms[n];
molecule = (struct MoleculeList *) calloc(no_molecules,
sizeof(struct MoleculeList));
if (molecule == NULL) {
printf("Unable to allocate memory for molecule structure\n");
exit(32);
}
molecule[0].start = 0;
molecule[0].end = no_atoms[0];
for (n=1; n < no_molecules; n++) {
molecule[n].start = molecule[n-1].end;
molecule[n].end = molecule[n].start + no_atoms[n];
}
/* Allocate space for atoms Atom structures */
atoms = (struct Atom *) calloc(total_no_atoms,sizeof(struct Atom));
if( atoms == NULL ) {
printf("Could not allocate memory for AtomList\n");
exit(32);
}
/* 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);
for(m=0; m < no_molecules; m++) {
for(k=molecule[m].start; k <
molecule[m].end; k++) {
atoms[k].molecule = m;
atoms[k].no = k;
fscanf(CarF,"%s %lf %lf %lf %*s %d %s %s %lf",
atoms[k].name,
&(atoms[k].x[0]),
&(atoms[k].x[1]),
&(atoms[k].x[2]),
&(atoms[k].molecule),
atoms[k].potential,
atoms[k].element,
&(atoms[k].q));
atoms[k].x[0] += shift[0];
atoms[k].x[1] += shift[1];
atoms[k].x[2] += shift[2];
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;
if (pflag > 1) {
printf(" There are %d atoms in %d molecules in this file\n",
total_no_atoms,no_molecules);
printf(" The total charge in the system is %7.3f.\n",total_q);
}
/* Search coordinates to find lowest and highest for x, y, and z */
if (periodic == 0) {
/* 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];
for ( m = 1; m < total_no_atoms; m++) {
if (atoms[m].x[k] < lowest) lowest = atoms[m].x[k];
if (atoms[m].x[k] > highest) highest = atoms[m].x[k];
}
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");
exit(32);
}
} else {
if (pflag > 2)
printf(" pbc[0] %f pbc[1] %f pbc[2] %f\n", pbc[0] ,pbc[1] ,pbc[2]);
if (TriclinicFlag == 0) {
for (k=0; k < 3; k++) {
box[0][k] = -0.5*pbc[k] + center[k] + shift[k];
box[1][k] = 0.5*pbc[k] + center[k] + shift[k];
box[2][k] = 0.0;
}
} else {
sq_c = pbc[2]*pbc[2];
cos_alpha = cos(pbc[3]*PI_180);
cos_gamma = cos(pbc[5]*PI_180);
sin_gamma = sin(pbc[5]*PI_180);
cos_beta = cos(pbc[4]*PI_180);
sin_beta = sin(pbc[4]*PI_180);
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];
box[0][0] = -0.5*A + center[0] + shift[0];
box[1][0] = 0.5*A + center[0] + shift[0];
box[0][1] = -0.5*B*sin_gamma + center[1] + shift[1];
box[1][1] = 0.5*B*sin_gamma + center[1] + shift[1];
box[0][2] = -0.5*sqrt(sq_c * sin_beta*sin_beta - C*(cos_alpha-cos_gamma*cos_beta)/sin_gamma) + center[2] + shift[2];
box[1][2] = 0.5*sqrt(sq_c * sin_beta*sin_beta - C*(cos_alpha-cos_gamma*cos_beta)/sin_gamma) + center[2] + shift[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 */
}
}
/* compute image flags */
set_box(box,hmat,hinv);
n = 0;
for (m = 0; m < total_no_atoms; m++) {
double tmp;
int w=0;
x2lamda(atoms[m].x,lamda,hinv,box[0]);
for (k = 0; k < 3; ++k) {
tmp = floor(lamda[k]);
atoms[m].image[k] = tmp;
lamda[k] -= tmp;
if (tmp != 0.0) ++w;
}
lamda2x(lamda, atoms[m].x,hmat,box[0]);
if (w > 0) ++n;
}
/* warn if atoms are outside the box */
if (n > 0) {
if (periodic) {
if (pflag > 1)
printf(" %d of %d atoms with nonzero image flags\n\n",n,total_no_atoms);
} else {
if (iflag == 0 || (pflag > 1))
printf(" %d of %d atoms outside the box\n\n",n,total_no_atoms);
condexit(32);
}
}
/* Close .car file */
if (fclose(CarF) !=0) {
printf("Error closing %s.car\n", rootname);
exit(31);
}
}
/* End ReadCarFile() */
Event Timeline
Log In to Comment