path: root/sci-chemistry/ambertools/files/ambertools-1.5-bugfix_1-21.patch

diff -urN amber11.orig/AmberTools/src/configure amber11/AmberTools/src/configure
--- amber11.orig/AmberTools/src/configure	2011-04-14 15:30:55.000000000 +0300
+++ amber11/AmberTools/src/configure	2011-10-25 15:01:28.222288537 +0300
@@ -650,7 +650,7 @@
     if [ "$cuda" = 'yes' -o "$cuda_SPSP" = 'yes' -o "$cuda_DPDP" = 'yes' ]; then
       pmemd_cu_includes='-I$(CUDA_HOME)/include -IB40C -IB40C/KernelCommon'
       pmemd_cu_defines='-DCUDA'
-      pmemd_cu_libs='-L$(CUDA_HOME)/lib64 -L$(CUDA_HOME)/lib -lcufft -lcudart ./cuda/cuda.a'
+      pmemd_cu_libs='-L$(CUDA_HOME)/lib64 -L$(CUDA_HOME)/lib -lcurand -lcufft -lcudart ./cuda/cuda.a'
       if [ "$optimise" = 'no' ]; then
         nvcc='$(CUDA_HOME)/bin/nvcc -use_fast_math -O0 -gencode arch=compute_13,code=sm_13 -gencode arch=compute_20,code=sm_20'
       else
@@ -909,7 +909,7 @@
 
         pmemd_cu_includes='-I$(CUDA_HOME)/include -IB40C -IB40C/KernelCommon'
         pmemd_cu_defines='-DCUDA'
-        pmemd_cu_libs='-L$(CUDA_HOME)/lib64 -L$(CUDA_HOME)/lib -lcufft -lcudart ./cuda/cuda.a'
+        pmemd_cu_libs='-L$(CUDA_HOME)/lib64 -L$(CUDA_HOME)/lib -lcurand -lcufft -lcudart ./cuda/cuda.a'
         if [ "$optimise" = 'yes' ]; then
             nvcc='$(CUDA_HOME)/bin/nvcc -use_fast_math -O3 -gencode arch=compute_13,code=sm_13 -gencode arch=compute_20,code=sm_20'
         else
@@ -2111,21 +2111,21 @@
     echo "Change to \$AMBERHOME/src/ and type 'make clean' followed by"
     echo "'make parallel_win'."
   elif [ "$cuda" = 'yes' -o "$cuda_SPSP" = 'yes' -o "$cuda_DPDP" = 'yes' ]; then
-    echo "The next step is to cd to ../../src and type 'make install'"
+    echo "The next step is to cd to ../../src and type 'make $installtype'"
   else
     echo "If you have amber11, the most common next step is to "
-    echo "'cd ../../src; make clean; make install'.  See the Users' Manual"
+    echo "'cd ../../src; make clean; make $installtype'.  See the Users' Manual"
     echo "for information on building a parallel version of AmberTools"
   fi
 elif [ "$cuda" = 'yes' -o "$cuda_SPSP" = 'yes' -o "$cuda_DPDP" = 'yes' ]; then
-  echo "The next step is to cd to ../../src and type 'make install'"
+  echo "The next step is to cd to ../../src and type 'make $installtype'"
 elif [ "$windows" = 'yes' ]; then
   echo "The next step is to switch to a dos shell with ifort enabled."
   echo "Change to \$AMBERHOME/src/ and type 'make serial_win'."
 elif [ "$openmp" = 'yes' ]; then
   echo "The next step is to type 'make nabonly'"
 else
-  echo "The next step is to type 'make install'"
+  echo "The next step is to type 'make $installtype'"
 fi
 echo " "
 
@@ -2149,7 +2149,7 @@
    echo " NOTE: after installing AmberTools, if you want to compile"
    echo "       Amber11 using AmberTools 1.5, you must run the "
    echo "       AT15_Amber11.py script first, e.g.:"
-   echo "       cd $AMBERHOME; ./AT15_Amber11.py; make $installtype"
+   echo "       cd $AMBERHOME; ./AT15_Amber11.py; cd src; make $installtype"
    echo "       (See the Amber11 Users' Manual for full instructions,"
    echo "       and be sure to have a backup, since the AT15_Amber11.py"
    echo "       script will modify your Amber11 files.)"
diff -urN amber11.orig/AmberTools/src/cpptraj/src/Action_Center.cpp amber11/AmberTools/src/cpptraj/src/Action_Center.cpp
--- amber11.orig/AmberTools/src/cpptraj/src/Action_Center.cpp	2011-04-14 15:30:16.000000000 +0300
+++ amber11/AmberTools/src/cpptraj/src/Action_Center.cpp	2011-10-25 15:01:28.108288642 +0300
@@ -72,7 +72,7 @@
 // else
 //    mass = NULL;
 
