Merge branch 'master' of deanwei:ansijing/prodata (a7285031) · Commits · DeAn Wei / prodata

include/MD.h

+0 −1

Original line number	Diff line number	Diff line
		@@ -28,7 +28,6 @@ typedef struct{

		typedef struct {
		int timestep;
		int atoms;
		real8 box[6];

		vector<string> bounds;

include/Util.h

+1 −1

Original line number	Diff line number	Diff line
		@@ -41,5 +41,5 @@ real8 LinearInterpolation(const Curve_t &curve, real8 x, real8 min = -1, real8 m
		void SwapTable(Table_t &table);
		void SwapLineList(LineList_t &list);
		void CleanMgData(MgData_t &mg);

		void NormalizeVec(real8 vec[3]);
		#endif

src/HandleExtendedDislocation.cpp

+247 −21

Original line number	Diff line number	Diff line
		@@ -140,14 +140,17 @@ void HandleExtendedDislocation_DDD(InArgs_t *inArgs)
		separation = 0.0;
		position = 0.0;
		for(i=0; i<seq.size(); i++){
		data[1][i] = LinearInterpolation(curve1, seq[i], boundMin[0], boundMax[0]);
		data[2][i] = LinearInterpolation(curve2, seq[i], boundMin[0], boundMax[0]);
		data[1][i] = linearinterpolation(curve1, seq[i], boundmin[0], boundmax[0]) - cubel;
		data[2][i] = linearinterpolation(curve2, seq[i], boundmin[0], boundmax[0]) - cubel;
		data[3][i] = 0.5 * (data[1][i] + data[2][i]);
		if(i<seq.size()-2){
		separation += fabs(data[1][i] - data[2][i]);
		position += data[3][i];
		}
		separation /= ((double)seq.size());
		position /= ((double)seq.size());
		}
		separation /= ((double)(seq.size()-1));
		position /= ((double)(seq.size()-1));
		position += cubel;

		for(i=0; i<seq.size(); i++){
		data[4][i] = 0.0;
		@@ -166,6 +169,9 @@ void HandleExtendedDislocation_DDD(InArgs_t *inArgs)

		for(i=0; i<seq.size(); i++){
		data[0][i] += 0.5*cubel;
		data[1][i] += cubel;
		data[2][i] += cubel;
		data[3][i] += cubel;
		data[4][i] /= ((double)seq.size());
		data[5][i] /= ((double)seq.size());
		data[6][i] /= ((double)seq.size());
		@@ -198,40 +204,193 @@ void HandleExtendedDislocation_DDD(InArgs_t *inArgs)
		return;
		}

		typedef struct {
		real8 x,y,z;
		vector<Atom_t *> nbr;
		}Probe_t;

		typedef struct {
		int timestep;
		double x,y,z;
		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){
		return(p.x<q.x);
		}else{
		return(p.z < q.z);
		}
		}else{
		return(p.y < q.y);
		}
		}

		bool com2(Probe_t p, Probe_t q)
		{
		return(p.y<p.y);
		}
		bool com3(EDState_t p, EDState_t q)
		{
		return(p.timestep<q.timestep);
		}

		void HandleExtendedDislocation_MD(InArgs_t *inArgs)
		{
		int i, j, file, index;
		double cutofflen = 2.556;
		int i, j, file, index, lastIndex, firstNbr, indexVar;
		int nums[2] = {18, 24}, logfile;
		int stepID;
		double cutofflen = 2.556, alpha = 0.005, beta = 1.0;
		double position[3], separation;
		MgData_t mg;
		LineList_t list;
		string cutoffName("cutoff"), dvarName("dvar"), dvar("c_vcna");
		double dir[3] = {0, 1, 0}, p0[3], dval=5;
		bool noP0 = 1;
		string paraName("para"), dvarName("dvar"), dvar("c_vcna"), p0Name("p0"), dirName("dir");
		string numsName("nums");

		if((index = GetValID(inArgs->priVars, cutoffName)) < inArgs->priVars.size()){
		vector<Atom_t>::iterator it;
		vector<EDState_t> states;
		EDState_t state;

		ofstream out;
		out.open(inArgs->outFiles[0].c_str(), ios::out);

		if((index = GetValID(inArgs->priVars, paraName)) < inArgs->priVars.size()){
		if(inArgs->priVars[index].vals.size() == 2){
		cutofflen = atof(inArgs->priVars[index].vals[0].c_str());
		alpha = atof(inArgs->priVars[index].vals[1].c_str());
		}
		printf("The cut-off (cutoff) length is %f\n", cutofflen);
		}
		printf("The parameters (para): cut-off length %f and resolution value %f\n", cutofflen, alpha);

		if((index = GetValID(inArgs->priVars, dvarName)) < inArgs->priVars.size()){
		dvar = inArgs->priVars[index].vals[0];
		if(inArgs->priVars[index].vals.size() == 2){
		dval = atof(inArgs->priVars[index].vals[1].c_str());
		}
		}
		printf("The determined var (dvar) is %s\n", dvar.c_str());
		printf("The determined var (dvar) is %s, value is %f\n", dvar.c_str(), dval);

		ReadDataFromMDLogFile(inArgs->auxFiles, list);
		if((index = GetValID(inArgs->priVars, p0Name)) < inArgs->priVars.size()){
		if(inArgs->priVars[index].vals.size() != 3){
		Fatal("You have to determine the initial point by -dp0 x y z");
		}
		p0[0] = atof(inArgs->priVars[index].vals[0].c_str());
		p0[1] = atof(inArgs->priVars[index].vals[1].c_str());
		p0[2] = atof(inArgs->priVars[index].vals[2].c_str());
		noP0 = 0;
		printf("The initial probe point (p0) is %f %f %f\n", p0[0], p0[1], p0[2]);
		}

		if((index = GetValID(inArgs->priVars, dirName)) < inArgs->priVars.size()){
		if(inArgs->priVars[index].vals.size() != 3){
		Fatal("at least 3 vals for %s", dirName.c_str());
		}
		dir[0] = atof(inArgs->priVars[index].vals[0].c_str());
		dir[1] = atof(inArgs->priVars[index].vals[1].c_str());
		dir[2] = atof(inArgs->priVars[index].vals[2].c_str());
		NormalizeVec(dir);
		}
		printf("The moving direction of probe is along %f, %f, %f\n", dir[0], dir[1], dir[2]);

		if((index = GetValID(inArgs->priVars, numsName)) < inArgs->priVars.size()){
		if(inArgs->priVars[index].vals.size() == 2){
		nums[0] = atoi(inArgs->priVars[index].vals[0].c_str());
		nums[1] = atoi(inArgs->priVars[index].vals[1].c_str());
		}
		}
		printf("The parameters (para): cut-off length %f and resolution value %f\n", cutofflen, alpha);
		logfile = ReadDataFromMDLogFile(inArgs->auxFiles, list);

		for(file=0; file<inArgs->inpFiles.size(); file++){
		ReadMGDataFile(inArgs->inpFiles[file], mg);

		for(i=0; i<mg.variables.size(); i++){
		if(dvar == mg.variables[i])break;
		}
		if((i=indexVar) == mg.variables.size()){
		Fatal("There is no %s in the mg 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){
		// vector<double>().swap(mg.atom[i].vars);
		mg.atom.erase(it);
		}
		}
		}

		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]);
		}

		sort(mg.atom.begin(), mg.atom.end(), com);

		vector<Probe_t> probes;
		Probe_t probe;
		double d = 0.0;

		probe.y = p0[1];
		probe.z = p0[2];

		lastIndex = 0;
		while(1){
		vector<Atom_t *>().swap(probe.nbr);
		probe.x += (dir[0]dcutofflen*alpha);
		probe.y += (dir[1]dcutofflen*alpha);
		probe.z += (dir[2]dcutofflen*alpha);

		firstNbr = 1;
		for(i=lastIndex; i<mg.atom.size(); i++){
		if((pow((mg.atom[i].y-probe.y), 2) +
		pow((mg.atom[i].z-probe.z), 2) < cutofflen*cutofflen)
		&& mg.atom[i].vars[indexVar] == dval){
		if(firstNbr){
		lastIndex = i;
		firstNbr = 0;
		}
		probe.nbr.push_back(&mg.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 0
		printf("Timestep %d, Atoms %d, Bouds %s %s %s, ",
		mg.timestep, mg.atoms, mg.bounds[0].c_str(),
		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("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++){
		@@ -239,8 +398,75 @@ void HandleExtendedDislocation_MD(InArgs_t *inArgs)
		}
		printf("\n");
		#endif
		for(i=0; i<probes.size(); i++){
		probes[i].x = 0.0;
		probes[i].y = 0.0;
		probes[i].z = 0.0;
		for(j=0; j<probes[i].nbr.size(); j++){
		probes[i].x += probes[i].nbr[j]->x;
		probes[i].y += probes[i].nbr[j]->y;
		probes[i].z += probes[i].nbr[j]->z;
		}
		probes[i].x /= (double)probes[i].nbr.size();
		probes[i].y /= (double)probes[i].nbr.size();
		probes[i].z /= (double)probes[i].nbr.size();

		}
		sort(probes.begin(), probes.end(), com2);

		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(logfile ==1){
		state.vars.resize(list.variables.size()-1);
		for(i=0; i<list.data[0].size(); i++){
		if(state.timestep == ((int)list.data[0][i]))break;
		}
		for(j=0; j<state.vars.size(); j++){
		state.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";
		if(logfile == 1){
		for(i=1; i<list.variables.size(); i++){
		out << ", " << list.variables[i];
		}
		}
		out << endl;

		for(i=0; i<states.size(); i++){
		out << states[i].timestep << " ";
		out << states[i].separation << " ";
		out << states[i].y << " ";
		for(j=0; j< states[i].vars.size(); j++){
		out << setprecision(10) << states[i].vars[j] << " ";
		}
		out << endl;

		// printf("%d %f %f %d\n", states[i].timestep, states[i].separation,
		// states[i].y, ((int)(states[i].vars.size())));
		}

		out.close();
		return;
		}

src/Main.cpp

+12 −8

Original line number	Diff line number	Diff line
		@@ -184,6 +184,7 @@ static void GetInArgs(int argc, char argv[], InArgs_t inArgs)
		* is wrong, so notify the user and terminate.
		*/
		if (argv[i][0] != '-') {
		printf("3\n");
		Usage(argv[0]);
		exit(1);
		}
		@@ -205,6 +206,7 @@ static void GetInArgs(int argc, char argv[], InArgs_t inArgs)

		if (j == OPT_MAX) {
		Usage(argv[0]);
		printf("1\n");
		exit(1);
		}

		@@ -214,6 +216,7 @@ static void GetInArgs(int argc, char argv[], InArgs_t inArgs)
		*/
		if (optList[j].optPaired) {
		if (i+1 >= argc) {
		printf("2\n");
		Usage(argv[0]);
		exit(1);
		} else {
		@@ -284,13 +287,6 @@ static void GetInArgs(int argc, char argv[], InArgs_t inArgs)
		swap(bakInps, inArgs->inpFiles);

		sort(inArgs->inpFiles.begin(), inArgs->inpFiles.end());
		if(inArgs->inpFiles.size()>1){
		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){
		strs.resize(inArgs->auxFiles.size());
		vector<string>().swap(bakInps);
		@@ -303,13 +299,21 @@ static void GetInArgs(int argc, char argv[], InArgs_t inArgs)
		}
		swap(bakInps, inArgs->auxFiles);
		}
		#if 0
		if(inArgs->inpFiles.size()>1){
		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){
		printf("Auxiliary file(s) :\n");
		for(i=0; i<inArgs->auxFiles.size(); i++){
		printf("%s \n", inArgs->auxFiles[i].c_str());
		}
		}

		#endif
		return;
		}

src/ReadData.cpp

+7 −4

Original line number	Diff line number	Diff line
		@@ -56,7 +56,7 @@ int ReadTecplotNormalData(string &file, Table_t &table, string &secLine)

		int ReadMGDataFile(const string &file, MgData_t &mgdata)
		{
		int i, j;
		int i, j, atoms;
		int idCol = 0, typeCol = 0;
		int xCol = 0, yCol = 0, zCol = 0;
		ifstream infile;
		@@ -87,7 +87,7 @@ int ReadMGDataFile(const string &file, MgData_t &mgdata)
		if(words.size()==4){
		if(words[2] == "OF" && words[3] == "ATOMS"){
		getline(infile,str);
		mgdata.atoms = atoi(str.c_str());
		atoms = atoi(str.c_str());
		}
		}
		}else if(words[1] == "BOX"){
		@@ -126,7 +126,7 @@ int ReadMGDataFile(const string &file, MgData_t &mgdata)
		j++;
		}

		mgdata.atom.resize(mgdata.atoms);
		mgdata.atom.resize(atoms);
		for(i=0; i<mgdata.atom.size(); i++){
		vector<string>().swap(subwords);
		getline(infile,str);
		@@ -203,6 +203,9 @@ int ReadDataFromMDLogFile(const vector<string> &files, LineList_t &list)
		printf("\n");
		}
		#endif
		if(list.data.size()==0){
		return(0);
		}
		return (1);
		}