diff --git a/CMakeLists.txt b/CMakeLists.txt
index d6025d112188b8cc8996e1454f813f06bcfadd6a..6a1dbbbca00307d2e84e4c96a3e5efae10340ea4 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -13,7 +13,7 @@ Set(CMAKE_CXX_COMPILER "g++")
 set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} ${OpenMP_C_FLAGS}")  
 set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}")  
 
-#FIND_PACKAGE( OpenMP REQUIRED)  
+FIND_PACKAGE( OpenMP REQUIRED)  
 if(OPENMP_FOUND)  
     set(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} ${OpenMP_EXE_LINKER_FLAGS}")  
     message(STATUS "OpenMP_C_FLAGS = " ${OpenMP_C_FLAGS})  
diff --git a/custom/Custom.cpp b/custom/Custom.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..f74ec5800fa14fd42626fde069c5adc72f5a7c03
--- /dev/null
+++ b/custom/Custom.cpp
@@ -0,0 +1,35 @@
+
+#include "Home.h"
+#include "Util.h"
+#include "ProDataIO.h"
+
+
+void SpecifyEquations1(Table_t &table)
+{
+        int     i;
+        real8   tau;
+ 
+        vector<vector<double> >::iterator   it;
+        int colSeparation =  GetColIDFromTable(table, "separation");
+
+        for(it = table.data.begin(); it != table.data.end(); ){
+
+            if((*it)[colSeparation] < 25){
+                table.data.erase(it);
+                continue;
+            }
+
+            it++;
+        }
+
+        return;
+}
+
+void SpecifyEquations(Table_t &table)
+{
+        int     i, j;
+
+        SpecifyEquations1(table);
+    
+        return;
+}
diff --git a/dev/reg.cpp b/dev/reg.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..b0951ff5f41f2c4daff2e3fe02b64b1485b33888
--- /dev/null
+++ b/dev/reg.cpp
@@ -0,0 +1,161 @@
+
+#include <iostream>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string>
+#include <regex.h>
+#include <assert.h>
+#include <string.h>
+
+using namespace std;
+
+int find_first(string input, string pattern, string &out){
+    regex_t reg;
+    regmatch_t pm[1];
+    int  iret = 0;
+    out = "";
+
+    /*编译正则表达式*/
+    iret = regcomp(&reg, pattern.c_str(), REG_EXTENDED|REG_NEWLINE);
+    if (iret != 0){
+        return -1;
+    }
+
+    iret = regexec(&reg, input.c_str(), 1, pm, 0);
+    if (iret == REG_NOMATCH){
+        out = "";
+        iret = input.length();
+    }else if (iret != 0) {
+        return -2;
+
+    }else{
+        out = input.substr(pm[0].rm_so,pm[0].rm_eo-pm[0].rm_so);
+        iret = pm[0].rm_eo;
+    }
+    regfree(&reg);
+    return iret;
+}
+
+int find_first(char *buff, char *pattern, char *outdata){
+    regex_t reg;
+    regmatch_t pm[1];
+    int  status = 0;
+    /*编译正则表达式*/
+    status = regcomp(&reg, pattern, REG_EXTENDED|REG_NEWLINE);  //扩展正则表达式和识别换行符
+    if (status != 0){    //成功返回0
+        return -1;
+    }
+    status = regexec(&reg, buff, 1, pm, 0);
+    if (status == REG_NOMATCH){
+        printf("no match!\n");
+        status = -1;
+    }
+    else if (status != 0) {
+        return -2;
+    }
+    else if (status == 0) {
+        int i, j;
+        for (i = pm[0].rm_so, j = 0; i < pm[0].rm_eo; i++, j++) {
+            outdata[j] = buff[i];
+        }
+        outdata[i] = '\0';
+    }
+    regfree(&reg);
+    return status;
+
+}
+
+int find_all(char *buff, char *pattern, char result[][20]){   //返回匹配个数
+    regex_t reg;
+    regmatch_t pm[1];
+    int  status = 0;
+    char * p = buff;
+    int count = 0;
+
+    /*编译正则表达式*/
+    status = regcomp(&reg, pattern, REG_EXTENDED|REG_NEWLINE);  //扩展正则表达式和识别换行符
+    if (status != 0){    //成功返回0
+        return -1;
+    }
+    int i = 0, j, k;
+    while((status = regexec(&reg, p, 1, pm, 0)) == 0) {
+        for(j = pm[0].rm_so, k = 0; j < pm[0].rm_eo; j++) {
+            result[i][k++] = p[j];
+        }
+        result[i][k] = '\0';
+        i++;
+        p += pm[0].rm_eo;
+        count++;
+        if (*p == '\0')  break;
+    }
+
+    regfree(&reg);
+    return count;
+
+}
+
+int print_file(const char *file_name, const char *pattern)
+
+{
+    regex_t reg;
+    regmatch_t pm[1];
+    int  status = 0;
+    int count = 0;
+    FILE *fp = fopen(file_name, "r+");
+    assert(fp);
+    char buff[1024] = {0};
+    char output[1024] = {0};
+
+    /*编译正则表达式*/
+    status = regcomp(&reg, pattern, REG_EXTENDED|REG_NEWLINE);  //扩展正则表达式和识别换行符
+    assert(status == 0);
+
+    while(fgets(buff, sizeof(buff), fp)) {     //循环读取文件
+
+        char * p = buff;
+
+        while(1) {
+            status = regexec(&reg, p, 1, pm, 0);
+            if (status == 0) {  //匹配成功
+                count++;
+                strncpy(output, p + pm[0].rm_so, pm[0].rm_eo - pm[0].rm_so);
+                cout<<"匹配:"<<output<<endl;
+                p += pm[0].rm_eo;
+            }
+            else {
+                break;
+            }
+        }
+    }
+    regfree(&reg);
+
+    return count;
+}
+
+int main(int args, char *argv[])
+
+{
+    char result[20][20] = {0};
+    char buf[] = "1231a4568b789c234";
+    char pattern[] = "23";
+
+    char resultfirst[20] = {0};
+
+    if(args != 3)exit(-1);
+
+    find_first(argv[1], argv[2], resultfirst);
+    cout<<strlen(resultfirst) << ":" << resultfirst<<endl;
+    cout << "***************************" <<endl;
+
+    int count = find_all(argv[1], argv[2], result);
+    for (int i = 0; i < count; i++) {
+        cout<<"result:"<<"i="<<i+1<<"----"<<result[i]<<endl;
+    }
+    cout << "***************************" <<endl;
+    count = print_file("test.txt", "[0-9]{5}");
+    cout<<"匹配的个数:"<<count<<endl;
+
+    return 0;
+
+}
+
diff --git a/include/Home.h b/include/Home.h
index aa7b5ef0387c1104cb430227523277c918e242af..10bf58e359fab5a19ba906a7d7b78f6ec7d3bd4d 100644
--- a/include/Home.h
+++ b/include/Home.h
@@ -14,6 +14,26 @@
 #include <string>
 #include <vector>
 #include <list>
+#include <map>
+
+#ifdef _OPENMP
+
+#include <omp.h>
+
+#define INIT_LOCK(a)    omp_init_lock((a))
+#define LOCK(a)         omp_set_lock((a))
+#define UNLOCK(a)       omp_unset_lock((a))
+#define DESTROY_LOCK(a) omp_destroy_lock((a))
+
+#else /* _OPENMP not defined */
+
+#define INIT_LOCK(a)
+#define LOCK(a)
+#define UNLOCK(a)
+#define DESTROY_LOCK(a)
+
+#endif /* end ifdef _OPENMP */
+
 
 #define real8 double
 using namespace std;
