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 |
|
|
}
|