-  if (!origin && P->ifbox==0) {
+  if (!origin && P->BoxType==0) {
     mprintf("    Error: Center::setup: Box center specified but no box information.\n");
     //fprintf(stdout,"                            Centering on origin.\n");
     return 1;
diff -urN amber11.orig/AmberTools/src/cpptraj/src/Action_Closest.cpp amber11/AmberTools/src/cpptraj/src/Action_Closest.cpp
--- amber11.orig/AmberTools/src/cpptraj/src/Action_Closest.cpp	2011-04-14 15:30:16.000000000 +0300
+++ amber11/AmberTools/src/cpptraj/src/Action_Closest.cpp	2011-10-25 15:01:28.108288642 +0300
@@ -181,8 +181,8 @@
   // NOTE: Should box be figured out from read-in coords?
   imageType = 0;
   if (!noimage) {
-    imageType = P->ifbox;
-    if (P->ifbox==0) {
+    imageType = P->BoxType;
+    if (P->BoxType==0) {
         mprintf("    Warning: Closest::setup: ");
         mprintf(" Imaging specified but no box information in prmtop %s\n",P->parmName);
         mprintf("             No imaging can occur..\n");
diff -urN amber11.orig/AmberTools/src/cpptraj/src/Action_Distance.cpp amber11/AmberTools/src/cpptraj/src/Action_Distance.cpp
--- amber11.orig/AmberTools/src/cpptraj/src/Action_Distance.cpp	2011-04-14 15:30:16.000000000 +0300
+++ amber11/AmberTools/src/cpptraj/src/Action_Distance.cpp	2011-10-25 15:01:28.194288563 +0300
@@ -57,8 +57,6 @@
   mprintf("    DISTANCE: %s to %s",Mask1.maskString, Mask2.maskString);
   if (noimage) 
     mprintf(", non-imaged");
-  else
-    mprintf(", imaged");
   if (useMass) 
     mprintf(", center of mass");
   else
@@ -68,6 +66,11 @@
   return 0;
 }
 
+/*
+ * Distance::setup()
+ * Determine what atoms each mask pertains to for the current parm file.
+ * Also determine whether imaging should be performed.
+ */
 int Distance::setup() {
 
   if ( Mask1.SetupMask(P,debug) ) return 1;
@@ -89,11 +92,19 @@
   // Check imaging - check box based on prmtop box
   imageType = 0;
   if (!noimage) {
-    imageType = P->ifbox;
-    if (P->ifbox==0 && debug>0) {
+    imageType = P->BoxType;
+    if (P->BoxType==0 && debug>0) {
       mprintf("    Warning: No box info in %s, disabling imaging.\n",P->parmName);
     }
   }
+
+  // Print imaging info for this parm
+  mprintf("    DISTANCE: %s to %s",Mask1.maskString, Mask2.maskString);
+  if (imageType > 0)
+    mprintf(", imaged");
+  else
+    mprintf(", imaging off");
+  mprintf(".\n");
         
   return 0;  
 }
diff -urN amber11.orig/AmberTools/src/cpptraj/src/Action_DSSP.cpp amber11/AmberTools/src/cpptraj/src/Action_DSSP.cpp
--- amber11.orig/AmberTools/src/cpptraj/src/Action_DSSP.cpp	2011-04-14 15:30:16.000000000 +0300
+++ amber11/AmberTools/src/cpptraj/src/Action_DSSP.cpp	2011-10-25 15:01:28.195288561 +0300
@@ -76,7 +76,7 @@
   int selected, atom, res;
   Residue RES;
 
-  // Set up mask
+  // Set up mask for this parm
   if ( Mask.SetupMask(P,debug) ) return 1;
   if ( Mask.None() ) {
     mprintf("      Error: DSSP::setup: Mask has no atoms.\n");
@@ -90,25 +90,25 @@
     Nres=P->nres;
   //mprintf("      DSSP: Setting up for %i residues.\n",Nres);
 
-  // Free up SecStruct if previously allocated.
-  // NOTE: In action setup, should check if Parm has really changed...
-  SecStruct.clear();
+  // Set up for each residue of the current Parm if not already set-up.
   for (res = 0; res < Nres; res++) {
-    RES.sstype=SECSTRUCT_NULL;
-    RES.isSelected=false;
-    RES.C=-1;
-    RES.O=-1;
-    RES.N=-1;
-    RES.H=-1;
-    RES.CO_HN_Hbond.assign( Nres, 0 );
-    RES.SSprob[0]=0.0;
-    RES.SSprob[1]=0.0;
-    RES.SSprob[2]=0.0;
-    RES.SSprob[3]=0.0;
-    RES.SSprob[4]=0.0;
-    RES.SSprob[5]=0.0;
-    RES.SSprob[6]=0.0;
-    SecStruct.push_back(RES);
+    if (res>=(int)SecStruct.size()) {
+      RES.sstype=SECSTRUCT_NULL;
+      RES.isSelected=false;
+      RES.C=-1;
+      RES.O=-1;
+      RES.N=-1;
+      RES.H=-1;
+      RES.CO_HN_Hbond.assign( Nres, 0 );
+      RES.SSprob[0]=0.0;
+      RES.SSprob[1]=0.0;
+      RES.SSprob[2]=0.0;
+      RES.SSprob[3]=0.0;
+      RES.SSprob[4]=0.0;
+      RES.SSprob[5]=0.0;
+      RES.SSprob[6]=0.0;
+      SecStruct.push_back(RES);
+    }
   }
 
   // Go through all atoms in mask. Set up a residue for each C, O, N, and H atom
diff -urN amber11.orig/AmberTools/src/cpptraj/src/Action_Image.cpp amber11/AmberTools/src/cpptraj/src/Action_Image.cpp
--- amber11.orig/AmberTools/src/cpptraj/src/Action_Image.cpp	2011-04-14 15:30:16.000000000 +0300
+++ amber11/AmberTools/src/cpptraj/src/Action_Image.cpp	2011-10-25 15:01:28.110288640 +0300
@@ -98,14 +98,14 @@
     useMass = false;
   }
 
-  if (P->ifbox==0) {
+  if (P->BoxType==0) {
     mprintf("    Error: Image::setup: Parm %s does not contain box information.\n",
             P->parmName);
     return 1;
   }
 
   ortho = false;  
-  if (P->ifbox==1 && triclinic==OFF) ortho=true;
+  if (P->BoxType==1 && triclinic==OFF) ortho=true;
 
   if (triclinic == FAMILIAR) {
     if (ComMask!=NULL) {
diff -urN amber11.orig/AmberTools/src/cpptraj/src/Action_Outtraj.cpp amber11/AmberTools/src/cpptraj/src/Action_Outtraj.cpp
--- amber11.orig/AmberTools/src/cpptraj/src/Action_Outtraj.cpp	2011-04-14 15:30:16.000000000 +0300
+++ amber11/AmberTools/src/cpptraj/src/Action_Outtraj.cpp	2011-10-25 15:01:28.129288623 +0300
@@ -6,6 +6,9 @@
 // CONSTRUCTOR
 Outtraj::Outtraj() {
   //fprintf(stderr,"Outtraj Con\n");
+  min=0.0;
+  max=0.0;
+  Dset=NULL;
 } 
 
 // DESTRUCTOR
@@ -13,11 +16,40 @@
 
 /*
  * Outtraj::init()
- * Action wrapper for trajout
+ * Action wrapper for trajout.
+ * Expected call: outtraj <filename> [ trajout args ] 
+ *                        [maxmin <dataset> min <min> max <max>
  */
 int Outtraj::init() {
+  char *datasetName;
+
+#ifdef MPI
+  mprintf("ERROR: OUTTRAJ currently not functional with MPI.\n");
+  return 1;
+#endif
+
+  mprintf("    OUTTRAJ: Will write to [%s]\n",A->Arg(1));
+  outtraj.SetDebug(debug);
+  // If maxmin, get the name of the dataset as well as the max and min values.
+  datasetName = A->getKeyString("maxmin",NULL);
+  if (datasetName!=NULL) {
+    Dset = DSL->Get(datasetName);
+    if (Dset==NULL) {
+      mprintf("Error: Outtraj maxmin: Could not get dataset %s\n",datasetName);
+      return 1;
+    } else {
+      // Currently dont allow for string datasets
+      if (Dset->Type()==STRING) {
+        mprintf("Error: Outtraj maxmin: String dataset (%s) not supported.\n",datasetName);
+        return 1;
+      }
+      max = A->getKeyDouble("max",0.0);
+      min = A->getKeyDouble("min",0.0);
+      mprintf("             maxmin: Printing trajectory frames based on %lf <= %s <= %lf\n",
+              min, datasetName, max);
+    }
+  }
 
-  mprintf("    OUTTRAJ: [%s]\n",A->ArgLine());
   return ( outtraj.Add(A,PFL,worldsize) );
 } 
 
@@ -33,8 +65,32 @@
  * Outtraj::action()
  */
 int Outtraj::action() {
+  double dVal;
+  int iVal;
 
-  return ( outtraj.Write(currentFrame, F, P) );
+  // If dataset defined, check if frame is within max/min
+  if (Dset!=NULL) {
+    if (Dset->Type() == DOUBLE) {
+      if (Dset->Get(&dVal, currentFrame)) return 1;
+    } else if (Dset->Type() == INT) {
+      if (Dset->Get(&iVal, currentFrame)) return 1;
+      dVal = (double) iVal;
+    } else
+      return 1;
+    //mprintf("DBG: maxmin: dVal = %lf\n",dVal);
+    // If value from dataset not within min/max, exit now.
+    if (dVal < min || dVal > max) return 0;
+  }
+  if ( outtraj.Write(currentFrame, F, P) != 0 ) return 1;
+  return 0;
 } 
 
+/*
+ * Outtraj::print()
+ * Close trajectory.
+ */
+void Outtraj::print() {
+  mprintf("  OUTTRAJ: [%s] Wrote %i frames.\n",A->Arg(1),outtraj.front()->CurrentFrame());
+  outtraj.Close();
+}
 
diff -urN amber11.orig/AmberTools/src/cpptraj/src/Action_Outtraj.h amber11/AmberTools/src/cpptraj/src/Action_Outtraj.h
--- amber11.orig/AmberTools/src/cpptraj/src/Action_Outtraj.h	2011-04-14 15:30:16.000000000 +0300
+++ amber11/AmberTools/src/cpptraj/src/Action_Outtraj.h	2011-10-25 15:01:28.129288623 +0300
@@ -6,6 +6,9 @@
 
 class Outtraj: public Action {
     TrajoutList outtraj;
+    double max;
+    double min;
+    DataSet *Dset;
   public:
     Outtraj();
     ~Outtraj();
@@ -13,5 +16,6 @@
     int init();
     //int setup();
     int action();
+    void print();
 };
 #endif
diff -urN amber11.orig/AmberTools/src/cpptraj/src/Action_Rms2d.cpp amber11/AmberTools/src/cpptraj/src/Action_Rms2d.cpp
--- amber11.orig/AmberTools/src/cpptraj/src/Action_Rms2d.cpp	1970-01-01 03:00:00.000000000 +0300
+++ amber11/AmberTools/src/cpptraj/src/Action_Rms2d.cpp	2011-10-25 15:01:28.074288674 +0300
@@ -0,0 +1,71 @@
+// Rms2d 
+#include "Action_Rms2d.h"
+#include "CpptrajStdio.h"
+
+// CONSTRUCTOR
+Rms2d::Rms2d() {
+  //fprintf(stderr,"Rms2d Con\n");
+  nofit=false;
+  useMass=false;
+} 
+
+// DESTRUCTOR
+Rms2d::~Rms2d() { }
+
+/*
+ * Rms2d::init()
+ * Expected call: rms2d <mask> <refmask> [rmsout filename] [mass] [nofit]
+ * Dataset name will be the last arg checked for. Check order is:
+ *    1) Keywords
+ *    2) Masks
+ *    3) Dataset name
+ */
+int Rms2d::init() {
+  char *mask0, *maskRef;
+
+  // Get keywords
+  nofit = A->hasKey("nofit");
+  useMass = A->hasKey("mass");
+  rmsdFile = A->getKeyString("rmsout",NULL);
+
+  // Get the RMS mask string for frames
+  mask0 = A->getNextMask();
+  FrameMask.SetMaskString(mask0);
+  // Get RMS mask string for reference
+  maskRef = A->getNextMask();
+  // If no reference mask specified, make same as RMS mask
+  if (maskRef==NULL) maskRef=mask0;
+  RefMask.SetMaskString(maskRef);
+
+  mprintf("    RMS2D: (%s) to (%s)",FrameMask.maskString,RefMask.maskString);
+  if (nofit)
+    mprintf(" (no fitting)");
+  if (useMass)
+    mprintf(" (mass-weighted)");
+  if (rmsdFile!=NULL) 
+    mprintf(" output to %s",rmsdFile);
+  mprintf("\n");
+
+  return 0;
+}
+
+/*
+ * Rms2d::setup()
+ * Not important for Rms2d, initial pass is only for storing frames.
+ */
+int Rms2d::setup() {
+  return 0;  
+}
+
+/*
+ * Rms2d::action()
+ * Store current frame as a reference frame.
+ */
+int Rms2d::action() {
+
+  ReferenceFrames.Add(F,P);
+  
+  return 0;
+} 
+
+
diff -urN amber11.orig/AmberTools/src/cpptraj/src/Action_Rms2d.h amber11/AmberTools/src/cpptraj/src/Action_Rms2d.h
--- amber11.orig/AmberTools/src/cpptraj/src/Action_Rms2d.h	1970-01-01 03:00:00.000000000 +0300
+++ amber11/AmberTools/src/cpptraj/src/Action_Rms2d.h	2011-10-25 15:01:28.075288673 +0300
@@ -0,0 +1,21 @@
+#ifndef INC_ACTION_RMS2D_H
+#define INC_ACTION_RMS2D_H
+// Rms2d
+#include "Action.h"
+
+class Rms2d: public Action {
+    FrameList ReferenceFrames;
+    bool nofit;
+    bool useMass;
+    AtomMask RefMask;
+    AtomMask FrameMask;
+    char *rmsdFile;
+  public:
+    Rms2d();
+    ~Rms2d();
+
+    int init();
+    int setup();
+    int action();
+};
+#endif
diff -urN amber11.orig/AmberTools/src/cpptraj/src/Action_Rmsd.cpp amber11/AmberTools/src/cpptraj/src/Action_Rmsd.cpp
--- amber11.orig/AmberTools/src/cpptraj/src/Action_Rmsd.cpp	2011-04-14 15:30:16.000000000 +0300
+++ amber11/AmberTools/src/cpptraj/src/Action_Rmsd.cpp	2011-10-25 15:01:28.229288531 +0300
@@ -256,6 +256,7 @@
     PerResRMSD=new DataSetList();
     resName[4]='\0';
     for (it=ResRange.begin(); it!=ResRange.end(); it++) {
+      //res = *it - 1; // res is the internal resnumber, *it the user resnumber
       // Get corresponding reference resnum - if none specified use current res
       if (RefRange.empty()) 
         refRes = (*it);
@@ -263,7 +264,18 @@
         refRes = RefRange.front();
         RefRange.pop_front();
       }
-      //res = *it - 1; // res is the internal resnumber, *it the user resnumber
+
+      // Check if either the residue num or the reference residue num out of range.
+      if ( (*it) < 1 || (*it) > P->nres) {
+        mprintf("    Warning: Rmsd: perres: Specified residue # %i is out of range.\n",*it);
+        continue;
+      }
+      if ( refRes < 1 || refRes > P->nres ) {
+        mprintf("    Warning: Rmsd: perres: Specified reference residue # %i is out of range.\n",
+                refRes);
+        continue;
+      }
+
       // Setup Dataset Name to be name of this residue 
       P->ResName(resName,(*it)-1);
       
diff -urN amber11.orig/AmberTools/src/cpptraj/src/AmberNetcdf.cpp amber11/AmberTools/src/cpptraj/src/AmberNetcdf.cpp
--- amber11.orig/AmberTools/src/cpptraj/src/AmberNetcdf.cpp	2011-04-14 15:30:16.000000000 +0300
+++ amber11/AmberTools/src/cpptraj/src/AmberNetcdf.cpp	2011-10-25 15:01:28.131288621 +0300
@@ -49,6 +49,7 @@
  * and close calls.
  */
 void AmberNetcdf::close() {
+  if (ncid<0) return;
   checkNCerr(nc_close(ncid),"Closing netcdf file.");
   if (debug>0) rprintf("Successfully closed ncid %i\n",ncid);
   ncid=-1;
@@ -92,8 +93,10 @@
  * Open the netcdf file, read all dimension and variable IDs, close. 
  */
 int AmberNetcdf::SetupRead() {
-  char *attrText; // For checking conventions and version 
-  int spatial; // For checking spatial dimensions
+  char *attrText;            // For checking conventions and version 
+  int spatial;               // For checking spatial dimensions
+  double box[6];             // For checking box type
+  size_t start[3], count[3]; // For checking box type
 
   if (open()) return 1;
 
@@ -147,15 +150,15 @@
   if ( nc_inq_varid(ncid,"cell_lengths",&cellLengthVID)==NC_NOERR ) {
     if (checkNCerr(nc_inq_varid(ncid,"cell_angles",&cellAngleVID),
       "Getting cell angles.")!=0) return 1;
-    if (debug>0) mprintf("  Netcdf Box information found.\n"); 
-    if (P->ifbox==0) {
-      mprintf("Warning: Netcdf file contains box info but no box info found\n");
-      mprintf("         in associated parmfile %s; defaulting to orthogonal.\n",
-              P->parmName);
-      isBox=1;
-    } else {
-      isBox=P->ifbox;
-    }
+    if (debug>0) mprintf("  Netcdf Box information found.\n");
+    // Determine box type from angles
+    start[0]=0; start[1]=0; start[2]=0; 
+    count[0]=1; count[1]=3; count[2]=0;
+    if ( checkNCerr(nc_get_vara_double(ncid, cellLengthVID, start, count, box),
+                    "Getting cell lengths.")!=0 ) return 1;
+    if ( checkNCerr(nc_get_vara_double(ncid, cellAngleVID, start, count, box+3),
+                    "Getting cell angles.")!=0 ) return 1;
+    CheckBoxType(box);
   } 
 
   // Replica Temperatures
@@ -258,7 +261,7 @@
     "Defining cell angular variable.")) return 1;
 
   // Box Info
-  if (isBox>0) {
+  if (BoxType!=0) {
     dimensionID[0]=frameDID;
     dimensionID[1]=cell_spatialDID;
     if (checkNCerr(nc_def_var(ncid,"cell_lengths",NC_DOUBLE,2,dimensionID,&cellLengthVID),
@@ -362,7 +365,7 @@
                   "Getting frame %i",set)!=0 ) return 1;
 
   // Read box info 
-  if (isBox!=0) {
+  if (BoxType!=0) {
     count [1]=3;
     count [2]=0;
     if ( checkNCerr(nc_get_vara_double(ncid, cellLengthVID, start, count, F->box),
@@ -399,7 +402,7 @@
   F->frameToFloat(Coord);
 
   // write coords
-  start[0]=set;
+  start[0]=currentFrame;
   start[1]=0;
   start[2]=0;
   count[0]=1;
@@ -409,7 +412,7 @@
     "Netcdf Writing frame %i",set)) return 1;
 
   // write box
-  if (isBox>0 && cellLengthVID!=-1) {
+  if (BoxType!=0 && cellLengthVID!=-1) {
     count[1]=3;
     count[2]=0;
     if (checkNCerr(nc_put_vara_double(ncid,cellLengthVID,start,count,F->box),
@@ -427,6 +430,8 @@
   
   nc_sync(ncid); // Necessary after every write??
 
+  currentFrame++;
+
   return 0;
 }  
 
diff -urN amber11.orig/AmberTools/src/cpptraj/src/AmberParm.cpp amber11/AmberTools/src/cpptraj/src/AmberParm.cpp
--- amber11.orig/AmberTools/src/cpptraj/src/AmberParm.cpp	2011-04-14 15:30:16.000000000 +0300
+++ amber11/AmberTools/src/cpptraj/src/AmberParm.cpp	2011-10-25 15:01:28.231288529 +0300
@@ -1,6 +1,12 @@
 /* AmberParm.cpp
- * Class that holds parameter information. Can be read in from Amber Topology
- * or PDB files.
+ * Class that holds parameter information. Can be read in from Amber Topology,
+ * PDB, or Mol2 files (implemented in the ReadParmXXX functions). The following
+ * parameters of AmberParm must always be set:
+ *   The NATOM, NRES, and IFBOX entries of the values array.
+ *   The names, resnames, resnums arrays.
+ *   The natom and nres variables.
+ * NOTES:
+ *   Eventually make the mol2 read parm function use the AddBond function.
  */
 #include <cstdlib>
 #include <cstring>
@@ -11,7 +17,7 @@
 #include "CpptrajStdio.h"
 
 #define AMBERPOINTERS 31
-
+#define TRUNCOCTBETA 109.4712206344906917365733534097672
 #define ELECTOAMBER 18.2223
 #define AMBERTOELEC 1/ELECTOAMBER
 // =============================================================
@@ -21,7 +27,7 @@
  * Given a residue number, set buffer with residue name. Replace blanks with _
  */
 void AmberParm::ResName(char *buffer, int res) {
-  if (res<0 || res>nres) return;
+  if (res<0 || res>=nres) return;
   strcpy(buffer, resnames[res]);
   if (buffer[3]==' ') buffer[3]='_';
 }
@@ -203,7 +209,9 @@
   bondsh=NULL;   
   types=NULL;   
   atomsPerMol=NULL; 
-  Box=NULL;    
+  Box[0]=0.0; Box[1]=0.0; Box[2]=0.0;
+  Box[3]=0.0; Box[4]=0.0; Box[5]=0.0;
+  BoxType=0; 
   pindex=0;      
   parmFrames=0; 
   outFrame=0;
@@ -217,7 +225,6 @@
   solventMoleculeStop=NULL;
   natom=0;      
   nres=0; 
-  ifbox=0; 
   parmName=NULL; 
   SurfaceInfo=NULL;
 }
@@ -234,7 +241,6 @@
   if (bonds!=NULL) free(bonds);
   if (bondsh!=NULL) free(bondsh);
   if (atomsPerMol!=NULL) free(atomsPerMol);
-  if (Box!=NULL) free(Box);
   if (solventMoleculeStart!=NULL) free(solventMoleculeStart);
   if (solventMoleculeStop!=NULL) free(solventMoleculeStop);
   if (solventMask!=NULL) free(solventMask);
@@ -271,54 +277,7 @@
   return NULL;
 }
 
-/* 
- * AmberParm::OpenParm()
- * Attempt to open file and read in parameters.
- */
-int AmberParm::OpenParm(char *filename) {
-  if ( File.SetupFile(filename,READ,UNKNOWN_FORMAT, UNKNOWN_TYPE,debug) ) return 1;
-
-  // Copy parm filename to parmName. Separate from File.filename in case of stripped parm
-  parmName=(char*) malloc( (strlen(File.basefilename)+1) * sizeof(char));
-  strcpy(parmName,File.basefilename);
-
-  if ( File.OpenFile() ) return 1;
-
-  switch (File.fileFormat) {
-    case AMBERPARM : if (ReadParmAmber()) return 1; break;
-    case PDBFILE   : if (ReadParmPDB()  ) return 1; break;
-    case MOL2FILE  : if (ReadParmMol2() ) return 1; break;
-    default: 
-      rprintf("Unknown parameter file type: %s\n",File.filename);
-      return 1;
-  }
-
-  File.CloseFile();
-
-  // Create a last dummy residue in resnums that holds natom, which would be
-  // the atom number of the next residue if it existed. Atom #s in resnums
-  // should correspond with cpptraj atom #s (start from 0) instead of Amber
-  // atom #s (start from 1). 
-  // Do this to be consistent with ptrajmask selection behavior - saves an 
-  // if-then statement.
-  resnums=(int*) realloc(resnums,(nres+1)*sizeof(int));
-  resnums[nres]=natom;
-  // DEBUG
-  //fprintf(stdout,"==== DEBUG ==== Resnums for %s:\n",File.filename);
-  //for (err=0; err<nres; err++) 
-  //  fprintf(stdout,"    %i: %i\n",err,resnums[err]);
-
-  // Set up solvent information
-  SetSolventInfo();
-
-  if (debug>0) {
-    mprintf("  Number of atoms= %i\n",natom);
-    mprintf("  Number of residues= %i\n",nres);
-  }
-
-  return 0;
-}
-
+// ---------========= ROUTINES PERTAINING TO SURFACE AREA =========---------
 /*
  * AmberParm::AssignLCPO()
  * Assign parameters for LCPO method. All radii are incremented by 1.4 Ang.
@@ -478,6 +437,7 @@
   return 0;
 }
 
+// ---------========= ROUTINES PERTAINING TO SOLVENT INFO =========---------
 /*
  * AmberParm::SetSolventInfo()
  * Assuming atomsPerMol has been read in, set solvent information.
@@ -536,7 +496,56 @@
 
   return 0; 
 }
-     
+    
+// --------========= ROUTINES PERTAINING TO READING PARAMETERS =========--------
+/* 
+ * AmberParm::OpenParm()
+ * Attempt to open file and read in parameters.
+ */
+int AmberParm::OpenParm(char *filename) {
+  if ( File.SetupFile(filename,READ,UNKNOWN_FORMAT, UNKNOWN_TYPE,debug) ) return 1;
+
+  // Copy parm filename to parmName. Separate from File.filename in case of stripped parm
+  parmName=(char*) malloc( (strlen(File.basefilename)+1) * sizeof(char));
+  strcpy(parmName,File.basefilename);
+
+  if ( File.OpenFile() ) return 1;
+
+  switch (File.fileFormat) {
+    case AMBERPARM : if (ReadParmAmber()) return 1; break;
+    case PDBFILE   : if (ReadParmPDB()  ) return 1; break;
+    case MOL2FILE  : if (ReadParmMol2() ) return 1; break;
+    default: 
+      rprintf("Unknown parameter file type: %s\n",File.filename);
+      return 1;
+  }
+
+  File.CloseFile();
+
+  // Create a last dummy residue in resnums that holds natom, which would be
+  // the atom number of the next residue if it existed. Atom #s in resnums
+  // should correspond with cpptraj atom #s (start from 0) instead of Amber
+  // atom #s (start from 1). 
+  // Do this to be consistent with ptrajmask selection behavior - saves an 
+  // if-then statement.
+  resnums=(int*) realloc(resnums,(nres+1)*sizeof(int));
+  resnums[nres]=natom;
+  // DEBUG
+  //fprintf(stdout,"==== DEBUG ==== Resnums for %s:\n",File.filename);
+  //for (err=0; err<nres; err++) 
+  //  fprintf(stdout,"    %i: %i\n",err,resnums[err]);
+
+  // Set up solvent information
+  SetSolventInfo();
+
+  if (debug>0) {
+    mprintf("  Number of atoms= %i\n",natom);
+    mprintf("  Number of residues= %i\n",nres);
+  }
+
+  return 0;
+}
+
 /* 
  * AmberParm::ReadParmAmber() 
  * Read parameters from Amber Topology file
@@ -544,6 +553,7 @@
 int AmberParm::ReadParmAmber() {
   int err, atom;
   int *solvent_pointer;
+  double *boxFromParm;
 
   if (debug>0) mprintf("Reading Amber Topology file %s\n",parmName);
 
@@ -555,7 +565,6 @@
 
   natom=values[NATOM];
   nres=values[NRES];
-  ifbox=values[IFBOX];
   if (debug>0)
     mprintf("    Amber top contains %i atoms, %i residues.\n",natom,nres);
 
@@ -583,7 +592,8 @@
   if (bonds==NULL) {mprintf("Error in bonds w/o H.\n"); err++;}
   bondsh=(int*) getFlagFileValues("BONDS_INC_HYDROGEN",values[NBONH]*3);
   if (bondsh==NULL) {mprintf("Error in bonds inc H.\n"); err++;}
-  if (ifbox>0) {
+  // Get solvent info if IFBOX>0
+  if (values[IFBOX]>0) {
     solvent_pointer=(int*) getFlagFileValues("SOLVENT_POINTERS",3);
     if (solvent_pointer==NULL) {
       mprintf("Error in solvent pointers.\n"); 
@@ -596,12 +606,20 @@
     }
     atomsPerMol=(int*) getFlagFileValues("ATOMS_PER_MOLECULE",molecules);
     if (atomsPerMol==NULL) {mprintf("Error in atoms per molecule.\n"); err++;}
-    Box=(double*) getFlagFileValues("BOX_DIMENSIONS",4);
-    if (Box==NULL) {mprintf("Error in Box information.\n"); err++;}
+    // boxFromParm = {OLDBETA, BOX(1), BOX(2), BOX(3)}
+    boxFromParm=(double*) getFlagFileValues("BOX_DIMENSIONS",4);
+    if (boxFromParm==NULL) {mprintf("Error in Box information.\n"); err++;}
+    // Determine box type: 1-Ortho, 2-Nonortho
+    SetBoxInfo(boxFromParm[0],boxFromParm[1],boxFromParm[2],boxFromParm[3]);
+    free(boxFromParm);
     if (debug>0) {
       mprintf("    %s contains box info: %i mols, first solvent mol is %i\n",
               parmName, molecules, firstSolvMol);
-      mprintf("    BOX: %lf %lf %lf %lf\n",Box[0],Box[1],Box[2],Box[3]);
+      mprintf("    BOX: %lf %lf %lf | %lf %lf %lf\n",Box[0],Box[1],Box[2],Box[3],Box[4],Box[5]);
+      if (BoxType==1)
+        mprintf("         Box is orthogonal.\n");
+      else
+        mprintf("         Box is non-orthogonal.\n");
     }
   }
 
@@ -686,6 +704,7 @@
   values = (int*) calloc(AMBERPOINTERS, sizeof(int));
   values[NATOM] = natom;
   values[NRES] = nres;
+  values[IFBOX] = 0;
 
   if (debug>0) 
     mprintf("    PDB contains %i atoms, %i residues, %i molecules.\n",
@@ -799,6 +818,7 @@
   values[NRES] = nres;
   values[NBONH] = numbondsh;
   values[MBONA] = numbonds;
+  values[IFBOX] = 0;
 
   mprintf("    Mol2 contains %i atoms, %i residues,\n", natom,nres);
   mprintf("    %i bonds to H, %i other bonds.\n", numbondsh,numbonds);
@@ -806,6 +826,61 @@
   return 0;
 }
 
+// ---------===========================================================---------
+/*
+ * AmberParm::SetBoxInfo()
+ * Given 3 box lengths and an angle determine the box type and set
+ * the box information. If called with negative beta, set no box.
+ * Currently recognized betas:
+ *   90.00 - Orthogonal
+ *  109.47 - Truncated octahedral
+ *   60.00 - Rhombic dodecahedron
+ * Any other beta just sets all angles to beta and a warning is printed.
+ */
+int AmberParm::SetBoxInfo(double beta, double bx, double by, double bz)  {
+  int ifbox=0;
+
+  // Determine box type from beta (none, ortho, non-ortho (truncated oct/triclinic)
+  if (beta<=0.0) {
+    if (BoxType>0) 
+      mprintf("    %s: Removing box information.\n",parmName);
+    BoxType=0;
+    ifbox=0;
+    Box[0]=0.0; Box[1]=0.0; Box[2]=0.0;
+    Box[3]=0.0; Box[4]=0.0; Box[5]=0.0;
+  } else if (beta == 90.0) {
+    BoxType=1;
+    ifbox=1;
+    Box[0]=bx; Box[1]=by; Box[2]=bz;
+    Box[3]=90.0; Box[4]=90.0; Box[5]=90.0;
+    if (debug>0) mprintf("    %s: Setting box to be orthogonal.\n",parmName);
+  } else if (beta > 109.47 && beta < 109.48) {
+    BoxType=2;
+    ifbox=2;
+    Box[0]=bx; Box[1]=by; Box[2]=bz;
+    //Box[3] = TRUNCOCTBETA;
+    Box[3] = beta;
+    Box[4]=Box[3]; Box[5]=Box[3];
+    if (debug>0) mprintf("    %s: Setting box to be a truncated octahedron, angle is %lf\n",
+                         parmName,Box[3]);
+  } else if (beta == 60.0) {
+    BoxType=2;
+    ifbox=1;
+    Box[0]=bx; Box[1]=by; Box[2]=bz;
+    Box[3]=60.0; Box[4]=90.0; Box[5]=60.0;
+    if (debug>0) 
+      mprintf("    %s: Setting box to be a rhombic dodecahedron, alpha=gamma=60.0, beta=90.0\n",
+              parmName);
+  } else {
+    BoxType=2;
+    ifbox=1;
+    Box[0]=bx; Box[1]=by; Box[2]=bz;
+    Box[3]=beta; Box[4]=beta; Box[5]=beta;
+    mprintf("    Warning: %s: Unrecognized box type, beta is %lf\n",beta);
+  }
+  return 0;  
+}
+
 /*
  * AmberParm::AtomInfo()
  * Print parm information for atom.
@@ -829,8 +904,8 @@
  */
 void AmberParm::Info(char *buffer) {
 
-  sprintf(buffer,"%i atoms, %i res, box %i, %i mol, %i solvent mol, %i frames",
-          natom,nres,ifbox,molecules,solventMolecules,parmFrames);
+  sprintf(buffer,"%i atoms, %i res, boxtype %i, %i mol, %i solvent mol, %i frames",
+          natom,nres,BoxType,molecules,solventMolecules,parmFrames);
 }
   
 // NOTE: The following atomToX functions do not do any memory checks!
@@ -963,6 +1038,7 @@
   return bonds;
 }
 
+// ---------===========================================================---------
 /*
  * AmberParm::modifyStateByMap()
  * Currently only intended for use with AtomMap.
@@ -1028,24 +1104,21 @@
   // Set up new parm information
   newParm->natom = this->natom;
   newParm->nres = this->nres;
-  newParm->ifbox = this->ifbox;
   newParm->parmFrames = this->parmFrames;
 
   // Give mapped parm the same pindex as original parm
   newParm->pindex = this->pindex;
 
   // Copy box information
-  if (this->Box!=NULL) {
-    newParm->Box=(double*) malloc(4*sizeof(double));
-    for (int i=0; i<4; i++)
-      newParm->Box[i] = this->Box[i];
-  }
+  for (int i=0; i<6; i++)
+    newParm->Box[i] = this->Box[i];
+  newParm->BoxType=this->BoxType;
 
   // Set values up
   // NOTE: Eventually set all pointers up?
   newParm->values[NATOM] = newParm->natom;
   newParm->values[NRES] = newParm->nres;
-  newParm->values[IFBOX] = newParm->ifbox;
+  newParm->values[IFBOX] = this->values[IFBOX];
 
   return newParm;
 }
@@ -1152,7 +1225,6 @@
   // Set up new parm information
   newParm->natom = j;
   newParm->nres = jres+1; 
-  newParm->ifbox = this->ifbox;
   newParm->parmFrames = this->parmFrames;
   if (this->molecules>0) 
     newParm->molecules = jmol+1;
@@ -1184,26 +1256,25 @@
   }
   
   // Copy box information
-  if (this->Box!=NULL) {
-    newParm->Box=(double*) malloc(4*sizeof(double));
-    for (i=0; i<4; i++)
-      newParm->Box[i] = this->Box[i];
-  }
+  for (i=0; i<6; i++)
+    newParm->Box[i] = this->Box[i];
+  newParm->BoxType=this->BoxType;
 
   // Set values up
   // NOTE: Eventually set all pointers up?
   newParm->values[NATOM] = newParm->natom;
   newParm->values[NRES] = newParm->nres;
-  newParm->values[IFBOX] = newParm->ifbox;
+  newParm->values[IFBOX] = this->values[IFBOX];
 
   mprintf("           New parmtop contains %i atoms.\n",newParm->natom);
   mprintf("                                %i residues.\n",newParm->nres);
   mprintf("                                %i molecules.\n",newParm->molecules);
-  mprintf("                                %i solvent molcules.\n",newParm->solventMolecules);
+  mprintf("                                %i solvent molecules.\n",newParm->solventMolecules);
 
   return newParm;
 }
 
+// ---------===========================================================---------
 /* 
  * AmberParm::DataToBuffer()
  * Return char buffer containing N data elements stored in I, D, or C with 
@@ -1267,6 +1338,7 @@
   char *buffer;
   int solvent_pointer[3];
   int atom;
+  double parmBox[4];
 
   if (parmName==NULL) return 1;
 
@@ -1346,7 +1418,7 @@
   }
 
   // SOLVENT POINTERS
-  if (ifbox>0) {
+  if (values[IFBOX]>0) {
     PrintFlagFormat(&outfile, "%FLAG SOLVENT_POINTERS", "%FORMAT(3I8)");
     solvent_pointer[0]=finalSoluteRes;
     solvent_pointer[1]=molecules;
@@ -1360,8 +1432,12 @@
     outfile.IO->Write(buffer, sizeof(char), BufferSize);
 
     // BOX DIMENSIONS
+    parmBox[0] = Box[4]; // beta
+    parmBox[1] = Box[0]; // boxX
+    parmBox[2] = Box[1]; // boxY
+    parmBox[3] = Box[2]; // boxZ
     PrintFlagFormat(&outfile, "%FLAG BOX_DIMENSIONS", "%FORMAT(5E16.8)");
-    buffer = DataToBuffer(buffer,"%FORMAT(5E16.8)", NULL, Box, NULL, 4);
+    buffer = DataToBuffer(buffer,"%FORMAT(5E16.8)", NULL, parmBox, NULL, 4);
     outfile.IO->Write(buffer, sizeof(char), BufferSize);
   }
 
diff -urN amber11.orig/AmberTools/src/cpptraj/src/AmberParm.h amber11/AmberTools/src/cpptraj/src/AmberParm.h
--- amber11.orig/AmberTools/src/cpptraj/src/AmberParm.h	2011-04-14 15:30:16.000000000 +0300
+++ amber11/AmberTools/src/cpptraj/src/AmberParm.h	2011-10-25 15:01:28.117288634 +0300
@@ -75,14 +75,14 @@
     int *resnums;         // IPRES(NRES) 
     int natom;            // NATOM
     int nres;             // NRES
-    int ifbox;            // IFBOX
     int finalSoluteRes;   // IPTRES
     int molecules;        // NSPM
     int firstSolvMol;     // NSPSOL
     int *atomsPerMol;     // NSP(NSPM)
     double *mass;         // AMASS(NATOM)
     double *charge;       // CHARGE(NATOM)
-    double *Box;          // OLDBETA, BOX(1), BOX(2), BOX(3)
+    double Box[6];        // X, Y, Z, alpha, beta, gamma 
+    int BoxType;          // 0: None, 1: Orthogonal, 2: Non-orthogonal
     // From Ptraj
     char *solventMask;         // T for atoms in the solvent
     int solventMolecules;      // number of solvent molecules
@@ -96,6 +96,7 @@
     ~AmberParm();
     void ResName(char *, int);
     int OpenParm(char *);
+    int SetBoxInfo(double,double,double,double);
     int SetSurfaceInfo();
     int SetSolventInfo();
     void AtomInfo(int);
diff -urN amber11.orig/AmberTools/src/cpptraj/src/AmberRestart.cpp amber11/AmberTools/src/cpptraj/src/AmberRestart.cpp
--- amber11.orig/AmberTools/src/cpptraj/src/AmberRestart.cpp	2011-04-14 15:30:16.000000000 +0300
+++ amber11/AmberTools/src/cpptraj/src/AmberRestart.cpp	2011-10-25 15:01:28.125288626 +0300
@@ -109,7 +109,7 @@
   frameSize+=frameSize;
 
   // If box coords are present, allocate extra space for them
-  if (isBox>0) {
+  if (BoxType!=0) {
     numBoxCoords=6;
     frameSize+=((numBoxCoords*12)+1);
   }
@@ -167,13 +167,12 @@
 
   // If 0 probably at EOF. No box or velo.
 //  } else if (lineSize==0) {
-    isBox=0;
+    BoxType=0;
     hasVelocity=0;
 
   // If 36 or 72 (+1 newline) box info.
   } else if (lineSize==37 || lineSize==73) {
-    isBox=1;
-    numBoxCoords = (lineSize-1) / 12;
+    getBoxType(frameBuffer,lineSize);
     hasVelocity=0;
 
   // If filled framebuffer again, has velocity info. Check for box after velocity.
@@ -182,15 +181,14 @@
     if (File->IO->Gets(buffer,82)==0) {
       lineSize=strlen(buffer);
       if (lineSize==37 || lineSize==73) {
-        isBox=1;
-        numBoxCoords = (lineSize-1) / 12;
+        getBoxType(buffer,lineSize);
       } else {
         mprintf("Error: AmberRestart::SetupRead():\n");
         mprintf("       Expect only 3 or 6 box coords in box coord line.\n");
         return 1;
       }
     } else
-      isBox=0;
+      BoxType=0;
 
   // Otherwise, who knows what was read?
   } else {
@@ -202,12 +200,12 @@
 
   // Recalculate the frame size
   if (hasVelocity) frameSize+=frameSize;
-  if (isBox) frameSize+=( (numBoxCoords*12) + 1 );
+  if (BoxType!=0) frameSize+=( (numBoxCoords*12) + 1 );
   frameBuffer=(char*) realloc(frameBuffer, frameSize*sizeof(char));
 
   if (debug > 0) {
-    mprintf("    Amber Restart isBox=%i hasVelocity=%i numBoxCoords=%i\n",
-            isBox,hasVelocity,numBoxCoords);
+    mprintf("    Amber Restart BoxType=%i hasVelocity=%i numBoxCoords=%i\n",
+            BoxType,hasVelocity,numBoxCoords);
     mprintf("    Amber Restart frameSize= %i\n",frameSize);
   }
   
@@ -221,6 +219,20 @@
 }
 
 /*
+ * AmberRestart::getBoxType()
+ * Based on input buffer, determine box type and num box coords.
+ */
+void AmberRestart::getBoxType(char *boxline, int boxlineSize) {
+  double box[6];
+  numBoxCoords = (boxlineSize-1) / 12;
+  if (numBoxCoords>3) {
+    sscanf(boxline, "%8lf%8lf%8lf%8lf%8lf%8lf",box,box+1,box+2,box+3,box+4,box+5);
+    CheckBoxType(box);
+  } else
+    BoxType = P->BoxType;
+}
+
+/*
  * AmberRestart::getFrame()
  * Get the restart file frame. If velocities are present, read those too.
  */
@@ -251,7 +263,7 @@
     //F->V->printAtomCoord(0);
   }
   // Convert box to Frame if present
-  if (isBox) {
+  if (BoxType!=0) {
     if ( (bufferPosition = F->BufferToBox(bufferPosition, numBoxCoords, 12))==NULL ) {
       mprintf("Error: AmberRestart::getFrame: * detected in box coordinates of %s\n",
               trajfilename);
@@ -287,7 +299,7 @@
   if (F->V!=NULL)  // NOTE: Use hasVelocity in addition/instead?
     bufferPosition = F->V->FrameToBuffer(bufferPosition,"%12.7lf",12,6);
   // Write box to buffer
-  if (isBox)
+  if (BoxType!=0)
     bufferPosition = F->BoxToBuffer(bufferPosition, numBoxCoords, "%12.7lf",12);
 
   //if (seekable) fseek(fp, titleSize+(set*frameSize),SEEK_SET);
@@ -297,6 +309,8 @@
 
   File->IO->Close();
 
+  currentFrame++;
+
   return 0;
 }
 
diff -urN amber11.orig/AmberTools/src/cpptraj/src/AmberRestart.h amber11/AmberTools/src/cpptraj/src/AmberRestart.h
--- amber11.orig/AmberTools/src/cpptraj/src/AmberRestart.h	2011-04-14 15:30:16.000000000 +0300
+++ amber11/AmberTools/src/cpptraj/src/AmberRestart.h	2011-10-25 15:01:28.118288633 +0300
@@ -14,6 +14,7 @@
 
   int SetupRead();
   int SetupWrite();
+  void getBoxType(char *, int);
 
   public:
 
diff -urN amber11.orig/AmberTools/src/cpptraj/src/AmberRestartNC.cpp amber11/AmberTools/src/cpptraj/src/AmberRestartNC.cpp
--- amber11.orig/AmberTools/src/cpptraj/src/AmberRestartNC.cpp	2011-04-14 15:30:16.000000000 +0300
+++ amber11/AmberTools/src/cpptraj/src/AmberRestartNC.cpp	2011-10-25 15:01:28.126288625 +0300
@@ -60,15 +60,8 @@
 /*
  * AmberRestartNC::open()
  * Open up Netcdf restart file and set all dimension and variable IDs.
- * This is done every time the file is opened up since Im not sure
- * the variable IDs stay the same throughout each opening.
- * Could eventually be separated.
- * NOTE: Replace attrText allocs with static buffer? 
  */
 int AmberRestartNC::open() {
-  char *attrText; // For checking conventions and version 
-  int spatial; // For checking spatial dimensions
-
   mprintf("DEBUG: AmberRestartNC::open() called for %s, ncid=%i\n",File->filename,ncid);
   // If already open, return
   if (ncid!=-1) return 0;
@@ -94,6 +87,23 @@
   // Netcdf files are always seekable
   seekable=1;
 
+  return 0;
+}
+
+/*
+ * AmberRestartNC::SetupRead()
+ * Set up netcdf restart file for reading, get all variable and dimension IDs. 
+ * Also check number of atoms against associated parmtop.
+ * NOTE: Replace attrText allocs with static buffer? 
+ */
+int AmberRestartNC::SetupRead() {
+  char *attrText; // For checking conventions and version 
+  int spatial; // For checking spatial dimensions
+  double box[6];
+  size_t start[2], count[2];
+
+  if (open()) return 1;
+
   // Get global attributes
   if (title==NULL) title = GetAttrText(ncid,NC_GLOBAL, "title");
   attrText = GetAttrText(ncid,NC_GLOBAL, "Conventions");
@@ -147,20 +157,18 @@
   mprintf("    Netcdf restart time= %lf\n",restartTime);
 
   // Box info
-  // NOTE: If no box info found in parm should really try to determine correct
-  //       box type from angles. 
   if ( nc_inq_varid(ncid,"cell_lengths",&cellLengthVID)==NC_NOERR ) {
     if (checkNCerr(nc_inq_varid(ncid,"cell_angles",&cellAngleVID),
       "Getting cell angles.")!=0) return 1;
-    if (debug>0) mprintf("  Netcdf restart Box information found.\n"); 
-    if (P->ifbox==0) {
-      mprintf("Warning: Netcdf restart file contains box info but no box info found\n");
-      mprintf("         in associated parmfile %s; defaulting to orthogonal.\n",
-              P->parmName);
-      isBox=1;
-    } else {
-      isBox=P->ifbox;
-    }
+    if (debug>0) mprintf("  Netcdf restart Box information found.\n");
+    // Determine box type from angles
+    start[0]=0; start[1]=0;
+    count[0]=3; count[1]=0; 
+    if ( checkNCerr(nc_get_vara_double(ncid, cellLengthVID, start, count, box),
+                    "Getting cell lengths.")!=0 ) return 1;
+    if ( checkNCerr(nc_get_vara_double(ncid, cellAngleVID, start, count, box+3),
+                    "Getting cell angles.")!=0 ) return 1;
+    CheckBoxType(box);
   } 
 
   // Replica Temperatures
@@ -175,15 +183,6 @@
   //int cell_spatialDID, cell_angularDID;
   //int spatialVID, cell_spatialVID, cell_angularVID;
 
-  return 0;
-}
-
-/*
- * AmberRestartNC::SetupRead()
- * Just a frontend to open for now. Also check number of atoms.
- */
-int AmberRestartNC::SetupRead() {
-  if (open()) return 1;
   if (ncatom!=P->natom) {
     mprintf("Warning: Number of atoms in NetCDF restart file %s (%i) does not\n",
             File->filename,ncatom);
@@ -210,7 +209,7 @@
  */
 int AmberRestartNC::setupWriteForSet(int set) {
   int dimensionID[NC_MAX_VAR_DIMS];
-  size_t start[3], count[3];
+  size_t start[2], count[2];
   char buffer[1024];
   char xyz[3];
   char abc[15] = { 'a', 'l', 'p', 'h', 'a', 
@@ -280,7 +279,7 @@
     "Defining cell angular variable.")) return 1;
 
   // Box Info
-  if (isBox>0) {
+  if (BoxType!=0) {
     dimensionID[0]=cell_spatialDID;
     if (checkNCerr(nc_def_var(ncid,"cell_lengths",NC_DOUBLE,1,dimensionID,&cellLengthVID),
       "Defining cell length variable.")) return 1;
@@ -388,9 +387,9 @@
   }
 
   // Read box info 
-  if (isBox!=0) {
-    count [0]=3;
-    count [1]=0;
+  if (BoxType!=0) {
+    count[0]=3;
+    count[1]=0;
     if ( checkNCerr(nc_get_vara_double(ncid, cellLengthVID, start, count, F->box),
                     "Getting cell lengths.")!=0 ) return 1;
     if ( checkNCerr(nc_get_vara_double(ncid, cellAngleVID, start, count, &(F->box[3])),
@@ -425,7 +424,7 @@
   }
 
   // write box
-  if (isBox>0 && cellLengthVID!=-1) {
+  if (BoxType!=0 && cellLengthVID!=-1) {
     count[0]=3;
     count[1]=0;
     if (checkNCerr(nc_put_vara_double(ncid,cellLengthVID,start,count,F->box),
@@ -439,6 +438,8 @@
   // Close file for this set
   close();
 
+  currentFrame++;
+
   return 0;
 }  
 
diff -urN amber11.orig/AmberTools/src/cpptraj/src/AmberTraj.cpp amber11/AmberTools/src/cpptraj/src/AmberTraj.cpp
--- amber11.orig/AmberTools/src/cpptraj/src/AmberTraj.cpp	2011-04-14 15:30:16.000000000 +0300
+++ amber11/AmberTools/src/cpptraj/src/AmberTraj.cpp	2011-10-25 15:01:28.226288534 +0300
@@ -104,16 +104,18 @@
     rprintf("Error: AmberTraj::getFrame: * detected in coordinates of %s\n",trajfilename);
     return 1;  
   } 
-  if (isBox) { 
+  if (BoxType!=0) { 
     if ( (bufferPosition = F->BufferToBox(bufferPosition,numBoxCoords,8))==NULL ) {
       rprintf("Error: AmberTraj::getFrame: * detected in box coordinates of %s\n",
               trajfilename);
       return 1;
     }
-    // Set box angles to parmtop default
-    F->box[3] = P->Box[0];
-    F->box[4] = P->Box[0];
-    F->box[5] = P->Box[0];
+    // Set box angles to parmtop default if not read in
+    if (numBoxCoords==3) {
+      F->box[3] = P->Box[3];
+      F->box[4] = P->Box[4];
+      F->box[5] = P->Box[5];
+    }
   }
   return 0;
 }
@@ -138,20 +140,22 @@
   }
 
   bufferPosition = F->FrameToBuffer(bufferPosition,"%8.3lf",8,10);
-  if (isBox) 
+  if (BoxType!=0) 
     bufferPosition = F->BoxToBuffer(bufferPosition,numBoxCoords,"%8.3lf",8);
 
   outFrameSize = (int) (bufferPosition - frameBuffer);
   
   //if (seekable) 
   // NOTE: Seek only needs to happen when traj file changes
-  offset = (off_t) set;
-  offset *= (off_t) outFrameSize;
-  offset += (off_t) titleSize;
-  File->IO->Seek( offset);
+  //offset = (off_t) currentFrame;
+  //offset *= (off_t) outFrameSize;
+  //offset += (off_t) titleSize;
+  //File->IO->Seek( offset);
 
   if (File->IO->Write(frameBuffer,sizeof(char),outFrameSize)) return 1;
 
+  currentFrame++;
+ 
   return 0;
 }
 
@@ -163,8 +167,10 @@
 int AmberTraj::SetupRead() {
   char buffer[BUFFER_SIZE];
   int frame_lines;
-  int lineSize;
-  long long int file_size, frame_size;
+  int lineSize, maxi;
+  bool sizeFound;
+  long long int file_size, frame_size, title_size, tmpfsize;
+  double box[6]; // For checking box coordinates
 
   // Attempt to open the file. open() sets the title and titleSize
   if (open()) return 1;
@@ -219,7 +225,7 @@
 
     if (strncmp(buffer,"REMD",4)==0 || strncmp(buffer,"HREMD",5)==0) {
       // REMD header - no box coords
-      isBox=0;
+      BoxType=0;
     } else if (lineSize<80) {
       /* Line is shorter than 80 chars, indicates box coords.
        * Length of the line HAS to be a multiple of 8, and probably could be
@@ -228,53 +234,108 @@
        */
       if (debug>0) mprintf("    Box line is %i chars.\n",lineSize);
       if ( ((lineSize-1)%24)!=0 ) {
-        rprintf("Error in box coord line.\nExpect only 3 or 6 box coords.\n");
+        mprintf("Error in box coord line. Expect only 3 or 6 box coords.\n");
         return 1;
       }
       numBoxCoords=(lineSize-1) / 8;
       if (debug>0) mprintf("    Detected %i box coords.\n",numBoxCoords);
-      frameSize+=lineSize;
-      // Reallocate frame buffer accordingly
-      frameBuffer=(char*) realloc(frameBuffer,frameSize * sizeof(char));
-      if (P->ifbox==0) { 
-        rprintf( "Warning: Box coords detected in trajectory but not defined in topology!\n");
-        rprintf("         Setting box type to rectangular.\n");
-        isBox=1;
+      // Determine box type based on angles. Angles are usually not printed 
+      // for orthogonal and truncated octahedral boxes, but check here just
+      // to be safe. If no angles present use parmtop Box Type.
+      if (numBoxCoords>3) {
+        sscanf(buffer, "%8lf%8lf%8lf%8lf%8lf%8lf",box,box+1,box+2,box+3,box+4,box+5);
+        CheckBoxType(box);
       } else {
-        isBox = P->ifbox;
+        BoxType = P->BoxType; 
+        // If box coords are present in traj but no box in parm, print warning
+        if (BoxType == 0)
+          mprintf("Warning: %s has box coordinates but no box info in %s\n",
+                  trajfilename,P->parmName);
       }
+      // Reallocate frame buffer accordingly
+      frameSize+=lineSize;
+      frameBuffer=(char*) realloc(frameBuffer,frameSize * sizeof(char));
     }
   }
 
-  /* Calculate number of frames. If not possible and this is not a
-   * compressed file the trajectory is probably corrupted. Frames will
-   * be read until EOF.
-   * NOTE: It is necessary to use the stat command to get the file size
-   * instead of fseek in case the file has been popen'd.
-   * NOTE: Need the uncompressed file size!
-   */
-  if (File->uncompressed_size>0)
-    file_size = File->uncompressed_size;
-  else
-    file_size=File->file_size;
+  // Calculate Frames and determine seekable. If not possible and this is not a
+  // compressed file the trajectory is probably corrupted. Frames will
+  // be read until EOF (Frames = -1).
   if (debug>0)
-    rprintf("Title offset=%i FrameSize=%i UncompressedFileSize=%lli\n",
-            titleSize,frameSize,file_size);
-  frame_size = (long long int) titleSize;
-  file_size = file_size - frame_size; // Subtract title size from file total size.
+    rprintf("Title offset=%i FrameSize=%i Size=%lu UncompressedFileSize=%lu\n",
+            titleSize,frameSize,File->file_size,File->uncompressed_size);
+  title_size = (long long int) titleSize;
   frame_size = (long long int) frameSize;
-  Frames = (int) (file_size / frame_size);
-
-  if ( (file_size % frame_size) == 0 ) {
-    seekable = 1;
-    stop = Frames;
+  // -----==== AMBER TRAJ COMPRESSED ====------
+  if (File->compressType!=NONE) {
+    // If the uncompressed size of compressed file is reported as <= 0,
+    // uncompressed size cannot be determined. Read coordinates until
+    // EOF.
+    if (File->uncompressed_size <= 0) {
+      mprintf("Warning: %s: Uncompressed size of trajectory could not be determined.\n",
+              File->filename);
+      if (File->compressType==BZIP2)
+        mprintf("         (This is normal for bzipped files)\n");
+      mprintf("         Number of frames could not be calculated.\n");
+      mprintf("         Frames will be read until EOF.\n");
+      Frames = -1;
+      seekable = false;
+    } else {
+      file_size = File->uncompressed_size;
+      file_size = file_size - title_size;
+      // Frame calculation for large gzip files
+      // If uncompressed size is less than compressed size, uncompressed
+      // size is likely > 4GB.
+      if (File->compressType == GZIP && file_size < (long long int) File->file_size) {
+        // Since this is gzip compressed, if the file_size % frame size != 0, 
+        // it could be that the uncompressed filesize > 4GB. Since 
+        //   ISIZE = uncompressed % 2^32, 
+        // try ((file_size + (2^32 * i)) % frame_size) and see if any are 0.
+        if ( (file_size % frame_size) != 0) {
+          // Determine the maximum number of iterations to try based on the
+          // fact that Amber trajectories typically compress about 3x with
+          // gzip.
+          tmpfsize = ((File->file_size * 4) - File->uncompressed_size) / 4294967296LL;
+          maxi = (int) tmpfsize;
+          maxi++;
+          if (debug>1)
+            mprintf("\tLooking for uncompressed gzip size > 4GB, %i iterations.\n",maxi);
+          tmpfsize = 0;
+          sizeFound=false;
+          for (int i = 0; i < maxi; i++ ) {
+            tmpfsize = (4294967296LL * i) + file_size;
+            if ( (tmpfsize % frame_size) == 0) {sizeFound=true; break;}
+          }
+          if (sizeFound) file_size = tmpfsize;
+        }
+      }
+      if ((file_size % frame_size) == 0) {
+        Frames = (int) (file_size / frame_size);
+        seekable = true;
+      } else {
+        mprintf("Warning: %s: Number of frames in compressed traj could not be determined.\n",
+                File->filename);
+        mprintf("         Frames will be read until EOF.\n");
+        Frames=-1;
+        seekable=false;
+      }
+    }
+  // ----==== AMBER TRAJ NOT COMPRESSED ====----
   } else {
-    stop = -1;
-    Frames = -1;
-    seekable = 0;
-    mprintf("  Error: Could not predict # frames in %s. This usually indicates \n",File->filename);
-    mprintf("         a corrupted trajectory. Frames will be read until EOF.\n");
+    file_size = File->file_size;
+    file_size = file_size - title_size;
+    Frames = (int) (file_size / frame_size);
+    if ( (file_size % frame_size) == 0 ) {
+      seekable = true;
+    } else {
+      mprintf("Warning: %s: Could not accurately predict # frames. This usually \n",
+              File->filename);
+      mprintf("         indicates a corrupted trajectory. Will attempt to read %i frames.\n",
+              Frames);
+      seekable=false;
+    }
   }
+  stop = Frames;
 
   if (debug>0)
     rprintf("Atoms: %i FrameSize: %i TitleSize: %i NumBox: %i Seekable: %i Frames: %i\n\n", 
@@ -315,7 +376,9 @@
   frameSize += hasREMD;
 
   // If box coords are present, allocate extra space for them
-  if (isBox>0) {
+  // NOTE: Currently only writing box lengths for all box types. This means
+  //       writing triclinic box type is currently not supported.
+  if (BoxType!=0) {
     numBoxCoords=3; // Only write out box lengths for trajectories
     frameSize+=((numBoxCoords*8)+1);
   }
diff -urN amber11.orig/AmberTools/src/cpptraj/src/ArgList.cpp amber11/AmberTools/src/cpptraj/src/ArgList.cpp
--- amber11.orig/AmberTools/src/cpptraj/src/ArgList.cpp	2011-04-14 15:30:16.000000000 +0300
+++ amber11/AmberTools/src/cpptraj/src/ArgList.cpp	2011-10-25 15:01:28.101288648 +0300
@@ -146,8 +146,16 @@
   return argline;
 }
 
-/*
- * ArgList::Command()
+/* ArgList::Arg()
+ * Return arg at specified position.
+ */
+char *ArgList::Arg(int pos) {
+  if (pos>-1 && pos<nargs) 
+    return arglist[pos];
+  return NULL;
+}
+
+/* ArgList::Command()
  * Check the first arg for command
  * Mark and return. Return even if marked.
  */
diff -urN amber11.orig/AmberTools/src/cpptraj/src/ArgList.h amber11/AmberTools/src/cpptraj/src/ArgList.h
--- amber11.orig/AmberTools/src/cpptraj/src/ArgList.h	2011-04-14 15:30:16.000000000 +0300
+++ amber11/AmberTools/src/cpptraj/src/ArgList.h	2011-10-25 15:01:28.102288647 +0300
@@ -15,6 +15,7 @@
     void Add(char *);
     void print();
     char *ArgLine();
+    char *Arg(int);
 
     char *Command();
     int CommandIs(const char *);
diff -urN amber11.orig/AmberTools/src/cpptraj/src/AtomMap.cpp amber11/AmberTools/src/cpptraj/src/AtomMap.cpp
--- amber11.orig/AmberTools/src/cpptraj/src/AtomMap.cpp	2011-04-14 15:30:16.000000000 +0300
+++ amber11/AmberTools/src/cpptraj/src/AtomMap.cpp	2011-10-25 15:01:28.081288667 +0300
@@ -6,6 +6,7 @@
 #include "TorsionRoutines.h"
 // DEBUG
 #include "Mol2File.h"
+#include <cstdio>
 
 //--------- PRIVATE ROUTINES ---------------------------------------
 /*
@@ -463,21 +464,27 @@
 // Write atommap out as a mol2 file, useful for checking bond info
 void atommap::WriteMol2(char *m2filename) {
   Mol2File outfile;
-  int *Selected;
+  AtomMask M1;
   // Temporary parm to play with
   AmberParm *tmpParm;
+  Frame *tmpFrame;
 
   // Create mask containing all atoms
-  Selected = (int*) malloc(natom*sizeof(int));
-  for (int atom=0; atom<natom; atom++) Selected[atom]=atom;
+  //for (int atom=0; atom<natom; atom++) Selected[atom]=atom;
   // Fake strip, just use as crap way to copy
-  tmpParm = P->modifyStateByMask(Selected,natom);
-  free(Selected);
+  //tmpParm = P->modifyStateByMask(Selected,natom);
+  //free(Selected);
   // Modify the bonds array to include this info
-  tmpParm->ResetBondInfo();
-  for (int atom=0; atom<natom; atom++) 
-    for (int bond=0; bond < M[atom].nbond; bond++) 
-      tmpParm->AddBond(atom, M[atom].bond[bond], 0);
+  //tmpParm->ResetBondInfo();
+  //for (int atom=0; atom<natom; atom++) 
+  //  for (int bond=0; bond < M[atom].nbond; bond++) 
+  //    tmpParm->AddBond(atom, M[atom].bond[bond], 0);
+  // Create mask with all mapped atoms
+  for (int atom=0; atom<natom; atom++) {if (M[atom].isMapped) M1.AddAtom(atom);}
+  // Strip so only mapped atoms remain
+  tmpParm = P->modifyStateByMask(M1.Selected,M1.Nselected);
+  tmpFrame = new Frame(M1.Nselected,NULL);
+  tmpFrame->SetFrameFromMask(F, &M1);
 
   // Trajectory Setup
   outfile.File=new PtrajFile();
@@ -485,13 +492,14 @@
   outfile.trajfilename = outfile.File->basefilename;
   outfile.debug=debug;
   outfile.SetTitle(m2filename);
-  outfile.P=P;
+  outfile.P=tmpParm;
   outfile.SetupWrite();
   outfile.open();
-  outfile.F=F;
+  outfile.F=tmpFrame;
   outfile.writeFrame(0);
   outfile.close();
   delete tmpParm;
+  delete tmpFrame;
 }
 // ============================================================================
 
@@ -1032,6 +1040,13 @@
   int iterations=0;
 
   numAtomsMapped=MapUniqueAtoms(Ref, Tgt);
+  // DEBUG
+  //char name[1024];
+  //sprintf(name,"Ref.%i.mol2",iterations);
+  //Ref->WriteMol2(name);
+  //sprintf(name,"Tgt.%i.mol2",iterations);
+  //Tgt->WriteMol2(name);
+  // END DEBUG
   if (debug>0)
     mprintf("*         MapUniqueAtoms: %i atoms mapped.\n",numAtomsMapped);
   if (numAtomsMapped==0) return 1;
@@ -1045,17 +1060,35 @@
     iterations++;
     // First assign based on bonds to unique (already mapped) atoms.
     numAtomsMapped=mapBondsToUnique(Ref,Tgt);
+    // DEBUG
+    //sprintf(name,"Ref.%i.u.mol2",iterations);
+    //Ref->WriteMol2(name);
+    //sprintf(name,"Tgt.%i.u.mol2",iterations);
+    //Tgt->WriteMol2(name);
+    // END DEBUG
     if (debug>0)
       mprintf("* [%3i] mapBondsToUnique: %i atoms mapped.\n",iterations,numAtomsMapped);
     if (numAtomsMapped<0) return 1;
     // Next assign based on chirality
     numAtomsMapped=mapChiral(Ref,Tgt);
+    // DEBUG
+    //sprintf(name,"Ref.%i.c.mol2",iterations);
+    //Ref->WriteMol2(name);
+    //sprintf(name,"Tgt.%i.c.mol2",iterations);
+    //Tgt->WriteMol2(name);
+    // END DEBUG
     if (debug>0)
       mprintf("* [%3i]        mapChiral: %i atoms mapped.\n",iterations,numAtomsMapped);
     if (numAtomsMapped<0) return 1;
     if (numAtomsMapped>0) continue;
     // Last assign based on index/element
     numAtomsMapped=mapByIndex(Ref,Tgt);
+    // DEBUG
+    //sprintf(name,"Ref.%i.i.mol2",iterations);
+    //Ref->WriteMol2(name);
+    //sprintf(name,"Tgt.%i.i.mol2",iterations);
+    //Tgt->WriteMol2(name);
+    // END DEBUG
     if (debug>0)
       mprintf("* [%3i]       mapByIndex: %i atoms mapped.\n",iterations,numAtomsMapped);
     if (numAtomsMapped<0) return 1;
diff -urN amber11.orig/AmberTools/src/cpptraj/src/CoordFileList.cpp amber11/AmberTools/src/cpptraj/src/CoordFileList.cpp
--- amber11.orig/AmberTools/src/cpptraj/src/CoordFileList.cpp	2011-04-14 15:30:16.000000000 +0300
+++ amber11/AmberTools/src/cpptraj/src/CoordFileList.cpp	2011-10-25 15:01:28.201288557 +0300
@@ -1,4 +1,5 @@
 // CoordFileList
+#include <cstddef>
 #include <cstring> // strcmp
 #include "CoordFileList.h"
 #include "CpptrajStdio.h"
diff -urN amber11.orig/AmberTools/src/cpptraj/src/cpptrajdepend amber11/AmberTools/src/cpptraj/src/cpptrajdepend
--- amber11.orig/AmberTools/src/cpptraj/src/cpptrajdepend	2011-04-14 15:30:16.000000000 +0300
+++ amber11/AmberTools/src/cpptraj/src/cpptrajdepend	2011-10-25 15:01:28.107288643 +0300
@@ -2,7 +2,7 @@
 AmberNetcdf.o : AmberNetcdf.cpp AmberNetcdf.h AmberParm.h ArgList.h AtomMask.h BaseFileIO.h CpptrajStdio.h Frame.h NetcdfRoutines.h PtrajFile.h Range.h TrajFile.h
 AmberParm.o : AmberParm.cpp AmberParm.h BaseFileIO.h CpptrajStdio.h Mol2FileRoutines.h PDBfileRoutines.h PtrajFile.h
 AmberTraj.o : AmberTraj.cpp AmberParm.h AmberTraj.h ArgList.h AtomMask.h BaseFileIO.h CpptrajStdio.h Frame.h PtrajFile.h Range.h TrajFile.h
-TrajFile.o : TrajFile.cpp AmberParm.h ArgList.h AtomMask.h BaseFileIO.h CpptrajStdio.h Frame.h PtrajFile.h PtrajMpi.h Range.h TrajFile.h
+TrajFile.o : TrajFile.cpp AmberParm.h ArgList.h AtomMask.h BaseFileIO.h CpptrajStdio.h Frame.h PtrajFile.h Range.h TrajFile.h
 Frame.o : Frame.cpp AmberParm.h AtomMask.h BaseFileIO.h CpptrajStdio.h DistRoutines.h Frame.h PtrajFile.h TorsionRoutines.h vectormath.h
 PtrajState.o : PtrajState.cpp Action.h AmberParm.h ArgList.h AtomMask.h BaseFileIO.h CoordFileList.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h Frame.h FrameList.h ParmFileList.h PtrajActionList.h PtrajFile.h PtrajMpi.h PtrajState.h Range.h ReferenceList.h TrajFile.h TrajinList.h TrajoutList.h
 ArgList.o : ArgList.cpp ArgList.h CpptrajStdio.h
@@ -13,7 +13,7 @@
 TrajoutList.o : TrajoutList.cpp AmberParm.h ArgList.h AtomMask.h BaseFileIO.h CoordFileList.h CpptrajStdio.h Frame.h ParmFileList.h PtrajFile.h Range.h TrajFile.h TrajoutList.h
 ReferenceList.o : ReferenceList.cpp AmberParm.h ArgList.h AtomMask.h BaseFileIO.h CoordFileList.h CpptrajStdio.h Frame.h FrameList.h ParmFileList.h PtrajFile.h Range.h ReferenceList.h TrajFile.h
 CoordFileList.o : CoordFileList.cpp AmberNetcdf.h AmberParm.h AmberRestart.h AmberRestartNC.h AmberTraj.h ArgList.h AtomMask.h BaseFileIO.h Conflib.h CoordFileList.h CpptrajStdio.h Frame.h Mol2File.h PDBfile.h ParmFileList.h PtrajFile.h Range.h RemdTraj.h TrajFile.h TrajinList.h
-PtrajActionList.o : PtrajActionList.cpp Action.h Action_Angle.h Action_Center.h Action_Closest.h Action_DSSP.h Action_Dihedral.h Action_Distance.h Action_Hbond.h Action_Image.h Action_Mask.h Action_NAstruct.h Action_Outtraj.h Action_Pucker.h Action_Radgyr.h Action_Rmsd.h Action_Strip.h Action_Surf.h AmberParm.h ArgList.h AtomMap.h AtomMask.h AxisType.h BaseFileIO.h CoordFileList.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h Frame.h FrameList.h ParmFileList.h PtrajActionList.h PtrajFile.h Range.h TrajFile.h TrajoutList.h
+PtrajActionList.o : PtrajActionList.cpp Action.h Action_Angle.h Action_Center.h Action_Closest.h Action_DSSP.h Action_Dihedral.h Action_Distance.h Action_Hbond.h Action_Image.h Action_Mask.h Action_NAstruct.h Action_Outtraj.h Action_Pucker.h Action_Radgyr.h Action_Rms2d.h Action_Rmsd.h Action_Strip.h Action_Surf.h AmberParm.h ArgList.h AtomMap.h AtomMask.h AxisType.h BaseFileIO.h CoordFileList.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h Frame.h FrameList.h ParmFileList.h PtrajActionList.h PtrajFile.h Range.h TrajFile.h TrajoutList.h
 DataSet.o : DataSet.cpp CpptrajStdio.h DataSet.h
 DataSetList.o : DataSetList.cpp CpptrajStdio.h DataSet.h DataSetList.h intDataSet.h mapDataSet.h stringDataSet.h
 vectormath.o : vectormath.cpp CpptrajStdio.h vectormath.h
@@ -23,9 +23,9 @@
 DataFile.o : DataFile.cpp BaseFileIO.h CpptrajStdio.h DataFile.h DataSet.h PtrajFile.h
 DataFileList.o : DataFileList.cpp BaseFileIO.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h PtrajFile.h
 PDBfile.o : PDBfile.cpp AmberParm.h ArgList.h AtomMask.h BaseFileIO.h CpptrajStdio.h Frame.h PDBfile.h PDBfileRoutines.h PtrajFile.h Range.h TrajFile.h
-PtrajFile.o : PtrajFile.cpp BaseFileIO.h Bzip2File.h CpptrajStdio.h GzipFile.h MpiFile.h NetcdfRoutines.h PDBfileRoutines.h PtrajFile.h PtrajMpi.h StdFile.h
+PtrajFile.o : PtrajFile.cpp BaseFileIO.h Bzip2File.h CpptrajStdio.h GzipFile.h Mol2FileRoutines.h MpiFile.h NetcdfRoutines.h PDBfileRoutines.h PtrajFile.h PtrajMpi.h StdFile.h
 PDBfileRoutines.o : PDBfileRoutines.cpp PDBfileRoutines.h
-AtomMap.o : AtomMap.cpp Action.h AmberParm.h ArgList.h AtomMap.h AtomMask.h BaseFileIO.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h Frame.h FrameList.h PDBfileRoutines.h ParmFileList.h PtrajFile.h TorsionRoutines.h
+AtomMap.o : AtomMap.cpp Action.h AmberParm.h ArgList.h AtomMap.h AtomMask.h BaseFileIO.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h Frame.h FrameList.h Mol2File.h ParmFileList.h PtrajFile.h Range.h TorsionRoutines.h TrajFile.h
 BaseFileIO.o : BaseFileIO.cpp BaseFileIO.h PtrajMpi.h
 StdFile.o : StdFile.cpp BaseFileIO.h StdFile.h
 GzipFile.o : GzipFile.cpp BaseFileIO.h CpptrajStdio.h GzipFile.h
@@ -62,3 +62,4 @@
 Action_Pucker.o : Action_Pucker.cpp Action.h Action_Pucker.h AmberParm.h ArgList.h AtomMask.h BaseFileIO.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h Frame.h FrameList.h ParmFileList.h PtrajFile.h
 Range.o : Range.cpp ArgList.h CpptrajStdio.h Range.h
 Action_Outtraj.o : Action_Outtraj.cpp Action.h Action_Outtraj.h AmberParm.h ArgList.h AtomMask.h BaseFileIO.h CoordFileList.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h Frame.h FrameList.h ParmFileList.h PtrajFile.h PtrajMpi.h Range.h TrajFile.h TrajoutList.h
+Action_Rms2d.o : Action_Rms2d.cpp Action.h Action_Rms2d.h AmberParm.h ArgList.h AtomMask.h BaseFileIO.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h Frame.h FrameList.h ParmFileList.h PtrajFile.h
diff -urN amber11.orig/AmberTools/src/cpptraj/src/DataFile.cpp amber11/AmberTools/src/cpptraj/src/DataFile.cpp
--- amber11.orig/AmberTools/src/cpptraj/src/DataFile.cpp	2011-04-14 15:30:16.000000000 +0300
+++ amber11/AmberTools/src/cpptraj/src/DataFile.cpp	2011-10-25 15:01:28.098288651 +0300
@@ -65,8 +65,49 @@
   isInverted=true;
 }
 
-/*
- * DataFile::AddSet()
+/* DataFile::SetPrecision()
+ * Set precision of the specified dataset to width.precision. If '*' specified 
+ * set for all datasets in file.
+ */
+void DataFile::SetPrecision(char *dsetName, int widthIn, int precisionIn) {
+  int precision, dset;
+  DataSet *Dset = NULL;
+
+  if (dsetName==NULL) {
+    mprintf("Error: SetPrecision must be called with dataset name or '*'.\n");
+    return;
+  }
+  if (widthIn<1) {
+    mprintf("Error: SetPrecision (%s): Cannot set width < 1.\n",filename);
+    return;
+  }
+  precision=precisionIn;
+  if (precisionIn<0) precision=0;
+  // If <dsetName>=='*' specified set precision for all data sets
+  if (dsetName[0]=='*') {
+    mprintf("    Setting width.precision for all sets in %s to %i.%i\n",
+            filename,widthIn,precision);
+    for (dset=0; dset<Nsets; dset++)
+      SetList[dset]->SetPrecision(widthIn,precision);
+
+  // Otherwise find dataset <dsetName> and set precision
+  } else {
+    mprintf("    Setting width.precision for dataset %s to %i.%i\n",
+            dsetName,widthIn,precision);
+    for (dset=0; dset<Nsets; dset++) {
+      if ( strcmp(SetList[dset]->Name(), dsetName)==0 ) {
+        Dset=SetList[dset];
+        break;
+      }
+    }
+    if (Dset!=NULL)
+      Dset->SetPrecision(widthIn,precision);
+    else
+      mprintf("Error: Dataset %s not found in datafile %s\n",dsetName,filename);
+  }
+}
+
+/* DataFile::AddSet()
  * Add given set to this datafile
  */
 int DataFile::AddSet(DataSet *D) {
diff -urN amber11.orig/AmberTools/src/cpptraj/src/DataFile.h amber11/AmberTools/src/cpptraj/src/DataFile.h
--- amber11.orig/AmberTools/src/cpptraj/src/DataFile.h	2011-04-14 15:30:16.000000000 +0300
+++ amber11/AmberTools/src/cpptraj/src/DataFile.h	2011-10-25 15:01:28.098288651 +0300
@@ -30,6 +30,7 @@
     void SetXlabel(char*);
     void SetInverted();
     void SetNoXcol();
+    void SetPrecision(char *, int, int);
     int AddSet(DataSet *);
     int NameIs(char *);
     void DataSetNames();
diff -urN amber11.orig/AmberTools/src/cpptraj/src/DataFileList.cpp amber11/AmberTools/src/cpptraj/src/DataFileList.cpp
--- amber11.orig/AmberTools/src/cpptraj/src/DataFileList.cpp	2011-04-14 15:30:16.000000000 +0300
+++ amber11/AmberTools/src/cpptraj/src/DataFileList.cpp	2011-10-25 15:01:28.201288557 +0300
@@ -1,4 +1,5 @@
 // DataFileList
+#include <cstddef>
 #include "DataFileList.h"
 #include "CpptrajStdio.h"
 
diff -urN amber11.orig/AmberTools/src/cpptraj/src/DataSet.cpp amber11/AmberTools/src/cpptraj/src/DataSet.cpp
--- amber11.orig/AmberTools/src/cpptraj/src/DataSet.cpp	2011-04-14 15:30:16.000000000 +0300
+++ amber11/AmberTools/src/cpptraj/src/DataSet.cpp	2011-10-25 15:01:28.094288655 +0300
@@ -1,6 +1,7 @@
 // DataSet
 #include <cstdlib>
 #include <cstring>
+#include <cstdio> // sprintf
 #include "DataSet.h"
 #include "CpptrajStdio.h"
 
@@ -11,12 +12,73 @@
   N=0;
   isDynamic=false;
   current=0;
+  width = 0;
+  precision = 0;
+  format = NULL;
+  dType = UNKNOWN_DATA;
 }
 
 // DESTRUCTOR
 DataSet::~DataSet() {
   //fprintf(stderr,"DataSet Destructor\n");
   if (name!=NULL) free(name);
+  if (format!=NULL) free(format);
+}
+
+/*
+ * DataSet::setFormatString()
+ * Set up the output format string for each data element based on the given 
+ * dataType and the current width, and precision.
+ */
+void DataSet::setFormatString() {
+  size_t stringWidth = 0;
+  int wWidth = 0;
+  int pWidth = 0;
+
+  if (format!=NULL) {free(format); format=NULL;}
+
+  // Calc num of chars necessary to hold width
+  wWidth = (width / 10) + 1;
+
+  switch (dType) {
+    case DOUBLE :
+      // Calc num of chars necessary to hold precision
+      pWidth = (precision / 10) + 1;
+      // String fmt: " %w.plf\0"
+      stringWidth = pWidth + wWidth + 6;
+      format = (char*) malloc( stringWidth * sizeof(char) );
+      sprintf(format, " %%%i.%ilf", width, precision);
+      break;
+    case STRING :
+      // String fmt: " %s"
+      format = (char*) malloc( 4 * sizeof(char) );
+      strcpy(format, " %s");
+      break;
+    case INT :
+      // String fmt: " %wi"
+      stringWidth = wWidth + 4;
+      format = (char*) malloc( stringWidth * sizeof(char) );
+      sprintf(format, " %%%ii", width);
+      break;
+    case UNKNOWN_DATA :
+      mprintf("Internal Error: setFormatString called with unknown data type.\n");
+  }
+
+  if (format==NULL) 
+    mprintf("Error: setFormatString: Could not allocate memory for string.\n");
+  // DEBUG
+  //else
+  //  mprintf("DEBUG: Format string: [%s]\n",format);
+}    
+
+/*
+ * DataSet::SetPrecision()
+ * Set dataset width and precision and recalc output format string.
+ */
+void DataSet::SetPrecision(int widthIn, int precisionIn) {
+  width=widthIn;
+  precision=precisionIn;
+  setFormatString();
 }
 
 /* 
diff -urN amber11.orig/AmberTools/src/cpptraj/src/DataSet.h amber11/AmberTools/src/cpptraj/src/DataSet.h
--- amber11.orig/AmberTools/src/cpptraj/src/DataSet.h	2011-04-14 15:30:16.000000000 +0300
+++ amber11/AmberTools/src/cpptraj/src/DataSet.h	2011-10-25 15:01:28.103288646 +0300
@@ -36,12 +36,18 @@
 class DataSet {
   protected:
     char *name;        // Name of the dataset
+    dataType dType;    // The dataset type
     int N;             // Number of data elements
     int current;       // The current data element
+    int width;         // The output width of a data element
+    int precision;     // The output precision of a data element (if applicable)
+    char *format;      // Format of output
+    
     bool isDynamic;    // True : N is not known, reallocate as N increases
                        // False: N is known, allocate for N
     // If not isDynamic, Allocate will reserve space for N data elements 
     virtual int Allocate( )      { return 0; }
+    void setFormatString();
 
   public:
 
@@ -51,13 +57,16 @@
     virtual int Xmax()               { return 0; }
     virtual int isEmpty(int)         { return 0; }
     virtual void Add( int, void * )  { return;   }
+    virtual int Get( void *, int )   { return 1; }
     virtual char *Write(char*, int)  { return 0; }
     virtual int Width()              { return 0; }
     virtual int Sync()               { return 0; }
 
+    void SetPrecision(int,int);
     int Setup(char*,int);
     void Info();
     char *Name() { return name; }
     int CheckSet();
+    dataType Type() {return dType;}
 };
 #endif 
diff -urN amber11.orig/AmberTools/src/cpptraj/src/Frame.cpp amber11/AmberTools/src/cpptraj/src/Frame.cpp
--- amber11.orig/AmberTools/src/cpptraj/src/Frame.cpp	2011-04-14 15:30:16.000000000 +0300
+++ amber11/AmberTools/src/cpptraj/src/Frame.cpp	2011-10-25 15:01:28.120288631 +0300
@@ -60,8 +60,38 @@
   if (Mass!=NULL) free(Mass);
 }
 
-/*
- * Frame::Copy()
+/* Frame::ZeroCoords()
+ * Set all coords to 0.0
+ */
+void Frame::ZeroCoords() {
+  for (int coord=0; coord < N; coord++)
+    X[coord]=0.0;
+}
+
+/* Frame::AddCoord()
+ * Add the coord values from the input frame to the coord values of 
+ * this frame.
+ */
+void Frame::AddCoord(Frame *FrameIn) {
+  if (FrameIn->N != this->N) {
+    mprintf("Error: Frame::AddCoord: Attempting to add %i coords to %i coords.\n",
+            FrameIn->N,this->N);
+  } else {
+    for (int coord=0; coord < N; coord++)
+      this->X[coord] += FrameIn->X[coord];
+  }
+}
+
+/* Frame::Divide()
+ * Divide all coord values by input. Dont do it if the number is too small.
+ */
+void Frame::Divide(double divisor) {
+  if (divisor < SMALL) return;
+  for (int coord=0; coord < N; coord++)
+    X[coord] /= divisor;
+}
+
+/* Frame::Copy()
  * Return a copy of the frame
  */
 Frame *Frame::Copy() {
@@ -423,21 +453,21 @@
  * Based on useMassIn, calculate geometric center (false) or center of mass 
  * (true) of the atoms in each mask.
  */
-double Frame::DIST2(AtomMask *Mask1, AtomMask *Mask2, bool useMassIn, int ifbox,
+double Frame::DIST2(AtomMask *Mask1, AtomMask *Mask2, bool useMassIn, int boxType,
                     double *ucell, double *recip) {
   double a1[3], a2[3];
 
   COM(Mask1, a1, useMassIn);
   COM(Mask2, a2, useMassIn);
 
-  if (ifbox == 0) 
+  if (boxType == 0) 
     return DIST2_NoImage(a1, a2);
-  else if (ifbox == 1) 
+  else if (boxType == 1) 
     return DIST2_ImageOrtho(a1, a2, this->box);
-  else if (ifbox == 2) 
+  else if (boxType == 2) 
     return DIST2_ImageNonOrtho(a1, a2, ucell, recip);
 
-  mprintf("    Error: Frame::DIST: Unrecognized box type (%i)\n.", ifbox);
+  mprintf("    Error: Frame::DIST: Unrecognized box type (%i)\n.", boxType);
 
   return (-1.0);
 }
@@ -445,8 +475,11 @@
 /*
  * Frame::DIST2()
  * Return the distance between atoms A1 and A2 with optional imaging.
+ *   0 = None
+ *   1 = Orthorhombic
+ *   2 = Non-orthorhombic
  */
-double Frame::DIST2(int A1, int A2, int ifbox, double *ucell, double *recip) {
+double Frame::DIST2(int A1, int A2, int boxType, double *ucell, double *recip) {
   int atom3;
   double a1[3], a2[3];
 
@@ -459,14 +492,14 @@
   a2[1] = X[atom3+1];
   a2[2] = X[atom3+2];
 
-  if (ifbox == 0)
+  if (boxType == 0)
     return DIST2_NoImage(a1, a2);
-  else if (ifbox == 1)
+  else if (boxType == 1)
     return DIST2_ImageOrtho(a1, a2, this->box);
-  else if (ifbox == 2) 
+  else if (boxType == 2) 
     return DIST2_ImageNonOrtho(a1, a2, ucell, recip);
 
-  mprintf("    Error: Frame::DIST: Unrecognized box type (%i)\n.", ifbox);
+  mprintf("    Error: Frame::DIST: Unrecognized box type (%i)\n.", boxType);
 
   return (-1.0);
 }
diff -urN amber11.orig/AmberTools/src/cpptraj/src/Frame.h amber11/AmberTools/src/cpptraj/src/Frame.h
--- amber11.orig/AmberTools/src/cpptraj/src/Frame.h	2011-04-14 15:30:16.000000000 +0300
+++ amber11/AmberTools/src/cpptraj/src/Frame.h	2011-10-25 15:01:28.092288656 +0300
@@ -22,6 +22,9 @@
     Frame(AtomMask *, double *);
     virtual ~Frame();             // Destructor is virtual since this class can be inherited
 
+    void ZeroCoords();
+    void AddCoord(Frame*);
+    void Divide(double);
     void printAtomCoord(int);
     void GetCoord(double *, int);
     void SetCoord(int, double *);
diff -urN amber11.orig/AmberTools/src/cpptraj/src/FrameList.cpp amber11/AmberTools/src/cpptraj/src/FrameList.cpp
--- amber11.orig/AmberTools/src/cpptraj/src/FrameList.cpp	2011-04-14 15:30:16.000000000 +0300
+++ amber11/AmberTools/src/cpptraj/src/FrameList.cpp	2011-10-25 15:01:28.073288675 +0300
@@ -13,7 +13,8 @@
     delete frameList[i];
 }
 
-/* FrameList::Add()
+/* 
+ * FrameList::Add()
  * Add given Frame to the FrameList. Store trajectory name that this frame
  * came from in frameNames. Store the associated parm in FrameParm. 
  */
@@ -29,6 +30,18 @@
   FrameParm.Add(P);
   Nframe++;
   return 0;
+}
+
+/*
+ * FrameList::Add()
+ * Add given Frame to the FrameList. Store the associated parm in FrameParm.
+ */
+int FrameList::Add(Frame *F, AmberParm *P) {
+  if (F==NULL || P==NULL) return 1;
+  frameList.push_back(F);
+  FrameParm.Add(P);
+  Nframe++;
+  return 0;
 }
 
 /* 
diff -urN amber11.orig/AmberTools/src/cpptraj/src/FrameList.h amber11/AmberTools/src/cpptraj/src/FrameList.h
--- amber11.orig/AmberTools/src/cpptraj/src/FrameList.h	2011-04-14 15:30:16.000000000 +0300
+++ amber11/AmberTools/src/cpptraj/src/FrameList.h	2011-10-25 15:01:28.073288675 +0300
@@ -18,6 +18,7 @@
     ~FrameList();
 
     int Add(Frame *, char *, AmberParm *,int);
+    int Add(Frame *, AmberParm *);
     AmberParm *GetFrameParm(int);
     int GetFrameIndex(char *);
     Frame *GetFrame(int idx);
diff -urN amber11.orig/AmberTools/src/cpptraj/src/GzipFile.cpp amber11/AmberTools/src/cpptraj/src/GzipFile.cpp
--- amber11.orig/AmberTools/src/cpptraj/src/GzipFile.cpp	2011-04-14 15:30:16.000000000 +0300
+++ amber11/AmberTools/src/cpptraj/src/GzipFile.cpp	2011-10-25 15:01:28.227288533 +0300
@@ -77,18 +77,24 @@
 
 /*
  * GzipFile::Read()
+ * NOTE: gzread returns 0 on EOF, -1 on error
  */
 int GzipFile::Read(void *buffer, size_t size, size_t count) {
   //size_t numread;
   int numread;
+  int expectedread;
+  
+  expectedread = (int)size;
+  expectedread *= (int)count;
   // Should never be able to call Read when fp is NULL.
   //if (fp==NULL) {
   //  fprintf(stdout,"Error: GzipFile::Read: Attempted to read NULL file pointer.\n");
   //  return 1;
   //}
-  numread = gzread(fp, buffer, size * count);
-  if (numread == -1) return -1;
-
+  numread = gzread(fp, buffer, expectedread);
+  if (numread != expectedread) return -1;
+  //if (numread < 1 ) return -1;
+  
   // NOTE: Check for errors here.
   return numread;
 }
diff -urN amber11.orig/AmberTools/src/cpptraj/src/intDataSet.cpp amber11/AmberTools/src/cpptraj/src/intDataSet.cpp
--- amber11.orig/AmberTools/src/cpptraj/src/intDataSet.cpp	2011-04-14 15:30:16.000000000 +0300
+++ amber11/AmberTools/src/cpptraj/src/intDataSet.cpp	2011-10-25 15:01:28.103288646 +0300
@@ -5,6 +5,14 @@
 #include "PtrajMpi.h"
 #include "CpptrajStdio.h"
 using namespace std;
+
+// CONSTRUCTOR
+intDataSet::intDataSet() {
+  width=12;
+  dType=INT;
+  setFormatString();
+}
+
 /*
  * intDataSet::Xmax(()
  * Return the maximum X value added to this set. By convention this is 
@@ -33,8 +41,21 @@
   current++;
 }
 
-/*
- * intDataSet::isEmpty()
+/* intDataSet::Get()
+ * Get data at frame, put into vOut. Return 1 if no data at frame.
+ */
+int intDataSet::Get(void *vOut, int frame) {
+  int *value;
+
+  if (vOut==NULL) return 1;
+  value = (int*) vOut;
+  it=Data.find( frame );
+  if (it == Data.end()) return 1;
+  *value = (*it).second;
+  return 0;
+}
+
+/* intDataSet::isEmpty()
  */
 int intDataSet::isEmpty(int frame) {
   it = Data.find( frame );
@@ -52,17 +73,17 @@
   it = Data.find( frame );
   if (it == Data.end()) 
     //sprintf(buffer," %12s","NoData");
-    sprintf(buffer," %12i", 0);
+    sprintf(buffer, format, 0);
   else 
-    sprintf(buffer," %12i",(*it).second);
-  return (buffer + 13);
+    sprintf(buffer, format, (*it).second);
+  return (buffer + width + 1);
 }
 
 /*
  * intDataSet::Width()
  */
 int intDataSet::Width() {
-  return 13;
+  return (width + 1);
 }
 
 /*
diff -urN amber11.orig/AmberTools/src/cpptraj/src/intDataSet.h amber11/AmberTools/src/cpptraj/src/intDataSet.h
--- amber11.orig/AmberTools/src/cpptraj/src/intDataSet.h	2011-04-14 15:30:16.000000000 +0300
+++ amber11/AmberTools/src/cpptraj/src/intDataSet.h	2011-10-25 15:01:28.104288645 +0300
@@ -11,9 +11,12 @@
     std::map<int,int> Data;
     std::map<int,int>::iterator it;
   public:
+    intDataSet();
+
     int Xmax();
     int isEmpty(int);
     void Add( int, void * );
+    int Get(void *, int);
     char *Write(char *, int);
     int Width();
     int Sync();
diff -urN amber11.orig/AmberTools/src/cpptraj/src/main.cpp amber11/AmberTools/src/cpptraj/src/main.cpp
--- amber11.orig/AmberTools/src/cpptraj/src/main.cpp	2011-04-14 15:30:16.000000000 +0300
+++ amber11/AmberTools/src/cpptraj/src/main.cpp	2011-10-25 15:01:28.228288532 +0300
@@ -6,7 +6,7 @@
 #include "PtrajMpi.h"
 #include "CpptrajStdio.h"
 #ifndef CPPTRAJ_VERSION_STRING
-#define CPPTRAJ_VERSION_STRING "V1.0.5"
+#define CPPTRAJ_VERSION_STRING "V1.1.1"
 #endif
 
 void Usage(char *programName) {
diff -urN amber11.orig/AmberTools/src/cpptraj/src/Makefile amber11/AmberTools/src/cpptraj/src/Makefile
--- amber11.orig/AmberTools/src/cpptraj/src/Makefile	2011-04-14 15:30:16.000000000 +0300
+++ amber11/AmberTools/src/cpptraj/src/Makefile	2011-10-25 15:01:28.075288673 +0300
@@ -16,7 +16,8 @@
        Action_Radgyr.cpp Conflib.cpp Action_Mask.cpp Action_Closest.cpp \
        NetcdfRoutines.cpp AmberRestartNC.cpp CpptrajStdio.cpp Mol2File.cpp \
        Mol2FileRoutines.cpp Action_NAstruct.cpp DistRoutines.cpp AxisType.cpp \
-       TorsionRoutines.cpp Action_Pucker.cpp Range.cpp Action_Outtraj.cpp
+       TorsionRoutines.cpp Action_Pucker.cpp Range.cpp Action_Outtraj.cpp \
+       Action_Rms2d.cpp
 
 OBJECTS=main.o AmberNetcdf.o AmberParm.o AmberTraj.o TrajFile.o \
         Frame.o PtrajState.o ArgList.o ptrajmask.o Action_Distance.o \
@@ -32,7 +33,8 @@
         Action_Radgyr.o Conflib.o Action_Mask.o Action_Closest.o \
         NetcdfRoutines.o AmberRestartNC.o CpptrajStdio.o Mol2File.o \
         Mol2FileRoutines.o Action_NAstruct.o DistRoutines.o AxisType.o \
-        TorsionRoutines.o Action_Pucker.o Range.o Action_Outtraj.o
+        TorsionRoutines.o Action_Pucker.o Range.o Action_Outtraj.o \
+        Action_Rms2d.o
 
 HEADERS=AmberNetcdf.h AmberParm.h AmberTraj.h TrajFile.h \
         Frame.h PtrajState.h ArgList.h ptrajmask.h Action.h Action_Distance.h \
@@ -48,7 +50,8 @@
         Action_Radgyr.h Conflib.h Action_Mask.h Action_Closest.h \
         NetcdfRoutines.h AmberRestartNC.h CpptrajStdio.h Mol2File.h \
         Mol2FileRoutines.h Action_NAstruct.h DistRoutines.h AxisType.h \
-        TorsionRoutines.h Action_Pucker.h Range.h Action_Outtraj.h
+        TorsionRoutines.h Action_Pucker.h Range.h Action_Outtraj.h \
+        Action_Rms2d.h
 
 all: cpptraj$(SFX)
 
diff -urN amber11.orig/AmberTools/src/cpptraj/src/Makefile_at amber11/AmberTools/src/cpptraj/src/Makefile_at
--- amber11.orig/AmberTools/src/cpptraj/src/Makefile_at	2011-04-14 15:30:16.000000000 +0300
+++ amber11/AmberTools/src/cpptraj/src/Makefile_at	2011-10-25 15:01:28.076288672 +0300
@@ -16,7 +16,8 @@
        Action_Radgyr.cpp Conflib.cpp Action_Mask.cpp Action_Closest.cpp \
        NetcdfRoutines.cpp AmberRestartNC.cpp CpptrajStdio.cpp Mol2File.cpp \
        Mol2FileRoutines.cpp Action_NAstruct.cpp DistRoutines.cpp AxisType.cpp \
-       TorsionRoutines.cpp Action_Pucker.cpp Range.cpp Action_Outtraj.cpp
+       TorsionRoutines.cpp Action_Pucker.cpp Range.cpp Action_Outtraj.cpp \
+       Action_Rms2d.cpp
 
 OBJECTS=main.o AmberNetcdf.o AmberParm.o AmberTraj.o TrajFile.o \
         Frame.o PtrajState.o ArgList.o ptrajmask.o Action_Distance.o \
@@ -32,7 +33,8 @@
         Action_Radgyr.o Conflib.o Action_Mask.o Action_Closest.o \
         NetcdfRoutines.o AmberRestartNC.o CpptrajStdio.o Mol2File.o \
         Mol2FileRoutines.o Action_NAstruct.o DistRoutines.o AxisType.o \
-        TorsionRoutines.o Action_Pucker.o Range.o Action_Outtraj.o
+        TorsionRoutines.o Action_Pucker.o Range.o Action_Outtraj.o \
+        Action_Rms2d.o
 
 HEADERS=AmberNetcdf.h AmberParm.h AmberTraj.h TrajFile.h \
         Frame.h PtrajState.h ArgList.h ptrajmask.h Action.h Action_Distance.h \
@@ -48,7 +50,8 @@
         Action_Radgyr.h Conflib.h Action_Mask.h Action_Closest.h \
         NetcdfRoutines.h AmberRestartNC.h CpptrajStdio.h Mol2File.h \
         Mol2FileRoutines.h Action_NAstruct.h DistRoutines.h AxisType.h \
-        TorsionRoutines.h Action_Pucker.h Range.h Action_Outtraj.h
+        TorsionRoutines.h Action_Pucker.h Range.h Action_Outtraj.h \
+        Action_Rms2d.h
 
 all: cpptraj$(SFX)
 
diff -urN amber11.orig/AmberTools/src/cpptraj/src/mapDataSet.cpp amber11/AmberTools/src/cpptraj/src/mapDataSet.cpp
--- amber11.orig/AmberTools/src/cpptraj/src/mapDataSet.cpp	2011-04-14 15:30:16.000000000 +0300
+++ amber11/AmberTools/src/cpptraj/src/mapDataSet.cpp	2011-10-25 15:01:28.104288645 +0300
@@ -6,6 +6,14 @@
 #include "CpptrajStdio.h"
 using namespace std;
 
+// CONSTRUCTOR
+mapDataSet::mapDataSet() {
+  width = 12;
+  precision = 4;
+  dType=DOUBLE;
+  setFormatString();
+}
+
 /*
  * mapDataSet::Xmax(()
  * Return the maximum X value added to this set. By convention this is 
@@ -34,8 +42,23 @@
   current++;
 }
 
-/*
- * mapDataSet::isEmpty()
+/* mapDataSet::Get()
+ * Get data at frame, put into vOut. Return 1 if no data at frame.
+ */
+int mapDataSet::Get(void *vOut, int frame) {
+  double *value;
+  
+  if (vOut==NULL) return 1;
+  //mprintf("DEBUG: Attempting to get double frame %i\n",frame);
+  value = (double*) vOut;
+  it=Data.find( frame );
+  if (it == Data.end()) return 1;
+  //mprintf("DEBUG: Double frame %i is %lf\n",frame,(*it).second);
+  *value = (*it).second;
+  return 0;
+}
+
+/* mapDataSet::isEmpty()
  */
 int mapDataSet::isEmpty(int frame) {
   it = Data.find( frame );
@@ -53,17 +76,17 @@
   it = Data.find( frame );
   if (it == Data.end()) 
     //sprintf(buffer," %12s","NoData");
-    sprintf(buffer," %12.4lf", 0.0);
+    sprintf(buffer, format, 0.0);
   else 
-    sprintf(buffer," %12.4lf",(*it).second);
-  return (buffer + 13);
+    sprintf(buffer, format,(*it).second);
+  return (buffer + width + 1);
 }
 
 /*
  * mapDataSet::Width()
  */
 int mapDataSet::Width() {
-  return 13;
+  return (width + 1);
 }
 
 /*
diff -urN amber11.orig/AmberTools/src/cpptraj/src/mapDataSet.h amber11/AmberTools/src/cpptraj/src/mapDataSet.h
--- amber11.orig/AmberTools/src/cpptraj/src/mapDataSet.h	2011-04-14 15:30:16.000000000 +0300
+++ amber11/AmberTools/src/cpptraj/src/mapDataSet.h	2011-10-25 15:01:28.105288645 +0300
@@ -7,14 +7,16 @@
  */
 #include <map>
 #include "DataSet.h"
-//using namespace std;
 class mapDataSet : public DataSet {
     std::map<int,double> Data;
     std::map<int,double>::iterator it;
   public:
+    mapDataSet();
+
     int Xmax();
     int isEmpty(int);
     void Add( int, void * );
+    int Get(void *, int);
     char *Write(char *, int);
     int Width();
     int Sync();
diff -urN amber11.orig/AmberTools/src/cpptraj/src/Mol2File.cpp amber11/AmberTools/src/cpptraj/src/Mol2File.cpp
--- amber11.orig/AmberTools/src/cpptraj/src/Mol2File.cpp	2011-04-14 15:30:16.000000000 +0300
+++ amber11/AmberTools/src/cpptraj/src/Mol2File.cpp	2011-10-25 15:01:28.127288624 +0300
@@ -239,6 +239,8 @@
   if (writeMode==2)
     File->IO->Close();
 
+  currentFrame++;
+
   return 0;
 }
  
diff -urN amber11.orig/AmberTools/src/cpptraj/src/PDBfile.cpp amber11/AmberTools/src/cpptraj/src/PDBfile.cpp
--- amber11.orig/AmberTools/src/cpptraj/src/PDBfile.cpp	2011-04-14 15:30:16.000000000 +0300
+++ amber11/AmberTools/src/cpptraj/src/PDBfile.cpp	2011-10-25 15:01:28.128288623 +0300
@@ -207,6 +207,9 @@
   } else if (writeMode==1) {
     File->IO->Printf("ENDMDL\n");
   }
+
+  currentFrame++;
+
   return 0;
 }
 
diff -urN amber11.orig/AmberTools/src/cpptraj/src/PtrajActionList.cpp amber11/AmberTools/src/cpptraj/src/PtrajActionList.cpp
--- amber11.orig/AmberTools/src/cpptraj/src/PtrajActionList.cpp	2011-04-14 15:30:16.000000000 +0300
+++ amber11/AmberTools/src/cpptraj/src/PtrajActionList.cpp	2011-10-25 15:01:28.105288645 +0300
@@ -20,6 +20,7 @@
 #include "Action_NAstruct.h"
 #include "Action_Pucker.h"
 #include "Action_Outtraj.h"
+#include "Action_Rms2d.h"
 
 // Constructor
 PtrajActionList::PtrajActionList() {
@@ -56,6 +57,7 @@
 
   // Decide what action this is based on the command.
   if      (A->CommandIs("distance")) {Act=new Distance;}
+  else if (A->CommandIs("rms2d"))    {Act=new Rms2d;   }
   else if (A->CommandIs("rmsd",3))   {Act=new Rmsd;    }
   else if (A->CommandIs("dihedral")) {Act=new Dihedral;}
   else if (A->CommandIs("atommap"))  {Act=new AtomMap; }
@@ -166,6 +168,8 @@
     err = ActionList[act]->DoAction(FrameAddress, frameIn);
     if (err==1) {
       // Treat actions that fail as if they could not be set up
+      mprintf("Warning: Action [%s] failed, frame %i.\n",ActionList[act]->CmdLine(),
+              frameIn);
       ActionList[act]->noSetup=1;
     } else if (err==2) {
       // Return value of 2 requests return to original frame
diff -urN amber11.orig/AmberTools/src/cpptraj/src/PtrajFile.cpp amber11/AmberTools/src/cpptraj/src/PtrajFile.cpp
--- amber11.orig/AmberTools/src/cpptraj/src/PtrajFile.cpp	2011-04-14 15:30:16.000000000 +0300
+++ amber11/AmberTools/src/cpptraj/src/PtrajFile.cpp	2011-10-25 15:01:28.196288560 +0300
@@ -113,6 +113,7 @@
 void PtrajFile::CloseFile() {
   if (isOpen) {
     IO->Close();
+    if (debug>0) rprintf("Closed %s.\n",filename);
     isOpen=0;
   }
 }
@@ -417,8 +418,8 @@
   // Standard file size is in the frame_stat struct
   uncompressed_size = IO->Size(filename);
 
-  // Determine format
-  // Read first 3 bytes again to determine format by magic number
+  // ========== Determine format ==========
+  // ---------- Read first 3 bytes again to determine format by magic number ----------
   IO->Open(filename,"rb"); // NOTE: Err Check
   memset(magic,0,3*sizeof(unsigned char));
   IO->Read(magic  ,1,1);
@@ -454,7 +455,7 @@
     return 0;
   }
 
-  // ID by file characteristics; read the first two lines
+  // ---------- ID by file characteristics; read the first two lines ----------
   // Initialize buffers to NULL
   buffer1[0]='\0';
   buffer2[0]='\0';
@@ -472,17 +473,6 @@
     }
   }
 
-  // Reopen and scan for Tripos mol2 molecule section
-  // 0 indicates section found.
-  IO->Open(filename,"r");
-  if (!Mol2ScanTo(this, MOLECULE)) {
-    if (debug>0) mprintf("  TRIPOS MOL2 file\n");
-    fileFormat=MOL2FILE;
-    IO->Close();
-    return 0;
-  }
-  IO->Close();
-
   // If both lines have PDB keywords, assume PDB
   if (isPDBkeyword(buffer1) && isPDBkeyword(buffer2)) {
     if (debug>0) mprintf("  PDB file\n");
@@ -491,7 +481,8 @@
   }
 
   // If either buffer contains a TRIPOS keyword assume Mol2
-  // NOTE: This will fail on tripos files with extensive header comments
+  // NOTE: This will fail on tripos files with extensive header comments.
+  //       A more expensive check for mol2 files is below.
   if (strncmp(buffer1,"@<TRIPOS>", 9)==0 ||
       strncmp(buffer2,"@<TRIPOS>", 9)==0) {
     if (debug>0) mprintf("  TRIPOS MOL2 file\n");
@@ -547,7 +538,19 @@
     }
   }
 
-  // NOTE: EXPERIMENTAL
+  // ---------- MORE EXPENSIVE CHECKS ----------
+  // Reopen and scan for Tripos mol2 molecule section
+  // 0 indicates section found.
+  IO->Open(filename,"r");
+  if (!Mol2ScanTo(this, MOLECULE)) {
+    if (debug>0) mprintf("  TRIPOS MOL2 file\n");
+    fileFormat=MOL2FILE;
+    IO->Close();
+    return 0;
+  }
+  IO->Close();
+
+  // ---------- EXPERIMENTAL ----------
   // If the file format is still undetermined and the file name is conflib.dat,
   // assume this is a conflib.dat file from LMOD. Cant think of a better way to
   // detect this since there is no magic number but the file is binary.
diff -urN amber11.orig/AmberTools/src/cpptraj/src/PtrajMpi.c amber11/AmberTools/src/cpptraj/src/PtrajMpi.c
--- amber11.orig/AmberTools/src/cpptraj/src/PtrajMpi.c	2011-04-14 15:30:16.000000000 +0300
+++ amber11/AmberTools/src/cpptraj/src/PtrajMpi.c	2011-10-25 15:01:28.100288649 +0300
@@ -460,7 +460,10 @@
   }
 
   err = MPI_Allreduce(input, Return, count, currentType, currentOp, MPI_COMM_WORLD);
-  if (err!=MPI_SUCCESS) printMPIerr(err, "Performing allreduce for %i elements\n",count);
+  if (err!=MPI_SUCCESS) {
+    printMPIerr(err, "Performing allreduce.\n");
+    rprintf("Error: allreduce failed for %i elements.\n",count);
+  }
 
   if (parallel_check_error(err)!=0) return 1;
   return 0;
diff -urN amber11.orig/AmberTools/src/cpptraj/src/PtrajState.cpp amber11/AmberTools/src/cpptraj/src/PtrajState.cpp
--- amber11.orig/AmberTools/src/cpptraj/src/PtrajState.cpp	2011-04-14 15:30:16.000000000 +0300
+++ amber11/AmberTools/src/cpptraj/src/PtrajState.cpp	2011-10-25 15:01:28.214288545 +0300
@@ -7,7 +7,8 @@
 #include "CpptrajStdio.h"
 
 // Constructor
-PtrajState::PtrajState() { 
+PtrajState::PtrajState() {
+  TotalErrors=0; 
   debug=0;
   showProgress=1;
 }
@@ -291,6 +292,7 @@
   char *df_cmd = NULL;
   char *name1 = NULL;
   char *name2 = NULL;
+  int width,precision;
   DataFile *df;
 
   if (DF_Args.empty()) return;
@@ -349,6 +351,19 @@
       }
       mprintf("    Not printing x column for datafile %s\n",name1);
       df->SetNoXcol();
+    
+    // datafile precision
+    // Usage: datafile precision <filename> <dataset> [<width>] [<precision>]
+    //        If width/precision not specified default to 12.4
+    } else if ( strcmp(df_cmd,"precision")==0 ) {
+      if (df==NULL) {
+        mprintf("Error: datafile precision: DataFile %s does not exist.\n",name1);
+        continue;
+      }
+      name2 = A->getNextString();
+      width = A->getNextInteger(12);
+      precision = A->getNextInteger(4);
+      df->SetPrecision(name2,width,precision);
     }
 
   } // END loop over datafile args
@@ -374,7 +389,7 @@
 
   // ========== S E T U P   P H A S E ========== 
   // Calculate frame division among trajectories
-  maxFrames=trajFileList.SetupFrames();
+  maxFrames=trajFileList.SetupFrames(worldrank,worldsize);
 
   // Parameter file information
   parmFileList.Print();
@@ -442,7 +457,7 @@
       // Perform Actions on Frame
       ptrajActionList.DoActions(&CurrentFrame, actionSet);
       // Do Output
-      outFileList.Write(outputSet, CurrentFrame, CurrentParm);
+      outFileList.Write(actionSet, CurrentFrame, CurrentParm);
 #ifdef DEBUG
       dbgprintf("\tDEBUG: %30s: %4i\n",CurrentParm->parmName,CurrentParm->outFrame);
 #endif
diff -urN amber11.orig/AmberTools/src/cpptraj/src/PtrajState.h amber11/AmberTools/src/cpptraj/src/PtrajState.h
--- amber11.orig/AmberTools/src/cpptraj/src/PtrajState.h	2011-04-14 15:30:16.000000000 +0300
+++ amber11/AmberTools/src/cpptraj/src/PtrajState.h	2011-10-25 15:01:28.093288656 +0300
@@ -20,6 +20,7 @@
     DataFileList DFL;                // List of datafiles that data sets will be written to
     ArgList *A;                      // Current argument list
     std::list<ArgList*> DF_Args;     // List of commands pertaining to datafile creation etc
+    int TotalErrors;                 // Sum of all returned error statuses
     int debug;
 
     void SetGlobalDebug(int);        // Set debug level for all components
diff -urN amber11.orig/AmberTools/src/cpptraj/src/ReferenceList.cpp amber11/AmberTools/src/cpptraj/src/ReferenceList.cpp
--- amber11.orig/AmberTools/src/cpptraj/src/ReferenceList.cpp	2011-04-14 15:30:16.000000000 +0300
+++ amber11/AmberTools/src/cpptraj/src/ReferenceList.cpp	2011-10-25 15:01:28.093288656 +0300
@@ -20,11 +20,15 @@
  */
 int ReferenceList::Add(ArgList *A, ParmFileList *parmFileList, int worldsize) {
   TrajFile *T;
-  int startArg;
+  int startArg,stopArg,offsetArg;
+  bool average = false;
 
   // Set up common arguments from arglist
   if (this->ProcessArgList(A,parmFileList)) return 1;
 
+  // Check if we want to obtain the average structure
+  average = A->hasKey("average");
+
   // Set up basic file to determine type and format
   T = this->SetupTrajectory(trajfilename, fileAccess, UNKNOWN_FORMAT, UNKNOWN_TYPE);
 
@@ -45,10 +49,18 @@
   // Get user-specified start arg
   // NOTE: For compatibility with ptraj start from 1
   startArg=A->getNextInteger(1);
-  T->SetArgs(startArg,startArg,1);
+  stopArg=startArg;
+  offsetArg=1;
+  // Get user-specified stop and offset only if getting avg structure
+  if (average) {
+    stopArg=A->getNextInteger(-1);
+    offsetArg=A->getNextInteger(1);
+  }
+  T->SetArgs(startArg,stopArg,offsetArg);
 
   // Add to trajectory file list
-  this->push_back(T); 
+  this->push_back(T);
+  Average.push_back(average); 
 
   return 0;
 }
@@ -60,9 +72,11 @@
  * place that frame in refFrames.
  */
 int ReferenceList::SetupRefFrames(FrameList *refFrames) {
-  int trajFrames;
-  Frame *F;
+  int trajFrames, global_set;
+  double Nframes;
+  Frame *F, *AvgFrame;
   int skipValue;
+  int refTrajNum = 0;
 
   mprintf("\nREFERENCE COORDS:\n");
   if (this->empty()) {
@@ -72,7 +86,13 @@
 
   for (it = this->begin(); it != this->end(); it++) {
     // Setup the reference traj for reading. Should only be 1 frame.
-    trajFrames=(*it)->setupFrameInfo(-1);
+    // NOTE: For MPI, calling setupFrameInfo with worldrank 0, worldsize 1 for 
+    //       all ranks. This is to ensure each thread has a copy of the ref 
+    //       struct.
+    //       Calling setupFrameInfo with -1 to ensure the Parm frame count is
+    //       not updated.
+
+    trajFrames=(*it)->setupFrameInfo(-1,0,1);
     if ((*it)->total_read_frames<1) {
       rprintf("Error: No frames could be read for reference %s, skipping\n",
       (*it)->trajfilename);
@@ -89,12 +109,34 @@
       skipValue=(*it)->skip;
       (*it)->skip=0;
     }
-    (*it)->Begin(&trajFrames, 0);
-    // Get and copy the 1 frame from Traj, then close
-    // NOTE: What happens when not seekable?
+    // Start trajectory read
+    global_set=0;
+    (*it)->Begin(&global_set, 0);
     (*it)->PrintInfo(1);
-    (*it)->NextFrame(&trajFrames);
-    F=(*it)->F->Copy();
+    // If averaging requested, loop over specified frames and avg coords.
+    if (Average[refTrajNum++]) {
+      mprintf("    Averaging over %i frames.\n",trajFrames);
+      AvgFrame = new Frame((*it)->P->natom, (*it)->P->mass);
+      AvgFrame->ZeroCoords();
+      global_set = 0;
+      Nframes = 0.0;
+      while ( (*it)->NextFrame(&global_set) ) {
+        AvgFrame->AddCoord( (*it)->F );
+        Nframes++;
+      }
+      if (Nframes < 1.0) { 
+        mprintf("Error: reference average: # frames read is less than 1.\n");
+        F=NULL;
+      } else {
+        AvgFrame->Divide( Nframes );
+        F=AvgFrame->Copy();
+      }
+      delete AvgFrame; 
+    // If no averaging, get and copy the 1 frame from Traj, then close
+    } else {
+      (*it)->NextFrame(&trajFrames);
+      F=(*it)->F->Copy();
+    }
     // DEBUG
     //fprintf(stdout,"DEBUG: Ref Coord Atom 0\n");
     //F->printAtomCoord(0);
diff -urN amber11.orig/AmberTools/src/cpptraj/src/ReferenceList.h amber11/AmberTools/src/cpptraj/src/ReferenceList.h
--- amber11.orig/AmberTools/src/cpptraj/src/ReferenceList.h	2011-04-14 15:30:16.000000000 +0300
+++ amber11/AmberTools/src/cpptraj/src/ReferenceList.h	2011-10-25 15:01:28.094288655 +0300
@@ -5,7 +5,7 @@
 #include "FrameList.h"
 
 class ReferenceList : public CoordFileList {
-
+    std::vector<bool> Average;
   public:
     
     ReferenceList();
diff -urN amber11.orig/AmberTools/src/cpptraj/src/RemdTraj.cpp amber11/AmberTools/src/cpptraj/src/RemdTraj.cpp
--- amber11.orig/AmberTools/src/cpptraj/src/RemdTraj.cpp	2011-04-14 15:30:16.000000000 +0300
+++ amber11/AmberTools/src/cpptraj/src/RemdTraj.cpp	2011-10-25 15:01:28.122288629 +0300
@@ -88,7 +88,7 @@
   // NOTE: Should check that this is the case for ALL frames.
   stop = T->Frames;
   Frames = T->Frames;
-  isBox = T->isBox;
+  BoxType = T->BoxType;
   trajfilename = T->File->basefilename;
   // Add it to the list
   REMDtraj.push_back(T);
@@ -190,12 +190,12 @@
       return 1;
     }
     // Check that #Frames and box info matches
-    if ( Frames!=T->Frames || isBox!=T->isBox ) {
+    if ( Frames!=T->Frames || BoxType!=T->BoxType ) {
       mprintf(
-              "    ERROR: REMDTRAJ: #Frames (%i) or box info (%i) in replica does not match\n",
-              T->Frames, T->isBox);
-      mprintf("                     values in lowest replica (Frames=%i, box=%i)\n",
-              Frames,isBox);
+              "    ERROR: REMDTRAJ: #Frames (%i) or box type (%i) in replica does not match\n",
+              T->Frames, T->BoxType);
+      mprintf("                     values in lowest replica (Frames=%i, boxtype=%i)\n",
+              Frames,BoxType);
       delete T;
       free(repFilename);
       free(Prefix);
diff -urN amber11.orig/AmberTools/src/cpptraj/src/stringDataSet.cpp amber11/AmberTools/src/cpptraj/src/stringDataSet.cpp
--- amber11.orig/AmberTools/src/cpptraj/src/stringDataSet.cpp	2011-04-14 15:30:16.000000000 +0300
+++ amber11/AmberTools/src/cpptraj/src/stringDataSet.cpp	2011-10-25 15:01:28.097288652 +0300
@@ -6,6 +6,13 @@
 #include "PtrajMpi.h"
 #include "CpptrajStdio.h"
 using namespace std;
+
+// CONSTRUCTOR
+stringDataSet::stringDataSet() {
+  dType=STRING;
+  setFormatString();
+}
+
 /*
  * stringDataSet::Xmax(()
  * Return the maximum X value added to this set. By convention this is 
@@ -54,10 +61,10 @@
 
   it = Data.find( frame );
   if (it == Data.end()) { 
-    sprintf(buffer," %s", "NoData");
+    sprintf(buffer, format, "NoData");
     return (buffer + 7);
   } else 
-    sprintf(buffer," %s",(*it).second.c_str());
+    sprintf(buffer, format, (*it).second.c_str());
 
   return (buffer + (*it).second.size() + 1);
 }
diff -urN amber11.orig/AmberTools/src/cpptraj/src/stringDataSet.h amber11/AmberTools/src/cpptraj/src/stringDataSet.h
--- amber11.orig/AmberTools/src/cpptraj/src/stringDataSet.h	2011-04-14 15:30:16.000000000 +0300
+++ amber11/AmberTools/src/cpptraj/src/stringDataSet.h	2011-10-25 15:01:28.097288652 +0300
@@ -12,6 +12,8 @@
     std::map<int,std::string> Data;
     std::map<int,std::string>::iterator it;
   public:
+    stringDataSet();
+
     int Xmax();
     int isEmpty(int);
     void Add( int, void * );
diff -urN amber11.orig/AmberTools/src/cpptraj/src/TrajFile.cpp amber11/AmberTools/src/cpptraj/src/TrajFile.cpp
--- amber11.orig/AmberTools/src/cpptraj/src/TrajFile.cpp	2011-04-14 15:30:16.000000000 +0300
+++ amber11/AmberTools/src/cpptraj/src/TrajFile.cpp	2011-10-25 15:01:28.122288629 +0300
@@ -2,7 +2,6 @@
 #include <cstdlib>
 #include <cstring>
 #include "TrajFile.h"
-#include "PtrajMpi.h" // worldrank and worldsize needed to calc frame division
 #include "CpptrajStdio.h"
 
 // CONSTRUCTOR
@@ -18,7 +17,7 @@
   start=0;
   stop=-1;
   offset=1;
-  isBox=0;
+  BoxType=0;
   title=NULL;
   P=NULL;
   frameskip=0;
@@ -45,12 +44,17 @@
 void TrajFile::SetTitle(char *titleIn) {
   size_t titleSize;
 
+  //mprintf("DEBUG: Attempting to set title for %s: [%s]\n",trajfilename,titleIn);
   if (titleIn==NULL) return;
   titleSize = strlen(titleIn);
-  if (titleSize==0) return;
+  //mprintf("       Title size is %i\n",titleSize);
+  if (titleSize==0) {
+    mprintf("Warning: TrajFile::SetTitle(): Title for %s is 0 length.\n",trajfilename);
+    return;
+  }
   this->title = (char*) malloc( (titleSize+1) * sizeof(char));
   if (this->title==NULL) {
-    mprintf("Error: TrajFile::SetTitle(): Could not allocate memory for title.\n");
+    mprintf("Error: SetTitle: Could not allocate memory for title of %s.\n",trajfilename);
     return;
   }
   strcpy(this->title, titleIn);
@@ -59,6 +63,30 @@
 }
 
 /*
+ * TrajFile::CheckBoxType()
+ * Set the trajectory box type (ortho/nonortho) based on box angles.
+ * Check the current box type against the associated parmfile box type.
+ * Print a warning if they are different.
+ */
+void TrajFile::CheckBoxType(double *box) {
+  // Determine orthogonal / non-orthogonal from angles
+  if (box[3]==0.0 || box[4]==0.0 || box[5]==0.0)
+    BoxType=0;
+  else if (box[3]==90.0 && box[4]==90.0 && box[5]==90.0)
+    BoxType=1;
+  else
+    BoxType=2;
+  if (P->BoxType != BoxType) {
+    mprintf("Warning: %s contains box info of type %i (beta %lf)\n",trajfilename,
+            BoxType,box[4]);
+    mprintf("         but associated parmfile %s has box type %i (beta %lf)\n",P->parmName, 
+            P->BoxType,P->Box[4]);
+    //mprintf("         Box information from trajectory will be used.\n");
+  }
+  if (debug>0) mprintf("    %s: Box type is %i (beta=%lf)\n",trajfilename,BoxType,box[4]);
+}
+
+/*
  * TrajFile::PrintInfo()
  * Print general trajectory information. Call TrajFile->Info for specific information.
  */
@@ -68,7 +96,7 @@
 
   mprintf(", Parm %i",P->pindex);
 
-  if (isBox) mprintf(" (with box info)");
+  if (BoxType>0) mprintf(" (with box info)");
 
   if (showExtended==0) {
     mprintf("\n");
@@ -85,7 +113,7 @@
     mprintf(": Writing %i frames", P->parmFrames);
     if (File->access==APPEND) mprintf(", appended"); // NOTE: Dangerous if REMD
   }
-  if (debug>0) mprintf(", %i atoms, Box %i, seekable %i",P->natom,isBox,seekable);
+  if (debug>0) mprintf(", %i atoms, Box %i, seekable %i",P->natom,BoxType,seekable);
   mprintf("\n");
 }
 
@@ -164,7 +192,7 @@
  * Note that the input frames start counting from 1, output starts counting from 0!
  * If called with maxFrames=-1 dont update the frame in parm file.
  */
-int TrajFile::setupFrameInfo(int maxFrames) {
+int TrajFile::setupFrameInfo(int maxFrames, int worldrank, int worldsize) {
   int Nframes;
   int ptraj_start_frame, ptraj_end_frame;
   int traj_start_frame, traj_end_frame;
diff -urN amber11.orig/AmberTools/src/cpptraj/src/TrajFile.h amber11/AmberTools/src/cpptraj/src/TrajFile.h
--- amber11.orig/AmberTools/src/cpptraj/src/TrajFile.h	2011-04-14 15:30:16.000000000 +0300
+++ amber11/AmberTools/src/cpptraj/src/TrajFile.h	2011-10-25 15:01:28.128288623 +0300
@@ -25,44 +25,44 @@
     int offset;         // Number of frames to skip while processing
 
     // --== Inherited by child classes ==--
-    virtual int open() { return 0; }        // Open the file, prepare for coord read/write
-    virtual void close() {}                 // Close the file
+    virtual int open() { return 0; } // Open the file, prepare for coord read/write
+    virtual void close() {}          // Close the file
+
+    void CheckBoxType(double *);     // Check that traj boxtype matches parm
 
   public:
     int debug;             // Level of debug information to print
     char *trajfilename;    // The base trajectory filename
-    // NOTE: I hate that the following are public. Only necessary for REMD processing!!
+    // NOTE: I hate that the 2 following are public. Only necessary for REMD processing!!
     int Frames;            // Total number of frames in trajectory
     int total_read_frames; // Total number of frames that will be read
-    int isBox;             // >0 means trajectory has box information
-
+    int BoxType;           // 0: None, 1: Ortho, 2: NonOrtho 
     Range *FrameRange;     // list of frames to be written out
     int hasTemperature;    // 1 means trajectory has temperature information
     PtrajFile *File;       // Class that handles basic file IO
     AmberParm *P;          // Memory address of the associated parmfile
     Frame *F;              // Hold coordinates of the current frame
-    int skip;              /* READ: If =1 do not process this input trajectory 
-                              WRITE: If =1 this traj has been set up for write */
+    int skip;              // READ: If =1 do not process this input trajectory 
+                           // WRITE: If =1 this traj has been set up for write
 
     TrajFile();            // Constructor
     virtual ~TrajFile();   // Destructor - virtual since this class is inherited.
 
-    int Start() { return start; }
+    int Start()        { return start;        }
+    int CurrentFrame() { return currentFrame; }
     void SetTitle(char *);   // Set trajectory title.
     void PrintInfo(int);     // Print trajectory Information
-    int setupFrameInfo(int); // Set actual start/stop based on total #frames and #threads 
-    int Begin(int *, int);   /* Prepare traj for processing. Set output start value, calcd in 
-                              * setupFrameInfo. Allocate memory for F. 
-                              */
-    int Begin();                 // Prepare trajectory for output
-    int NextFrame(int*);         // Put the next target frame into F.
-    void End();                  // Close trajectory and free F memory
-    void progressBar();          // Display trajectory progress to screen
-//    void progressBar2();         // Display trajectory progress to screen
-   
-    void SetArgs(int,int,int);   // Set the stop, start, and offset args from user input
+    int Begin(int *, int);   // Prepare traj for processing. Set output start value, calcd in 
+                             // setupFrameInfo. Allocate memory for F. 
+    int Begin();             // Prepare trajectory for output
+    int NextFrame(int*);     // Put the next target frame into F.
+    void End();              // Close trajectory and free F memory
+    void progressBar();      // Display trajectory progress to screen
+//    void progressBar2();   // Display trajectory progress to screen
+    int setupFrameInfo(int,int,int); // Set actual start/stop based on total #frames and #threads 
+    void SetArgs(int,int,int);       // Set the stop, start, and offset args from user input
     // --== Inherited by child classes ==--
-    virtual int getFrame(int)      { return 1; } // Read the next coord frame into F
+    virtual int getFrame(int)      { return 1; } // Read specified frame into F
     virtual int SetupRead()        { return 1; } // Set file up for reading
     virtual int WriteArgs(ArgList*){ return 0; } // (Opt.) Process any args related to writing
     virtual int SetupWrite()       { return 1; } // Set file up for writing
diff -urN amber11.orig/AmberTools/src/cpptraj/src/TrajinList.cpp amber11/AmberTools/src/cpptraj/src/TrajinList.cpp
--- amber11.orig/AmberTools/src/cpptraj/src/TrajinList.cpp	2011-04-14 15:30:16.000000000 +0300
+++ amber11/AmberTools/src/cpptraj/src/TrajinList.cpp	2011-10-25 15:01:28.201288557 +0300
@@ -1,4 +1,5 @@
 // TrajinList
+#include <cstddef>
 #include "TrajinList.h"
 #include "RemdTraj.h"
 #include "CpptrajStdio.h"
@@ -70,7 +71,7 @@
  * actual start and stop and how many frames total will be processed.
  * Return the number of frames to be processed.
  */
-int TrajinList::SetupFrames() {
+int TrajinList::SetupFrames(int worldrank, int worldsize) {
   int maxFrames, trajFrames;
 
   mprintf("\nTRAJECTORIES:\n");
@@ -78,7 +79,7 @@
   maxFrames=0;
 
   for (it = this->begin(); it != this->end(); it++) {
-    trajFrames = (*it)->setupFrameInfo(maxFrames);
+    trajFrames = (*it)->setupFrameInfo(maxFrames,worldrank,worldsize);
     if (trajFrames==-1) {
       maxFrames=-1;
     }
diff -urN amber11.orig/AmberTools/src/cpptraj/src/TrajinList.h amber11/AmberTools/src/cpptraj/src/TrajinList.h
--- amber11.orig/AmberTools/src/cpptraj/src/TrajinList.h	2011-04-14 15:30:16.000000000 +0300
+++ amber11/AmberTools/src/cpptraj/src/TrajinList.h	2011-10-25 15:01:28.090288658 +0300
@@ -13,7 +13,7 @@
     // NOTE: worldsize is passed in as last arg to avoid include of PtrajMpi
     int Add(ArgList *A, ParmFileList *, int);
     // TRAJIN: Set up frames to be processed 
-    int SetupFrames();
+    int SetupFrames(int,int);
 };
 #endif
 
diff -urN amber11.orig/AmberTools/src/cpptraj/src/TrajoutList.cpp amber11/AmberTools/src/cpptraj/src/TrajoutList.cpp
--- amber11.orig/AmberTools/src/cpptraj/src/TrajoutList.cpp	2011-04-14 15:30:16.000000000 +0300
+++ amber11/AmberTools/src/cpptraj/src/TrajoutList.cpp	2011-10-25 15:01:28.202288556 +0300
@@ -1,4 +1,5 @@
 // TrajoutList
+#include <cstddef>
 #include "TrajoutList.h"
 #include "CpptrajStdio.h"
 
@@ -20,7 +21,6 @@
  */
 int TrajoutList::Add(ArgList *A, ParmFileList *parmFileList, int worldsize) {
   TrajFile *T;
-  int boxInfo;
   FileFormat writeFormat;
   FileType writeType;
   char *onlyframes; 
@@ -29,7 +29,6 @@
   if (this->ProcessArgList(A,parmFileList)) return 1;
 
   // Init variables
-  boxInfo=0;
   writeFormat=AMBERTRAJ; 
   writeType=UNKNOWN_TYPE;
  
@@ -56,10 +55,6 @@
   if (worldsize>1 && writeFormat!=AMBERRESTART) 
     writeType=MPIFILE;
 
-  // Set box info from parm file unless nobox is set.
-  boxInfo=P->ifbox;
-  if (A->hasKey("nobox")) boxInfo=0;
-
   // Set up basic file for given type and format
   // If type is unknown it will be determined from extension or will be standard (default)
   T = this->SetupTrajectory(trajfilename, fileAccess, writeFormat, writeType);
@@ -90,10 +85,11 @@
   // Set parameter file
   T->P=P;
 
-  // Set box information (only needed for write)
-  T->isBox=boxInfo;
+  // Set box type from parm file unless "nobox" specified  
+  T->BoxType=P->BoxType;
+  if (A->hasKey("nobox")) T->BoxType=0;
 
-  // No setup here; Write is set up after first frame read in PtrajState::Run
+  // No more setup here; Write is set up when first frame written.
   // Add to trajectory file list
   this->push_back(T); 
 
diff -urN amber11.orig/AmberTools/src/etc/chemistry/amber/readparm.py amber11/AmberTools/src/etc/chemistry/amber/readparm.py
--- amber11.orig/AmberTools/src/etc/chemistry/amber/readparm.py	2011-04-14 15:30:17.000000000 +0300
+++ amber11/AmberTools/src/etc/chemistry/amber/readparm.py	2011-10-25 15:01:28.200288558 +0300
@@ -118,22 +118,44 @@
    """ Parses the fortran format statement. Recognizes ints, exponents, and strings.
        Returns the number of items/line, size of each item, and type of data """
 
+   # Get rid of ( and ) specifiers in Fortran format strings. This is a hack, but
+   # should work for existing chamber prmtop files
+
+   format_string = format_string.replace('(','').replace(')','')
+
+   # Fix case for E, I, and F
+
+   format_string = format_string.replace('e','E')
+   format_string = format_string.replace('i','I')
+   format_string = format_string.replace('f','F')
+
    if 'a' in format_string: # this is a string
       format_parts = format_string.split('a')
-      return int(format_parts[0]), int(format_parts[1]), 'str'
+      try:
+         return int(format_parts[0]), int(format_parts[1]), 'str', None
+      except:
+         return 1, 80, 'str', None
 
    elif 'I' in format_string: # this is an integer
       format_parts = format_string.split('I')
-      return int(format_parts[0]), int(format_parts[1]), 'int'
+      if len(format_parts[0].strip()) == 0: format_parts[0] = 1
+      return int(format_parts[0]), int(format_parts[1]), 'int', None
 
    elif 'E' in format_string: # this is a floating point decimal
       format_parts = format_string.split('E')
       decimal_parts = format_parts[1].split('.')
-      return int(format_parts[0]), int(decimal_parts[0]), 'dec'
+      if len(format_parts[0].strip()) == 0: format_parts[0] = 1
+      return int(format_parts[0]), int(decimal_parts[0]), 'dec', int(decimal_parts[1])
+   
+   elif 'F' in format_string: # this is also a floating point decimal
+      format_parts = format_string.split('F')
+      decimal_parts = format_parts[1].split('.')
+      if len(format_parts[0].strip()) == 0: format_parts[0] = 1
+      return int(format_parts[0]), int(decimal_parts[0]), 'dec', int(decimal_parts[1])
 
    else:
       print >> stderr, 'Error: Unrecognized format "%s"!' % format_string
-      return 1, 80, 'str'
+      return 1, 80, 'str', None
 
 # +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
@@ -275,7 +297,7 @@
 
          elif prmlines[i][0:7] == '%FORMAT':
             self.formats[current_flag] = prmlines[i][8:len(prmlines[i].strip())-1]
-            number_items_perline, size_item, dat_type = _parseFormat(self.formats[current_flag])
+            number_items_perline, size_item, dat_type, junk = _parseFormat(self.formats[current_flag])
             gathering_data = True
 
          elif gathering_data:
@@ -339,17 +361,17 @@
          flag = self.flag_list[i]
          new_prm.write('%%FLAG %s\n' % flag)
          new_prm.write('%%FORMAT(%s)\n' % self.formats[flag])
-         number_items_perline, size_item, dat_type = _parseFormat(self.formats[flag])
-         if dat_type == 'dec':
-            decnum = int(self.formats[flag].split('E')[1].split('.')[1])
+         number_items_perline, size_item, dat_type, decnum = _parseFormat(self.formats[flag])
          line = ''
          num_items = 0
          if len(self.parm_data[flag]) == 0: # empty field...
             new_prm.write('\n')
             continue
          for j in range(len(self.parm_data[flag])): # write data in new_prm
-            if dat_type == 'dec':
+            if dat_type == 'dec' and 'E' in self.formats[flag].upper():
                line += ('%%%s.%sE' % (size_item, decnum)) % self.parm_data[flag][j] 
+            elif dat_type == 'dec' and 'F' in self.formats[flag].upper():
+               line += ('%%%s.%sF' % (size_item, decnum)) % self.parm_data[flag][j] 
             elif dat_type == 'int':
                line += ('%%%sd' % size_item) % self.parm_data[flag][j] 
             else:
@@ -764,7 +786,8 @@
          self.LJ_types[self.parm_data["AMBER_ATOM_TYPE"][i]] = self.parm_data["ATOM_TYPE_INDEX"][i]
          
       for i in range(self.pointers["NTYPES"]):
-         lj_index = (i + 1) * (i + 2) / 2 - 1 # n(n+1)/2 th position adjusted for indexing from 0
+         lj_index = self.parm_data["NONBONDED_PARM_INDEX"][
+                     self.pointers["NTYPES"] * i + i - 1] - 1
          if self.parm_data["LENNARD_JONES_BCOEF"][lj_index] < 1.0e-6:
             self.LJ_radius.append(0)
             self.LJ_depth.append(0)
diff -urN amber11.orig/AmberTools/src/leap/src/leap/amber.c amber11/AmberTools/src/leap/src/leap/amber.c
--- amber11.orig/AmberTools/src/leap/src/leap/amber.c	2011-03-10 20:12:06.000000000 +0200
+++ amber11/AmberTools/src/leap/src/leap/amber.c	2011-10-25 15:01:28.139288613 +0300
@@ -36,6 +36,10 @@
  *              UNITs and PARMSETs.
  */
 
+        /*      Arunima Singh (UGA)     */
+        /*         20110420             */
+        /*   	 Added SCEE/SCNB   	*/
+
 #include        "basics.h"
 #include        "vector.h"
 #include        "classes.h"
@@ -240,23 +244,47 @@
  *
  *      Read the proper torsion parameter terms.
  */
+
+/* Arunima Singh added reading in of scee and scnb fields for proper torsions 20110420 */
 static void
 zAmberReadParmSetPropers( PARMSET psParms, FILE *fIn )
 {
 STRING          sLine;
-int             iRead, iN;
+int             iRead, iN, iRead2, iRead3;
 STRING          saStr[10];
-double          dDivisions, dKp, dP0, dN;
+double          dDivisions, dKp, dP0, dN, dScee, dScnb;
+char		*cScee, *cScnb;
 
     memset(saStr, 0, sizeof(saStr));                    /* for Purify */
     while (1) {
         FGETS( sLine, fIn );
         NODASHES(sLine);
-        iRead = sscanf( sLine, "%s %s %s %s %lf %lf %lf %lf", 
+
+	iRead = sscanf( sLine, "%s %s %s %s %lf %lf %lf %lf ", 
                                 saStr[0], saStr[1], saStr[2], saStr[3], 
                                 &dDivisions, &dKp, &dP0, &dN );
-        if ( iRead <= 0 ) 
+
+	cScee = strstr(sLine, "SCEE");
+	if(cScee!=NULL) {
+		iRead2 = sscanf( cScee, "SCEE=%lf", &dScee);
+		iRead++;
+		}
+
+        cScnb = strstr(sLine, "SCNB");	
+        if(cScnb!=NULL) {
+                iRead3 = sscanf( cScnb, "SCNB=%lf", &dScnb);
+                iRead++;
+                }
+
+	if ( iRead <= 0 ) 
                 break;
+
+	if ( iRead == 8 )	/*Arunima Singh*/
+	{
+		dScee = 1.2;
+		dScnb = 2.0;
+	}
+	
         MESSAGE(( "Read: %s\n", sLine ));
 
         if ( sLine[0] == ' ' && sLine[1] == ' ') {
@@ -277,14 +305,37 @@
         zAmberConvertWildCard( saStr[3] );
         iParmSetAddProperTerm( psParms, 
                                 saStr[0], saStr[1], saStr[2], saStr[3],
-                                abs(iN), dKp, dP0*DEGTORAD, "" );
+                                abs(iN), dKp, dP0*DEGTORAD, dScee, dScnb, "" );
+
+
         while( iN < 0 ) {
                 FGETS( sLine, fIn );
                 NODASHES(sLine);
                 MESSAGE(( "Read extra term: %s\n", sLine ));
-                iRead = sscanf( &sLine[11], "%lf %lf %lf %lf",
+                iRead = sscanf( &sLine[11], "%lf %lf %lf %lf ",
                                 &dDivisions, &dKp, &dP0, &dN );
-                if ( iRead<=0 ) break;
+               
+	        cScee = strstr(sLine, "SCEE");
+        	if(cScee!=NULL) {
+                	iRead2 = sscanf( cScee, "SCEE=%lf", &dScee);
+	                iRead++;
+        	        }
+
+	        cScnb = strstr(sLine, "SCNB");
+        	if(cScnb!=NULL) {
+	                iRead3 = sscanf( cScnb, "SCNB=%lf", &dScnb);
+        	        iRead++;
+                	}
+ 
+
+	        if ( iRead == 4 )       /*Arunima Singh*/
+	        {
+	                dScee = 1.2;
+                	dScnb = 2.0;
+        	}
+
+
+		if ( iRead<=0 ) break;
 
                 if ( dDivisions == 0.0 ) 
                         dDivisions = 1.0;
@@ -292,7 +343,7 @@
                 iN = (int)floor(dN+0.5);
                 iParmSetAddProperTerm( psParms,
                                 saStr[0], saStr[1], saStr[2], saStr[3],
-                                abs(iN), dKp, dP0*DEGTORAD, "" );
+                                abs(iN), dKp, dP0*DEGTORAD, dScee, dScnb, "" );
         }
         if ( iRead <= 0 )
                 break;
@@ -311,7 +362,7 @@
 STRING          sLine;
 int             iRead, iN;
 STRING          saStr[10];
-double          dKp, dP0, dN;
+double          dKp, dP0, dN, dScee, dScnb;
 BOOL            bPrintLine;
 
     memset(saStr, 0, sizeof(saStr));                    /* for Purify */
@@ -335,6 +386,8 @@
         zAmberConvertWildCard( saStr[2] );
         zAmberConvertWildCard( saStr[3] );
         iN = (int)dN;
+	dScee = 0.0;
+	dScnb = 0.0;
 
         /*
          *  check everything in case a format or other user error
@@ -362,7 +415,7 @@
 
         iParmSetAddImproperTerm( psParms, 
                                 saStr[0], saStr[1], saStr[2], saStr[3],
-                                iN, dKp, dP0*DEGTORAD, "" );
+                                iN, dKp, dP0*DEGTORAD, dScee, dScnb, "" );
     }
     if ( iRead > 0 )
         VP0(( "WARNING: incomplete Improper Torsion line:\n%s", sLine ));
diff -urN amber11.orig/AmberTools/src/leap/src/leap/build.c amber11/AmberTools/src/leap/src/leap/build.c
--- amber11.orig/AmberTools/src/leap/src/leap/build.c	2011-03-10 20:12:06.000000000 +0200
+++ amber11/AmberTools/src/leap/src/leap/build.c	2011-10-25 15:01:28.140288613 +0300
@@ -42,8 +42,9 @@
  *		
  */
 
-
-
+        /*      Arunima Singh (UGA)     */
+        /*         20110420             */
+        /*      Added SCEE/SCNB   	*/
 
 #include	"basics.h"
 
@@ -1567,7 +1568,7 @@
 LOOP            lAtoms, lTemp;
 ATOM            aAtom, aAtom1, aAtom2, aAtom3, aAtom4;
 BOOL            bM1, bM2, bM3, bM4, bOneMinimizedAtom;
-double          dKb, dR0, dKt, dT0, dTkub, dRkub, dKp, dP0;
+double          dKb, dR0, dKt, dT0, dTkub, dRkub, dKp, dP0, dScee, dScnb;
 STRING		sAtom1, sAtom2, sAtom3, sAtom4, sDesc;
 PARMSET		psTemp;
 TORSION		tTorsion;
@@ -1757,7 +1758,7 @@
 	    ParmSetTORSIONTerm( tTorsion, i,
 	    		&iIndex,
 			sAtom1, sAtom2, sAtom3, sAtom4,
-			&iN, &dKp, &dP0, sDesc );
+			&iN, &dKp, &dP0, &dScee, &dScnb, sDesc );
 	    if ( !bMinimizerAddTorsion( mStrain, 
 	    				aAtom1, aAtom2, aAtom3, aAtom4,
 					(double)iN, dKp, dP0 )) {
diff -urN amber11.orig/AmberTools/src/leap/src/leap/model.c amber11/AmberTools/src/leap/src/leap/model.c
--- amber11.orig/AmberTools/src/leap/src/leap/model.c	2011-04-14 15:30:17.000000000 +0300
+++ amber11/AmberTools/src/leap/src/leap/model.c	2011-10-25 15:01:28.141288612 +0300
@@ -39,7 +39,9 @@
  *              are determined from atom types.
  */
 
-
+        /*      Arunima Singh (UGA)     */
+        /*         20110420             */
+        /*      Added SCEE/SCNB         */
 
 
 #include	"basics.h"
@@ -781,16 +783,6 @@
 	    maPAtom++;
 	}
     }
-    if ( iAtomCoordination(aX) != 
-          1 + maPAtom - &(mtTorsions.maaXBonds[0]) ) {
-	VP0(( "Error: Atom %s has force field coordination %i\n"
-	      "       but only %i bonded neighbors.\n"
-	      "       The cause may be an incorrect atom type, and\n"
-	      "       the effect may be a crash very soon.\n",
-	      sContainerFullDescriptor((CONTAINER)aX,s1),
-	      iAtomCoordination(aX),
-	      1 + maPAtom - &(mtTorsions.maaXBonds[0]) ));
-    }
     for ( i=iAtomCoordination(aX); i<MAXBONDS; i++ ) {
 	maPAtom->aAtom = NULL;
 	maPAtom++;
@@ -824,16 +816,6 @@
 	    maPAtom++;
 	}
     }
-    if ( iAtomCoordination(aY) != 
-          1 + maPAtom - &(mtTorsions.maaYBonds[0]) ) {
-	VP0(( "Error: Atom %s has force field coordination %i\n"
-	      "       but only %i bonded neighbors.\n"
-	      "       The cause may be an incorrect atom type, and\n"
-	      "       the effect may be a crash very soon.\n",
-	      sContainerFullDescriptor((CONTAINER)aY,s1),
-	      iAtomCoordination(aY),
-	      1 + maPAtom - &(mtTorsions.maaYBonds[0]) ));
-    }
     for ( i=iAtomCoordination(aY); i<MAXBONDS; i++ ) {
 	maPAtom->aAtom = NULL;
 	maPAtom++;
@@ -1205,6 +1187,8 @@
 		int	iN;
 		double	dK;
 		double	dE;
+		double  dSce;	/*Arunima Singh*/
+		double  dScn;	/*Arunima Singh*/
 		} H_PROPERPARMt;
 
 typedef	struct	{
@@ -1222,12 +1206,13 @@
 
 	/* Keep iHybrid2 <= iHybrid3 */
 
+/* Arunima Singh 20110413. Added Scee and Scnb AMBER defaults of 1.2 and 2.0 */
 #define	TFORCE	20.0
 static	H_PROPERPARMt	SppaPropers[] = {
-{ 	HSP3,	HSP3,	3,	1.0,		0.0 },	/* Non bond */
-{ 	HSP2,	HSP3,	6,	-2.0,		0.0 },	/* Non bond */
-{	HSP2,	HSP2,	2,	-4.0,		0.0 },	/* Pi bond overlap */
-{	HSP1,	HSP1,	1,	0.0,		0.0 }	/* Not interesting */
+{ 	HSP3,	HSP3,	3,	1.0,		0.0,	1.2,	2.0 },	/* Non bond */
+{ 	HSP2,	HSP3,	6,	-2.0,		0.0,    1.2,    2.0 },	/* Non bond */
+{	HSP2,	HSP2,	2,	-4.0,		0.0,    1.2,    2.0 },	/* Pi bond overlap */
+{	HSP1,	HSP1,	1,	0.0,		0.0,    1.2,    2.0 }	/* Not interesting */
 };
 
 #define	AFORCE	100.0
@@ -1290,6 +1275,8 @@
 					SppaPropers[i].iN,
 					SppaPropers[i].dK,
 					SppaPropers[i].dE,
+					SppaPropers[i].dSce,
+					SppaPropers[i].dScn,
 					sDesc );
 	}
     }
diff -urN amber11.orig/AmberTools/src/leap/src/leap/parmSet.c amber11/AmberTools/src/leap/src/leap/parmSet.c
--- amber11.orig/AmberTools/src/leap/src/leap/parmSet.c	2011-03-10 20:12:06.000000000 +0200
+++ amber11/AmberTools/src/leap/src/leap/parmSet.c	2011-10-25 15:01:28.143288610 +0300
@@ -54,7 +54,9 @@
  *		search routines REQUIRE pre-ordering.
  */
 
-
+        /*      Arunima Singh (UGA)     */
+        /*         20110420             */
+        /*      Added SCEE/SCNB         */
 
 #include	"basics.h"
 
@@ -965,6 +967,8 @@
 	/*
 	 *  copy into 1 vararray & set iType
 	 */
+	
+	/*Arunima Singh added scee and scnb into the database */
 	vaTorsTypes = vaVarArrayCopy2( psLib->vaTorsions, psLib->vaImpropers );
 	tP = PVAI(vaTorsTypes,TORSIONPARMt,0);
 	for (i=0; i<iVarArrayElementCount(psLib->vaTorsions); i++, tP++)
@@ -980,8 +984,12 @@
                 7, "n",
 		    (char *)&(PVAI(vaTorsTypes,TORSIONPARMt,0)->iN),
 		    iVarArrayElementSize(vaTorsTypes),
-                0, NULL, NULL, 0,
-                0, NULL, NULL, 0,
+                9, "scee", 
+		    (char *)&(PVAI(vaTorsTypes,TORSIONPARMt,0)->dScee), 
+		    iVarArrayElementSize(vaTorsTypes),
+                10, "scnb", 
+		    (char *)&(PVAI(vaTorsTypes,TORSIONPARMt,0)->dScnb), 
+		    iVarArrayElementSize(vaTorsTypes),
                 0, NULL, NULL, 0,
                 0, NULL, NULL, 0,
                 0, NULL, NULL, 0,
@@ -1008,7 +1016,7 @@
                 4, "type4",
 		    (char *)&(PVAI(vaTorsTypes,TORSIONPARMt,0)->sType4),
                     iVarArrayElementSize(vaTorsTypes),
-                9, "desc",
+                11, "desc",
 		    (char *)&(PVAI(vaTorsTypes,TORSIONPARMt,0)->sDesc),
                     iVarArrayElementSize(vaTorsTypes)
                 );
@@ -1586,10 +1594,15 @@
  *      Add a torsion parameter to the PARMSET.
  *      Return the index. 
  */
+
+        /*      Arunima Singh (UGA)     */
+        /*         20110420             */
+        /*      Added SCEE/SCNB         */
+
 int
 iParmSetAddProperTerm( PARMSET psLib, 
 	char *sType1, char *sType2, char *sType3, char *sType4, 
-	int iN, double dKp, double dP0, char *sDesc )
+	int iN, double dKp, double dP0, double dScee, double dScnb, char *sDesc )
 {
 TORSIONPARMt    tpTorsion;
 
@@ -1605,6 +1618,8 @@
     tpTorsion.dKp = dKp;
     tpTorsion.iN  = iN;
     tpTorsion.dP0 = dP0;
+    tpTorsion.dScee = dScee;
+    tpTorsion.dScnb = dScnb;
     strcpy( tpTorsion.sOrder, "0123" );
     if ( sDesc != NULL )
     	strcpy( tpTorsion.sDesc, sDesc );
@@ -1633,7 +1648,7 @@
 int
 iParmSetAddImproperTerm( PARMSET psLib, 
 	char *sType1, char *sType2, char *sType3, char *sType4, 
-	int iN, double dKp, double dP0, char *sDesc )
+	int iN, double dKp, double dP0, double dScee, double dScnb, char *sDesc )
 {
 TORSIONPARMt    tpImproper;
 orderStr	sOrder;
@@ -2022,7 +2037,7 @@
 void
 ParmSetTORSIONTerm( TORSION tTorsion, int iTorsionIndex, int *iPParmSetIndex,
 		   char *cPType1, char *cPType2, char *cPType3, char *cPType4,
-		   int *iPN, double *dPKp, double *dPP0, char *sDesc )
+		   int *iPN, double *dPKp, double *dPP0, double *dPScee, double *dPScnb, char *sDesc )
 {
 TORSION_MATCHt	*tmPCur;
 
@@ -2035,6 +2050,8 @@
     *iPN = tmPCur->tpTorsion.iN;
     *dPKp = tmPCur->tpTorsion.dKp;
     *dPP0 = tmPCur->tpTorsion.dP0;
+    *dPScee = tmPCur->tpTorsion.dScee;
+    *dPScnb = tmPCur->tpTorsion.dScnb;
     strcpy(sDesc, tmPCur->tpTorsion.sDesc);
 }
 
@@ -2050,7 +2067,7 @@
 BOOL
 bParmSetTORSIONAddProperTerm( TORSION tTorsion,
 	char *cPType1, char *cPType2, char *cPType3, char *cPType4,
-	int iN, double dKp, double dP0, char *sDesc )
+	int iN, double dKp, double dP0, double dScee, double dScnb, char *sDesc )
 {
 TORSIONPARMt	tpTorsion;
 
@@ -2066,6 +2083,8 @@
     tpTorsion.iN = iN;
     tpTorsion.dKp = dKp;
     tpTorsion.dP0 = dP0;
+    tpTorsion.dScee = dScee;
+    tpTorsion.dScnb = dScnb;
     strcpy(tpTorsion.sDesc, sDesc);
     strcpy( tpTorsion.sOrder, "0123" );
     
@@ -2087,7 +2106,7 @@
 BOOL
 bParmSetTORSIONAddImproperTerm( TORSION tTorsion,
 	char *cPType1, char *cPType2, char *cPType3, char *cPType4,
-	int iN, double dKp, double dP0, char *sDesc )
+	int iN, double dKp, double dP0, double dScee, double dScnb, char *sDesc )
 {
 TORSIONPARMt	tpTorsion;
 orderStr	sOrder;
@@ -2104,6 +2123,8 @@
     tpTorsion.iN = iN;
     tpTorsion.dKp = dKp;
     tpTorsion.dP0 = dP0;
+    tpTorsion.dScee = dScee;
+    tpTorsion.dScnb = dScnb;
     strcpy(tpTorsion.sDesc, sDesc);
     strcpy( tpTorsion.sOrder, sOrder );
     
@@ -2516,7 +2537,7 @@
 void
 ParmSetTorsion( PARMSET psLib, int i, 
 	char *sType1, char *sType2, char *sType3, char *sType4,
-	int *iPN, double *dPKp, double *dPP0, char *sDesc)
+	int *iPN, double *dPKp, double *dPP0, double *dPScee, double *dPScnb, char *sDesc)
 {
 TORSIONPARMt   *tpPTorsion;
 
@@ -2531,6 +2552,8 @@
 	strcpy( sType4, WILD_CARD_TYPE );
 	*iPN  = 0;
 	*dPKp = 0;
+	*dPScee = 0;
+        *dPScnb = 0;
 	*dPP0 = 0;
 	strcpy( sDesc, "??" );
 	return;
@@ -2543,7 +2566,10 @@
     *iPN  = tpPTorsion->iN;
     *dPKp = tpPTorsion->dKp;
     *dPP0 = tpPTorsion->dP0;
+    *dPScee = tpPTorsion->dScee;
+    *dPScnb = tpPTorsion->dScnb; 
     strcpy( sDesc, tpPTorsion->sDesc );
+
 }
 
 
@@ -2557,7 +2583,7 @@
 void    
 ParmSetImproper( PARMSET psLib, int i, 
 	char *sType1, char *sType2, char *sType3, char *sType4,
-	int *iPN, double *dPKp, double *dPP0, char *sDesc)
+	int *iPN, double *dPKp, double *dPP0, double *dPScee, double *dPScnb, char *sDesc)
 {
 TORSIONPARMt   *tpPImproper;
 
@@ -2573,6 +2599,8 @@
 	*iPN  = 0;
 	*dPKp = 0;
 	*dPP0 = 0;
+	*dPScee = 0;
+	*dPScnb = 0;
 	strcpy( sDesc, "??" );
 	return;
     }
@@ -2584,6 +2612,8 @@
     *iPN  = tpPImproper->iN;
     *dPKp = tpPImproper->dKp;
     *dPP0 = tpPImproper->dP0;
+    *dPScee = tpPImproper->dScee;
+    *dPScnb = tpPImproper->dScnb;
     strcpy( sDesc, tpPImproper->sDesc );
 }
 
@@ -2740,7 +2770,7 @@
 void
 ParmSetUpdateTorsion( PARMSET psLib, int i, 
 	char *sType1, char *sType2, char *sType3, char *sType4,
-	int *iPN, double *dPKp, double *dPP0, char *sDescription)
+	int *iPN, double *dPKp, double *dPP0, double *dPScee, double *dPScnb, char *sDescription)
 {
 TORSIONPARMt	*tpPTorsion;
 orderStr	sOrder;
@@ -2754,6 +2784,8 @@
     if (         iPN != (int*)NULL   ) tpPTorsion->iN = *iPN;
     if (        dPKp != (double*)NULL) tpPTorsion->dKp = *dPKp;
     if (        dPP0 != (double*)NULL) tpPTorsion->dP0 = *dPP0;
+    if (      dPScee != (double*)NULL) tpPTorsion->dScee = *dPScee;
+    if (      dPScnb != (double*)NULL) tpPTorsion->dScnb = *dPScnb;
     if (sDescription != (char*)NULL  ) strcpy(tpPTorsion->sDesc, sDescription);
 
     strcpy( sOrder, "0123" );
@@ -2777,7 +2809,7 @@
 void
 ParmSetUpdateImproper( PARMSET psLib, int i, 
 	char *sType1, char *sType2, char *sType3, char *sType4,
-	int *iPN, double *dPKp, double *dPP0, char *sDescription)
+	int *iPN, double *dPKp, double *dPP0, double *dPScee, double *dPScnb, char *sDescription)
 {
 TORSIONPARMt   *tpPTorsion;
 orderStr	sOrder;
@@ -2790,6 +2822,8 @@
     if (         iPN != (int*)NULL   ) tpPTorsion->iN = *iPN;
     if (        dPKp != (double*)NULL) tpPTorsion->dKp = *dPKp;
     if (        dPP0 != (double*)NULL) tpPTorsion->dP0 = *dPP0;
+    if (      dPScee != (double*)NULL) tpPTorsion->dScee = *dPScee;
+    if (      dPScnb != (double*)NULL) tpPTorsion->dScnb = *dPScnb;
     if (sDescription != (char*)NULL  ) strcpy(tpPTorsion->sDesc, sDescription);
 
     strcpy( sOrder, "0123" );
@@ -2979,6 +3013,7 @@
 			tpPCur->sType1, tpPCur->sType2, 
 			tpPCur->sType3, tpPCur->sType4,
 			tpPCur->iN, tpPCur->dKp, tpPCur->dP0/DEGTORAD,
+			tpPCur->dScee, tpPCur->dScnb,
 			tpPCur->sOrder, tpPCur->sDesc );
 	}
 }
@@ -2996,6 +3031,7 @@
 			tpPCur->sType1, tpPCur->sType2, 
 			tpPCur->sType3, tpPCur->sType4,
 			tpPCur->iN, tpPCur->dKp, tpPCur->dP0/DEGTORAD, 
+			tpPCur->dScee, tpPCur->dScnb,
 			tpPCur->sDesc );
 	}
 }
diff -urN amber11.orig/AmberTools/src/leap/src/leap/parmSet.h amber11/AmberTools/src/leap/src/leap/parmSet.h
--- amber11.orig/AmberTools/src/leap/src/leap/parmSet.h	2011-03-10 20:12:06.000000000 +0200
+++ amber11/AmberTools/src/leap/src/leap/parmSet.h	2011-10-25 15:01:28.144288609 +0300
@@ -42,6 +42,9 @@
  *
  *
  */
+        /*      Arunima Singh (UGA)     */
+        /*         20110420             */
+        /*      Added SCEE/SCNB         */
  
 #ifndef PARMSET_H
 #define PARMSET_H
@@ -96,6 +99,8 @@
 	int		iN;
 	double          dKp;
 	double          dP0;
+	double		dScee;		/* for propers */ /* Arunima Singh */ 
+	double          dScnb;		/* for propers */ /* Arunima Singh */
 	orderStr	sOrder;		/* for impropers */
 	DESCRIPTION     sDesc;
 } TORSIONPARMt;
@@ -158,10 +163,10 @@
 			double dKt, double dT0, double dTkub, double dRkub, char *sDesc);
 extern int	iParmSetAddProperTerm(PARMSET psLib,
 			char *sType1, char *sType2, char *sType3, char *sType4,
-			int iN, double dKp, double dP0, char *sDesc);
+			int iN, double dKp, double dP0, double dScee, double dScnb, char *sDesc);	/* for proper dihedrals */ /* Arunima Singh */
 extern int	iParmSetAddImproperTerm(PARMSET psLib,
 			char *sType1, char *sType2, char *sType3, char *sType4,
-			int iN, double dKp, double dP0, char *sDesc);
+			int iN, double dKp, double dP0, double dScee, double dScnb, char *sDesc);	/* Arunima Singh : added for proper dihedrals, but has to be included in impropers to maintain common data structure usage */ 
 extern int	iParmSetAddHBond(PARMSET psLib, char *sType1, char *sType2,
 			double dA, double dB, char *sDesc);
 
@@ -216,15 +221,15 @@
 extern void	ParmSetTORSIONTerm(TORSION tTorsion, int iTorsionIndex, 
 			int *iPParmSetIndex,
 			char *cPTyp1, char *cPTyp2, char *cPTyp3, char *cPTyp4,
-			int *iPN, double *dPKp, double *dPP0, char *sDesc );
+			int *iPN, double *dPKp, double *dPP0, double *dPScee, double *dPScnb, char *sDesc );
 extern BOOL	bParmSetTORSIONAddProperTerm(TORSION tTorsion,
 			char *cPType1, char *cPType2, 
 			char *cPType3, char *cPType4,
-			int iN, double dKp, double dP0, char *sDesc);
+			int iN, double dKp, double dP0, double dScee, double dScnb, char *sDesc);
 extern BOOL	bParmSetTORSIONAddImproperTerm(TORSION tTorsion,
 			char *cPType1, char *cPType2, 
 			char *cPType3, char *cPType4,
-			int iN, double dKp, double dP0, char *sDesc);
+			int iN, double dKp, double dP0, double dScee, double dScnb, char *sDesc);
 extern void	ParmSetTORSIONOrderAtoms();	
 extern void	ParmSetImproperOrderAtoms( TORSION tTorsion, int iTorsionIndex,
 			char *cPaTypes[4], int iaIndexes[4] );
@@ -268,10 +273,10 @@
 			char *sDesc);
 extern void	ParmSetTorsion(PARMSET psLib, int i,
 			char *sType1, char *sType2, char *sType3, char *sType4,
-			int *iPN, double *dPKp, double *dPP0, char *sDesc);
+			int *iPN, double *dPKp, double *dPP0, double *dPScee, double *dPScnb, char *sDesc);
 extern void	ParmSetImproper(PARMSET psLib, int i,
 			char *sType1, char *sType2, char *sType3, char *sType4,
-			int *iPN, double *dPKp, double *dPP0, char *sDesc);
+			int *iPN, double *dPKp, double *dPP0, double *dPScee, double *dPScnb, char *sDesc);
 extern void	ParmSetHBond(PARMSET psLib, int i, char *sType1, char *sType2,
 			double *dPA, double *dPB, char *sDesc);
 
@@ -296,11 +301,11 @@
 			double *dPKt, double *dPT0, char *sDescription);
 extern void	ParmSetUpdateTorsion(PARMSET psLib, int i,
 			char *sType1, char *sType2, char *sType3, char *sType4,
-			int *iPN, double *dPKp, double *dPP0, 
+			int *iPN, double *dPKp, double *dPP0, double *dPScee, double *dPScnb,
 			char *sDescription);
 extern void	ParmSetUpdateImproper(PARMSET psLib, int i,
 			char *sType1, char *sType2, char *sType3, char *sType4,
-			int *iPN, double *dPKp, double *dPP0, 
+			int *iPN, double *dPKp, double *dPP0, double *dPScee, double *dPScnb,
 			char *sDescription);
 extern void	ParmSetUpdateHBond(PARMSET psLib, int i, 
 			char *sType1, char *sType2,
diff -urN amber11.orig/AmberTools/src/leap/src/leap/unitio.c amber11/AmberTools/src/leap/src/leap/unitio.c
--- amber11.orig/AmberTools/src/leap/src/leap/unitio.c	2011-03-10 20:12:06.000000000 +0200
+++ amber11/AmberTools/src/leap/src/leap/unitio.c	2011-10-25 15:01:28.192288568 +0300
@@ -70,6 +70,10 @@
 *       zbUnitIOIndexBondParameters and zUnitDoAtoms are now "extern functions" 
 */ 
 
+        /*      Arunima Singh (UGA)     */
+        /*         20110420             */
+        /*      Added SCEE/SCNB         */
+
 #include <time.h>
 
 #include        "basics.h"
@@ -1585,7 +1589,9 @@
         tC.dKp = tpA->dKp;
         tC.iN = tpA->iN;
         tC.dP0 = tpA->dP0;
-        VarArrayAdd(vaB, (GENP) & tC);
+        tC.dScee = tpA->dScee;
+	tC.dScnb = tpA->dScnb;
+	VarArrayAdd(vaB, (GENP) & tC);
         iIndex++;
         iOldIndex = iParmOffset + iA + 1;
 
@@ -1615,6 +1621,10 @@
                 continue;
             if (tpB->dP0 != tpA->dP0)
                 continue;
+	    if (tpB->dScee != tpA->dScee)
+                continue;
+            if (tpB->dScnb != tpA->dScnb)
+                continue;
 
             /*
              *  B is a duplicate of A
@@ -1698,12 +1708,12 @@
     int iTerm, iPertTerm;
     BOOL bDone, bUse, bUsePert, bCopy, bCopyPert, bEnd, bPertEnd;
     int iN, iPertIndex, iPertN, iLastN, iLastPertN;
-    double dKp, dP0, dPertKp, dPertP0;
+    double dKp, dP0, dScee, dScnb, dPertKp, dPertP0, dPertScee, dPertScnb;
     BOOL bCalc14, bCalcPert14;
 #ifdef  DEBUG2
     STRING s1, s2, s3, s4;
     int iTParm, iTmp;
-    double dTK, dTP;
+    double dTK, dTP, dTScee, dTScnb;
     STRING sT1, sT2, sT3, sT4, sTemp;
 #endif
         STRING sDesc;
@@ -2016,9 +2026,9 @@
             ParmSetTORSIONTerm(tTorsion, i,
                                &iTParm,
                                sT1, sT2, sT3, sT4,
-                               &iTmp, &dTK, &dTP, sTemp);
-            MESSAGE(("Term %3d  %d %s-%s-%s-%s  %d  %lf  %lf\n",
-                     i, iTParm, sT1, sT2, sT3, sT4, iTmp, dTK, dTP));
+                               &iTmp, &dTK, &dTP, &dTScee, &dTScnb, sTemp);
+            MESSAGE(("Term %3d  %d %s-%s-%s-%s  %d  %lf  %lf %lf %lf\n",
+                     i, iTParm, sT1, sT2, sT3, sT4, iTmp, dTK, dTP, dTScee, dTScnb));
         }
         if (bPerturbTorsion) {
             MESSAGE(("Pert%s %s-%s-%s-%s found %d terms\n",
@@ -2029,9 +2039,9 @@
                 ParmSetTORSIONTerm(tPertTorsion, i,
                                    &iTParm,
                                    sT1, sT2, sT3, sT4,
-                                   &iTmp, &dTK, &dTP, sTemp);
-                MESSAGE(("Term %3d  %d %s-%s-%s-%s  %d  %lf  %lf\n",
-                         i, iTParm, sT1, sT2, sT3, sT4, iTmp, dTK, dTP));
+                                   &iTmp, &dTK, &dTP, &dTScee, &dTScnb, sTemp);
+                MESSAGE(("Term %3d  %d %s-%s-%s-%s  %d  %lf  %lf %lf %lf\n",
+                         i, iTParm, sT1, sT2, sT3, sT4, iTmp, dTK, dTP dTScee, dTScnb));
             }
         }
 #endif
@@ -2061,7 +2071,7 @@
             ParmSetTORSIONTerm(tTorsion, iTerm,
                                &iIndex,
                                sAtom1, sAtom2, sAtom3, sAtom4,
-                               &iN, &dKp, &dP0, sDesc);
+                               &iN, &dKp, &dP0, &dScee, &dScnb, sDesc);
             MESSAGE(("First non-perturbed multiplicity: %d\n", iN));
         } else {
             if (bProper) {
@@ -2082,7 +2092,7 @@
                 ParmSetTORSIONTerm(tPertTorsion, iPertTerm,
                                    &iPertIndex,
                                    sPert1, sPert2, sPert3, sPert4,
-                                   &iPertN, &dPertKp, &dPertP0, sDesc);
+                                   &iPertN, &dPertKp, &dPertP0, &dPertScee, &dPertScnb, sDesc);
                 MESSAGE(("First perturbed multiplicity: %d\n", iPertN));
             } else
                 bPertEnd = TRUE;
@@ -2175,14 +2185,14 @@
                 if (bProper)
                     iIndex = iParmSetAddProperTerm(uUnit->psParameters,
                                                    sAtom1, sAtom2, sAtom3,
-                                                   sAtom4, iN, dKp, dP0,
-                                                   sDesc);
+                                                   sAtom4, iN, dKp, dP0, 
+                                                   dScee, dScnb, sDesc);
 /*                else if ( !GDefaults.iCharmm )    ???---should I do this????     */
                 else
                     iIndex = iParmSetAddImproperTerm(uUnit->psParameters,
                                                      sAtom1, sAtom2,
                                                      sAtom3, sAtom4, iN,
-                                                     dKp, dP0, sDesc);
+                                                     dKp, dP0, dScee, dScnb, sDesc);
             }
             if (bCopyPert) {
                 if (bProper) {
@@ -2190,13 +2200,13 @@
                                                        sPert1, sPert2,
                                                        sPert3, sPert4,
                                                        iPertN, dPertKp,
-                                                       dPertP0, sDesc);
+                                                       dPertP0, dScee, dScnb, sDesc);
                 } else {
                     iPertIndex =
                         iParmSetAddImproperTerm(uUnit->psParameters,
                                                 sPert1, sPert2, sPert3,
                                                 sPert4, iPertN, dPertKp,
-                                                dPertP0, sDesc);
+                                                dPertP0, dScee, dScnb, sDesc);
                 }
                 MESSAGE(("iPertIndex = %d\n", iPertIndex));
             }
