Changes in src/atom_trajectoryparticle.cpp [d74077:952f38]
- File:
-
- 1 edited
-
src/atom_trajectoryparticle.cpp (modified) (13 diffs)
Legend:
- Unmodified
- Added
- Removed
-
src/atom_trajectoryparticle.cpp
rd74077 r952f38 12 12 #include "config.hpp" 13 13 #include "element.hpp" 14 #include " info.hpp"15 #include " log.hpp"14 #include "Helpers/Info.hpp" 15 #include "Helpers/Log.hpp" 16 16 #include "parser.hpp" 17 17 #include "ThermoStatContainer.hpp" 18 #include " verbose.hpp"18 #include "Helpers/Verbose.hpp" 19 19 20 20 /** Constructor of class TrajectoryParticle. … … 38 38 { 39 39 for (int i=NDIM;i--;) 40 *temperature += getType()->mass * Trajectory.U.at(step)[i]* Trajectory.U.at(step)[i];40 *temperature += type->mass * Trajectory.U.at(step)[i]* Trajectory.U.at(step)[i]; 41 41 }; 42 42 … … 65 65 for(int d=0;d<NDIM;d++) { 66 66 Trajectory.U.at(Step)[d] -= CoGVelocity->at(d); 67 *ActualTemp += 0.5 * getType()->mass * Trajectory.U.at(Step)[d] * Trajectory.U.at(Step)[d];67 *ActualTemp += 0.5 * type->mass * Trajectory.U.at(Step)[d] * Trajectory.U.at(Step)[d]; 68 68 } 69 69 }; … … 113 113 Trajectory.R.at(NextStep)[d] = Trajectory.R.at(NextStep-1)[d]; 114 114 Trajectory.R.at(NextStep)[d] += configuration->Deltat*(Trajectory.U.at(NextStep-1)[d]); // s(t) = s(0) + v * deltat + 1/2 a * deltat^2 115 Trajectory.R.at(NextStep)[d] += 0.5*configuration->Deltat*configuration->Deltat*(Trajectory.F.at(NextStep)[d]/ getType()->mass); // F = m * a and s =115 Trajectory.R.at(NextStep)[d] += 0.5*configuration->Deltat*configuration->Deltat*(Trajectory.F.at(NextStep)[d]/type->mass); // F = m * a and s = 116 116 } 117 117 // Update U 118 118 for (int d=0; d<NDIM; d++) { 119 119 Trajectory.U.at(NextStep)[d] = Trajectory.U.at(NextStep-1)[d]; 120 Trajectory.U.at(NextStep)[d] += configuration->Deltat * (Trajectory.F.at(NextStep)[d]+Trajectory.F.at(NextStep-1)[d]/ getType()->mass); // v = F/m * t120 Trajectory.U.at(NextStep)[d] += configuration->Deltat * (Trajectory.F.at(NextStep)[d]+Trajectory.F.at(NextStep-1)[d]/type->mass); // v = F/m * t 121 121 } 122 122 // Update R (and F) … … 137 137 void TrajectoryParticle::SumUpKineticEnergy( int Step, double *TotalMass, Vector *TotalVelocity ) const 138 138 { 139 *TotalMass += getType()->mass; // sum up total mass139 *TotalMass += type->mass; // sum up total mass 140 140 for(int d=0;d<NDIM;d++) { 141 TotalVelocity->at(d) += Trajectory.U.at(Step)[d]* getType()->mass;141 TotalVelocity->at(d) += Trajectory.U.at(Step)[d]*type->mass; 142 142 } 143 143 }; … … 154 154 for (int d=0; d<NDIM; d++) { 155 155 U[d] *= ScaleTempFactor; 156 *ekin += 0.5* getType()->mass * U[d]*U[d];156 *ekin += 0.5*type->mass * U[d]*U[d]; 157 157 } 158 158 }; … … 170 170 for (int d=0; d<NDIM; d++) { 171 171 *G += U[d] * F[d]; 172 *E += U[d]*U[d]* getType()->mass;172 *E += U[d]*U[d]*type->mass; 173 173 } 174 174 }; … … 185 185 if (FixedIon == 0) // even FixedIon moves, only not by other's forces 186 186 for (int d=0; d<NDIM; d++) { 187 U[d] += configuration->Deltat/ getType()->mass * ( (G_over_E) * (U[d]*getType()->mass) );188 *ekin += getType()->mass * U[d]*U[d];187 U[d] += configuration->Deltat/type->mass * ( (G_over_E) * (U[d]*type->mass) ); 188 *ekin += type->mass * U[d]*U[d]; 189 189 } 190 190 }; … … 198 198 void TrajectoryParticle::Thermostat_Langevin(int Step, gsl_rng * r, double *ekin, config *configuration) 199 199 { 200 double sigma = sqrt(configuration->Thermostats->TargetTemp/ getType()->mass); // sigma = (k_b T)/m (Hartree/atomicmass = atomiclength/atomictime)200 double sigma = sqrt(configuration->Thermostats->TargetTemp/type->mass); // sigma = (k_b T)/m (Hartree/atomicmass = atomiclength/atomictime) 201 201 Vector &U = Trajectory.U.at(Step); 202 202 if (FixedIon == 0) { // even FixedIon moves, only not by other's forces … … 211 211 } 212 212 for (int d=0; d<NDIM; d++) 213 *ekin += 0.5* getType()->mass * U[d]*U[d];213 *ekin += 0.5*type->mass * U[d]*U[d]; 214 214 } 215 215 }; … … 227 227 for (int d=0; d<NDIM; d++) { 228 228 U[d] *= sqrt(1+(configuration->Deltat/configuration->Thermostats->TempFrequency)*(ScaleTempFactor-1)); 229 *ekin += 0.5* getType()->mass * U[d]*U[d];229 *ekin += 0.5*type->mass * U[d]*U[d]; 230 230 } 231 231 } … … 241 241 if (FixedIon == 0) { // even FixedIon moves, only not by other's forces 242 242 for (int d=0; d<NDIM; d++) { 243 *delta_alpha += U[d]*U[d]* getType()->mass;243 *delta_alpha += U[d]*U[d]*type->mass; 244 244 } 245 245 } … … 256 256 if (FixedIon == 0) { // even FixedIon moves, only not by other's forces 257 257 for (int d=0; d<NDIM; d++) { 258 U[d] += configuration->Deltat/ getType()->mass * (configuration->Thermostats->alpha * (U[d] * getType()->mass));259 *ekin += (0.5* getType()->mass) * U[d]*U[d];258 U[d] += configuration->Deltat/type->mass * (configuration->Thermostats->alpha * (U[d] * type->mass)); 259 *ekin += (0.5*type->mass) * U[d]*U[d]; 260 260 } 261 261 } 262 262 }; 263 264 265 std::ostream & TrajectoryParticle::operator << (std::ostream &ost) const266 {267 ParticleInfo::operator<<(ost);268 ost << "," << getPosition();269 return ost;270 }271 272 std::ostream & operator << (std::ostream &ost, const TrajectoryParticle &a)273 {274 a.ParticleInfo::operator<<(ost);275 ost << "," << a.getPosition();276 return ost;277 }278
Note:
See TracChangeset
for help on using the changeset viewer.
