Changes in src/atom_trajectoryparticle.cpp [a3fded:112b09]
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
src/atom_trajectoryparticle.cpp
ra3fded r112b09 5 5 * Author: heber 6 6 */ 7 8 #include "Helpers/MemDebug.hpp" 7 9 8 10 #include "atom.hpp" … … 13 15 #include "log.hpp" 14 16 #include "parser.hpp" 15 #include "ThermoStatContainer.hpp"16 17 #include "verbose.hpp" 17 18 … … 196 197 void TrajectoryParticle::Thermostat_Langevin(int Step, gsl_rng * r, double *ekin, config *configuration) 197 198 { 198 double sigma = sqrt(configuration->T hermostats->TargetTemp/type->mass); // sigma = (k_b T)/m (Hartree/atomicmass = atomiclength/atomictime)199 double sigma = sqrt(configuration->TargetTemp/type->mass); // sigma = (k_b T)/m (Hartree/atomicmass = atomiclength/atomictime) 199 200 Vector &U = Trajectory.U.at(Step); 200 201 if (FixedIon == 0) { // even FixedIon moves, only not by other's forces 201 202 // throw a dice to determine whether it gets hit by a heat bath particle 202 if (((((rand()/(double)RAND_MAX))*configuration->T hermostats->TempFrequency) < 1.)) {203 if (((((rand()/(double)RAND_MAX))*configuration->TempFrequency) < 1.)) { 203 204 DoLog(3) && (Log() << Verbose(3) << "Particle " << *this << " was hit (sigma " << sigma << "): " << sqrt(U[0]*U[0]+U[1]*U[1]+U[2]*U[2]) << " -> "); 204 205 // pick three random numbers from a Boltzmann distribution around the desired temperature T for each momenta axis … … 224 225 if (FixedIon == 0) { // even FixedIon moves, only not by other's forces 225 226 for (int d=0; d<NDIM; d++) { 226 U[d] *= sqrt(1+(configuration->Deltat/configuration->T hermostats->TempFrequency)*(ScaleTempFactor-1));227 U[d] *= sqrt(1+(configuration->Deltat/configuration->TempFrequency)*(ScaleTempFactor-1)); 227 228 *ekin += 0.5*type->mass * U[d]*U[d]; 228 229 } … … 254 255 if (FixedIon == 0) { // even FixedIon moves, only not by other's forces 255 256 for (int d=0; d<NDIM; d++) { 256 U[d] += configuration->Deltat/type->mass * (configuration-> Thermostats->alpha * (U[d] * type->mass));257 U[d] += configuration->Deltat/type->mass * (configuration->alpha * (U[d] * type->mass)); 257 258 *ekin += (0.5*type->mass) * U[d]*U[d]; 258 259 }
Note:
See TracChangeset
for help on using the changeset viewer.