OOPSE/utils/zsub.cpp

#include <cstdio>
#include <cmath>
#include <cstring>
#include <vector>

#include "CmdlineZsub.h"
#include "SimInfo.hpp"
#include "SimSetup.hpp"
#include "ZConsReader.hpp"

using namespace std;

/**
 * @author: Teng Lin
 * @date: Feb 25, 2004
 * @version: 1.0
 * @email: tlin@nd.edu
 * 
 * This program is a substitue of quickLate. It not only implemnet all of the functions in
 * quickLate, but also provide the functions of replacing atom types for z-constraint.
 * This implement is kind of wierd!!!
 */
vector<int> getAtomToZMol(SimInfo* info, vector<int>& zmolIndex); 
void replaceAtomType(SimInfo* info, vector<int>& zmolIndex, vector<double>& zmolPos, 
                                           vector<int> atomToZMol,double zconsTol);
void convertSSD(char* buffer, char* prefix, double pos[3], double RotMat[3][3], int dipole_given);
void dumpToXYZ(SimInfo* info, ofstream* xyzStream, int repeatX, int repeatY, int repeatZ,
                                 int periodicBox_given, int dipole_given, int water_given);
void matVecMul3(double a[3][3], double b[3]);

int main(int argc, char* argv[]){
  gengetopt_args_info args_info;
  SimSetup startMe;
  SimInfo* info;
  string bassFileName;
  string dumpFileName;
  string xyzFileName;
  int nFrame;
  ZConsReader* zReader;
  DumpReader* dumpReader;
  ofstream* xyzStream;
  int nZmol;
  double zconsTol;
  vector<int> zmolIndex;
  vector<double> zmolPos;
  vector<int> atomToZMol; 
  int havezcons;
  //parse the command line option
    if (cmdline_parser (argc, argv, &args_info) != 0)
      exit(1) ;


  //get the dumpfile name and bass file name
  if (args_info.input_given){
    dumpFileName = args_info.input_arg;
  }
  else{
    cerr << "Does not have input file name" << endl;
    exit(1);
  }
  bassFileName = dumpFileName;
  bassFileName = bassFileName.substr(0, bassFileName.rfind(".")) + ".bass";

  if (args_info.output_given){
    xyzFileName = args_info.output_arg;
  }
  else{
    xyzFileName = dumpFileName;
    xyzFileName = xyzFileName.substr(0, xyzFileName.rfind(".")) + ".xyz";
  }
  
  info = new SimInfo();
  startMe.setSimInfo(info );

  startMe.parseFile( bassFileName.c_str());

  startMe.createSim();

  dumpReader = new DumpReader(dumpFileName.c_str());
 
  nFrame = dumpReader->getNframes();

  if (args_info.replace_given && info->nZconstraints > 0){    
    
    zReader = new ZConsReader(info);
    nZmol = zReader->getNumZMol();
    zmolIndex = zReader->getZConsIndex();
    zmolPos = zReader->getZConsPos();
    
    atomToZMol = getAtomToZMol(info, zmolIndex);

    //get tolerance
    GenericData* data;
    DoubleData* toleranceData;
    data = info->getProperty(ZCONSTOL_ID);

    if(!data) {        
      cerr << "Zsub error: can not get tolerance \n" << endl;
      exit(1);
    }
    else{

      toleranceData = dynamic_cast<DoubleData*>(data);
       
      if(!toleranceData){
        cerr << "Zsub error: can not get tolerance \n" << endl;
        exit(1);
      }
      else{
        zconsTol = toleranceData->getData();
      }
    }    


  }

  xyzStream = new ofstream(xyzFileName.c_str());
  for (int i = 0; i < nFrame; i += args_info.frame_arg){
    dumpReader->readFrame(info, i);
    if (args_info.replace_given && info->nZconstraints > 0){    
      replaceAtomType(info, zmolIndex, zmolPos, atomToZMol, zconsTol);
    }
    dumpToXYZ(info, xyzStream,args_info.repeatX_arg,
                        args_info.repeatY_arg, args_info.repeatZ_arg, 
                        args_info.periodicBox_given, args_info.dipole_given,
                        args_info.water_given);
     
  }

  delete xyzStream;
  delete zReader;
  delete dumpReader;
}

