| 1 | /** \file atom.cpp | 
|---|
| 2 | * | 
|---|
| 3 | * Function implementations for the class atom. | 
|---|
| 4 | * | 
|---|
| 5 | */ | 
|---|
| 6 |  | 
|---|
| 7 | #include "atom.hpp" | 
|---|
| 8 | #include "bond.hpp" | 
|---|
| 9 | #include "element.hpp" | 
|---|
| 10 | #include "memoryallocator.hpp" | 
|---|
| 11 | #include "vector.hpp" | 
|---|
| 12 |  | 
|---|
| 13 | /************************************* Functions for class atom *************************************/ | 
|---|
| 14 |  | 
|---|
| 15 | /** Constructor of class atom. | 
|---|
| 16 | */ | 
|---|
| 17 | atom::atom() | 
|---|
| 18 | { | 
|---|
| 19 | father = this;  // generally, father is itself | 
|---|
| 20 | previous = NULL; | 
|---|
| 21 | next = NULL; | 
|---|
| 22 | Ancestor = NULL; | 
|---|
| 23 | type = NULL; | 
|---|
| 24 | sort = NULL; | 
|---|
| 25 | FixedIon = 0; | 
|---|
| 26 | GraphNr = -1; | 
|---|
| 27 | ComponentNr = NULL; | 
|---|
| 28 | IsCyclic = false; | 
|---|
| 29 | SeparationVertex = false; | 
|---|
| 30 | LowpointNr = -1; | 
|---|
| 31 | AdaptiveOrder = 0; | 
|---|
| 32 | MaxOrder = false; | 
|---|
| 33 | // set LCNode::Vector to our Vector | 
|---|
| 34 | node = &x; | 
|---|
| 35 | }; | 
|---|
| 36 |  | 
|---|
| 37 | /** Constructor of class atom. | 
|---|
| 38 | */ | 
|---|
| 39 | atom::atom(atom *pointer) | 
|---|
| 40 | { | 
|---|
| 41 | Name = NULL; | 
|---|
| 42 | previous = NULL; | 
|---|
| 43 | next = NULL; | 
|---|
| 44 | father = pointer;  // generally, father is itself | 
|---|
| 45 | Ancestor = NULL; | 
|---|
| 46 | GraphNr = -1; | 
|---|
| 47 | ComponentNr = NULL; | 
|---|
| 48 | IsCyclic = false; | 
|---|
| 49 | SeparationVertex = false; | 
|---|
| 50 | LowpointNr = -1; | 
|---|
| 51 | AdaptiveOrder = 0; | 
|---|
| 52 | MaxOrder = false; | 
|---|
| 53 | type = pointer->type;  // copy element of atom | 
|---|
| 54 | x.CopyVector(&pointer->x); // copy coordination | 
|---|
| 55 | v.CopyVector(&pointer->v); // copy velocity | 
|---|
| 56 | FixedIon = pointer->FixedIon; | 
|---|
| 57 | nr = -1; | 
|---|
| 58 | sort = &nr; | 
|---|
| 59 | node = &x; | 
|---|
| 60 | } | 
|---|
| 61 |  | 
|---|
| 62 |  | 
|---|
| 63 | /** Destructor of class atom. | 
|---|
| 64 | */ | 
|---|
| 65 | atom::~atom() | 
|---|
| 66 | { | 
|---|
| 67 | Free<int>(&ComponentNr, "atom::~atom: *ComponentNr"); | 
|---|
| 68 | Free<char>(&Name, "atom::~atom: *Name"); | 
|---|
| 69 | Trajectory.R.clear(); | 
|---|
| 70 | Trajectory.U.clear(); | 
|---|
| 71 | Trajectory.F.clear(); | 
|---|
| 72 | }; | 
|---|
| 73 |  | 
|---|
| 74 |  | 
|---|
| 75 | /** Climbs up the father list until NULL, last is returned. | 
|---|
| 76 | * \return true father, i.e. whose father points to itself, NULL if it could not be found or has none (added hydrogen) | 
|---|
| 77 | */ | 
|---|
| 78 | atom *atom::GetTrueFather() | 
|---|
| 79 | { | 
|---|
| 80 | atom *walker = this; | 
|---|
| 81 | do { | 
|---|
| 82 | if (walker == walker->father) // top most father is the one that points on itself | 
|---|
| 83 | break; | 
|---|
| 84 | walker = walker->father; | 
|---|
| 85 | } while (walker != NULL); | 
|---|
| 86 | return walker; | 
|---|
| 87 | }; | 
|---|
| 88 |  | 
|---|
| 89 | /** Sets father to itself or its father in case of copying a molecule. | 
|---|
| 90 | */ | 
|---|
| 91 | void atom::CorrectFather() | 
|---|
| 92 | { | 
|---|
| 93 | if (father->father == father)   // same atom in copy's father points to itself | 
|---|
| 94 | father = this;  // set father to itself (copy of a whole molecule) | 
|---|
| 95 | else | 
|---|
| 96 | father = father->father;  // set father to original's father | 
|---|
| 97 |  | 
|---|
| 98 | }; | 
|---|
| 99 |  | 
|---|
| 100 | /** Check whether father is equal to given atom. | 
|---|
| 101 | * \param *ptr atom to compare father to | 
|---|
| 102 | * \param **res return value (only set if atom::father is equal to \a *ptr) | 
|---|
| 103 | */ | 
|---|
| 104 | void atom::EqualsFather ( atom *ptr, atom **res ) | 
|---|
| 105 | { | 
|---|
| 106 | if ( ptr == father ) | 
|---|
| 107 | *res = this; | 
|---|
| 108 | }; | 
|---|
| 109 |  | 
|---|
| 110 | /** Checks whether atom is within the given box. | 
|---|
| 111 | * \param offset offset to box origin | 
|---|
| 112 | * \param *parallelepiped box matrix | 
|---|
| 113 | * \return true - is inside, false - is not | 
|---|
| 114 | */ | 
|---|
| 115 | bool atom::IsInParallelepiped(Vector offset, double *parallelepiped) | 
|---|
| 116 | { | 
|---|
| 117 | return (node->IsInParallelepiped(offset, parallelepiped)); | 
|---|
| 118 | }; | 
|---|
| 119 |  | 
|---|
| 120 | /** Output of a single atom. | 
|---|
| 121 | * \param ElementNo cardinal number of the element | 
|---|
| 122 | * \param AtomNo cardinal number among these atoms of the same element | 
|---|
| 123 | * \param *out stream to output to | 
|---|
| 124 | * \param *comment commentary after '#' sign | 
|---|
| 125 | * \return true - \a *out present, false - \a *out is NULL | 
|---|
| 126 | */ | 
|---|
| 127 | bool atom::Output(ofstream *out, int ElementNo, int AtomNo, const char *comment) const | 
|---|
| 128 | { | 
|---|
| 129 | if (out != NULL) { | 
|---|
| 130 | *out << "Ion_Type" << ElementNo << "_" << AtomNo << "\t"  << fixed << setprecision(9) << showpoint; | 
|---|
| 131 | *out << x.x[0] << "\t" << x.x[1] << "\t" << x.x[2]; | 
|---|
| 132 | *out << "\t" << FixedIon; | 
|---|
| 133 | if (v.Norm() > MYEPSILON) | 
|---|
| 134 | *out << "\t" << scientific << setprecision(6) << v.x[0] << "\t" << v.x[1] << "\t" << v.x[2] << "\t"; | 
|---|
| 135 | if (comment != NULL) | 
|---|
| 136 | *out << " # " << comment << endl; | 
|---|
| 137 | else | 
|---|
| 138 | *out << " # molecule nr " << nr << endl; | 
|---|
| 139 | return true; | 
|---|
| 140 | } else | 
|---|
| 141 | return false; | 
|---|
| 142 | }; | 
|---|
| 143 | bool atom::Output(ofstream *out, int *ElementNo, int *AtomNo, const char *comment) | 
|---|
| 144 | { | 
|---|
| 145 | AtomNo[type->Z]++;  // increment number | 
|---|
| 146 | if (out != NULL) { | 
|---|
| 147 | *out << "Ion_Type" << ElementNo[type->Z] << "_" << AtomNo[type->Z] << "\t"  << fixed << setprecision(9) << showpoint; | 
|---|
| 148 | *out << x.x[0] << "\t" << x.x[1] << "\t" << x.x[2]; | 
|---|
| 149 | *out << "\t" << FixedIon; | 
|---|
| 150 | if (v.Norm() > MYEPSILON) | 
|---|
| 151 | *out << "\t" << scientific << setprecision(6) << v.x[0] << "\t" << v.x[1] << "\t" << v.x[2] << "\t"; | 
|---|
| 152 | if (comment != NULL) | 
|---|
| 153 | *out << " # " << comment << endl; | 
|---|
| 154 | else | 
|---|
| 155 | *out << " # molecule nr " << nr << endl; | 
|---|
| 156 | return true; | 
|---|
| 157 | } else | 
|---|
| 158 | return false; | 
|---|
| 159 | }; | 
|---|
| 160 |  | 
|---|
| 161 | /** Output of a single atom as one lin in xyz file. | 
|---|
| 162 | * \param *out stream to output to | 
|---|
| 163 | * \return true - \a *out present, false - \a *out is NULL | 
|---|
| 164 | */ | 
|---|
| 165 | bool atom::OutputXYZLine(ofstream *out) const | 
|---|
| 166 | { | 
|---|
| 167 | if (out != NULL) { | 
|---|
| 168 | *out << type->symbol << "\t" << x.x[0] << "\t" << x.x[1] << "\t" << x.x[2] << "\t" << endl; | 
|---|
| 169 | return true; | 
|---|
| 170 | } else | 
|---|
| 171 | return false; | 
|---|
| 172 | }; | 
|---|
| 173 |  | 
|---|
| 174 | /** Output of a single atom as one lin in xyz file. | 
|---|
| 175 | * \param *out stream to output to | 
|---|
| 176 | * \param *ElementNo array with ion type number in the config file this atom's element shall have | 
|---|
| 177 | * \param *AtomNo array with atom number in the config file this atom shall have, is increase by one automatically | 
|---|
| 178 | * \param step Trajectory time step to output | 
|---|
| 179 | * \return true - \a *out present, false - \a *out is NULL | 
|---|
| 180 | */ | 
|---|
| 181 | bool atom::OutputTrajectory(ofstream *out, int *ElementNo, int *AtomNo, int step) const | 
|---|
| 182 | { | 
|---|
| 183 | AtomNo[type->Z]++; | 
|---|
| 184 | if (out != NULL) { | 
|---|
| 185 | *out << "Ion_Type" << ElementNo[type->Z] << "_" << AtomNo[type->Z] << "\t"  << fixed << setprecision(9) << showpoint; | 
|---|
| 186 | *out << Trajectory.R.at(step).x[0] << "\t" << Trajectory.R.at(step).x[1] << "\t" << Trajectory.R.at(step).x[2]; | 
|---|
| 187 | *out << "\t" << FixedIon; | 
|---|
| 188 | if (Trajectory.U.at(step).Norm() > MYEPSILON) | 
|---|
| 189 | *out << "\t" << scientific << setprecision(6) << Trajectory.U.at(step).x[0] << "\t" << Trajectory.U.at(step).x[1] << "\t" << Trajectory.U.at(step).x[2] << "\t"; | 
|---|
| 190 | if (Trajectory.F.at(step).Norm() > MYEPSILON) | 
|---|
| 191 | *out << "\t" << scientific << setprecision(6) << Trajectory.F.at(step).x[0] << "\t" << Trajectory.F.at(step).x[1] << "\t" << Trajectory.F.at(step).x[2] << "\t"; | 
|---|
| 192 | *out << "\t# Number in molecule " << nr << endl; | 
|---|
| 193 | return true; | 
|---|
| 194 | } else | 
|---|
| 195 | return false; | 
|---|
| 196 | }; | 
|---|
| 197 |  | 
|---|
| 198 | /** Output of a single atom as one lin in xyz file. | 
|---|
| 199 | * \param *out stream to output to | 
|---|
| 200 | * \param step Trajectory time step to output | 
|---|
| 201 | * \return true - \a *out present, false - \a *out is NULL | 
|---|
| 202 | */ | 
|---|
| 203 | bool atom::OutputTrajectoryXYZ(ofstream *out, int step) const | 
|---|
| 204 | { | 
|---|
| 205 | if (out != NULL) { | 
|---|
| 206 | *out << type->symbol << "\t"; | 
|---|
| 207 | *out << Trajectory.R.at(step).x[0] << "\t"; | 
|---|
| 208 | *out << Trajectory.R.at(step).x[1] << "\t"; | 
|---|
| 209 | *out << Trajectory.R.at(step).x[2] << endl; | 
|---|
| 210 | return true; | 
|---|
| 211 | } else | 
|---|
| 212 | return false; | 
|---|
| 213 | }; | 
|---|
| 214 |  | 
|---|
| 215 | /** Prints all bonds of this atom from given global lists. | 
|---|
| 216 | * \param *out stream to output to | 
|---|
| 217 | * \param *NumberOfBondsPerAtom array with number of bonds per atomic index | 
|---|
| 218 | * \param ***ListOfBondsPerAtom array per atomic index of array with pointer to bond | 
|---|
| 219 | * \return true - \a *out present, false - \a *out is NULL | 
|---|
| 220 | */ | 
|---|
| 221 | bool atom::OutputBondOfAtom(ofstream *out, int *NumberOfBondsPerAtom, bond ***ListOfBondsPerAtom) const | 
|---|
| 222 | { | 
|---|
| 223 | if (out != NULL) { | 
|---|
| 224 | *out << Verbose(4) << "Atom " << Name << "/" << nr << " with " << NumberOfBondsPerAtom[nr] << " bonds: "; | 
|---|
| 225 | int TotalDegree = 0; | 
|---|
| 226 | for (int j=0;j<NumberOfBondsPerAtom[nr];j++) { | 
|---|
| 227 | *out << *ListOfBondsPerAtom[nr][j] << "\t"; | 
|---|
| 228 | TotalDegree += ListOfBondsPerAtom[nr][j]->BondDegree; | 
|---|
| 229 | } | 
|---|
| 230 | *out << " -- TotalDegree: " << TotalDegree << endl; | 
|---|
| 231 | return true; | 
|---|
| 232 | } else | 
|---|
| 233 | return false; | 
|---|
| 234 | }; | 
|---|
| 235 |  | 
|---|
| 236 | ostream & operator << (ostream &ost, const atom &a) | 
|---|
| 237 | { | 
|---|
| 238 | ost << "[" << a.Name << "|" << &a << "]"; | 
|---|
| 239 | return ost; | 
|---|
| 240 | }; | 
|---|
| 241 |  | 
|---|
| 242 | ostream & atom::operator << (ostream &ost) | 
|---|
| 243 | { | 
|---|
| 244 | ost << "[" << Name << "|" << this << "]"; | 
|---|
| 245 | return ost; | 
|---|
| 246 | }; | 
|---|
| 247 |  | 
|---|
| 248 | /** Compares the indices of \a this atom with a given \a ptr. | 
|---|
| 249 | * \param ptr atom to compare index against | 
|---|
| 250 | * \return true - this one's is smaller, false - not | 
|---|
| 251 | */ | 
|---|
| 252 | bool atom::Compare(const atom &ptr) | 
|---|
| 253 | { | 
|---|
| 254 | if (nr < ptr.nr) | 
|---|
| 255 | return true; | 
|---|
| 256 | else | 
|---|
| 257 | return false; | 
|---|
| 258 | }; | 
|---|
| 259 |  | 
|---|
| 260 | bool operator < (atom &a, atom &b) | 
|---|
| 261 | { | 
|---|
| 262 | return a.Compare(b); | 
|---|
| 263 | }; | 
|---|
| 264 |  | 
|---|