Permalink
Browse files

core: Use generic rotation code in thermostat.

  • Loading branch information...
fweik committed Nov 7, 2018
1 parent 8ef71c2 commit 9cdfe22f02e33cfde93b2caa3e5445f34f8d4915
Showing with 6 additions and 67 deletions.
  1. +2 −1 src/core/electrostatics_magnetostatics/maggs.cpp
  2. +4 −66 src/core/thermostat.hpp
@@ -54,7 +54,8 @@
#include "integrate.hpp"
#include "nonbonded_interactions/nonbonded_interaction_data.hpp"
#include "particle_data.hpp"
#include "thermostat.hpp"
#include "debug.hpp"
#include <cstdio>
#include <cstdlib>
#include <cstring>
@@ -30,6 +30,7 @@
#include "integrate.hpp"
#include "particle_data.hpp"
#include "random.hpp"
#include "rotation.hpp"
#include "Vector.hpp"
@@ -114,66 +115,6 @@ void thermo_heat_up();
/** pendant to \ref thermo_heat_up */
void thermo_cool_down();
#ifdef ROTATION
inline void thermo_define_rotation_matrix(Particle *p, double A[9]) {
double q0q0 = p->r.quat[0];
q0q0 *= q0q0;
double q1q1 = p->r.quat[1];
q1q1 *= q1q1;
double q2q2 = p->r.quat[2];
q2q2 *= q2q2;
double q3q3 = p->r.quat[3];
q3q3 *= q3q3;
A[0 + 3 * 0] = q0q0 + q1q1 - q2q2 - q3q3;
A[1 + 3 * 1] = q0q0 - q1q1 + q2q2 - q3q3;
A[2 + 3 * 2] = q0q0 - q1q1 - q2q2 + q3q3;
A[0 + 3 * 1] =
2 * (p->r.quat[1] * p->r.quat[2] + p->r.quat[0] * p->r.quat[3]);
A[0 + 3 * 2] =
2 * (p->r.quat[1] * p->r.quat[3] - p->r.quat[0] * p->r.quat[2]);
A[1 + 3 * 0] =
2 * (p->r.quat[1] * p->r.quat[2] - p->r.quat[0] * p->r.quat[3]);
A[1 + 3 * 2] =
2 * (p->r.quat[2] * p->r.quat[3] + p->r.quat[0] * p->r.quat[1]);
A[2 + 3 * 0] =
2 * (p->r.quat[1] * p->r.quat[3] + p->r.quat[0] * p->r.quat[2]);
A[2 + 3 * 1] =
2 * (p->r.quat[2] * p->r.quat[3] - p->r.quat[0] * p->r.quat[1]);
}
inline void thermo_convert_forces_body_to_space(Particle *p, double *force) {
double A[9];
thermo_define_rotation_matrix(p, A);
force[0] = A[0 + 3 * 0] * p->f.f[0] + A[1 + 3 * 0] * p->f.f[1] +
A[2 + 3 * 0] * p->f.f[2];
force[1] = A[0 + 3 * 1] * p->f.f[0] + A[1 + 3 * 1] * p->f.f[1] +
A[2 + 3 * 1] * p->f.f[2];
force[2] = A[0 + 3 * 2] * p->f.f[0] + A[1 + 3 * 2] * p->f.f[1] +
A[2 + 3 * 2] * p->f.f[2];
}
inline void thermo_convert_vel_space_to_body(Particle *p,
const Vector3d &vel_space,
Vector3d &vel_body) {
double A[9];
thermo_define_rotation_matrix(p, A);
vel_body[0] = A[0 + 3 * 0] * vel_space[0] + A[0 + 3 * 1] * vel_space[1] +
A[0 + 3 * 2] * vel_space[2];
vel_body[1] = A[1 + 3 * 0] * vel_space[0] + A[1 + 3 * 1] * vel_space[1] +
A[1 + 3 * 2] * vel_space[2];
vel_body[2] = A[2 + 3 * 0] * vel_space[0] + A[2 + 3 * 1] * vel_space[1] +
A[2 + 3 * 2] * vel_space[2];
}
#endif // ROTATION
#ifdef NPT
/** add velocity-dependent noise and friction for NpT-sims to the particle's
velocity
@@ -279,10 +220,8 @@ inline void friction_thermo_langevin(Particle *p) {
(langevin_pref1_temp[1] != langevin_pref1_temp[2]) ||
(langevin_pref2_temp[0] != langevin_pref2_temp[1]) ||
(langevin_pref2_temp[1] != langevin_pref2_temp[2]);
Vector3d velocity_body = {0.0, 0.0, 0.0};
if (aniso_flag) {
thermo_convert_vel_space_to_body(p, velocity, velocity_body);
}
auto const velocity_body =
(aniso_flag) ? convert_vector_space_to_body(*p, velocity) : Vector3d{};
#endif
// Do the actual thermostatting
@@ -326,9 +265,8 @@ inline void friction_thermo_langevin(Particle *p) {
#ifdef PARTICLE_ANISOTROPY
if (aniso_flag) {
double particle_force[3] = {0.0, 0.0, 0.0};
auto const particle_force = convert_vector_body_to_space(*p, p->f.f);
thermo_convert_forces_body_to_space(p, particle_force);
for (int j = 0; j < 3; j++) {
#ifdef EXTERNAL_FORCES
if (!(p->p.ext_flag & COORD_FIXED(j)))

0 comments on commit 9cdfe22

Please sign in to comment.