@@ -2215,7 +2225,7 @@
                     ParmSetTORSIONTerm(tTorsion, iTerm,
                                        &iIndex,
                                        sAtom1, sAtom2, sAtom3, sAtom4,
-                                       &iN, &dKp, &dP0, sDesc);
+                                       &iN, &dKp, &dP0, &dScee, &dScnb, sDesc);
                 }
                 MESSAGE(
                         ("Advancing non-perturbed multiplicity to %d\n",
@@ -2231,7 +2241,7 @@
                     ParmSetTORSIONTerm(tPertTorsion, iPertTerm,
                                        &iPertIndex,
                                        sPert1, sPert2, sPert3, sPert4,
-                                       &iPertN, &dPertKp, &dPertP0, sDesc);
+                                       &iPertN, &dPertKp, &dPertP0, &dPertScee, &dPertScnb, sDesc);
                 }
                 MESSAGE(
                         ("Advancing perturbed multiplicity to %d\n",
@@ -3988,7 +3998,7 @@
     SAVEATOMt *saPAtom;
     SAVETORSIONt *stPTorsion;
     SAVERESTRAINTt *srPRestraint;
-    double dMass, dPolar, dR, dKb, dR0, dKt, dT0, dTkub, dRkub, dKp, dP0,
+    double dMass, dPolar, dR, dKb, dR0, dKt, dT0, dTkub, dRkub, dKp, dP0, dScee, dScnb,
         dC, dD, dTemp;
     STRING sAtom1, sAtom2, sAtom3, sAtom4, sType1, sType2;
     int iN, iAtoms, iMaxAtoms, iTemp, iAtom, iCalc14, iProper;
@@ -4585,16 +4595,16 @@
              iParmSetTotalImproperParms(uUnit->psParameters)));
     for (i = 0; i < iParmSetTotalTorsionParms(uUnit->psParameters); i++) {
         ParmSetTorsion(uUnit->psParameters, i, sAtom1, sAtom2,
-                       sAtom3, sAtom4, &iN, &dKp, &dP0, sDesc);
-        MESSAGE(("Torsion %d  %s-%s-%s-%s %d %lf %lf\n",
-                 i, sAtom1, sAtom2, sAtom3, sAtom4, iN, dKp, dP0));
+                       sAtom3, sAtom4, &iN, &dKp, &dP0, &dScee, &dScnb, sDesc);
+        MESSAGE(("Torsion %d  %s-%s-%s-%s %d %lf %lf %lf %lf\n",
+                 i, sAtom1, sAtom2, sAtom3, sAtom4, iN, dKp, dP0, dScee, dScnb));
         FortranWriteDouble(dKp);
     }
     for (i = 0; i < iParmSetTotalImproperParms(uUnit->psParameters); i++) {
         ParmSetImproper(uUnit->psParameters, i, sAtom1, sAtom2,
-                        sAtom3, sAtom4, &iN, &dKp, &dP0, sDesc);
-        MESSAGE(("Improper %d  %s-%s-%s-%s %d %lf %lf\n",
-                 i, sAtom1, sAtom2, sAtom3, sAtom4, iN, dKp, dP0));
+                        sAtom3, sAtom4, &iN, &dKp, &dP0, &dScee, &dScnb, sDesc);
+        MESSAGE(("Improper %d  %s-%s-%s-%s %d %lf %lf %lf %lf\n",
+                 i, sAtom1, sAtom2, sAtom3, sAtom4, iN, dKp, dP0, dScee, dScnb));
         FortranWriteDouble(dKp);
     }
     /* Write the torsion RESTRAINT constants AND set the index */
