| 1 |
#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 |
if (zReader->hasNextFrame())
|
| 114 |
zReader->readNextFrame();
|
| 115 |
|
| 116 |
}
|
| 117 |
|
| 118 |
|
| 119 |
|
| 120 |
}
|
| 121 |
|
| 122 |
|
| 123 |
xyzStream = new ofstream(xyzFileName.c_str());
|
| 124 |
for (int i = 0; i < nFrame; i += args_info.frame_arg){
|
| 125 |
dumpReader->readFrame(info, i);
|
| 126 |
|
| 127 |
if (args_info.replace_given && info->nZconstraints > 0){
|
| 128 |
|
| 129 |
while (zReader->getCurTime() < info->currentTime){
|
| 130 |
// (zReader->hasNextFrame())
|
| 131 |
zReader->readNextFrame();
|
| 132 |
}
|
| 133 |
|
| 134 |
zmolPos = zReader->getZConsPos();
|
| 135 |
|
| 136 |
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 |
char* atomType;
|
| 208 |
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 |
atomType = atoms[i]->getType();
|
| 220 |
if (!strcmp(atomType, "SSD") || !strcmp(atomType+2, "SSD"))
|
| 221 |
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 |
atomType = atoms[l]->getType();
|
| 249 |
|
| 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 |
sprintf(buffer,"%s\t%lf\t%lf\t%lf\n",
|
| 294 |
atomType,pos[0], pos [1], pos[2]);
|
| 295 |
}
|
| 296 |
}
|
| 297 |
else{
|
| 298 |
if (dipole_given)
|
| 299 |
sprintf(buffer,"%s\t%lf\t%lf\t%lf\t0.0\t0.0\t0.0\n",
|
| 300 |
atomType,pos[0], pos [1], pos[2]);
|
| 301 |
|
| 302 |
else
|
| 303 |
sprintf(buffer,"%s\t%lf\t%lf\t%lf\n",
|
| 304 |
atomType,pos[0], pos [1], pos[2]);
|
| 305 |
}//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 |
isZmol = curIndex;
|
| 399 |
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 |
}
|
