tools/reax/mol_fra.c

   1 /*** modified on 12/09/2010 by Aidan Thompson               ***/
   2 /*** modified on 10/21/2010 by avv/cfdrc                    ***/
   3 /*** this code was given to LAMMPS by Sergey Zybin, CalTech ***/
   4
   5 #include <stdio.h>
   6 #include <stdlib.h>
   7 #include <math.h>
   8
   9 #define Cutoff "Cutoff.dic"
  10 #define Snapshots "bonds.reax"
  11 #define Output1 "fra.out"
  12 #define Output2 "mol.out"
  13 #define Output3 "num.out"
  14 #define Output4 "dipol.out"
  15 #define Output5 "vibtemp.out"
  16 #define Output6 "conc.out"
  17
  18 #define MaxNumBonds 10 //Max Number of bonds per atom
  19 #define MaxMolTypes 80000 //Max Number of molecule types
  20 #define Natomtypes 4 //C,H,O,N
  21
  22 #define Carbon 0
  23 #define Hydrogen 1
  24 #define Nitrogen 3
  25 #define Oxygen 2
  26
  27 // Prototypes for functions
  28
  29 void FileOpenError(FILE *fp,char *name);
  30 void Allocation();
  31 int AtomIndicator(char *symbol);
  32 void SkipLine(FILE *fp);
  33
  34 int Natom,Nmolecule;//total numbers of atoms and molecules
  35 int Nmoltype=0;//number of species
  36 int *Atomdic;//this array stores atom type accessible by atom id-1.
  37 int *MolName;//an array that stores number of atoms of each type for each molecule by using counters of Natomtypes*imoleculetypes+atomtype so if H2 is the second molecule entry the array will have 4=0, 5=2, 6=0, 7=0.
  38 int *MolType;//array with the same info as MolName but accessed by Nmoltype, not molecule id.
  39 int *Mol;//an array that stores the molecule id for each atom
  40 int *NMol;//an array that stores the number of molecules of this species for this frame
  41 double *Netc;//an array that stores charge
  42 double Mass[Natomtypes];
  43 double *rx,*ry,*rz,*Charge;
  44 double *vx,*vy,*vz;
  45 double COdic[Natomtypes][Natomtypes];
  46 double latc_x[3],latc_y[3],latc_z[3];//lattice parameters of the box read from config.out
  47 double Ilatc_x[3],Ilatc_y[3],Ilatc_z[3];//normalized box lengths
  48 double Origin[3];//origin coordinates of the box
  49 FILE *trj,*dpol,*vtemp,*conc;
  50
  51 main() {
  52
  53   void GetAtomlist();
  54   void GetAtomlist_mod();
  55   void GetCOdict();
  56   void AnalyzeTrajectory();
  57   void PostProcess();
  58
  59   GetAtomlist_mod(Snapshots);
  60   GetCOdict(Cutoff);
  61   AnalyzeTrajectory(Snapshots,Output3);
  62   PostProcess(Output1,Output2,Output3);
  63
  64 }
  65
  66 void GetAtomlist_mod(char *name) {
  67 /*This subroutine reads the bonds.petn file to set the number of atoms.  It also checks that a full frame exists because it loops over the first frame to populate the Atomdic with the proper type for each atom id.*/
  68   FILE *fp;
  69   char buffer[1000];
  70   int  i,j, Nlink, iatom, itemp, id;
  71   double BO,dtemp;//BO is bond order
  72
  73   fp=fopen(name,"r");
  74   FileOpenError(fp,name);
  75
  76  while(fscanf(fp,"%s",buffer)!=EOF) {
  77
  78    if(strcmp(buffer,"particles")==0) {
  79       fscanf(fp,"%d",&Natom);
  80       printf("# of atoms : %d\n",Natom);
  81       Allocation();
  82     }
  83     if(strcmp(buffer,"q")==0) {
  84
  85      for(i=0;i<Natom;i++) {
  86         fscanf(fp,"%d",&iatom);
  87         iatom--;
  88         fscanf(fp,"%d",&itemp);
  89         itemp--;
  90         if (itemp == Carbon){
  91            sprintf(buffer, "C");
  92         }else if(itemp == Hydrogen) {
  93           sprintf(buffer, "H");
  94         }else if(itemp == Oxygen) {
  95           sprintf(buffer, "O");
  96         }else if(itemp == Nitrogen) {
  97           sprintf(buffer, "N");
  98         }
  99
 100         fscanf(fp,"%d",&Nlink);
 101         for(j=0;j<Nlink;j++) {
 102           fscanf(fp,"%d",&itemp);
 103           itemp--;
 104         }
 105         fscanf(fp,"%d",&itemp);
 106
 107         for(j=0;j<Nlink;j++) {
 108           fscanf(fp,"%lf",&BO);
 109         }
 110         fscanf(fp,"%lf",&dtemp);
 111         fscanf(fp,"%lf",&dtemp);
 112         fscanf(fp,"%lf",&dtemp);
 113
 114         id=AtomIndicator(buffer);
 115         Atomdic[iatom]=id;
 116       }
 117       if (i == Natom) goto Ending;
 118   }
 119 }
 120 Ending:;
 121 }
 122
 123
 124
 125 void GetAtomlist(char *name) {
 126 /*This subroutine reads the head of the configuration file to set the lattice parameters and the number of atoms.  It also checks that a full frame exists because it loops over the first frame to populate the Atomdic with the proper type for each atom id.*/
 127
 128   FILE *fp;
 129   int i,iatom,count=0,itemp,id,flag;
 130   double dtemp;
 131   double a,b,c,d,e,f,g,h,ii,det;
 132   char buffer[1000];
 133
 134   FILE *temp;
 135   fp=fopen(name,"r");
 136   FileOpenError(fp,name);
 137
 138   while(fscanf(fp,"%s",buffer)!=EOF) {
 139     if(strcmp(buffer,"vectors:")==0) {
 140       fscanf(fp,"%s",buffer);
 141       for(i=0;i<3;i++) fscanf(fp,"%lf",&latc_x[i]);
 142       fscanf(fp,"%s",buffer);
 143       for(i=0;i<3;i++) fscanf(fp,"%lf",&latc_y[i]);
 144       fscanf(fp,"%s",buffer);
 145       for(i=0;i<3;i++) fscanf(fp,"%lf",&latc_z[i]);
 146     }
 147     else if(strcmp(buffer,"origin:")==0) {
 148       fscanf(fp,"%s",buffer);
 149       for(i=0;i<3;i++) fscanf(fp,"%lf",&Origin[i]);
 150     }
 151     else if(strcmp(buffer,"particles")==0) {
 152       fscanf(fp,"%d",&Natom);
 153       printf("# of atoms : %d\n",Natom);
 154       Allocation();
 155     }
 156     else if(strcmp(buffer,"positions")==0) {
 157       for(i=0,flag=0;i<Natom;i++) {
 158         fscanf(fp,"%d",&id);
 159         id--;
 160         fscanf(fp,"%s",buffer);
 161         iatom=AtomIndicator(buffer);
 162         Atomdic[id]=iatom;
 163         count++;
 164 //        fscanf(fp,"%d",&itemp);
 165 //        fprintf(temp,"\t%d",itemp);
 166 //        fscanf(fp,"%lf",&dtemp);
 167 //        fprintf(temp,"\t%lf",dtemp);
 168 //        fscanf(fp,"%lf",&dtemp);
 169 //        fprintf(temp,"\t%lf",dtemp);
 170 //        fscanf(fp,"%lf",&dtemp);
 171 //        fprintf(temp,"\t%lf",dtemp);
 172 //        fscanf(fp,"%lf",&dtemp);
 173 //        fprintf(temp,"\t%lf\n",dtemp);
 174         fgets(buffer,1000,fp);
 175         flag=1;
 176 /*         fscanf(fp,"%lf",&dtemp); */
 177 /*         fscanf(fp,"%lf",&dtemp); */
 178 /*         fscanf(fp,"%lf",&dtemp); */
 179 /*         fscanf(fp,"%lf",&dtemp); */
 180       }
 181       if(flag==1) break;
 182     }
 183   }
 184   if(Natom!=count) {
 185     printf("Can't read whole # of atoms.\n");
 186     exit(0);
 187   }
 188
 189
 190   a=latc_x[0]; b=latc_y[0]; c=latc_z[0];
 191   d=latc_x[1]; e=latc_y[1]; f=latc_z[1];
 192   g=latc_x[2]; h=latc_y[2]; ii=latc_z[2];
 193   det=a*e*ii-a*f*h-d*b*ii+d*c*h+g*b*f-g*c*e;
 194   if (det != 0) {
 195   Ilatc_x[0]=(e*ii-f*h)/det;
 196   Ilatc_x[1]=-(b*ii-c*h)/det;
 197   Ilatc_x[2]=(b*f-c*e)/det;
 198
 199   Ilatc_y[0]=-(d*ii-f*g)/det;
 200   Ilatc_y[1]=(a*ii-c*g)/det;
 201   Ilatc_y[2]=-(a*f-c*d)/det;
 202
 203   Ilatc_z[0]=(e*g-d*h)/det;
 204   Ilatc_z[1]=-(b*g-a*h)/det;
 205   Ilatc_z[2]=(a*e-b*d)/det;
 206   }
 207
 208   close(fp);
 209 }
 210
 211 void GetCOdict(char *name) {
 212 /* This subroutine populates the cut-off dictionary for each bond and throws an error if any possible combination of atoms do not have a specified bond*/
 213   FILE *fp;
 214   int i,j,iatom,jatom;
 215   double cut;
 216 #define BUFLEN 1000
 217   char buffer[BUFLEN];
 218   fp=fopen(name,"r");
 219   FileOpenError(fp,name);
 220
 221   for(i=0;i<Natomtypes;i++) {
 222     for(j=0;j<Natomtypes;j++)  COdic[i][j]=-1;
 223   }
 224
 225   while(fscanf(fp,"%s",buffer)!=EOF) {
 226     if (strcmp(buffer,"#") == 0) {
 227       SkipLine(fp);
 228       continue;
 229     }
 230     iatom=AtomIndicator(buffer);
 231     fscanf(fp,"%s",buffer);
 232     jatom=AtomIndicator(buffer);
 233     fscanf(fp,"%lf",&cut);
 234     COdic[iatom][jatom]=cut;
 235     COdic[jatom][iatom]=cut;
 236   }
 237   for(i=0;i<Natomtypes;i++) {
 238     for(j=0;j<Natomtypes;j++) {
 239       if(COdic[i][j]==-1) {
 240         printf("Can't fine a certain pair for cut off dictionary!\n");
 241         exit(0);
 242       }
 243     }
 244   }
 245
 246
 247   close(fp);
 248 }
 249 void AnalyzeTrajectory(char *name,char *name2) {
 250 /***************************************************************************/
 251 /*  This subroutine loops through the data both with the bond lists and the position files.
 252     Step 1: It checks the timestep for the bond list and then reads in the position data ensuring that the timesteps match.
 253     Step 2: It loops over the bonds to identify molecules
 254     Step 3: It loops over the molecules to determine name and species.
 255     Step 4: It calculates and writes dipole.out and vibtemp.out.
 256     Step 5: It writes num.out */
 257 /**************************************************************************/
 258
 259   FILE *fp,*out3;//fp is the file with the bond lists.  out3 refers to a scratch file to keep track of how many of each type of molecules exist.
 260   int i,j,k,l,iatom,jatom,Nlink,count,iatomtype,jatomtype;
 261 /* jatom is the neighbor, iatom is the self, Nlink is the number of neighbors, count is a running index to identify order of atoms (e.g. atom 5 is on line 38)*/
 262   int frame;
 263   int itemp,Molid;
 264 /*itemp is a generic variable as needed, Molid is the line of the atom (e.g. atom 5 is on line 38 means Molid[4]=37 because of the zero).*/
 265
 266   double BO,dtemp;//BO is bond order
 267   char buffer[1000];
 268   int Cnct[MaxNumBonds];//connectivity table for a given atom
 269
 270   int Nspec,flag_mol,flag_spec,flag_identity,nd,timestep;
 271 /*  Nspec is the number of species during this frame.  flag_mol is a binary flag with 1 indicating that this atom belongs to this molecule.  flag_spec is a binary flag with 1 indicating that a new species has been found.  I have no idea what purpose flag_identity serves, but a 1 indicates that flag_spec is 1.
 272 nd is a return value from CheckExistence.*/
 273
 274   int Name[Natomtypes],MolTemp[Natom];
 275 /*Name is a counter for how many atoms of a given type are in this particular molecule. MolTemp is an array that keeps track of molecule id by atom id.  For example, if atom 5 belongs to molecule 3 MolTemp[4]=3*/
 276
 277   int CheckExistence();
 278   void GetAtomdata();
 279   void GetCOMDipole();
 280
 281   fp=fopen(name,"r");
 282   FileOpenError(fp,name);
 283   out3=fopen(name2,"w");
 284   FileOpenError(out3,name2);
 285
 286   frame=0;
 287
 288   while(fscanf(fp,"%s",buffer)!=EOF) {
 289 /* Step 1 */
 290     if(strcmp(buffer,"Timestep")==0) {
 291       fscanf(fp,"%d",&timestep);
 292       //      GetAtomdata(timestep);
 293     }
 294
 295 /* Step 2 */
 296     if(strcmp(buffer,"q")==0) {
 297       frame++;
 298       count=0;
 299 /*             zeroing out the counters        */
 300       for(i=0;i<Natom;i++) {
 301         MolTemp[i]=-1;
 302         NMol[i]=1;
 303       }
 304 /*          Step 2A reads the atomid as iatom, the first itemp is a discard of the atom type, Nlink is the number of neighbors, and then the loop over all the bonded neighbors is stored in Cnct, and the last itemp is the ignored 0 to mark the end.*/
 305       for(i=0;i<Natom;i++) {
 306         fscanf(fp,"%d",&iatom);
 307         iatom--;
 308         if(MolTemp[iatom]==-1) {
 309           MolTemp[iatom]=count;
 310           count++;
 311         }
 312         fscanf(fp,"%d",&itemp);
 313         fscanf(fp,"%d",&Nlink);
 314         for(j=0;j<Nlink;j++) {
 315           fscanf(fp,"%d",&itemp);
 316           itemp--;
 317           Cnct[j]=itemp;
 318         }
 319         fscanf(fp,"%d",&itemp);
 320
 321 /*         Step 2B reads the bond orders from the bond file and assigns atoms to molecules.  The first check is if a bond exists.  The second check is if i and j already belong to the same molecule.  If they don't and j has no molecule, assign it to i's molecule.  If they do, then j was assigned to a different number so loop through the molecule list and set all the atoms on that molecule to the new i molecule.  Then read and discard the remainder of the line.*/
 322
 323         for(j=0;j<Nlink;j++) {
 324           fscanf(fp,"%lf",&BO);
 325           jatom=Cnct[j];
 326           Molid=MolTemp[jatom];
 327           iatomtype=Atomdic[iatom];
 328           jatomtype=Atomdic[jatom];
 329           if(BO>COdic[iatomtype][jatomtype]) {
 330             if(Molid!=MolTemp[iatom]) {
 331               if(MolTemp[jatom]==-1) MolTemp[jatom]=MolTemp[iatom];
 332               else {
 333                 for(k=0;k<Natom;k++) {
 334                   if(MolTemp[k]==Molid) MolTemp[k]=MolTemp[iatom];
 335                 }
 336               }
 337             }
 338           }
 339         }
 340         fscanf(fp,"%lf",&dtemp);
 341         fscanf(fp,"%lf",&dtemp);
 342         fscanf(fp,"%lf",&dtemp);
 343       }
 344
 345 /*       Step 3 builds the chemical formula for each molecule.  For each molecule (the for loop over i with count as the maximum), loop over all the atoms.*/
 346
 347       for(i=0;i<MaxMolTypes;i++) Netc[i]=0.0;
 348       for(i=0,Nspec=0,Nmolecule=0;i<count;i++) {
 349
 350  /*     Step 3A: if the atom belongs to this molecule, increment the appropriate type counter (e.g. we have another carbon), flip flag_mol to 1, set the Mol value to this molecule, and add the atomic charge to the running charge total. */
 351
 352         for(j=0;j<Natomtypes;j++) Name[j]=0;
 353         for(j=0,flag_mol=0;j<Natom;j++) {
 354           if(MolTemp[j]==i) {
 355             jatom=Atomdic[j];
 356             Name[jatom]++;
 357             flag_mol=1;
 358             Mol[j]=Nmolecule;
 359             Netc[Nmolecule]+=Charge[j];
 360           }
 361         }
 362
 363 /*     Step 3B: If this molecule contains any atoms (flag_mol==1), then determine if it is a new species.  MolName is populated as new species are found.  In one of the least efficient matching mechanisms ever, we walk through all the known species names.  If no match is found, the flag_spec is set to 1.  A zero flag_spec means that we have a match, increment the appropriate species counter.  Otherwise, if we have a zero number of species or a flipped flag_spec, we have a new species and need to populate its name.  Then we increment the number of species as appropriate and always add a new molecule.*/
 364
 365         if(flag_mol==1) {
 366           flag_identity=1;
 367           for(k=0;k<Nspec;k++) {
 368             flag_spec=0;
 369             for(l=0;l<Natomtypes;l++) {
 370               if(MolName[Natomtypes*k+l]!=Name[l]) flag_spec=1;
 371             }
 372             if(flag_spec==0) NMol[k]++;
 373             flag_identity*=flag_spec;
 374           }
 375           if(Nspec==0 || flag_identity==1) {
 376             for(l=0;l<Natomtypes;l++) {
 377               MolName[Natomtypes*Nspec+l]=Name[l];
 378             }
 379             Nspec++;
 380           }
 381           Nmolecule++;
 382         }//ends the if(flag_mol==1) loop
 383
 384       }//ends the for i up to count loop
 385
 386 /*  Step 4 does the end of frame cleanup including the writes for COM */
 387       fprintf(out3,"%d\t",timestep);
 388       printf("Timestep = %d\n",timestep);
 389       fprintf(dpol,"%d\n",Nmolecule);
 390       fprintf(dpol,"TSTEP %d\n",timestep);
 391       fprintf(vtemp,"%d\n",Nmolecule);
 392       fprintf(vtemp,"TSTEP %d\n",timestep);
 393       fprintf(conc,"TSTEP %d\n",timestep);
 394       fprintf(conc,"%d\n",Nmolecule);
 395       GetCOMDipole(Nmolecule,timestep,frame);
 396
 397 /*  Step 5 prints the list of species and their quantities   */
 398       fprintf(out3,"%d\t",Nspec);
 399       for(i=0;i<Nspec;i++) {
 400         nd=CheckExistence(i);
 401         fprintf(out3,"%d ",nd);
 402         fprintf(out3,"%d\t",NMol[i]);
 403       }
 404       fprintf(out3,"\n");
 405     }
 406   }
 407   fflush(out3);
 408   fclose(out3);
 409   fclose(fp);
 410 }
 411
 412 int CheckExistence(int id) {
 413 /* The purpose of this function is to have a list of only the species that occur in this frame. */
 414
 415   int i,j,num,molid,itemp;
 416   int flag,flag_identity=1;
 417
 418   for(i=0;i<Nmoltype;i++) {
 419     flag=0;
 420     for(j=0;j<Natomtypes;j++) {
 421       molid=MolType[Natomtypes*i+j];
 422       if(molid!=MolName[Natomtypes*id+j]) flag=1;
 423     }
 424     if(flag==0) return i;
 425     flag_identity*=flag;
 426   }
 427
 428   for(i=0;i<Natomtypes;i++) {
 429     MolType[Natomtypes*Nmoltype+i]=MolName[Natomtypes*id+i];
 430   }
 431   Nmoltype++;
 432
 433   return Nmoltype-1;
 434 }
 435
 436 void FileOpenError(FILE *fp,char *name) {
 437   if(fp==NULL) {
 438     printf("file open error(%s).\n",name);
 439     exit(0);
 440   }
 441 }
 442
 443 void Allocation() {
 444 /* This subroutine sets up most of the important adjustable arrays based on the number of atoms.  It also opens the dipole and vibtemp files for writing.*/
 445   int flag=0;
 446   char name[100];
 447   Atomdic=(int *)malloc(Natom*sizeof(int));
 448   Mol=(int *)malloc(Natom*sizeof(int));
 449   MolName=(int *)malloc(Natom*Natomtypes*sizeof(int));
 450 //20Jun07 JLB: I have made the absurdly huge array much smaller
 451 //  MolType=(int *)malloc(Natom*MaxMolTypes*sizeof(int));
 452   MolType=(int *)malloc((Natomtypes+20)*MaxMolTypes*sizeof(int));
 453   NMol=(int *)malloc(Natom*sizeof(int));
 454   rx=(double *)malloc(Natom*sizeof(double));
 455   ry=(double *)malloc(Natom*sizeof(double));
 456   rz=(double *)malloc(Natom*sizeof(double));
 457   vx=(double *)malloc(Natom*sizeof(double));
 458   vy=(double *)malloc(Natom*sizeof(double));
 459   vz=(double *)malloc(Natom*sizeof(double));
 460   Charge=(double *)malloc(Natom*sizeof(double));
 461   Netc=(double *)malloc(MaxMolTypes*sizeof(double));
 462
 463   if(Atomdic==NULL || NMol==NULL || Mol==NULL || MolType==NULL) flag=1;
 464   if(rx==NULL || ry==NULL || rz==NULL || Charge==NULL) flag=1;
 465   if(vx==NULL || vy==NULL || vz==NULL) flag=1;
 466   if(Netc==NULL) flag=1;
 467
 468   if(flag==1) {
 469     printf("malloc error.\n");
 470     exit(0);
 471   }
 472   //  trj=fopen(Atomlist,"r");
 473   //  FileOpenError(trj,Atomlist);
 474   dpol=fopen(Output4,"w");
 475   FileOpenError(dpol,Output4);
 476   vtemp=fopen(Output5,"w");
 477   FileOpenError(vtemp,Output5);
 478   conc=fopen(Output6,"w");
 479   FileOpenError(conc,Output6);
 480 }
 481
 482 int AtomIndicator(char *symbol) {
 483 /* This subroutine identifies the atom type from the config file and sets the masses.  Why set the masses here, I don't know.*/
 484
 485   int value;
 486   Mass[Hydrogen]=1.0079;
 487   Mass[Carbon]=12.0107;
 488   Mass[Oxygen]=15.9994;
 489   Mass[Nitrogen]=14.0067;
 490   if(strcmp(symbol,"H")==0) value=Hydrogen;
 491   else if(strcmp(symbol,"C")==0) value=Carbon;
 492   else if(strcmp(symbol,"O")==0) value=Oxygen;
 493   else if(strcmp(symbol,"N")==0) value=Nitrogen;
 494   else {
 495     printf("Can't recognize the atom type : %s.\n",symbol);
 496     exit(0);
 497   }
 498
 499   return value;
 500 }
 501 void PostProcess(char *name,char *name2,char *name3) {
 502
 503 /*  This subroutine does the global lists after all the frames have been read.  out is the fra.out with the list of species and number of each per frame.  out2 is mol.out with the chemical identities of each species. out3 is num.out which is the scratch file with number of species per frame in a non-plottable format.*/
 504
 505   int i,j;
 506   FILE *out,*out2,*out3;
 507   int step,molnum,id,num,itemp;
 508   int numcount[Nmoltype];
 509
 510   out=fopen(name,"w");
 511   FileOpenError(out,name);
 512   out2=fopen(name2,"w");
 513   FileOpenError(out2,name2);
 514   out3=fopen(name3,"r");
 515   FileOpenError(out3,name3);
 516
 517   fprintf(out,"Timestep\t");
 518   fprintf(out2,"%d\n\n",Nmoltype);
 519
 520 /*  To get the chemical identity of a species, this process loops over each atom type and writes the number of atoms of that type in the formula.    */
 521   for(i=0;i<Nmoltype;i++) {
 522     fprintf(out2,"%d\t",i);
 523     for(j=0;j<Natomtypes;j++) {
 524       itemp=MolType[Natomtypes*i+j];
 525       if(itemp!=0) {
 526         if(j==Carbon) {
 527           fprintf(out,"C");
 528           fprintf(out2,"C");
 529         }
 530         else if(j==Hydrogen) {
 531           fprintf(out,"H");
 532           fprintf(out2,"H");
 533         }
 534         else if(j==Oxygen) {
 535           fprintf(out,"O");
 536           fprintf(out2,"O");
 537         }
 538         else if(j==Nitrogen) {
 539           fprintf(out,"N");
 540           fprintf(out2,"N");
 541         }
 542         if(itemp!=1) {
 543           fprintf(out,"%d",itemp);
 544           fprintf(out2,"%d",itemp);
 545         }
 546       }
 547     }
 548     fprintf(out,"\t");
 549     fprintf(out2,"\n");
 550   }
 551   fprintf(out,"\n");
 552
 553 /*    for each molecule type, read the num.out scratch file and find the number of those molecules.  Then write it out in a nice column format suitable for plotting. */
 554   while(fscanf(out3,"%d",&step)!=EOF) {
 555     fscanf(out3,"%d",&num);
 556     for(i=0;i<Nmoltype;i++) numcount[i]=0;
 557     for(i=0;i<num;i++) {
 558       fscanf(out3,"%d",&id);
 559       fscanf(out3,"%d",&molnum);
 560       numcount[id]=molnum;
 561     }
 562     fprintf(out,"%d\t",step);
 563     for(i=0;i<Nmoltype;i++) {
 564       fprintf(out,"%d\t",numcount[i]);
 565     }
 566     fprintf(out,"\n");
 567   }
 568 }
 569 void GetAtomdata(int tstep) {
 570 /* This subroutine reads in the position data for a frame*/
 571
 572   int i,iatom,count=0,itemp,id,flag;
 573   double dtemp;
 574   int step;
 575   char buffer[1000];
 576
 577   while(fscanf(trj,"%s",buffer)!=EOF) {
 578     if(strcmp(buffer,"Timestep")==0) {
 579       fscanf(trj,"%d",&step);
 580       if(step!=tstep) {
 581         printf("The sequences are not matched between bond data file and coordinate file\n");
 582         printf("step in config file : %d\tStep in bond file : %d\n",tstep,step);
 583         exit(0);
 584       }
 585     }
 586     else if(strcmp(buffer,"positions")==0) {
 587       for(i=0,flag=0;i<Natom;i++) {
 588         fscanf(trj,"%d",&id);
 589         id--;
 590         fscanf(trj,"%s",buffer);
 591         iatom=AtomIndicator(buffer);
 592         Atomdic[id]=iatom;
 593         fscanf(trj,"%d",&itemp);
 594         fscanf(trj,"%lf",&dtemp);
 595         Charge[id]=dtemp;
 596         fscanf(trj,"%lf",&dtemp);
 597         rx[id]=dtemp;
 598         fscanf(trj,"%lf",&dtemp);
 599         ry[id]=dtemp;
 600         fscanf(trj,"%lf",&dtemp);
 601         rz[id]=dtemp;
 602         fscanf(trj,"%lf",&dtemp);
 603         vx[id]=dtemp;
 604         fscanf(trj,"%lf",&dtemp);
 605         vy[id]=dtemp;
 606         fscanf(trj,"%lf",&dtemp);
 607         vz[id]=dtemp;
 608         fscanf(trj,"%lf",&dtemp);
 609         flag=1;
 610       }
 611       if(flag==1) break;
 612     }
 613   }
 614 }
 615 void GetCOMDipole(int max,int tstep,int frame) {
 616 /* max is the number of molecules */
 617   int i;
 618   void AnalCOMDipole();
 619
 620   for(i=0;i<max;i++) AnalCOMDipole(i);
 621
 622 }
 623 void AnalCOMDipole(int molid) {
 624
 625   int i,j,ix,iy,iz,count,itemp;
 626   int ntot;
 627   double mass,chg;
 628   int atomid;
 629   double totmass;
 630   double x0,y0,z0,xtemp,ytemp,ztemp,x,y,z,r2,minr2,dx,dy,dz;
 631   double COMx,COMy,COMz;
 632   double vCOMx,vCOMy,vCOMz;
 633   double Dpolx,Dpoly,Dpolz,Dpol;
 634   double VibTempx,VibTempy,VibTempz,VibTemp;
 635   int Name[Natomtypes];
 636
 637   double tempx,tempy,tempz;
 638
 639   totmass=0.0;
 640   vCOMx=vCOMy=vCOMz=0.0;
 641
 642 /* Again in an inefficient way, we loop over all the atoms to find the ones that are members of this molecule (molid). Then calculate the desired quantities, including the fact that the first atom sets our standard for measurement.  Then, we write the quantities to the COM files.
 643
 644 JB 26Sep07: to facilitate my write out of the chemical formula, redo the molecule determination here. */
 645   for(j=0;j<Natomtypes;j++) Name[j]=0;
 646   for(i=0,count=0;i<Natom;i++) {
 647     if(Mol[i]==molid) {
 648       atomid=Atomdic[i];
 649       Name[atomid]++;
 650       mass=Mass[atomid];
 651       chg=Charge[i];
 652       totmass+=mass;
 653       vCOMx+=mass*vx[i];
 654       vCOMy+=mass*vy[i];
 655       vCOMz+=mass*vz[i];
 656       if(count==0) {
 657 //        printf("i,  rx[i], ry[i], rz[i] , mass %d %e %e %e %e\n", i,  rx[i], ry[i], rz[i] , mass);
 658         x0=rx[i];
 659         y0=ry[i];
 660         z0=rz[i];
 661         COMx=mass*x0;
 662         COMy=mass*y0;
 663         COMz=mass*z0;
 664       }
 665       if(count!=0) {
 666         minr2=1.0e8;
 667         for(ix=-1;ix<2;ix++) {
 668           for(iy=-1;iy<2;iy++) {
 669             for(iz=-1;iz<2;iz++) {
 670               xtemp=rx[i]+ix*latc_x[0]+iy*latc_y[0]+iz*latc_z[0];
 671               ytemp=ry[i]+ix*latc_x[1]+iy*latc_y[1]+iz*latc_z[1];
 672               ztemp=rz[i]+ix*latc_x[2]+iy*latc_y[2]+iz*latc_z[2];
 673               dx=xtemp-x0;
 674               dy=ytemp-y0;
 675               dz=ztemp-z0;
 676               r2=dx*dx+dy*dy+dz*dz;
 677               if(r2<minr2) {
 678                 minr2=r2;
 679                 x=xtemp;
 680                 y=ytemp;
 681                 z=ztemp;
 682               }
 683             }
 684           }
 685         }
 686         rx[i]=x;
 687         ry[i]=y;
 688         rz[i]=z;
 689 //        printf("AAA i,  rx[i], ry[i], rz[i] , mass %d %e %e %e %e\n", i,  rx[i], ry[i], rz[i] , mass);
 690         COMx+=mass*x;
 691         COMy+=mass*y;
 692         COMz+=mass*z;
 693       }
 694       count++;
 695     }
 696   }
 697 //  printf ("AAA totmass %e\n",totmass);
 698   COMx/=totmass;
 699   COMy/=totmass;
 700   COMz/=totmass;
 701
 702   vCOMx/=totmass;
 703   vCOMy/=totmass;
 704   vCOMz/=totmass;
 705
 706   VibTempx=VibTempy=VibTempz=0.0;
 707   for(i=0,ntot=0,Dpolx=0.0,Dpoly=0.0,Dpolz=0.0;i<Natom;i++) {
 708     if(Mol[i]==molid) {
 709       atomid=Atomdic[i];
 710       mass=Mass[atomid];
 711       chg=Charge[i];
 712       totmass+=mass;
 713       x=rx[i]-COMx;
 714       y=ry[i]-COMy;
 715       z=rz[i]-COMz;
 716       Dpolx+=chg*x;
 717       Dpoly+=chg*y;
 718       Dpolz+=chg*z;
 719       VibTempx+=mass*(vx[i]-vCOMx)*(vx[i]-vCOMx);
 720       VibTempy+=mass*(vy[i]-vCOMy)*(vy[i]-vCOMy);
 721       VibTempz+=mass*(vz[i]-vCOMz)*(vz[i]-vCOMz);
 722       ntot++;
 723     }
 724   }
 725   VibTempx/=(double)ntot;
 726   VibTempy/=(double)ntot;
 727   VibTempz/=(double)ntot;
 728
 729   VibTemp=(VibTempx+VibTempy+VibTempz)/3.0;
 730
 731   VibTemp*=1.2027E+06;
 732
 733
 734   COMx-=Origin[0]; COMy-=Origin[1]; COMz-=Origin[2];
 735 //  printf("VVVV %e %e %e\n",Ilatc_x[0],Ilatc_x[1],Ilatc_x[2]);
 736   tempx=Ilatc_x[0]*COMx;
 737   tempx+=Ilatc_x[1]*COMy;
 738   tempx+=Ilatc_x[2]*COMz;
 739
 740   tempy=Ilatc_y[0]*COMx;
 741   tempy+=Ilatc_y[1]*COMy;
 742   tempy+=Ilatc_y[2]*COMz;
 743
 744   tempz=Ilatc_z[0]*COMx;
 745   tempz+=Ilatc_z[1]*COMy;
 746   tempz+=Ilatc_z[2]*COMz;
 747
 748   tempx-=floor(tempx);
 749   tempy-=floor(tempy);
 750   tempz-=floor(tempz);
 751
 752   COMx=latc_x[0]*tempx+latc_y[0]*tempy+latc_z[0]*tempz;
 753   COMy=latc_x[1]*tempx+latc_y[1]*tempy+latc_z[1]*tempz;
 754   COMz=latc_x[2]*tempx+latc_y[2]*tempy+latc_z[2]*tempz;
 755 //  printf ("LATICE %e %e %e %e %e %e \n",latc_x[0], tempx, latc_y[0], tempy,latc_z[0],tempz );
 756   COMx+=Origin[0]; COMy+=Origin[1]; COMz+=Origin[2];
 757
 758   Dpolx/=0.2082267;
 759   Dpoly/=0.2082267;
 760   Dpolz/=0.2082267;
 761
 762   Dpol=Dpolx*Dpolx+Dpoly*Dpoly+Dpolz*Dpolz;
 763   Dpol=sqrt(Dpol);
 764 //  printf("BBB %e %e %e\n",COMx,COMy,COMz);
 765   fprintf(dpol,"He %e %e %e ",COMx,COMy,COMz);
 766   fprintf(dpol,"%e %e %e %e\n",Dpol,Dpolx,Dpoly,Dpolz);
 767   fprintf(vtemp,"He %e %e %e",COMx,COMy,COMz);
 768   fprintf(vtemp,"%e %e %e %e\n",VibTemp,vCOMx,vCOMy,vCOMz);
 769
 770 /* JB 26Sep07: I am adding the routines that will write the chemical formula and the COM to a separate file for use with the composition profile.  Step 1: determine the chemical formula for this molecule.*/
 771     for(j=0;j<Natomtypes;j++) {
 772       itemp=Name[j];
 773       if(itemp!=0) {
 774         if(j==Carbon) {
 775           fprintf(conc,"C");
 776         }
 777         else if(j==Hydrogen) {
 778           fprintf(conc,"H");
 779         }
 780         else if(j==Oxygen) {
 781           fprintf(conc,"O");
 782         }
 783         else if(j==Nitrogen) {
 784           fprintf(conc,"N");
 785         }
 786         if(itemp!=1) {
 787           fprintf(conc,"%d",itemp);
 788         }
 789       }
 790     }
 791   fprintf(conc,"\t");
 792   fprintf(conc,"%e %e %e\n",COMx,COMy,COMz);
 793
 794 }
 795
 796 void SkipLine(FILE *fp) {
 797   int c;
 798   while ((c = getc(fp)) != EOF)
 799     if (c == '\n') break;
 800 }
 801