[14de469] | 1 | /** \file atom.cpp
|
---|
[1907a7] | 2 | *
|
---|
[14de469] | 3 | * Function implementations for the class atom.
|
---|
[1907a7] | 4 | *
|
---|
[14de469] | 5 | */
|
---|
| 6 |
|
---|
[357fba] | 7 | #include "atom.hpp"
|
---|
[e41951] | 8 | #include "bond.hpp"
|
---|
[4a7776a] | 9 | #include "config.hpp"
|
---|
[f66195] | 10 | #include "element.hpp"
|
---|
[266237] | 11 | #include "lists.hpp"
|
---|
[29812d] | 12 | #include "memoryallocator.hpp"
|
---|
[ccd9f5] | 13 | #include "parser.hpp"
|
---|
[f66195] | 14 | #include "vector.hpp"
|
---|
[1907a7] | 15 |
|
---|
[14de469] | 16 | /************************************* Functions for class atom *************************************/
|
---|
| 17 |
|
---|
| 18 | /** Constructor of class atom.
|
---|
| 19 | */
|
---|
[1907a7] | 20 | atom::atom()
|
---|
[14de469] | 21 | {
|
---|
[357fba] | 22 | father = this; // generally, father is itself
|
---|
[14de469] | 23 | previous = NULL;
|
---|
| 24 | next = NULL;
|
---|
| 25 | Ancestor = NULL;
|
---|
| 26 | type = NULL;
|
---|
| 27 | sort = NULL;
|
---|
[943d02] | 28 | FixedIon = 0;
|
---|
[14de469] | 29 | GraphNr = -1;
|
---|
| 30 | ComponentNr = NULL;
|
---|
[683914] | 31 | IsCyclic = false;
|
---|
[14de469] | 32 | SeparationVertex = false;
|
---|
| 33 | LowpointNr = -1;
|
---|
[db942e] | 34 | AdaptiveOrder = 0;
|
---|
[362b0e] | 35 | MaxOrder = false;
|
---|
[357fba] | 36 | // set LCNode::Vector to our Vector
|
---|
| 37 | node = &x;
|
---|
[14de469] | 38 | };
|
---|
| 39 |
|
---|
[2319ed] | 40 | /** Constructor of class atom.
|
---|
| 41 | */
|
---|
| 42 | atom::atom(atom *pointer)
|
---|
| 43 | {
|
---|
| 44 | Name = NULL;
|
---|
| 45 | previous = NULL;
|
---|
| 46 | next = NULL;
|
---|
[89c8b2] | 47 | father = pointer; // generally, father is itself
|
---|
[2319ed] | 48 | Ancestor = NULL;
|
---|
[357fba] | 49 | GraphNr = -1;
|
---|
| 50 | ComponentNr = NULL;
|
---|
| 51 | IsCyclic = false;
|
---|
| 52 | SeparationVertex = false;
|
---|
| 53 | LowpointNr = -1;
|
---|
| 54 | AdaptiveOrder = 0;
|
---|
| 55 | MaxOrder = false;
|
---|
[2319ed] | 56 | type = pointer->type; // copy element of atom
|
---|
| 57 | x.CopyVector(&pointer->x); // copy coordination
|
---|
| 58 | v.CopyVector(&pointer->v); // copy velocity
|
---|
| 59 | FixedIon = pointer->FixedIon;
|
---|
| 60 | nr = -1;
|
---|
| 61 | sort = &nr;
|
---|
[89c8b2] | 62 | node = &x;
|
---|
[2319ed] | 63 | }
|
---|
| 64 |
|
---|
| 65 |
|
---|
[14de469] | 66 | /** Destructor of class atom.
|
---|
| 67 | */
|
---|
[1907a7] | 68 | atom::~atom()
|
---|
[14de469] | 69 | {
|
---|
[266237] | 70 | BondList::const_iterator Runner;
|
---|
| 71 | while (!ListOfBonds.empty()) {
|
---|
| 72 | Runner = ListOfBonds.begin();
|
---|
| 73 | removewithoutcheck(*Runner);
|
---|
| 74 | }
|
---|
| 75 | unlink(this);
|
---|
[a33931] | 76 | Free<int>(&ComponentNr, "atom::~atom: *ComponentNr");
|
---|
[c26f44] | 77 | Free<char>(&Name, "atom::~atom: *Name");
|
---|
[fcd7b6] | 78 | Trajectory.R.clear();
|
---|
| 79 | Trajectory.U.clear();
|
---|
| 80 | Trajectory.F.clear();
|
---|
[14de469] | 81 | };
|
---|
| 82 |
|
---|
| 83 |
|
---|
| 84 | /** Climbs up the father list until NULL, last is returned.
|
---|
| 85 | * \return true father, i.e. whose father points to itself, NULL if it could not be found or has none (added hydrogen)
|
---|
| 86 | */
|
---|
| 87 | atom *atom::GetTrueFather()
|
---|
| 88 | {
|
---|
| 89 | atom *walker = this;
|
---|
| 90 | do {
|
---|
| 91 | if (walker == walker->father) // top most father is the one that points on itself
|
---|
| 92 | break;
|
---|
| 93 | walker = walker->father;
|
---|
| 94 | } while (walker != NULL);
|
---|
| 95 | return walker;
|
---|
| 96 | };
|
---|
| 97 |
|
---|
[e65246] | 98 | /** Sets father to itself or its father in case of copying a molecule.
|
---|
| 99 | */
|
---|
| 100 | void atom::CorrectFather()
|
---|
| 101 | {
|
---|
| 102 | if (father->father == father) // same atom in copy's father points to itself
|
---|
| 103 | father = this; // set father to itself (copy of a whole molecule)
|
---|
| 104 | else
|
---|
| 105 | father = father->father; // set father to original's father
|
---|
| 106 |
|
---|
| 107 | };
|
---|
| 108 |
|
---|
| 109 | /** Check whether father is equal to given atom.
|
---|
| 110 | * \param *ptr atom to compare father to
|
---|
| 111 | * \param **res return value (only set if atom::father is equal to \a *ptr)
|
---|
| 112 | */
|
---|
| 113 | void atom::EqualsFather ( atom *ptr, atom **res )
|
---|
| 114 | {
|
---|
| 115 | if ( ptr == father )
|
---|
| 116 | *res = this;
|
---|
| 117 | };
|
---|
| 118 |
|
---|
[e9f8f9] | 119 | /** Checks whether atom is within the given box.
|
---|
| 120 | * \param offset offset to box origin
|
---|
| 121 | * \param *parallelepiped box matrix
|
---|
| 122 | * \return true - is inside, false - is not
|
---|
| 123 | */
|
---|
| 124 | bool atom::IsInParallelepiped(Vector offset, double *parallelepiped)
|
---|
| 125 | {
|
---|
| 126 | return (node->IsInParallelepiped(offset, parallelepiped));
|
---|
| 127 | };
|
---|
| 128 |
|
---|
[266237] | 129 | /** Counts the number of bonds weighted by bond::BondDegree.
|
---|
| 130 | * \param bonds times bond::BondDegree
|
---|
| 131 | */
|
---|
| 132 | int atom::CountBonds() const
|
---|
| 133 | {
|
---|
| 134 | int NoBonds = 0;
|
---|
| 135 | for (BondList::const_iterator Runner = ListOfBonds.begin(); Runner != ListOfBonds.end(); (++Runner))
|
---|
| 136 | NoBonds += (*Runner)->BondDegree;
|
---|
| 137 | return NoBonds;
|
---|
| 138 | };
|
---|
| 139 |
|
---|
[14de469] | 140 | /** Output of a single atom.
|
---|
| 141 | * \param ElementNo cardinal number of the element
|
---|
| 142 | * \param AtomNo cardinal number among these atoms of the same element
|
---|
| 143 | * \param *out stream to output to
|
---|
[1907a7] | 144 | * \param *comment commentary after '#' sign
|
---|
[e41951] | 145 | * \return true - \a *out present, false - \a *out is NULL
|
---|
[14de469] | 146 | */
|
---|
[fcd7b6] | 147 | bool atom::Output(ofstream *out, int ElementNo, int AtomNo, const char *comment) const
|
---|
[14de469] | 148 | {
|
---|
| 149 | if (out != NULL) {
|
---|
| 150 | *out << "Ion_Type" << ElementNo << "_" << AtomNo << "\t" << fixed << setprecision(9) << showpoint;
|
---|
[943d02] | 151 | *out << x.x[0] << "\t" << x.x[1] << "\t" << x.x[2];
|
---|
| 152 | *out << "\t" << FixedIon;
|
---|
| 153 | if (v.Norm() > MYEPSILON)
|
---|
| 154 | *out << "\t" << scientific << setprecision(6) << v.x[0] << "\t" << v.x[1] << "\t" << v.x[2] << "\t";
|
---|
[437922] | 155 | if (comment != NULL)
|
---|
| 156 | *out << " # " << comment << endl;
|
---|
[e9f8f9] | 157 | else
|
---|
| 158 | *out << " # molecule nr " << nr << endl;
|
---|
| 159 | return true;
|
---|
| 160 | } else
|
---|
| 161 | return false;
|
---|
| 162 | };
|
---|
[fcd7b6] | 163 | bool atom::Output(ofstream *out, int *ElementNo, int *AtomNo, const char *comment)
|
---|
[e9f8f9] | 164 | {
|
---|
| 165 | AtomNo[type->Z]++; // increment number
|
---|
| 166 | if (out != NULL) {
|
---|
| 167 | *out << "Ion_Type" << ElementNo[type->Z] << "_" << AtomNo[type->Z] << "\t" << fixed << setprecision(9) << showpoint;
|
---|
| 168 | *out << x.x[0] << "\t" << x.x[1] << "\t" << x.x[2];
|
---|
| 169 | *out << "\t" << FixedIon;
|
---|
| 170 | if (v.Norm() > MYEPSILON)
|
---|
| 171 | *out << "\t" << scientific << setprecision(6) << v.x[0] << "\t" << v.x[1] << "\t" << v.x[2] << "\t";
|
---|
| 172 | if (comment != NULL)
|
---|
| 173 | *out << " # " << comment << endl;
|
---|
[437922] | 174 | else
|
---|
| 175 | *out << " # molecule nr " << nr << endl;
|
---|
[14de469] | 176 | return true;
|
---|
| 177 | } else
|
---|
| 178 | return false;
|
---|
| 179 | };
|
---|
| 180 |
|
---|
| 181 | /** Output of a single atom as one lin in xyz file.
|
---|
| 182 | * \param *out stream to output to
|
---|
[e41951] | 183 | * \return true - \a *out present, false - \a *out is NULL
|
---|
[14de469] | 184 | */
|
---|
| 185 | bool atom::OutputXYZLine(ofstream *out) const
|
---|
| 186 | {
|
---|
| 187 | if (out != NULL) {
|
---|
| 188 | *out << type->symbol << "\t" << x.x[0] << "\t" << x.x[1] << "\t" << x.x[2] << "\t" << endl;
|
---|
| 189 | return true;
|
---|
| 190 | } else
|
---|
| 191 | return false;
|
---|
| 192 | };
|
---|
| 193 |
|
---|
[fcd7b6] | 194 | /** Output of a single atom as one lin in xyz file.
|
---|
| 195 | * \param *out stream to output to
|
---|
[e41951] | 196 | * \param *ElementNo array with ion type number in the config file this atom's element shall have
|
---|
| 197 | * \param *AtomNo array with atom number in the config file this atom shall have, is increase by one automatically
|
---|
| 198 | * \param step Trajectory time step to output
|
---|
| 199 | * \return true - \a *out present, false - \a *out is NULL
|
---|
[fcd7b6] | 200 | */
|
---|
| 201 | bool atom::OutputTrajectory(ofstream *out, int *ElementNo, int *AtomNo, int step) const
|
---|
| 202 | {
|
---|
| 203 | AtomNo[type->Z]++;
|
---|
| 204 | if (out != NULL) {
|
---|
| 205 | *out << "Ion_Type" << ElementNo[type->Z] << "_" << AtomNo[type->Z] << "\t" << fixed << setprecision(9) << showpoint;
|
---|
| 206 | *out << Trajectory.R.at(step).x[0] << "\t" << Trajectory.R.at(step).x[1] << "\t" << Trajectory.R.at(step).x[2];
|
---|
| 207 | *out << "\t" << FixedIon;
|
---|
| 208 | if (Trajectory.U.at(step).Norm() > MYEPSILON)
|
---|
| 209 | *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";
|
---|
| 210 | if (Trajectory.F.at(step).Norm() > MYEPSILON)
|
---|
| 211 | *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";
|
---|
| 212 | *out << "\t# Number in molecule " << nr << endl;
|
---|
| 213 | return true;
|
---|
| 214 | } else
|
---|
| 215 | return false;
|
---|
| 216 | };
|
---|
| 217 |
|
---|
[681a8a] | 218 | /** Output of a single atom as one lin in xyz file.
|
---|
| 219 | * \param *out stream to output to
|
---|
[e41951] | 220 | * \param step Trajectory time step to output
|
---|
| 221 | * \return true - \a *out present, false - \a *out is NULL
|
---|
[681a8a] | 222 | */
|
---|
| 223 | bool atom::OutputTrajectoryXYZ(ofstream *out, int step) const
|
---|
| 224 | {
|
---|
| 225 | if (out != NULL) {
|
---|
| 226 | *out << type->symbol << "\t";
|
---|
| 227 | *out << Trajectory.R.at(step).x[0] << "\t";
|
---|
| 228 | *out << Trajectory.R.at(step).x[1] << "\t";
|
---|
| 229 | *out << Trajectory.R.at(step).x[2] << endl;
|
---|
| 230 | return true;
|
---|
| 231 | } else
|
---|
| 232 | return false;
|
---|
| 233 | };
|
---|
| 234 |
|
---|
[266237] | 235 | /** Output graph info of this atom.
|
---|
| 236 | * \param *out output stream
|
---|
| 237 | */
|
---|
| 238 | void atom::OutputGraphInfo(ofstream *out) const
|
---|
| 239 | {
|
---|
| 240 | *out << Verbose(2) << "Atom " << Name << " is " << ((SeparationVertex) ? "a" : "not a") << " separation vertex, components are ";
|
---|
| 241 | OutputComponentNumber(out);
|
---|
| 242 | *out << " with Lowpoint " << LowpointNr << " and Graph Nr. " << GraphNr << "." << endl;
|
---|
| 243 | };
|
---|
| 244 |
|
---|
| 245 | /** Output a list of flags, stating whether the bond was visited or not.
|
---|
| 246 | * Note, we make use of the last entry of the ComponentNr always being -1 if allocated.
|
---|
| 247 | * \param *out output stream for debugging
|
---|
| 248 | */
|
---|
| 249 | void atom::OutputComponentNumber(ofstream *out) const
|
---|
| 250 | {
|
---|
| 251 | if (ComponentNr != NULL) {
|
---|
| 252 | for (int i=0; ComponentNr[i] != -1; i++)
|
---|
| 253 | *out << ComponentNr[i] << " ";
|
---|
| 254 | }
|
---|
| 255 | };
|
---|
| 256 |
|
---|
| 257 |
|
---|
| 258 | /** Prints all bonds of this atom with total degree.
|
---|
[e41951] | 259 | * \param *out stream to output to
|
---|
| 260 | * \return true - \a *out present, false - \a *out is NULL
|
---|
| 261 | */
|
---|
[266237] | 262 | bool atom::OutputBondOfAtom(ofstream *out) const
|
---|
[e41951] | 263 | {
|
---|
| 264 | if (out != NULL) {
|
---|
[49f802c] | 265 | #ifdef ADDHYDROGEN
|
---|
| 266 | if (type->Z != 1) { // regard only non-hydrogen
|
---|
| 267 | #endif
|
---|
[266237] | 268 | *out << Verbose(4) << "Atom " << Name << "/" << nr << " with " << ListOfBonds.size() << " bonds: ";
|
---|
[49f802c] | 269 | int TotalDegree = 0;
|
---|
[266237] | 270 | for (BondList::const_iterator Runner = ListOfBonds.begin(); Runner != ListOfBonds.end(); ++Runner) {
|
---|
| 271 | *out << **Runner << "\t";
|
---|
| 272 | TotalDegree += (*Runner)->BondDegree;
|
---|
[49f802c] | 273 | }
|
---|
| 274 | *out << " -- TotalDegree: " << TotalDegree << endl;
|
---|
| 275 | #ifdef ADDHYDROGEN
|
---|
[e41951] | 276 | }
|
---|
[49f802c] | 277 | #endif
|
---|
[e41951] | 278 | return true;
|
---|
| 279 | } else
|
---|
| 280 | return false;
|
---|
| 281 | };
|
---|
| 282 |
|
---|
[266237] | 283 | /** Output of atom::nr along with all bond partners.
|
---|
| 284 | * \param *AdjacencyFile output stream
|
---|
| 285 | */
|
---|
| 286 | void atom::OutputAdjacency(ofstream *AdjacencyFile) const
|
---|
| 287 | {
|
---|
| 288 | *AdjacencyFile << nr << "\t";
|
---|
| 289 | for (BondList::const_iterator Runner = ListOfBonds.begin(); Runner != ListOfBonds.end(); (++Runner))
|
---|
| 290 | *AdjacencyFile << (*Runner)->GetOtherAtom(this)->nr << "\t";
|
---|
| 291 | *AdjacencyFile << endl;
|
---|
| 292 | };
|
---|
| 293 |
|
---|
| 294 | /** Puts a given bond into atom::ListOfBonds.
|
---|
| 295 | * \param *Binder bond to insert
|
---|
| 296 | */
|
---|
| 297 | bool atom::RegisterBond(bond *Binder)
|
---|
| 298 | {
|
---|
| 299 | bool status = false;
|
---|
| 300 | if (Binder != NULL) {
|
---|
| 301 | if (Binder->Contains(this)) {
|
---|
| 302 | ListOfBonds.push_back(Binder);
|
---|
| 303 | status = true;
|
---|
| 304 | } else {
|
---|
| 305 | cout << Verbose(1) << "ERROR: " << *Binder << " does not contain " << *this << "." << endl;
|
---|
| 306 | }
|
---|
| 307 | } else {
|
---|
| 308 | cout << Verbose(1) << "ERROR: Binder is " << Binder << "." << endl;
|
---|
| 309 | }
|
---|
| 310 | return status;
|
---|
| 311 | };
|
---|
| 312 |
|
---|
| 313 | /** Removes a given bond from atom::ListOfBonds.
|
---|
| 314 | * \param *Binder bond to remove
|
---|
| 315 | */
|
---|
| 316 | bool atom::UnregisterBond(bond *Binder)
|
---|
| 317 | {
|
---|
| 318 | bool status = false;
|
---|
| 319 | if (Binder != NULL) {
|
---|
| 320 | if (Binder->Contains(this)) {
|
---|
| 321 | ListOfBonds.remove(Binder);
|
---|
| 322 | status = true;
|
---|
| 323 | } else {
|
---|
| 324 | cout << Verbose(1) << "ERROR: " << *Binder << " does not contain " << *this << "." << endl;
|
---|
| 325 | }
|
---|
| 326 | } else {
|
---|
| 327 | cout << Verbose(1) << "ERROR: Binder is " << Binder << "." << endl;
|
---|
| 328 | }
|
---|
| 329 | return status;
|
---|
| 330 | };
|
---|
| 331 |
|
---|
| 332 | /** Removes all bonds from atom::ListOfBonds.
|
---|
| 333 | * \note Does not do any memory de-allocation.
|
---|
| 334 | */
|
---|
| 335 | void atom::UnregisterAllBond()
|
---|
| 336 | {
|
---|
| 337 | ListOfBonds.clear();
|
---|
| 338 | };
|
---|
| 339 |
|
---|
| 340 | /** Corrects the bond degree by one at most if necessary.
|
---|
| 341 | * \param *out output stream for debugging
|
---|
| 342 | */
|
---|
| 343 | int atom::CorrectBondDegree(ofstream *out)
|
---|
| 344 | {
|
---|
| 345 | int NoBonds = 0;
|
---|
| 346 | int OtherNoBonds = 0;
|
---|
| 347 | int FalseBondDegree = 0;
|
---|
| 348 | atom *OtherWalker = NULL;
|
---|
| 349 | bond *CandidateBond = NULL;
|
---|
| 350 |
|
---|
| 351 | *out << Verbose(3) << "Walker " << *this << ": " << (int)this->type->NoValenceOrbitals << " > " << NoBonds << "?" << endl;
|
---|
| 352 | NoBonds = CountBonds();
|
---|
| 353 | if ((int)(type->NoValenceOrbitals) > NoBonds) { // we have a mismatch, check all bonding partners for mismatch
|
---|
| 354 | for (BondList::const_iterator Runner = ListOfBonds.begin(); Runner != ListOfBonds.end(); (++Runner)) {
|
---|
| 355 | OtherWalker = (*Runner)->GetOtherAtom(this);
|
---|
| 356 | OtherNoBonds = OtherWalker->CountBonds();
|
---|
| 357 | *out << Verbose(3) << "OtherWalker " << *OtherWalker << ": " << (int)OtherWalker->type->NoValenceOrbitals << " > " << NoBonds << "?" << endl;
|
---|
| 358 | if ((int)(OtherWalker->type->NoValenceOrbitals) > NoBonds) { // check if possible candidate
|
---|
| 359 | if ((CandidateBond == NULL) || (ListOfBonds.size() > OtherWalker->ListOfBonds.size())) { // pick the one with fewer number of bonds first
|
---|
| 360 | CandidateBond = (*Runner);
|
---|
| 361 | *out << Verbose(3) << "New candidate is " << *CandidateBond << "." << endl;
|
---|
| 362 | }
|
---|
| 363 | }
|
---|
| 364 | }
|
---|
| 365 | if ((CandidateBond != NULL)) {
|
---|
| 366 | CandidateBond->BondDegree++;
|
---|
| 367 | *out << Verbose(2) << "Increased bond degree for bond " << *CandidateBond << "." << endl;
|
---|
| 368 | } else {
|
---|
| 369 | *out << Verbose(2) << "Could not find correct degree for atom " << *this << "." << endl;
|
---|
| 370 | FalseBondDegree++;
|
---|
| 371 | }
|
---|
| 372 | }
|
---|
| 373 | return FalseBondDegree;
|
---|
| 374 | };
|
---|
| 375 |
|
---|
[321a11] | 376 | ostream & operator << (ostream &ost, const atom &a)
|
---|
[14de469] | 377 | {
|
---|
| 378 | ost << "[" << a.Name << "|" << &a << "]";
|
---|
| 379 | return ost;
|
---|
| 380 | };
|
---|
| 381 |
|
---|
[055861] | 382 | ostream & atom::operator << (ostream &ost)
|
---|
| 383 | {
|
---|
| 384 | ost << "[" << Name << "|" << this << "]";
|
---|
| 385 | return ost;
|
---|
| 386 | };
|
---|
| 387 |
|
---|
[14de469] | 388 | /** Compares the indices of \a this atom with a given \a ptr.
|
---|
| 389 | * \param ptr atom to compare index against
|
---|
| 390 | * \return true - this one's is smaller, false - not
|
---|
[1907a7] | 391 | */
|
---|
[321a11] | 392 | bool atom::Compare(const atom &ptr)
|
---|
[14de469] | 393 | {
|
---|
| 394 | if (nr < ptr.nr)
|
---|
| 395 | return true;
|
---|
| 396 | else
|
---|
| 397 | return false;
|
---|
| 398 | };
|
---|
| 399 |
|
---|
[4a7776a] | 400 | /** Extends the trajectory STL vector to the new size.
|
---|
| 401 | * Does nothing if \a MaxSteps is smaller than current size.
|
---|
| 402 | * \param MaxSteps
|
---|
| 403 | */
|
---|
| 404 | void atom::ResizeTrajectory(int MaxSteps)
|
---|
| 405 | {
|
---|
| 406 | if (Trajectory.R.size() <= (unsigned int)(MaxSteps)) {
|
---|
| 407 | //cout << "Increasing size for trajectory array of " << keyword << " to " << (MaxSteps+1) << "." << endl;
|
---|
| 408 | Trajectory.R.resize(MaxSteps+1);
|
---|
| 409 | Trajectory.U.resize(MaxSteps+1);
|
---|
| 410 | Trajectory.F.resize(MaxSteps+1);
|
---|
| 411 | }
|
---|
| 412 | };
|
---|
| 413 |
|
---|
| 414 | /** Copies a given trajectory step \a src onto another \a dest
|
---|
| 415 | * \param dest index of destination step
|
---|
| 416 | * \param src index of source step
|
---|
| 417 | */
|
---|
| 418 | void atom::CopyStepOnStep(int dest, int src)
|
---|
| 419 | {
|
---|
| 420 | if (dest == src) // self assignment check
|
---|
| 421 | return;
|
---|
| 422 |
|
---|
| 423 | for (int n=NDIM;n--;) {
|
---|
| 424 | Trajectory.R.at(dest).x[n] = Trajectory.R.at(src).x[n];
|
---|
| 425 | Trajectory.U.at(dest).x[n] = Trajectory.U.at(src).x[n];
|
---|
| 426 | Trajectory.F.at(dest).x[n] = Trajectory.F.at(src).x[n];
|
---|
| 427 | }
|
---|
| 428 | };
|
---|
| 429 |
|
---|
| 430 | /** Performs a velocity verlet update of the trajectory.
|
---|
| 431 | * Parameters are according to those in configuration class.
|
---|
| 432 | * \param NextStep index of sequential step to set
|
---|
| 433 | * \param *configuration pointer to configuration with parameters
|
---|
| 434 | * \param *Force matrix with forces
|
---|
| 435 | */
|
---|
| 436 | void atom::VelocityVerletUpdate(int NextStep, config *configuration, ForceMatrix *Force)
|
---|
| 437 | {
|
---|
| 438 | //a = configuration.Deltat*0.5/walker->type->mass; // (F+F_old)/2m = a and thus: v = (F+F_old)/2m * t = (F + F_old) * a
|
---|
| 439 | for (int d=0; d<NDIM; d++) {
|
---|
| 440 | Trajectory.F.at(NextStep).x[d] = -Force->Matrix[0][nr][d+5]*(configuration->GetIsAngstroem() ? AtomicLengthToAngstroem : 1.);
|
---|
| 441 | Trajectory.R.at(NextStep).x[d] = Trajectory.R.at(NextStep-1).x[d];
|
---|
| 442 | Trajectory.R.at(NextStep).x[d] += configuration->Deltat*(Trajectory.U.at(NextStep-1).x[d]); // s(t) = s(0) + v * deltat + 1/2 a * deltat^2
|
---|
| 443 | Trajectory.R.at(NextStep).x[d] += 0.5*configuration->Deltat*configuration->Deltat*(Trajectory.F.at(NextStep).x[d]/type->mass); // F = m * a and s = 0.5 * F/m * t^2 = F * a * t
|
---|
| 444 | }
|
---|
| 445 | // Update U
|
---|
| 446 | for (int d=0; d<NDIM; d++) {
|
---|
| 447 | Trajectory.U.at(NextStep).x[d] = Trajectory.U.at(NextStep-1).x[d];
|
---|
| 448 | Trajectory.U.at(NextStep).x[d] += configuration->Deltat * (Trajectory.F.at(NextStep).x[d]+Trajectory.F.at(NextStep-1).x[d]/type->mass); // v = F/m * t
|
---|
| 449 | }
|
---|
| 450 | // Update R (and F)
|
---|
| 451 | // out << "Integrated position&velocity of step " << (NextStep) << ": (";
|
---|
| 452 | // for (int d=0;d<NDIM;d++)
|
---|
| 453 | // out << Trajectory.R.at(NextStep).x[d] << " "; // next step
|
---|
| 454 | // out << ")\t(";
|
---|
| 455 | // for (int d=0;d<NDIM;d++)
|
---|
| 456 | // cout << Trajectory.U.at(NextStep).x[d] << " "; // next step
|
---|
| 457 | // out << ")" << endl;
|
---|
| 458 | };
|
---|
| 459 |
|
---|
| 460 | /** Sums up mass and kinetics.
|
---|
| 461 | * \param Step step to sum for
|
---|
| 462 | * \param *TotalMass pointer to total mass sum
|
---|
| 463 | * \param *TotalVelocity pointer to tota velocity sum
|
---|
| 464 | */
|
---|
| 465 | void atom::SumUpKineticEnergy( int Step, double *TotalMass, Vector *TotalVelocity )
|
---|
| 466 | {
|
---|
| 467 | *TotalMass += type->mass; // sum up total mass
|
---|
| 468 | for(int d=0;d<NDIM;d++) {
|
---|
| 469 | TotalVelocity->x[d] += Trajectory.U.at(Step).x[d]*type->mass;
|
---|
| 470 | }
|
---|
| 471 | };
|
---|
| 472 |
|
---|
[5034e1] | 473 | /** Outputs the current atom::AdaptiveOrder and atom::MaxOrder to \a *file.
|
---|
| 474 | * \param *file output stream
|
---|
| 475 | */
|
---|
| 476 | void atom::OutputOrder(ofstream *file)
|
---|
| 477 | {
|
---|
| 478 | *file << nr << "\t" << (int)AdaptiveOrder << "\t" << (int)MaxOrder << endl;
|
---|
[266237] | 479 | //cout << Verbose(2) << "Storing: " << nr << "\t" << (int)AdaptiveOrder << "\t" << (int)MaxOrder << "." << endl;
|
---|
[5034e1] | 480 | }
|
---|
| 481 |
|
---|
[49f802c] | 482 | /** Returns squared distance to a given vector.
|
---|
| 483 | * \param origin vector to calculate distance to
|
---|
| 484 | * \return distance squared
|
---|
| 485 | */
|
---|
| 486 | double atom::DistanceSquaredToVector(Vector &origin)
|
---|
| 487 | {
|
---|
| 488 | return origin.DistanceSquared(&x);
|
---|
| 489 | };
|
---|
| 490 |
|
---|
| 491 | /** Adds kinetic energy of this atom to given temperature value.
|
---|
| 492 | * \param *temperature add on this value
|
---|
| 493 | * \param step given step of trajectory to add
|
---|
| 494 | */
|
---|
| 495 | void atom::AddKineticToTemperature(double *temperature, int step) const
|
---|
| 496 | {
|
---|
| 497 | for (int i=NDIM;i--;)
|
---|
| 498 | *temperature += type->mass * Trajectory.U.at(step).x[i]* Trajectory.U.at(step).x[i];
|
---|
| 499 | };
|
---|
| 500 |
|
---|
| 501 | /** Returns distance to a given vector.
|
---|
| 502 | * \param origin vector to calculate distance to
|
---|
| 503 | * \return distance
|
---|
| 504 | */
|
---|
| 505 | double atom::DistanceToVector(Vector &origin)
|
---|
| 506 | {
|
---|
| 507 | return origin.Distance(&x);
|
---|
| 508 | };
|
---|
| 509 |
|
---|
[1907a7] | 510 | bool operator < (atom &a, atom &b)
|
---|
[14de469] | 511 | {
|
---|
| 512 | return a.Compare(b);
|
---|
| 513 | };
|
---|
| 514 |
|
---|
[ccd9f5] | 515 | /** Evaluates some constraint potential if atom moves from \a startstep at once to \endstep in trajectory.
|
---|
| 516 | * \param startstep trajectory begins at
|
---|
| 517 | * \param endstep trajectory ends at
|
---|
| 518 | * \param **PermutationMap if atom switches places with some other atom, there is no translation but a permutaton noted here (not in the trajectories of each).
|
---|
| 519 | * \param *Force Force matrix to store result in
|
---|
| 520 | */
|
---|
| 521 | void atom::EvaluateConstrainedForce(int startstep, int endstep, atom **PermutationMap, ForceMatrix *Force)
|
---|
| 522 | {
|
---|
| 523 | double constant = 10.;
|
---|
| 524 | atom *Sprinter = PermutationMap[nr];
|
---|
| 525 | // set forces
|
---|
| 526 | for (int i=NDIM;i++;)
|
---|
| 527 | Force->Matrix[0][nr][5+i] += 2.*constant*sqrt(Trajectory.R.at(startstep).Distance(&Sprinter->Trajectory.R.at(endstep)));
|
---|
| 528 | };
|
---|
[4a7776a] | 529 |
|
---|
| 530 | /** Correct velocity against the summed \a CoGVelocity for \a step.
|
---|
| 531 | * \param *ActualTemp sum up actual temperature meanwhile
|
---|
| 532 | * \param Step MD step in atom::Tracjetory
|
---|
| 533 | * \param *CoGVelocity remnant velocity (i.e. vector sum of all atom velocities)
|
---|
| 534 | */
|
---|
| 535 | void atom::CorrectVelocity(double *ActualTemp, int Step, Vector *CoGVelocity)
|
---|
| 536 | {
|
---|
| 537 | for(int d=0;d<NDIM;d++) {
|
---|
| 538 | Trajectory.U.at(Step).x[d] -= CoGVelocity->x[d];
|
---|
| 539 | *ActualTemp += 0.5 * type->mass * Trajectory.U.at(Step).x[d] * Trajectory.U.at(Step).x[d];
|
---|
| 540 | }
|
---|
| 541 | };
|
---|
| 542 |
|
---|
| 543 | /** Scales velocity of atom according to Woodcock thermostat.
|
---|
| 544 | * \param ScaleTempFactor factor to scale the velocities with (i.e. sqrt of energy scale factor)
|
---|
| 545 | * \param Step MD step to scale
|
---|
| 546 | * \param *ekin sum of kinetic energy
|
---|
| 547 | */
|
---|
| 548 | void atom::Thermostat_Woodcock(double ScaleTempFactor, int Step, double *ekin)
|
---|
| 549 | {
|
---|
| 550 | double *U = Trajectory.U.at(Step).x;
|
---|
| 551 | if (FixedIon == 0) // even FixedIon moves, only not by other's forces
|
---|
| 552 | for (int d=0; d<NDIM; d++) {
|
---|
| 553 | U[d] *= ScaleTempFactor;
|
---|
| 554 | *ekin += 0.5*type->mass * U[d]*U[d];
|
---|
| 555 | }
|
---|
| 556 | };
|
---|
| 557 |
|
---|
| 558 | /** Scales velocity of atom according to Gaussian thermostat.
|
---|
| 559 | * \param Step MD step to scale
|
---|
| 560 | * \param *G
|
---|
| 561 | * \param *E
|
---|
| 562 | */
|
---|
| 563 | void atom::Thermostat_Gaussian_init(int Step, double *G, double *E)
|
---|
| 564 | {
|
---|
| 565 | double *U = Trajectory.U.at(Step).x;
|
---|
| 566 | double *F = Trajectory.F.at(Step).x;
|
---|
| 567 | if (FixedIon == 0) // even FixedIon moves, only not by other's forces
|
---|
| 568 | for (int d=0; d<NDIM; d++) {
|
---|
| 569 | *G += U[d] * F[d];
|
---|
| 570 | *E += U[d]*U[d]*type->mass;
|
---|
| 571 | }
|
---|
| 572 | };
|
---|
| 573 |
|
---|
| 574 | /** Determines scale factors according to Gaussian thermostat.
|
---|
| 575 | * \param Step MD step to scale
|
---|
| 576 | * \param GE G over E ratio
|
---|
| 577 | * \param *ekin sum of kinetic energy
|
---|
| 578 | * \param *configuration configuration class with TempFrequency and TargetTemp
|
---|
| 579 | */
|
---|
| 580 | void atom::Thermostat_Gaussian_least_constraint(int Step, double G_over_E, double *ekin, config *configuration)
|
---|
| 581 | {
|
---|
| 582 | double *U = Trajectory.U.at(Step).x;
|
---|
| 583 | if (FixedIon == 0) // even FixedIon moves, only not by other's forces
|
---|
| 584 | for (int d=0; d<NDIM; d++) {
|
---|
| 585 | U[d] += configuration->Deltat/type->mass * ( (G_over_E) * (U[d]*type->mass) );
|
---|
| 586 | *ekin += type->mass * U[d]*U[d];
|
---|
| 587 | }
|
---|
| 588 | };
|
---|
| 589 |
|
---|
| 590 | /** Scales velocity of atom according to Langevin thermostat.
|
---|
| 591 | * \param Step MD step to scale
|
---|
| 592 | * \param *r random number generator
|
---|
| 593 | * \param *ekin sum of kinetic energy
|
---|
| 594 | * \param *configuration configuration class with TempFrequency and TargetTemp
|
---|
| 595 | */
|
---|
| 596 | void atom::Thermostat_Langevin(int Step, gsl_rng * r, double *ekin, config *configuration)
|
---|
| 597 | {
|
---|
| 598 | double sigma = sqrt(configuration->TargetTemp/type->mass); // sigma = (k_b T)/m (Hartree/atomicmass = atomiclength/atomictime)
|
---|
| 599 | double *U = Trajectory.U.at(Step).x;
|
---|
| 600 | if (FixedIon == 0) { // even FixedIon moves, only not by other's forces
|
---|
| 601 | // throw a dice to determine whether it gets hit by a heat bath particle
|
---|
| 602 | if (((((rand()/(double)RAND_MAX))*configuration->TempFrequency) < 1.)) {
|
---|
| 603 | cout << Verbose(3) << "Particle " << *this << " was hit (sigma " << sigma << "): " << sqrt(U[0]*U[0]+U[1]*U[1]+U[2]*U[2]) << " -> ";
|
---|
| 604 | // pick three random numbers from a Boltzmann distribution around the desired temperature T for each momenta axis
|
---|
| 605 | for (int d=0; d<NDIM; d++) {
|
---|
| 606 | U[d] = gsl_ran_gaussian (r, sigma);
|
---|
| 607 | }
|
---|
| 608 | cout << sqrt(U[0]*U[0]+U[1]*U[1]+U[2]*U[2]) << endl;
|
---|
| 609 | }
|
---|
| 610 | for (int d=0; d<NDIM; d++)
|
---|
| 611 | *ekin += 0.5*type->mass * U[d]*U[d];
|
---|
| 612 | }
|
---|
| 613 | };
|
---|
| 614 |
|
---|
| 615 | /** Scales velocity of atom according to Berendsen thermostat.
|
---|
| 616 | * \param Step MD step to scale
|
---|
| 617 | * \param ScaleTempFactor factor to scale energy (not velocity!) with
|
---|
| 618 | * \param *ekin sum of kinetic energy
|
---|
| 619 | * \param *configuration configuration class with TempFrequency and Deltat
|
---|
| 620 | */
|
---|
| 621 | void atom::Thermostat_Berendsen(int Step, double ScaleTempFactor, double *ekin, config *configuration)
|
---|
| 622 | {
|
---|
| 623 | double *U = Trajectory.U.at(Step).x;
|
---|
| 624 | if (FixedIon == 0) { // even FixedIon moves, only not by other's forces
|
---|
| 625 | for (int d=0; d<NDIM; d++) {
|
---|
| 626 | U[d] *= sqrt(1+(configuration->Deltat/configuration->TempFrequency)*(ScaleTempFactor-1));
|
---|
| 627 | *ekin += 0.5*type->mass * U[d]*U[d];
|
---|
| 628 | }
|
---|
| 629 | }
|
---|
| 630 | };
|
---|
| 631 |
|
---|
| 632 | /** Initializes current run of NoseHoover thermostat.
|
---|
| 633 | * \param Step MD step to scale
|
---|
| 634 | * \param *delta_alpha additional sum of kinetic energy on return
|
---|
| 635 | */
|
---|
| 636 | void atom::Thermostat_NoseHoover_init(int Step, double *delta_alpha)
|
---|
| 637 | {
|
---|
| 638 | double *U = Trajectory.U.at(Step).x;
|
---|
| 639 | if (FixedIon == 0) { // even FixedIon moves, only not by other's forces
|
---|
| 640 | for (int d=0; d<NDIM; d++) {
|
---|
| 641 | *delta_alpha += U[d]*U[d]*type->mass;
|
---|
| 642 | }
|
---|
| 643 | }
|
---|
| 644 | };
|
---|
| 645 |
|
---|
| 646 | /** Initializes current run of NoseHoover thermostat.
|
---|
| 647 | * \param Step MD step to scale
|
---|
| 648 | * \param *ekin sum of kinetic energy
|
---|
| 649 | * \param *configuration configuration class with TempFrequency and Deltat
|
---|
| 650 | */
|
---|
| 651 | void atom::Thermostat_NoseHoover_scale(int Step, double *ekin, config *configuration)
|
---|
| 652 | {
|
---|
| 653 | double *U = Trajectory.U.at(Step).x;
|
---|
| 654 | if (FixedIon == 0) { // even FixedIon moves, only not by other's forces
|
---|
| 655 | for (int d=0; d<NDIM; d++) {
|
---|
| 656 | U[d] += configuration->Deltat/type->mass * (configuration->alpha * (U[d] * type->mass));
|
---|
| 657 | *ekin += (0.5*type->mass) * U[d]*U[d];
|
---|
| 658 | }
|
---|
| 659 | }
|
---|
| 660 | };
|
---|