@@ -4612,13 +4622,13 @@
     FortranFormat(5, DBLFORMAT);
     for (i = 0; i < iParmSetTotalTorsionParms(uUnit->psParameters); i++) {
         ParmSetTorsion(uUnit->psParameters, i, sAtom1, sAtom2,
-                       sAtom3, sAtom4, &iN, &dKp, &dP0, sDesc);
+                       sAtom3, sAtom4, &iN, &dKp, &dP0, &dScee, &dScnb, sDesc);
         dTemp = iN;
         FortranWriteDouble(dTemp);
     }
     for (i = 0; i < iParmSetTotalImproperParms(uUnit->psParameters); i++) {
         ParmSetImproper(uUnit->psParameters, i, sAtom1, sAtom2,
-                        sAtom3, sAtom4, &iN, &dKp, &dP0, sDesc);
+                        sAtom3, sAtom4, &iN, &dKp, &dP0, &dScee, &dScnb, sDesc);
         dTemp = iN;
         FortranWriteDouble(dTemp);
     }
@@ -4637,12 +4647,12 @@
     FortranFormat(5, DBLFORMAT);
     for (i = 0; i < iParmSetTotalTorsionParms(uUnit->psParameters); i++) {
         ParmSetTorsion(uUnit->psParameters, i, sAtom1, sAtom2,
-                       sAtom3, sAtom4, &iN, &dKp, &dP0, sDesc);
+                       sAtom3, sAtom4, &iN, &dKp, &dP0, &dScee, &dScnb, sDesc);
         FortranWriteDouble(dP0);
     }
     for (i = 0; i < iParmSetTotalImproperParms(uUnit->psParameters); i++) {
         ParmSetImproper(uUnit->psParameters, i, sAtom1, sAtom2,
-                        sAtom3, sAtom4, &iN, &dKp, &dP0, sDesc);
+                        sAtom3, sAtom4, &iN, &dKp, &dP0, &dScee, &dScnb, sDesc);
         FortranWriteDouble(dP0);
     }
     /* Write the torsion RESTRAINT constants AND set the index */
