# Temperature of a Lennard-Jones system [closed]

#ifndef LJ_HPP
#define LJ_HPP

#include <vector>

#include "common.hpp"
#include "vec3.hpp"

class LennardJones
{
public:
real kB;
real epsilon;
real sigma;

public:
LennardJones()
{
kB = 0;
epsilon = 0;
sigma = 0;
}

LennardJones(double epsilon, double sigma, real kB) : epsilon(epsilon), sigma(sigma), kB(kB) {}

real getPotential(const Vec3& distance) const
{
real r_mag = std::sqrt(distance.x * distance.x + distance.y * distance.y + distance.z * distance.z);
real s_over_r = sigma / r_mag;
real s_over_r6 = pow(s_over_r, 6);
return epsilon * (s_over_r6 * s_over_r6 - 2.0 * s_over_r6);
}

real getPotentialAttractive(const Vec3& distance) const
{
real r_mag = std::sqrt(distance.x * distance.x + distance.y * distance.y + distance.z * distance.z);
real s_over_r = sigma / r_mag;
real s_over_r6 = pow(s_over_r, 6);
real attrPotential = (-2.0) * epsilon * s_over_r6;
return attrPotential;
}

real getPotentialRepulsive(const Vec3& distance) const
{
real r_mag = std::sqrt(distance.x * distance.x + distance.y * distance.y + distance.z * distance.z);
real s_over_r = sigma / r_mag;
real s_over_r12 = pow(s_over_r, 12);
real repPotential = epsilon * s_over_r12;
return repPotential;
}

Vec3 getForceApprox(const Vec3& distance) const
{
// Define a small distance for the derivative approximation
real dr = 1e-6;
Vec3 force;
std::vector<Vec3> r_plus_dr = { distance, distance, distance };
r_plus_dr[0].x += dr;
r_plus_dr[1].y += dr;
r_plus_dr[2].z += dr;

// The force is the negative derivative of the potential energy
force.x = -(getPotential(r_plus_dr[0]) - getPotential(distance)) / dr;
force.y = -(getPotential(r_plus_dr[1]) - getPotential(distance)) / dr;
force.z = -(getPotential(r_plus_dr[2]) - getPotential(distance)) / dr;
return force;
}

Vec3 getAcceleration(const Vec3& distance, double mass) const
{
double rSquared = distance.x * distance.x + distance.y * distance.y + distance.z * distance.z;
double r6 = std::pow(sigma * sigma / rSquared, 3);
double r12 = r6 * r6;

double magnitude = 24.0 * epsilon / mass * (2.0 * r12 - r6) / rSquared;

Vec3 acceleration(distance.x * magnitude, distance.y * magnitude, distance.z * magnitude);
return acceleration;
}

Vec3 getForce(Vec3 position) const
{
real r_mag = std::sqrt(position.x * position.x + position.y * position.y + position.z * position.z);
real s_over_r = sigma / r_mag;
real s_over_r6 = s_over_r * s_over_r * s_over_r * s_over_r * s_over_r * s_over_r;
real s_over_r12 = s_over_r6 * s_over_r6;
real factor = 24.0 * epsilon * (2.0 * s_over_r12 - s_over_r6) / (r_mag * r_mag * r_mag);

Vec3 force;
force.x = factor * position.x;
force.y = factor * position.y;
force.z = factor * position.z;
return force;
}

real getKineticEnergy(real mass, const Vec3& velocity) const
{
real vx = velocity.x;
real vy = velocity.y;
real vz = velocity.z;
return 0.5 * mass * (vx * vx + vy * vy + vz * vz);
}

real getTotalEnergy(real mass, const Vec3& distance, const Vec3& velocity) const
{
real potentialEnergy = getPotential(distance);
real kineticEnergy = getKineticEnergy(mass, velocity);
return potentialEnergy + kineticEnergy;
}

double getTemperature(real mass, const Vec3& velocity) const
{
double kineticEnergy = getKineticEnergy(mass, velocity);
double temperature = (2.0 * kineticEnergy) / (3.0 * kB);
return temperature;
}

void setTemperature(Vec3& velocity, real currentTemperature, real targetTemperature) const
{
real scalingFactor = std::sqrt(targetTemperature / currentTemperature);
velocity *= scalingFactor;
}
};
#endif // !LJ_HPP


Are the formulas used in getTemperature() and setTemperature() correct?

EDIT: My professor told me that these formulas are incorrect as the temperature is a feature of the system not individual particles.

So, I am now trying the following inside System class that holds a list of particles.

    double getTemperature()
{
int n_particles = particles_.size();
real currentTemperature = 0.0;
for (int i = 0; i < n_particles; i++)
{
real lenSq = particles_[i].velocity.magnitudeSquared();

currentTemperature += 0.5 * mass * lenSq / n_particles;
}
return currentTemperature;
}

void setTemperature(real targetTemperature)
{
int n_particles = particles_.size();

real currentTemperature = getTemperature();

real scalingFactor = std::sqrt(targetTemperature / currentTemperature);

for (int i = 0; i < n_particles; i++)
{
particles_[i].velocity *= scalingFactor;
}
}
$$$$
`
• What have you tried so far? Sep 18 at 14:32
• @MPIchael, Check the edit. Sep 18 at 15:12