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;
  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:	1075
Committed:	Mon Mar 1 21:17:27 2004 UTC (20 years, 4 months ago) by tim
File size:	11841 byte(s)
Log Message:	Fix a couple of bugs in zsub
#	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	tim	1075	char* atomType;
194	tim	1074	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	tim	1075	atomType = atoms[i]->getType();
206			if (!strcmp(atomType, "SSD") \|\| !strcmp(atomType+2, "SSD"))
207	tim	1074	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	tim	1075	atomType = atoms[l]->getType();
235	tim	1074
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	tim	1075	sprintf(buffer,"%s\t%lf\t%lf\t%lf\n",
280	tim	1074	atomType,pos[0], pos [1], pos[2]);
281			}
282			}
283			else{
284	tim	1075	if (dipole_given)
285			sprintf(buffer,"%s\t%lf\t%lf\t%lf\t0.0\t0.0\t0.0\n",
286	tim	1074	atomType,pos[0], pos [1], pos[2]);
287	tim	1075
288			else
289			sprintf(buffer,"%s\t%lf\t%lf\t%lf\n",
290			atomType,pos[0], pos [1], pos[2]);
291	tim	1074	}//end if atoms[i]->isDirectional()
292
293			(*xyzStream) << buffer;
294			}//end for int l = 0
295
296			}//end for int k =0
297			}// end for int j = 0
298			}//end for int i = 0
299
300			}
301
302			void convertSSD(char* buffer, char* prefix, double pos[3], double RotMat[3][3], int dipole_given){
303			double h1[3] = {0.0, -0.75695, 0.5206};
304			double h2[3] = {0.0, 0.75695, 0.5206};
305			double ox[3] = {0.0, 0.0, -0.0654};
306			double u[3] = {0, 0, 1};
307
308			matVecMul3(RotMat, u);
309			matVecMul3(RotMat, h1);
310			matVecMul3(RotMat, h2);
311			matVecMul3(RotMat, ox);
312
313			if (dipole_given){
314			sprintf(buffer, "%sX\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\n",
315			prefix,
316			pos[0], pos[1], pos[2],
317			u[0], u[1], u[2]);
318
319			sprintf(buffer+strlen(buffer), "%sO\t%lf\t%lf\t%lf\t0.0\t0.0\t0.0\n",
320			prefix,
321			pos[0] + ox[0],
322			pos[1] + ox[1],
323			pos[2] + ox[2] );
324
325			sprintf(buffer+strlen(buffer), "%sH\t%lf\t%lf\t%lf\t0.0\t0.0\t0.0\n",
326			prefix,
327			pos[0] + h1[0],
328			pos[1] + h1[1],
329			pos[2] + h1[2] );
330
331			sprintf(buffer+strlen(buffer), "%sH\t%lf\t%lf\t%lf\t0.0\t0.0\t0.0\n",
332			prefix,
333			pos[0] + h2[0],
334			pos[1] + h2[1],
335			pos[2] + h2[2] );
336			}
337			else{
338			sprintf(buffer, "%sX\t%lf\t%lf\t%lf\n",
339			prefix,
340			pos[0], pos[1], pos[2]);
341
342			sprintf(buffer+strlen(buffer), "%sO\t%lf\t%lf\t%lf\n",
343			prefix,
344			pos[0] + ox[0],
345			pos[1] + ox[1],
346			pos[2] + ox[2] );
347
348			sprintf(buffer+strlen(buffer), "%sH\t%lf\t%lf\t%lf\n",
349			prefix,
350			pos[0] + h1[0],
351			pos[1] + h1[1],
352			pos[2] + h1[2] );
353
354			sprintf(buffer+strlen(buffer), "%sH\t%lf\t%lf\t%lf\n",
355			prefix,
356			pos[0] + h2[0],
357			pos[1] + h2[1],
358			pos[2] + h2[2] );
359			}
360			}
361
362			vector<int> getAtomToZMol(SimInfo* info, vector<int>& zmolIndex){
363			Molecule* mols;
364			Atom** atoms;
365			int nmol;
366			int natom;
367			vector<int> atomToZMol;
368			int isZmol;
369			size_t curIndex;
370
371			atomToZMol.resize(info->n_atoms);
372			mols = info->molecules;
373			nmol = info->n_mol;
374			curIndex = 0;
375
376			for(int i = 0; i < nmol; i++){
377			atoms =mols[i].getMyAtoms();
378
379			if (curIndex >= zmolIndex.size()){
380			isZmol = -1;
381			}
382			else{
383			if(mols[i].getGlobalIndex() == zmolIndex[curIndex]){
384	tim	1075	isZmol = curIndex;
385	tim	1074	curIndex ++;
386			}
387			else
388			isZmol = -1;
389			}
390
391
392
393			natom = mols[i].getNAtoms();
394			for(int j = 0; j < natom; j++)
395			#ifndef IS_MPI
396			atomToZMol[atoms[j]->getIndex()] = isZmol;
397			#else
398			atomToZMol[atoms[j]->getGlobalIndex()] = isZmol;
399			#endif
400			}
401
402			return atomToZMol;
403			}
404
405			void matVecMul3(double a[3][3], double b[3]){
406			double inVect[3];
407
408			inVect[0] = b[0];
409			inVect[1] = b[1];
410			inVect[2] = b[2];
411
412			b[0] = a[0][0]inVect[0] + a[0][1]inVect[1] + a[0][2]*inVect[2];
413			b[1] = a[1][0]inVect[0] + a[1][1]inVect[1] + a[1][2]*inVect[2];
414			b[2] = a[2][0]inVect[0] + a[2][1]inVect[1] + a[2][2]*inVect[2];
415			}