@@ -4650,6 +4660,60 @@
     RESTRAINTLOOP(RESTRAINTTORSION, dN, i + 1);
     FortranEndLine();
 
+	/*	Arunima Singh (UGA)	*/
+	/*	   20110420		*/
+	/*   Added Printing SCEE/SCNB	*/
+
+
+    /* -17.5A- SCEE for torsions */
+    FortranDebug("-17.5A-");
+    
+    MESSAGE(("Writing SCEE scale factor for torsion interaction\n"));
+    FortranFormat(1, "%-80s");
+    FortranWriteString("%FLAG SCEE_SCALE_FACTOR");
+    FortranWriteString("%FORMAT(5E16.8)");
+    FortranFormat(5, DBLFORMAT);
+    for (i = 0; i < iParmSetTotalTorsionParms(uUnit->psParameters); i++) {
+        ParmSetTorsion(uUnit->psParameters, i, sAtom1, sAtom2,
+                       sAtom3, sAtom4, &iN, &dKp, &dP0, &dScee, &dScnb, sDesc);
+        FortranWriteDouble(dScee);
+    }
+    for (i = 0; i < iParmSetTotalImproperParms(uUnit->psParameters); i++) {
+        ParmSetImproper(uUnit->psParameters, i, sAtom1, sAtom2,
+                        sAtom3, sAtom4, &iN, &dKp, &dP0, &dScee, &dScnb, sDesc);
+        FortranWriteDouble(dScee);
+    }
+    /* Write the torsion RESTRAINT constants AND set the index */
+    /* for where the interaction can find its constants */
+    RESTRAINTLOOP(RESTRAINTTORSION, dX0, i + 1);
+    FortranEndLine();
+
+    /* -17.5B- SCNB for torsions */
+    FortranDebug("-17.5B-");
+
+    MESSAGE(("Writing SCNB scale factor for torsion interaction\n"));
+    FortranFormat(1, "%-80s");
+    FortranWriteString("%FLAG SCNB_SCALE_FACTOR");
+    FortranWriteString("%FORMAT(5E16.8)");
+    FortranFormat(5, DBLFORMAT);
+    for (i = 0; i < iParmSetTotalTorsionParms(uUnit->psParameters); i++) {
+        ParmSetTorsion(uUnit->psParameters, i, sAtom1, sAtom2,
+                       sAtom3, sAtom4, &iN, &dKp, &dP0, &dScee, &dScnb, sDesc);
+        FortranWriteDouble(dScnb);
+    }
+    for (i = 0; i < iParmSetTotalImproperParms(uUnit->psParameters); i++) {
+        ParmSetImproper(uUnit->psParameters, i, sAtom1, sAtom2,
+                        sAtom3, sAtom4, &iN, &dKp, &dP0, &dScee, &dScnb, sDesc);
+        FortranWriteDouble(dScnb);
+    }
+    /* Write the torsion RESTRAINT constants AND set the index */
+    /* for where the interaction can find its constants */
+    RESTRAINTLOOP(RESTRAINTTORSION, dX0, i + 1);
+    FortranEndLine();
+	
+	/*     	    Arunima Singh	*/
+	/*    End of adding SCEE/SCNB	*/
+
     /* -18- Not used, reserved for future use, uses NATYP */
     /* Corresponds to the AMBER SOLTY array */
     FortranDebug("-18-");
