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 implemnets all of the functions in
 * quickLate, but also provides 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();
      }

      if (zReader->hasNextFrame())
        zReader->readNextFrame();

    }    


  }


  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){    

      while (zReader->getCurTime() < info->currentTime){
        // (zReader->hasNextFrame())
          zReader->readNextFrame();
      }

      zmolPos = zReader->getZConsPos();

      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;
  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++){
    atomType = 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);
          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\n",
                          atomType,pos[0], pos [1], pos[2]);
            }
          }
          else{
              if (dipole_given)
                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 = 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:	1116
Committed:	Thu Apr 15 22:15:21 2004 UTC (20 years, 2 months ago) by tim
File size:	12117 byte(s)
Log Message:	fix a bug in setting exclude list
#	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	tim	1116	* This program is a substitue of quickLate. It not only implemnets all of the functions in
20			* quickLate, but also provides the functions of replacing atom types for z-constraint.
21	tim	1074	* 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	tim	1108
92	tim	1074	//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	tim	1108
113			if (zReader->hasNextFrame())
114			zReader->readNextFrame();
115
116	tim	1074	}
117
118
119
120			}
121
122	tim	1108
123	tim	1074	xyzStream = new ofstream(xyzFileName.c_str());
124			for (int i = 0; i < nFrame; i += args_info.frame_arg){
125			dumpReader->readFrame(info, i);
126	tim	1108
127	tim	1074	if (args_info.replace_given && info->nZconstraints > 0){
128	tim	1108
129			while (zReader->getCurTime() < info->currentTime){
130			// (zReader->hasNextFrame())
131			zReader->readNextFrame();
132			}
133
134			zmolPos = zReader->getZConsPos();
135
136	tim	1074	replaceAtomType(info, zmolIndex, zmolPos, atomToZMol, zconsTol);
137			}
138			dumpToXYZ(info, xyzStream,args_info.repeatX_arg,
139			args_info.repeatY_arg, args_info.repeatZ_arg,
140			args_info.periodicBox_given, args_info.dipole_given,
141			args_info.water_given);
142
143			}
144
145			delete xyzStream;
146			delete zReader;
147			delete dumpReader;
148			}
149
150			void replaceAtomType(SimInfo* info, vector<int>& zmolIndex, vector<double>& zmolPos,
151			vector<int> atomToZMol,double zconsTol){
152			vector<int> zconsStatus;
153			Molecule* mols;
154			Atom** atoms;
155			double com[3];
156			size_t curIndex;
157			char newAtomType[2000];
158			int molIndex;
159
160			mols = info->molecules;
161			atoms = info->atoms;
162			curIndex = 0;
163
164			//determine the status of zconstraint molecules
165			for(int i = 0; i< info->n_mol && curIndex < zmolIndex.size(); i++){
166
167			if (mols[i].getGlobalIndex() == zmolIndex[curIndex]){
168			mols[i].getCOM(com);
169			if (fabs(com[2] - zmolPos[curIndex]) < zconsTol)
170			zconsStatus.push_back(1);
171			else
172			zconsStatus.push_back(0);
173
174			curIndex++;
175			}
176
177			}
178
179			for(int j = 0; j < info->n_atoms; j++){
180			molIndex = atomToZMol[j];
181			if (molIndex >= 0){
182			switch(zconsStatus[molIndex]){
183			case 1 :
184			strcpy(newAtomType, "ZF");
185			break;
186			case 0 :
187			strcpy(newAtomType, "ZM");
188			break;
189			default:
190			cerr << "error in zconsStatus" <<endl;
191			}
192
193			strcat(newAtomType, atoms[j]->getType());
194			atoms[j]->setType(newAtomType);
195			}
196			}
197
198			}
199
200
201			void dumpToXYZ(SimInfo* info, ofstream* xyzStream, int repeatX, int repeatY, int repeatZ,
202			int periodicBox_given, int dipole_given, int water_given) {
203			char buffer[2000];
204			int newN;
205			Atom** atoms;
206			double pos[3];
207	tim	1075	char* atomType;
208	tim	1074	double HMat[3][3];
209			int numNonSSD;
210			int numSSD;
211
212			numSSD = 0;
213			numNonSSD = 0;
214
215			atoms = info->atoms;
216
217			//Determine the number of SSD molecules and Non-SSD molecules
218			for (int i = 0; i < info->n_atoms; i++){
219	tim	1075	atomType = atoms[i]->getType();
220			if (!strcmp(atomType, "SSD") \|\| !strcmp(atomType+2, "SSD"))
221	tim	1074	numSSD++;
222			else
223			numNonSSD++;
224			}
225
226			if(water_given)
227			newN = numNonSSD + numSSD * 4;
228			else
229			newN = numNonSSD;
230
231			newN = newN * (repeatX + 1) * (repeatY + 1) * (repeatZ + 1);
232
233			sprintf(buffer, "%d\n%lf;\t%lf\t%lf\t%lf;\t%lf\t%lf\t%lf;\t%lf\t%lf\t%lf;\n",
234			newN, info->getTime(),
235			info->Hmat[0][0], info->Hmat[1][0],info->Hmat[2][0],
236			info->Hmat[0][1], info->Hmat[1][1],info->Hmat[2][1],
237			info->Hmat[0][2], info->Hmat[1][2], info->Hmat[2][2] );
238			(*xyzStream) << buffer;
239
240			for(int i = 0; i < repeatX+1; i++){
241			for(int j = 0; j < repeatY+1; j++){
242			for(int k = 0; k < repeatY+1; k++){
243
244			info->getBoxM(HMat);
245
246			for (int l = 0; l < info->n_atoms; l++){
247			atoms[l]->getPos(pos);
248	tim	1075	atomType = atoms[l]->getType();
249	tim	1074
250			if (periodicBox_given)
251			info->wrapVector(pos);
252
253			for(int m = 0; m < 3; m++) {
254			pos[m] += i * HMat[m][0] + j * HMat[m][1] + k * HMat[m][2];
255			}
256
257			if (atoms[l]->isDirectional()){
258			double RotMat[3][3];
259			DirectionalAtom* dAtom;
260
261			dAtom = (DirectionalAtom*) atoms[l];
262			dAtom->getA(RotMat);
263
264			if (strcmp(atomType, "SSD") == 0){
265
266			if (water_given)
267			convertSSD(buffer, "",pos, RotMat, dipole_given);
268			else
269			continue;
270			}
271			else if (strcmp(atomType+2, "SSD") == 0){
272			if (water_given){
273			//just a hard code here
274			char prefix[3];
275			prefix[0] = atomType[0];
276			prefix[1] = atomType[1];
277			prefix[2] = '\0';
278			convertSSD(buffer, prefix,pos, RotMat, dipole_given);
279			}
280			else
281			continue;
282			}
283
284			else{
285			if (dipole_given){
286			double u[3] = {0, 0, 1};
287
288			matVecMul3(RotMat, u);
289			sprintf(buffer,"%s\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\n",
290			atomType, pos[0], pos [1], pos[2], u[0], u[1], u[2] );
291			}
292			else
293	tim	1075	sprintf(buffer,"%s\t%lf\t%lf\t%lf\n",
294	tim	1074	atomType,pos[0], pos [1], pos[2]);
295			}
296			}
297			else{
298	tim	1075	if (dipole_given)
299			sprintf(buffer,"%s\t%lf\t%lf\t%lf\t0.0\t0.0\t0.0\n",
300	tim	1074	atomType,pos[0], pos [1], pos[2]);
301	tim	1075
302			else
303			sprintf(buffer,"%s\t%lf\t%lf\t%lf\n",
304			atomType,pos[0], pos [1], pos[2]);
305	tim	1074	}//end if atoms[i]->isDirectional()
306
307			(*xyzStream) << buffer;
308			}//end for int l = 0
309
310			}//end for int k =0
311			}// end for int j = 0
312			}//end for int i = 0
313
314			}
315
316			void convertSSD(char* buffer, char* prefix, double pos[3], double RotMat[3][3], int dipole_given){
317			double h1[3] = {0.0, -0.75695, 0.5206};
318			double h2[3] = {0.0, 0.75695, 0.5206};
319			double ox[3] = {0.0, 0.0, -0.0654};
320			double u[3] = {0, 0, 1};
321
322			matVecMul3(RotMat, u);
323			matVecMul3(RotMat, h1);
324			matVecMul3(RotMat, h2);
325			matVecMul3(RotMat, ox);
326
327			if (dipole_given){
328			sprintf(buffer, "%sX\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\n",
329			prefix,
330			pos[0], pos[1], pos[2],
331			u[0], u[1], u[2]);
332
333			sprintf(buffer+strlen(buffer), "%sO\t%lf\t%lf\t%lf\t0.0\t0.0\t0.0\n",
334			prefix,
335			pos[0] + ox[0],
336			pos[1] + ox[1],
337			pos[2] + ox[2] );
338
339			sprintf(buffer+strlen(buffer), "%sH\t%lf\t%lf\t%lf\t0.0\t0.0\t0.0\n",
340			prefix,
341			pos[0] + h1[0],
342			pos[1] + h1[1],
343			pos[2] + h1[2] );
344
345			sprintf(buffer+strlen(buffer), "%sH\t%lf\t%lf\t%lf\t0.0\t0.0\t0.0\n",
346			prefix,
347			pos[0] + h2[0],
348			pos[1] + h2[1],
349			pos[2] + h2[2] );
350			}
351			else{
352			sprintf(buffer, "%sX\t%lf\t%lf\t%lf\n",
353			prefix,
354			pos[0], pos[1], pos[2]);
355
356			sprintf(buffer+strlen(buffer), "%sO\t%lf\t%lf\t%lf\n",
357			prefix,
358			pos[0] + ox[0],
359			pos[1] + ox[1],
360			pos[2] + ox[2] );
361
362			sprintf(buffer+strlen(buffer), "%sH\t%lf\t%lf\t%lf\n",
363			prefix,
364			pos[0] + h1[0],
365			pos[1] + h1[1],
366			pos[2] + h1[2] );
367
368			sprintf(buffer+strlen(buffer), "%sH\t%lf\t%lf\t%lf\n",
369			prefix,
370			pos[0] + h2[0],
371			pos[1] + h2[1],
372			pos[2] + h2[2] );
373			}
374			}
375
376			vector<int> getAtomToZMol(SimInfo* info, vector<int>& zmolIndex){
377			Molecule* mols;
378			Atom** atoms;
379			int nmol;
380			int natom;
381			vector<int> atomToZMol;
382			int isZmol;
383			size_t curIndex;
384
385			atomToZMol.resize(info->n_atoms);
386			mols = info->molecules;
387			nmol = info->n_mol;
388			curIndex = 0;
389
390			for(int i = 0; i < nmol; i++){
391			atoms =mols[i].getMyAtoms();
392
393			if (curIndex >= zmolIndex.size()){
394			isZmol = -1;
395			}
396			else{
397			if(mols[i].getGlobalIndex() == zmolIndex[curIndex]){
398	tim	1075	isZmol = curIndex;
399	tim	1074	curIndex ++;
400			}
401			else
402			isZmol = -1;
403			}
404
405
406
407			natom = mols[i].getNAtoms();
408			for(int j = 0; j < natom; j++)
409			#ifndef IS_MPI
410			atomToZMol[atoms[j]->getIndex()] = isZmol;
411			#else
412			atomToZMol[atoms[j]->getGlobalIndex()] = isZmol;
413			#endif
414			}
415
416			return atomToZMol;
417			}
418
419			void matVecMul3(double a[3][3], double b[3]){
420			double inVect[3];
421
422			inVect[0] = b[0];
423			inVect[1] = b[1];
424			inVect[2] = b[2];
425
426			b[0] = a[0][0]inVect[0] + a[0][1]inVect[1] + a[0][2]*inVect[2];
427			b[1] = a[1][0]inVect[0] + a[1][1]inVect[1] + a[1][2]*inVect[2];
428			b[2] = a[2][0]inVect[0] + a[2][1]inVect[1] + a[2][2]*inVect[2];
429			}