void replaceAtomType(SimInfo* info, vector<int>& zmolIndex, vector<double>& zmolPos, 
                                           vector<int> atomToZMol,double zconsTol){
      vector<int> zconsStatus;
      Molecule* mols;
      Atom** atoms;
      double com[3];
      size_t curIndex;
      char newAtomType[2000];
      int molIndex;

      mols = info->molecules;
      atoms = info->atoms;
      curIndex = 0;

      //determine the status of zconstraint molecules
      for(int i = 0; i< info->n_mol && curIndex < zmolIndex.size(); i++){

        if (mols[i].getGlobalIndex() == zmolIndex[curIndex]){
          mols[i].getCOM(com);
          if (fabs(com[2] - zmolPos[curIndex]) < zconsTol)
            zconsStatus.push_back(1);
          else
            zconsStatus.push_back(0);

          curIndex++;
        }
        
      }
      
      for(int j = 0; j < info->n_atoms; j++){
        molIndex = atomToZMol[j];
        if (molIndex >= 0){
          switch(zconsStatus[molIndex]){
            case 1 :
              strcpy(newAtomType, "ZF");
              break;
            case 0 :
              strcpy(newAtomType, "ZM");
              break;
            default:
              cerr << "error in zconsStatus" <<endl;
          }
              
           strcat(newAtomType, atoms[j]->getType());
           atoms[j]->setType(newAtomType);
        }
      }

}


void dumpToXYZ(SimInfo* info, ofstream* xyzStream, int repeatX, int repeatY, int repeatZ,
                                 int periodicBox_given, int dipole_given, int water_given) {
  char buffer[2000];
  int newN;
  Atom** atoms;
  double pos[3];
  char atomType[100];
  double HMat[3][3];
  int numNonSSD;
  int numSSD;

  numSSD = 0;
  numNonSSD = 0;

  atoms = info->atoms;

  //Determine the number of SSD molecules and Non-SSD molecules
  for (int i = 0; i < info->n_atoms; i++){
    atoms[i]->getType();
    if (strcmp(atomType, "SSD") || strcmp(atomType+2, "SSD"))
      numSSD++;
    else
      numNonSSD++;
  }
  
  if(water_given)
     newN = numNonSSD + numSSD * 4;
  else
    newN = numNonSSD;
  
  newN =  newN * (repeatX + 1) * (repeatY + 1) * (repeatZ + 1); 
  
  sprintf(buffer, "%d\n%lf;\t%lf\t%lf\t%lf;\t%lf\t%lf\t%lf;\t%lf\t%lf\t%lf;\n",
            newN, info->getTime(),
            info->Hmat[0][0], info->Hmat[1][0],info->Hmat[2][0],
            info->Hmat[0][1], info->Hmat[1][1],info->Hmat[2][1],
            info->Hmat[0][2], info->Hmat[1][2], info->Hmat[2][2] );
 (*xyzStream) << buffer;
 
  for(int i = 0; i < repeatX+1; i++){
    for(int j = 0; j < repeatY+1; j++){    
      for(int k = 0; k < repeatY+1; k++){    

         info->getBoxM(HMat);

        for (int l = 0; l < info->n_atoms; l++){
          atoms[l]->getPos(pos);
          strcpy(atomType, atoms[l]->getType());

          if (periodicBox_given)
              info->wrapVector(pos);
              
          for(int m = 0; m < 3; m++) {
            pos[m] += i * HMat[m][0] + j * HMat[m][1] + k * HMat[m][2];
          }

          if (atoms[l]->isDirectional()){     
            double RotMat[3][3];
            DirectionalAtom* dAtom;
            
            dAtom = (DirectionalAtom*) atoms[l];
            dAtom->getA(RotMat);
            
            if (strcmp(atomType, "SSD") == 0){

              if (water_given)
                convertSSD(buffer, "",pos, RotMat, dipole_given);
              else
                continue;
            }
            else if (strcmp(atomType+2, "SSD") == 0){
              if (water_given){
                //just a hard code here
                char prefix[3];
                prefix[0] = atomType[0];
                prefix[1] = atomType[1];
                prefix[2] = '\0';
                convertSSD(buffer, prefix,pos, RotMat, dipole_given);              
              }
              else
                continue;
            }
            
            else{
              if (dipole_given){
                double u[3] = {0, 0, 1};
                
                matVecMul3(RotMat, u);
                sprintf(buffer,"%s\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\n",
                          atomType, pos[0], pos [1], pos[2], u[0], u[1], u[2] );
              }
              else
                sprintf(buffer,"%s\t%lf\t%lf\t%lf\t0.0\t0.0\t0.0\n",
                          atomType,pos[0], pos [1], pos[2]);
            }
          }
          else{
              sprintf(buffer,"%s\t%lf\t%lf\t%lf\n",
                        atomType,pos[0], pos [1], pos[2]);
          }//end if atoms[i]->isDirectional()

          (*xyzStream) << buffer;
        }//end for int l = 0
        
      }//end for int k =0
    }// end for int j = 0
  }//end for int i = 0
  
}