@@ -32,7 +52,7 @@ typedef struct {
 typedef struct {
         double  burgMag;
         int     seed;
-        int     type;
+        int     type, nThreads;
         bool    help;
 
         vector<string>  inpFiles, outFiles;
@@ -67,6 +87,7 @@ typedef	enum{
     OPT_TYPE,
     OPT_PRIVATEVALS,
     OPT_AUXFILE,
+    OPT_THREADS,
 	OPT_MAX
 }OPT_t;
 
@@ -75,9 +96,15 @@ typedef	enum{
     FTYPE_PROC_EXTEND_DIS_DDD,
     FTYPE_PROC_EXTEND_DIS_MD,
     FTYPE_SPECIFY_EQUATIONS,
+    FTYPE_GENERATE_DISLOCATION,
 	FTYPE_MAX
 }FTYPE_t;
 
+typedef enum{
+    INT_DATA = 0,
+    DOUBLE_DATA,
+    STRING_DATA
+}DATA_TYPE;
 /*
  *      Define a structure to hold a command line option's id (type),
  *      name, the shortest possible unique abbreviation of the option
@@ -91,11 +118,25 @@ typedef struct {
 } Option_t;
 
 typedef struct {
+        int     type;
+        string  name;
+        string  val;
+}Variable_t;
+
+typedef struct {
+    string                  T, F;
+    int                     i,j,k;
+    double                  solutionTime;
     vector<string>          variables;
+    map<string, string>     aux;
     vector<vector<double> > data;
 } Table_t;
 
 typedef struct {
+    string                  T, F;
+    int                     i,j,k;
+    double                  solutionTime;
+    map<string, string>     aux;
     vector<string>          variables;
     vector<vector<double> > data;
 }LineList_t;
diff --git a/include/MD.h b/include/MD.h
index 32fa5599ed9f47d1cf7459f57c36d3fba9aab27f..59359b09f40c93a5860abc387d4c3d75b5651b88 100644
--- a/include/MD.h
+++ b/include/MD.h
@@ -28,13 +28,14 @@ typedef struct{
 
 typedef struct {
     int     timestep;
-    real8   box[6];
 
+    vector<vector<real8> > box;
     vector<string> bounds;
     vector<string> variables;
     vector<Atom_t> atom;
-}MgData_t;
+}Dump_t;
 
 void HandleExtendedDislocation_MD(InArgs_t *inArgs);
+void GenerateDislocation(InArgs_t *inArgs);
 
 #endif
diff --git a/include/Math.h b/include/Math.h
index 60c835af37a46d22e1b4960f96566dbe967717c4..7c356474e0d639b452d53cc5068bdd73e3113bd5 100644
--- a/include/Math.h
+++ b/include/Math.h
@@ -18,6 +18,38 @@
 #include <algorithm>
 #include <numeric>
 
+#ifndef MAX
+#define MAX(a,b) ((a)>(b)?(a):(b))
+#endif
+
+#ifndef MIN
+#define MIN(a,b) ((a)<(b)?(a):(b))
+#endif
+
+#define SQUARE(a) ((a)*(a))
+#define VECTOR_ADD(a,b)  {(a)[0] += (b)[0]; (a)[1] += (b)[1]; (a)[2] += (b)[2];}
+#define VECTOR_COPY(a,b) {(a)[0] = (b)[0]; (a)[1] = (b)[1]; (a)[2] = (b)[2];}
+#define VECTOR_ZERO(a)   {(a)[0] = 0; (a)[1]  = 0; (a)[2] = 0;}
+#define VECTOR_MULTIPLY(a,x)	{(a)[0] *= (x);    (a)[1] *= (x);    (a)[2] *= (x);}
+
+#define NO_INTERSECTION     0
+#define POINT_INTERSECTION  1
+#define LINE_INTERSECTION   2
+#define PLANE_INTERSECTION  3
+#define EPS1 1.0E-3
+#define EPS2 1.0E-1
+
+int PlanePlaneIntersection(double *n1, double *p1, double *n2, double*p2,
+                           double *vec1, double *vec2);
+int LinePlaneIntersection(double *d, double *p1, double *n, double *p2, 
+                          double *vec1, double *vec2);
+int LineLineIntersection(double *d1, double *p1, double *d2, double *p2,
+                         double *vec1, double *vec2, double *dis);
+int PointPlaneIntersection(double *p1, double *n, double *p2, double *vec, double *dis);
+int PointLineIntersection(double *p1, double *d, double *p2, double *vec, double *dis);
+int PointPointIntersection(double *p1, double *p2, double *vec, double *dis);
+void cross(real8 a[3], real8 b[3], real8 c[3]);
+real8 Normal(real8 a[3]);
 
 void AverageLines(InArgs_t *inArgs);
 
diff --git a/include/ProDataIO.h b/include/ProDataIO.h
index b6eba257cb5d07fab0dafbeda0f7c34378017e34..214215f4a2c42166ee15111fb397c55b0bec6f5b 100644
--- a/include/ProDataIO.h
+++ b/include/ProDataIO.h
@@ -27,7 +27,10 @@ int ReadTecplotNormalData(string &file, Table_t &table, string &secLine);
 void WriteTecplotNormalData(const Table_t &table, const string &file, double precision, string secLine = ""); 
 void    WriteTecplotNormalData(const LineList_t &list, const string &file,
                                double precision = 6);
-int ReadMGDataFile(const string &file, MgData_t &mgdata);
+int ReadDumpFile(const string &file, Dump_t &dum);
 int ReadDataFromMDLogFile(const vector<string> &files, LineList_t &list);
+int WriteDumpFile(const string &file, Dump_t &dum, int precision = 6); 
+int MGToLMPDataFile(const string &file, Dump_t &dum, int precision = 10);
+bool ReadLMPFile(const string file, Dump_t &dum);
 
 #endif
diff --git a/include/Util.h b/include/Util.h
index fd1ed19180995f6494be390ee6a350367bd1df0f..10eab95d08e7c32807ea42f28af45b4c87b1498b 100644
--- a/include/Util.h
+++ b/include/Util.h
@@ -21,6 +21,11 @@
 
 #include "MD.h"
 
+#define DotProduct(vec1,vec2)       \
+            ((vec1[0])*(vec2[0]) +  \
+             (vec1[1])*(vec2[1]) +  \
+             (vec1[2])*(vec2[2]))
+
 #define  NUMBER_DATA    0x0001
 #define  CHAR_DATA      0x0002
 
@@ -43,12 +48,15 @@ void ZImage(real8 boundMin[3], real8 boundMax[3], real8 *x, real8 *y, real8 *z);
 real8 LinearInterpolation(const Curve_t &curve, real8 x, real8 min = -1, real8 max = -1);
 void SwapTable(Table_t &table);
 void SwapLineList(LineList_t &list);
-void CleanMgData(MgData_t &mg);
+void CleanDump(Dump_t &dum);
 void NormalizeVec(real8 vec[3]);
 void StitchTecplotData(vector<Table_t> &tables, Table_t &table, int eigenID = 0);
 int DataType(const string &str);
 void SpecifyEquations_PLTDATA(InArgs_t *inArgs);
 void SpecifyEquations(Table_t &table);
+void FormatVector(real8 vec[3], const char *msg);
 
+void InitList(LineList_t &list);
+void InitTable(Table_t &table);
 
 #endif
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index c146e29b1489922770a44acf2d388b875820c600..bb3283c8aea34d508c52f0969b5927d1ced71e58 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -1,10 +1,15 @@
 Set (SRCS 
     Main.cpp  ReadData.cpp  Util.cpp Math.cpp HandleExtendedDislocation.cpp
-    WriteData.cpp Parse.cpp)
+    WriteData.cpp Parse.cpp GenerateDislocation.cpp 
+    ../custom/Custom.cpp)
+
+set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -O3 -Wno-unknown-pragmas -std=c++11" )
+set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O3 -std=c++11" )
+
+if(OPENMP_FOUND)  
+set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -openmp")
+endif(OPENMP_FOUND)
 
-set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -O3 -Wno-unknown-pragmas -DPARALLEL=1 " )
-set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O3 " )
-set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ")
 add_executable(prodata ${SRCS})
 add_definitions( -DPARALLEL=1)
 
diff --git a/src/GenerateDislocation.cpp b/src/GenerateDislocation.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..0f5ca20f07308e27f87d9660d63919bbf06df676
--- /dev/null
+++ b/src/GenerateDislocation.cpp
@@ -0,0 +1,326 @@
+#include <algorithm>
+#include <functional>
+
+#include "Home.h"
+#include "Util.h"
+#include "ProDataIO.h"
+#include "DDD.h"
+#include "MD.h"
+#include "Math.h"
+
+using namespace std;
+
+void GenerateScrewDislocation(InArgs_t *inArgs)
+{
+    int         index, i;    
+    Dump_t    dum;
+    string      posName("pos"), lineName("line"), gDirName("gdir"), burgMagName("burgMag");
+    string      screwTypeName("stype"), multiDisName("multi");
+    int         screwType = 0, nPairs = 0;
+    double      rsize = 0, delta;
+
+    real8       pos[3] = {0.0, 0.0, 0.0};
+    real8       line[3] = {1.0, 0.0, 0.0};
+    real8       gDir[3] = {0.0, 1.0, 0.0};
+    real8       normal[3];
+    real8       burgMag = 2.556, boundMin[3], boundMax[3], argument;
+
+    if((index = GetValID(inArgs->priVars, posName)) < inArgs->priVars.size()){
+        if(inArgs->priVars[index].vals.size() == 3){
+            pos[0] = atof(inArgs->priVars[index].vals[0].c_str());
+            pos[1] = atof(inArgs->priVars[index].vals[1].c_str());
+            pos[2] = atof(inArgs->priVars[index].vals[2].c_str());
+        }
+    }
+    printf("The position of dislocation (pos): {%f, %f, %f}\n", pos[0], pos[1], pos[2]);
+
+    if((index = GetValID(inArgs->priVars, lineName)) < inArgs->priVars.size()){
+        if(inArgs->priVars[index].vals.size() == 3){
+            line[0] = atof(inArgs->priVars[index].vals[0].c_str());
+            line[1] = atof(inArgs->priVars[index].vals[1].c_str());
+            line[2] = atof(inArgs->priVars[index].vals[2].c_str());
+        }
+    }
+    printf("The line direction of dislocation (line): {%f, %f, %f}\n", line[0], line[1], line[2]);
+
+    if((index = GetValID(inArgs->priVars, gDirName)) < inArgs->priVars.size()){
+        if(inArgs->priVars[index].vals.size() == 3){
+            gDir[0] = atof(inArgs->priVars[index].vals[0].c_str());
+            gDir[1] = atof(inArgs->priVars[index].vals[1].c_str());
+            gDir[2] = atof(inArgs->priVars[index].vals[2].c_str());
+        }
+    }
+    printf("The glide direction of dislocation (gDir): {%f, %f, %f}\n", gDir[0], gDir[1], gDir[2]);
+
+    if((index = GetValID(inArgs->priVars, burgMagName)) < inArgs->priVars.size()){
+         burgMag = atof(inArgs->priVars[index].vals[0].c_str());
+    }
+    printf("The magnitude of dislocation (burgMag): %f\n", burgMag);
+
+    if((index = GetValID(inArgs->priVars, screwTypeName)) < inArgs->priVars.size()){
+         screwType = atoi(inArgs->priVars[index].vals[0].c_str());
+    }
+    printf("The type of screw dislocation (stype): %d\n", screwType);
+
+    if((index = GetValID(inArgs->priVars,  multiDisName)) < inArgs->priVars.size()){
+        if(inArgs->priVars[index].vals.size() == 3){
+            nPairs = atoi(inArgs->priVars[index].vals[0].c_str());
+            rsize = atof(inArgs->priVars[index].vals[1].c_str());
+            delta = atof(inArgs->priVars[index].vals[2].c_str()); 
+        }
+    }
+    printf("Multi-pairs (multi): %d pairs with a mesh size %f\n", nPairs, rsize);
+
+    NormalizeVec(gDir);
+    NormalizeVec(line);
+    cross(line, gDir, normal);
+    NormalizeVec(normal);
+    FormatVector(normal, "normal");
+
+    if(fabs(DotProduct(line, gDir)) > 1.0E-5)Fatal("the glide direction should be prependicular to the line direction");
+
+    if(inArgs->help)return;
+
+    if(strstr(inArgs->inpFiles[0].c_str(), ".lmp") != NULL){
+        ReadLMPFile(inArgs->inpFiles[0], dum);
+    }else{
+        ReadDumpFile(inArgs->inpFiles[0], dum);
+    }
+
+    for(i=0; i<3; i++){
+        boundMin[i] = dum.box[i][0];
+        boundMax[i] = dum.box[i][1];
+    } 
+
+    real8  point[3], dvec[3], dis, vec[3], a, b, shift;
+
+    if(nPairs > 0){
+    while(nPairs > 0){
+
+#pragma omp parallel for private(point, dvec, dis, vec, a, b, shift) shared(dum)
+        for(i=0; i<dum.atom.size(); i++){
+            point[0] = dum.atom[i].x;
+            point[1] = dum.atom[i].y;
+            point[2] = dum.atom[i].z;
+//            printf("%f %f %f %d\n", dum.atom[i].x, dum.atom[i].y, dum.atom[i].z, dum.atom[i].id);
+        
+            if(POINT_INTERSECTION == PointLineIntersection(point, line, pos, dvec, &dis)){
+//                printf("%d: %f %f %f - %f %f %f  @@ %f\n", i, pos[0], pos[1], pos[2], point[0], point[1], point[2], dis);
+                continue;
+            }
+        
+            vec[0] = dvec[0] - point[0];
+            vec[1] = dvec[1] - point[1];
+            vec[2] = dvec[2] - point[2];
+        
+            a = DotProduct(vec, gDir);
+            b = DotProduct(vec, normal); 
+        
+            if(screwType == 0){
+                shift = (0.5*burgMag*atan2(b, a)/M_PI);
+            }else{
+                shift = (0.5*burgMag*(M_PI + atan2(b, a)/M_PI));
+            }
+        
+            dum.atom[i].x += (shift*line[0]);
+            dum.atom[i].y += (shift*line[1]);
+            dum.atom[i].z += (shift*line[2]);
+        
+            FoldBox(boundMin, boundMax, &(dum.atom[i].x), &(dum.atom[i].y), &(dum.atom[i].z));
+        }
+
+        printf("pair %d: burgMag %f, position {%f,%f,%f}, rsize %f\n", nPairs, burgMag, pos[0], pos[1], pos[2], rsize);
+        pos[0] += (rsize*gDir[0]);
+        pos[1] += (rsize*gDir[1]);
+        pos[2] += (rsize*gDir[2]);
+
+        burgMag = -burgMag;
+
+        rsize = delta;
+        nPairs--;
+    }    
+    }else{
+        for(i=0; i<dum.atom.size(); i++){
+            FoldBox(boundMin, boundMax, &(dum.atom[i].x), &(dum.atom[i].y), &(dum.atom[i].z));
+        }
+    }
+
+//    WriteDumpFile(inArgs->outFiles[0], dum); 
+    MGToLMPDataFile(inArgs->outFiles[0], dum);
+    return;
+}
+
+
+#define DisDisPlacement(r) (0.5 + sqrt(M_PI)*(r)*(15.0+10.0*M_PI*(r)*(r) + 2.0*M_PI*M_PI*(r)*(r)*(r)*(r)) \
+                              / (4.0*pow((2.0 + M_PI*(r)*(r)), 2.5)) )
+
+void GenerateExtendedDislocation(InArgs_t *inArgs)
+{
+    int     index, nDiss, i;
+    string  posName("pos"), lineName("line"), separationName("separation"), gDirName("gdir");
+    string  burgName("burg");
+    real8   boundMax[3], boundMin[3], pos[3];
+    real8   line[3] = {1,0,0}, separation = 10.0, gDir[3] = {0,1,0}, normal[3];
+    vector<vector<real8> > burgList, posList;
+
+    if((index = GetValID(inArgs->priVars, posName)) < inArgs->priVars.size()){
+        if(inArgs->priVars[index].vals.size() == 3){
+            pos[0] = atof(inArgs->priVars[index].vals[0].c_str());
+            pos[1] = atof(inArgs->priVars[index].vals[1].c_str());
+            pos[2] = atof(inArgs->priVars[index].vals[2].c_str());
+        }
+    }
+    printf("The drag dislocation position (pos): {%f, %f, %f}\n", pos[0], pos[1], pos[2]);
+
+    if((index = GetValID(inArgs->priVars, separationName)) < inArgs->priVars.size()){
+        separation = atof(inArgs->priVars[index].vals[0].c_str());
+    }
+    printf("The separation of dislocations (separation): %f\n", separation);
+
+    if((index = GetValID(inArgs->priVars, lineName)) < inArgs->priVars.size()){
+        if(inArgs->priVars[index].vals.size() == 3){
+            line[0] = atof(inArgs->priVars[index].vals[0].c_str());
+            line[1] = atof(inArgs->priVars[index].vals[1].c_str());
+            line[2] = atof(inArgs->priVars[index].vals[2].c_str());
+        }
+        NormalizeVec(line);
+    }
+    printf("The line direction of dislocation (line): {%f, %f, %f}\n", line[0], line[1], line[2]);
+
+    if((index = GetValID(inArgs->priVars, gDirName)) < inArgs->priVars.size()){
+        if(inArgs->priVars[index].vals.size() == 3){
+            gDir[0] = atof(inArgs->priVars[index].vals[0].c_str());
+            gDir[1] = atof(inArgs->priVars[index].vals[1].c_str());
+            gDir[2] = atof(inArgs->priVars[index].vals[2].c_str());
+        }
+        NormalizeVec(gDir);
+    }
+    printf("The glide direction of dislocation (gDir): {%f, %f, %f}\n", gDir[0], gDir[1], gDir[2]);
+
+    if((index = GetValID(inArgs->priVars, burgName)) < inArgs->priVars.size()){
+
+        nDiss = (int)(inArgs->priVars[index].vals.size());
+        nDiss /= 3;
+        burgList.resize(nDiss);
+        posList.resize(nDiss);
+
+        printf("%d dislocations: \n", nDiss);
+        for(i=0; i<burgList.size(); i++){
+            burgList[i].resize(3);
+            posList[i].resize(3);
+            burgList[i][0] = atof(inArgs->priVars[index].vals[3*i].c_str());
+            burgList[i][1] = atof(inArgs->priVars[index].vals[3*i+1].c_str());
+            burgList[i][2] = atof(inArgs->priVars[index].vals[3*i+2].c_str());
+
+            posList[i][0] = pos[0] + 2*gDir[0]*separation*((real8)i);
+            posList[i][1] = pos[1] + 2*gDir[1]*separation*((real8)i);
+            posList[i][2] = pos[2] + 2*gDir[2]*separation*((real8)i);
+            
+            printf("%d: position {%f,%f,%f}, burg {%f,%f,%f}\n", i, posList[i][0], posList[i][1], posList[i][2],
+                    burgList[i][0], burgList[i][1], burgList[i][2]);
+        }
+    }else{
+        Fatal("use -dburg to define dislocation lists");
+    }
+
+    if(inArgs->help)return;
+    printf("Reading atoms file %s ...\n", inArgs->inpFiles[0].c_str());
+
+    Dump_t    dum;
+    if(strstr(inArgs->inpFiles[0].c_str(), ".lmp") != NULL){
+        ReadLMPFile(inArgs->inpFiles[0], dum);
+    }else{
+        ReadDumpFile(inArgs->inpFiles[0], dum);
+    }
+
+    for(i=0; i<3; i++){
+        boundMin[i] = dum.box[i][0];
+        boundMax[i] = dum.box[i][1];
+    } 
+    printf("%d atoms have been read\n", (int)dum.atom.size());
+
+    if(fabs(DotProduct(line, gDir)) > 1.0E-2){
+        Fatal("glide direction should be perpendicular to line direction.");
+    }
+
+    cross(line, gDir, normal);
+    NormalizeVec(normal);
+    printf("The glide plane is {%f,%f,%f}\n", normal[0], normal[1], normal[2]);
+
+    int     j;
+    real8   point[3], vec1[3], vec2[3], displace, vec[3], dis1, dis2;
+    real8   a, b, angle;
+
+
+#pragma omp parallel for shared(dum) private (j,point,vec,vec1,vec2,dis1,dis2,displace,a,b, angle)
+    for(i=0; i<dum.atom.size(); i++){
+
+        point[0] = dum.atom[i].x;
+        point[1] = dum.atom[i].y;
+        point[2] = dum.atom[i].z;
+        
+        for(j=0; j<nDiss; j++){
+            VECTOR_COPY(vec, posList[j]);
+
+            PointLineIntersection(point, line, vec, vec1, &dis1);
+            vec1[0] -= point[0];
+            vec1[1] -= point[1];
+            vec1[2] -= point[2];
+
+            NormalizeVec(vec1);
+
+            a = DotProduct(vec1, gDir);
+            b = DotProduct(vec1, normal); 
+
+            angle = atan2(b, a)/M_PI;
+            if(angle < 0)angle += 1;
+            angle /= 2.0;
+        
+            displace = angle*DisDisPlacement(dis1);
+            if(std::isnan(displace)){
+                displace = 0;
+                printf("0");
+            }
+            dum.atom[i].x += burgList[j][0]*displace;
+            dum.atom[i].y += burgList[j][1]*displace;
+            dum.atom[i].z += burgList[j][2]*displace;
+            
+        }
+//        FoldBox(boundMin, boundMax, &(dum.atom[i].x), &(dum.atom[i].y), &(dum.atom[i].z));
+    }
+    
+    MGToLMPDataFile(inArgs->outFiles[0], dum);
+
+    if(inArgs->outFiles.size() == 2){
+        WriteDumpFile(inArgs->outFiles[1], dum); 
+    }
+    return;
+}
+
+void GenerateDislocation(InArgs_t *inArgs){
+    int     index;
+    string  disTypeName("type");
+    int     type = 0;
+    
+    if((index = GetValID(inArgs->priVars, disTypeName)) < inArgs->priVars.size()){
+        if(inArgs->priVars[index].vals.size() == 1){
+            type = atoi(inArgs->priVars[index].vals[0].c_str());
+        }
+    }
+
+    switch(type){
+        case 0:
+            GenerateScrewDislocation(inArgs);
+            break;
+
+        case 1:
+            GenerateExtendedDislocation(inArgs);
+            break;
+
+        default:
+            GenerateExtendedDislocation(inArgs);
+            break;
+    }
+
+    return;
+}
diff --git a/src/HandleExtendedDislocation.cpp b/src/HandleExtendedDislocation.cpp
index 1814f7e936ce9a1b85278d6383c3832f6d230f94..220a8f2885bc55a5ee49c8939be016225ee9ec7a 100644
--- a/src/HandleExtendedDislocation.cpp
+++ b/src/HandleExtendedDislocation.cpp
@@ -127,8 +127,13 @@ void HandleExtendedDislocation_DDD(InArgs_t *inArgs)
         vector<string>().swap(table.variables);
 
         if(!ReadTecplotNormalData(inArgs->inpFiles[file], table, secLine))continue;
+
         vector<string>().swap(words);
         words = split(secLine, "\"");
+
+        if(words.size() == 0){
+            Fatal("can not find the timesolution in file %s", inArgs->inpFiles[file].c_str());
+        }
  
         if(table.data.size() < 2)Fatal("the data size is %d", table.data.size());
  
@@ -385,15 +390,17 @@ typedef struct {
 }Probe_t;
 
 typedef struct {
-    int timestep;
+    long int timestep;
     double x,y,z;
+    double absy;
+    double vy;
     double separation;
     vector<double> vars;
 }EDState_t;
 
 bool com(Atom_t p, Atom_t q){
-    if(fabs(p.y - q.y) < 0.1){
-        if(fabs(p.z - q.z)< 0.1){
+    if(fabs(p.y - q.y) < 1.0E-20){
+        if(fabs(p.z - q.z)< 1.0E-20){
             return(p.x<q.x);
         }else{
             return(p.z < q.z);
@@ -405,30 +412,35 @@ bool com(Atom_t p, Atom_t q){
 
 bool com2(Probe_t p, Probe_t q)
 {
-    return(p.y<p.y);
+    return(p.y<q.y);
 }
 bool com3(EDState_t p, EDState_t q)
 {
     return(p.timestep<q.timestep);
 }
+bool com4(Probe_t p, Probe_t q)
+{
+    return(p.nbr.size() > q.nbr.size());
+}
 
 void HandleExtendedDislocation_MD(InArgs_t *inArgs)
 {
-    int         i, j, file, index, lastIndex, firstNbr, indexVar;
+    int         i, j, file, index, firstNbr, indexVar;
     int         nums[2] = {18, 24}, logfile;
     int         stepID;
-    double      cutofflen = 2.556, alpha = 0.005, beta = 1.0;
-    double      position[3], separation, dis, maxDis;
-    MgData_t    mg;
+    double      cutofflen = 2.556, alpha = 0.1, beta = 1.0;
+    double      position[3];
+    Dump_t    dum;
     LineList_t  list;
-    double      dir[3] = {0, 1, 0}, p0[3], dval=5; 
+    double      dir[3] = {0, 1, 0}, p0[3], dval=5, effeSepRange[2] = {1.0,40.0}; 
     bool        noP0 = 1;
     string      paraName("para"), dvarName("dvar"), dvar("c_vcna"), p0Name("p0"), dirName("dir");
-    string      numsName("nums");
+    string      numsName("nums"), ppstring("pp"), separationName("separation");
+    int     lastIndex = 0;
+    double  dis, minDis;
 
     vector<Atom_t>::iterator     it; 
     vector<EDState_t>   states;
-    EDState_t   state;
 
     ofstream out;
     out.open(inArgs->outFiles[0].c_str(), ios::out);
@@ -477,105 +489,112 @@ void HandleExtendedDislocation_MD(InArgs_t *inArgs)
             nums[1] = atoi(inArgs->priVars[index].vals[1].c_str());
         }
     }
-    printf("The range of effective atoms of a paritial (nums) is [%d,%d]\n", nums[0], nums[1]);
+    printf("The range of effective number of a paritial (nums) is [%d,%d]\n", nums[0], nums[1]);
+
+    if((index = GetValID(inArgs->priVars, separationName)) < inArgs->priVars.size()){
+        if(inArgs->priVars[index].vals.size() == 2){
+            effeSepRange[0] = atof(inArgs->priVars[index].vals[0].c_str());
+            effeSepRange[1] = atof(inArgs->priVars[index].vals[1].c_str());
+        }
+    }
+    printf("The range of effective separtion of the extended dislocation (separation) is [%f,%f]\n", effeSepRange[0], effeSepRange[1]);
 
     if(inArgs->help)return;
     logfile = ReadDataFromMDLogFile(inArgs->auxFiles, list);
 
+    states.resize(inArgs->inpFiles.size());
     for(file=0; file<inArgs->inpFiles.size(); file++){
-        ReadMGDataFile(inArgs->inpFiles[file], mg);
+        states[file].timestep = 1E10;
+    }
 
-        for(i=0; i<mg.variables.size(); i++){
-            if(dvar == mg.variables[i])break;
+#pragma omp parallel for private(dum, it, i, j, indexVar, firstNbr, lastIndex, dis, minDis) shared(states) 
+    for(file=0; file<inArgs->inpFiles.size(); file++){
+        #pragma omp critical
+        {
+            ReadDumpFile(inArgs->inpFiles[file], dum);
+        }
+
+        for(i=0; i<dum.variables.size(); i++){
+            if(dvar == dum.variables[i])break;
         }
-        if((i=indexVar) == mg.variables.size()){
-            Fatal("There is no %s in the mg file %s", dvar.c_str(), 
+        if((indexVar = i) == dum.variables.size()){
+            Fatal("There is no %s in the dum file %s", dvar.c_str(), 
                     inArgs->inpFiles[file].c_str());
         }
         
-        for(it = mg.atom.end(), i=mg.atom.size()-1; it != mg.atom.begin(); it--, i--){
-            if((*it).x < mg.box[0] + (mg.box[1]-mg.box[0])*0.05 ||
-               (*it).y < mg.box[2] + (mg.box[3]-mg.box[2])*0.05 ||
-               (*it).z < mg.box[4] + (mg.box[5]-mg.box[4])*0.05 ||
-               (*it).x > mg.box[0] + (mg.box[1]-mg.box[0])*0.95 ||
-               (*it).y > mg.box[2] + (mg.box[3]-mg.box[2])*0.95 ||
-               (*it).z > mg.box[4] + (mg.box[5]-mg.box[4])*0.95){
-        
-                if(mg.atom[i].vars[indexVar] == dval){
-                    mg.atom.erase(it); 
-                }
+        for(i=0; i<dum.atom.size(); ){
+            if(dum.atom[i].x < dum.box[0][0] + 5 ||
+               dum.atom[i].y < dum.box[1][0] + 5 ||
+               dum.atom[i].z < dum.box[2][0] + 5 ||
+               dum.atom[i].x > dum.box[0][1] - 5 ||
+               dum.atom[i].y > dum.box[1][1] - 5 ||
+               dum.atom[i].z > dum.box[2][1] - 5){
+                dum.atom.erase(dum.atom.begin()+i);
+                continue;
             }
+            if(dum.atom[i].vars[indexVar] != dval){
+                dum.atom.erase(dum.atom.begin()+i);
+                continue;
+            }
+
+
+            i++;
         }
 
+        if(dum.atom.size() == 0){
+            continue;
+        }
+        
         if(noP0 ==1 && file == 0){
-            j = 0;
-            for(i=0; i<mg.atom.size(); i++){
-                 if(mg.atom[i].vars[indexVar] != dval)continue;
-                 j++;
-                p0[0] += mg.atom[i].x;
-                p0[2] += mg.atom[i].z;
-            }
-            if(j==0){
-                Fatal("The");
-            }
-            p0[1]  = mg.box[2];
-            p0[0] /= ((double)j);
-            p0[2] /= ((double)j);
-            printf("The initial probe point (p0) is %f %f %f\n", p0[0], p0[1], p0[2]);
+            Fatal("The initial probe point (p0) is %f %f %f\n", p0[0], p0[1], p0[2]);
         }
 
-        sort(mg.atom.begin(), mg.atom.end(), com);
+        sort(dum.atom.begin(), dum.atom.end(), com);
 
         vector<Probe_t> probes;
         Probe_t         probe;
-        double          d = 0.0;
 
-        probe.y = p0[1];
+        probe.y = dum.atom[0].y;
         probe.z = p0[2];
 
-        lastIndex = 0;
+        lastIndex = 0;  
         while(1){
             vector<Atom_t *>().swap(probe.nbr);
 
-            probe.x += (dir[0]*d*cutofflen*alpha);
-            probe.y += (dir[1]*d*cutofflen*alpha);
-            probe.z += (dir[2]*d*cutofflen*alpha);
-            
             firstNbr = 1;
-            maxDis = -1E10;
-            for(i=lastIndex; i<mg.atom.size(); i++){
-                dis = sqrt(pow((mg.atom[i].y-probe.y), 2) +
-                      pow((mg.atom[i].z-probe.z), 2) );
-                if(dis < cutofflen && mg.atom[i].vars[indexVar] == dval){
+            minDis = 1E10;
+            for(i=lastIndex; i<dum.atom.size(); i++){
+                dis = sqrt(pow((dum.atom[i].y-probe.y), 2) +
+                      pow((dum.atom[i].z-probe.z), 2) );
+                if(dis < minDis){
+                    minDis = dis;
+                }
+
+                if(dis < cutofflen){
                     if(firstNbr){
                         lastIndex = i;
                         firstNbr = 0;
                     }
-                    probe.nbr.push_back(&mg.atom[i]);
+                    probe.nbr.push_back(&(dum.atom[i]));
                 }
             }
             if(probe.nbr.size() > nums[0] && probe.nbr.size()<nums[1]){
                 probes.push_back(probe);
                 vector<Atom_t *>().swap(probe.nbr);
             }
-            if(probe.y > mg.atom.back().y ||
-               probe.z > mg.atom.back().z)break;
-            d++;
+            if(probe.y > dum.atom.back().y)break;
+
+            if(minDis < alpha){
+                probe.x += (dir[0]*alpha);
+                probe.y += (dir[1]*alpha);
+                probe.z += (dir[2]*alpha);
+            }else{
+                probe.x += (dir[0]*minDis);
+                probe.y += (dir[1]*minDis);
+                probe.z += (dir[2]*minDis);
+            }
         }
 
-#if 0
-        printf("Timestep %d, Atoms %d, Bouds %s %s %s, ", 
-               mg.timestep, (int)mg.atom.size(), mg.bounds[0].c_str(), 
-               mg.bounds[1].c_str(), mg.bounds[2].c_str());
-//        printf("Box %e %e %e %e %e %e\n", mg.box[0], mg.box[1],
-//               mg.box[2], mg.box[3], mg.box[4], mg.box[5]);
-
-        printf("variables are ");
-        for(i=0; i<mg.variables.size(); i++){
-            printf("%s(%d) ", mg.variables[i].c_str(), (int)mg.atom.size());
-        }
-        printf("\n");
-#endif
         for(i=0; i<probes.size(); i++){
             probes[i].x = 0.0;
             probes[i].y = 0.0;
@@ -590,40 +609,57 @@ void HandleExtendedDislocation_MD(InArgs_t *inArgs)
             probes[i].z /= (double)probes[i].nbr.size();
 
         }
-        sort(probes.begin(), probes.end(), com2);
 
-        separation = 0;
+        double separation = 0;
 
         if(probes.size() > 0){
             separation = fabs(probes[0].y-probes.back().y);
-            if(separation > 0.1 && separation < 20){
-                state.x = (probes[0].x+probes.back().x)*0.5;
-                state.y = (probes[0].y+probes.back().y)*0.5;
-                state.z = (probes[0].z+probes.back().z)*0.5;
-                state.timestep = mg.timestep;
-                state.separation = separation;
+            if(effeSepRange[0] > 0 && effeSepRange[1] >0 ){
+                 sort(probes.begin(), probes.end(), com2);
+                if(separation >  effeSepRange[0] && separation <  effeSepRange[1] ){
+                    states[file].x = (probes[0].x+probes.back().x)*0.5;
+                    states[file].y = (probes[0].y+probes.back().y)*0.5;
+                    states[file].z = (probes[0].z+probes.back().z)*0.5;
+                    states[file].timestep = dum.timestep;
+                    states[file].separation = separation;
+
+                    if(logfile ==1){
+                        states[file].vars.resize(list.variables.size()-1);
+                        for(i=0; i<list.data[0].size(); i++){
+                            if(states[file].timestep == ((int)list.data[0][i]))break;
+                        }
+                        for(j=0; j<states[file].vars.size(); j++){
+                            states[file].vars[j] = list.data[j+1][i];
+                        }
+                    }
+                }
+            }else{
+                sort(probes.begin(), probes.end(), com4);
+                
+                states[file].x = probes[0].x;
+                states[file].y = probes[0].y;
+                states[file].z = probes[0].z;
+                states[file].timestep = dum.timestep;
+                states[file].separation = 0;
+                printf("file %d, nbrsize %d, timestep %d %f\n",file, (int)probes[0].nbr.size(), (int)states[file].timestep, states[file].y);
                 if(logfile ==1){
-                    state.vars.resize(list.variables.size()-1);
+                    states[file].vars.resize(list.variables.size()-1);
                     for(i=0; i<list.data[0].size(); i++){
-                        if(state.timestep == ((int)list.data[0][i]))break;
+                        if(states[file].timestep == ((int)list.data[0][i]))break;
                     }
-                    for(j=0; j<state.vars.size(); j++){
-                        state.vars[j] = list.data[j+1][i];
+                    for(j=0; j<states[file].vars.size(); j++){
+                        states[file].vars[j] = list.data[j+1][i];
                     }
                 }
-                states.push_back(state);
+
             }
         }
         
     }
 
-    if(states.size() == 0){
-        Fatal("No out data can be dumped, please check parameters");
-    }
-
     sort(states.begin(), states.end(), com3);
 
-    out << "variables = timestep, separation, y";
+    out << "variables = timestep, separation, y, absy, vy";
     if(logfile == 1){
         for(i=1; i<list.variables.size(); i++){
             out << ", " << list.variables[i];
@@ -631,10 +667,28 @@ void HandleExtendedDislocation_MD(InArgs_t *inArgs)
     }
     out << endl;
 
+    real8    lastAbsy, vy, lasty, lastTimestep;
+    ReadDumpFile(inArgs->inpFiles[0], dum);
+    real8    ybox = dum.box[1][1];
+    printf("effective states %d, boxy is %f\n", (int)states.size(), dum.box[1][1]);
     for(i=0; i<states.size(); i++){
+        if(states[i].timestep > 1E9)break;
+        
+        lastAbsy = ((i==0) ? states[i].y : states[i-1].absy);
+        lasty = ((i==0) ? states[i].y : states[i-1].y);
+        lastTimestep = ((i==0) ? states[i].timestep : states[i-1].timestep);
         out << states[i].timestep << " ";
         out << states[i].separation << " ";
         out << states[i].y  << " ";
+
+        vy = states[i].y - lasty;
+        vy -= rint(vy * 1.0/ybox) * ybox;
+        states[i].absy = lastAbsy+vy;
+        out << states[i].absy  << " ";
+
+        states[i].vy = ((i==0) ? 0 : (vy)/(states[i].timestep - lastTimestep));
+        out << states[i].vy  << " ";
+
         for(j=0; j< states[i].vars.size(); j++){
             out << setprecision(10) << states[i].vars[j] << " "; 
         }
diff --git a/src/Main.cpp b/src/Main.cpp
index 1d656fd01ad02bf257621bb644e4f5d2bf3a6943..1678f0900cdf057ab81030bd826bce698305b516 100644
--- a/src/Main.cpp
+++ b/src/Main.cpp
@@ -50,7 +50,8 @@ Option_t optList[OPT_MAX] = {
         {OPT_SEED,          "seed",     	3, 1},
         {OPT_TYPE,          "type",     	1, 1},
         {OPT_PRIVATEVALS,   "d",            1, 1},
-        {OPT_AUXFILE,       "auxfile",      3, 1}     
+        {OPT_AUXFILE,       "auxfile",      3, 1},     
+        {OPT_THREADS,       "nthreads",     1, 1}     
 };
 
 
@@ -104,7 +105,11 @@ static void InitDefaultValues(InArgs_t *inArgs)
         inArgs->seed       = time(0) + getpid();
         inArgs->type       = 0;
         inArgs->help       = 0;
-        	
+#ifdef _OPENMP
+        inArgs->nThreads   = 4;
+#else
+        inArgs->nThreads   = 1;
+#endif        	
         return;
 }
 
@@ -252,6 +257,9 @@ static void GetInArgs(int argc, char *argv[], InArgs_t *inArgs)
                 case OPT_AUXFILE:
                     swap(inArgs->auxFiles, argValues);
                     break;
+                case OPT_THREADS:
+                    inArgs->nThreads = atoi(argValue.c_str());
+                    break;
                 default:
                     break;
             }
@@ -294,15 +302,15 @@ static void GetInArgs(int argc, char *argv[], InArgs_t *inArgs)
             swap(bakInps, inArgs->auxFiles);
         }
         sort(inArgs->auxFiles.begin(), inArgs->auxFiles.end());
-#if 0
-        if(inArgs->inpFiles.size()>1){
+#if 1
+        if(inArgs->inpFiles.size()<5){
             printf("Input files are:\n");
             for(i=0; i<inArgs->inpFiles.size(); i++){
                 printf("%s \n", inArgs->inpFiles[i].c_str());
             }
         }
     
-        if(inArgs->auxFiles.size() > 0){
+        if(inArgs->auxFiles.size() < 5){
             printf("Auxiliary file(s) :\n");
             for(i=0; i<inArgs->auxFiles.size(); i++){
                 printf("%s \n", inArgs->auxFiles[i].c_str());
@@ -324,6 +332,11 @@ int main(int argc, char *argv[])
         CheckArgSanity(&inArgs);
         PrintArgs(&inArgs); 
 
+#ifdef _OPENMP
+        omp_set_num_threads(inArgs.nThreads); 
+
+        printf("Threads count: %d \n", omp_get_max_threads());
+#endif
 /*
  *      All that's left is to invoke the proper function to
  *      handle data 
@@ -340,6 +353,10 @@ int main(int argc, char *argv[])
                 break;
             case FTYPE_SPECIFY_EQUATIONS:
                 SpecifyEquations_PLTDATA(&inArgs);
+                break;
+            case FTYPE_GENERATE_DISLOCATION:
+                GenerateDislocation(&inArgs);
+                break;
         }
 
         exit(0);
diff --git a/src/Math.cpp b/src/Math.cpp
index 9ff5679ecf1fe51cfd4b62d6fcdb20b6b274fd74..0693623c2009ed922a7dadd3426835d25bfe8a4c 100644
--- a/src/Math.cpp
+++ b/src/Math.cpp
@@ -12,34 +12,50 @@ void AverageLines(InArgs_t *inArgs)
 {
     int     index, i, j, k, colX = -1, colY;
     int     readState;
-    bool    firstFile = 1, calTau = 0, specifyEqu = 0;
-    real8   rsize = 20, min, max, effNums = 0, value;
+    bool    firstFile = 1, specifyEqu = 0;
+    real8   rsize = 0, min, max, effNums = 0, value;
     string  rsizeName("rsize"), varsName("vars"), overName("over"), specifyEquName("spe");
-    string  str("Ave_"), secLine, tauName("tau"), overVar;
+    string  str("Ave_"), secLine, tauName("tau"), overVar, weighName("weigh"), weightCoeff;
 
+    bool        calTau = 0;
     LineList_t  list;
-    Curve_t     curve;
+    InitList(list);
 
+    Curve_t     curve;
     vector<int>         varID;
     vector<real8>       seq;
     vector<Table_t>     tables(inArgs->inpFiles.size());
     vector<string>      s1, s2;
     vector<vector<vector<real8> > > array;
 
+    bool weigh = false;
+    if((index = GetValID(inArgs->priVars, weighName)) < inArgs->priVars.size()){
+        weigh = true;
+        weightCoeff = inArgs->priVars[index].vals[0];
+    }
+    if(weigh){
+        printf("Weight coffeicient (weigh) %s is defined\n", weightCoeff.c_str());
+    }else{
+        printf("No weight coffeicient (weigh) is defied\n");
+    }
+    
     if((index = GetValID(inArgs->priVars, rsizeName)) < inArgs->priVars.size()){
         rsize = atof(inArgs->priVars[index].vals[0].c_str());
     }
-    printf("The remesh size (rsize) is %f\n", rsize);
-
+    if(rsize ==  0){
+        printf("Not change the original remesh size (rsize)\n");
+    }else{
+        printf("The remesh size (rsize) is %f\n", rsize);
+    }
     if((index = GetValID(inArgs->priVars, tauName)) < inArgs->priVars.size()){
-        calTau = atoi(inArgs->priVars[index].vals[0].c_str());
+        calTau = true;
     }
     if(calTau){
         printf("The variance (tau) will been caculated\n");
     }else{
         printf("The variance (tau) will not been caculated\n");
     }
-
+    
     if((index = GetValID(inArgs->priVars, specifyEquName)) < inArgs->priVars.size()){
         specifyEqu = atoi(inArgs->priVars[index].vals[0].c_str());
     }
@@ -90,10 +106,41 @@ void AverageLines(InArgs_t *inArgs)
         Fatal("There is no input file.");
     }
 
-    printf("Ranges of files:\n ");
+    map<string, string>::iterator   iter;
+    vector<real8>   weightList(inArgs->inpFiles.size());
+    real8   totalWeight = 0.0;       
+    int     firstReadFile = -1;        
+    Table_t auxTable;
+    InitTable(auxTable);
+
     for(i=0; i<inArgs->inpFiles.size(); i++){
         readState = ReadTecplotNormalData(inArgs->inpFiles[i], tables[i], secLine);
-        if(!readState)continue;
+        if(!readState){
+            weightList[i] = 0.0;
+            Fatal("can not read file %s", inArgs->inpFiles[i].c_str());
+            continue;
+        }else{
+            if(firstReadFile < 0){
+                firstReadFile = i;
+            }
+        }
+
+        if(weigh){
+//            for(iter=tables[i].aux.begin(); iter != tables[i].aux.end(); iter++){
+//                printf("%s = %s\n", iter->first.c_str(), iter->second.c_str());
+//            }
+//            Fatal("T");
+    
+            if((iter = tables[i].aux.find(weightCoeff)) != tables[i].aux.end()){
+                weightList[i] = atof(iter->second.c_str());
+            }else{
+                printf("waring: no weight variale %s in file %s\n", weightCoeff.c_str(), inArgs->inpFiles[i].c_str());
+                weightList[i] = 1.0;
+            }
+        }else{
+            weightList[i] = 1.0;
+        }
+        totalWeight += weightList[i];
 
         if(specifyEqu){
             SpecifyEquations(tables[i]);
@@ -138,19 +185,55 @@ void AverageLines(InArgs_t *inArgs)
             min = tables[i].data[0][colX];
             max = tables[i].data[tables[i].data.size()-1][colX];
             firstFile = 0;
+
+            if(!tables[i].aux.empty()){
+                auxTable.data.resize(inArgs->inpFiles.size());
+                auxTable.i = (int)inArgs->inpFiles.size();
+                printf("Auxiliary variables: ");
+
+                auxTable.variables.push_back("file");
+                auxTable.data[0].push_back((real8)i);
+
+                for(const auto &pair : tables[i].aux){
+                    printf("%s ", pair.first.c_str());
+                    auxTable.variables.push_back(pair.first);
+                    auxTable.data[0].push_back(atof(pair.second.c_str()));
+                }
+                printf("\n");
+            }
+
         }else{
             if(min < tables[i].data[0][colX]) min = tables[i].data[0][colX];
             if(max > tables[i].data[tables[i].data.size()-1][colX]){
                 max = tables[i].data[tables[i].data.size()-1][colX];
             }
+
+            if(auxTable.variables.size()>0){
+                auxTable.data[i].resize(auxTable.variables.size());
+                j = 0;
+                auxTable.data[i][j++] = (real8)i;
+
+                for(const auto &pair : tables[i].aux){
+                    auxTable.data[i][j++] = (atof(pair.second.c_str()));
+                }
+            }
         }
 
-        printf("[%f,%f]\n", tables[i].data[0][colX], tables[i].data[tables[i].data.size()-1][colX]);
+        printf("File %s(%d): Range is [%f,%f], %d points\n", 
+                inArgs->inpFiles[i].c_str(), i, tables[i].data[0][colX], tables[i].data[tables[i].data.size()-1][colX],
+                (int)(tables[i].data.size()));
     }
     printf("The effective range of %s is [%f,%f]\n", overVar.c_str(), min, max);
-    printf("%d files was read\n", (int)effNums);
+    printf("%d files have been read, the total weight is %f\n", (int)effNums, totalWeight);
 
-    seq = GenerateSequence(min, max, rsize);
+    if(rsize == 0){
+        seq.resize(tables[firstReadFile].data.size());
+        for(i=0; i<seq.size(); i++){
+            seq[i] = tables[firstReadFile].data[i][varID[colX]];
+        }
+    }else{
+        seq = GenerateSequence(min, max, rsize);
+    }
     if(calTau){
         array.resize((int)(effNums));
         for(i=0; i<array.size(); i++){
@@ -190,13 +273,27 @@ void AverageLines(InArgs_t *inArgs)
             }else{
                 for(k=0; k<curve.ax.size(); k++){
                     curve.ay[k] = tables[i].data[k][varID[j]];
+//                    printf("%d %e %e\n", k, curve.ay[k], tables[i].data[k][varID[j]]);
                 }
 
-                for(k=0; k<seq.size(); k++){
-                    value = LinearInterpolation(curve, seq[k], min, max);
-                    list.data[j][k] += (value/effNums);
-                    if(calTau){
-                        array[i][j-1][k] = value; 
+                if(rsize == 0){
+                    if(seq.size() != curve.ay.size()){
+                        Fatal("the %d table has a differnt size (%d) from the basic one (%d)",
+                               i, (int)seq.size(), (int)curve.ay.size());
+                    }
+                    for(k=0; k<seq.size(); k++){
+                        list.data[j][k] += (curve.ay[k]*weightList[i]/totalWeight);
+                    }
+//                    printf("%e %e %e %e %d\n", list.data[j][k],curve.ay[k],weightList[i],totalWeight, (int)curve.ay.size()); 
+                    if(calTau)array[i][j-1][k] = curve.ay[k];
+
+                }else{
+                    for(k=0; k<seq.size(); k++){
+                        value = LinearInterpolation(curve, seq[k], min, max);
+                        list.data[j][k] += (value*weightList[i]/totalWeight);
+                        if(calTau){
+                            array[i][j-1][k] = value; 
+                        }
                     }
                 }
             }
@@ -220,7 +317,7 @@ void AverageLines(InArgs_t *inArgs)
         for(i=0; i<array.size(); i++){
             for(j=0; j<array[i].size(); j++){
                 for(k=0; k<array[i][j].size(); k++){
-                    list.data[j+oriVals][k] += (pow((array[i][j][k] - list.data[j+1][k]),2)/((real8)array.size()));
+                    list.data[j+oriVals][k] += (pow((array[i][j][k] - list.data[j+1][k]),2)*weightList[i]/totalWeight);
                 }
             }
         }
@@ -232,7 +329,469 @@ void AverageLines(InArgs_t *inArgs)
         }
     }
 
+    if(auxTable.data.size() > 0){
+
+        real8 aveAuxVal;
+        char  val[50];
+        for(i=1; i<auxTable.variables.size(); i++){
+            if(auxTable.variables[i] ==  weightCoeff)continue;
+            
+            aveAuxVal = 0.0;
+            for(j=0; j<auxTable.data.size(); j++){
+                aveAuxVal += (auxTable.data[j][i]*weightList[j]/totalWeight);
+            }
+            
+            snprintf(val, sizeof(val), "%e", aveAuxVal);
+            auxTable.aux[std::string("Ave_")+auxTable.variables[i]] = val;
+            list.aux[std::string("Ave_")+auxTable.variables[i]] = val;
+        }
+
+        if(inArgs->outFiles.size() < 2){
+//            WriteTecplotNormalData(auxTable, std::string("auxTable.plt"), 10);
+        }else{
+            WriteTecplotNormalData(auxTable,  inArgs->outFiles[1], 10);
+        }
+    }
+
+    list.i = ((int)list.data[0].size());
     WriteTecplotNormalData(list, inArgs->outFiles[0], 10);
 
     return;
 }
+
+
+int PlanePlaneIntersection(double *n1, double *p1, double *n2, double*p2,
+                           double *vec1, double *vec2)
+{
+    double  n1x = n1[0], n1y = n1[1], n1z = n1[2];
+    double  n2x = n2[0], n2y = n2[1], n2z = n2[2];
+    double  p1x = p1[0], p1y = p1[1], p1z = p1[2];
+    double  p2x = p2[0], p2y = p2[1], p2z = p2[2];
+    double  dir[3], dir1[3], dir2[3], dot, dis, vec3[3], vec4[4], vec5[3], vec6[3];
+    int     type;
+
+    VECTOR_ZERO(vec1);
+    VECTOR_ZERO(vec2);
+
+    if(sqrt(SQUARE(fabs(-(n1y*n2x) + n1x*n2y)) + 
+       SQUARE(fabs(n1z*n2x - n1x*n2z)) + 
+       SQUARE(fabs(-(n1z*n2y) + n1y*n2z))) < EPS1){
+        if(fabs(n1x*(p1x - p2x) + n1y*(p1y - p2y) + n1z*(p1z - p2z)) < EPS2){
+            VECTOR_COPY(vec1,n2);
+            VECTOR_COPY(vec2,p2);
+            return(PLANE_INTERSECTION);
+        }
+
+        VECTOR_COPY(vec1,n2);
+        VECTOR_COPY(vec2,p2);
+        return(NO_INTERSECTION);
+    } 
+
+    vec1[0] = -(n1z*n2y) + n1y*n2z;
+    vec1[1] = n1z*n2x - n1x*n2z;
+    vec1[2] = -(n1y*n2x) + n1x*n2y;
+    NormalizeVec(vec1);
+
+    if(sqrt(SQUARE(fabs(p1x - p2x)) + SQUARE(fabs(p1y - p2y)) + 
+       SQUARE(fabs(p1z - p2z))) < EPS2){
+        VECTOR_COPY(vec2,p2);
+        return(LINE_INTERSECTION);
+    }
+
+    cross(vec1, n1, vec3);
+    if(Normal(vec3)< 1.0E-4){
+        Fatal("PlanePlaneIntersection in %s at %d\n", __FILE__, __LINE__);
+    }
+    NormalizeVec(vec3);
+
+    type = LinePlaneIntersection(vec3, p1, n2, p2, vec4, vec2);
+    if(type == NO_INTERSECTION){
+        Fatal("PlanePlaneIntersection in %s at %d\n", __FILE__, __LINE__);
+    }
+#if 0
+    dir[0] = p2[0] - p1[0];
+    dir[1] = p2[1] - p1[1];
+    dir[2] = p2[2] - p1[2];
+    if(Normal(dir) < 1.0E-2){
+        VECTOR_COPY(vec2, p2);
+        return(LINE_INTERSECTION);
+    }
+    NormalizeVec(dir);
+
+    dot = DotProduct(dir, n1);
+    dir1[0] = dir[0] - dot*n1[0];
+    dir1[1] = dir[1] - dot*n1[1];
+    dir1[2] = dir[2] - dot*n1[2];
+    if(Normal(dir1) < 1.0E-4){
+        VECTOR_COPY(vec2, p1);
+        return(LINE_INTERSECTION);
+    }
+
+    NormalizeVec(dir1);
+
+    dot = DotProduct(dir, n2);
+    dir2[0] = dir[0] - dot*n2[0];
+    dir2[1] = dir[1] - dot*n2[1];
+    dir2[2] = dir[2] - dot*n2[2];
+    if(Normal(dir2) < 1.0E-4){
+        VECTOR_COPY(vec2, p2);
+        return(LINE_INTERSECTION);
+    }
+
+    NormalizeVec(dir2);
+
+    type = LineLineIntersection(dir1, p1, dir2, p2, vec5, vec2, &dis);
+
+    if(fabs(dis) > 2){
+        FormatVector(vec3, "p1");
+        FormatVector(vec4, "p2");
+        printf("Warning: PlanePlaneIntersection dis is %f\n", dis);
+    }
+#endif
+    return(LINE_INTERSECTION);
+} 
+
+int LinePlaneIntersection(double *d, double *p1, double *n, double *p2, 
+                          double *vec1, double *vec2)
+{
+    double  dx = d[0], dy = d[1], dz = d[2];
+    double  nx = n[0], ny = n[1], nz = n[2];
+    double  p1x = p1[0], p1y = p1[1], p1z = p1[2];
+    double  p2x = p2[0], p2y = p2[1], p2z = p2[2];
+    double  dis;
+
+    VECTOR_ZERO(vec1);
+    VECTOR_ZERO(vec2);
+
+    if(fabs(dx*nx + dy*ny + dz*nz) < EPS1){
+        VECTOR_COPY(vec1, d);
+        PointPlaneIntersection(p1, n, p2, vec2, &dis);
+        if(fabs(dis) < EPS2){
+            VECTOR_COPY(vec2, p1);
+            return(LINE_INTERSECTION);
+        }else{ 
+            return(NO_INTERSECTION);
+        }
+    }
+
+    double  t = (-(nx*p1x) - ny*p1y - nz*p1z + nx*p2x + ny*p2y + nz*p2z)/
+                (dx*nx + dy*ny + dz*nz);
+
+    vec2[0] = p1x + t*dx;    
+    vec2[1] = p1y + t*dy;    
+    vec2[2] = p1z + t*dz;
+
+    if(std::isnan(vec2[0]) || std::isinf(vec2[0]) ||
+       std::isnan(vec2[1]) || std::isinf(vec2[1]) ||
+       std::isnan(vec2[2]) || std::isinf(vec2[2])){
+        printf(" LinePlaneIntersection: %f,%f,%f\n", 
+                vec2[0], vec2[1], vec2[2]);
+    }
+
+    if(std::isnan(vec2[0]))vec2[0] = 0.0;
+    if(std::isnan(vec2[1]))vec2[1] = 0.0;
+    if(std::isnan(vec2[2]))vec2[2] = 0.0;
+    
+    if(std::isinf(vec2[0]) || std::isinf(vec2[1]) || std::isinf(vec2[2]) ){
+        printf("d: {%f,%f,%f}\n", d[0], d[1], d[2]);
+        printf("p1: {%f,%f,%f}\n", p1[0], p1[1], p1[2]);
+        printf("n: {%f,%f,%f}\n", n[0], n[1], n[2]);
+        printf("p2: {%f,%f,%f}\n", p2[0], p2[1], p2[2]);
+        exit(0);
+    }
+#if 0
+    if(sqrt(SQUARE(vec2[2]) + SQUARE(vec2[2]) + SQUARE(vec2[2])) > 10000){
+        printf("vec2: {%f,%f,%f}\n", vec2[0], vec2[1], vec2[2]);
+        printf("d: {%f,%f,%f}\n", d[0], d[1], d[2]);
+        printf("p1: {%f,%f,%f}\n", p1[0], p1[1], p1[2]);
+        printf("n: {%f,%f,%f}\n", n[0], n[1], n[2]);
+        printf("p2: {%f,%f,%f}\n", p2[0], p2[1], p2[2]);
+        printf("PointPlaneIntersection\n");
+    }
+#endif
+
+    return(POINT_INTERSECTION);    
+}
+
+int PointLineIntersection(double *p1, double *d, double *p2,
+                         double *vec, double *dis){
+
+    double  p1x = p1[0], p1y=p1[1], p1z = p1[2];
+    double  p2x = p2[0], p2y=p2[1], p2z = p2[2];
+    double  dx = d[0], dy=d[1], dz = d[2];
+    double  norm, norm2, invNorm, t;
+
+    VECTOR_ZERO(vec);
+    norm = sqrt(SQUARE(dx) + SQUARE(dy) + SQUARE(dz));
+    invNorm = 1.0/norm;
+
+    if(norm < EPS1){
+        Fatal("PointLineIntersection: line direction is zero\n");
+        return(NO_INTERSECTION);
+    }
+
+    norm2 = sqrt(SQUARE(p1x-p2x) + SQUARE(p1y-p2y) + SQUARE(p1z-p2z));
+    if(norm2 < EPS1){
+        vec[0] = p2[0];
+        vec[1] = p2[1];
+        vec[2] = p2[2];
+        *dis = 0.0;
+        printf("0");
+        return(POINT_INTERSECTION);
+    }
+
+    if(fabs(dx*(p1x-p2x) + dy*(p1y-p2y) + dz*(p1z-p2z)) > norm*norm2*(1-EPS1)){
+        vec[0] = p1[0];
+        vec[1] = p1[1];
+        vec[2] = p1[2];
+        *dis = 0.0;
+        printf("0");
+        return(POINT_INTERSECTION);
+    }
+
+    t = (dx*(p1x - p2x) + dy*(p1y - p2y) + dz*(p1z - p2z))*(invNorm*invNorm);
+    vec[0] = p2x + dx*t;
+    vec[1] = p2y + dy*t;
+    vec[2] = p2z + dz*t;
+
+    *dis = sqrt(SQUARE(vec[0]-p1[0]) + SQUARE(vec[1]-p1[1]) + SQUARE(vec[2]-p1[2]) );
+
+    if(std::isnan(*dis) || std::isinf(*dis)){
+        printf(" PointLineIntersection: %f {%f,%f,%f}\n", 
+               *dis, vec[0], vec[1], vec[2]);
+        vec[0] = p1[0];
+        vec[1] = p1[1];
+        vec[2] = p1[2];
+    }
+
+    if(std::isnan(*dis))*dis = 0.0;
+    if(std::isinf(*dis)) *dis = 0.0;
+#if 0
+    if(*dis < EPS2 && sqrt(SQUARE(vec[0]-p1[0]) + SQUARE(vec[0]-p1[0]) + SQUARE(vec[0]-p1[0])) > 10000){
+        printf("vec: {%f,%f,%f}\n", vec[0], vec[1], vec[2]);
+        printf("p1: {%f,%f,%f}\n", p1[0], p1[1], p1[2]);
+        printf("d: {%f,%f,%f}\n", d[0], d[1], d[2]);
+        printf("p2: {%f,%f,%f}\n", p2[0], p2[1], p2[2]);
+        printf("PointLineIntersection");
+    }
+#endif
+    if(*dis < EPS2){
+        return(POINT_INTERSECTION);
+    }else{
+        return(NO_INTERSECTION);
+    }
+}
+
+int LineLineIntersection(double *d1, double *p1, double *d2, double *p2,
+                         double *vec1, double *vec2, double *dis){
+    double  d1x = d1[0], d1y = d1[1], d1z = d1[2];
+    double  d2x = d2[0], d2y = d2[1], d2z = d2[2];
+    double  p1x = p1[0], p1y = p1[1], p1z = p1[2];
+    double  p2x = p2[0], p2y = p2[1], p2z = p2[2];
+
+    VECTOR_ZERO(vec1);
+    VECTOR_ZERO(vec2);
+
+    if(SQUARE(d1x) + SQUARE(d1y) + SQUARE(d1z) < EPS1*EPS1){
+        Fatal("LineLineIntersection: line direction is zero L1\n");
+        return(NO_INTERSECTION);
+    }
+    if(SQUARE(d2x) + SQUARE(d2y) + SQUARE(d2z) < EPS1*EPS1){
+        Fatal("LineLineIntersection: line direction is zero L2\n");
+        return(NO_INTERSECTION);
+    }
+
+    if(fabs(1.0-fabs(d1x*d2x+d1y*d2y+d1z*d2z)) < EPS1){
+        if(PointLineIntersection(p1, d2, p2, vec2, dis) == POINT_INTERSECTION){
+            VECTOR_COPY(vec1, d1);
+            *dis = fabs(*dis);
+            return(LINE_INTERSECTION);
+        }else{
+            *dis = fabs(*dis);
+            VECTOR_COPY(vec1, d2);
+            VECTOR_COPY(vec2, p2);
+            return(NO_INTERSECTION);
+        }
+    }
+
+/*
+ *  return the closed point on the first line
+ */
+    vec1[0] = p1x - (d1x*(-4*(SQUARE(d2x) + SQUARE(d2y) + SQUARE(d2z))*
+              (d1x*(p1x - p2x) + d1y*(p1y - p2y) + d1z*(p1z - p2z)) + 
+              4*(d1x*d2x + d1y*d2y + d1z*d2z)*
+              (d2x*(p1x - p2x) + d2y*(p1y - p2y) + d2z*(p1z - p2z))))/
+              (4*SQUARE(d1x*d2x + d1y*d2y + d1z*d2z) - 
+              4*(SQUARE(d1x) + SQUARE(d1y) + SQUARE(d1z))*
+              (SQUARE(d2x) + SQUARE(d2y) + SQUARE(d2z)));        
+    vec1[1] = p1y - (d1y*(-4*(SQUARE(d2x) + SQUARE(d2y) + SQUARE(d2z))*
+              (d1x*(p1x - p2x) + d1y*(p1y - p2y) + d1z*(p1z - p2z)) + 
+              4*(d1x*d2x + d1y*d2y + d1z*d2z)*
+              (d2x*(p1x - p2x) + d2y*(p1y - p2y) + d2z*(p1z - p2z))))/
+              (4*SQUARE(d1x*d2x + d1y*d2y + d1z*d2z) - 
+              4*(SQUARE(d1x) + SQUARE(d1y) + SQUARE(d1z))*
+              (SQUARE(d2x) + SQUARE(d2y) + SQUARE(d2z)));
+    vec1[2] = p1z - (d1z*(-4*(SQUARE(d2x) + SQUARE(d2y) + SQUARE(d2z))*
+              (d1x*(p1x - p2x) + d1y*(p1y - p2y) + d1z*(p1z - p2z)) + 
+              4*(d1x*d2x + d1y*d2y + d1z*d2z)*
+              (d2x*(p1x - p2x) + d2y*(p1y - p2y) + d2z*(p1z - p2z))))/
+              (4*SQUARE(d1x*d2x + d1y*d2y + d1z*d2z) - 
+              4*(SQUARE(d1x) + SQUARE(d1y) + SQUARE(d1z))*
+              (SQUARE(d2x) + SQUARE(d2y) + SQUARE(d2z))); 
+
+/*
+ *  return the closed point on the second line
+ */
+    vec2[0] = p2x + (d2x*(SQUARE(d1z)*
+              (d2x*p1x + d2y*p1y - d2x*p2x - d2y*p2y) + 
+              SQUARE(d1y)*(d2x*(p1x - p2x) + d2z*(p1z - p2z)) + 
+              d1x*d1z*(-(d2z*p1x) - d2x*p1z + d2z*p2x + d2x*p2z) + 
+              SQUARE(d1x)*(d2y*p1y + d2z*p1z - d2y*p2y - d2z*p2z) + 
+              d1y*(d1x*(-(d2y*p1x) - d2x*p1y + d2y*p2x + d2x*p2y) + 
+              d1z*(-(d2z*p1y) - d2y*p1z + d2z*p2y + d2y*p2z))))/
+              (SQUARE(d1z)*(SQUARE(d2x) + SQUARE(d2y)) - 
+              2*d1x*d1z*d2x*d2z - 2*d1y*d2y*(d1x*d2x + d1z*d2z) + 
+              SQUARE(d1y)*(SQUARE(d2x) + SQUARE(d2z)) + 
+              SQUARE(d1x)*(SQUARE(d2y) + SQUARE(d2z)));
+
+    vec2[1] = p2y + (d2y*(SQUARE(d1z)*
+              (d2x*p1x + d2y*p1y - d2x*p2x - d2y*p2y) + 
+              SQUARE(d1y)*(d2x*(p1x - p2x) + d2z*(p1z - p2z)) + 
+              d1x*d1z*(-(d2z*p1x) - d2x*p1z + d2z*p2x + d2x*p2z) + 
+              SQUARE(d1x)*(d2y*p1y + d2z*p1z - d2y*p2y - d2z*p2z) + 
+              d1y*(d1x*(-(d2y*p1x) - d2x*p1y + d2y*p2x + d2x*p2y) + 
+              d1z*(-(d2z*p1y) - d2y*p1z + d2z*p2y + d2y*p2z))))/
+              (SQUARE(d1z)*(SQUARE(d2x) + SQUARE(d2y)) - 
+              2*d1x*d1z*d2x*d2z - 2*d1y*d2y*(d1x*d2x + d1z*d2z) + 
+              SQUARE(d1y)*(SQUARE(d2x) + SQUARE(d2z)) + 
+              SQUARE(d1x)*(SQUARE(d2y) + SQUARE(d2z)));
+
+    vec2[2] = p2z + (d2z*(SQUARE(d1z)*
+              (d2x*p1x + d2y*p1y - d2x*p2x - d2y*p2y) + 
+              SQUARE(d1y)*(d2x*(p1x - p2x) + d2z*(p1z - p2z)) + 
+              d1x*d1z*(-(d2z*p1x) - d2x*p1z + d2z*p2x + d2x*p2z) + 
+              SQUARE(d1x)*(d2y*p1y + d2z*p1z - d2y*p2y - d2z*p2z) + 
+              d1y*(d1x*(-(d2y*p1x) - d2x*p1y + d2y*p2x + d2x*p2y) + 
+              d1z*(-(d2z*p1y) - d2y*p1z + d2z*p2y + d2y*p2z))))/
+              (SQUARE(d1z)*(SQUARE(d2x) + SQUARE(d2y)) - 
+              2*d1x*d1z*d2x*d2z - 2*d1y*d2y*(d1x*d2x + d1z*d2z) + 
+              SQUARE(d1y)*(SQUARE(d2x) + SQUARE(d2z)) + 
+              SQUARE(d1x)*(SQUARE(d2y) + SQUARE(d2z)));
+
+    *dis = sqrt(SQUARE(d1z*(d2y*(p1x - p2x) + d2x*(-p1y + p2y)) + 
+           d1y*(-(d2z*p1x) + d2x*p1z + d2z*p2x - d2x*p2z) + 
+           d1x*(d2z*p1y - d2y*p1z - d2z*p2y + d2y*p2z))/
+           (SQUARE(d1z)*(SQUARE(d2x) + SQUARE(d2y)) - 
+           2*d1x*d1z*d2x*d2z - 2*d1y*d2y*(d1x*d2x + d1z*d2z) + 
+           SQUARE(d1y)*(SQUARE(d2x) + SQUARE(d2z)) + 
+           SQUARE(d1x)*(SQUARE(d2y) + SQUARE(d2z))));
+
+
+    if(std::isnan(*dis) || std::isinf(*dis)){
+        printf("LineLineIntersection: %f\n", *dis);
+    }
+    if(std::isnan(*dis)) *dis = 0.0;
+    if(std::isinf(*dis)) *dis = 0.0;
+
+    if(fabs(*dis) < EPS2 && sqrt(SQUARE(vec2[0]-p1[0]) + SQUARE(vec2[0]-p1[0]) + SQUARE(vec2[0]-p1[0])) > 10000){
+        printf("vec2: {%f,%f,%f}\n", vec2[0], vec2[1], vec2[2]);
+        printf("d1: {%f,%f,%f}\n", d1[0], d1[1], d1[2]);
+        printf("p1: {%f,%f,%f}\n", p1[0], p1[1], p1[2]);
+        printf("d2: {%f,%f,%f}\n", d2[0], d2[1], d2[2]);
+        printf("p2: {%f,%f,%f}\n", p2[0], p2[1], p2[2]);
+        printf("LineLineIntersection\n");
+    }
+
+    if(*dis < EPS2){
+        return(POINT_INTERSECTION);
+    }else{
+        return(NO_INTERSECTION);
+    }
+}
+
+
+int PointPlaneIntersection(double *p1, double *n, double *p2, double *vec, double *dis)
+{
+    double  p1x = p1[0], p1y = p1[1], p1z = p1[2];
+    double  p2x = p2[0], p2y = p2[1], p2z = p2[2];
+    double  nx = n[0], ny = n[1], nz = n[2];
+    double  temVec[3];
+
+    if(nx*nx + ny*ny + nz*nz < EPS1*EPS1){
+        printf("PointPlaneIntersection: normal is close to zero.\n");
+        return(0);
+    }
+
+    *dis = sqrt(SQUARE(nx*p1x + ny*p1y + nz*p1z - nx*p2x - ny*p2y - nz*p2z)/
+           (SQUARE(nx) + SQUARE(ny) + SQUARE(nz)));
+
+    if(std::isnan(*dis) || std::isinf(*dis)){
+        printf(" PointPlaneIntersection: %f\n", *dis);
+    }
+
+    if(std::isnan(*dis)) *dis = 0.0;
+    if(std::isinf(*dis)) *dis = 0.0;
+
+    vec[0] = (SQUARE(ny)*p1x + SQUARE(nz)*p1x + SQUARE(nx)*p2x + 
+             nx*ny*(-p1y + p2y) + nx*nz*(-p1z + p2z))/
+             (SQUARE(nx) + SQUARE(ny) + SQUARE(nz));
+    vec[1] = (SQUARE(nx)*p1y + SQUARE(nz)*p1y + nx*ny*(-p1x + p2x) + 
+             SQUARE(ny)*p2y + ny*nz*(-p1z + p2z))/
+             (SQUARE(nx) + SQUARE(ny) + SQUARE(nz));
+    vec[2] = (SQUARE(nx)*p1z + SQUARE(ny)*p1z + nx*nz*(-p1x + p2x) + 
+             ny*nz*(-p1y + p2y) + SQUARE(nz)*p2z)/
+             (SQUARE(nx) + SQUARE(ny) + SQUARE(nz));
+
+    temVec[0] = p1[0] - vec[0];
+    temVec[1] = p1[1] - vec[1];
+    temVec[2] = p1[2] - vec[2];
+    NormalizeVec(temVec);
+    if(DotProduct(temVec, n) < 0){
+        *dis *= -1.0;
+    }
+    
+    if(fabs(*dis) < EPS2 && sqrt(SQUARE(vec[0]-p1[0]) + SQUARE(vec[0]-p1[0]) + SQUARE(vec[0]-p1[0])) > 10000){
+        printf("vec: {%f,%f,%f}\n", vec[0], vec[1], vec[2]);
+        printf("p1: {%f,%f,%f}\n", p1[0], p1[1], p1[2]);
+        printf("n: {%f,%f,%f}\n", n[0], n[1], n[2]);
+        printf("p2: {%f,%f,%f}\n", p2[0], p2[1], p2[2]);
+        printf("PointPlaneIntersection\n");
+    }
+    if(fabs(*dis) < EPS2){
+        return(POINT_INTERSECTION);
+    }else{
+        return(NO_INTERSECTION);
+    }
+}
+
+int PointPointIntersection(double *p1, double *p2, double *vec, double *dis){
+    vec[0] = p2[0] - p1[0];
+    vec[1] = p2[1] - p1[1];
+    vec[2] = p2[2] - p1[2];
+
+    *dis = sqrt(vec[0]*vec[0] + vec[1]*vec[1] + vec[2]*vec[2]);
+    if(std::isnan(*dis) || std::isinf(*dis)){
+        printf(" PointPointIntersection: %f\n", *dis);
+    }
+    if(std::isnan(*dis)) *dis = 0.0;
+    if(std::isinf(*dis)) *dis = 0.0;
+    
+    if(*dis < EPS2){
+        return(POINT_INTERSECTION);
+    }else{
+        return(NO_INTERSECTION);
+    }
+}
+
+
+void cross(real8 a[3], real8 b[3], real8 c[3])
+{
+    c[0]=a[1]*b[2]-a[2]*b[1];
+    c[1]=a[2]*b[0]-a[0]*b[2];
+    c[2]=a[0]*b[1]-a[1]*b[0];
+}
+
+real8 Normal(real8 a[3])
+{
+        return( sqrt(a[0]*a[0]+a[1]*a[1]+a[2]*a[2]) );
+}
diff --git a/src/Parse.cpp b/src/Parse.cpp
index 8fc563b0b1c7dcba132141ae86ed1a540849ac7a..d4af27eae048f603ad241895c16687860e5e24b5 100644
--- a/src/Parse.cpp
+++ b/src/Parse.cpp
@@ -2,6 +2,7 @@
 #include <dirent.h>
 #include <unistd.h>
 #include <sys/stat.h>
+#include <regex>
 
 #include "Home.h"
 #include "Util.h"
@@ -22,7 +23,6 @@ vector<string> GetFiles(const string &file)
     
     memset(basePath, '\0', sizeof(basePath));
     getcwdReturn = getcwd(basePath, sizeof(basePath));
-//    bpath = basePath;
 
     iPos = f.find_last_of('/');
     if(iPos == string::npos){
@@ -32,7 +32,7 @@ vector<string> GetFiles(const string &file)
         fdir = curDir + f.substr(0,iPos);
         fname = f.substr(iPos+1);
     }
-
+    std::regex regex(fname);
     
     if((dir = opendir(fdir.c_str())) == NULL){
         Fatal("Can not open dir %s", fdir.c_str());
@@ -44,8 +44,7 @@ vector<string> GetFiles(const string &file)
           continue;
         }else if(ptr->d_type == 8){   //file
             dname = ptr->d_name;
-            int loc = dname.find(fname, 0); 
-            if(loc != string::npos){
+            if(std::regex_match(dname, regex)){
                 strs.push_back(fdir + "/" + dname);
             }
         }else if(ptr->d_type == 10){    // link file
diff --git a/src/ReadData.cpp b/src/ReadData.cpp
index e56240a0bd95beed64484672cbfdb91fb8dfa6a6..dac84087ecd4e57ed2b54217f60ab7f4ac5e80da 100644
--- a/src/ReadData.cpp
+++ b/src/ReadData.cpp
@@ -1,4 +1,6 @@
 
+#include <strings.h>
+#include <regex>
 #include "Home.h"
 #include "Util.h"
 #include "ProDataIO.h"
@@ -6,86 +8,183 @@
 #include "DDD.h"
 #include "MD.h"
 
+#define     MAXLINELENGTH 1024
+
 using namespace std;
+const char* strstri(const char* str, const char* subStr)
+{
+    int len = strlen(subStr);
+    if(len == 0)
+    {
+        return NULL;
+    }
+    while(*str)
+    {
+        if(strncasecmp(str, subStr, len) == 0)
+        {
+            return str;
+        }
+        ++str;
+    }
+    return NULL;
+}
 
 int ReadTecplotNormalData(string &file, Table_t &table, string &secLine)
 {
-    int         lineNo = 1, i, j, type;
-    string      str, v("v");
-    ifstream    infile;
+    int  i, numPoints = 1, j;
+    char str[MAXLINELENGTH], delim[] = " \t\",\n", *p, *token;
+    char quotation[]="\"\n";
+    char name[256], *rt;
+    FILE *fp;
+    if((fp = fopen(file.c_str(), "r+")) == NULL){
+        printf("Warning: can not open fie %s\n", file.c_str());
+        return(0);
+    }else{
+//        printf("Reading tecplot file %s ...\n", file.c_str());
+    }
 
-    vector<string>  line;
-    vector<double>  col;
-    infile.open(file.c_str());
+    InitTable(table);
+    secLine = "";
 
-    if(!infile){
-        printf("Warning: cant not the file %s", file.c_str());
-        return (0);
-    }
+    std::regex aux_equation("AUXDATA\\s+\\S+\\s*=\\s*[-\"\'+.a-zA-Z0-9]+");
+    std::regex equation("\\S+\\s*=\\s*[-\"\'+.a-zA-Z0-9]+");
+    std::regex var_name("[a-zA-Z]+[a-zA-z0-9_\\-]*");
+    std::regex value("[a-zA-z0-9_\\.\\+\\-]+");
+    std::smatch equation_match;
+    std::smatch v_match;
 
-    vector<string>().swap(table.variables);
-    vector<vector<real8> >().swap(table.data);
-/*
- *  Read the first line.
- */
-    getline(infile, str);
-    if(str.empty()){
-        printf("Warning: there is noting in the file %s", file.c_str());
-        return (0);
-    }
+    std::string buff0, buff1, buff2;
 
-    line = split(str, " ");
-    for(i=0; i<line.size(); i++){
-        if(atof(line[i].c_str()) == 0)break;
-    }
-    vector<string>().swap(line);
+    Variable_t  auxVar;
 
-    if(i == 0){
-        line = split(str, "=");
-        table.variables = split(line[1], ",");
-        for(i=0; i<table.variables.size(); i++){
-            WashString(table.variables[i]);
-        }
-        vector<string>().swap(line);
-        col.resize(table.variables.size());
+    bool    firstPoint = 1;
+    int     currSize;
+    char    *p2;
+    while(1){
+        rt = fgets(str, MAXLINELENGTH, fp);
+        if(feof(fp))break;
+        if(str == "\n" || str == NULL)continue;
 
-    }else{
-        line = split(str, " ");
-        table.variables.resize(line.size());
-        col.resize(table.variables.size());
-        for(i=0;i<line.size();i++){
-            table.variables[i] = v + (char)i;
-            col[i] = atof(line[i].c_str());
-        }
-        table.data.push_back(col);
-    }
+        if(firstPoint){
+            if(strstr(str, "variables") != NULL){
+                token = strtok(str, "=");
+                while(1){
+                    if((token = strtok(NULL, delim)) != NULL){
+                        table.variables.push_back(token);
+                    }else{
+                        break;
+                    }
+//                    table.variables.push_back(strtok(NULL, delim));
+                }
+//                printf("variables: ");
+//                for(i=0;i<table.variables.size();i++)printf("%s ", table.variables[i].c_str());
+//                printf("\n");
+                continue;
+            }
 
+            buff0 = str;
+            if(std::regex_search(buff0, equation_match, aux_equation)){
+                buff1 = equation_match[0].str();
+                if(std::regex_search(buff1, v_match, var_name)){
+                    buff2 = v_match.suffix().str();
+                    if(std::regex_search(buff2, v_match, var_name)){
+                        buff1 = v_match[0].str();
+                        buff2 = v_match.suffix().str();
+                        if(std::regex_search(buff2, v_match, value)){
+                            table.aux[buff1] = v_match[0].str();
+                        }
+                    }
+                }
+                continue;
+            }else if(std::regex_search(buff0,equation_match, equation)){
+                 buff1 = buff0;
+                 while(std::regex_search(buff1,equation_match, equation)){
+                    buff2 = equation_match[0].str();
+                    if(std::regex_search(buff2,v_match, var_name)){
+                        buff2 = v_match[0].str();
+                        buff1 = v_match.suffix().str();
+                        if(std::regex_search(buff1, v_match, value)){
+                            if(buff2 == "T"){
+                                table.T = v_match[0].str(); 
+                            }else if (buff2 == "F"){
+                                table.F = v_match[0].str(); 
+                            }else if (buff2 == "I" || buff2 == "i"){
+                                table.i = atoi(v_match[0].str().c_str());
+                            }else if (buff2 == "J" || buff2 == "j"){
+                                table.j = atoi(v_match[0].str().c_str());
+                            }else if (buff2 == "K" || buff2 == "k"){
+                                table.k = atoi(v_match[0].str().c_str());
+                            }else if (buff2 == "SOLUTIONTIME" || buff2 == "solutiontime"){
+                                table.solutionTime = atof(v_match[0].str().c_str());
+                            }else{
+                                printf("Warning: can not identify %s", buff2.c_str());
+                            }
+                        }
+                    }
+                    buff1 = equation_match.suffix().str();
+                }
+                continue;
+            }
 
-    j = 0;
-    while(getline(infile,str))
-    {   
-        lineNo++;
-        line = split(str, " ");
-        if(lineNo == 2)secLine = str;
-       
-        if(line.size() != table.variables.size()){
-            continue; 
+            firstPoint = 0;
+            if(table.variables.size() == 0){
+                char strBak[MAXLINELENGTH];
+                sprintf(strBak, "%s", str);
+                i = 0;
+                token = strtok(str, " \n");
+                if(token != NULL){
+                    snprintf(name, sizeof(name), "v%d", i); 
+                    table.variables.push_back(name);
+                    
+                    while((token = strtok(NULL, " \n")) != NULL){
+                        i++;
+                        snprintf(name, sizeof(name), "v%d", i); 
+                        table.variables.push_back(name);
+                    } 
+                    printf("Undifined variables: ");
+                    for(i=0;i<table.variables.size();i++)printf("%s ", table.variables[i].c_str());
+                    printf("\n");
+                }
+                sprintf(str, "%s", strBak);
+            }
+            numPoints = table.i*table.j*table.k;
         }
-        for(i=0; i<col.size(); i++){
-            col[i] = atof(line[i].c_str());
+
+
+//        if(strtok(str, " ") == NULL)continue;
+        if(numPoints > 1){
+            currSize = (int)table.data.size();
+            table.data.resize(currSize+numPoints);
+            for(i=currSize; i<table.data.size() && !feof(fp); i++){
+                table.data[i].resize(table.variables.size());
+                table.data[i][0] = atof(strtok(str, " ")); 
+                for(j=1; j<table.variables.size(); j++){
+                    table.data[i][j] = atof(strtok(NULL, " \n"));
+                }
+                rt = fgets(str, MAXLINELENGTH, fp);
+            }
+        }else{
+            currSize = (int)table.data.size();
+            table.data.resize(currSize+numPoints);
+            table.data[currSize].resize(table.variables.size());
+            table.data[currSize][0] = atof(strtok(str, " ")); 
+            for(j=1; j<table.variables.size(); j++){
+                table.data[currSize][j] = atof(strtok(NULL, " \n"));
+            }
+//            printf("%d: ", currSize);
+//            for(i=0;i<table.data[currSize].size(); i++)printf("%f ", table.data[currSize][i]);
+//            printf("\n");
         }
-        table.data.push_back(col);
-        j++;
-        vector<string>().swap(line);
     }
-    infile.close();
+    fclose(fp);
 
-//    printf("Finsish reading input file %s\n", file.c_str());    
+    if((int)table.data.size() != table.i*table.j*table.k)table.i = (int)table.data.size();
+    printf("Finish reading input file %s, %d points\n", file.c_str(), currSize+1);    
     return 1;            
 }
 
 
-int ReadMGDataFile(const string &file, MgData_t &mgdata)
+int ReadDumpFile(const string &file, Dump_t &dum)
 {
     int         i, j, atoms;
     int         idCol = 0, typeCol = 0;
@@ -93,7 +192,7 @@ int ReadMGDataFile(const string &file, MgData_t &mgdata)
     ifstream    infile;
     string      str;
 
-    CleanMgData(mgdata);
+    CleanDump(dum);
 
     infile.open(file.c_str());
 
@@ -114,7 +213,7 @@ int ReadMGDataFile(const string &file, MgData_t &mgdata)
         if(words[0] == "ITEM:"){
             if(words[1] == "TIMESTEP"){
                 getline(infile,str);
-                mgdata.timestep = atoi(str.c_str());
+                dum.timestep = atoi(str.c_str());
             }else if(words[1] == "NUMBER"){
                 if(words.size()==4){
                     if(words[2] == "OF" && words[3] == "ATOMS"){
@@ -125,10 +224,12 @@ int ReadMGDataFile(const string &file, MgData_t &mgdata)
             }else if(words[1] == "BOX"){
                 if(words.size()==6){
                     if(words[2] == "BOUNDS"){
-                        mgdata.bounds.resize(3);
-                        mgdata.bounds[0] = words[3];
-                        mgdata.bounds[1] = words[4];
-                        mgdata.bounds[2] = words[5];
+                        dum.bounds.resize(3);
+                        dum.bounds[0] = words[3];
+                        dum.bounds[1] = words[4];
+                        dum.bounds[2] = words[5];
+                        
+                        dum.box.resize(3);
                         for(i=0; i<3; i++){
                             vector<string>().swap(subwords);
                             getline(infile,str);
@@ -136,19 +237,40 @@ int ReadMGDataFile(const string &file, MgData_t &mgdata)
                             if(subwords.size() != 2){
                                 Fatal("in file %s, can not read %s", file.c_str(), str.c_str());
                             }
-                            mgdata.box[i*2] = atof(subwords[0].c_str());
-                            mgdata.box[1+i*2] = atof(subwords[1].c_str());
+                            dum.box[i].resize(2);
+                            dum.box[i][0] = atof(subwords[0].c_str());
+                            dum.box[i][1] = atof(subwords[1].c_str());
+                        }
+                    }
+                }else if(words.size() == 9){
+                    dum.bounds.resize(6);
+                    for(i=0; i<6; i++){
+                        dum.bounds[i] = words[3+i];
+                    }
+                    dum.box.resize(3);
+
+                    for(i=0; i<3; i++){
+                        vector<string>().swap(subwords);
+                        getline(infile,str);
+                        subwords = split(str, " ");
+                        if(subwords.size() != 3){
+                            Fatal("in file %s, can not read %s", file.c_str(), str.c_str());
                         }
+                        dum.box[i].resize(3);
+                        dum.box[i][0] = atof(subwords[0].c_str());
+                        dum.box[i][1] = atof(subwords[1].c_str());
+                        dum.box[i][2] = atof(subwords[2].c_str());
+
                     }
                 } 
              }else if(words[1] == "ATOMS"){
-                if(words.size()<=7){
-                    printf("Waring: in file %s, can not read %s", file.c_str(), str.c_str());
+                if(words.size()<7){
+                    printf("Waring: in file %s, can not read %s - 1\n", file.c_str(), str.c_str());
                     return(0);
                 }
 
                 varCol.resize(words.size()-7);
-                mgdata.variables.resize(words.size()-7);
+                dum.variables.resize(words.size()-7);
                 j = 0;
                 for(i=2; i<words.size(); i++){
                     if(words[i] == "x"){xCol = i-2;        continue;}
@@ -156,30 +278,30 @@ int ReadMGDataFile(const string &file, MgData_t &mgdata)
                     if(words[i] == "z"){zCol = i-2;        continue;}
                     if(words[i] == "id"){idCol = i-2;      continue;}
                     if(words[i] == "type"){typeCol = i-2;  continue;}
-                    mgdata.variables[j] = words[i];
+                    dum.variables[j] = words[i];
                     varCol[j] = i-2;
                     j++;
                 } 
 
-                mgdata.atom.resize(atoms);
-                for(i=0; i<mgdata.atom.size(); i++){
+                dum.atom.resize(atoms);
+                for(i=0; i<dum.atom.size(); i++){
                     vector<string>().swap(subwords);
                     getline(infile,str);
                     subwords = split(str, " ");
-                    if(subwords.size() != mgdata.variables.size()+5){
-                        printf("Warning: in file %s, can not read %s", file.c_str(), str.c_str());
+                    if(subwords.size() != dum.variables.size()+5){
+                        printf("Warning: in file %s, can not read %s - 2\n", file.c_str(), str.c_str());
                         return (0);
                     }
 
-                    mgdata.atom[i].x = atof(subwords[xCol].c_str());
-                    mgdata.atom[i].y = atof(subwords[yCol].c_str());
-                    mgdata.atom[i].z = atof(subwords[zCol].c_str());
-                    mgdata.atom[i].id = atoi(subwords[idCol].c_str());
-                    mgdata.atom[i].type = atoi(subwords[typeCol].c_str());
+                    dum.atom[i].x = atof(subwords[xCol].c_str());
+                    dum.atom[i].y = atof(subwords[yCol].c_str());
+                    dum.atom[i].z = atof(subwords[zCol].c_str());
+                    dum.atom[i].id = atoi(subwords[idCol].c_str());
+                    dum.atom[i].type = atoi(subwords[typeCol].c_str());
 
-                    mgdata.atom[i].vars.resize(mgdata.variables.size());
-                    for(j=0; j<mgdata.variables.size(); j++){
-                        mgdata.atom[i].vars[j] = atof(subwords[varCol[j]].c_str());
+                    dum.atom[i].vars.resize(dum.variables.size());
+                    for(j=0; j<dum.variables.size(); j++){
+                        dum.atom[i].vars[j] = atof(subwords[varCol[j]].c_str());
                     }
                 }
             }else{
@@ -247,3 +369,83 @@ int ReadDataFromMDLogFile(const vector<string> &files, LineList_t &list)
 }
 
 
+bool ReadLMPFile(const string file, Dump_t &dum)
+{
+    CleanDump(dum);
+
+    FILE *fp;
+    if((fp = fopen(file.c_str(), "r+")) == NULL){
+        return(0);
+    }
+
+    char str[MAXLINELENGTH], delim[] = " \t", *p, *rt;
+    dum.box.resize(3);
+    dum.bounds.resize(3);
+    for(int i = 0; i<3; i++)dum.box[i].resize(2);
+
+    while(!feof(fp)){
+        rt = fgets(str, MAXLINELENGTH, fp);
+        
+        if((p = strstr(str, "atoms")) != NULL){
+            dum.atom.resize(atoi(str));
+            printf("atoms: %d\n",(int)dum.atom.size());
+            continue;
+        }
+
+        if((p = strstr(str, "xlo xhi")) != NULL){
+            dum.box[0][0] = atof(strtok(str, delim));
+            dum.box[0][1] = atof(strtok(NULL, delim));
+            printf("box x: [%f, %f]\n", dum.box[0][0], dum.box[0][1]);
+            dum.bounds[0] = "xx";
+            continue;
+        }
+
+        if((p = strstr(str, "ylo yhi")) != NULL){
+            dum.box[1][0] = atof(strtok(str, delim));
+            dum.box[1][1] = atof(strtok(NULL, delim));
+            dum.bounds[1] = "yy";
+            printf("box y: [%f, %f]\n", dum.box[1][0], dum.box[1][1]);
+            continue;
+        }
+
+        if((p = strstr(str, "zlo zhi")) != NULL){
+            dum.box[2][0] = atof(strtok(str, delim));
+            dum.box[2][1] = atof(strtok(NULL, delim));
+            dum.bounds[2] = "zz";
+            printf("box z: [%f, %f]\n", dum.box[2][0], dum.box[2][1]);
+            continue;
+        }
+
+        if((p = strstr(str, "xy xz yz")) != NULL){
+            dum.box[0].resize(3);
+            dum.box[0].resize(3);
+            dum.box[0].resize(3);
+            dum.bounds.resize(6);
+
+            dum.box[0][2] = atof(strtok(str, delim));
+            dum.box[1][2] = atof(strtok(NULL, delim));
+            dum.box[2][2] = atof(strtok(NULL, delim));
+            dum.bounds[3] = "xy";
+            dum.bounds[4] = "xz";
+            dum.bounds[5] = "yz";
+
+            printf("triclinic box: [%f, %f %f]\n", dum.box[0][2], dum.box[1][2], dum.box[2][2]);
+            continue;
+        }
+
+        if((p = strstr(str, "Atoms")) != NULL){
+            rt = fgets(str, MAXLINELENGTH, fp);
+            for(int i=0; i<dum.atom.size(); i++){
+                rt = fgets(str, MAXLINELENGTH, fp);
+                dum.atom[i].id = atoi(strtok(str, delim));
+                dum.atom[i].type = atoi(strtok(NULL, delim));
+                dum.atom[i].x = atof(strtok(NULL, delim));
+                dum.atom[i].y = atof(strtok(NULL, delim));
+                dum.atom[i].z = atof(strtok(NULL, delim));
+            }
+        }
+    }
+
+    fclose(fp);
+    return 1;
+}
diff --git a/src/Util.cpp b/src/Util.cpp
index 75ab6360ab763a2a0243e7fc9ebd5b55d66f89e4..7108913c193aee1b4c497adb19178a74731a39f3 100644
--- a/src/Util.cpp
+++ b/src/Util.cpp
@@ -121,7 +121,6 @@ vector<double>  GenerateSequence(double from, double to, double meshSize)
     return(seq);
 }
 
-
 void FoldBox(real8 boundMin[3], real8 boundMax[3], real8 *x, real8 *y, real8 *z)
 {
         real8   xc, yc, zc;
@@ -257,21 +256,19 @@ void SwapLineList(LineList_t &list){
 }
 
 
-void CleanMgData(MgData_t &mg)
+void CleanDump(Dump_t &dum)
 {
     int i;
-    mg.timestep = 0;
-    for(i=0; i<6; i++){
-        mg.box[i] = 0.0;
-    }
+    dum.timestep = 0;
 
-    vector<string>().swap(mg.bounds);
-    vector<string>().swap(mg.variables);
+    vector<vector<double> >().swap(dum.box);
+    vector<string>().swap(dum.bounds);
+    vector<string>().swap(dum.variables);
     
-    for(i=0; i<mg.atom.size(); i++){
-        vector<double>().swap(mg.atom[i].vars);
+    for(i=0; i<dum.atom.size(); i++){
+        vector<double>().swap(dum.atom[i].vars);
     }
-    vector<Atom_t>().swap(mg.atom);
+    vector<Atom_t>().swap(dum.atom);
 
     return;
 }
@@ -334,100 +331,26 @@ void StitchTecplotData(vector<Table_t> &tables, Table_t &table, int eigenID)
     return;            
 }
 
-void SpecifyEquations(Table_t &table)
-{
-#if 1
-    int     i, j;
-
-    int colBurgID =  GetColIDFromTable(table, "burgID");
-    int colIndex =  GetColIDFromTable(table, "index");
-
-    for(i=0; i<table.data.size(); i++){
-        real8 &burgID = table.data[i][colBurgID];
-        real8 &index = table.data[i][colIndex];
-
-        if(burgID > 6.5){
-            if(burgID == 7 || burgID == 12 || burgID == 17){
-                burgID = -1;
-            }else{
-                burgID = 7;
-            }
-        }
-    }
-
-#else
-    int     i, j;
-    bool    changed;
-    real8   tauAve, diff;
-    
-    vector<real8>   vec(3);
-    table.variables.push_back("tau_averaged_single");
-    table.variables.push_back("diff");
-
-    int colSigmaS1 = GetColIDFromTable(table, "sigma_sin1");
-    int colSS1 = GetColIDFromTable(table, "s_sin1");
-    int colTauS1 = GetColIDFromTable(table, "tau_sin1");
-
-    int colSigmaS2 = GetColIDFromTable(table, "sigma_sin2");
-    int colSS2 = GetColIDFromTable(table, "s_sin2");
-    int colTauS2 = GetColIDFromTable(table, "tau_sin2");
-
-    int colTauBic = GetColIDFromTable(table, "tau_bic");
-
-    for(i=0; i<table.data.size(); i++){
-
-        real8 &sigmaS1 = table.data[i][colSigmaS1];
-        real8 &sS1 = table.data[i][colSS1];
-        real8 &tauS1 = table.data[i][colTauS1];
-
-        real8 &sigmaS2 = table.data[i][colSigmaS2];
-        real8 &sS2 = table.data[i][colSS2];
-        real8 &tauS2 = table.data[i][colTauS2];
-        real8 tauBic = table.data[i][colTauBic];
-
-        changed = ((sigmaS1 > sigmaS2) ? 1 : 0);
-        if(changed){
-            vec[0] = sigmaS1;
-            vec[1] = sS1;
-            vec[2] = tauS1;
-
-            sigmaS1 = sigmaS2;
-            sS1 = sS2;
-            tauS1 = tauS2;
-
-            sigmaS2 = vec[0];
-            sS2 = vec[1];
-            tauS2 = vec[2];
-        }
-
-        if(sS1 == sS2){
-            tauAve = 0.5*(tauS1 + tauS2);
-        }else{
-            if(sS1 > sS2){
-                tauAve = tauS1;
-            }else{
-                tauAve = tauS2;
-            }
-        }
-
-        diff = tauBic - tauAve; 
-
-        table.data[i].push_back(tauAve);
-        table.data[i].push_back(diff);
-    }
-#endif
-    return;
-}
-
 
 void SpecifyEquations_PLTDATA(InArgs_t *inArgs)
 {
-    int     i;
+    int     i, index;
     bool    readState;
-    string  secLine;
+    string  secLine, noBackupName = ("nobackup"), backupFile, sufBack = (".bak");
+    char    name[256];
+    bool    backup = 1;
 
     vector<Table_t> tables;
 
+    if((index = GetValID(inArgs->priVars, noBackupName)) < inArgs->priVars.size()){
+        backup = 0;
+    }
+    if(backup){
+        printf("Backup: Yes.\n");
+    }else{
+        printf("Backup: NO! (nobackup)\n");
+    }
+
     if(inArgs->help)return;
 
     if(inArgs->inpFiles.size() == 0){
@@ -438,12 +361,58 @@ void SpecifyEquations_PLTDATA(InArgs_t *inArgs)
     for(i=0; i<inArgs->inpFiles.size(); i++){
         readState = ReadTecplotNormalData(inArgs->inpFiles[i], tables[i], secLine);
         if(!readState)continue;
-
+        if(backup){
+            backupFile = inArgs->inpFiles[i] + sufBack;
+            WriteTecplotNormalData(tables[i], backupFile, 10, secLine); 
+        }
         SpecifyEquations(tables[i]);
     
-        WriteTecplotNormalData(tables[i], inArgs->outFiles[i], 10, secLine); 
+        WriteTecplotNormalData(tables[i], inArgs->inpFiles[i], 10, secLine); 
     }
 
     return;
 
 }
+
+void FormatVector(real8 vec[3], const char *msg){
+	printf("%s vector:\n",msg);
+	printf("{%.15f,%.15f,%.15f}\n", vec[0], vec[1], vec[2]);
+}
+
+void InitList(LineList_t &list){
+    vector<string>().swap(list.variables);
+    vector<vector<real8> >().swap(list.data);
+    list.aux.clear();
+    list.i=1;
+    list.j=1;
+    list.k=1;
+    list.solutionTime = -1;
+    list.T="";
+    list.F = "Point";
+    return;
+}
+
+void InitTable(Table_t &table){
+    vector<string>().swap(table.variables);
+    vector<vector<real8> >().swap(table.data);
+    table.aux.clear();
+    table.i=1;
+    table.j=1;
+    table.k=1;
+    table.solutionTime = -1;
+    table.T="";
+    table.F = "Point";
+    return;
+}
+
+bool FindLinearPart(real8 (*line)[3], const int nums, int range[2])
+{
+    if(nums < 3)return 0;
+
+    int i, j;
+
+    return 0;
+} 
+
+
+
diff --git a/src/WriteData.cpp b/src/WriteData.cpp
index 134c02da303906e0c5fa739b752c08ef0ea139bd..8a39034829ace786e751e8fbf0336f008cfb5a8f 100644
--- a/src/WriteData.cpp
+++ b/src/WriteData.cpp
@@ -1,9 +1,11 @@
+#include <ctime>
 
 #include "Home.h"
 #include "Util.h"
 #include "ProDataIO.h"
 #include "MD.h"
 #include "DDD.h"
+#include "Math.h"
 
 using namespace std;
 
@@ -12,7 +14,6 @@ void WriteTecplotNormalData(const LineList_t &list, const string &file, double p
     int     i, j;
     string  fn(file), plt(".plt");
 
-
     if(fn.length() > 4){
         int loc = fn.find(plt, fn.length()-5); 
         if(loc == string::npos){
@@ -34,6 +35,22 @@ void WriteTecplotNormalData(const LineList_t &list, const string &file, double p
         }
     }
 
+    if(((int)list.data[0].size()) != list.i*list.j*list.k){
+        Fatal("somthing wrong with the list size %d != %dx%dx%d",
+              (int)list.data[0].size(), list.i, list.j, list.k);
+    }
+
+    out << "zone i = " << list.i << ", j = " << list.j << ", k = " << list.k;
+    if(list.T.length() != 0) out << ", T = \"" << list.T << "\"";
+    if(list.solutionTime >=0) out << ", SOLUTIONTIME = " << list.solutionTime;
+    out << ", F = " << list.F << "\n";
+
+    if(!list.aux.empty()){
+        for(const auto &pair : list.aux){
+            out << "AUXDATA " << pair.first << " = \"" << pair.second << "\"\n";
+        }
+    }
+
     for(i=0; i<list.data[0].size(); i++){
         for(j=0; j<list.variables.size(); j++){
             out << setprecision(precision) <<list.data[j][i] << " ";
@@ -72,8 +89,19 @@ void WriteTecplotNormalData(const Table_t &table, const string &file, double pre
         }
     }
 
-    if(secLine != ""){
-        out << secLine << endl;
+    if(((int)table.data.size()) != table.i*table.j*table.k){
+        Fatal("somthing wrong with the table size %d != %dx%dx%d",
+              (int)table.data.size(), table.i, table.j, table.k);
+    }
+    out << "zone i = " << table.i << ", j = " << table.j << ", k = " << table.k;
+    if(table.T.length() != 0) out << ", T = \"" << table.T << "\"";
+    if(table.solutionTime >=0) out << ", SOLUTIONTIME = " << table.solutionTime;
+    out << ", F = " << table.F << "\n";
+
+    if(!table.aux.empty()){
+        for(const auto &pair : table.aux){
+            out << "AUXDATA " << pair.first << " = \"" << pair.second << "\"\n";
+        }
     }
 
     for(i=0; i<table.data.size(); i++){
@@ -87,3 +115,111 @@ void WriteTecplotNormalData(const Table_t &table, const string &file, double pre
 //    printf("Finsish writing output file %s\n", fn.c_str());    
     return;
 }
+
+
+int WriteDumpFile(const string &file, Dump_t &dum, int precision) 
+{
+    int     i, j;
+    string  fn(file), plt(".dum");
+
+    if(fn.length() > 4){
+        int loc = fn.find(plt, fn.length()-5); 
+        if(loc == string::npos){
+            fn += ".dum";
+        }
+    }else{
+        fn += ".dum";
+    }
+
+    ofstream out;
+    out.open(fn.c_str(), ios::out);
+    
+    out << "ITEM: TIMESTEP" << endl;
+    out << dum.timestep << endl;
+
+    out << "ITEM: NUMBER OF ATOMS" << endl;
+    out << dum.atom.size() << endl;
+    
+    out << "ITEM: BOX BOUNDS";
+    for(i=0; i<dum.bounds.size(); i++){
+        out << " " << dum.bounds[i];
+    }
+    out << endl;
+
+    for(i=0; i<dum.box.size(); i++){
+        for(j=0; j<dum.box[i].size(); j++){
+            out <<  setprecision(precision) << dum.box[i][j] << " ";
+        }
+        out << endl;
+    }
+
+    out << "ITEM: ATOMS id type x y z";
+    for(i=0; i<dum.variables.size(); i++){
+        out << " " << dum.variables[i];
+    }
+    out << endl;
+
+    for(i=0; i<dum.atom.size(); i++){
+        out << dum.atom[i].id << " ";
+        out << dum.atom[i].type << " ";
+        out << setprecision(precision) << dum.atom[i].x << " ";
+        out << setprecision(precision) << dum.atom[i].y << " ";
+        out << setprecision(precision) << dum.atom[i].z << " ";
+        
+        for(j=0; j<dum.atom[i].vars.size(); j++){
+            out << setprecision(precision) << dum.atom[i].vars[j] << " "; 
+        }
+        out << endl;
+    }
+
+    out.close();
+    return(0);
+}
+
+
+
+int MGToLMPDataFile(const string &file, Dump_t &dum, int precision)
+{
+    int     i, j, nTypes = 0;
+    string  fn(file), plt(".lmp");
+
+    if(fn.length() > 4){
+        int loc = fn.find(plt, fn.length()-5); 
+        if(loc == string::npos){
+            fn += ".lmp";
+        }
+    }else{
+        fn += ".lmp";
+    }
+
+    ofstream out;
+    out.open(fn.c_str(), ios::out);
+
+    time_t  now = time(0);
+    char* dt = ctime(&now);
+
+    out << "LMMPS: " << dt << endl;
+
+    out << dum.atom.size() << " atoms" << endl << endl;
+
+    for(i=0; i<dum.atom.size(); i++){
+        nTypes = MAX(nTypes, dum.atom[i].type);
+    }
+    out << nTypes << " atom types" << endl << endl;
+
+    out << setprecision(precision) << dum.box[0][0] << " " << dum.box[0][1] << " xlo xhi" << endl;
+    out << setprecision(precision) << dum.box[1][0] << " " << dum.box[1][1] << " ylo yhi" << endl;
+    out << setprecision(precision) << dum.box[2][0] << " " << dum.box[2][1] << " zlo zhi" << endl;
+
+    out << endl << "Atoms" << endl << endl;
+
+    for(i=0; i<dum.atom.size(); i++){
+        out << dum.atom[i].id << " "<< dum.atom[i].type << " ";
+        out << setprecision(precision) << dum.atom[i].x << " ";
+        out << setprecision(precision) << dum.atom[i].y << " ";
+        out << setprecision(precision) << dum.atom[i].z << endl;
+    }
+
+    out.close();
+    return 0;
+}