From 74e31cccb7554c3ed4836adbc9fb4726fb3d33ea Mon Sep 17 00:00:00 2001 From: HanatoK Date: Thu, 30 May 2024 13:54:10 -0500 Subject: [PATCH] Initial implementation of the Cremer-Pople parameters --- namd/colvars/src/Makefile.namd | 3 +- src/colvar.cpp | 4 + src/colvar.h | 5 + src/colvarcomp_pucker.cpp | 235 +++++++++++++++++++++++++++++++++ src/colvarcomp_pucker.h | 62 +++++++++ 5 files changed, 308 insertions(+), 1 deletion(-) create mode 100644 src/colvarcomp_pucker.cpp create mode 100644 src/colvarcomp_pucker.h diff --git a/namd/colvars/src/Makefile.namd b/namd/colvars/src/Makefile.namd index 636f29681..b53f5ac10 100644 --- a/namd/colvars/src/Makefile.namd +++ b/namd/colvars/src/Makefile.namd @@ -39,5 +39,6 @@ COLVARSLIB = \ $(DSTDIR)/colvars_memstream.o \ $(DSTDIR)/colvartypes.o \ $(DSTDIR)/colvarvalue.o \ - $(DSTDIR)/nr_jacobi.o + $(DSTDIR)/nr_jacobi.o \ + $(DSTDIR)/colvarcomp_pucker.o diff --git a/src/colvar.cpp b/src/colvar.cpp index 58eb87fd0..70c1e8a87 100644 --- a/src/colvar.cpp +++ b/src/colvar.cpp @@ -17,6 +17,7 @@ #include "colvarvalue.h" #include "colvarparse.h" #include "colvarcomp.h" +#include "colvarcomp_pucker.h" #include "colvar.h" #include "colvarbias.h" #include "colvars_memstream.h" @@ -881,6 +882,9 @@ void colvar::define_component_types() add_component_type("euler phi angle of the optimal orientation", "eulerPhi"); add_component_type("euler psi angle of the optimal orientation", "eulerPsi"); add_component_type("euler theta angle of the optimal orientation", "eulerTheta"); + add_component_type("Cremer-Pople parameter Q", "cpQ"); + add_component_type("Cremer-Pople parameter theta", "cptheta"); + add_component_type("Cremer-Pople parameter phi", "cpphi"); #ifdef LEPTON add_component_type("CV with support of the Lepton custom function", "customColvar"); diff --git a/src/colvar.h b/src/colvar.h index 443e1e4bd..d16d960b5 100644 --- a/src/colvar.h +++ b/src/colvar.h @@ -641,6 +641,11 @@ class colvar : public colvarparse, public colvardeps { class cartesian; class orientation; + // Cremer-Pople parameters + class cpQ; + class cptheta; + class cpphi; + // components that do not handle any atoms directly class map_total; diff --git a/src/colvarcomp_pucker.cpp b/src/colvarcomp_pucker.cpp new file mode 100644 index 000000000..3abe309c4 --- /dev/null +++ b/src/colvarcomp_pucker.cpp @@ -0,0 +1,235 @@ +#include "colvarcomp_pucker.h" +#include "colvarmodule.h" +#include + +struct cpABC { + cvm::real A; + cvm::real B; + cvm::real C; +}; + +struct cpABC_grad { + std::vector& dA_dr; + std::vector& dB_dr; + std::vector& dC_dr; +}; + +cpABC calc_cpABC(const cvm::atom_group& r, cpABC_grad* grad = nullptr) { + // Calculate the center of geometry + cvm::rvector cog{0, 0, 0}; + for (auto it = r.begin(); it != r.end(); ++it) cog += (it->pos); + cog /= r.size(); + // Calculate R + std::vector R; + for (auto it = r.begin(); it != r.end(); ++it) R.push_back(it->pos - cog); + // Calculate R' and R'' + cvm::rvector Rp{0, 0, 0}; + cvm::rvector Rpp{0, 0, 0}; + // Save the sine and cosine factors for derivative calculations + std::vector sin_f1, cos_f1; + for (size_t i = 0; i < r.size(); ++i) { + const cvm::real factor = 2.0 * M_PI * i / r.size(); + const cvm::real sin_f = std::sin(factor); + const cvm::real cos_f = std::cos(factor); + Rp += sin_f * R[i]; + Rpp += cos_f * R[i]; + if (grad) { + sin_f1.push_back(sin_f); + cos_f1.push_back(cos_f); + } + } + // Calculate the normal vector + const cvm::rvector cross = cvm::rvector::outer(Rp, Rpp); + const cvm::real denominator = cross.norm(); + const cvm::rvector n_hat = cross / denominator; + // Project R_j onto the normal vector + std::vector z; + for (size_t i = 0; i < r.size(); ++i) { + z.push_back(R[i] * n_hat); + } + // Calculate A, B and C from z_j + cpABC result{0, 0, 0}; + // Save the sine and cosine factors for derivative calculations + std::vector sin_f2, cos_f2; + for (size_t i = 0; i < r.size(); ++i) { + const cvm::real factor = 2.0 * M_PI * 2.0 * i / r.size(); + const cvm::real sin_f = std::sin(factor); + const cvm::real cos_f = std::cos(factor); + result.A += z[i] * sin_f; + result.B += z[i] * cos_f; + result.C += z[i] * ((-2.0) * (i % 2) + 1.0); + if (grad) { + sin_f2.push_back(sin_f); + cos_f2.push_back(cos_f); + } + } + // Gradients of ABC with respect to r + if (grad) { + grad->dA_dr.assign(r.size(), cvm::rvector{0, 0, 0}); + grad->dB_dr.assign(r.size(), cvm::rvector{0, 0, 0}); + grad->dC_dr.assign(r.size(), cvm::rvector{0, 0, 0}); + const cvm::real tmp2 = 1.0 - 1.0 / r.size(); + const cvm::real tmp3 = -1.0 / r.size(); + const cvm::real one_denom = 1.0 / denominator; + const cvm::real one_denom_sq = one_denom * one_denom; + std::vector triples; + for (size_t k = 0; k < r.size(); ++k) { + triples.push_back(R[k] * cross); + } + for (size_t j = 0; j < r.size(); ++j) { + // ∂R'/∂rj + cvm::real dRp_drj_f = 0; + cvm::real dRpp_drj_f = 0; + for (size_t k = 0; k < r.size(); ++k) { + if (j == k) { + dRp_drj_f += sin_f1[k] * tmp2; + dRpp_drj_f += cos_f1[k] * tmp2; + } else { + dRp_drj_f += sin_f1[k] * tmp3; + dRpp_drj_f += cos_f1[k] * tmp3; + } + } + // ∂R'/∂rj × R'' + const cvm::rvector dRp_drjx_times_Rpp = cvm::rvector::outer({dRp_drj_f, 0, 0}, Rpp); + const cvm::rvector dRp_drjy_times_Rpp = cvm::rvector::outer({0, dRp_drj_f, 0}, Rpp); + const cvm::rvector dRp_drjz_times_Rpp = cvm::rvector::outer({0, 0, dRp_drj_f}, Rpp); + // R' × ∂R''/∂rj + const cvm::rvector Rp_times_dRpp_drjx = cvm::rvector::outer(Rp, {dRpp_drj_f, 0, 0}); + const cvm::rvector Rp_times_dRpp_drjy = cvm::rvector::outer(Rp, {0, dRpp_drj_f, 0}); + const cvm::rvector Rp_times_dRpp_drjz = cvm::rvector::outer(Rp, {0, 0, dRpp_drj_f}); + // Derivative of the norm + const cvm::real dnorm_dx = one_denom * ((cross * dRp_drjx_times_Rpp) + (Rp_times_dRpp_drjx * cross)); + const cvm::real dnorm_dy = one_denom * ((cross * dRp_drjy_times_Rpp) + (Rp_times_dRpp_drjy * cross)); + const cvm::real dnorm_dz = one_denom * ((cross * dRp_drjz_times_Rpp) + (Rp_times_dRpp_drjz * cross)); + // ABC derivatives with respect to r + for (size_t k = 0; k < r.size(); ++k) { + cvm::real dtriple_dx = 0; + cvm::real dtriple_dy = 0; + cvm::real dtriple_dz = 0; + if (j == k) { + dtriple_dx = tmp2 * cross.x + R[k] * dRp_drjx_times_Rpp + R[k] * Rp_times_dRpp_drjx; + dtriple_dy = tmp2 * cross.y + R[k] * dRp_drjy_times_Rpp + R[k] * Rp_times_dRpp_drjy; + dtriple_dz = tmp2 * cross.z + R[k] * dRp_drjz_times_Rpp + R[k] * Rp_times_dRpp_drjz; + } else { + dtriple_dx = tmp3 * cross.x + R[k] * dRp_drjx_times_Rpp + R[k] * Rp_times_dRpp_drjx; + dtriple_dy = tmp3 * cross.y + R[k] * dRp_drjy_times_Rpp + R[k] * Rp_times_dRpp_drjy; + dtriple_dz = tmp3 * cross.z + R[k] * dRp_drjz_times_Rpp + R[k] * Rp_times_dRpp_drjz; + } + // Derivative of z_j wrt r_j + const cvm::rvector dzk_drj{one_denom_sq * (dtriple_dx * denominator - dnorm_dx * triples[k]), + one_denom_sq * (dtriple_dy * denominator - dnorm_dy * triples[k]), + one_denom_sq * (dtriple_dz * denominator - dnorm_dz * triples[k])}; + grad->dA_dr[j] += sin_f2[k] * dzk_drj; + grad->dB_dr[j] += cos_f2[k] * dzk_drj; + grad->dB_dr[j] += ((-2.0) * (k % 2) + 1.0) * dzk_drj; + } + } + } + return result; +} + +colvar::cpQ::cpQ() { + set_function_type("cpQ"); + x.type(colvarvalue::type_scalar); + enable(f_cvc_explicit_gradient); +} + +int colvar::cpQ::init(const std::string& conf) { + int error_code = cvc::init(conf); + + atoms = parse_group(conf, "atoms"); + if (!atoms || atoms->size() != 6) { + return error_code | COLVARS_INPUT_ERROR; + } + return error_code; +} + +void colvar::cpQ::calc_value() { + cpABC_grad grad{dA_dr, dB_dr, dC_dr}; + cpABC result = calc_cpABC(*atoms, &grad); + A = result.A; + B = result.B; + C = result.C; + x.real_value = std::sqrt((2.0 * A * A + 2.0 * B * B + C * C) / 6.0); +} + +void colvar::cpQ::calc_gradients() { + for (size_t ia = 0; ia < atoms->size(); ia++) { + (*atoms)[ia].grad = (2.0 * A * dA_dr[ia] + 2.0 * B * dB_dr[ia] + C * dC_dr[ia]) / (6.0 * x.real_value); + } +} + +colvar::cptheta::cptheta() { + set_function_type("cptheta"); + init_as_angle(); + enable(f_cvc_explicit_gradient); +} + +int colvar::cptheta::init(const std::string& conf) { + int error_code = cvc::init(conf); + + atoms = parse_group(conf, "atoms"); + if (!atoms || atoms->size() != 6) { + return error_code | COLVARS_INPUT_ERROR; + } + return error_code; +} + +void colvar::cptheta::calc_value() { + cpABC_grad grad{dA_dr, dB_dr, dC_dr}; + cpABC result = calc_cpABC(*atoms, &grad); + A = result.A; + B = result.B; + C = result.C; + x.real_value = 180.0 / M_PI * std::acos(std::sqrt(C * C / (2.0 * A * A + 2.0 * B * B + C * C))); +} + +void colvar::cptheta::calc_gradients() { + const cvm::real tmp1 = 2.0 * (A * A + B * B); + const cvm::real tmp2 = tmp1 + C * C; + const cvm::real tmp3 = -180.0 / M_PI * std::sqrt(tmp2 / tmp1) * std::sqrt(tmp2 / (4.0 * C * C)); + for (size_t ia = 0; ia < atoms->size(); ia++) { + const cvm::rvector tmp4 = 2.0 * C * tmp2 * dC_dr[ia]; + const cvm::rvector tmp5 = C * C * (4.0 * A * dA_dr[ia] + 4.0 * B * dB_dr[ia] + 2.0 * C * dC_dr[ia]); + (*atoms)[ia].grad = tmp3 * (tmp4 - tmp5) / (tmp2 * tmp2); + } +} + +colvar::cpphi::cpphi() { + set_function_type("cpphi"); + x.type(colvarvalue::type_scalar); + provide(f_cvc_periodic); + enable(f_cvc_periodic); + period = 360.0; + init_scalar_boundaries(0, 360.0); + enable(f_cvc_explicit_gradient); +} + +int colvar::cpphi::init(const std::string& conf) { + int error_code = cvc::init(conf); + + atoms = parse_group(conf, "atoms"); + if (!atoms || atoms->size() != 6) { + return error_code | COLVARS_INPUT_ERROR; + } + return error_code; +} + +void colvar::cpphi::calc_value() { + cpABC_grad grad{dA_dr, dB_dr, dC_dr}; + cpABC result = calc_cpABC(*atoms, &grad); + A = result.A; + B = result.B; + C = result.C; + x.real_value = 180.0 / M_PI * std::atan2(-A, B); + if (x.real_value < 0) { + x.real_value += 360.0; + } +} + +void colvar::cpphi::calc_gradients() { + const cvm::real factor = 180.0 / M_PI / (A * A + B * B); + for (size_t ia = 0; ia < atoms->size(); ia++) { + (*atoms)[ia].grad = factor * (-dA_dr[ia] * B + dB_dr[ia] * A); + } +} diff --git a/src/colvarcomp_pucker.h b/src/colvarcomp_pucker.h new file mode 100644 index 000000000..e4e8586e6 --- /dev/null +++ b/src/colvarcomp_pucker.h @@ -0,0 +1,62 @@ +// -*- c++ -*- +#ifndef COLVARCOMP_PUCKER_H +#define COLVARCOMP_PUCKER_H + +#include "colvarcomp.h" + +class colvar::cpQ : public colvar::cvc +{ +protected: + /// Atom group + cvm::atom_group *atoms = nullptr; + cvm::real A; + cvm::real B; + cvm::real C; + std::vector dA_dr; + std::vector dB_dr; + std::vector dC_dr; +public: + cpQ(); + virtual int init(std::string const &conf); + virtual void calc_value(); + virtual void calc_gradients(); +}; + +class colvar::cptheta : public colvar::cvc +{ +protected: + /// Atom group + cvm::atom_group *atoms = nullptr; + cvm::real A; + cvm::real B; + cvm::real C; + std::vector dA_dr; + std::vector dB_dr; + std::vector dC_dr; +public: + cptheta(); + virtual int init(std::string const &conf); + virtual void calc_value(); + virtual void calc_gradients(); +}; + +class colvar::cpphi : public colvar::cvc +{ +protected: + /// Atom group + cvm::atom_group *atoms = nullptr; + cvm::real A; + cvm::real B; + cvm::real C; + std::vector dA_dr; + std::vector dB_dr; + std::vector dC_dr; +public: + cpphi(); + virtual int init(std::string const &conf); + virtual void calc_value(); + virtual void calc_gradients(); +}; + +#endif // COLVARCOMP_PUCKER_H +