void convertSSD(char* buffer, char* prefix, double pos[3], double RotMat[3][3], int dipole_given){
  double h1[3] = {0.0, -0.75695, 0.5206};
  double h2[3] = {0.0, 0.75695, 0.5206};
  double ox[3] = {0.0, 0.0, -0.0654};
  double u[3] = {0, 0, 1};

  matVecMul3(RotMat, u);
  matVecMul3(RotMat, h1);
  matVecMul3(RotMat, h2);
  matVecMul3(RotMat, ox);
  
  if (dipole_given){
    sprintf(buffer, "%sX\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\n",
              prefix,
              pos[0], pos[1], pos[2],
              u[0], u[1], u[2]);

    sprintf(buffer+strlen(buffer), "%sO\t%lf\t%lf\t%lf\t0.0\t0.0\t0.0\n",
              prefix,
              pos[0] + ox[0],
              pos[1] + ox[1],
              pos[2] + ox[2] );

    sprintf(buffer+strlen(buffer), "%sH\t%lf\t%lf\t%lf\t0.0\t0.0\t0.0\n",
              prefix,
              pos[0] + h1[0],
              pos[1] + h1[1],
              pos[2] + h1[2] );

    sprintf(buffer+strlen(buffer), "%sH\t%lf\t%lf\t%lf\t0.0\t0.0\t0.0\n",
              prefix,
              pos[0] + h2[0],
              pos[1] + h2[1],
              pos[2] + h2[2] );
  }
  else{
    sprintf(buffer, "%sX\t%lf\t%lf\t%lf\n",
              prefix,
              pos[0], pos[1], pos[2]);

    sprintf(buffer+strlen(buffer), "%sO\t%lf\t%lf\t%lf\n",
              prefix,
              pos[0] + ox[0],
              pos[1] + ox[1],
              pos[2] + ox[2] );

    sprintf(buffer+strlen(buffer), "%sH\t%lf\t%lf\t%lf\n",
              prefix,
              pos[0] + h1[0],
              pos[1] + h1[1],
              pos[2] + h1[2] );

    sprintf(buffer+strlen(buffer), "%sH\t%lf\t%lf\t%lf\n",
              prefix,
              pos[0] + h2[0],
              pos[1] + h2[1],
              pos[2] + h2[2] );
  }
}