@@ -6009,7 +6073,7 @@
 SAVEATOMt       *saPAtom;
 SAVETORSIONt    *stPTorsion;
 SAVERESTRAINTt  *srPRestraint;
-double          dMass, dPolar, dR, dKb, dR0, dKt, dT0, dTkub, dRkub, dKp, dP0, dC, dD, dTemp;
+double          dMass, dPolar, dR, dKb, dR0, dKt, dT0, dTkub, dRkub, dKp, dP0, dScee, dScnb, dC, dD, dTemp;
 STRING          sAtom1, sAtom2, sAtom3, sAtom4, sType1, sType2;
 int             iN, iAtoms, iMaxAtoms, iTemp, iAtom, iCalc14, iProper;
 int             iElement, iHybridization, iStart, iFirstSolvent;
@@ -6547,19 +6611,19 @@
     for ( i=0; i<iParmSetTotalTorsionParms(uUnit->psParameters); i++ ) {
         ParmSetTorsion( uUnit->psParameters, i, sAtom1, sAtom2, 
                         sAtom3, sAtom4,
-                        &iN, &dKp, &dP0, sDesc );
-        MESSAGE(( "Torsion %d  %s-%s-%s-%s %d %lf %lf\n",
+                        &iN, &dKp, &dP0, &dScee, &dScnb, sDesc );
+        MESSAGE(( "Torsion %d  %s-%s-%s-%s %d %lf %lf %lf %lf\n",
                         i, sAtom1, sAtom2, sAtom3, sAtom4,
-                        iN, dKp, dP0 ));
+                        iN, dKp, dP0, dScee, dScnb ));
         FortranWriteDouble( dKp );
     }
     for ( i=0; i<iParmSetTotalImproperParms(uUnit->psParameters); i++ ) {
         ParmSetImproper( uUnit->psParameters, i, sAtom1, sAtom2, 
                         sAtom3, sAtom4,
-                        &iN, &dKp, &dP0, sDesc );
-        MESSAGE(( "Improper %d  %s-%s-%s-%s %d %lf %lf\n",
+                        &iN, &dKp, &dP0, &dScee, &dScnb, sDesc );
+        MESSAGE(( "Improper %d  %s-%s-%s-%s %d %lf %lf %lf %lf\n",
                         i, sAtom1, sAtom2, sAtom3, sAtom4,
-                        iN, dKp, dP0 ));
+                        iN, dKp, dP0, dScee, dScnb ));
         FortranWriteDouble( dKp );
     }
                 /* Write the torsion RESTRAINT constants AND set the index */
@@ -6575,14 +6639,14 @@
     for ( i=0; i<iParmSetTotalTorsionParms(uUnit->psParameters); i++ ) {
         ParmSetTorsion( uUnit->psParameters, i, sAtom1, sAtom2, 
                         sAtom3, sAtom4,
-                        &iN, &dKp, &dP0, sDesc );
+                        &iN, &dKp, &dP0, &dScee, &dScnb, sDesc );
         dTemp = iN;
         FortranWriteDouble( dTemp );
     }
     for ( i=0; i<iParmSetTotalImproperParms(uUnit->psParameters); i++ ) {
         ParmSetImproper( uUnit->psParameters, i, sAtom1, sAtom2, 
                         sAtom3, sAtom4,
-                        &iN, &dKp, &dP0, sDesc );
+                        &iN, &dKp, &dP0, &dScee, &dScnb, sDesc );
         dTemp = iN;
         FortranWriteDouble( dTemp );
     }
@@ -6599,13 +6663,13 @@
     for ( i=0; i<iParmSetTotalTorsionParms(uUnit->psParameters); i++ ) {
         ParmSetTorsion( uUnit->psParameters, i, sAtom1, sAtom2, 
                         sAtom3, sAtom4,
-                        &iN, &dKp, &dP0, sDesc );
+                        &iN, &dKp, &dP0, &dScee, &dScnb, sDesc );
         FortranWriteDouble( dP0 );
     }
     for ( i=0; i<iParmSetTotalImproperParms(uUnit->psParameters); i++ ) {
         ParmSetImproper( uUnit->psParameters, i, sAtom1, sAtom2, 
                         sAtom3, sAtom4,
-                        &iN, &dKp, &dP0, sDesc );
+                        &iN, &dKp, &dP0, &dScee, &dScnb, sDesc );
         FortranWriteDouble( dP0 );
     }
                 /* Write the torsion RESTRAINT constants AND set the index */
diff -urN amber11.orig/AmberTools/src/leap/src/leap/xaImproperParmTable.c amber11/AmberTools/src/leap/src/leap/xaImproperParmTable.c
--- amber11.orig/AmberTools/src/leap/src/leap/xaImproperParmTable.c	2011-03-10 20:12:06.000000000 +0200
+++ amber11/AmberTools/src/leap/src/leap/xaImproperParmTable.c	2011-10-25 15:01:28.193288566 +0300
@@ -37,7 +37,9 @@
  *		Handle editing of parameters in a table format.
  */
 
-
+        /*      Arunima Singh (UGA)     */
+        /*         20110420             */
+        /*      Added SCEE/SCNB         */
 
 #include	<X11/IntrinsicP.h>
 #include	<X11/StringDefs.h>
@@ -59,7 +61,10 @@
 #define	NC		4
 #define	KPC		5
 #define P0C		6
-#define DESCC		7
+#define SCEE            7       /*Arunima Singh*/
+#define SCNB            8       /*Arunima Singh*/
+#define DESCC           9       /*Arunima Singh: Changed value from 7 to 9 to include scee and scnb*/
+ 
 
 #define	MAXTYPELEN	5
 #define DESCLEN		32
@@ -129,13 +134,15 @@
 int			iN;
 double			dKp;
 double			dP0;
+double			dScee;
+double			dScnb;
 char			sDesc[DESCLEN];
 
 
     iptPCur = (IMPROPERPARMTABLEt*)PXATClientPointer(tTable);
     
     ParmSetImproper( iptPCur->psParmSet, iRow, 
-	    sType1, sType2, sType3, sType4, &iN, &dKp, &dP0, sDesc );
+	    sType1, sType2, sType3, sType4, &iN, &dKp, &dP0, &dScee, &dScnb, sDesc );
 
     switch ( iCol ) {
 	    case TYPE1C:
@@ -171,6 +178,14 @@
 	    	}
 		return(SsBuffer);
 		break;
+            case SCEE:
+                sprintf( SsBuffer, DBLFMT, dScee );
+                return(SsBuffer);
+                break;
+            case SCNB:
+                sprintf( SsBuffer, DBLFMT, dScnb );
+                return(SsBuffer);
+                break;
 	    case DESCC:
 	    	strcpy( SsBuffer, sDesc );
 		return( SsBuffer );
@@ -259,6 +274,16 @@
 		if ( !strcmp( cPData, "0" ))  break;
 		return("P0 field must be '0' or 'Pi'.");
 		break;
+            case SCEE:
+                if ( !bStringToDouble( cPData, &dValue ) ) {
+                    return("Invalid character in Scee field.");
+                }
+                break;
+            case SCNB:
+                if ( !bStringToDouble( cPData, &dValue ) ) {
+                    return("Invalid character in Scnb field.");
+                }
+                break;
 	    case DESCC:
 		if ( strlen(cPData)>DESCLEN-1 ) {
 		    sprintf( SsError, "%s %d characters.",
@@ -384,11 +409,33 @@
 		return("P0 field must be '0' or 'Pi'.");
 	}
 
+        /*
+         *  Scee
+         */
+        *iPErrCol = 7;
+        cPData = col[7];
+        if ( !bStringToDouble( cPData, &dValue ) ) {
+                return("Invalid character in Scee field.");
+        } else if ( dValue < 0.0 ) {
+                return("Scee cannot be negative.");
+        }
+
+        /*
+         *  Scnb
+         */
+        *iPErrCol = 8;
+        cPData = col[8];
+        if ( !bStringToInt( cPData, &iValue ) ) {
+                return("Invalid character in Scnb field.");
+        } else if ( dValue < 0.0 ) {
+                return("Scnb cannot be negative.");
+        }
+
 	/*
 	 *  Desc
 	 */
-	*iPErrCol = 7;
-	cPData = col[7];
+	*iPErrCol = 9;
+	cPData = col[9];
 	if ( strlen(cPData) > DESCLEN-1 ) {
 		sprintf( SsError, "%s %d characters.",
 			"Parameter Description cannot be longer than",
@@ -413,7 +460,8 @@
 int			iN;
 double			dKp;
 double			dP0;
-
+double			dScee;
+double			dScnb;
 
 	iptPCur = (IMPROPERPARMTABLEt*) PXATClientPointer(tTable);
 	psParmSet = iptPCur->psParmSet;
@@ -444,7 +492,7 @@
 					col[0], col[1], col[2], col[3] ));
         	if ( iRow != iParmSetAddImproperTerm( psParmSet,
 					col[0], col[1], col[2], col[3],
-					iN, dKp, dP0, col[7] ) )
+					iN, dKp, dP0, dScee, dScnb, col[9] ) )
 			DFATAL(( "programming err 2 in zXAIPTAcceptRow\n" ));
 	} else {
 		/*
@@ -452,7 +500,7 @@
 		 */
 		ParmSetUpdateImproper( psParmSet, iRow, 
 				col[0], col[1], col[2], col[3],
-				&iN, &dKp, &dP0, col[7] );
+				&iN, &dKp, &dP0, &dScee, &dScnb, col[9] );
 	}
 }
 
