1 // MODEL.H : the "model" class that stores the "atom" and "bond" objects. 2 3 // Copyright (C) 1998 Tommi Hassinen. 4 5 // This package is free software; you can redistribute it and/or modify 6 // it under the terms of the GNU General Public License as published by 7 // the Free Software Foundation; either version 2 of the License, or 8 // (at your option) any later version. 9 10 // This package is distributed in the hope that it will be useful, 11 // but WITHOUT ANY WARRANTY; without even the implied warranty of 12 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 13 // GNU General Public License for more details. 14 15 // You should have received a copy of the GNU General Public License 16 // along with this package; if not, write to the Free Software 17 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA 18 19 /*################################################################################################*/ 20 21 #ifndef MODEL_H 22 #define MODEL_H 23 24 #include "libghemicaldefine.h" 25 26 class model; 27 class crd_set; 28 29 // the readpdb classes here are not fundamental, but are located here for convenience... 30 // the readpdb classes here are not fundamental, but are located here for convenience... 31 // the readpdb classes here are not fundamental, but are located here for convenience... 32 33 class readpdb_mdata; 34 struct readpdb_mdata_chain; 35 struct readpdb_data_atom; 36 struct readpdb_data_ssbond; 37 38 /*################################################################################################*/ 39 40 class sequencebuilder; // seqbuild.h 41 class resonance_structures; // resonance.h 42 43 class setup; // engine.h 44 class engine; // engine.h 45 46 class geomopt_param; // geomopt.h 47 48 class moldyn; // moldyn.h 49 class moldyn_param; // moldyn.h 50 51 #include "atom.h" 52 #include "bond.h" 53 54 #include "constraint.h" 55 56 #include "chn_info.h" 57 58 #include <list> 59 #include <vector> 60 #include <fstream> 61 #include <iostream> 62 #include <algorithm> 63 using namespace std; 64 65 #define NOT_FOUND 0x7fffffff // numeric_limits<i32s>::max()?!?!?! 66 67 #define SF_NUM_TYPES 37 68 69 /*################################################################################################*/ 70 71 /// The model class contains information about all atoms and bonds in a model. 72 73 class model 74 { 75 protected: 76 77 setup * current_setup; ///< Must always be a valid pointer!!! 78 resonance_structures * rs; ///< This pointer may be either NULL or a valid one. 79 80 list<atom> atom_list; 81 list<bond> bond_list; 82 83 list<constraint_dst> const_D_list; 84 85 vector<crd_set *> cs_vector; ///< This determines how many crd_sets there are in the model. 86 i32u crd_table_size_glob; ///< This determines how big the crd_table arrays are; always >= than above!!! 87 88 bool is_index_clean; ///< Some flags that show which information is up-to-date... 89 bool is_groups_clean; ///< Some flags that show which information is up-to-date... 90 bool is_groups_sorted; ///< Some flags that show which information is up-to-date... 91 92 i32s qm_total_charge; 93 i32s qm_current_orbital; 94 95 bool use_boundary_potential; 96 f64 saved_boundary_potential_rad_solute; 97 f64 saved_boundary_potential_rad_solvent; 98 99 bool use_periodic_boundary_conditions; 100 f64 saved_periodic_box_HALFdim[3]; 101 102 i32s nmol; 103 vector<chn_info> * ref_civ; // vector<chn_info *> ?!?!?! 104 105 ifstream * trajfile; // trajectory files... 106 i32s traj_num_atoms; // trajectory files... 107 i32s total_traj_frames; // trajectory files... 108 i32s current_traj_frame; // trajectory files... 109 110 public: 111 112 // the verbosity levels control the amount of log messages: 113 // 0 = silent ; 1 = errors ; 2 = warnings ; 3 = notices ; 4 = debug... 114 115 int verbosity; 116 117 // the ecomp stuff is here: 118 // ^^^^^^^^^^^^^^^^^^^^^^^^ 119 120 public: 121 122 bool ecomp_enabled; 123 124 protected: 125 126 vector<const char *> ecomp_grp_names; 127 128 /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// 129 /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// 130 // for thermodynamics stuff, several "topologies" (and therefore also several elements, atom types etc...) are needed. 131 // plan1) define two different elements, bondtypes for the topologies; make it possible to create MM "engine" objects for 132 // either of the topologies, so that there are exactly the same terms in both objects. then make a third "engine" object with 133 // energy term parameters intermediate of those two topologies -> one can softly change 1st topology to the 2nd one. 134 // plan2) do not really store/maintain two sets of elements/bondtypes/etc but instead make two separate systems, only slightly 135 // modified. later, it's possible to compare them and find the differences (at least if some "anchor" atom pairs are given). 136 // then, add dummy atoms to both so that exactly same terms appear in engine classes. make this even in many steps?!?!?! 137 // still many practical questions remain; do dummy atoms interfere in typerules? how to handle dummy atoms in calculations 138 // (for example, in a case O-H -> Cl-du an atom will disappear??? do not accept structures like du-du-...???) 139 /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// 140 /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// 141 142 public: 143 144 static const char libversion[16]; 145 static char libdata_path[256]; // libghemical_init() will modify this... 146 147 static const char sf_symbols[20 + 1]; 148 static const i32s sf_types[SF_NUM_TYPES]; 149 static const bool sf_is_polar[SF_NUM_TYPES]; 150 151 static sequencebuilder * amino_builder; 152 static sequencebuilder * nucleic_builder; 153 154 friend void CopyCRD(model *, engine *, i32u); 155 friend void CopyCRD(engine *, model *, i32u); 156 157 friend fGL plot_GetESPValue(fGL *, model *, fGL *); 158 friend fGL plot_GetVDWSValue(fGL *, model *, fGL *); 159 160 friend void DefineSecondaryStructure(model *); 161 friend f64 HBondEnergy(model *, i32s *, i32s *); 162 163 friend class sequencebuilder; 164 friend class mfinder; 165 166 friend class engine; 167 friend class engine_bp; 168 friend class engine_pbc; 169 170 friend class setup1_qm; 171 friend class eng1_qm_mopac; 172 173 friend class setup1_mm; 174 friend class eng1_mm_pbc; 175 176 friend class eng1_mm_tripos52_nbt_bp; 177 friend class eng1_mm_tripos52_nbt_mim; 178 179 friend class eng1_mm_default_nbt_bp; 180 friend class eng1_mm_default_nbt_mim; 181 friend class default_tables; 182 183 friend class setup1_sf; 184 185 friend class moldyn; 186 187 public: 188 189 model(void); 190 virtual ~model(void); 191 192 virtual void ThreadLock(void); 193 virtual void ThreadUnlock(void); 194 195 virtual void NoThreadsIterate(void); 196 197 virtual bool SetProgress(double, double *); 198 199 virtual void Message(const char *); 200 virtual void WarningMessage(const char *); 201 virtual void ErrorMessage(const char *); 202 203 setup * GetCurrentSetup(void); 204 void ReplaceCurrentSetup(setup *); 205 206 resonance_structures * GetRS(void); 207 208 void CreateRS(void); 209 void DestroyRS(void); 210 211 static void OpenLibDataFile(ifstream &, bool, const char *); 212 213 // what to do for this one??? 214 // what to do for this one??? 215 // what to do for this one??? 216 /// The project-class will override this function... 217 virtual void UpdateAllGraphicsViews(bool = false) { } 218 219 // what to do for this one??? 220 // what to do for this one??? 221 // what to do for this one??? 222 /// Add a message string to the logfile. This is just a default for console... PrintToLog(const char * p1)223 virtual void PrintToLog(const char * p1) { cout << "PrintToLog: " << p1 << endl; } 224 225 /// This will return the number of coordinate sets. 226 /** It is supposed that at least one coordinate set exists all the time!!! */ 227 i32u GetCRDSetCount(void); 228 229 // what to do for this one??? 230 // what to do for this one??? 231 // what to do for this one??? 232 bool GetCRDSetVisible(i32u); 233 void SetCRDSetVisible(i32u, bool); 234 235 void PushCRDSets(i32u); 236 void PopCRDSets(i32u); 237 238 void CopyCRDSet(i32u, i32u); 239 void SwapCRDSets(i32u, i32u); 240 241 void CenterCRDSet(i32u, bool); 242 void OrientCRDSet(i32u, bool, fGL *); 243 244 void ReserveCRDSets(i32u); 245 246 virtual void DiscardCurrEng(void); 247 248 // rename this!!! 249 // rename this!!! 250 // rename this!!! 251 void SetupPlotting(void); 252 253 // methods for adding new atoms and bonds: 254 // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 255 // 2008-07-28 : AddAtom() was renamed into AddAtom_lg() in order to avoid problems caused by 256 // an another AddAtom()-something included in the Win32 API ; it is hard to understand why this is 257 // so bad problem in win32 but it is ; not even limited to MSVC -> perhaps about name mangling etc... 258 259 virtual void AddAtom_lg(atom &); ///< This will just push new atom to the atom list. 260 virtual void RemoveAtom(iter_al); ///< This will delete all bonds associated with this atom, and erase atom from the list... 261 262 virtual void AddBond(bond &); ///< This will add neighbor infos for both atoms and add the new bond into the bond list. 263 virtual void RemoveBond(iter_bl); ///< This will remove infos from the atoms and erase bond from the bond list. 264 265 virtual void AddConstraint(constraint_dst &); 266 virtual void RemoveConstraint(iter_CDl); 267 268 void SystemWasModified(void); ///< This should be called ALWAYS if ANY modification is done to the system. Automatically called by Add/Remove/Atom/Bond. GUI should change if change in element etc... 269 270 void ClearModel(void); ///< This will remove all atoms and bonds. 271 272 // methods for accessing atom/bond lists: 273 // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 274 // ALL THESE ITERATORS SHOULD BE READ-ONLY!!! is it possible to modify the containers using iterators only??? 275 GetAtomsBegin(void)276 iter_al GetAtomsBegin(void) { return atom_list.begin(); } //const; GetAtomsEnd(void)277 iter_al GetAtomsEnd(void) { return atom_list.end(); } //const; GetLastAtom(void)278 atom & GetLastAtom(void) { return atom_list.back(); } GetAtomCount(void)279 int GetAtomCount(void) { return (int) atom_list.size(); } ///< This will return atom_list.size(). 280 GetBondsBegin(void)281 iter_bl GetBondsBegin(void) { return bond_list.begin(); } //const; GetBondsEnd(void)282 iter_bl GetBondsEnd(void) { return bond_list.end(); } //const; GetLastBond(void)283 bond & GetLastBond(void) { return bond_list.back(); } GetBondCount(void)284 int GetBondCount(void) { return (int) bond_list.size(); } ///< This will return bond_list.size(). 285 GetConstD_begin(void)286 iter_CDl GetConstD_begin(void) { return const_D_list.begin(); } //const; GetConstD_end(void)287 iter_CDl GetConstD_end(void) { return const_D_list.end(); } //const; GetLastConstD(void)288 constraint_dst & GetLastConstD(void) { return const_D_list.back(); } GetConstD_count(void)289 int GetConstD_count(void) { return (int) const_D_list.size(); } ///< This will return const_D_list.size(). 290 291 // methods for ??? 292 // ^^^^^^^^^^^^^^^ 293 GetQMTotalCharge(void)294 i32s GetQMTotalCharge(void) { return qm_total_charge; } SetQMTotalCharge(i32s tc)295 void SetQMTotalCharge(i32s tc) { qm_total_charge = tc; } 296 GetQMCurrentOrbital(void)297 i32s GetQMCurrentOrbital(void) { return qm_current_orbital; } SetQMCurrentOrbital(i32s orb)298 void SetQMCurrentOrbital(i32s orb) { qm_current_orbital = orb; } 299 GetMoleculeCount(void)300 i32s GetMoleculeCount(void) { return nmol; } ///< This will return nmol (see UpdateGroups() for more info). 301 IsEmpty(void)302 bool IsEmpty(void) { return (bond_list.empty() && atom_list.empty()); } ///< This will return whether the system is empty. 303 IsIndexClean(void)304 bool IsIndexClean(void) { return is_index_clean; } IsGroupsClean(void)305 bool IsGroupsClean(void) { return is_groups_clean; } IsGroupsSorted(void)306 bool IsGroupsSorted(void) { return is_groups_sorted; } 307 308 iter_al FindAtomByIndex(i32s); 309 iter_CDl FindAtomConstraint(atom &); 310 311 void GetRange(i32s, i32s, iter_al *); ///< This is just a default version of GetRange() using the full range of atom list iterators... 312 313 /// GetRange is used to get a range of atoms that form molecules, residues etc... 314 /** This will reduce the initial range of two atom list iterators to some subrange 315 with certain values in certain atom::id[]-fields. Before using this you MUST call 316 model::UpdateGroups() to arrange the atom list! What the above explanation really tries 317 to say, is that using this function you can pick up a certain part of the model; for 318 example a molecule, or a chain in a macromolecule, or a range of residue in a chain. */ 319 void GetRange(i32s, iter_al *, i32s, iter_al *); 320 321 void GetRange(i32s, iter_bl *); ///< This GetRange() gives a range of bond iterators for a molecule... 322 323 i32s FindPath(atom *, atom *, i32s, i32s, i32s = 0); 324 vector<bond *> * FindPathV(atom *, atom *, i32s, i32s, i32s = 0); 325 bool FindRing(atom *, atom *, signed char *, i32s, i32s, i32s = 0); 326 327 /** Adding or removing atoms or bonds will generally scramble the grouping, 328 and when this happens one should call this function to discard all grouping information. */ 329 virtual void InvalidateGroups(void); 330 331 void UpdateIndex(void); 332 333 /** This will set the atom ID numbers so that molecules form groups. Will also validate model::nmol 334 but will not yet permit the use of model::GetRange()-functions since atom list is not sorted. */ 335 void UpdateGroups(void); 336 337 /** This will group the atom list so that molecules (and chains/residues if defined) will form 338 continuous groups. Will permit the use of model::GetRange()-functions. */ 339 void SortGroups(void); 340 341 /** This will validate ref_civ and fill it with chain descriptions using all sequencebuilder objects. 342 Finally it will call model::SortGroups() so that chains and residues in the atom list will be ordered 343 sequentially in contiguous blocks. */ 344 virtual void UpdateChains(bool = false); 345 GetCI(void)346 vector<chn_info> * GetCI(void) { return ref_civ; } 347 348 private: 349 350 /** This will set molecule numbers quickly using a recursive search algorithm. 351 This is private because only model::UpdateGroups() should use this... */ 352 void GatherAtoms(atom *, i32s); 353 354 atom * cp_FindAtom(iter_al *, i32s); ///< only for CheckProtonation()... 355 356 /** This will update atom::formal_charge using ref_civ->p_state information. 357 This is private because only model::AddHydrogens() should use this... */ 358 void CheckProtonation(void); 359 360 public: 361 362 void AddHydrogens(void); 363 void AddHydrogens(atom *); 364 365 void RemoveHydrogens(void); 366 367 void SolvateBox(fGL, fGL, fGL, fGL = 1.0, model * = NULL, const char * = NULL); 368 void SolvateSphere(fGL, fGL, fGL = 1.0, model * = NULL); 369 fGL S_Initialize(fGL, model **); 370 371 // here we have a set of Do???()-functions. the idea is that there is a set 372 // of features that we wish to behave EXACTLY same way in each target/platform. 373 // we then create a Do???()-function for the feature, and hide the details of 374 // the user interface in a set of virtual functions. 375 376 /** This will perform an energy calculation, and report the result. */ 377 void DoEnergy(void); 378 379 /** This will perform geometry optimization. */ 380 void DoGeomOpt(geomopt_param &, bool); 381 382 /** This will perform molecular dynamics. */ 383 void DoMolDyn(moldyn_param &, bool); 384 385 /** This will perform a random search using torsions as variables. */ 386 void DoRandomSearch(i32s, i32s, bool); 387 388 /** This will perform a systematic search using torsions as variables. */ 389 void DoSystematicSearch(i32s, i32s, bool); 390 391 /** This will perform a MonteCarlo search using torsions as variables. */ 392 void DoMonteCarloSearch(i32s, i32s, i32s, bool); 393 394 public: 395 396 /// A function for reading in Brookhaven PDB/ENT files. 397 /** The readpdb functions here are used to import PDB files as correctly as possible. 398 The PDB files represent experimental results, and in many cases the structures in files 399 have gapped and/or incomplete sequences, incomplete residues, and so on. 400 401 The readpdb functions do the import in two stages: in first stage read in the "metadata" 402 (all headers and remarks about the data including the original sequence), and in second stage 403 read in the data as correctly as possible. later, results from these two can be compared, for 404 example to evaluate quality of the data or to match the data with records in other databases. */ 405 readpdb_mdata * readpdb_ReadMData(const char *); 406 407 void readpdb_ReadData(const char *, readpdb_mdata *, i32s); 408 i32s readpdb_ReadData_sub1(vector<readpdb_data_atom> &, i32s *, const char *, bool); 409 void readpdb_ReadData_sub2(vector<readpdb_data_atom> &, i32s *, const char *, const char *, char); 410 411 // methods related to MD trajectories... 412 // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 413 void OpenTrajectory(const char *); 414 void ReadTrajectoryFrame(void); 415 void CloseTrajectory(void); 416 417 i32s GetTotalFrames(void); 418 ifstream * GetTrajectoryFile(void); 419 420 i32s GetCurrentFrame(void); 421 void SetCurrentFrame(i32s); 422 423 void WriteTrajectoryHeader(ofstream &, int); 424 void WriteTrajectoryFrame(ofstream &, moldyn *); 425 426 void EvaluateBFact(void); 427 void EvaluateDiffConst(double); 428 429 // methods related to E-groups (implemented in eng-classes). 430 // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 431 int ecomp_AddGroup(const char *); 432 void ecomp_DeleteGroups(void); 433 }; 434 435 /*################################################################################################*/ 436 437 /// A coordinate-set class. 438 439 class crd_set 440 { 441 protected: 442 443 char * description; 444 445 // bring also the coloring information here??? -> would allow different colors for different crdsets!!! 446 // bring also the coloring information here??? -> would allow different colors for different crdsets!!! 447 // bring also the coloring information here??? -> would allow different colors for different crdsets!!! 448 449 public: 450 451 // these are tricky to protect, since these are still developing and are not yet used much. 452 // so, these are public as long as some reasonable strategy is found... 453 454 fGL accum_weight; 455 fGL accum_value; 456 bool visible; 457 458 public: 459 460 crd_set(void); 461 crd_set(const crd_set &); 462 ~crd_set(void); 463 464 void SetDescription(const char *); 465 }; 466 467 /*################################################################################################*/ 468 469 // define struct readpdb_mdata_chain before class readpdb_mdata, since the latter uses 470 // former in some inline functions... 471 472 // how to best relate readpdb_mdata_chain and chn_info !??!?!?!?!? 473 // maybe just by storing the alpha-carbon pointers here... 474 475 struct readpdb_mdata_chain 476 { 477 char chn_id; 478 char * seqres; 479 480 vector<i32s> missing_residues; 481 vector<atom *> alpha_carbons; 482 }; 483 484 // class readpdb_mdata is a class just to make the memory management easier. the data members in 485 // the class are filled in model::readpdb_ReadMData(), and at end the object is just to be deleted. 486 487 class readpdb_mdata 488 { 489 public: 490 491 vector<readpdb_mdata_chain *> chn_vector; 492 493 public: 494 readpdb_mdata(void)495 readpdb_mdata(void) 496 { 497 } 498 ~readpdb_mdata(void)499 ~readpdb_mdata(void) 500 { 501 for (i32u n1 = 0;n1 < chn_vector.size();n1++) 502 { 503 delete[] chn_vector[n1]->seqres; // delete the sequence... 504 delete chn_vector[n1]; // delete the whole record... 505 } 506 } 507 }; 508 509 // READPDB_MAX_CRDSETS is relevant only if READPDB_ENABLE_MULTIPLE_CRDSETS is defined... 510 // READPDB_MAX_CRDSETS is relevant only if READPDB_ENABLE_MULTIPLE_CRDSETS is defined... 511 // READPDB_MAX_CRDSETS is relevant only if READPDB_ENABLE_MULTIPLE_CRDSETS is defined... 512 513 #define READPDB_MAX_CRDSETS 10 514 515 struct readpdb_data_atom 516 { 517 char chn_id; 518 519 i32s res_num; 520 char res_name[5]; 521 char atm_name[5]; 522 fGL crd[READPDB_MAX_CRDSETS][3]; 523 524 atom * ref; 525 }; 526 527 struct readpdb_data_ssbond 528 { 529 char chn_id; 530 i32s res_num; 531 532 atom * ref; 533 }; 534 535 /*################################################################################################*/ 536 537 void libghemical_init(const char *); 538 539 /*################################################################################################*/ 540 541 #endif // MODEL_H 542 543 // eof 544