vector<int> getAtomToZMol(SimInfo* info, vector<int>& zmolIndex){
  Molecule* mols;
  Atom** atoms;
  int nmol;
  int natom;
  vector<int> atomToZMol;
  int isZmol;
  size_t curIndex;
  
  atomToZMol.resize(info->n_atoms);
  mols = info->molecules;
  nmol = info->n_mol;
  curIndex = 0;

  for(int i = 0; i < nmol; i++){
    atoms =mols[i].getMyAtoms();

    if (curIndex >= zmolIndex.size()){
      isZmol = -1; 
    }
    else{
      if(mols[i].getGlobalIndex() == zmolIndex[curIndex]){
      isZmol = zmolIndex[curIndex];
      curIndex ++;
      }
      else
        isZmol = -1; 
    }

       
    natom = mols[i].getNAtoms();
    for(int j = 0; j < natom; j++)
#ifndef IS_MPI
      atomToZMol[atoms[j]->getIndex()] = isZmol;  
#else
      atomToZMol[atoms[j]->getGlobalIndex()] = isZmol;  
#endif 
  }

  return atomToZMol;
}

void matVecMul3(double a[3][3], double b[3]){
  double inVect[3];
  
  inVect[0] = b[0];
  inVect[1] = b[1];
  inVect[2] = b[2];

  b[0] = a[0][0]*inVect[0] + a[0][1]*inVect[1] + a[0][2]*inVect[2];
  b[1] = a[1][0]*inVect[0] + a[1][1]*inVect[1] + a[1][2]*inVect[2];
  b[2] = a[2][0]*inVect[0] + a[2][1]*inVect[1] + a[2][2]*inVect[2];
}
Revision:	1074
Committed:	Mon Mar 1 20:01:50 2004 UTC (21 years, 8 months ago) by tim
File size:	11678 byte(s)
Log Message:	Adding zsub, a program which can be used to replace atom type for zconstraint into OOPSE
#	User	Rev	Content
1	tim	1074	#include <cstdio>
2			#include <cmath>
3			#include <cstring>
4			#include <vector>
5
6			#include "CmdlineZsub.h"
7			#include "SimInfo.hpp"
8			#include "SimSetup.hpp"
9			#include "ZConsReader.hpp"
10
11			using namespace std;
12
13			/**
14			* @author: Teng Lin
15			* @date: Feb 25, 2004
16			* @version: 1.0
17			* @email: tlin@nd.edu
18			*
19			* This program is a substitue of quickLate. It not only implemnet all of the functions in
20			* quickLate, but also provide the functions of replacing atom types for z-constraint.
21			* This implement is kind of wierd!!!
22			*/
23			vector<int> getAtomToZMol(SimInfo* info, vector<int>& zmolIndex);
24			void replaceAtomType(SimInfo* info, vector<int>& zmolIndex, vector<double>& zmolPos,
25			vector<int> atomToZMol,double zconsTol);
26			void convertSSD(char* buffer, char* prefix, double pos[3], double RotMat[3][3], int dipole_given);
27			void dumpToXYZ(SimInfo* info, ofstream* xyzStream, int repeatX, int repeatY, int repeatZ,
28			int periodicBox_given, int dipole_given, int water_given);
29			void matVecMul3(double a[3][3], double b[3]);
30
31			int main(int argc, char* argv[]){
32			gengetopt_args_info args_info;
33			SimSetup startMe;
34			SimInfo* info;
35			string bassFileName;
36			string dumpFileName;
37			string xyzFileName;
38			int nFrame;
39			ZConsReader* zReader;
40			DumpReader* dumpReader;
41			ofstream* xyzStream;
42			int nZmol;
43			double zconsTol;
44			vector<int> zmolIndex;
45			vector<double> zmolPos;
46			vector<int> atomToZMol;
47			int havezcons;
48			//parse the command line option
49			if (cmdline_parser (argc, argv, &args_info) != 0)
50			exit(1) ;
51
52
53			//get the dumpfile name and bass file name
54			if (args_info.input_given){
55			dumpFileName = args_info.input_arg;
56			}
57			else{
58			cerr << "Does not have input file name" << endl;
59			exit(1);
60			}
61			bassFileName = dumpFileName;
62			bassFileName = bassFileName.substr(0, bassFileName.rfind(".")) + ".bass";
63
64			if (args_info.output_given){
65			xyzFileName = args_info.output_arg;
66			}
67			else{
68			xyzFileName = dumpFileName;
69			xyzFileName = xyzFileName.substr(0, xyzFileName.rfind(".")) + ".xyz";
70			}
71
72			info = new SimInfo();
73			startMe.setSimInfo(info );
74
75			startMe.parseFile( bassFileName.c_str());
76
77			startMe.createSim();
78
79			dumpReader = new DumpReader(dumpFileName.c_str());
80
81			nFrame = dumpReader->getNframes();
82
83			if (args_info.replace_given && info->nZconstraints > 0){
84
85			zReader = new ZConsReader(info);
86			nZmol = zReader->getNumZMol();
87			zmolIndex = zReader->getZConsIndex();
88			zmolPos = zReader->getZConsPos();
89
90			atomToZMol = getAtomToZMol(info, zmolIndex);
91
92			//get tolerance
93			GenericData* data;
94			DoubleData* toleranceData;
95			data = info->getProperty(ZCONSTOL_ID);
96
97			if(!data) {
98			cerr << "Zsub error: can not get tolerance \n" << endl;
99			exit(1);
100			}
101			else{
102
103			toleranceData = dynamic_cast<DoubleData*>(data);
104
105			if(!toleranceData){
106			cerr << "Zsub error: can not get tolerance \n" << endl;
107			exit(1);
108			}
109			else{
110			zconsTol = toleranceData->getData();
111			}
112			}
113
114
115
116			}
117
118			xyzStream = new ofstream(xyzFileName.c_str());
119			for (int i = 0; i < nFrame; i += args_info.frame_arg){
120			dumpReader->readFrame(info, i);
121			if (args_info.replace_given && info->nZconstraints > 0){
122			replaceAtomType(info, zmolIndex, zmolPos, atomToZMol, zconsTol);
123			}
124			dumpToXYZ(info, xyzStream,args_info.repeatX_arg,
125			args_info.repeatY_arg, args_info.repeatZ_arg,
126			args_info.periodicBox_given, args_info.dipole_given,
127			args_info.water_given);
128
129			}
130
131			delete xyzStream;
132			delete zReader;
133			delete dumpReader;
134			}
135
136			void replaceAtomType(SimInfo* info, vector<int>& zmolIndex, vector<double>& zmolPos,
137			vector<int> atomToZMol,double zconsTol){
138			vector<int> zconsStatus;
139			Molecule* mols;
140			Atom** atoms;
141			double com[3];
142			size_t curIndex;
143			char newAtomType[2000];
144			int molIndex;
145
146			mols = info->molecules;
147			atoms = info->atoms;
148			curIndex = 0;
149
150			//determine the status of zconstraint molecules
151			for(int i = 0; i< info->n_mol && curIndex < zmolIndex.size(); i++){
152
153			if (mols[i].getGlobalIndex() == zmolIndex[curIndex]){
154			mols[i].getCOM(com);
155			if (fabs(com[2] - zmolPos[curIndex]) < zconsTol)
156			zconsStatus.push_back(1);
157			else
158			zconsStatus.push_back(0);
159
160			curIndex++;
161			}
162
163			}
164
165			for(int j = 0; j < info->n_atoms; j++){
166			molIndex = atomToZMol[j];
167			if (molIndex >= 0){
168			switch(zconsStatus[molIndex]){
169			case 1 :
170			strcpy(newAtomType, "ZF");
171			break;
172			case 0 :
173			strcpy(newAtomType, "ZM");
174			break;
175			default:
176			cerr << "error in zconsStatus" <<endl;
177			}
178
179			strcat(newAtomType, atoms[j]->getType());
180			atoms[j]->setType(newAtomType);
181			}
182			}
183
184			}
185
186
187			void dumpToXYZ(SimInfo* info, ofstream* xyzStream, int repeatX, int repeatY, int repeatZ,
188			int periodicBox_given, int dipole_given, int water_given) {
189			char buffer[2000];
190			int newN;
191			Atom** atoms;
192			double pos[3];
193			char atomType[100];
194			double HMat[3][3];
195			int numNonSSD;
196			int numSSD;
197
198			numSSD = 0;
199			numNonSSD = 0;
200
201			atoms = info->atoms;
202
203			//Determine the number of SSD molecules and Non-SSD molecules
204			for (int i = 0; i < info->n_atoms; i++){
205			atoms[i]->getType();
206			if (strcmp(atomType, "SSD") \|\| strcmp(atomType+2, "SSD"))
207			numSSD++;
208			else
209			numNonSSD++;
210			}
211
212			if(water_given)
213			newN = numNonSSD + numSSD * 4;
214			else
215			newN = numNonSSD;
216
217			newN = newN * (repeatX + 1) * (repeatY + 1) * (repeatZ + 1);
218
219			sprintf(buffer, "%d\n%lf;\t%lf\t%lf\t%lf;\t%lf\t%lf\t%lf;\t%lf\t%lf\t%lf;\n",
220			newN, info->getTime(),
221			info->Hmat[0][0], info->Hmat[1][0],info->Hmat[2][0],
222			info->Hmat[0][1], info->Hmat[1][1],info->Hmat[2][1],
223			info->Hmat[0][2], info->Hmat[1][2], info->Hmat[2][2] );
224			(*xyzStream) << buffer;
225
226			for(int i = 0; i < repeatX+1; i++){
227			for(int j = 0; j < repeatY+1; j++){
228			for(int k = 0; k < repeatY+1; k++){
229
230			info->getBoxM(HMat);
231
232			for (int l = 0; l < info->n_atoms; l++){
233			atoms[l]->getPos(pos);
234			strcpy(atomType, atoms[l]->getType());
235
236			if (periodicBox_given)
237			info->wrapVector(pos);
238
239			for(int m = 0; m < 3; m++) {
240			pos[m] += i * HMat[m][0] + j * HMat[m][1] + k * HMat[m][2];
241			}
242
243			if (atoms[l]->isDirectional()){
244			double RotMat[3][3];
245			DirectionalAtom* dAtom;
246
247			dAtom = (DirectionalAtom*) atoms[l];
248			dAtom->getA(RotMat);
249
250			if (strcmp(atomType, "SSD") == 0){
251
252			if (water_given)
253			convertSSD(buffer, "",pos, RotMat, dipole_given);
254			else
255			continue;
256			}
257			else if (strcmp(atomType+2, "SSD") == 0){
258			if (water_given){
259			//just a hard code here
260			char prefix[3];
261			prefix[0] = atomType[0];
262			prefix[1] = atomType[1];
263			prefix[2] = '\0';
264			convertSSD(buffer, prefix,pos, RotMat, dipole_given);
265			}
266			else
267			continue;
268			}
269
270			else{
271			if (dipole_given){
272			double u[3] = {0, 0, 1};
273
274			matVecMul3(RotMat, u);
275			sprintf(buffer,"%s\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\n",
276			atomType, pos[0], pos [1], pos[2], u[0], u[1], u[2] );
277			}
278			else
279			sprintf(buffer,"%s\t%lf\t%lf\t%lf\t0.0\t0.0\t0.0\n",
280			atomType,pos[0], pos [1], pos[2]);
281			}
282			}
283			else{
284			sprintf(buffer,"%s\t%lf\t%lf\t%lf\n",
285			atomType,pos[0], pos [1], pos[2]);
286			}//end if atoms[i]->isDirectional()
287
288			(*xyzStream) << buffer;
289			}//end for int l = 0
290
291			}//end for int k =0
292			}// end for int j = 0
293			}//end for int i = 0
294
295			}
296
297			void convertSSD(char* buffer, char* prefix, double pos[3], double RotMat[3][3], int dipole_given){
298			double h1[3] = {0.0, -0.75695, 0.5206};
299			double h2[3] = {0.0, 0.75695, 0.5206};
300			double ox[3] = {0.0, 0.0, -0.0654};
301			double u[3] = {0, 0, 1};
302
303			matVecMul3(RotMat, u);
304			matVecMul3(RotMat, h1);
305			matVecMul3(RotMat, h2);
306			matVecMul3(RotMat, ox);
307
308			if (dipole_given){
309			sprintf(buffer, "%sX\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\n",
310			prefix,
311			pos[0], pos[1], pos[2],
312			u[0], u[1], u[2]);
313
314			sprintf(buffer+strlen(buffer), "%sO\t%lf\t%lf\t%lf\t0.0\t0.0\t0.0\n",
315			prefix,
316			pos[0] + ox[0],
317			pos[1] + ox[1],
318			pos[2] + ox[2] );
319
320			sprintf(buffer+strlen(buffer), "%sH\t%lf\t%lf\t%lf\t0.0\t0.0\t0.0\n",
321			prefix,
322			pos[0] + h1[0],
323			pos[1] + h1[1],
324			pos[2] + h1[2] );
325
326			sprintf(buffer+strlen(buffer), "%sH\t%lf\t%lf\t%lf\t0.0\t0.0\t0.0\n",
327			prefix,
328			pos[0] + h2[0],
329			pos[1] + h2[1],
330			pos[2] + h2[2] );
331			}
332			else{
333			sprintf(buffer, "%sX\t%lf\t%lf\t%lf\n",
334			prefix,
335			pos[0], pos[1], pos[2]);
336
337			sprintf(buffer+strlen(buffer), "%sO\t%lf\t%lf\t%lf\n",
338			prefix,
339			pos[0] + ox[0],
340			pos[1] + ox[1],
341			pos[2] + ox[2] );
342
343			sprintf(buffer+strlen(buffer), "%sH\t%lf\t%lf\t%lf\n",
344			prefix,
345			pos[0] + h1[0],
346			pos[1] + h1[1],
347			pos[2] + h1[2] );
348
349			sprintf(buffer+strlen(buffer), "%sH\t%lf\t%lf\t%lf\n",
350			prefix,
351			pos[0] + h2[0],
352			pos[1] + h2[1],
353			pos[2] + h2[2] );
354			}
355			}
356
357			vector<int> getAtomToZMol(SimInfo* info, vector<int>& zmolIndex){
358			Molecule* mols;
359			Atom** atoms;
360			int nmol;
361			int natom;
362			vector<int> atomToZMol;
363			int isZmol;
364			size_t curIndex;
365
366			atomToZMol.resize(info->n_atoms);
367			mols = info->molecules;
368			nmol = info->n_mol;
369			curIndex = 0;
370
371			for(int i = 0; i < nmol; i++){
372			atoms =mols[i].getMyAtoms();
373
374			if (curIndex >= zmolIndex.size()){
375			isZmol = -1;
376			}
377			else{
378			if(mols[i].getGlobalIndex() == zmolIndex[curIndex]){
379			isZmol = zmolIndex[curIndex];
380			curIndex ++;
381			}
382			else
383			isZmol = -1;
384			}
385
386
387
388			natom = mols[i].getNAtoms();
389			for(int j = 0; j < natom; j++)
390			#ifndef IS_MPI
391			atomToZMol[atoms[j]->getIndex()] = isZmol;
392			#else
393			atomToZMol[atoms[j]->getGlobalIndex()] = isZmol;
394			#endif
395			}
396
397			return atomToZMol;
398			}
399
400			void matVecMul3(double a[3][3], double b[3]){
401			double inVect[3];
402
403			inVect[0] = b[0];
404			inVect[1] = b[1];
405			inVect[2] = b[2];
406
407			b[0] = a[0][0]inVect[0] + a[0][1]inVect[1] + a[0][2]*inVect[2];
408			b[1] = a[1][0]inVect[0] + a[1][1]inVect[1] + a[1][2]*inVect[2];
409			b[2] = a[2][0]inVect[0] + a[2][1]inVect[1] + a[2][2]*inVect[2];
410			}