diff -urN amber11.orig/AmberTools/src/leap/src/leap/xaTorsionParmTable.c amber11/AmberTools/src/leap/src/leap/xaTorsionParmTable.c
--- amber11.orig/AmberTools/src/leap/src/leap/xaTorsionParmTable.c	2011-03-10 20:12:06.000000000 +0200
+++ amber11/AmberTools/src/leap/src/leap/xaTorsionParmTable.c	2011-10-25 15:01:28.194288563 +0300
@@ -37,6 +37,9 @@
  *		Handle editing of parameters in a table format.
  */
 
+        /*      Arunima Singh (UGA)     */
+        /*         20110420             */
+        /*      Added SCEE/SCNB         */
 
 #include	<X11/IntrinsicP.h>
 #include	<X11/StringDefs.h>
@@ -58,7 +61,9 @@
 #define	NC		4
 #define	KPC		5
 #define P0C		6
-#define DESCC		7
+#define SCEE		7	/*Arunima Singh*/
+#define SCNB		8	/*Arunima Singh*/
+#define DESCC		9	/*Arunima Singh: Changed value from 7 to 9 to include scee and scnb*/
 
 #define	MAXTYPELEN	5
 #define DESCLEN		32
@@ -112,7 +117,7 @@
  *	zcPXATPTGetElement
  *
  *	Get the values for the elements of the TABLE from the
- *	particular Torsion Parmeter Entry.
+ *	particular Torsion Parameter Entry.
  */
 static char *
 zcPXATPTGetElement( TABLE tTable, int iCol, int iRow )
@@ -126,13 +131,15 @@
 int			iN;
 double			dKp;
 double			dP0;
+double			dScee;
+double			dScnb;
 char			sDesc[DESCLEN];
 
 
     tptPCur = (TORSIONPARMTABLEt*)PXATClientPointer(tTable);
     
     ParmSetTorsion( tptPCur->psParmSet, iRow, 
-	    sType1, sType2, sType3, sType4, &iN, &dKp, &dP0, sDesc );
+	    sType1, sType2, sType3, sType4, &iN, &dKp, &dP0, &dScee, &dScnb, sDesc );
     
     switch ( iCol ) {
 	    case TYPE1C:
@@ -168,6 +175,14 @@
 	    	}
 		return(SsBuffer);
 		break;
+	    case SCEE:
+		sprintf( SsBuffer, DBLFMT, dScee );
+		return(SsBuffer);
+                break;
+            case SCNB:
+                sprintf( SsBuffer, DBLFMT, dScnb );
+                return(SsBuffer);
+                break; 
 	    case DESCC:
 	    	strcpy( SsBuffer, sDesc );
 		return( SsBuffer );
@@ -256,6 +271,16 @@
 		if ( !strcmp( cPData, "0" ))  break;
 		return("P0 field must be '0' or 'Pi'.");
 		break;
+            case SCEE:
+                if ( !bStringToDouble( cPData, &dValue ) ) {
+                    return("Invalid character in Scee field.");
+                }
+                break;
+            case SCNB:
+                if ( !bStringToDouble( cPData, &dValue ) ) {
+                    return("Invalid character in Scnb field.");
+                }
+                break;
 	    case DESCC:
 		if ( strlen(cPData)>DESCLEN-1 ) {
 		    sprintf( SsError, "%s %d characters.",
@@ -379,11 +404,33 @@
 		return("P0 field must be '0' or 'Pi'.");
 	}
 
+        /*
+         *  Scee
+         */
+        *iPErrCol = 7;
+        cPData = col[7];
+        if ( !bStringToDouble( cPData, &dValue ) ) {
+                return("Invalid character in Scee field.");
+        } else if ( dValue < 0.0 ) {
+                return("Scee cannot be negative.");
+        }
+
+        /*
+         *  Scnb
+         */
+        *iPErrCol = 8;
+        cPData = col[8];
+        if ( !bStringToInt( cPData, &iValue ) ) {
+                return("Invalid character in Scnb field.");
+        } else if ( dValue < 0.0 ) {
+                return("Scnb cannot be negative.");
+        }
+
 	/*
 	 *  Desc
 	 */
-	*iPErrCol = 7;
-	cPData = col[7];
+	*iPErrCol = 9;
+	cPData = col[9];
 	if ( strlen(cPData) > DESCLEN-1 ) {
 		sprintf( SsError, "%s %d characters.",
 			"Parameter Description cannot be longer than",
@@ -408,6 +455,8 @@
 int			iN;
 double			dKp;
 double			dP0;
+double			dScee;
+double			dScnb;
 
 	tptPCur = (TORSIONPARMTABLEt*) PXATClientPointer(tTable);
 	psParmSet = tptPCur->psParmSet;
@@ -438,15 +487,15 @@
 					col[0], col[1], col[2], col[3] ));
         	if ( iRow != iParmSetAddProperTerm( psParmSet,
 					col[0], col[1], col[2], col[3],
-					iN, dKp, dP0, col[7] ) )
-			DFATAL(( "programming err 2 in zXATPTAcceptRow\n" ));
+					iN, dKp, dP0, dScee, dScnb, col[9] ) )	/*Arunima singh*/
+			DFATAL(( "programming err 2 in zXATPTAcceptRow\n" ));	
 	} else {
 		/*
 		 *  update row in place
 		 */
 		ParmSetUpdateTorsion( psParmSet, iRow, 
 				col[0], col[1], col[2], col[3],
-				&iN, &dKp, &dP0, col[7] );
+				&iN, &dKp, &dP0, &dScee, &dScnb, col[9] );	/*Arunima singh*/
 	}
 }
 
diff -urN amber11.orig/AmberTools/src/mmpbsa_py/MMPBSA_mods/alamdcrd.py amber11/AmberTools/src/mmpbsa_py/MMPBSA_mods/alamdcrd.py
--- amber11.orig/AmberTools/src/mmpbsa_py/MMPBSA_mods/alamdcrd.py	2011-04-14 15:30:17.000000000 +0300
+++ amber11/AmberTools/src/mmpbsa_py/MMPBSA_mods/alamdcrd.py	2011-10-25 15:01:28.132288620 +0300
@@ -371,7 +371,7 @@
          coords_received = _scaledistance(coords_tosend, chdist)
 
          for i in range(3):
-            new_coords.append(coords_received[x+3])
+            new_coords.append(coords_received[i+3])
 
          coords_tosend = []
          coords_received = []
diff -urN amber11.orig/AmberTools/src/mmpbsa_py/MMPBSA.pypp amber11/AmberTools/src/mmpbsa_py/MMPBSA.pypp
--- amber11.orig/AmberTools/src/mmpbsa_py/MMPBSA.pypp	2011-04-14 15:30:17.000000000 +0300
+++ amber11/AmberTools/src/mmpbsa_py/MMPBSA.pypp	2011-10-25 15:01:28.238288523 +0300
@@ -412,6 +412,7 @@
    INPUT['solvcut']        = float(INPUT['solvcut'])
    INPUT['tolerance']      = float(INPUT['tolerance'])
    INPUT['rism_verbose']   =   int(INPUT['rism_verbose'])
+   INPUT['inp']            =   int(INPUT['inp'])
 except ValueError, err:
    print >> sys.stderr, 'Error: Invalid input data types! Check input file for proper float/int/string arguments.'
    print >> sys.stderr, '       %s' % err
@@ -998,7 +999,7 @@
          if master: MMPBSA_timer.StartTimer('ptraj')
 
          # create the dummy inpcrd files
-         if isinerr != 1 and (INPUT['gbrun'] or INPUT['pbrun']):
+         if isinerr != 1 and (INPUT['gbrun'] or INPUT['pbrun'] or INPUT['rismrun']):
             os.system('%s %s _MMPBSA_mutant_complexinpcrd.in > _MMPBSA_ptraj9.out 2>&1' % (ptraj, FILES['mutant_complex_prmtop']))
             if not stability and FILES['mutant_receptor_prmtop'] == FILES['receptor_prmtop']:
                os.system('%s %s _MMPBSA_mutant_ligandinpcrd.in > _MMPBSA_ptraj11.out 2>&1' % (ptraj, FILES['mutant_ligand_prmtop']))
diff -urN amber11.orig/AmberTools/src/mmpbsa_py/setup.sh amber11/AmberTools/src/mmpbsa_py/setup.sh
--- amber11.orig/AmberTools/src/mmpbsa_py/setup.sh	2011-04-14 15:30:18.000000000 +0300
+++ amber11/AmberTools/src/mmpbsa_py/setup.sh	2011-10-25 15:01:28.133288619 +0300
@@ -10,7 +10,8 @@
 
 # require AMBERHOME
 if [ -z $AMBERHOME ]; then
-   export AMBERHOME=`dirname \`dirname $PWD\``
+   AMBERHOME=`dirname \`dirname $PWD\``
+   export AMBERHOME=`dirname $AMBERHOME`
    echo "AMBERHOME is not set.  Assuming it is $AMBERHOME"
 fi
 
diff -urN amber11.orig/AmberTools/src/pbsa/amg1r5.f amber11/AmberTools/src/pbsa/amg1r5.f
--- amber11.orig/AmberTools/src/pbsa/amg1r5.f	2011-04-14 15:30:19.000000000 +0300
+++ amber11/AmberTools/src/pbsa/amg1r5.f	2011-10-25 15:01:28.203288555 +0300
@@ -3707,7 +3707,7 @@
 
 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 !+ [Enter a one-line description of real function cgeps here]
