1 /********************************************************************** 2 Copyright (C) 2011 by Albert DeFusco 3 4 This program is free software; you can redistribute it and/or modify 5 it under the terms of the GNU General Public License as published by 6 the Free Software Foundation version 2 of the License. 7 8 This program is distributed in the hope that it will be useful, 9 but WITHOUT ANY WARRANTY; without even the implied warranty of 10 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 11 GNU General Public License for more details. 12 ***********************************************************************/ 13 14 #include <openbabel/babelconfig.h> 15 #include <openbabel/obmolecformat.h> 16 #include <openbabel/mol.h> 17 #include <openbabel/atom.h> 18 #include <openbabel/elements.h> 19 #include <openbabel/bond.h> 20 21 #include <openbabel/obiter.h> 22 23 #include <sstream> 24 #include <map> 25 #include <cstdlib> 26 27 using namespace std; 28 namespace OpenBabel 29 { 30 31 class LMPDATFormat : public OBMoleculeFormat 32 { 33 public: 34 //Register this format type ID LMPDATFormat()35 LMPDATFormat() 36 { 37 OBConversion::RegisterFormat("lmpdat", this, "chemical/x-lmpdat"); 38 OBConversion::RegisterOptionParam("q", nullptr, 1, OBConversion::OUTOPTIONS); 39 OBConversion::RegisterOptionParam("d", nullptr, 1, OBConversion::OUTOPTIONS); 40 } 41 Description()42 virtual const char* Description() //required 43 { 44 return 45 "The LAMMPS data format\n\n" 46 47 "LAMMPS is a classical molecular dynamics code, and an acronym for\n" 48 "Large-scale Atomic/Molecular Massively Parallel Simulator.\n\n" 49 50 "Write Options, e.g. -xq\n" 51 " q \"water-model\" Set atomic charges for water.\n" 52 " There are two options: SPC (default) or SPCE\n" 53 " d \"length\" Set the length of the boundary box around the molecule.\n" 54 " The default is to make a cube around the molecule\n" 55 " adding 50% to the most positive and negative\n" 56 " cartesian coordinate.\n" 57 ; 58 }; 59 60 GetMIMEType()61 virtual const char* GetMIMEType() 62 { return "chemical/x-lmpdat"; }; 63 Flags()64 virtual unsigned int Flags() 65 { 66 return NOTREADABLE | WRITEONEONLY; 67 }; 68 69 virtual bool WriteMolecule(OBBase* pOb, OBConversion* pConv); 70 }; 71 //*** 72 73 //Make an instance of the format class 74 LMPDATFormat theLMPDATFormat; 75 76 //////////////////////////////////////////////////////////////// 77 WriteMolecule(OBBase * pOb,OBConversion * pConv)78 bool LMPDATFormat::WriteMolecule(OBBase* pOb, OBConversion* pConv) 79 { 80 OBMol* pmol = dynamic_cast<OBMol*>(pOb); 81 if (pmol == nullptr) 82 return false; 83 84 //Define some references so we can use the old parameter names 85 ostream &ofs = *pConv->GetOutStream(); 86 OBMol &mol = *pmol; 87 OBAtom *a; 88 OBAtom *b; 89 OBAtom *c; 90 OBAtom *d; 91 92 char buffer[BUFF_SIZE]; 93 94 95 //Very basic atom typing by element only 96 //TODO: The maps should become smarts strings instead of element names 97 double ThisMass; 98 string ThisAtom; 99 map<string, double> AtomMass; 100 FOR_ATOMS_OF_MOL(atom, mol) 101 { 102 ThisMass=atom->GetAtomicMass(); 103 ThisAtom=OBElements::GetSymbol(atom->GetAtomicNum()); 104 AtomMass[ThisAtom] = ThisMass; 105 } 106 map<string, int> AtomType; 107 int AtomIndex=0; 108 map<string, double>::iterator itr; 109 //Set AtomType integer 110 for(itr=AtomMass.begin();itr!=AtomMass.end();++itr) 111 { 112 AtomIndex++; 113 AtomType[itr->first] = AtomIndex; 114 } 115 116 //Determine unique bonds 117 mol.ConnectTheDots(); 118 char BondString[BUFF_SIZE]; 119 map<string, int> BondType; 120 FOR_BONDS_OF_MOL(bond,mol) 121 { 122 a=bond->GetBeginAtom(); 123 b=bond->GetEndAtom(); 124 125 //To let the unordered map determine unique keys, 126 //always put the heavier atom first 127 if(b->GetAtomicNum() > a->GetAtomicNum() ) 128 { 129 snprintf(BondString,BUFF_SIZE, 130 "%2s:%2s", 131 OBElements::GetSymbol(b->GetAtomicNum()), 132 OBElements::GetSymbol(a->GetAtomicNum())); 133 } 134 else 135 { 136 snprintf(BondString,BUFF_SIZE, 137 "%2s:%2s", 138 OBElements::GetSymbol(a->GetAtomicNum()), 139 OBElements::GetSymbol(b->GetAtomicNum())); 140 } 141 BondType[BondString] = 0; 142 } 143 map<string, int>::iterator intitr; 144 int BondIndex=0; 145 //Set the BondType integer 146 for(intitr=BondType.begin();intitr!=BondType.end();++intitr) 147 { 148 BondIndex++; 149 BondType[intitr->first] = BondIndex; 150 } 151 152 //Determine unique angles 153 mol.FindAngles(); 154 int anglecount=0; 155 char AngleString[BUFF_SIZE]; 156 map<string, int> AngleType; 157 FOR_ANGLES_OF_MOL(angle,mol) 158 { 159 anglecount++; 160 a = mol.GetAtom((*angle)[0]+1); 161 b = mol.GetAtom((*angle)[1]+1); 162 c = mol.GetAtom((*angle)[2]+1); 163 //Origanization trick: 164 //1. The atom "a" is acutally the middle 165 //2. Always write the heavy element first 166 if(b->GetAtomicNum() > c->GetAtomicNum() ) 167 { 168 snprintf(AngleString,BUFF_SIZE, 169 "%2s:%2s:%2s", 170 OBElements::GetSymbol(b->GetAtomicNum()), 171 OBElements::GetSymbol(a->GetAtomicNum()), 172 OBElements::GetSymbol(c->GetAtomicNum())); 173 } 174 else 175 { 176 snprintf(AngleString,BUFF_SIZE, 177 "%2s:%2s:%2s", 178 OBElements::GetSymbol(c->GetAtomicNum()), 179 OBElements::GetSymbol(a->GetAtomicNum()), 180 OBElements::GetSymbol(b->GetAtomicNum())); 181 } 182 AngleType[AngleString]=0; 183 } 184 int AngleIndex=0; 185 //Set the AtomType integer 186 for(intitr=AngleType.begin();intitr!=AngleType.end();++intitr) 187 { 188 AngleIndex++; 189 AngleType[intitr->first] = AngleIndex; 190 } 191 192 //dihedrals 193 mol.FindTorsions(); 194 int dihedralcount=0; 195 char DihedralString[BUFF_SIZE]; 196 map<string, int>DihedralType; 197 FOR_TORSIONS_OF_MOL(dihedral, mol) 198 { 199 dihedralcount++; 200 a = mol.GetAtom((*dihedral)[0]+1); 201 b = mol.GetAtom((*dihedral)[1]+1); 202 c = mol.GetAtom((*dihedral)[2]+1); 203 d = mol.GetAtom((*dihedral)[3]+1); 204 //place the heavier element first 205 //the same may have to be done with the inner two as well 206 if(a->GetAtomicNum() > d->GetAtomicNum() ) 207 { 208 snprintf(DihedralString,BUFF_SIZE, 209 "%2s:%2s:%2s:%2s", 210 OBElements::GetSymbol(a->GetAtomicNum()), 211 OBElements::GetSymbol(b->GetAtomicNum()), 212 OBElements::GetSymbol(c->GetAtomicNum()), 213 OBElements::GetSymbol(d->GetAtomicNum())); 214 } 215 else 216 { 217 snprintf(DihedralString,BUFF_SIZE, 218 "%2s:%2s:%2s:%2s", 219 OBElements::GetSymbol(d->GetAtomicNum()), 220 OBElements::GetSymbol(b->GetAtomicNum()), 221 OBElements::GetSymbol(c->GetAtomicNum()), 222 OBElements::GetSymbol(a->GetAtomicNum())); 223 } 224 DihedralType[DihedralString]=0; 225 } 226 int DihedralIndex=0; 227 //Set DihedralType integer 228 for(intitr=DihedralType.begin();intitr!=DihedralType.end();++intitr) 229 { 230 DihedralIndex++; 231 DihedralType[intitr->first] = DihedralIndex; 232 } 233 234 //The box lengths 235 vector3 vmin,vmax; 236 vmax.Set(-10E10,-10E10,-10E10); 237 vmin.Set( 10E10, 10E10, 10E10); 238 FOR_ATOMS_OF_MOL(a,mol) 239 { 240 if (a->x() < vmin.x()) 241 vmin.SetX(a->x()); 242 if (a->y() < vmin.y()) 243 vmin.SetY(a->y()); 244 if (a->z() < vmin.z()) 245 vmin.SetZ(a->z()); 246 247 if (a->x() > vmax.x()) 248 vmax.SetX(a->x()); 249 if (a->y() > vmax.y()) 250 vmax.SetY(a->y()); 251 if (a->z() > vmax.z()) 252 vmax.SetZ(a->z()); 253 } 254 255 //For now, a simple cube may be the best way to go. 256 //It may be necessary to set the boxlength to enforce 257 //a density. 258 const char *boxLn = pConv->IsOption("d",OBConversion::OUTOPTIONS); 259 double xlo,xhi; 260 char BoxString[BUFF_SIZE]; 261 if(boxLn) 262 { 263 xlo=-atof(boxLn); 264 xhi=atof(boxLn); 265 } 266 else 267 { 268 double cmin=10e-10; 269 double cmax=-10e10; 270 if ( vmax.x() > cmax ) cmax=vmax.x(); 271 if ( vmax.y() > cmax ) cmax=vmax.y(); 272 if ( vmax.z() > cmax ) cmax=vmax.z(); 273 if ( vmin.x() < cmin ) cmin=vmin.x(); 274 if ( vmin.y() < cmin ) cmin=vmin.y(); 275 if ( vmin.z() < cmin ) cmin=vmin.z(); 276 277 double length=cmax-cmin; 278 xlo = cmin -0.5;//- 0.01*length; 279 xhi = cmax +0.5;//+ 0.01*length; 280 } 281 snprintf(BoxString,BUFF_SIZE, 282 "%10.5f %10.5f xlo xhi\n%10.5f %10.5f ylo yhi\n%10.5f %10.5f zlo zhi\n", 283 xlo,xhi, 284 xlo,xhi, 285 xlo,xhi); 286 287 288 //%%%%%%%%%%%%%% Now writting the data file %%%%%%%%%%%%% 289 290 //The LAMMPS header stars here 291 ofs << "LAMMPS data file generated by OpenBabel" << endl; 292 ofs << mol.NumAtoms() << " atoms" << endl; 293 ofs << mol.NumBonds() << " bonds" << endl; 294 ofs << anglecount << " angles" << endl; // no NumAngles()? 295 ofs << dihedralcount << " dihedrals" << endl; 296 ofs << 0 << " impropers" << endl; 297 ofs << AtomType.size() << " atom types" << endl; 298 ofs << BondType.size() << " bond types" << endl; 299 ofs << AngleType.size() << " angle types" << endl; 300 ofs << DihedralType.size() << " dihedral types" << endl; 301 ofs << 0 << " improper types" << endl; 302 ofs << BoxString << endl; 303 304 305 //Write the atom types 306 ofs << endl << endl << "Masses" << endl << endl; 307 for(itr=AtomMass.begin();itr!=AtomMass.end();++itr) 308 { 309 double mass=itr->second; 310 ofs << AtomType[itr->first] << " " << mass << " # " << itr->first << endl; 311 } 312 ofs << endl; 313 314 315 316 //Set atomic charges for atom_style=full 317 //These are charges for the SPC water model 318 const char *selectCharges = pConv->IsOption("q",OBConversion::OUTOPTIONS); 319 map<string, double> AtomCharge; 320 for(itr=AtomMass.begin();itr!=AtomMass.end();++itr) 321 { 322 if(selectCharges) 323 { 324 if(strcmp(selectCharges,"spce")==0) 325 { 326 if(itr->second == 15.9994) 327 AtomCharge[itr->first] = -0.8472; 328 if(itr->second == 1.00794) 329 AtomCharge[itr->first] = 0.4236; 330 } 331 else if(strcmp(selectCharges,"spc")==0) 332 { 333 if(itr->second == 15.9994) 334 AtomCharge[itr->first] = -0.820; 335 if(itr->second == 1.00794) 336 AtomCharge[itr->first] = 0.410; 337 } 338 } 339 else 340 { 341 if(itr->second == 15.9994) 342 AtomCharge[itr->first] = -0.820; 343 if(itr->second == 1.00794) 344 AtomCharge[itr->first] = 0.410; 345 } 346 347 } 348 349 350 351 //Write atomic coordinates 352 vector<OBMol> mols; 353 vector<OBMol>::iterator molitr; 354 mols=mol.Separate(); 355 int atomcount=0; 356 int molcount=0; 357 ofs << endl; 358 ofs << "Atoms" << endl << endl; 359 snprintf(buffer,BUFF_SIZE,"#%3s %4s %4s %10s %10s %10s %10s\n", 360 "idx","mol","type","charge","x","y","z"); 361 //ofs << buffer; 362 for(molitr=mols.begin();molitr!=mols.end();++molitr) 363 { 364 molcount++; 365 FOR_ATOMS_OF_MOL(atom,*molitr) 366 { 367 int atomid=5; 368 double charge=0.5; 369 atomcount++; 370 ThisAtom=OBElements::GetSymbol(atom->GetAtomicNum()); 371 snprintf(buffer,BUFF_SIZE,"%-4d %4d %4d %10.5f %10.5f %10.5f %10.5f # %3s\n", 372 atomcount,molcount, 373 AtomType[ThisAtom], 374 AtomCharge[ThisAtom], 375 atom->GetX(), 376 atom->GetY(), 377 atom->GetZ(), 378 ThisAtom.c_str()); 379 ofs << buffer; 380 } 381 } 382 383 384 //Write Bonds 385 BondIndex=0; 386 ofs << endl << endl; 387 ofs << "Bonds" << endl; 388 ofs << endl; 389 snprintf(buffer,BUFF_SIZE, 390 "#%3s %4s %4s %4s\n", 391 "idx","type","atm1","atom2"); 392 //ofs << buffer; 393 FOR_BONDS_OF_MOL(bond,mol) 394 { 395 BondIndex++; 396 a = bond->GetBeginAtom(); 397 b = bond->GetEndAtom(); 398 //To let the unordered map determine unique keys, 399 //always put the heavier atom first 400 if(b->GetAtomicNum() > a->GetAtomicNum() ) 401 { 402 snprintf(BondString,BUFF_SIZE, 403 "%2s:%2s", 404 OBElements::GetSymbol(b->GetAtomicNum()), 405 OBElements::GetSymbol(a->GetAtomicNum())); 406 snprintf(buffer,BUFF_SIZE, 407 "%-4d %4d %4d %4d # %5s\n", 408 BondIndex, 409 BondType[BondString], 410 bond->GetEndAtomIdx(), 411 bond->GetBeginAtomIdx(), 412 BondString); 413 } 414 else 415 { 416 snprintf(BondString,BUFF_SIZE, 417 "%2s:%2s", 418 OBElements::GetSymbol(a->GetAtomicNum()), 419 OBElements::GetSymbol(b->GetAtomicNum())); 420 snprintf(buffer,BUFF_SIZE, 421 "%-4d %4d %4d %4d # %5s\n", 422 BondIndex, 423 BondType[BondString], 424 bond->GetBeginAtomIdx(), 425 bond->GetEndAtomIdx(), 426 BondString); 427 } 428 ofs << buffer; 429 } 430 ofs << endl; 431 432 //Write Angles 433 AngleIndex=0; 434 ofs << endl; 435 ofs << "Angles" << endl; 436 ofs << endl; 437 snprintf(buffer,BUFF_SIZE, 438 "#%3s %4s %4s %4s %4s\n", 439 "idx","type","atm1","atm2","atm3"); 440 //ofs << buffer; 441 FOR_ANGLES_OF_MOL(angle,mol) 442 { 443 AngleIndex++; 444 a = mol.GetAtom((*angle)[0]+1); 445 b = mol.GetAtom((*angle)[1]+1); 446 c = mol.GetAtom((*angle)[2]+1); 447 //Origanization trick: 448 //1. The atom "a" is acutally the middle 449 //2. Always write the heavy element first 450 if(b->GetAtomicNum() > c->GetAtomicNum() ) 451 { 452 snprintf(AngleString,BUFF_SIZE, 453 "%2s:%2s:%2s", 454 OBElements::GetSymbol(b->GetAtomicNum()), 455 OBElements::GetSymbol(a->GetAtomicNum()), 456 OBElements::GetSymbol(c->GetAtomicNum())); 457 snprintf(buffer,BUFF_SIZE, 458 "%-4d %4d %4d %4d %4d # %8s\n", 459 AngleIndex, 460 AngleType[AngleString], 461 b->GetIdx(), 462 a->GetIdx(), 463 c->GetIdx(), 464 AngleString); 465 } 466 else 467 { 468 snprintf(AngleString,BUFF_SIZE, 469 "%2s:%2s:%2s", 470 OBElements::GetSymbol(c->GetAtomicNum()), 471 OBElements::GetSymbol(a->GetAtomicNum()), 472 OBElements::GetSymbol(b->GetAtomicNum())); 473 snprintf(buffer,BUFF_SIZE, 474 "%-4d %4d %4d %4d %4d # %8s\n", 475 AngleIndex, 476 AngleType[AngleString], 477 c->GetIdx(), 478 a->GetIdx(), 479 b->GetIdx(), 480 AngleString); 481 } 482 ofs << buffer; 483 484 } 485 ofs << endl; 486 487 //Write Dihedrals 488 if(dihedralcount>0) 489 { 490 ofs << endl; 491 ofs << "Dihedrals" << endl; 492 ofs << endl; 493 snprintf(buffer,BUFF_SIZE, 494 "#%3s %4s %4s %4s %4s %4s\n", 495 "idx","type","atm1","atm2","atm3","atm4"); 496 //ofs << buffer; 497 DihedralIndex=0; 498 FOR_TORSIONS_OF_MOL(dihedral, mol) 499 { 500 DihedralIndex++; 501 a = mol.GetAtom((*dihedral)[0]+1); 502 b = mol.GetAtom((*dihedral)[1]+1); 503 c = mol.GetAtom((*dihedral)[2]+1); 504 d = mol.GetAtom((*dihedral)[3]+1); 505 //place the heavier element first 506 //the same may have to be done with the inner two as well 507 if(a->GetAtomicNum() > d->GetAtomicNum() ) 508 { 509 snprintf(DihedralString,BUFF_SIZE, 510 "%2s:%2s:%2s:%2s", 511 OBElements::GetSymbol(a->GetAtomicNum()), 512 OBElements::GetSymbol(b->GetAtomicNum()), 513 OBElements::GetSymbol(c->GetAtomicNum()), 514 OBElements::GetSymbol(d->GetAtomicNum())); 515 snprintf(buffer,BUFF_SIZE, 516 "%-4d %4d %4d %4d %4d %4d # %11s\n", 517 DihedralIndex, 518 DihedralType[DihedralString], 519 a->GetIdx(), 520 b->GetIdx(), 521 c->GetIdx(), 522 d->GetIdx(), 523 DihedralString); 524 } 525 else 526 { 527 snprintf(DihedralString,BUFF_SIZE, 528 "%2s:%2s:%2s:%2s", 529 OBElements::GetSymbol(d->GetAtomicNum()), 530 OBElements::GetSymbol(b->GetAtomicNum()), 531 OBElements::GetSymbol(c->GetAtomicNum()), 532 OBElements::GetSymbol(a->GetAtomicNum())); 533 snprintf(buffer,BUFF_SIZE, 534 "%-4d %4d %4d %4d %4d %4d # %11s\n", 535 DihedralIndex, 536 DihedralType[DihedralString], 537 d->GetIdx(), 538 b->GetIdx(), 539 c->GetIdx(), 540 a->GetIdx(), 541 DihedralString); 542 } 543 ofs << buffer; 544 } 545 } 546 547 548 return(true); 549 } 550 551 } //namespace OpenBabel 552