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

tools/msi2lmp/README

+Stephanie Teich-McGoldrick (Sandai) is the current maintainer
+of the msi2lmp tool.  She can be contacted at steichm at sandia.gov
+
+18 Jul 2013 Axel Kohlmeyer <akohlmey@gmail.com>
+
+Added support for writing out image flags
+Improved accuracy of atom masses
+Added flag for shifting the entire system
+Fixed some minor logic bugs and prepared
+for supporting other force fields and morse style bonds.
+
+12 Jul 2013 Axel Kohlmeyer <akohlmey@gmail.com>
+
+Fixed the bug that caused improper coefficients to be wrong
+Cleaned up the handling of box parameters and center the box
+by default around the system/molecule. Added a flag to make
+this step optional and center the box around the origin instead.
+Added a regression test script with examples.
+
+1 Jul 2013 Axel Kohlmeyer <akohlmey@gmail.com>
+
+Cleanup and improved port to windows.
+Removed some more static string limits.
+Added print level 3 for additional output.
+Make code stop at missing force field parameters
+and added -i flag to override this.
+Safer argument checking.
+Provide short versions for all flags.
+
+23 Sep 2011
+
+added support for triclinic boxes
+see msi2lmp/TriclinicModification.pdf doc for details
+
+-----------------------------
+
+ msi2lmp V3.6 4/10/2005
+
+ This program uses the .car and .mdf files from MSI/Biosyms's INSIGHT
+ program to produce a LAMMPS data file.
+
+ 1. Building msi2lmp3
+
+    Use the Makefile in the src directory. It is
+    currently set up for gcc. One will have to modify
+    it to use a different compiler.
+
+ 2. Testing the program
+
+    There are three pairs (.car and .mdf) files in the
+    test directory: crambin, nylon and phen3_cff97. The
+    atom types in crambin and nylon are cvff (Class I) atom
+    types and those in phen3_cff97 are cff9x (Class II) atom types.
+    Two forcefield files, cvff.frc and cff91.frc, are needed
+    generate lammps data files for these three test files. To
+    run it you would:
+
+   % setenv BIOSYM_LIBRARY ../biosym_frc_files
+   % ../src/msi2lmp.exe nylon -class I -frc cvff > data.nylon
+   % ../src/msi2lmp.exe crambin -class I -frc cvff > data.crambin
+   % ../src/msi2lmp.exe phen3_cff97 -class II -frc cff91 > data.phen3_cff97
+
+   Three files should be generated: data.nylon, data.crambin
+   and data.phen3_cff97. These can be compared against 
+   data.x in the directory correct. If there are differences, 
+   first recompile the program with no optimization and try again.
+   If there are still differences, send email to jec@mayo.edu
+
+   Note: you will see many "Unable to find..." parameters messages
+         in the phen3_cff97 test case. Most of those parameters 
+         exist in cff95.frc, but not in cff91.frc
+
+ 3. To run the program
+
+   The program is started by supplying information at the command prompt
+   according to the usage described below.  
+
+   USAGE: msi2lmp.exe ROOTNAME {-2001} {-print #} {-class #} {-frc FRC_FILE} 
+
+   -- msi2lmp.exe is the name of the executable
+   -- ROOTNAME is the base name of the .car and .mdf files
+   -- -2001
+         Output lammps files for LAMMPS version 2001 (F90 version)
+         Default is to write output for the C++ version of LAMMPS
+
+   -- -print (or -p)
+	 # is the print level  0 - silent except for error messages
+	                       1 - minimal (default)
+                               2 - verbose (usual for developing and
+                                   checking new data files for consistency)
+                               3 - even more verbose (additional debug info)
+
+   -- -ignore (or -i)   ignore errors about missing force field parameters
+                        and treat them as warnings instead.
+
+   -- -class  (or -c)
+        # is the class of forcefield to use (I  or 1 = Class I e.g., CVFF)
+                                             (II or 2 = Class II e.g., CFFx)
+        default is -class I
+
+   -- -frc    (or -f) specifies name of the forcefield file (e.g., cff91)
+ 
+     If the file name includes a directory component (or drive letter on Windows),
+     then the name is used as is. Otherwise, the program looks for the forcefield
+     file in $BIOSYM_LIBRARY (or %BIOSYM_LIBRARY% on Windows).
+     If $BIOSYM_LIBRARY is not set, ../biosym_frc_files is used (for testing).
+     If the file name does not end in .frc, then .frc is appended to the name.
+
+     For example,  -frc cvff (assumes cvff.frc is in $BIOSYM_LIBRARY
+                              or ../biosym_frc_files)
+
+                   -frc cff/cff91 (assumes cff91.frc is in cff)
+
+                   -frc /usr/local/biosym/forcefields/cff95
+                       (assumes cff95.frc is in /usr/local/biosym/forcefields/)
+
+     By default, the program uses $BIOSYM_LIBRARY/cvff.frc or ../biosym_frc_files/cvff.frc
+
+  -- the LAMMPS data file is written to ROOTNAME.lammps{01/05},
+     protocol and error information is written to the screen.
+
+****************************************************************
+*
+* Msi2lmp3
+*
+* This is the third version of a program that generates a LAMMPS
+* data file based on the information in MSI .car (atom
+* coordinates), .mdf (molecular topology) and .frc (forcefield) 
+* files. The .car and .mdf files are specific to a molecular
+* system while the .frc file is specific to a forcefield version.
+* The only coherency needed between .frc and .car/.mdf files are
+* the atom types. 
+*
+* The first version was written by Steve Lustig at Dupont, but
+* required using Discover to derive internal coordinates and
+* forcefield parameters
+*
+* The second version was written by Michael Peachey while an
+* intern in the Cray Chemistry Applications Group managed
+* by John Carpenter. This version derived internal coordinates
+* from the mdf file and looked up parameters in the frc file
+* thus eliminating the need for Discover.
+*
+* The third version was written by John Carpenter to optimize
+* the performance of the program for large molecular systems
+* (the original  code for deriving atom numbers was quadratic in time)
+* and to make the program fully dynamic. The second version used
+* fixed dimension arrays for the internal coordinates.
+*
+* The current maintainer is only reluctantly doing so because John Mayo no longer
+* needs this code.
+*
+* V3.2 corresponds to adding code to MakeLists.c to gracefully deal with
+* systems that may only be molecules of 1 to 3 atoms. In V3.1, the values
+* for number_of_dihedrals, etc. could be unpredictable in these systems.
+*
+* V3.3 was generated in response to a strange error reading a MDF file generated by
+* Accelys' Materials Studio GUI. Simply rewriting the input part of ReadMdfFile.c 
+* seems to have fixed the problem.
+*
+* V3.4 and V3.5 are minor upgrades to fix bugs associated mostly with .car and .mdf files
+* written by Accelys' Materials Studio GUI.
+*
+* V3.6 outputs to LAMMPS 2005 (C++ version).
+*
+* Contact: Kelly L. Anderson, kelly.anderson@cantab.net
+* 
+* April 2005

tools/msi2lmp/biosym_frc_files/README

+
+Copy cvff.frc and cff9*.frc or any other
+.frc files you have licensed to this 
+directory for use with msi2lmp. 
+

tools/msi2lmp/biosym_frc_files/clayff.frc

+!CLAYFF forcefield
+
+#atom_types	cvff
+
+!Ver  Ref  Type    Mass      Element  Connections   Comment
+!---- ---  ----  ----------  -------  -----------------------------------------
+ 1.0   1    st     28.08550    Si          4  
+ 1.0   1    ao     26.98154    Al          6     
+ 1.0   1    at     26.98154    Al          4   
+ 1.0   1   mgo     24.30500    Mg          6   
+ 1.0   1   cao     40.08000    Ca          6
+ 1.0   1   feo     55.84700    Fe          6
+ 1.0   1   lio     6.941000    Li          6
+ 1.0   1    ob    15.99940      O          2
+ 1.0   1  obss     15.99940     O          3   
+ 1.0   1  obts     15.99940     O          2   
+ 1.0   1  obos     15.99940     O          2   
+ 1.0   1   ohs     15.99940     O          2   
+ 1.0   1    oh     15.99940     O          2
+ 1.0   1   oh-     15.99940     O          1
+ 1.0   1    o*     15.99940     O          2
+ 1.0   1    ho     1.007970     H          1   
+ 1.0   1    h*     1.007970     H          1
+ 1.0   1    Na     22.99000    Na          0
+ 1.0   1     K     39.10        K          0 
+ 1.0   1    Cs     132.9100    Cs          0
+ 1.0   1    Ca     40.07980    Ca          0
+ 1.0   1    Ba     137.3300    Ba          0 
+ 1.0   1    Mg      24.3050    Mg          0
+ 1.0   1    Sr      87.6200    Sr          0
+ 1.0   1    Pb     207.2000    Pb          0
+ 1.0   1    Cl     35.45300    Cl          0
+
+
+#equivalence	cvff 
+
+> Equivalence table for any variant of cvff 
+
+!		         	  Equivalences
+!                 -----------------------------------------
+!Ver  Ref   Type  NonB     Bond    Angle    Torsion    OOP
+!---- ---   ----  ----     ----    -----    -------    ----
+ 1.0   1    h     h        h        h        h         h   
+
+
+
+#auto_equivalence	cvff_auto
+
+!		         	  Equivalences
+!                 -----------------------------------------                       
+!Ver  Ref   Type  NonB Bond   Bond     Angle    Angle     Torsion   Torsion      OOP      OOP 
+!                      Inct           End atom Apex atom End Atoms Center Atoms End Atom Center Atom
+!---- ---   ----  ---- ------ ----  ---------- --------- --------- -----------  -------- ----------- 
+ 2.0  18    h     h     h     h_       h_       h_        h_        h_           h_       h_
+
+
+
+#hbond_definition	cvff 
+
+
+#morse_bond	cvff 
+
+> E = D * (1 - exp(-ALPHA*(R - R0)))^2
+
+!Ver  Ref     I     J          R0         D           ALPHA
+!---- ---    ----  ----     -------    --------      -------
+ 2.3  23     no    o-        1.2178    140.2486      2.0000
+
+
+
+#quadratic_bond	cvff 
+
+> E = K2 * (R - R0)^2
+
+!Ver  Ref     I     J          R0         K2    
+!---- ---    ----  ----     -------    -------- 
+ 2.1  28      oh   ho        1.0000    553.9350
+ 2.1  28     ohs   ho        1.0000    553.9350
+
+
+#quadratic_angle	cvff 
+
+> E = K2 * (Theta - Theta0)^2
+
+!Ver  Ref     I     J     K       Theta0         K2        
+!---- ---    ----  ----  ----    --------     -------
+ 2.3  23     cp    cp    c'      120.0000     34.6799
+
+
+
+#torsion_1      cvff_auto
+
+> E = Kphi * [ 1 + cos(n*Phi - Phi0) ]
+
+!Ver  Ref     I     J      K      L           Kphi        n           Phi0
+!---- ---    ----  ----   ----   ----      -------      ------     -------
+ 2.0  18     *     c_     n3n_   *          0.0500        3         0.
+                                                                      
+                                                                      
+
+#out_of_plane	cvff_auto
+
+> E = Kchi * [ 1 + cos(n*Chi - Chi0) ]
+
+!Ver  Ref     I     J     K     L           Kchi        n           Chi0
+!---- ---    ----  ----  ----  ----      -------      ------     -------
+ 2.0  18     *     c'_   *      *        10.0000        2       180.0000
+
+
+#nonbond(12-6)	cvff 
+
+@type A-B
+@combination geometric
+
+> E = Aij/r^12 - Bij/r^6
+> where  Aij = sqrt( Ai * Aj )
+>        Bij = sqrt( Bi * Bj )
+
+!Ver  Ref     I           A             B 
+!---- ---    ----    -----------   -----------
+ 1.0   1    st           12.3645       0.00954
+ 1.0   1    ao          196.1446       0.03230
+ 1.0   1    at           12.3645       0.00954
+ 1.0   1   mgo         1636.3265       0.07688
+ 1.0   1   cao         17814.73	       0.5987 
+ 1.0   1   feo           702.54	       0.0504
+ 1.0   1   lio           112.01	       0.0201
+ 1.0   1    ob       629358.0000     625.50000
+ 1.0   1  obss      629358.0000     625.50000
+ 1.0   1  obts      629358.0000     625.50000
+ 1.0   1  obos      629358.0000     625.50000
+ 1.0   1   ohs      629358.0000     625.50000
+ 1.0   1    oh      629358.0000     625.50000
+ 1.0   1   oh-      629358.0000     625.50000
+ 1.0   1    o*      629358.0000     625.50000
+ 1.0   1    ho        0.00000001       0.00000
+ 1.0   1    h*        0.00000001       0.00000
+ 1.0   1    Na       14763.1719      87.65132
+ 1.0   1    K        754506.86       549.37
+ 1.0   1    Cs       3998193.96	    1264.63
+ 1.0   1    Ca      125966.6068     224.46969
+ 1.0   1    Ba      1799606.56	    582.25
+ 1.0   1    Mg      1369.00	     69.22
+ 1.0   1    Sr      1185860.37	    688.73
+ 1.0   1    Pb      861150.71	    638.08
+ 1.0   1    Cl      21081006.97     2905.31
+
+
+#bond_increments        cvff
+
+!Ver  Ref     I     J       DeltaIJ     DeltaJI
+!---- ---    ----  ----     -------     -------
+ 2.3  23     no    o-        0.1684     -0.1684             
+

tools/msi2lmp/biosym_frc_files/pcff.rlb

+VERSION
+elib

tools/msi2lmp/src/CheckLists.c

+#include "msi2lmp.h"
+
+void CheckLists() {
+  int i;
+
+  for (i=0; i < total_no_bonds; i++) {
+    if ((atoms[bonds[i].members[0]].type != bondtypes[bonds[i].type].types[0])
+        || (atoms[bonds[i].members[1]].type != bondtypes[bonds[i].type].types[1])) {
+      fprintf(stderr,"Warning atom types in bond %d are inconsistent with bond type %d\n",i,bonds[i].type);
+    }
+  }
+
+  for (i=0; i < total_no_angles;i++) {
+    if ((atoms[angles[i].members[0]].type != angletypes[angles[i].type].types[0])
+        || (atoms[angles[i].members[1]].type != angletypes[angles[i].type].types[1])
+        || (atoms[angles[i].members[2]].type != angletypes[angles[i].type].types[2])) {
+      fprintf(stderr,"Warning atom types in angle %d are inconsistent with angle type %d\n", i,angles[i].type);
+    }
+  }
+
+  for (i=0; i < total_no_dihedrals; i++) {
+    if ((atoms[dihedrals[i].members[0]].type != dihedraltypes[dihedrals[i].type].types[0])
+        || (atoms[dihedrals[i].members[1]].type != dihedraltypes[dihedrals[i].type].types[1])
+        || (atoms[dihedrals[i].members[2]].type != dihedraltypes[dihedrals[i].type].types[2])
+        || (atoms[dihedrals[i].members[3]].type != dihedraltypes[dihedrals[i].type].types[3])) {
+      fprintf(stderr,"Warning atom types in dihedral %d are inconsistent with dihedral type %d\n",i,dihedrals[i].type);
+    }
+  }
+
+  for (i=0; i < total_no_oops; i++) {
+
+    if ((atoms[oops[i].members[0]].type != ooptypes[oops[i].type].types[0])
+        || (atoms[oops[i].members[1]].type != ooptypes[oops[i].type].types[1])
+        || (atoms[oops[i].members[2]].type != ooptypes[oops[i].type].types[2])
+        || (atoms[oops[i].members[3]].type != ooptypes[oops[i].type].types[3])) {
+      fprintf(stderr,"Warning atom types in oop %d are inconsistent with oop type %d\n",i,oops[i].type);
+    }
+  }
+}

tools/msi2lmp/src/Forcefield.h

+/******************************
+*
+*   This is the header file for the routine that reads the forcefield file
+*   into memory in order to speed up searching.
+*
+*   It defines the data structures used to store the force field in memory
+*/
+
+#define MAX_NO_MEMS 6
+#define MAX_NO_PARAMS 8
+
+struct FrcFieldData {
+  float  ver;                /* Version number of forcefield entry */
+  int    ref;                /* Reference within forcefield */
+  char   ff_types[MAX_NO_MEMS][5];
+  double ff_param[MAX_NO_PARAMS];
+};
+
+struct FrcFieldItem {
+  char keyword[25];
+  int  number_of_members;    /* number of members of item */
+  int  number_of_parameters; /* number of parameters of item */
+  int  entries;              /* number of entries in item list */
+  struct FrcFieldData *data; /* contains all eqiuv and param data */
+};
+
+extern struct FrcFieldItem ff_atomtypes, equivalence, ff_vdw, ff_bond, ff_morse, ff_ang, ff_tor, ff_oop,
+  ff_bonbon, ff_bonang, ff_angtor, ff_angangtor, ff_endbontor, ff_midbontor, ff_angang, ff_bonbon13;
+
+/* prototypes */
+extern void InitializeItems(void);
+extern void SearchAndFill(struct FrcFieldItem *item);
+
+

tools/msi2lmp/src/InitializeItems.c

+/*
+*   This function fills in the keyword field, the number of members for each
+*   item and the number of parameters for each item
+*
+*/
+
+#include "msi2lmp.h"
+#include "Forcefield.h"
+
+#include <string.h>
+
+void InitializeItems(void)
+{
+  /* ATOM TYPES */
+  strcpy(ff_atomtypes.keyword,"#atom_types");
+  ff_atomtypes.number_of_members = 1;
+  ff_atomtypes.number_of_parameters = 1;
+
+  /* EQUIVALENCE */
+
+  strcpy(equivalence.keyword,"#equivalence");
+  equivalence.number_of_members = 6;
+  equivalence.number_of_parameters = 0;
+
+  /* NON-BOND */
+
+  strcpy(ff_vdw.keyword,"#nonbond");
+  ff_vdw.number_of_members = 1;
+  ff_vdw.number_of_parameters = 2;
+
+  /* BOND */
+
+  ff_bond.number_of_members = 2;
+  if (forcefield & FF_TYPE_CLASS1) {
+    strcpy(ff_bond.keyword,"#quadratic_bond");
+    ff_bond.number_of_parameters = 2;
+  }
+
+  if (forcefield & FF_TYPE_CLASS2) {
+    strcpy(ff_bond.keyword,"#quartic_bond");
+    ff_bond.number_of_parameters = 4;
+  }
+
+  /* MORSE */
+
+  if (forcefield & FF_TYPE_CLASS1) {
+    ff_morse.number_of_members = 2;
+    strcpy(ff_morse.keyword,"#morse_bond");
+    ff_morse.number_of_parameters = 3;
+  }
+
+  /* ANGLE */
+
+  ff_ang.number_of_members = 3;
+  if (forcefield & FF_TYPE_CLASS1) {
+    strcpy(ff_ang.keyword,"#quadratic_angle");
+    ff_ang.number_of_parameters = 2;
+  }
+
+  if (forcefield & FF_TYPE_CLASS2) {
+    strcpy(ff_ang.keyword,"#quartic_angle");
+    ff_ang.number_of_parameters = 4;
+  }
+
+  /* TORSION */
+
+  ff_tor.number_of_members = 4;
+  if (forcefield & FF_TYPE_CLASS1) {
+    strcpy(ff_tor.keyword,"#torsion_1");
+    ff_tor.number_of_parameters = 3;
+  } 
+
+  if (forcefield & FF_TYPE_CLASS2) {
+    strcpy(ff_tor.keyword,"#torsion_3");
+    ff_tor.number_of_parameters = 6;
+  }
+
+  /* OOP */
+
+  ff_oop.number_of_members = 4;
+  if (forcefield & FF_TYPE_CLASS1) {
+    strcpy(ff_oop.keyword,"#out_of_plane");
+    ff_oop.number_of_parameters = 3;
+  }
+
+  if (forcefield & FF_TYPE_CLASS2) {
+    strcpy(ff_oop.keyword,"#wilson_out_of_plane");
+    ff_oop.number_of_parameters = 2;
+  }
+
+  if (forcefield & FF_TYPE_CLASS2) {
+    /* BOND-BOND */
+
+    strcpy(ff_bonbon.keyword,"#bond-bond");
+    ff_bonbon.number_of_members = 3;
+    ff_bonbon.number_of_parameters = 1;
+
+    /* BOND-ANGLE */
+
+    strcpy(ff_bonang.keyword,"#bond-angle");
+    ff_bonang.number_of_members = 3;
+    ff_bonang.number_of_parameters = 2;
+
+    /* ANGLE-TORSION */
+
+    strcpy(ff_angtor.keyword,"#angle-torsion_3");
+    ff_angtor.number_of_members = 4;
+    ff_angtor.number_of_parameters = 6;
+
+    /* ANGLE-ANGLE-TORSION */
+
+    strcpy(ff_angangtor.keyword,"#angle-angle-torsion_1");
+    ff_angangtor.number_of_members = 4;
+    ff_angangtor.number_of_parameters = 1;
+
+    /* END-BOND-TORSION */
+
+    strcpy(ff_endbontor.keyword,"#end_bond-torsion_3");
+    ff_endbontor.number_of_members = 4;
+    ff_endbontor.number_of_parameters = 6;
+
+    /* MID-BOND-TORSION */
+
+    strcpy(ff_midbontor.keyword,"#middle_bond-torsion_3");
+    ff_midbontor.number_of_members = 4;
+    ff_midbontor.number_of_parameters = 3;
+
+    /* ANGLE-ANGLE */
+
+    strcpy(ff_angang.keyword,"#angle-angle");
+    ff_angang.number_of_members = 4;
+    ff_angang.number_of_parameters = 1;
+
+    /* BOND-BOND-1-3 */
+
+    strcpy(ff_bonbon13.keyword,"#bond-bond_1_3");
+    ff_bonbon13.number_of_members = 4;
+    ff_bonbon13.number_of_parameters = 1;
+  }
+}

tools/msi2lmp/src/Makefile

+
+TARGET	 = msi2lmp.exe
+
+SRCS     = msi2lmp.c \
+           ReadCarFile.c \
+           ReadMdfFile.c \
+	   MakeLists.c \
+           ReadFrcFile.c \
+           InitializeItems.c \
+	   SearchAndFill.c \
+           GetParameters.c \
+	   CheckLists.c \
+           WriteDataFile.c
+
+OBJS     = $(SRCS:.c=.o)
+
+HEADERS  = msi2lmp.h Forcefield.h
+
+CC       = gcc
+CFLAGS   = -O -Wall -W -g 
+FRCFILE  = cvff.frc
+FRCFILE2 = cff91.frc
+README   = README
+MKFILE   = Makefile
+
+$(TARGET) : $(OBJS)
+	$(CC) $(CFLAGS) -o $(TARGET) $(OBJS) -lm
+
+.c.o:
+	$(CC) $(CFLAGS) -c $<
+clean:
+	rm -f $(OBJS) $(TARGET)
+
+# dependencies
+CheckLists.o: CheckLists.c msi2lmp.h
+GetParameters.o: GetParameters.c msi2lmp.h Forcefield.h
+InitializeItems.o: InitializeItems.c msi2lmp.h Forcefield.h
+MakeLists.o: MakeLists.c msi2lmp.h
+msi2lmp.o: msi2lmp.c msi2lmp.h
+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
+WriteDataFile.o: WriteDataFile.c msi2lmp.h Forcefield.h

tools/msi2lmp/src/ReadCarFile.c

+/*
+*  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][0];
+  h[4] = box[2][1];
+  h[5] = box[2][2];
+  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 %f",
+             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 (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() */