-real function cgeps(k,s2,a,u,f,ia,ja,iw, &
+_REAL_ function cgeps(k,s2,a,u,f,ia,ja,iw, &
       imin,imax,iminw,m,ierr,ium)
    
    implicit _REAL_ (a-h,o-z)
@@ -3753,7 +3753,7 @@
 
 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 !+ [Enter a one-line description of real function cgalf here]
-real function cgalf(k,s2,a,u,f,ia,ja,iw, &
+_REAL_ function cgalf(k,s2,a,u,f,ia,ja,iw, &
       imin,imax,iminw,m)
    
    implicit _REAL_ (a-h,o-z)
@@ -3971,7 +3971,7 @@
 
 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 !+ [Enter a one-line description of real function random here]
-real function random(s)
+_REAL_ function random(s)
    
    !     FUNCTION TO CREATE "RANDOM" SEQUENCE OF NUMBERS BETWEEN 0 AND 0.1
    
diff -urN amber11.orig/AmberTools/src/pbsa/Makefile amber11/AmberTools/src/pbsa/Makefile
--- amber11.orig/AmberTools/src/pbsa/Makefile	2011-04-14 15:30:19.000000000 +0300
+++ amber11/AmberTools/src/pbsa/Makefile	2011-10-25 15:01:28.228288532 +0300
@@ -204,6 +204,14 @@
 	   cd ../c9x-complex && $(MAKE) libmc.a; \
 	 fi
 
+pb_init.o: pb_init.f
+	$(FPP) $(FPPFLAGS) $(AMBERFPPFLAGS) $< > _$<
+	$(FC) -c -O0 $(FFLAGS) $(AMBERFFLAGS) -o $@ _$<
+
+pb_init.SANDER.o: pb_init.f
+	$(FPP) $(FPPFLAGS) $(AMBERFPPFLAGS) -DSANDER $< > _$<
+	$(FC) -c -O0 $(FFLAGS) $(AMBERFFLAGS) -o $@ _$<
+
 %.LIBPBSA.o: %.f
 	$(FPP) $(FPPFLAGS) $(AMBERFPPFLAGS) -DLIBPBSA $< > _$<
 	$(FC) -c $(FOPTFLAGS) $(FFLAGS) $(AMBERFFLAGS) -o $@ _$<
diff -urN amber11.orig/AmberTools/src/ptraj/actions.c amber11/AmberTools/src/ptraj/actions.c
--- amber11.orig/AmberTools/src/ptraj/actions.c	2011-04-14 15:30:19.000000000 +0300
+++ amber11/AmberTools/src/ptraj/actions.c	2011-10-25 15:01:28.235288526 +0300
@@ -46,6 +46,7 @@
 #include "ptraj.h"
 #include <string.h>
 #include <stdarg.h>
+#include <float.h>
 
 /*
  *  The code in this file implements the various "actions" of ptraj to
@@ -4809,8 +4810,14 @@
    *  set the minimum distance to the solute to be larger than 
    *  any possible (imaged!!!) distance
    */
-  min = box[0] + box[1] + box[2];
-  min = min * min;
+  if (box[0]!=0.0 && box[1]!=0.0 && box[2]!=0.0) {
+    min = box[0] + box[1] + box[2];
+    min = min * min;
+  } else {
+    // No box information, set to arbitrarily large max
+    min=DBL_MAX;
+  }
+
   for (i=0; i < oldstate->solventMolecules; i++) {
     minDistance[i] = min;
     sortindex[i] = i;
@@ -16752,9 +16759,12 @@
      */
 
     structInfo = (transformSecondaryStructInfo *) action->carg1;
-    if (structInfo->filename != NULL)
+    if (structInfo->filename != NULL) {
+      i = (int) strlen(structInfo->filename);
+      i+=5;
+      structInfo->filename = (char*) realloc(structInfo->filename, i*sizeof(char));
       outFile = safe_fopen(strcat(structInfo->filename, ".sum"), "w");
-    else
+    } else
       outFile = stdout;
 
     fprintf(outFile, "#ResNum\t3-10-Helix\talpha-helix\tPI-Helix\tparallel-Sheet\tantip.-Sheet\tTurn\n");
@@ -18123,17 +18133,26 @@
     /*
      *  process the remaining arguments
      */
-    if ( (buffer = argumentStackKeyToString(argumentStackPointer, "principal", NULL)) != NULL ) {
+    if ( argumentStackContains(argumentStackPointer, "principal")  ) {
       vectorInfo->mode = VECTOR_PRINCIPAL_X;
-      if (strcmp(buffer, "x") == 0) 
-	vectorInfo->mode = VECTOR_PRINCIPAL_X;
-      else if (strcmp(buffer, "y") == 0) 
-	vectorInfo->mode = VECTOR_PRINCIPAL_Y;
-      else if (strcmp(buffer, "z") == 0) 
-	vectorInfo->mode = VECTOR_PRINCIPAL_Z;
-      safe_free(buffer);
-    } 
-    else if (argumentStackContains(argumentStackPointer, "dipole")) 
+      /* DRR - Check the next argument for x, y, or z. If none of these 
+       *       are present, put the argument back on the stack.
+       */
+      buffer = getArgumentString(argumentStackPointer, NULL);
+      if (buffer!=NULL) {
+        if (strcmp(buffer, "x") == 0) { 
+          vectorInfo->mode = VECTOR_PRINCIPAL_X; 
+          safe_free(buffer);
+        } else if (strcmp(buffer, "y") == 0) {
+	  vectorInfo->mode = VECTOR_PRINCIPAL_Y;
+          safe_free(buffer);
+        } else if (strcmp(buffer, "z") == 0) {
+	  vectorInfo->mode = VECTOR_PRINCIPAL_Z;
+          safe_free(buffer);
+        } else
+          pushBottomStack(argumentStackPointer, buffer); 
+      }
+    } else if (argumentStackContains(argumentStackPointer, "dipole")) 
       vectorInfo->mode = VECTOR_DIPOLE;
     else if (argumentStackContains(argumentStackPointer, "box")) 
       vectorInfo->mode = VECTOR_BOX;
@@ -18263,6 +18282,11 @@
        vectorInfo->mode == VECTOR_CORRIRED ||
        vectorInfo->mode == VECTOR_MASK){
       buffer = getArgumentString(argumentStackPointer, NULL);
+      if (buffer==NULL) {
+        fprintf(stdout,"Error: vector: specified vector mode requires a second mask.\n");
+        freeTransformVectorMemory(action);
+        return -1;
+      }
       vectorInfo->mask2 = processAtomMask(buffer, action->state);
       safe_free(buffer);
     }
@@ -18885,7 +18909,7 @@
   /*
    *  USAGE:
    *
-   *    watershell mask [lower <lower cut>] [upper <upper cut>] [noimage]
+   *    watershell mask filename [lower <lower cut>] [upper <upper cut>] [noimage]
    *
    *  action argument usage:
    *
@@ -18917,11 +18941,25 @@
     action->iarg1 = argumentStackContains(argumentStackPointer, "noimage");
 
     buffer = getArgumentString(argumentStackPointer, NULL);
+    if (buffer==NULL) {
+      fprintf(stdout,"ERROR: WATERSHELL: Solute mask must be specified.\n");
+      return -1;
+    }
+
     info->soluteMask = processAtomMask(buffer, action->state);
+    if (info->soluteMask==NULL) {
+      fprintf(stdout,"ERROR: WATERSHELL: Solute mask %s corresponds to 0 atoms.\n",buffer);
+      safe_free(buffer);
+      return -1;
+    }
     safe_free(buffer);
 
     info->filename = getArgumentString(argumentStackPointer, NULL);
-    
+    if (info->filename==NULL) {
+      fprintf(stdout,"ERROR: WATERSHELL: Output filename must be specified.\n");
+      return -1;
+    }
+
     info->lowerCutoff = 
       argumentStackKeyToDouble(argumentStackPointer, "lower", info->lowerCutoff);
     info->upperCutoff = 
@@ -18932,6 +18970,17 @@
       info->solventMask = processAtomMask(buffer, action->state);
     else
       info->solventMask = processAtomMask(":WAT", action->state);
+    if (info->solventMask==NULL) {
+      if (buffer!=NULL)
+        fprintf(stdout,"ERROR: WATERSHELL: Solvent mask %s corresponds to 0 atoms.\n",buffer);
+      else {
+        fprintf(stdout,"ERROR: WATERSHELL: Default solvent mask :WAT corresponds to 0 atoms.\n");
+        fprintf(stdout,
+          "                   Solvent mask can be specified as the third argument.\n");
+      }
+      safe_free(buffer);
+      return -1;
+    }
     safe_free(buffer);
 
     action->carg1 = (void *) info;
diff -urN amber11.orig/AmberTools/src/ptraj/ptraj.c amber11/AmberTools/src/ptraj/ptraj.c
--- amber11.orig/AmberTools/src/ptraj/ptraj.c	2011-04-14 15:30:19.000000000 +0300
+++ amber11/AmberTools/src/ptraj/ptraj.c	2011-10-25 15:01:28.224288536 +0300
@@ -2059,7 +2059,8 @@
   int start = 1;
   int stop = 1;
   int frame_lines, title_size, seekable;
-  long long int file_size, frame_size;
+  int maxi, sizeFound; // For large gzip file calc, Amber Traj
+  long long int file_size, frame_size, tmpfsize;
   long int endoffset;
   float *binposScratch;
   FILE *fp;
@@ -2439,16 +2440,54 @@
       file_size = file_size - frame_size; // Subtract title size from file total size.
       frame_size = (long long int) trajInfo->frameSize;
       trajInfo->Nframes = (int) (file_size / frame_size);
+
+      // Frame calculation for large gzip files
+      // If uncompressed size is less than compressed size, uncompressed
+      // size is likely > 4GB.
+      if (trajInfo->compressType == 1 && file_size < (long long int)frame_stat.st_size) {
+        // Since this is gzip compressed, if the file_size % frame size != 0, 
+        // it could be that the uncompressed filesize > 4GB. Since 
+        //   ISIZE = uncompressed % 2^32, 
+        // try ((file_size + (2^32 * i)) % frame_size) and see if any are 0.
+        if ( (file_size % frame_size) != 0) {
+          // Determine the maximum number of iterations to try based on the
+          // fact that Amber trajectories typically compress about 3x with
+          // gzip. If the number of frames cannot accurately be calculated 
+          // use the max estimated file size to estimate # frames so that
+          // ptraj actions allocate enough memory.
+          tmpfsize = (long long int) frame_stat.st_size;
+          tmpfsize *= 4;
+          tmpfsize = (tmpfsize - file_size) / 4294967296LL;
+          maxi = (int) tmpfsize;
+          maxi++;
+          if (prnlev>1)
+            printf("\tLooking for uncompressed gzip size > 4GB, %i iterations.\n",maxi);
+          tmpfsize = 0;
+          sizeFound=0;
+          for (i = 0; i < maxi; i++ ) {
+            tmpfsize = (4294967296LL * i) + file_size;
+            if ( (tmpfsize % frame_size) == 0) {sizeFound=1; break;}
+          }
+          if (sizeFound) {
+            printf("Warning: Cannot accurately determine # of frames in gzipped trajectory %s.\n",
+                   filename);
+            printf("         This usually indicates the trajectory is corrupted.\n");
+            printf("         Ptraj will attempt to estimate the correct number of frames.\n");
+          }
+          file_size = tmpfsize; 
+          trajInfo->Nframes = (int) (file_size / frame_size);
+        }
+      }
+
       if (prnlev>0) fprintf(stdout,"    File has %i frames.\n",trajInfo->Nframes);
       if ( (file_size % frame_size) == 0 ) {
 	seekable = 1;
-	stop = trajInfo->Nframes;
       } else {
-	stop = -1;
 	seekable = 0;
 	fprintf(stderr, "%s: Could not predict number of frames for AMBER trajectory file: %s\n", ROUTINE, filename);
 	fprintf(stderr, "\tIf this is not a compressed file then there is a problem\n");
       }
+    stop = trajInfo->Nframes;
 
     actualAtoms = totalAtoms;
     //trajInfo->frameSize = frame_size;
@@ -2569,6 +2608,8 @@
   trajInfo->isBox = isBox;
   trajInfo->isVelocity = isVelocity;
   trajInfo->type = type;
+  // If compressed, no seek possible
+  if (trajInfo->compressType > 0) trajInfo->seekable=0;
 
   // This should eventually get its own structure
   if (type==COORD_CHARMM_TRAJECTORY) 
diff -urN amber11.orig/AmberTools/src/rism/amber_rism_interface.f amber11/AmberTools/src/rism/amber_rism_interface.f
--- amber11.orig/AmberTools/src/rism/amber_rism_interface.f	2011-04-14 15:30:19.000000000 +0300
+++ amber11/AmberTools/src/rism/amber_rism_interface.f	2011-10-25 15:01:28.241288520 +0300
@@ -815,7 +815,6 @@
         !redo calculation with charges off
         call rism3d_unsetCharges(rism_3d)
         call rism3d_solve(rism_3d,rismprm%saveprogress,rismprm%progress,rismprm%maxstep,rismprm%tolerance)
-        call rism3d_resetCharges(rism_3d)
 
         !setup memory space
         apol_exchem => apol_mpi_buffer(1:rism_3d%solv%natom) 
@@ -825,6 +824,7 @@
         apol_exchem = rism3d_exchem(rism_3d,rismprm%asympCorr)*KB*rism_3d%solv%temperature
         apol_exchemGF = rism3d_exchemGF(rism_3d,rismprm%asympCorr)*KB*rism_3d%solv%temperature
         
+        call rism3d_resetCharges(rism_3d)
         !parallel communication
 #ifdef MPI
 #  ifdef USE_MPI_IN_PLACE
diff -urN amber11.orig/AmberTools/src/rism/rism1d_potential_c.f amber11/AmberTools/src/rism/rism1d_potential_c.f
--- amber11.orig/AmberTools/src/rism/rism1d_potential_c.f	2011-04-14 15:30:19.000000000 +0300
+++ amber11/AmberTools/src/rism/rism1d_potential_c.f	2011-10-25 15:01:28.242288519 +0300
@@ -124,7 +124,7 @@
     type(solvMDL), intent(in) :: mdl
     _REAL_, intent(in) :: density
     integer :: oldnv
-    integer :: iv1, iv2, ivv, iv, iat, imlt
+    integer :: iv1, iv2, ivv, iv, iat, imlt, isp
 
     !increment scalar variables
     oldnv = this%nv
@@ -139,7 +139,10 @@
     this%nat(this%nsp) = mdl%ntype
     this%mta => safemem_realloc(this%mta, max(ubound(this%mta,1),mdl%ntype), this%nsp)
     this%mta(1:mdl%ntype,this%nsp) = mdl%multi 
-
+    !zero out multiplicity for unused indicies in previously add species
+    do isp=1, this%nsp
+       this%mta(this%nat(isp)+1:,isp)=0
+    end do
     this%mtv => safemem_realloc(this%mtv,this%nv)
     this%mtv(oldnv+1:this%nv) = mdl%multi
 
@@ -178,7 +181,7 @@
     this%rminv => safemem_realloc(this%rminv,this%nv)
     this%rminv(oldnv+1:this%nv) = mdl%rmin/2d0
 
-    this%rma => safemem_realloc(this%rma,3,ubound(this%mta,1),maxval(this%nat),this%nsp)
+    this%rma => safemem_realloc(this%rma,3,maxval(this%mta),maxval(this%nat),this%nsp,.true.,.false.)
     iv = 0
     do iat=1,this%nat(this%nsp)
        do imlt = 1, this%mta(iat,this%nsp)
@@ -186,7 +189,7 @@
           this%rma(:,imlt,iat,this%nsp) = mdl%coord(:,iv)
        end do
     end do
-    this%wlmvv => safemem_realloc(this%wlmvv,ubound(this%mta,1),this%nv,this%nv)
+    this%wlmvv => safemem_realloc(this%wlmvv,maxval(this%mta),this%nv,this%nv)
     call intramolecular_dist(this)
 
     !data arrays in the original 1D-RISM code had an offset of 0 in the first index (like C).
@@ -269,7 +272,7 @@
          r,rs, rs6,ri6, usr,usra
 
     _REAL_ ::  d0x(this%nv), d0y(this%nv), d1z(this%nv), &
-         wkvv(this%nv,this%nv)
+         wkvv(this%nv,this%nv), erfc_test
     _REAL_, external :: erfc
 
 
@@ -478,10 +481,21 @@
              r = (ir-1)*this%dr
              k = (ir-1)*this%dk
              qvv = charge*this%qspv(iv1)*this%qspv(iv2)
-             this%hlrvv(ir,ivv) = -qvv/this%dielconst &
-                  * 0.5d0 * exp((this%smear*this%kappa/2.d0)**2) &
-                  * ( exp(-this%kappa*r)*erfc(this%kappa*this%smear/2.d0 - r/this%smear) &
-                  - exp(this%kappa*r)*erfc(this%kappa*this%smear/2.d0 + r/this%smear) )
+             !for large grids with short Debye lengths the positive
+             !exponent can overflow.  erfc is already zero at this
+             !point.  So we test erfc to see if we can avoid the
+             !exponent
+             erfc_test = erfc(this%kappa*this%smear/2.d0 + r/this%smear)
+             if(erfc_test > sqrt(tiny(1d0)))then
+                this%hlrvv(ir,ivv) = -qvv/this%dielconst &
+                     * 0.5d0 * exp((this%smear*this%kappa/2.d0)**2) &
+                     * ( exp(-this%kappa*r)*erfc(this%kappa*this%smear/2.d0 - r/this%smear) &
+                     - exp(this%kappa*r)*erfc_test )
+             else
+                this%hlrvv(ir,ivv) = -qvv/this%dielconst &
+                     * 0.5d0 * exp((this%smear*this%kappa/2.d0)**2) &
+                     * ( exp(-this%kappa*r)*erfc(this%kappa*this%smear/2.d0 - r/this%smear))
+             end if
              this%hlkvv(ir,ivv) = -qvv/this%dielconst &
                   * 4.d0*pi*exp(-(0.5d0*this%smear*k)**2) * k/(k**2+this%kappa**2)
           enddo
diff -urN amber11.orig/AmberTools/src/rism/rism3d_c.f amber11/AmberTools/src/rism/rism3d_c.f
--- amber11.orig/AmberTools/src/rism/rism3d_c.f	2011-04-14 15:30:19.000000000 +0300
+++ amber11/AmberTools/src/rism/rism3d_c.f	2011-10-25 15:01:28.245288516 +0300
@@ -1086,26 +1086,26 @@
   !reallocate arrays that require preservation of their contents
 #if defined(MPI)
   this%cuvWRK => safemem_realloc(this%cuvWRK,this%grid%nr(1),this%grid%nr(2),this%grid%nr(3),&
-       this%solv%natom,this%NVec,.true.)
+       this%solv%natom,this%NVec,.true.,.true.)
   if(rism3d_solu_charged(this%solu))then
      this%oldcuvChg => safemem_realloc(this%oldcuvChg,this%grid%nr(1),this%grid%nr(2),&
-          this%grid%nr(3),this%solv%natom,this%ncuvsteps,.true.)
+          this%grid%nr(3),this%solv%natom,this%ncuvsteps,.true.,.true.)
      this%oldcuv => this%oldcuvChg
   else
      this%oldcuvNoChg => safemem_realloc(this%oldcuvNoChg,this%grid%nr(1),this%grid%nr(2),&
-          this%grid%nr(3),this%solv%natom,this%ncuvsteps,.true.)
+          this%grid%nr(3),this%solv%natom,this%ncuvsteps,.true.,.true.)
      this%oldcuv => this%oldcuvNoChg
   end if
 #else
   this%cuvWRK => safemem_realloc(this%cuvWRK,this%grid%nr(1),this%grid%nr(2),this%grid%nr(3),&
-       this%solv%natom,this%NVec,.true.)
+       this%solv%natom,this%NVec,.true.,.true.)
   if(rism3d_solu_charged(this%solu))then
      this%oldcuvChg => safemem_realloc(this%oldcuvChg,this%grid%nr(1),this%grid%nr(2),&
-          this%grid%nr(3),this%solv%natom,this%ncuvsteps,.true.)
+          this%grid%nr(3),this%solv%natom,this%ncuvsteps,.true.,.true.)
      this%oldcuv => this%oldcuvChg
   else
      this%oldcuvNoChg => safemem_realloc(this%oldcuvNoChg,this%grid%nr(1),this%grid%nr(2),&
-          this%grid%nr(3),this%solv%natom,this%ncuvsteps,.true.)
+          this%grid%nr(3),this%solv%natom,this%ncuvsteps,.true.,.true.)
      this%oldcuv => this%oldcuvNoChg
   end if
 !!$  this%oldcuv => safemem_realloc(this%oldcuv,this%grid%nr(1),this%grid%nr(2),&
diff -urN amber11.orig/AmberTools/src/rism/rism3d_opendx.f amber11/AmberTools/src/rism/rism3d_opendx.f
--- amber11.orig/AmberTools/src/rism/rism3d_opendx.f	2011-04-14 15:30:19.000000000 +0300
+++ amber11/AmberTools/src/rism/rism3d_opendx.f	2011-10-25 15:01:28.205288553 +0300
@@ -49,9 +49,11 @@
      module procedure readDXHeader_file, readDXHeader_unit
   end interface readDXHeader
 
+#if 0
   interface readDX
      module procedure readDX_3D,readDX_3D_data, readDX_1D, readDX_1D_data
   end interface readDX
+#endif
 contains
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !!!Reads in a OpenDX file header.  Expects
@@ -160,6 +162,7 @@
     end if
   end subroutine readDXHeader_unit
 
+#if 0
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !!!Reads in a OpenDX file, returning data in a pointer, allocating necessary memory.  Expects
 !!!a 3D, regularly spaced grid in ASCII format.
@@ -311,6 +314,7 @@
     call readDX_data(unit,data,ndata,npos)
     close(unit)
   end subroutine readDX_1D_data
+#endif
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !!!write to file in open DX format for use in VMD.  When writing in
diff -urN amber11.orig/AmberTools/src/rism/rism3d.snglpnt.nab amber11/AmberTools/src/rism/rism3d.snglpnt.nab
--- amber11.orig/AmberTools/src/rism/rism3d.snglpnt.nab	2011-04-14 15:30:19.000000000 +0300
+++ amber11/AmberTools/src/rism/rism3d.snglpnt.nab	2011-10-25 15:01:28.243288518 +0300
@@ -632,6 +632,7 @@
 mm_options(sprintf("centering=%d, zerofrc=%d, apply_rism_force=%d",rismOpt.centering,rismOpt.zerofrc,rismOpt.apply_rism_force));
 mm_options(sprintf("polarDecomp=%d",rismOpt.polarDecomp));
 mm_options(sprintf("ntwrism=%d, verbose=%d, progress=%d",rismOpt.ntwrism,rismOpt.verbose,rismOpt.progress));
+mm_options(sprintf("asympCorr=%d",rismOpt.asympcorr));
 mm_options(sprintf("saveprogress=%d",rismOpt.saveprogress));
 mme_init( m, NULL, "::Z", p_xyz, NULL);
 
diff -urN amber11.orig/AmberTools/src/rism/safemem.f amber11/AmberTools/src/rism/safemem.f
--- amber11.orig/AmberTools/src/rism/safemem.f	2011-04-14 15:30:19.000000000 +0300
+++ amber11/AmberTools/src/rism/safemem.f	2011-10-25 15:01:28.247288515 +0300
@@ -46,13 +46,13 @@
        !logical :: Current number of bytes of logical memory allocated
        !char    :: Current number of bytes of character memory allocated
        !total   :: Current number of bytes of total memory allocated
-       integer*8 :: int,real,logical,char,total
+       integer*8 :: int=0,real=0,logical=0,char=0,total=0
        !maxint     :: Current number of bytes of integer memory allocated
        !maxreal    :: Current number of bytes of real memory allocated
        !maxlogical :: Current number of bytes of logical memory allocated
        !maxchar    :: Current number of bytes of character memory allocated
        !max        :: Current number of bytes of total memory allocated
-       integer*8 :: maxint,maxreal,maxlogical,maxchar,max
+       integer*8 :: maxint=0,maxreal=0,maxlogical=0,maxchar=0,max=0
     end type memTracker
 
     !BYTES_PER_GIGABYTES :: used to convert between bytes and GB
@@ -352,6 +352,9 @@
 !!!   n4    :: size of the fourth dimension
 !!!   preserve :: boolean indicating whether or not to preserve the contents of 
 !!!           the array (optional)
+!!!   center :: when preserving the contents of the array, center the 
+!!!             data in the array instead of preserving the value at each
+!!!             index 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     FUNCTION safemem_realloc_4d_real(p, n1, n2, n3, n4, preserve, center)
       implicit none
@@ -369,7 +372,7 @@
 #endif /*RISM_DEBUG*/
       prsrv = .true.
       if(present(preserve)) prsrv=preserve
-      cntr = .true.
+      cntr = .false.
       if(present(center)) cntr=center
 
       if(.not.prsrv .and. associated(p)) then
@@ -499,6 +502,9 @@
 !!!   n5    :: size of the fifth dimension
 !!!   preserve :: boolean indicating whether or not to preserve the contents of 
 !!!           the array (optional)
+!!!   center :: when preserving the contents of the array, center the 
+!!!             data in the array instead of preserving the value at each
+!!!             index 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     FUNCTION safemem_realloc_5d_real(p, n1, n2, n3, n4, n5, preserve, center)
       implicit none
@@ -516,7 +522,7 @@
 #endif /*RISM_DEBUG*/
       prsrv = .true.
       if(present(preserve)) prsrv = preserve
-      cntr=.true.
+      cntr=.false.
       if(present(center)) cntr = center
 
       if(.not.prsrv .and. associated(p)) then
@@ -1006,7 +1012,7 @@
       !work around for old gfortran bug
       !  if(associated(p)) deallocate(p,STAT=safemem_dealloc_pointer_1d_logical)
       if(associated(p))then
-         call memadd_c(totalmem,-product(ubound(p))*kind(p))
+         call memadd_l(totalmem,-product(ubound(p))*kind(p))
          deallocate(p,STAT=temp)
       end if
       safemem_dealloc_pointer_1d_logical = temp
diff -urN amber11.orig/AmberTools/src/sqm/qm2_energy.f amber11/AmberTools/src/sqm/qm2_energy.f
--- amber11.orig/AmberTools/src/sqm/qm2_energy.f	2011-03-11 20:43:13.000000000 +0200
+++ amber11/AmberTools/src/sqm/qm2_energy.f	2011-10-25 15:01:28.206288552 +0300
@@ -277,7 +277,7 @@
 
       ! Add PM6 corrections to Heat of Formation
       if (qmmm_nml%qmtheory%PM6) then
-         escf = escf + pm6_correction()
+         escf = escf + hofCorrection()
       end if
 
       call timer_stop(TIME_QMMMENERGYSCF)
diff -urN amber11.orig/AmberTools/src/sqm/qm2_get_qm_forces.f amber11/AmberTools/src/sqm/qm2_get_qm_forces.f
--- amber11.orig/AmberTools/src/sqm/qm2_get_qm_forces.f	2011-03-11 20:43:13.000000000 +0200
+++ amber11/AmberTools/src/sqm/qm2_get_qm_forces.f	2011-10-25 15:01:28.208288550 +0300
@@ -250,7 +250,7 @@
    ! --------------------------------------------
    if (qmmm_nml%qmtheory%PM6) then
       natom = qmmm_struct%nquant_nlink
-      call pm6_correction(natom, dxyzqm)
+      call hofCorrectionGradient(natom, dxyzqm)
    end if
 
    if(qmmm_nml%peptide_corr) then
diff -urN amber11.orig/AmberTools/src/sqm/qm2_pm6_hof_module.f amber11/AmberTools/src/sqm/qm2_pm6_hof_module.f
--- amber11.orig/AmberTools/src/sqm/qm2_pm6_hof_module.f	2011-03-11 20:43:13.000000000 +0200
+++ amber11/AmberTools/src/sqm/qm2_pm6_hof_module.f	2011-10-25 15:01:28.209288549 +0300
@@ -21,8 +21,7 @@
   implicit none
 
   private
-  public :: pm6_correction
-  public :: corInfoType, print
+  public :: corInfoType, print, hofCorrection, hofCorrectionGradient
   public :: cct, nsp2
   public :: strlen
 
@@ -30,11 +29,6 @@
      module procedure printCorInfoType
   end interface
 
-  interface pm6_correction
-     module procedure hofCorrection
-     module procedure hofCorrectionGradient
-  end interface
-
   ! Data type collecting information on PM6 HOF corrections
   type corInfoType
      logical :: inUse  ! correction in use ?
diff -urN amber11.orig/AmberTools/test/cpptraj/Test_General/RunTest.sh amber11/AmberTools/test/cpptraj/Test_General/RunTest.sh
--- amber11.orig/AmberTools/test/cpptraj/Test_General/RunTest.sh	2011-04-14 15:30:21.000000000 +0300
+++ amber11/AmberTools/test/cpptraj/Test_General/RunTest.sh	2011-10-25 15:01:28.215288544 +0300
@@ -25,7 +25,7 @@
 DoTest PhiPsi.dat.save PhiPsi.dat
 DoTest test.crd.save test.crd
 DoTest a1.dat.save a1.dat
-DoTest test.rst7.103.save Restart/test.rst7.103
+DoTest test.rst7.103.save Restart/test.rst7.123
 DoTest test.nc.save test.nc
 DoTest r4.dat.save r4.dat
 DoTest a2.dat.gz.save a2.dat.gz
diff -urN amber11.orig/AmberTools/test/mmpbsa_py/11_3D-RISM/FINAL_RESULTS_MMPBSA.dat.save amber11/AmberTools/test/mmpbsa_py/11_3D-RISM/FINAL_RESULTS_MMPBSA.dat.save
--- amber11.orig/AmberTools/test/mmpbsa_py/11_3D-RISM/FINAL_RESULTS_MMPBSA.dat.save	2011-04-14 15:30:22.000000000 +0300
+++ amber11/AmberTools/test/mmpbsa_py/11_3D-RISM/FINAL_RESULTS_MMPBSA.dat.save	2011-10-25 15:01:28.247288515 +0300
@@ -1,4 +1,4 @@
-| Run on Tue Mar  8 15:33:18 PST 2011
+| Run on Fri Oct  7 17:12:46 EDT 2011
 
 |Input file:
 |--------------------------------------------------------------
@@ -32,14 +32,14 @@
 -------------------------------------------------------------------------------
 VDWAALS                     -4.6380                0.3291              0.2327
 EEL                       -425.3449                4.5554              3.2212
-ERISM                     -164.0501                3.9611              2.8009
-   POLAR                  -219.8772                4.0588              2.8700
-   APOLAR                   55.8271                0.0977              0.0691
+ERISM                     -164.0499                3.9609              2.8008
+   POLAR                  -235.0781                4.4994              3.1816
+   APOLAR                   71.0282                0.5386              0.3808
 
 G gas                     -429.9829                4.2263              2.9885
-G solv                    -164.0501                3.9611              2.8009
+G solv                    -164.0499                3.9609              2.8008
 
-TOTAL                     -594.0331                0.2653              0.1876
+TOTAL                     -594.0329                0.2655              0.1877
 
 
 Receptor:
@@ -48,8 +48,8 @@
 VDWAALS                     -3.9302                0.2702              0.1910
 EEL                       -254.1833                4.5196              3.1959
 ERISM                      -79.2939                3.0131              2.1306
-   POLAR                  -118.9799                3.1556              2.2314
-   APOLAR                   39.6860                0.1425              0.1008
+   POLAR                  -126.9784                3.4713              2.4546
+   APOLAR                   47.6845                0.4582              0.3240
 
 G gas                     -258.1136                4.2494              3.0048
 G solv                     -79.2939                3.0131              2.1306
@@ -62,12 +62,12 @@
 -------------------------------------------------------------------------------
 VDWAALS                     -0.7056                0.0606              0.0428
 EEL                       -170.5448                1.1864              0.8389
-ERISM                      -85.2510                0.0719              0.0509
-   POLAR                  -101.3724                0.0489              0.0346
-   APOLAR                   16.1214                0.0230              0.0163
+ERISM                      -85.2510                0.0719              0.0508
+   POLAR                  -108.6574                0.0351              0.0249
+   APOLAR                   23.4064                0.1070              0.0757
 
 G gas                     -171.2505                1.2470              0.8818
-G solv                     -85.2510                0.0719              0.0509
+G solv                     -85.2510                0.0719              0.0508
 
 TOTAL                     -256.5015                1.1751              0.8309
 
@@ -77,15 +77,15 @@
 -------------------------------------------------------------------------------
 VDWAALS                     -0.0021                0.0017              0.0012
 EEL                         -0.6167                1.2222              0.8642
-ERISM                        0.4947                1.0198              0.7211
-   POLAR                     0.4750                0.9520              0.6732
-   APOLAR                    0.0198                0.0678              0.0480
+ERISM                        0.4950                1.0196              0.7210
+   POLAR                     0.5577                0.9929              0.7021
+   APOLAR                   -0.0627                0.0267              0.0189
 
 DELTA G gas                 -0.6189                1.2239              0.8654
-DELTA G solv                 0.4947                1.0198              0.7211
+DELTA G solv                 0.4950                1.0196              0.7210
 
 
- DELTA G binding =             -0.1241 +/-      0.2041                 0.1443
+ DELTA G binding =             -0.1239 +/-      0.2043                 0.1445
 
 
 -------------------------------------------------------------------------------
@@ -99,8 +99,8 @@
 VDWAALS                     -4.6380                0.3291              0.2327
 EEL                       -425.3449                4.5554              3.2212
 ERISM                     -187.9280                4.0576              2.8691
-   POLAR                  -221.3631                4.0723              2.8795
-   APOLAR                   33.4351                0.0147              0.0104
+   POLAR                  -236.5648                4.5137              3.1916
+   APOLAR                   48.6367                0.4561              0.3225
 
 G gas                     -429.9829                4.2263              2.9885
 G solv                    -187.9280                4.0576              2.8691
@@ -113,12 +113,12 @@
 -------------------------------------------------------------------------------
 VDWAALS                     -3.9302                0.2702              0.1910
 EEL                       -254.1833                4.5196              3.1959
-ERISM                      -94.4232                3.0803              2.1781
-   POLAR                  -119.8020                3.1352              2.2169
-   APOLAR                   25.3788                0.0549              0.0388
+ERISM                      -94.4233                3.0804              2.1781
+   POLAR                  -127.8003                3.4511              2.4403
+   APOLAR                   33.3771                0.3707              0.2621
 
 G gas                     -258.1136                4.2494              3.0048
-G solv                     -94.4232                3.0803              2.1781
+G solv                     -94.4233                3.0804              2.1781
 
 TOTAL                     -352.5368                1.1691              0.8267
 
@@ -128,14 +128,14 @@
 -------------------------------------------------------------------------------
 VDWAALS                     -0.7056                0.0606              0.0428
 EEL                       -170.5448                1.1864              0.8389
-ERISM                      -94.0188                0.0538              0.0381
-   POLAR                  -102.0695                0.0274              0.0194
-   APOLAR                    8.0506                0.0264              0.0187
+ERISM                      -94.0185                0.0541              0.0383
+   POLAR                  -109.3544                0.0565              0.0399
+   APOLAR                   15.3358                0.1106              0.0782
 
 G gas                     -171.2505                1.2470              0.8818
-G solv                     -94.0188                0.0538              0.0381
+G solv                     -94.0185                0.0541              0.0383
 
-TOTAL                     -265.2693                1.1932              0.8437
+TOTAL                     -265.2690                1.1929              0.8435
 
 
 Differences (Complex - Receptor - Ligand):
@@ -143,15 +143,15 @@
 -------------------------------------------------------------------------------
 VDWAALS                     -0.0021                0.0017              0.0012
 EEL                         -0.6167                1.2222              0.8642
-ERISM                        0.5140                1.0311              0.7291
-   POLAR                     0.5083                0.9645              0.6820
-   APOLAR                    0.0057                0.0666              0.0471
+ERISM                        0.5137                1.0314              0.7293
+   POLAR                     0.5900                1.0061              0.7114
+   APOLAR                   -0.0762                0.0253              0.0179
 
 DELTA G gas                 -0.6189                1.2239              0.8654
-DELTA G solv                 0.5140                1.0311              0.7291
+DELTA G solv                 0.5137                1.0314              0.7293
 
 
- DELTA G binding =             -0.1049 +/-      0.1928                 0.1363
+ DELTA G binding =             -0.1051 +/-      0.1926                 0.1362
 
 
 -------------------------------------------------------------------------------
diff -urN amber11.orig/AmberTools/test/nab/rism3d.snglpnt.2.out.check amber11/AmberTools/test/nab/rism3d.snglpnt.2.out.check
--- amber11.orig/AmberTools/test/nab/rism3d.snglpnt.2.out.check	2011-04-14 15:30:23.000000000 +0300
+++ amber11/AmberTools/test/nab/rism3d.snglpnt.2.out.check	2011-10-25 15:01:28.248288514 +0300
@@ -32,6 +32,7 @@
 	mm_options:  ntwrism=1
 	mm_options:  verbose=0
 	mm_options:  progress=0
+	mm_options:  asympCorr=1
 	mm_options:  saveprogress=0
 |3D-RISM thermodynamic data key
 |solute_epot [kcal/mol] :            Total               LJ          Coulomb             Bond            Angle         Dihedral           H-Bond            LJ-14       Coulomb-14       Restraints          3D-RISM
@@ -74,13 +75,13 @@
 3D-RISM processing complete.
 
 |Timing summary:
-|              Initialize       0.042
+|              Initialize       0.046
 |             Molec. Dyn.       0.000
 |             Normal Mode       0.000
 |             Conj. Grad.       0.000
 |                  Newton       0.000
 |-------------------------
-|                   Total       0.042
+|                   Total       0.046
 
 |1st derivative timing summary:
 |             constraints       0.000
@@ -91,7 +92,7 @@
 |                 nonbond       0.000
 |               gen. Born       0.000
 |       Poisson Boltzmann       0.000
-|                 3D-RISM      10.167
+|                 3D-RISM       8.211
 |                   Other       0.000
 |-------------------------
-|                   Total      10.167
+|                   Total       8.212
diff -urN amber11.orig/AmberTools/test/nab/rism3d.snglpnt.out.check amber11/AmberTools/test/nab/rism3d.snglpnt.out.check
--- amber11.orig/AmberTools/test/nab/rism3d.snglpnt.out.check	2011-04-14 15:30:23.000000000 +0300
+++ amber11/AmberTools/test/nab/rism3d.snglpnt.out.check	2011-10-25 15:01:28.248288514 +0300
@@ -31,6 +31,7 @@
 	mm_options:  ntwrism=0
 	mm_options:  verbose=0
 	mm_options:  progress=0
+	mm_options:  asympCorr=1
 	mm_options:  saveprogress=0
 |3D-RISM thermodynamic data key
 |solute_epot [kcal/mol] :            Total               LJ          Coulomb             Bond            Angle         Dihedral           H-Bond            LJ-14       Coulomb-14       Restraints          3D-RISM
@@ -55,10 +56,10 @@
 rism_volume:  1.93595981E+002
 rism_exNumb:                  -6.38184158E+000 -1.27637588E+001
 rism_exChrg: -3.20562831E-005  5.40924893E+000 -5.40928099E+000
-rism_polar_: -1.31661477E+001 -4.40204657E+001  3.08543179E+001
-rism_apolar:  2.96543934E+001  7.56776450E+001 -4.60232516E+001
-rism_polGF_: -1.34197324E+001 -4.39508092E+001  3.05310768E+001
-rism_apolGF:  1.97398606E+001  7.11960080E+001 -5.14561474E+001
+rism_polar_: -1.47850626E+001 -1.09412958E+001 -3.84376682E+000
+rism_apolar:  3.12733083E+001  4.25984751E+001 -1.13251668E+001
+rism_polGF_: -1.50386472E+001 -1.08716393E+001 -4.16700794E+000
+rism_apolGF:  2.13587755E+001  3.81168382E+001 -1.67580627E+001
 
 Frame: 2 of 3
 solute_epot: -3.41386230E-001  2.56854615E+000 -8.23476129E+001  3.38717466E-002  3.59735615E-001  7.49048012E+000  0.00000000E+000  4.97044290E+000  5.01391318E+001  0.00000000E+000  1.64440184E+001
@@ -68,10 +69,10 @@
 rism_volume:  1.93535504E+002
 rism_exNumb:                  -6.37993323E+000 -1.27596468E+001
 rism_exChrg:  9.31177385E-005  5.40763142E+000 -5.40753830E+000
-rism_polar_: -1.31967073E+001 -4.40467251E+001  3.08500178E+001
-rism_apolar:  2.96407256E+001  7.56851573E+001 -4.60444317E+001
-rism_polGF_: -1.34505678E+001 -4.39769309E+001  3.05263631E+001
-rism_apolGF:  1.97272712E+001  7.12038497E+001 -5.14765785E+001
+rism_polar_: -1.48191496E+001 -1.09491800E+001 -3.86996959E+000
+rism_apolar:  3.12631679E+001  4.25876122E+001 -1.13244443E+001
+rism_polGF_: -1.50730101E+001 -1.08793858E+001 -4.19362428E+000
+rism_apolGF:  2.13497135E+001  3.81063046E+001 -1.67565911E+001
 
 Frame: 3 of 3
 solute_epot: -6.77945271E-001  2.31243593E+000 -8.22812685E+001  9.08980189E-002  3.84286510E-001  7.48894041E+000  0.00000000E+000  4.91668330E+000  5.00274224E+001  0.00000000E+000  1.63826567E+001
@@ -81,21 +82,21 @@
 rism_volume:  1.93468855E+002
 rism_exNumb:                  -6.37771162E+000 -1.27552065E+001
 rism_exChrg:  9.18518609E-005  5.40574837E+000 -5.40565652E+000
-rism_polar_: -1.32423609E+001 -4.40992117E+001  3.08568508E+001
-rism_apolar:  2.96250176E+001  7.57232232E+001 -4.60982056E+001
-rism_polGF_: -1.34968274E+001 -4.40292606E+001  3.05324332E+001
-rism_apolGF:  1.97126459E+001  7.12421943E+001 -5.15295484E+001
+rism_polar_: -1.48691739E+001 -1.09507135E+001 -3.91846037E+000
+rism_apolar:  3.12518306E+001  4.25747250E+001 -1.13228944E+001
+rism_polGF_: -1.51236404E+001 -1.08807624E+001 -4.24287800E+000
+rism_apolGF:  2.13394589E+001  3.80936961E+001 -1.67542372E+001
 
 3D-RISM processing complete.
 
 |Timing summary:
-|              Initialize       0.043
+|              Initialize       0.047
 |             Molec. Dyn.       0.000
 |             Normal Mode       0.000
 |             Conj. Grad.       0.000
 |                  Newton       0.000
 |-------------------------
-|                   Total       0.043
+|                   Total       0.047
 
 |1st derivative timing summary:
 |             constraints       0.000
@@ -106,7 +107,7 @@
 |                 nonbond       0.000
 |               gen. Born       0.000
 |       Poisson Boltzmann       0.000
-|                 3D-RISM      32.639
+|                 3D-RISM      30.135
 |                   Other       0.000
 |-------------------------
-|                   Total      32.639
+|                   Total      30.135
diff -urN amber11.orig/AmberTools/test/nab/Run.prm amber11/AmberTools/test/nab/Run.prm
--- amber11.orig/AmberTools/test/nab/Run.prm	2011-04-14 15:30:23.000000000 +0300
+++ amber11/AmberTools/test/nab/Run.prm	2011-10-25 15:01:28.133288619 +0300
@@ -17,7 +17,7 @@
 
 echo "checking the prmtop file:"
 echo ""
-tail -3209 tprmtop > tprmtop1
+tail -3228 tprmtop > tprmtop1
 ../dacdif tprmtop.check tprmtop1
 
 rm -f prm prm.c tleap.out leap.log tprmcrd prm.out1 tprmtop
diff -urN amber11.orig/AmberTools/test/nab/tprmtop.check amber11/AmberTools/test/nab/tprmtop.check
--- amber11.orig/AmberTools/test/nab/tprmtop.check	2011-03-10 20:12:18.000000000 +0200
+++ amber11/AmberTools/test/nab/tprmtop.check	2011-10-25 15:01:28.135288617 +0300
@@ -1,3 +1,4 @@
+%VERSION  VERSION_STAMP = V0001.000  DATE = 05/09/11  12:59:19                  
 %FLAG TITLE                                                                     
 %FORMAT(20a4)                                                                   
                                                                                 
@@ -541,6 +542,24 @@
   3.14159400E+00  3.14159400E+00  3.14159400E+00  3.14159400E+00  3.14159400E+00
   3.14159400E+00  3.14159400E+00  3.14159400E+00  3.14159400E+00  3.14159400E+00
   3.14159400E+00  3.14159400E+00  3.14159400E+00
+%FLAG SCEE_SCALE_FACTOR                                                         
+%FORMAT(5E16.8)                                                                 
+  1.20000000E+00  1.20000000E+00  1.20000000E+00  1.20000000E+00  1.20000000E+00
+  1.20000000E+00  1.20000000E+00  1.20000000E+00  1.20000000E+00  1.20000000E+00
+  1.20000000E+00  1.20000000E+00  1.20000000E+00  1.20000000E+00  1.20000000E+00
+  1.20000000E+00  1.20000000E+00  1.20000000E+00  1.20000000E+00  1.20000000E+00
+  1.20000000E+00  1.20000000E+00  1.20000000E+00  1.20000000E+00  1.20000000E+00
+  1.20000000E+00  1.20000000E+00  1.20000000E+00  1.20000000E+00  1.20000000E+00
+  0.00000000E+00  0.00000000E+00  0.00000000E+00
+%FLAG SCNB_SCALE_FACTOR                                                         
+%FORMAT(5E16.8)                                                                 
+  2.00000000E+00  2.00000000E+00  2.00000000E+00  2.00000000E+00  2.00000000E+00
+  2.00000000E+00  2.00000000E+00  2.00000000E+00  2.00000000E+00  2.00000000E+00
+  2.00000000E+00  2.00000000E+00  2.00000000E+00  2.00000000E+00  2.00000000E+00
+  2.00000000E+00  2.00000000E+00  2.00000000E+00  2.00000000E+00  2.00000000E+00
+  2.00000000E+00  2.00000000E+00  2.00000000E+00  2.00000000E+00  2.00000000E+00
+  2.00000000E+00  2.00000000E+00  2.00000000E+00  2.00000000E+00  2.00000000E+00
+  0.00000000E+00  0.00000000E+00  0.00000000E+00
 %FLAG SOLTY                                                                     
 %FORMAT(5E16.8)                                                                 
   0.00000000E+00  0.00000000E+00  0.00000000E+00  0.00000000E+00  0.00000000E+00
diff -urN amber11.orig/AmberTools/test/rism1d/check1d amber11/AmberTools/test/rism1d/check1d
--- amber11.orig/AmberTools/test/rism1d/check1d	2011-04-14 15:30:29.000000000 +0300
+++ amber11/AmberTools/test/rism1d/check1d	2011-10-25 15:01:28.249288513 +0300
@@ -85,7 +85,7 @@
 #It is unrealistic to achieve low relative error for these very small number so we only 
 #use the abolute error
 if [ -s $1.bvv.save ]; then
-   $AMBERHOME/AmberTools/test/dacdif -a 1.e-7 $1.bvv.save $1.bvv
+   $AMBERHOME/AmberTools/test/dacdif -a 2.e-7 $1.bvv.save $1.bvv
 fi
 
 /bin/rm -f $1.xvv.delhv0 $1.xvv.delhv0.save $1.xvv.xvv $1.xvv.xvv.save\
diff -urN amber11.orig/AT15_Amber11.py amber11/AT15_Amber11.py
--- amber11.orig/AT15_Amber11.py	2011-04-14 15:30:11.000000000 +0300
+++ amber11/AT15_Amber11.py	2011-10-25 15:01:28.221288538 +0300
@@ -145,4 +145,29 @@
 makefile.close()
       
 # Now we are done
-print "\nAmber 11 patched for AmberTools 1.5. Run\nmake %s\nto build amber11\n" % buildtype
+print "\nAmber 11 patched for AmberTools 1.5. Run\ncd src; make %s\nto build amber11\n" % buildtype
+
+# Now print out details about expected test FAILUREs and
+# errors:
+
+warning = """NOTE: Because PBSA has changed since Amber 11 was released, some
+tests are known to fail and others are known to quit in error. These
+can be safely ignored.
+
+Tests that error: Tests in $AMBERHOME/test/sander_pbsa_frc
+   Run.argasp.min    Run.dadt.min      Run.dgdc.min
+   Run.lysasp.min    Run.polyALA.min   Run.polyAT.min
+   Run.argasp.min    Run.dadt.min      Run.dgdc.min
+   Run.lysasp.min    Run.polyALA.min   Run.polyAT.min
+   Run.argasp.min    Run.dadt.min      Run.dgdc.min
+   Run.lysasp.min    Run.polyALA.min   Run.polyAT.min
+
+Tests that produce possible FAILUREs:
+   cd sander_pbsa_ipb2   && ./Run.110D.min
+   cd sander_pbsa_lpb    && ./Run.lsolver.min (only some of them fail here)
+   cd sander_pbsa_tsr    && ./Run.tsrb.min
+   cd sander_pbsa_decres && ./Run.pbsa_decres
+   mm_pbsa.pl tests 02, 03, and 05
+"""
+
+print warning
diff -urN amber11.orig/dat/antechamber/ATOMTYPE_AMBER.DEF amber11/dat/antechamber/ATOMTYPE_AMBER.DEF
--- amber11.orig/dat/antechamber/ATOMTYPE_AMBER.DEF	2011-04-14 15:30:40.000000000 +0300
+++ amber11/dat/antechamber/ATOMTYPE_AMBER.DEF	2011-10-25 15:01:28.198288560 +0300
@@ -49,7 +49,7 @@
 ATD  F     *   9   1   &
 ATD  Cl    *   17  1   &
 ATD  Br    *   35  1   &
-ATD  F     *   53  1   &
+ATD  I     *   53  1   &
 ATD  P     *   15  &
 ATD  N1    *   7   1   &   
 ATD  NB    *   7   2   *   *   [RG5,AR1.AR2.AR3]         & 
diff -urN amber11.orig/dat/leap/cmd/leaprc.ff99SBnmr amber11/dat/leap/cmd/leaprc.ff99SBnmr
--- amber11.orig/dat/leap/cmd/leaprc.ff99SBnmr	2011-04-14 15:30:40.000000000 +0300
+++ amber11/dat/leap/cmd/leaprc.ff99SBnmr	2011-10-25 15:01:28.239288522 +0300
@@ -74,13 +74,6 @@
 	{ "OL"  "O" "sp3" }
 	{ "AC"  "C" "sp3" }
 	{ "EC"  "C" "sp3" }
-#ildn
-	{ "C3"  "C" "sp3" }
-	{ "C4"  "C" "sp3" }
-	{ "C5"   "C" "sp2" }
-	{ "C6"   "C" "sp2" }
-	{ "NP"   "N" "sp2" }
-	{ "OM"   "O" "sp2" } # guess!!!
 }
 #
 #	Load the main parameter set.
@@ -95,9 +88,9 @@
 #	Load main chain and terminating 
 #	amino acid libraries (i.e. ff94 libs)
 #
-loadOff all_amino94ildn.lib
-loadOff all_aminoct94ildn.lib
-loadOff all_aminont94ildn.lib
+loadOff all_amino94.lib
+loadOff all_aminoct94.lib
+loadOff all_aminont94.lib
 #
 #       Load water and ions
 # 
diff -urN amber11.orig/dat/leap/parm/parm10.dat amber11/dat/leap/parm/parm10.dat
--- amber11.orig/dat/leap/parm/parm10.dat	2011-04-14 15:30:40.000000000 +0300
+++ amber11/dat/leap/parm/parm10.dat	2011-10-25 15:01:28.218288541 +0300
@@ -783,8 +783,10 @@
 OS-P -OS-CT   1    0.25          0.0            -3.         JCC,7,(1986),230
 OS-P -OS-CT   1    1.20          0.0             2.         gg&gt ene.631g*/mp2
 H1-CT-C -O    1    0.80          0.0            -1.         Junmei et al, 1999
+H1-CT-C -O    1    0.00          0.0            -2.         Explicit of wild card X-C-CT-X
 H1-CT-C -O    1    0.08        180.0             3.         Junmei et al, 1999
 H1-CX-C -O    1    0.80          0.0            -1.         Junmei et al, 1999 (was H1-CT-C -O )
+H1-CX-C -O    1    0.00          0.0            -2.         Explicit of wild card X-C-CT-X
 H1-CX-C -O    1    0.08        180.0             3.         Junmei et al, 1999 (was H1-CT-C -O )
 HC-CT-C -O    1    0.80          0.0            -1.         Junmei et al, 1999
 HC-CT-C -O    1    0.00          0.0            -2.         Explicit of wild card X-C-CT-X
diff -urN amber11.orig/dat/leap/parm/parm99.dat amber11/dat/leap/parm/parm99.dat
--- amber11.orig/dat/leap/parm/parm99.dat	2011-04-14 15:30:40.000000000 +0300
+++ amber11/dat/leap/parm/parm99.dat	2011-10-25 15:01:28.221288538 +0300
@@ -546,7 +546,7 @@
 OS-P -OS-CT   1    0.25          0.0            -3.         JCC,7,(1986),230
 OS-P -OS-CT   1    1.20          0.0             2.         gg&gt ene.631g*/mp2
 H1-CT-C -O    1    0.80          0.0            -1.         Junmei et al, 1999
-HC-CT-C -O    1    0.00          0.0            -2.         Explicit of wild card X-C-CT-X
+H1-CT-C -O    1    0.00          0.0            -2.         Explicit of wild card X-C-CT-X
 H1-CT-C -O    1    0.08        180.0             3.         Junmei et al, 1999
 HC-CT-C -O    1    0.80          0.0            -1.         Junmei et al, 1999
 HC-CT-C -O    1    0.00          0.0            -2.         Explicit of wild card X-C-CT-X